CrossSection.hpp

Go to the documentation of this file.

#ifndef CROSS_SECTION_HPP
#define CROSS_SECTION_HPP
 
#include <Rtypes.h>
#include <TMatrixD.h>
#include <TString.h>
#include <map>
#include <utility>
#include <vector>
 
class TCutG;
class TFile;
class TGraph;
class TH2F;
class TCanvas;
class TPad;
class TObject;
class TGraphErrors;
class TGraphAsymmErrors;
struct CrossSectionChannel;

class CrossSection {
public:
  Bool_t Run();

  static Bool_t ExitResidue(const TString &exit, Int_t z_beam, Int_t a_beam,
                            Int_t &z, Int_t &a);
  static TString Label(const CrossSectionChannel &ch);
  static TString SubtractedLabel(const CrossSectionChannel &ch);
 
private:
  struct TalysCurve {
    TString label;
    TGraph *axn;
  };
  struct Point {
    Int_t reac;
    Double_t e_in, e_out, e_eff, e_eff_lo, e_eff_hi; // E_cm [MeV]
    Double_t n_reac, n_denom;
    Double_t sub_frac, sub_frac_err;
    Double_t sigma, sys; // [mb]
    Double_t stat_lo, stat_hi;
    Double_t ErrLo() const; 
    Double_t ErrHi() const; 
  };
  struct ChannelResult {
    const CrossSectionChannel *ch;
    std::vector<TalysCurve> talys;
    std::vector<TalysCurve> talys_sub;
    std::vector<Point> points;
  };
  struct FitComponents {
    std::vector<TGraph *> graph;
    std::vector<TString> label;
    std::vector<TString> tag;
    Int_t n_keep = 0;
    Int_t Size() const { return Int_t(graph.size()); }
  };
  struct ScaleFit {
    std::vector<Double_t> scale;
    std::vector<Double_t> scale_err;
    std::vector<Double_t> scale_raw;
    std::vector<Double_t> scale_raw_err;
    std::vector<Int_t> free_idx;
    std::vector<Int_t> pinned;
    TMatrixD cov;
    Double_t chi2 = 0.0;
    Int_t nf = 0;
  };
  struct FitFrame {
    Double_t fx_lo = 1.0e9, fx_hi = -1.0e9, fy_lo = 1.0e9, fy_hi = -1.0e9;
    Double_t pull_max = 0.0;
  };
  struct FitPanels {
    TCanvas *c = nullptr;
    TPad *plot = nullptr;
    TPad *side = nullptr;
    TPad *top = nullptr;
    TPad *bottom = nullptr;
    Int_t plot_right_margin_px = 0;
    Int_t extra_width_px = 0;
    std::vector<std::pair<TObject *, TString>> entries;
  };
 
  Bool_t LoadCache();
  Bool_t LoadBeam();
 
  // Every residual-production graph of every model in
  // root_files/talys/talys_xs.root, keyed by model then (Z, A).
  void LoadTalys();
 
  // Sum the named exits for each TALYS model.
  std::vector<TalysCurve> ExitCurves(const CrossSectionChannel &ch,
                                     const std::vector<TString> &exits,
                                     const char *what) const;
  Bool_t RunChannel(const CrossSectionChannel &ch, ChannelResult &out);
 
  // One strip's point, or kFALSE (with a printed reason) when it has none.
  Bool_t Strip(const CrossSectionChannel &ch,
               const std::vector<TalysCurve> &talys,
               const std::vector<TalysCurve> &talys_sub, Int_t reac, Point &pt);
  void CompareReference(const ChannelResult &r) const;
 
