ICrossSectionAnalysis.cxx
Go to the documentation of this file.
3 #include "CAFAna/Core/Spectrum.h"
4 
9 
10 #include "TH2.h"
11 #include "TH3.h"
12 
13 #include <iostream>
14 
15 namespace ana {
16  ///////////////////////////////////////////////////////////////////////
17  void
19  {
20  fPOT = POT;
21  }
22 
23  ///////////////////////////////////////////////////////////////////////
24  void
26  {
27  std::for_each(fXSecSpectra.begin(), fXSecSpectra.end(),
28  [this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*> xsec)
29  {
30  this->fXSec
31  .insert(std::pair<XSecSysts::Syst_t, TH1 *>
32  (xsec.first,
33  CrossSection(xsec.first)));
34  });
35  }
36 
37  ///////////////////////////////////////////////////////////////////////
38  Cut *
40  {
41  if(ptr)
42  return new Cut(*ptr);
43  else
44  return NULL;
45  }
46 
47  ///////////////////////////////////////////////////////////////////////
48  NuTruthCut *
50  {
51  if(ptr)
52  return new NuTruthCut(*ptr);
53  else
54  return NULL;
55  }
56 
57  ///////////////////////////////////////////////////////////////////////
58  HistAxis *
60  {
61  if(ptr)
62  return new HistAxis(*ptr);
63  else
64  return NULL;
65  }
66 
67  ///////////////////////////////////////////////////////////////////////
70  {
71  if(ptr)
72  return new NuTruthHistAxis(*ptr);
73  else
74  return NULL;
75  }
76 
77  ///////////////////////////////////////////////////////////////////////
78  Var *
80  {
81  if(ptr)
82  return new Var(*ptr);
83  else
84  return NULL;
85  }
86 
87  ///////////////////////////////////////////////////////////////////////
88  NuTruthVar *
90  {
91  if(ptr)
92  return new NuTruthVar(*ptr);
93  else
94  return NULL;
95  }
96 
97  ///////////////////////////////////////////////////////////////////////
98  SystShifts *
100  {
101  if(ptr)
102  return new SystShifts(*ptr);
103  else
104  return NULL;
105  }
106 
107  ///////////////////////////////////////////////////////////////////////
108  bool
110  {
111  return
112  HasRecoAxis() &&
113  HasTruthAxis() &&
114  HasSelectionCut() &&
115  HasSignalCut() &&
116  HasFiducialVolume() &&
117  HasWeights() &&
118  HasNuTruthWeights() &&
119  HasSystShifts();
120  }
121 
122  ///////////////////////////////////////////////////////////////////////
123  bool
125  {
126  return
127  HasSpectra() &&
128  HasUnfoldMethod () &&
129  HasUnfoldReg() &&
130  HasDimensions() &&
132  }
133 
134  ///////////////////////////////////////////////////////////////////////
135  const Spectrum *
137  {
138  return fXSecSpectra[syst]->GetSignalEstimator()->Signal(fData);
139  }
140 
141  ///////////////////////////////////////////////////////////////////////
142  void
144  {
145  assert( HasRecoAxis() &&
146  HasSelectionCut() &&
147  "Analysis not ready to fill data spectrum" );
148 
149  fData = new Spectrum(*loader,
150  *fRecoAxis,
151  *fSelectionCut,
152  kNoShift,
153  kUnweighted);
154  }
155 
156  ///////////////////////////////////////////////////////////////////////
157  void
159  {
160  fData = data;
161  }
162  ///////////////////////////////////////////////////////////////////////
163  bool
165  {
166  if(fENuBinning)
167  return true;
168  else {
169  std::cout << "Binning not initialized" << std::endl;
170  return false;
171  }
172  }
173 
174  ///////////////////////////////////////////////////////////////////////
175  bool
177  {
178  if(fPDG != 0)
179  return true;
180  else {
181  std::cout << "Particle PDG not set" << std::endl;
182  return false;
183  }
184  }
185 
186  ///////////////////////////////////////////////////////////////////////
187  bool
189  {
190  if(fData)
191  return true;
192  else {
193  std::cout << "CrossSectionAnalysis object does not have data" << std::endl;
194  return false;
195  }
196  }
197 
198  ///////////////////////////////////////////////////////////////////////
199  bool
201  {
202  if(!fXSecSpectra.empty())
203  return true;
204  else {
205  std::cout << "CrossSectionAnalysis object does not have MC spectra" << std::endl;
206  return false;
207  }
208  }
209 
210  ///////////////////////////////////////////////////////////////////////
211  bool
213  {
214  if(!fWeights.empty())
215  return true;
216  else {
217  std::cout << "Weights not initialized" << std::endl;
218  return false;
219  }
220  }
221 
222  ///////////////////////////////////////////////////////////////////////
223  bool
225  {
226  if(!fNuTruthWeights.empty())
227  return true;
228  else {
229  std::cout << "NuTruthWeights not initialized" << std::endl;
230  return false;
231  }
232  }
233 
234  ///////////////////////////////////////////////////////////////////////
235  bool
237  {
238  if(!fSystShifts.empty())
239  return true;
240  else {
241  std::cout << "SystShifts not initialized" << std::endl;
242  return false;
243  }
244  }
245 
246  ///////////////////////////////////////////////////////////////////////
247  bool
249  {
250  if(!fXSec.empty())
251  return true;
252  else {
253  std::cout << "CrossSectionAnalysis object has not caluclated cross sections" << std::endl;
254  return false;
255  }
256  }
257 
258  ///////////////////////////////////////////////////////////////////////
259  bool
261  {
262  if(fUnfoldMethod != 0)
263  return true;
264  else {
265  std::cout << "Unfolding method not set: " << fUnfoldMethod << std::endl;
266  return false;
267  }
268  }
269 
270  ///////////////////////////////////////////////////////////////////////
271  bool
273  {
274  if(fUnfoldReg != 0)
275  return true;
276  else {
277  std::cout << "Unfolding reg not set" << std::endl;
278  return false;
279  }
280  }
281 
282  ///////////////////////////////////////////////////////////////////////
283  bool
285  {
286  if(fNDims != 0)
287  return true;
288  else {
289  std::cout << "Measurement dimensions not set" << std::endl;
290  return false;
291  }
292  }
293 
294  ///////////////////////////////////////////////////////////////////////
295  bool
297  {
298  if(fSignalCut)
299  return true;
300  else {
301  std::cout << "Signal cut not defined" << std::endl;
302  return false;
303  }
304  }
305 
306  ///////////////////////////////////////////////////////////////////////
307  bool
309  {
310  if(fSelectionCut)
311  return true;
312  else {
313  std::cout << "Selection Cut not defined" << std::endl;
314  return false;
315  }
316  }
317 
318  ///////////////////////////////////////////////////////////////////////
319  bool
321  {
322  if(fRecoAxis)
323  return true;
324  else {
325  std::cout << "Reco axis not initialized" << std::endl;
326  return false;
327  }
328  }
329 
330  ///////////////////////////////////////////////////////////////////////
331  bool
333  {
334  if(fTruthAxis)
335  return true;
336  else {
337  std::cout << "Truth axis not initialized" << std::endl;
338  return false;
339  }
340  }
341 
342  ///////////////////////////////////////////////////////////////////////
343  bool
345  {
346  if((fFidMax.X() != 0 &&
347  fFidMax.Y() != 0 &&
348  fFidMax.Z() != 0) ||
349  (fFidMin.X() != 0 &&
350  fFidMin.Y() != 0 &&
351  fFidMin.Z() != 0))
352  return true;
353  else {
354  std::cout << "Fiducial volume not defined" << std::endl;
355  return false;
356  }
357  }
358 
359  ///////////////////////////////////////////////////////////////////////
360  void
363  {
364  fXSecSpectra.insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *>
365  (syst, xsec));
366  }
367 
368  ///////////////////////////////////////////////////////////////////////
369  void
372  ISignalEstimator * signal_est)
373  {
374 
375  assert( CanFillSpectra() &&
376  HasENuBinning() &&
377  HasPDG() &&
378  "CrossSectionAnalysis object not yet ready to setup systematic spectra." );
379  fXSecSpectra.insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *>
380  (syst,
381  new CrossSectionSpectra(*loader,
382  *fRecoAxis,
383  *fTruthAxis,
384  *fSelectionCut,
385  *fSignalCut,
386  signal_est,
387  fFidMax,
388  fFidMin,
389  *fENuBinning,
390  fPDG,
391  *fWeights[syst],
392  *fNuTruthWeights[syst],
393  *fSystShifts[syst])));
394  }
395 
396  ///////////////////////////////////////////////////////////////////////
397  void
399  ReweightableSpectrum * reco_true,
400  ISignalEstimator * signal_est,
401  Spectrum * flux,
402  Spectrum * mc_selected,
403  Spectrum * sig_selected_truth,
404  Spectrum * sig_selected_reco,
405  Spectrum * sig_truth)
406  {
407  fXSecSpectra.insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *>
408  (syst,
409  new CrossSectionSpectra(reco_true,
410  signal_est,
411  flux,
412  mc_selected,
413  sig_selected_truth,
414  sig_selected_reco,
415  sig_truth)));
416  }
417 
418  ///////////////////////////////////////////////////////////////////////
419  const CrossSectionSpectra*
421  {
422  return fXSecSpectra[syst];
423  }
424 
425  ///////////////////////////////////////////////////////////////////////
426  const TargetCount*
428  {
429  return &fTarget;
430  }
431 
432  ///////////////////////////////////////////////////////////////////////
433  const TVector3*
435  {
436  return &fFidMax;
437  }
438 
439  ///////////////////////////////////////////////////////////////////////
440  const TVector3*
442  {
443  return &fFidMin;
444  }
445 
446  ///////////////////////////////////////////////////////////////////////
447  void
449  const Var * var)
450  {
451  fWeights.insert(std::pair<XSecSysts::Syst_t, const Var *>
452  (syst, new Var(*var)));
453  }
454 
455  ///////////////////////////////////////////////////////////////////////
456  void
458  const NuTruthVar * nt_var)
459  {
460  fNuTruthWeights.insert(std::pair<XSecSysts::Syst_t, const NuTruthVar *>
461  (syst, new NuTruthVar(*nt_var)));
462  }
463 
464  ///////////////////////////////////////////////////////////////////////
465  void
467  const SystShifts * shift)
468  {
469  fSystShifts.insert(std::pair<XSecSysts::Syst_t, const SystShifts *>
470  (syst, new SystShifts(*shift)));
471  }
472 
473  ///////////////////////////////////////////////////////////////////////
474  void
476  {
477  fSelectionCut = cut;
478  }
479 
480  ///////////////////////////////////////////////////////////////////////
481  void
483  {
484  fSignalCut = cut;
485  }
486 
487  ///////////////////////////////////////////////////////////////////////
488  void
490  {
491  fRecoAxis = axis;
492  }
493 
494  ///////////////////////////////////////////////////////////////////////
495  void
497  {
498  fTruthAxis = axis;
499  }
500 
501  ///////////////////////////////////////////////////////////////////////
502  void
504  {
505  fTarget = target;
506  }
507 
508  ///////////////////////////////////////////////////////////////////////
509  void
511  {
512  fFidMax = fidmax;
513  }
514 
515  ///////////////////////////////////////////////////////////////////////
516  void
518  {
519  fFidMin = fidmin;
520  }
521 
522  ///////////////////////////////////////////////////////////////////////
523  TH1 *
525  {
526  return (TH1 *) fXSec[syst]->Clone();
527  }
528 
529  ///////////////////////////////////////////////////////////////////////
530  TH1 *
532  {
534  "CrossSectionAnalysis not ready for cross section calculation" );
535 
536  const TH1 * unfolded_signal_est = UnfoldedSignal(syst);
537  const TH1 * efficiency = Efficiency(syst);
538  TH1 * flux = Flux(syst);
539  flux->Scale(1e-4); // Convert nu/m^2 to nu/cm^2
540  double ntargets = NTargets();
541  TH1 * xsec = (TH1*) unfolded_signal_est->Clone();
542  xsec->Divide(efficiency);
543  xsec->Divide(flux);
544  xsec->Scale(1./ntargets);
545  if(fIsDifferential)
546  DivideByBinWidth(xsec);
547  return xsec;
548  }
549 
550  ///////////////////////////////////////////////////////////////////////
551  const TH1 *
553  {
554  const ReweightableSpectrum * recotrue = fXSecSpectra[syst]->GetRecoTrue();
555  const Spectrum * signal = SignalEstimate(syst);
556 
557  switch(fUnfoldMethod) {
558  case(UnfoldMethod_t::kTikhonov) : {
559  UnfoldTikhonov unfold(*recotrue, fUnfoldReg);
560  return ToHist(unfold.Truth(*signal));
561 
562  }
563  case(UnfoldMethod_t::kSVD) : {
564  UnfoldSVD unfold(*recotrue, fUnfoldReg);
565  return ToHist(unfold.Truth(*signal));
566 
567  }
568  case(UnfoldMethod_t::kMaxEnt) : {
569  UnfoldMaxEnt unfold(*recotrue, fUnfoldReg);
570  return ToHist(unfold.Truth(*signal));
571 
572  }
573  default : {
574  UnfoldIterative unfold(*recotrue, fUnfoldReg);
575  return ToHist(unfold.Truth(*signal));
576  }
577  } // end switch-case
578  }
579 
580  ///////////////////////////////////////////////////////////////////////
581  double
583  {
584  return fTarget.NNucleons();
585  }
586 
587  ///////////////////////////////////////////////////////////////////////
588  int
590  {
591  return fNDims;
592  }
593 
594  ///////////////////////////////////////////////////////////////////////
595  TH1 *
597  double POT)
598  {
599  if(POT < 0)
600  POT = spec.POT();
601 
602  switch(fNDims) {
603  case(1) : {
604  return spec.ToTH1(POT);
605  }
606  case(2) : {
607  return spec.ToTH2(POT);
608  }
609  case(3) : {
610  return spec.ToTH3(POT);
611  }
612  default : {
613  std::cout << "Dimension of measurement undetermined." << std::endl;
614  abort();
615  }
616  } // end switch case
617  }
618 
619  ///////////////////////////////////////////////////////////////////////
620  void
622  {
623  TH1 * binning;
624  switch(fNDims) {
625  case(1) : {
626  binning = (TH1D*) xsec->Clone();
627  DX(binning);
628  break;
629  }
630  case(2) : {
631  binning = (TH2D*) xsec->Clone();
632  DXDY(binning);
633  break;
634  }
635  case(3) : {
636  binning = (TH3D*) xsec->Clone();
637  DXDYDZ(binning);
638  break;
639  }
640  default : {
641  std::cout << "Dimension of measurement undetermined." << std::endl;
642  abort();
643  }
644  } // end of switch-case
645  }
646 
647  ///////////////////////////////////////////////////////////////////////
648  void
650  {
651  for(int i = 1; i <= (int) binning->GetNbinsX(); i++) {
652  binning->SetBinContent(i, binning->GetXaxis()->GetBinWidth(i));
653  }
654  }
655 
656  ///////////////////////////////////////////////////////////////////////
657  void
659  {
660  for(int i = 1; i <= (int) binning->GetNbinsX(); i++) {
661  double dx = binning->GetXaxis()->GetBinWidth(i);
662  for(int j = 1; j <= (int) binning->GetNbinsY(); j++) {
663  double dy = binning->GetYaxis()->GetBinWidth(j);
664  binning->SetBinContent(i, j, dx * dy);
665  }
666  }
667  }
668 
669  ///////////////////////////////////////////////////////////////////////
670  void
672  {
673  for(int i = 1; i <= (int) binning->GetNbinsX(); i++) {
674  double dx = binning->GetXaxis()->GetBinWidth(i);
675  for(int j = 1; j <= (int) binning->GetNbinsY(); j++) {
676  double dy = binning->GetYaxis()->GetBinWidth(j);
677  for(int k = 1; k <= (int) binning->GetNbinsZ(); k++) {
678  double dz = binning->GetZaxis()->GetBinWidth(k);
679  binning->SetBinContent(i, j, k, dx * dy * dz);
680  }
681  }
682  }
683  }
684 
685  ///////////////////////////////////////////////////////////////////////
686  void
688  {
689  this->fUnfoldMethod = method;
690  }
691 
692  ///////////////////////////////////////////////////////////////////////
693  void
695  {
696  this->fUnfoldReg = reg;
697  }
698 
699  ///////////////////////////////////////////////////////////////////////
701  {
702  std::for_each(fXSecSpectra.begin(), fXSecSpectra.end(),
703  [](std::pair<XSecSysts::Syst_t, CrossSectionSpectra *> xsec)
704  {
705  delete xsec.second;
706  });
707  std::for_each(fWeights.begin(), fWeights.end(),
708  [](std::pair<XSecSysts::Syst_t, const Var *> weight)
709  {
710  delete weight.second;
711  });
712  std::for_each(fSystShifts.begin(), fSystShifts.end(),
713  [](std::pair<XSecSysts::Syst_t, const SystShifts *> shift)
714  {
715  delete shift.second;
716  });
717  delete fSelectionCut;
718  delete fSignalCut;
719  delete fRecoAxis;
720  delete fTruthAxis;
721  }
722 }
Maximum entropy unfolding.
Definition: UnfoldMaxEnt.h:8
virtual TH1 * Flux(XSecSysts::Syst_t)=0
http://arxiv.org/abs/hep-ph/9509307
Definition: UnfoldSVD.h:10
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
TH1 * GetCrossSection(XSecSysts::Syst_t)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const XML_Char * target
Definition: expat.h:268
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 Var weight
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Spectrum with the value of a second variable, allowing for reweighting
Loaders::FluxType flux
_Cut< caf::SRNeutrinoProxy > NuTruthCut
Cut designed to be used over the nuTree, ie all neutrinos, not just those that got slices...
Definition: Cut.h:104
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const TH1 * UnfoldedSignal(XSecSysts::Syst_t)
double NNucleons() const
Number of nucleons (mass * avogadro&#39;s number)
Definition: TargetCount.h:31
const XML_Char const XML_Char * data
Definition: expat.h:268
virtual TH1 * Efficiency(XSecSysts::Syst_t)=0
void SetWeight(XSecSysts::Syst_t, const Var *)
double dy[NP][NC]
double dx[NP][NC]
void InitDataSpectrum(SpectrumLoaderBase *)
double dz[NP][NC]
std::map< XSecSysts::Syst_t, CrossSectionSpectra * > fXSecSpectra
virtual Spectrum Truth(const Spectrum &reco) override
Definition: UnfoldSVD.cxx:24
_Var< caf::SRNeutrinoProxy > NuTruthVar
Var designed to be used over the nuTree, ie all neutrinos, not just those that got slices...
Definition: Var.h:82
TH1 * CrossSection(XSecSysts::Syst_t)
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
void SetSystShifts(XSecSysts::Syst_t, const SystShifts *)
std::map< XSecSysts::Syst_t, const NuTruthVar * > fNuTruthWeights
void SetAxisNT(const NuTruthHistAxis *)
const double j
Definition: BetheBloch.cxx:29
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:610
double POT() const
Definition: Spectrum.h:227
void SetNuTruthWeight(XSecSysts::Syst_t, const NuTruthVar *)
Double_t xsec[nknots]
Definition: testXsec.C:47
virtual Spectrum Truth(const Spectrum &reco) override
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
void SetUnfoldingMethod(UnfoldMethod_t)
Base class for the various types of spectrum loader.
void SetSignalCut(const NuTruthCut *)
virtual Spectrum Truth(const Spectrum &reco) override
const Cut cut
Definition: exporter_fd.C:30
std::map< XSecSysts::Syst_t, const Var * > fWeights
TH3 * ToTH3(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 3D to obtain TH3.
Definition: Spectrum.cxx:187
void InitSystematicSpectra(XSecSysts::Syst_t, CrossSectionSpectra *)
std::map< XSecSysts::Syst_t, const SystShifts * > fSystShifts
assert(nhit_max >=nhit_nbins)
UnfoldMethod_t
enum for unfolding methods
const NuTruthHistAxis * fTruthAxis
std::map< XSecSysts::Syst_t, TH1 * > fXSec
TH1 * ToHist(const Spectrum &, double=-1)
const CrossSectionSpectra * GetCrossSectionSpectra(XSecSysts::Syst_t)
_HistAxis< NuTruthVar > NuTruthHistAxis
Definition: HistAxis.h:106
Float_t e
Definition: plot.C:35
void efficiency()
Definition: efficiency.C:58
Count number of targets within the main part of the ND (vSuperBlockND)
Definition: TargetCount.h:10
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
const Spectrum * SignalEstimate(XSecSysts::Syst_t)