SingleNucAnalysis.cxx
Go to the documentation of this file.
1 #include "CAFAna/Core/Ratio.h"
2 #include "CAFAna/Core/Spectrum.h"
7 
9 
10 #include "TH3.h"
11 #include "TFile.h"
12 #include "TCanvas.h"
13 #include "TVectorD.h"
14 
15 #include <cassert>
16 #include <vector>
17 #include <string> // std::string
18 #include <iostream> // std::cout
19 #include <sstream> // std::stringstream
20 #include <iomanip> // std::setprecision
21 
22 namespace ana
23 {
25  {
26  Spectrum bkgd = fBkgdEst->Background();
27  Spectrum sig = SignalEst(bkgd);
28  Spectrum sig_unfolded = UnfoldedSignal(sig);
29  Ratio eff = Efficiency();
30  double ntargs = NTargets();
31  //double flux = fFlux->Integral(fData.POT());
32  TH1* hflux = fFlux->ToTH1(fData.POT());
33  double flux = hflux->Integral();
34  TH2* xsec = (sig_unfolded/eff).ToTH2(fData.POT());
35  for(int iBinX = 1; iBinX<=xsec->GetXaxis()->GetNbins(); iBinX++)
36  {
37  for(int iBinY = 1; iBinY<=xsec->GetYaxis()->GetNbins(); iBinY++)
38  {
39  double binwidth = xsec->GetXaxis()->GetBinWidth(iBinX)*xsec->GetYaxis()->GetBinWidth(iBinY);
40  double newcontent = xsec->GetBinContent(iBinX,iBinY)/binwidth;
41  xsec->SetBinContent(iBinX,iBinY,newcontent);
42  }
43  }
44  hflux->Scale(1e-4); // Convert nu/m^2 to nu/cm^2
45  //xsec->Divide(hflux);
46  xsec->Scale(1./ntargs);
47  xsec->Scale(1./flux);
48  //xsec->Scale(1./(nnuc*flux));
49 
50  return xsec;
51  }
52 
54  {
55  return Efficiency().ToTH2();
56  }
57 
58  //----------------------------------------------------------------------
60  {
62  }
63 
64  //----------------------------------------------------------------------
65 
67  {
68  Spectrum bkgd = fBkgdEst->Background();
69  return UnfoldedSignal(SignalEst(bkgd)).ToTH2(fData.POT());
70  }
71 
72  //----------------------------------------------------------------------
73 
75  {
76  TH1 *h = fFlux->ToTH1(1);
77  //h->Scale(1e-4);
78  h->GetYaxis()->SetTitle("Flux / POT / m^{2}");
79  return h;
80  }
81 
82  //----------------------------------------------------------------------
83 
85  {
86  Spectrum bkgd = fBkgdEst->Background();
87  return bkgd.ToTH2(fData.POT());
88  }
89 
90  //----------------------------------------------------------------------
91 
93  {
94  return fData.ToTH2(fData.POT());
95  }
96 
97  //----------------------------------------------------------------------
98 
100  {
101  return fRecoTrue.ToTH2(fData.POT());
102  }
103 
104  //----------------------------------------------------------------------
105  std::unique_ptr<SingleNucAnalysis> SingleNucAnalysis::LoadFrom(TDirectory* dir){
106  return LoadFrom(dir, kIterative);
107  }
108 
109  //----------------------------------------------------------------------
110  std::unique_ptr<SingleNucAnalysis> SingleNucAnalysis::LoadFrom(TDirectory* dir,
111  UnfoldMethod_t unfoldmethod)
112  {
113  Spectrum* data = ana::LoadFrom<Spectrum>
114  (dir->GetDirectory("fData")).release();
115  Spectrum* mc = ana::LoadFrom<Spectrum>
116  (dir->GetDirectory("fMC")).release();
117  Spectrum* mcsig= ana::LoadFrom<Spectrum>
118  (dir->GetDirectory("fMCSig")).release();
119  Spectrum* mcsignutree = ana::LoadFrom<Spectrum>
120  (dir->GetDirectory("fMCSigNuTree")).release();
121  ReweightableSpectrum* rt = ana::LoadFrom<ReweightableSpectrum>
122  (dir->GetDirectory("fRecoTrue")).release();
123 
125  (dir->GetDirectory("fBkgdEst")).release();
126 
127  TVectorD* unfoldreg = (TVectorD*)dir->Get("fUnfoldReg");
128  TVector3* min = (TVector3*)dir->Get("fFidMin");
129  TVector3* max = (TVector3*)dir->Get("fFidMax");
130 
131  Spectrum* flux = ana::LoadFrom<Spectrum>
132  (dir->GetDirectory("fFlux")).release();
133 
134  return std::unique_ptr<SingleNucAnalysis>(
135  new SingleNucAnalysis(*data, *mc, mcsig, mcsignutree, *rt, bkgdest,
136  (*unfoldreg)[0], unfoldmethod, *min, *max, flux));
137  }
138 
139 
140 }
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:166
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TH2 * ToTH2() const
Definition: Ratio.cxx:78
Spectrum UnfoldedSignal(Spectrum signal, ReweightableSpectrum *rt=NULL)
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:149
Spectrum with the value of a second variable, allowing for reweighting
virtual Spectrum Background() const =0
static std::unique_ptr< SingleNucAnalysis > LoadFrom(TDirectory *dir)
Loaders::FluxType flux
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const XML_Char const XML_Char * data
Definition: expat.h:268
TH1F * hflux
Definition: Xdiff_gwt.C:121
double NTargets()
Default implementation returns number of nucleons. Override if needed.
Spectrum SignalEst(Spectrum bkgd)
const double ntargs
double POT() const
Definition: Spectrum.h:227
Double_t xsec[nknots]
Definition: testXsec.C:47
Represent the ratio between two spectra.
Definition: Ratio.h:8
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
UnfoldMethod_t
Enumerator for unfolding methods. TODO: Allow swapping to RooUnfold, TUnfold, etc?
TDirectory * dir
Definition: macro.C:5
std::unique_ptr< IBkgdEstimator > LoadFrom< IBkgdEstimator >(TDirectory *dir, const std::string &label)
ReweightableSpectrum fRecoTrue
Float_t e
Definition: plot.C:35
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:66
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
For use with Var2D.
Definition: UtilsExt.cxx:115
TH2D * ToTH2(double pot) const