  // The figure for these channels; name is the file's basename.
  void Draw(const std::vector<const ChannelResult *> &rs,
            const TString &name) const;
  // The second figure: one model's sum and per-exit curves, the sum with a
  // scale per exit fitted to the points, and a deviation panel.
  void DrawFit(const ChannelResult &r, const TString &name) const;
  // The model's per-exit curves for the fit figure, or kFALSE (with the
  // reason printed) when it has none of the channel's exits.
  Bool_t CollectFitComponents(const ChannelResult &r, Int_t m,
                              FitComponents &comp) const;
  // The weighted least-squares scale per exit, exits that go negative
  // Negative scales are pinned at zero; point_scale fits pre-subtraction data.
  void FitExitScales(const ChannelResult &r, const FitComponents &comp,
                     ScaleFit &fs,
                     const std::vector<Double_t> *point_scale = nullptr) const;
  // Compare subtracted-exit scales before and after correction.
  void PrintSubtractionCheck(const ChannelResult &r, const FitComponents &comp,
                             const ScaleFit &fs) const;
  // The "scaled to this work's points" block: the per-exit lines, the
  // free-scale correlations, chi2, and the pull table header.
  void PrintScaleFit(const CrossSectionChannel &ch, Int_t m,
                     const FitComponents &comp, const ScaleFit &fs,
                     Int_t np) const;
  // The scaled sum and its 3-sigma band on the model's grid, and the
  // unscaled sum; kFALSE (the caller returns) when the grid is empty.
  Bool_t BuildScaledCurves(const ChannelResult &r, const FitComponents &comp,
                           const ScaleFit &fs, TGraph *&sum, TGraph *&fit,
                           TGraphErrors *&band) const;
  // The measured and deviation graphs and each point's pull (printed as a
  // row); extends the frame.
  void BuildMeasuredGraphs(const ChannelResult &r, TGraph *fit,
                           TGraphAsymmErrors *&measured,
                           TGraphAsymmErrors *&deviation,
                           FitFrame &frame) const;
  // The published reference points, nullptr when the channel has none;
  // extends the frame.
  TGraphAsymmErrors *BuildPublishedGraph(const CrossSectionChannel &ch,
                                         FitFrame &frame) const;
  // The "#times" scale text and the legend's width in pixels from the
  // full label list.
  void FitLegendWidthPx(const ScaleFit &fs, const CrossSectionChannel &ch,
                        const FitComponents &comp, const TString &model_label,
                        Bool_t has_published, TString &scale_text,
                        Double_t &legend_width_px) const;
  // The canvas with its side-legend pad and the top panel: the frame, the
  // 3-sigma band, the scaled and unscaled sums, the per-exit components,
  // the published and measured points, and the preliminary stamp.
  void DrawFitTopPanel(const CrossSectionChannel &ch, const FitFrame &frame,
                       const TString &scale_text, Double_t legend_width_px,
                       TGraph *sum, TGraph *fit, TGraphErrors *band,
                       const FitComponents &comp, TGraphAsymmErrors *published,
                       TGraphAsymmErrors *measured, FitPanels &pads) const;
  // The pull panel: the zero and +-3 lines and the deviation points.
  void DrawPullPanel(const FitFrame &frame, TGraphAsymmErrors *deviation,
                     FitPanels &pads) const;
  // The "PRELIMINARY" stamp: red, bold, in the top-left corner just inside
  // the frame's margins; `pad` must be the current pad.
  static void StampPreliminary(const TPad *pad);
 
  static Long64_t ReadCount(TFile &f, const char *name, Bool_t &ok);
  static TGraph *Clipped(TGraph *g, Double_t e_lo, Double_t e_hi);
  static Double_t Interpolated(TGraph *g, Double_t e);
  static Double_t EffectiveEnergy(TGraph *axn, Double_t e_out, Double_t e_in);
  static TCutG *ScaledCut(TCutG *cut, Double_t scale);
  static Double_t CountInCut(TH2F *scatter, TCutG *cut, Double_t scale);
  static Double_t Enclosed(Double_t nsigma);
  static void PoissonInterval(Double_t count, Double_t &below, Double_t &above);
  // The cut-variation systematic on a strip's cross section [mb]: per
  // condition the larger change of tagged / denominator under its two
  // shifts, from the per-variant counts the cache carries, in quadrature;
  // `detail` lists them. `per_count` converts a count ratio to mb.
  Double_t CutVariation(Int_t reac, Double_t per_count, TString &detail) const;
  // The gas-pressure uncertainty on a point [mb]: relative, 1 / pressure.
  Double_t GasPressureSys(Double_t sigma) const;
  // Return the mean subtracted branching fraction and half its model spread.
  static Double_t SubtractedFraction(const std::vector<TalysCurve> &talys,
                                     const std::vector<TalysCurve> &talys_sub,
                                     Double_t e_out, Double_t e_in,
                                     Double_t &spread);
  // Average a TALYS curve across a strip; covered reports grid coverage.
  static Double_t StripAverage(TGraph *g, Double_t e_out, Double_t e_in,
                               Bool_t &covered);
 
  TFile *cache_ = nullptr;
  Long64_t n_seen_ = 0, n_beam_ = 0;
  Double_t n_gas_ = 0.0, areal_ = 0.0;
  Double_t dE_[18], e_strip0_ = 0.0, cm_frac_ = 0.0, e_mid_[17];
  std::vector<TString> talys_labels_;
  std::vector<std::map<std::pair<Int_t, Int_t>, TGraph *>> talys_raw_;
};
 
#endif