CrossSectionAnalysis.cxx
Go to the documentation of this file.
1 #include "CAFAna/Core/Ratio.h"
2 #include "CAFAna/Core/Spectrum.h"
5 
7 
12 
13 
14 #include "TObjString.h"
15 #include "TVectorD.h"
16 
17 #include <cassert>
18 #include <vector>
19 #include <string>
20 #include <iostream>
21 
22 namespace ana
23 {
25  SpectrumLoaderBase& lDA,
26  HistAxis axisreco,
27  NuTruthHistAxis axistrue,
28  Binning EBins,
29  Cut sel,
30  NuTruthCut isSig,
31  IBkgdEstimator* bkgdest,
32  TVector3 min,
33  TVector3 max,
34  double unfoldReg,
35  UnfoldMethod_t unfoldmetod,
36  const SystShifts& shiftMC,
37  const SystShifts& shiftDA,
38  const NuTruthVar& wei,
39  const Var& weiDA)
40 
41  : fData(lDA, axisreco, sel, shiftDA, weiDA),
42  fMC(lMC, axisreco, sel, shiftMC, VarFromNuTruthVar(wei, 1) ),
43  fRecoTrue(lMC,
44  axisreco, HistAxisFromNuTruthHistAxis(axistrue),
45  sel && CutFromNuTruthCut(isSig),
46  shiftMC, VarFromNuTruthVar(wei)),
47  fBkgdEst(bkgdest),
48  fNTargets(0),
49  fUnfoldReg(unfoldReg),
50  fUnfoldMethod(unfoldmetod),
51  fFidMin(min),
52  fFidMax(max),
53  fTarget(min,max,1e6)
54  {
55  std::pair<Spectrum*,Spectrum*> eff = ana::Efficiency(lMC, axistrue,
56  sel, isSig, shiftMC, wei);
57  fMCSig = eff.first;
58  fMCSigNuTree = eff.second;
59  fFlux = DeriveFlux(lMC, EBins, 14, &fFidMin, &fFidMax, shiftMC, wei);
60  }
61 
63  SpectrumLoaderBase& lDA,
64  HistAxis axisreco,
65  NuTruthHistAxis axistrue,
66  Binning EBins,
67  Cut sel,
68  NuTruthCut isSig,
69  IBkgdEstimator* bkgdest,
70  TVector3 min,
71  TVector3 max,
72  double nTargets,
73  double unfoldReg,
74  UnfoldMethod_t unfoldmetod,
75  const SystShifts& shiftMC,
76  const SystShifts& shiftDA,
77  const NuTruthVar& wei,
78  const Var& weiDA)
79 
80  : fData(lDA, axisreco, sel, shiftDA, weiDA),
81  fMC(lMC, axisreco, sel, shiftMC, VarFromNuTruthVar(wei, 1) ),
82  fRecoTrue(lMC,
83  axisreco, HistAxisFromNuTruthHistAxis(axistrue),
84  sel && CutFromNuTruthCut(isSig),
85  shiftMC, VarFromNuTruthVar(wei)),
86  fBkgdEst(bkgdest),
87  fNTargets(nTargets),
88  fUnfoldReg(unfoldReg),
89  fUnfoldMethod(unfoldmetod),
90  fFidMin(min),
91  fFidMax(max)
92  {
93  std::pair<Spectrum*,Spectrum*> eff = ana::Efficiency(lMC, axistrue,
94  sel, isSig, shiftMC, wei);
95  fMCSig = eff.first;
96  fMCSigNuTree = eff.second;
97  fFlux = DeriveFlux(lMC, EBins, 14, &fFidMin, &fFidMax, shiftMC, wei);
98  }
99 
101  {
102  Spectrum bkgd = fBkgdEst->Background();
103  Spectrum sig = SignalEst(bkgd);
104  Spectrum sig_unfolded = UnfoldedSignal(sig);
105  Ratio eff = Efficiency();
106  double ntargs = NTargets();
107  //double flux = fFlux->Integral(fData.POT());
108  TH1* hflux = fFlux->ToTH1(fData.POT());
109 
110  TH1* xsec = (sig_unfolded/eff).ToTH1(fData.POT());
111  hflux->Scale(1e-4); // Convert nu/m^2 to nu/cm^2
112  xsec->Divide(hflux);
113  xsec->Scale(1./ntargs);
114  //xsec->Scale(1./(nnuc*flux));
115 
116  return xsec;
117  }
118 
120  {
121  return fData - fBkgdEst->Background();
122  }
123 
126  {
127 
128  if(rt == NULL)
129  rt = &fRecoTrue;
130 
131  switch(fUnfoldMethod){
132  case(UnfoldMethod_t::kTikhonov) : {
133  UnfoldTikhonov unfold(*rt, fUnfoldReg);
134  return unfold.Truth(sig);
135  }
136  case(UnfoldMethod_t::kSVD) : {
137  UnfoldSVD unfold(*rt, fUnfoldReg);
138  return unfold.Truth(sig);
139  }
140  case(UnfoldMethod_t::kMaxEnt) : {
141  UnfoldMaxEnt unfold(*rt, fUnfoldReg);
142  return unfold.Truth(sig);
143  }
144  default : {
145  UnfoldIterative unfold(*rt, fUnfoldReg);
146  return unfold.Truth(sig);
147  }
148 
149  }// end of switch-case
150 
151  }
152 
154  {
155  return (*fMCSig) / (*fMCSigNuTree);
156  }
157 
158  /// Default implementation returns number of nucleons. Override if needed
160  {
161  if (fNTargets!=0) return fNTargets;
162  else return NucleonCount();
163  }
164 
166  {
167  return fTarget.NNucleons();
168  }
169 
171  {
172  return fTarget.NNuclei(Z, A);
173  }
174 
175  //----------------------------------------------------------------------
176  // Public functions allowing user to print histograms of internal spectra
177  //----------------------------------------------------------------------
178 
180  {
181  return Efficiency().ToTH1();
182  }
183 
184  //----------------------------------------------------------------------
185 
187  {
188  return SignalEst(fBkgdEst->Background()).ToTH1(fData.POT());
189  }
190 
191  //----------------------------------------------------------------------
192 
194  {
195  Spectrum bkgd = fBkgdEst->Background();
196  return UnfoldedSignal(SignalEst(bkgd)).ToTH1(fData.POT());
197  }
198 
199  //----------------------------------------------------------------------
200 
202  {
203  TH1 *h = fFlux->ToTH1(1);
204  //h->Scale(1e-4);
205  h->GetYaxis()->SetTitle("Flux / POT / m^{2}");
206  return h;
207  }
208 
209  //----------------------------------------------------------------------
210 
212  {
213  Spectrum bkgd = fBkgdEst->Background();
214  return bkgd.ToTH1(fData.POT());
215  }
216 
217  //----------------------------------------------------------------------
218 
220  {
221  return fData.ToTH1(fData.POT());
222  }
223 
224  //----------------------------------------------------------------------
225 
227  {
228  return fRecoTrue.ToTH2(fData.POT());
229  }
230 
231  //----------------------------------------------------------------------
232  void CrossSectionAnalysis::SaveTo(TDirectory* dir, const std::string& name) const {
233  TDirectory *tmp = gDirectory;
234 
235  dir = dir->mkdir(name.c_str()); // switch to subdir
236  dir->cd();
237 
238  TObjString("CrossSectionAnlaysis").Write("type");
239 
240  fData.SaveTo(dir, "fData");
241  fMC.SaveTo(dir, "fMC");
242  fMCSig->SaveTo(dir, "fMCSig");
243  fMCSigNuTree->SaveTo(dir, "fMCSigNuTree");
244  fRecoTrue.SaveTo(dir, "fRecoTrue");
245 
246  fBkgdEst->SaveTo(dir, "fBkgdEst");
247 
248  TVectorD unfoldreg(1);
249  unfoldreg[0] = fUnfoldReg;
250  unfoldreg.Write("fUnfoldReg");
251 
252  fFidMin.Write("fFidMin");
253  fFidMax.Write("fFidMax");
254  fFlux->SaveTo(dir, "fFlux");
255 
256  dir->Write();
257  delete dir;
258 
259  tmp->cd();
260  }
261 
262  //----------------------------------------------------------------------
263  std::unique_ptr<CrossSectionAnalysis> CrossSectionAnalysis::LoadFrom(TDirectory* dir, const std::string& name, UnfoldMethod_t unfoldmethod)
264  {
265  dir = dir->GetDirectory(name.c_str()); // switch to subdir
266  assert(dir);
267 
268  Spectrum* data = ana::LoadFrom<Spectrum>(dir, "fData").release();
269  Spectrum* mc = ana::LoadFrom<Spectrum>(dir, "fMC").release();
270  Spectrum* mcsig= ana::LoadFrom<Spectrum>(dir, "fMCSig").release();
271  Spectrum* mcsignutree = ana::LoadFrom<Spectrum>(dir, "fMCSigNuTree").release();
272  ReweightableSpectrum* rt = ana::LoadFrom<ReweightableSpectrum>(dir, "fRecoTrue").release();
273 
274  IBkgdEstimator* bkgdest = ana::LoadFrom<IBkgdEstimator>(dir, "fBkgdEst").release();
275 
276  TVectorD* unfoldreg = (TVectorD*)dir->Get("fUnfoldReg");
277  TVector3* min = (TVector3*)dir->Get("fFidMin");
278  TVector3* max = (TVector3*)dir->Get("fFidMax");
279 
280  Spectrum* flux = ana::LoadFrom<Spectrum>(dir, "fFlux").release();
281 
282  delete dir;
283 
284  return std::unique_ptr<CrossSectionAnalysis>(
285  new CrossSectionAnalysis(*data, *mc, mcsig, mcsignutree, *rt, bkgdest,
286  (*unfoldreg)[0], unfoldmethod, *min, *max, flux));
287  }
288 
289  //------------------------------------------------------------------------------
290  std::unique_ptr<CrossSectionAnalysis> CrossSectionAnalysis::LoadFrom(TDirectory* dir,
291  const std::string& name,
292  UnfoldMethod_t unfoldmethod,
293  double nTargets)
294  {
295  dir = dir->GetDirectory(name.c_str()); // switch to subdir
296  assert(dir);
297 
298  Spectrum* data = ana::LoadFrom<Spectrum>(dir, "fData").release();
299  Spectrum* mc = ana::LoadFrom<Spectrum>(dir, "fMC").release();
300  Spectrum* mcsig= ana::LoadFrom<Spectrum>(dir, "fMCSig").release();
301  Spectrum* mcsignutree = ana::LoadFrom<Spectrum>(dir, "fMCSigNuTree").release();
302  ReweightableSpectrum* rt = ana::LoadFrom<ReweightableSpectrum>(dir, "fRecoTrue").release();
303 
304  IBkgdEstimator* bkgdest = ana::LoadFrom<IBkgdEstimator>(dir, "fBkgdEst").release();
305 
306  TVectorD* unfoldreg = (TVectorD*)dir->Get("fUnfoldReg");
307  TVector3* min = (TVector3*)dir->Get("fFidMin");
308  TVector3* max = (TVector3*)dir->Get("fFidMax");
309 
310  Spectrum* flux = ana::LoadFrom<Spectrum>(dir, "fFlux").release();
311 
312  delete dir;
313 
314  return std::unique_ptr<CrossSectionAnalysis>(
315  new CrossSectionAnalysis(*data, *mc, mcsig, mcsignutree, *rt, bkgdest, nTargets,
316  (*unfoldreg)[0], unfoldmethod, *min, *max, flux));
317  }
318 
319 }
Maximum entropy unfolding.
Definition: UnfoldMaxEnt.h:8
double NucleusCount(int Z=-1, int A=-1)
const XML_Char * name
Definition: expat.h:151
http://arxiv.org/abs/hep-ph/9509307
Definition: UnfoldSVD.h:10
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
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
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
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void SaveTo(TDirectory *dir, const std::string &name) const
Spectrum with the value of a second variable, allowing for reweighting
Float_t tmp
Definition: plot.C:36
virtual Spectrum Background() const =0
Loaders::FluxType flux
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
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
Float_t Z
Definition: plot.C:38
virtual double NTargets()
Default implementation returns number of nucleons. Override if needed.
double NNuclei(int Z=-1, int A=-1) const
Number of nuclei, optionally restricted to a particular element or isotope.
TH1F * hflux
Definition: Xdiff_gwt.C:121
virtual Spectrum Truth(const Spectrum &reco) override
Definition: UnfoldSVD.cxx:24
static std::unique_ptr< CrossSectionAnalysis > LoadFrom(TDirectory *dir, const std::string &name, UnfoldMethod_t unfoldmethod=kIterative)
virtual void SaveTo(TDirectory *dir, const std::string &name) const =0
HistAxis HistAxisFromNuTruthHistAxis(NuTruthHistAxis ntha, double _default)
Definition: HistAxis.cxx:9
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:517
Spectrum SignalEst(Spectrum bkgd)
const double ntargs
double POT() const
Definition: Spectrum.h:219
Double_t xsec[nknots]
Definition: testXsec.C:47
Represent the ratio between two spectra.
Definition: Ratio.h:8
virtual Spectrum Truth(const Spectrum &reco) override
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
Base class for the various types of spectrum loader.
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
UnfoldMethod_t
Enumerator for unfolding methods. TODO: Allow swapping to RooUnfold, TUnfold, etc?
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
virtual Spectrum Truth(const Spectrum &reco) override
static const double A
Definition: Units.h:82
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:7
Spectrum * DeriveFlux(SpectrumLoaderBase &loader, const Binning &bins, int pdg, const TVector3 *min, const TVector3 *max, const SystShifts &shift, const NuTruthVar weight)
Definition: Flux.cxx:56
TDirectory * dir
Definition: macro.C:5
void SaveTo(TDirectory *dir, const std::string &name) const
assert(nhit_max >=nhit_nbins)
std::unique_ptr< IBkgdEstimator > LoadFrom< IBkgdEstimator >(TDirectory *dir, const std::string &label)
ReweightableSpectrum fRecoTrue
CrossSectionAnalysis(SpectrumLoaderBase &lMC, SpectrumLoaderBase &lDA, HistAxis axisreco, NuTruthHistAxis axistrue, Binning EBins, Cut sel, NuTruthCut isSig, IBkgdEstimator *bkgdest, TVector3 min, TVector3 max, double unfoldReg, UnfoldMethod_t unfoldmethod=kIterative, const SystShifts &shiftMC=kNoShift, const SystShifts &shiftDA=kNoShift, const NuTruthVar &wei=kNuTruthUnweighted, const Var &weiDA=kUnweighted)
Template for Var and SpillVar.
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:68
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:7
TH2D * ToTH2(double pot) const