NumuCC2p2hAnalysis.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  hflux->Scale(1e-4); // Convert nu/m^2 to nu/cm^2
34  double flux = hflux->Integral();
35  TH2* xsec = (sig_unfolded/eff).ToTH2(fData.POT());
36  for(int iBinX = 1; iBinX<=xsec->GetXaxis()->GetNbins(); iBinX++)
37  {
38  for(int iBinY = 1; iBinY<=xsec->GetYaxis()->GetNbins(); iBinY++)
39  {
40  double binwidth = xsec->GetXaxis()->GetBinWidth(iBinX)*xsec->GetYaxis()->GetBinWidth(iBinY);
41  double newcontent = xsec->GetBinContent(iBinX,iBinY)/binwidth;
42  xsec->SetBinContent(iBinX,iBinY,newcontent);
43  }
44  }
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  Spectrum bkgd = fBkgdEst->Background();
56  Spectrum sig = SignalEst(bkgd);
57  Spectrum sig_unfolded = UnfoldedSignal(sig);
58  Ratio eff = Efficiency();
59  double ntargs = NTargets();
60  //double flux = fFlux->Integral(fData.POT());
61  TH1* hflux = fFlux->ToTH1(fData.POT());
62  hflux->Scale(1e-4); // Convert nu/m^2 to nu/cm^2
63  double flux = hflux->Integral();
64  TH1* xsec = (sig_unfolded/eff).ToTH2(fData.POT());
65  for(int iBinX = 1; iBinX<=xsec->GetXaxis()->GetNbins(); iBinX++)
66  {
67  double binwidth = xsec->GetXaxis()->GetBinWidth(iBinX);
68  double newcontent = xsec->GetBinContent(iBinX)/binwidth;
69  xsec->SetBinContent(iBinX,newcontent);
70  }
71 
72  //xsec->Divide(hflux);
73  xsec->Scale(1./ntargs);
74  xsec->Scale(1./flux);
75  //xsec->Scale(1./(nnuc*flux));
76 
77  return xsec;
78  }
79 
80 
82  {
83  return Efficiency().ToTH2();
84  }
85 
86  //----------------------------------------------------------------------
88  {
90  }
91 
92  //----------------------------------------------------------------------
93 
95  {
96  Spectrum bkgd = fBkgdEst->Background();
97  return UnfoldedSignal(SignalEst(bkgd)).ToTH2(fData.POT());
98  }
99 
100  //----------------------------------------------------------------------
101 
103  {
104  TH1 *h = fFlux->ToTH1(1);
105  //h->Scale(1e-4);
106  h->GetYaxis()->SetTitle("Flux / POT / m^{2}");
107  return h;
108  }
109 
110  //----------------------------------------------------------------------
111 
113  {
114  Spectrum bkgd = fBkgdEst->Background();
115  return bkgd.ToTH2(fData.POT());
116  }
117 
118  //----------------------------------------------------------------------
119 
121  {
122  return fData.ToTH2(fData.POT());
123  }
124 
125  //----------------------------------------------------------------------
126 
128  {
129  return fRecoTrue.ToTH2(fData.POT());
130  }
131  //----------------------------------------------------------------------
132  std::unique_ptr<NumuCC2p2hAnalysis> NumuCC2p2hAnalysis::LoadFrom(TDirectory* dir){
133  return LoadFrom(dir, kIterative);
134  }
135 
136  //----------------------------------------------------------------------
137  std::unique_ptr<NumuCC2p2hAnalysis> NumuCC2p2hAnalysis::LoadFrom(TDirectory* dir,
138  UnfoldMethod_t unfoldmethod)
139  {
140  Spectrum* data = ana::LoadFrom<Spectrum>
141  (dir->GetDirectory("fData")).release();
142  Spectrum* mc = ana::LoadFrom<Spectrum>
143  (dir->GetDirectory("fMC")).release();
144  Spectrum* mcsig= ana::LoadFrom<Spectrum>
145  (dir->GetDirectory("fMCSig")).release();
146  Spectrum* mcsignutree = ana::LoadFrom<Spectrum>
147  (dir->GetDirectory("fMCSigNuTree")).release();
148  ReweightableSpectrum* rt = ana::LoadFrom<ReweightableSpectrum>
149  (dir->GetDirectory("fRecoTrue")).release();
150 
152  (dir->GetDirectory("fBkgdEst")).release();
153 
154  TVectorD* unfoldreg = (TVectorD*)dir->Get("fUnfoldReg");
155  TVector3* min = (TVector3*)dir->Get("fFidMin");
156  TVector3* max = (TVector3*)dir->Get("fFidMax");
157 
158  Spectrum* flux = ana::LoadFrom<Spectrum>
159  (dir->GetDirectory("fFlux")).release();
160 
161  return std::unique_ptr<NumuCC2p2hAnalysis>(new NumuCC2p2hAnalysis(*data, *mc, mcsig, mcsignutree, *rt, bkgdest,
162  (*unfoldreg)[0], unfoldmethod, *min, *max, flux));
163  }
164  //-----------------------------------------------------------------------
165  std::unique_ptr<NumuCC2p2hAnalysis> NumuCC2p2hAnalysis::LoadFrom(TDirectory* dir,
166  UnfoldMethod_t unfoldmethod,
167  double nTargets)
168  {
169  Spectrum* data = ana::LoadFrom<Spectrum>
170  (dir->GetDirectory("fData")).release();
171  Spectrum* mc = ana::LoadFrom<Spectrum>
172  (dir->GetDirectory("fMC")).release();
173  Spectrum* mcsig= ana::LoadFrom<Spectrum>
174  (dir->GetDirectory("fMCSig")).release();
175  Spectrum* mcsignutree = ana::LoadFrom<Spectrum>
176  (dir->GetDirectory("fMCSigNuTree")).release();
177  ReweightableSpectrum* rt = ana::LoadFrom<ReweightableSpectrum>
178  (dir->GetDirectory("fRecoTrue")).release();
179 
181  (dir->GetDirectory("fBkgdEst")).release();
182 
183  TVectorD* unfoldreg = (TVectorD*)dir->Get("fUnfoldReg");
184  TVector3* min = (TVector3*)dir->Get("fFidMin");
185  TVector3* max = (TVector3*)dir->Get("fFidMax");
186 
187  Spectrum* flux = ana::LoadFrom<Spectrum>
188  (dir->GetDirectory("fFlux")).release();
189 
190  return std::unique_ptr<NumuCC2p2hAnalysis>(new NumuCC2p2hAnalysis(*data, *mc, mcsig, mcsignutree, *rt, bkgdest, nTargets,
191  (*unfoldreg)[0], unfoldmethod, *min, *max, flux));
192  }
193 }
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
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
virtual double NTargets()
Default implementation returns number of nucleons. Override if needed.
TH1F * hflux
Definition: Xdiff_gwt.C:121
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)
static std::unique_ptr< NumuCC2p2hAnalysis > LoadFrom(TDirectory *dir)
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