NumuCCIncAnalysis.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 
18 #include <string> // std::string
19 #include <iostream> // std::cout
20 #include <sstream> // std::stringstream
21 #include <iomanip> // std::setprecision
22 
23 
24 //-----------------------------------------------------------------
25 namespace ana
26 {
28  SpectrumLoader &lDA,
29  TrivialBkgdEstimator* bkgdest,
30  TrivialBkgdEstimator* bkgdest3D,
31  XSecType_t xsectype,
32  double unfoldReg,
33  UnfoldMethod_t unfoldmethod,
34  const SystShifts &shiftMC,
35  const SystShifts &shiftDA,
36  const NuTruthVar wei,
37  const Var weiDA,
38  const bool efficiencyIn3D)
39 
40  : CrossSectionAnalysis(lMC, lDA,
41  xsectype == kVsTCos?
44  xsectype == kVsTCos?
47  enubins,
50  bkgdest,
51  vtxmin, vtxmax,
52  unfoldReg,
53  unfoldmethod,
54  shiftMC, shiftDA, wei, weiDA),
55  fData3D(lDA, kRecoMuKEVsCosVsEnuStandardAxis, kAllNumuCCCuts, shiftDA, weiDA),
59  shiftMC, VarFromNuTruthVar(wei))
60  {
61 
62  fUnfoldMethod = unfoldmethod;
63  fXSecType = xsectype;
64 
65  // for the 3D efficiency, use the 3D axis and extract the true and true-selected signal
66  // from it
67  std::pair<Spectrum*,Spectrum*> eff = ana::Efficiency(lMC,
70  kIsTrueSigST, shiftMC, wei);
71  fMCSig3D = eff.first;
72  fMCSigNuTree3D = eff.second;
73  fEfficiencyIn3D = efficiencyIn3D;
74  fBkgdEst3D = bkgdest3D;
75 
76  }// End of constructor
77 
78  //-----------------------------------------------------------------
79 
81  TrivialBkgdEstimator* bkgdest,
82  TrivialBkgdEstimator* bkgdest3D,
83  XSecType_t xsectype,
84  double unfoldReg,
85  UnfoldMethod_t unfoldmethod,
86  const SystShifts &shiftMC,
87  const SystShifts &shiftDA,
88  const NuTruthVar wei,
89  const Var weiDA)
90 
91  : CrossSectionAnalysis(lMC, lMC,
92  xsectype == kVsTCos?
95  xsectype == kVsTCos?
98  enubins,
99  kIsTrueSig,
100  kIsTrueSigST,
101  bkgdest,
102  vtxmin, vtxmax,
103  unfoldReg,
104  unfoldmethod,
105  shiftMC, shiftDA, wei, weiDA),
108  fRecoTrue3D(lMC,
111  shiftMC, VarFromNuTruthVar(wei))
112  {
113 
114  fUnfoldMethod = unfoldmethod;
115  fXSecType = xsectype;
116 
117  std::pair<Spectrum*,Spectrum*> eff = ana::Efficiency(lMC, kTrueMuKEVsCosVsEnuStandardAxisST,
118  kAllNumuCCCuts, kIsTrueSigST, shiftMC, wei);
119 
120  fMCSig3D = eff.first;
121  fMCSigNuTree3D = eff.second;
122  fBkgdEst3D = bkgdest3D;
123 
124  fEfficiencyIn3D = false;
125  }// End of constructor
126 
127  //-----------------------------------------------------------------
128 
130  {
131  if(!fEfficiencyIn3D)
132  return (*fMCSig) / (*fMCSigNuTree);
133  else
134  return (*fMCSig3D) / (*fMCSigNuTree3D);
135  }// end of Efficiency
136 
137 
138  //-----------------------------------------------------------------
139 
141  {
142  if(!fEfficiencyIn3D)
143  return fData - fBkgdEst->Background();
144 
145  else
146  return fData3D - fBkgdEst3D->Background();
147  }
148 
149  //-----------------------------------------------------------------
150 
152  {
153  if(!fEfficiencyIn3D)
155  else
157  }
158 
159  //-----------------------------------------------------------------
160 
161 
163  {
164  if(fXSecType == kVsE){
165  std::cerr<<"NumuCCIncAnalysis class has been setup for 1D cross-section vs Energy\n"
166  <<"Can not retrieve 2D results. Returning NULL histogram\n";
167  return NULL;
168  }
169 
170  Ratio eff = Efficiency();
171  Spectrum sig_unfolded = UnfoldedSignal(SignalEst());
172 
173  Spectrum xsec_spec = (sig_unfolded/eff);
174 
175  TH1* hflux = fFlux->ToTH1(fData.POT());
176  hflux->Scale(1e-4); // Convert nu/m^2 to nu/cm^2
177  double flux = hflux->Integral();
178 
179  xsec_spec.Scale(1./flux);
180  xsec_spec.Scale(1./fNNucleons);
181 
182  TH2 *xsec = ToTH2(xsec_spec, fData.POT(), kPOT,
184  return xsec;
185 
186  }// End of Result
187 
188  //-----------------------------------------------------------------
189 
190  std::vector<TH1*> NumuCCIncAnalysis::Result()
191  {
192 
193  Ratio eff = Efficiency();
194  Spectrum sig_unfolded = UnfoldedSignal(SignalEst());
195 
196  Spectrum xsec_spec = (sig_unfolded/eff);
197 
198  TH1* hflux = fFlux->ToTH1(fData.POT());
199  hflux->Scale(1e-4); // Convert nu/m^2 to nu/cm^2
200  double flux = hflux->Integral();
201 
202  if(fXSecType == kVsE){
203  TH1* xsec;
204  if(!fEfficiencyIn3D)
205  xsec = xsec_spec.ToTH1(fData.POT(), kPOT);
206  else
207  xsec = Project(xsec_spec, fData.POT());
208 
209  xsec->Divide(hflux);
210  xsec->Scale(1./fNNucleons);
211  xsec->GetYaxis()->SetTitle("#sigma (cm^{2}/nucleon)");
212  xsec->GetXaxis()->SetTitle("Neutrino Energy, E_{#nu} (GeV)");
213  CenterTitles(xsec);
214  return {xsec};
215  }
216  else{
217  xsec_spec.Scale(1./flux);
218  xsec_spec.Scale(1./fNNucleons);
219  return ReturnHists(xsec_spec);
220  }
221  }// End of Result
222 
223 
224  //-----------------------------------------------------------------
225 
227  {
228 
229  if(fXSecType != kVsE){
230  TH2* temp = fEfficiencyIn3D ? (TH2*)Project(Efficiency()):
232  return PlotHelper( temp);
233  }
234  else{
235  if(fEfficiencyIn3D)
236  return {Project(Efficiency())};
237  return {Efficiency().ToTH1()};
238  }
239  }
240 
241  //-----------------------------------------------------------------
242 
244  {
245  return ReturnHists(SignalEst());
246  }
247 
248  //-----------------------------------------------------------------
249 
251  {
253  }
254 
255  //-----------------------------------------------------------------
256 
258  {
260  }
261 
262  //-----------------------------------------------------------------
263 
265  {
267  }
268 
269  //-----------------------------------------------------------------
270 
272  {
274  }
275 
276  //-----------------------------------------------------------------
277 
279  {
280  Spectrum bkgd = fBkgdEst->Background();
281 
282  if(fXSecType != kVsE)
283  return PlotHelper( ToTH2( bkgd, fData.POT(), kPOT,
285  else
286  return {bkgd.ToTH1(fData.POT())};
287  }
288 
289  //-----------------------------------------------------------------
290 
291  std::vector<TH1*> NumuCCIncAnalysis::PlotData()
292  {
293  if(fXSecType != kVsE)
294  return PlotHelper( ToTH2( fData, fData.POT(), kPOT,
296  else
297  return {fData.ToTH1(fData.POT())};
298  }
299 
300  //-----------------------------------------------------------------
301 
302  std::vector<TH1*> NumuCCIncAnalysis::PlotMCReco()
303  {
305  ReturnHists(fMC);
306  }
307 
308  //-----------------------------------------------------------------
309 
310  std::vector<TH1*> NumuCCIncAnalysis::ReturnHists( Spectrum spec)
311  {
312 
313  if(fXSecType != kVsE){
314  TH2* temp = fEfficiencyIn3D ? (TH2*)Project(spec, fData.POT()):
315  ToTH2( spec, fData.POT(), kPOT,
317  return PlotHelper( temp);
318  }
319  else{
320  if(fEfficiencyIn3D)
321  return {Project(spec, fData.POT())};
322  return {spec.ToTH1(fData.POT())};
323  }
324  }
325 
326  //-----------------------------------------------------------------
327 
329  {
330 
331  if(fXSecType == kVsE){
332  TH1* h = temp->Project3D("z");
333  return h;
334  }
335  else {
336  TH2* h = (TH2*)temp->Project3D("yx");
337  h->Scale(1, "width"); // normalize by bin-width
338  return h;
339  }
340  }
341 
342  //-----------------------------------------------------------------
343 
345  const double pot)
346  {
347  TH3* temp = ana::ToTH3(spec, pot, kPOT, angbinsCustom,
348  mukebins, enubins);
349 
350  return ProjectHist(temp);
351  }
352 
353  //-----------------------------------------------------------------
354 
355 
357  {
358  TH3* temp = ana::ToTH3(rat, angbinsCustom,
359  mukebins, enubins);
360 
361  return ProjectHist(temp);
362  }
363 
364  //-----------------------------------------------------------------
365 
366 
367  std::vector<TH1*> NumuCCIncAnalysis::PlotHelper(const TH2* spec)
368  {
369  std::vector<TH1*> ret;
370 
371  int nbins = angbinsCustom.NBins();
372  for(int ibin = 1; ibin < nbins+1; ibin++){
373 
374  std::stringstream low;
375  low<<std::fixed<<std::setprecision(3)<<
376  spec->GetXaxis()->GetBinLowEdge(ibin);
377  std::stringstream high;
378  high<<std::fixed<<std::setprecision(3)<<
379  spec->GetXaxis()->GetBinUpEdge(ibin);
380 
381  std::string caption = low.str() + " #leq cos #theta_{#mu} < " + high.str();
382 
383  TH1* h = spec->ProjectionY(Form("bin%d",ibin),ibin,ibin);
384  h->SetTitle(caption.c_str());
385  h->GetXaxis()->SetTitle("Muon Kinetic Energy, T_{#mu} (GeV)");
386  h->GetYaxis()->SetTitle("#frac{d^{2} #sigma}{dcos#theta dT_{#mu}} #left(#frac{10^{-38} cm^{2}}{nucleon GeV}#right)");
387  h->Scale(1e+38);
388  h->SetName(UniqueName().c_str());
389  CenterTitles(h);
390  ret.push_back(h);
391  }// end loop over cos theta bins
392 
393  return ret;
394  }// End of PlotHelper
395 
396  //-----------------------------------------------------------------
397 
398  void NumuCCIncAnalysis::SaveTo(TDirectory* dir) const {
399 
401 
402  dir->cd();
403  fData3D.SaveTo(dir->mkdir("fData3D"));
404  fMC3D.SaveTo(dir->mkdir("fMC3D"));
405  fMCSig3D->SaveTo(dir->mkdir("fMCSig3D"));
406  fMCSigNuTree3D->SaveTo(dir->mkdir("fMCSigNuTree3D"));
407  fBkgdEst3D->SaveTo(dir->mkdir("fBkgdEst3D"));
408  fRecoTrue3D.SaveTo(dir->mkdir("fRecoTrue3D"));
409  }
410 
411  //-----------------------------------------------------------------
412 
413  std::unique_ptr<NumuCCIncAnalysis> NumuCCIncAnalysis::LoadFrom(TDirectory* dir,
414  UnfoldMethod_t unfoldmethod,
415  bool eff3D, bool muonkine,
416  double regStrength)
417  {
418  Spectrum* data = ana::LoadFrom<Spectrum>
419  (dir->GetDirectory("fData")).release();
420  Spectrum* mc = ana::LoadFrom<Spectrum>
421  (dir->GetDirectory("fMC")).release();
422  Spectrum* mcsig= ana::LoadFrom<Spectrum>
423  (dir->GetDirectory("fMCSig")).release();
424  Spectrum* mcsignutree = ana::LoadFrom<Spectrum>
425  (dir->GetDirectory("fMCSigNuTree")).release();
426 
427  Spectrum* data3D = ana::LoadFrom<Spectrum>
428  (dir->GetDirectory("fData3D")).release();
429  Spectrum* mc3D = ana::LoadFrom<Spectrum>
430  (dir->GetDirectory("fMC3D")).release();
431  Spectrum* mcsig3D= ana::LoadFrom<Spectrum>
432  (dir->GetDirectory("fMCSig3D")).release();
433  Spectrum* mcsignutree3D = ana::LoadFrom<Spectrum>
434  (dir->GetDirectory("fMCSigNuTree3D")).release();
435 
436 
437  ReweightableSpectrum* rt = ana::LoadFrom<ReweightableSpectrum>
438  (dir->GetDirectory("fRecoTrue")).release();
439  ReweightableSpectrum* rt3D = ana::LoadFrom<ReweightableSpectrum>
440  (dir->GetDirectory("fRecoTrue3D")).release();
441 
442 
444  (dir->GetDirectory("fBkgdEst")).release();
446  (dir->GetDirectory("fBkgdEst3D")).release();
447 
448 
449  double unfoldreg = regStrength;
450  if(regStrength == -5 ){
451  TVectorD* temp = (TVectorD*)dir->Get("fUnfoldReg");
452  unfoldreg = (*temp)[0];
453  }
454 
455  Spectrum* flux = ana::LoadFrom<Spectrum>
456  (dir->GetDirectory("fFlux")).release();
457 
458 
460  if(data->GetBinnings()[0].NBins() == enubins.NBins() && !(eff3D && muonkine))
461  xsectype = ana::NumuCCIncAnalysis::kVsE;
462 
463  return std::unique_ptr<NumuCCIncAnalysis>(new NumuCCIncAnalysis(*data, *data3D, *mc, *mc3D,
464  mcsig, mcsig3D, mcsignutree, mcsignutree3D,
465  *rt, *rt3D,
466  bkgdest, bkgdest3D,
467  unfoldreg, unfoldmethod, xsectype, flux,
468  eff3D));
469 
470  }
471 
472 
473 }
const NuTruthHistAxis kTrueMuKEVsCosVsEnuStandardAxisST("True T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)", angvsmukevsebins, kTrueMuKEVsCosVsEnuST)
const std::vector< Binning > & GetBinnings() const
Definition: Spectrum.h:264
Ratio Efficiency()
Return the efficiency Ratio.
const HistAxis kTrueEStandardAxis
const NuTruthCut kIsTrueSigST
const NuTruthHistAxis kTrueMuKEVsCosStandardAxisST("True T_{#mu} vs cos #{theta};T_{#mu} [GeV];cos #{theta}", angvsmukebins, kTrueMuKEVsCosST)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Spectrum UnfoldedSignal(Spectrum signal, ReweightableSpectrum *rt=NULL)
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:67
XSecType_t fXSecType
Which xsec are we measureing? vs Energy or vs Muon Kinematics?
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:148
const TVector3 vtxmin(-130,-176, 225)
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void SaveTo(TDirectory *dir, const std::string &name) const
IBkgdEstimator * fBkgdEst3D
Spectrum with the value of a second variable, allowing for reweighting
Divide bin contents by bin widths.
Definition: UtilsExt.h:31
NumuCCIncAnalysis(SpectrumLoader &lMC, SpectrumLoader &lDA, TrivialBkgdEstimator *bkgdest, TrivialBkgdEstimator *bkgdest3D, XSecType_t xsectype=kVsE, double unfoldReg=1, UnfoldMethod_t unfoldmethod=kIterative, const SystShifts &shiftMC=kNoShift, const SystShifts &shiftDA=kNoShift, const NuTruthVar wei=kXSecCVWgt2017ST *kPPFXFluxCVWgtST, const Var weiDA=kXSecCVWgt2017 *kPPFXFluxCVWgt, const bool efficiencyIn3D=false)
Constructor to create a NumuCCIncAnalysis object. lMC is the MC loader lDA is the data loader bkgdest...
OStream cerr
Definition: OStream.cxx:7
std::vector< TH1 * > PlotHelper(const TH2 *spec)
Helper function to convert a 2D plot of the muon cos theta vs muon kinetic energy into a vector of 1D...
ReweightableSpectrum fRecoTrue3D
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
virtual Spectrum Background() const =0
const HistAxis kTrueMuKEVsCosVsEnuStandardAxis
TH1 * ProjectHist(const TH3 *temp)
Helper function to project from 3D to 2D or 1D
const HistAxis kRecoMuKEVsCosStandardAxis("Reconstructed T_{#mu} vs cos #{theta};T_{#mu} [GeV];cos #{theta}", angvsmukebins, kRecoMuKEVsCos)
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:237
Loaders::FluxType flux
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
std::vector< TH1 * > ReturnHists(Spectrum spec)
Makes the vector of 1D hists by calling Project and PlotHelper.
std::vector< TH1 * > PlotEfficiencyCorrectedSignal()
Returns a vector of histograms with the best estimate of signal, after background subtraction...
const XML_Char const XML_Char * data
Definition: expat.h:268
std::vector< TH1 * > PlotSignalEstimate()
Returns a vector of histograms containing signal estimate that is, (data - background) in reco space...
std::vector< TH1 * > PlotTrueSelectedSignal()
Returns a vector of histograms of MC in true space, with true and reco signal selection cuts applied...
std::vector< TH1 * > Result()
Returns a vector of histograms of the cross-section results.
const int nbins
Definition: cellShifts.C:15
const Binning enubins
std::vector< TH1 * > PlotUnfoldedSignal()
Returns a vector of histograms containing unfolded signal.
TH1 * Project(const Spectrum spec, const double pot)
TH1F * hflux
Definition: Xdiff_gwt.C:121
std::vector< TH1 * > PlotData()
Returns a vector of histograms of data in reco space.
#define pot
virtual void SaveTo(TDirectory *dir, const std::string &name) const =0
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:506
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
std::pair< Spectrum *, Spectrum * > Efficiency(SpectrumLoaderBase &loader, const NuTruthHistAxis &axis, const Cut &seln, const NuTruthCut &truthSeln, const SystShifts &shift, const NuTruthVar &weight)
Definition: Efficiency.cxx:10
const Cut kIsTrueSig
UnfoldMethod_t
Enumerator for unfolding methods. TODO: Allow swapping to RooUnfold, TUnfold, etc?
bool fEfficiencyIn3D
Are we measuring the cross-section in 3D?
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:7
const Binning mukebins
Definition: NumuCCIncBins.h:90
static std::unique_ptr< NumuCCIncAnalysis > LoadFrom(TDirectory *dir, UnfoldMethod_t unfoldmethod, bool eff3D=false, bool muonkine=false, double regStrength=-5)
Load the spectra etc from the file and create a NumuCCAnalysis object.
std::vector< TH1 * > PlotMCReco()
Returns a vector of histograms of MC in the reco space For comparison with the output of ana::NumuCCI...
const HistAxis kRecoMuKEVsCosVsEnuStandardAxis("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevsebins, kRecoMuKEVsCosVsEnu)
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
TDirectory * dir
Definition: macro.C:5
XSecType_t
Enumeration of the type of cross-section we are measuring.
int NBins() const
Definition: Binning.h:29
void SaveTo(TDirectory *dir, const std::string &name) const
TH2 * Result2D()
Represent the double-differential result in 2D, rather than in 1D slices.
std::unique_ptr< IBkgdEstimator > LoadFrom< IBkgdEstimator >(TDirectory *dir, const std::string &label)
std::vector< TH1 * > PlotEfficiency()
Generic organizational class for a cross section analysis.
const HistAxis kTrueMuKEVsCosStandardAxis
const TVector3 vtxmax(160, 160, 1000)
Template for Var and SpillVar.
Float_t e
Definition: plot.C:35
const Binning angbinsCustom
Definition: NumuCCIncBins.h:72
void SaveTo(TDirectory *dir) const
Save the NumuCCAnalysis cross-section object to file so it can be used later.
const HistAxis kRecoEStandardAxis("Reconstructed Neutrino Energy (GeV)", enubins, kRecoE)
const NuTruthHistAxis kTrueEStandardAxisST("True Neutrino Energy (GeV)", enubins, kTrueEST)
std::vector< TH1 * > PlotTrueSignal()
Returns a vector of histograms with MC signal in true space drawn from nuTree. For comparison with th...
std::vector< TH1 * > PlotBackgroundEstimate()
Returns a vector of histograms with background estimate.
Spectrum SignalEst()
Return estimated signal, ie data-bg.
TH3 * ToTH3(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Same as ToTH2, but with 3 dimensions.
Definition: UtilsExt.cxx:162
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
Just return the MC prediction for the background.
const Cut kAllNumuCCCuts
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
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:7
Spectrum UnfoldedSignal(Spectrum sig)
Function to return the unfolded spectrum. It checks if the cross-section is being computed in 3D or 1...
const double fNNucleons
Pre-computed number of nucleon targets in our fiducial volume.
enum BeamMode string