FluxReweight.cxx
Go to the documentation of this file.
1 #include "TH2.h"
2 #include "TH1.h"
3 #include "TF1.h"
4 #include "TLegend.h"
5 #include "TCanvas.h"
6 #include "TObjString.h"
8 #include "CAFAna/Core/Ratio.h"
11 
12 #include <TDirectory.h>
13 
14 #include "CAFAna/Cuts/BeamCuts.h"
18 
19 // Designed for FluxDecomp, can't be used in anything else for now.
20 //---------------------------------------------------------------------------
21 
22 namespace ana
23 {
25  SpectrumLoaderBase& ndloader,
26  const HistAxis& axis,
27  const Cut& fdcut,
28  const SystShifts& shiftMC,
29  const Var& weight,
32  const Cut& ndcut,
33  const FluxDecomp& decomposition,
34  SpectrumLoaderBase& fdloader,
35  const Cut& flavors,
36  SpectrumLoaderBase& extrafdloaderswap,
37  SpectrumLoaderBase& extrafdloadertau
38  )
39  : fRecoNDPi(ndloader, axis, ndcut && flavors && kIsAllFromPions, shiftMC, weight),
40  fRecoNDOther(ndloader, axis, ndcut && flavors && !kIsAllFromPions, shiftMC, weight),
41  fTrueToRecoFDPi(fdloader, axis, fdcut && flavors && kIsAllFromPions, shiftMC, weight),
42  fTrueToRecoFDOther(fdloader, axis, fdcut && flavors && !kIsAllFromPions, shiftMC, weight),
43  // can remove beamnue and antinu cut
44 // fNueCCFromPiPtPz(fdloader, ktpzAxis, ktpTAxis, axis, fdcut && kIsBeamNue && kIsAllFromPions&& kIsAntiNu, shiftMC, weight),
45 // fNueCCFromOther(fdloader, ktpzAxis, ktpTAxis, axis, fdcut && kIsBeamNue && !kIsAllFromPions&& kIsAntiNu, shiftMC, weight),
46  fDecomp(&decomposition),
47  fOwnDecomp(false),
48  fLabel(label),
49  fLatex(latex)
50  {
51 /* extrafdloaderswap.AddReweightableSpectrum(
52  fNueCCFromPiPtPz, axis.GetVar1D(), fdcut && fdflavors && kIsAllFromPions, shiftMC, weight);
53  extrafdloadertau.AddReweightableSpectrum(
54  fNueCCFromPiPtPz, axis.GetVar1D(), fdcut && fdflavors && kIsAllFromPions, shiftMC, weight);
55 
56  extrafdloaderswap.AddReweightableSpectrum(
57  fNueCCFromOther, axis.GetVar1D(), fdcut && fdflavors && !kIsAllFromPions, shiftMC, weight);
58  extrafdloadertau.AddReweightableSpectrum(
59  fNueCCFromOther, axis.GetVar1D(), fdcut && fdflavors && !kIsAllFromPions, shiftMC, weight);
60 */
61  extrafdloaderswap.AddReweightableSpectrum(
62  fTrueToRecoFDPi, axis.GetVar1D(), kTrueE, fdcut && flavors && kIsAllFromPions, shiftMC, weight);
63  extrafdloadertau.AddReweightableSpectrum(
64  fTrueToRecoFDPi, axis.GetVar1D(), kTrueE, fdcut && flavors && kIsAllFromPions, shiftMC, weight);
65  extrafdloaderswap.AddReweightableSpectrum(
66  fTrueToRecoFDOther, axis.GetVar1D(), kTrueE, fdcut && flavors && !kIsAllFromPions, shiftMC, weight);
67  extrafdloadertau.AddReweightableSpectrum(
68  fTrueToRecoFDOther, axis.GetVar1D(), kTrueE, fdcut && flavors && !kIsAllFromPions, shiftMC, weight);
69 
70  }
71 
73  {
74  if(fOwnDecomp) delete fDecomp;
75  }
76 
78  {
81 // ReweightableSpectrum PiPtPz(fNueCCFromPiPtPz);
82 
83  // Pi weight ratio as a TH2D
84  Spectrum pidecompresult(fDecomp->NuePiEstimate());
85  // reweight fd pipz
86 // PiPtPz.ReweightByRecoRatio(pidecompresult);
87  Spectrum otherdecompresult(fDecomp->NueKaEstimate());
88 
89 // TH1* piweight = (TH1*) ((TH3*)PiPtPz.ToTH3(1e20))->ProjectionZ();
90 
91  Ratio pidataMC = FormSmartRatio(
92  pidecompresult, fRecoNDPi,
93  fLabel, "MC ND Pi Reco",
95 
96  resultPi.ReweightToRecoSpectrum( fTrueToRecoFDPi.Unoscillated() * pidataMC );
97 
98  Ratio otherdataMC = FormSmartRatio(
99  otherdecompresult, fRecoNDOther,
100  fLabel, "MC ND Other Reco",
102 
103  resultOther.ReweightToRecoSpectrum( fTrueToRecoFDOther.Unoscillated() * otherdataMC );
104 
105 /*
106  resultPi.ReweightByRecoRatio(piweight);
107  // normalization for now, later use ptpz weight
108  Spectrum kadecompresult(fDecomp->NueKaEstimate());
109  Ratio dataMC = FormSmartRatio(
110  kadecompresult, fRecoND,
111  fLabel, "MC ND Reco",
112  fTrueToRecoFDOther.Unoscillated() );
113 
114  resultOther.ReweightToRecoSpectrum( fTrueToRecoFDOther.Unoscillated() * dataMC );
115 */
116  fReturnFDPi = resultPi;
117  fReturnFDOther = resultOther;
118 
119  return (resultPi + resultOther);
120 
121  }
122 
123  void FluxReweight::SaveTo(TDirectory* dir, const std::string& name) const
124  {
125  TDirectory* tmp = gDirectory;
126 
127  dir = dir->mkdir(name.c_str()); // switch to subdir
128  dir->cd();
129 
130  TObjString("FluxReweight").Write("type");
131  fRecoNDPi.SaveTo(dir, "RecoNDPi");
132  fRecoNDOther.SaveTo(dir, "RecoNDOther");
133  fTrueToRecoFDPi.SaveTo(dir, "TrueToRecoFDPi");
134  fTrueToRecoFDOther.SaveTo(dir, "TrueToRecoFDOther");
135 
136 // fNueCCFromPiPtPz.SaveTo(, "NueCCFromPiPtPz"));
137 // fNueCCFromOther.SaveTo(, "NueCCFromOther"));
138  fDecomp->SaveTo(dir, "Decomp");
139  TObjString(fLabel.c_str()).Write("Label");
140  TObjString(fLatex.c_str()).Write("Latex");
141 
142  dir->Write();
143  delete dir;
144 
145  tmp->cd();
146  }
147 
148  std::unique_ptr<FluxReweight>
150  {
151  dir = dir->GetDirectory(name.c_str()); // switch to subdir
152  assert(dir);
153 
154 // TObjString* dr = (TObjString*)dir->Get("DecompRes");
155 // assert(dr);
156  TObjString* label = (TObjString*)dir->Get("Label");
157  TObjString* latex = (TObjString*)dir->Get("Latex");
158  assert(label);
159  assert(latex);
160 
162  *(Spectrum::LoadFrom(dir, "RecoNDPi")),
163  *(Spectrum::LoadFrom(dir, "RecoNDOther")),
164  *(OscillatableSpectrum::LoadFrom(dir, "TrueToRecoFDPi")),
165  *(OscillatableSpectrum::LoadFrom(dir, "TrueToRecoFDOther")),
166 // *(ReweightableSpectrum::LoadFrom(dir, "NueCCFromPiPtPz")),
167 // *(ReweightableSpectrum::LoadFrom(dir, "NueCCFromOther")),
168  *(ana::LoadFrom<FluxDecomp>(dir, "Decomp").release()),
169  label->GetString().Data(),
170  latex->GetString().Data()
171  );
172 
173  ret->fOwnDecomp = true;
174 // delete dr;
175  delete dir;
176  delete label;
177  delete latex;
178  return std::unique_ptr<FluxReweight>(ret);
179  }
180 
181  void FluxReweight::SavePlots(TDirectory* dir, double potFD) const
182  {
183  TDirectory* tmp = gDirectory;
184  dir->cd();
185 
186  Spectrum fdExtReco( Return().Unoscillated() );
187 
188  Spectrum fdMCRecoPi( fTrueToRecoFDPi.Unoscillated() );
189  Spectrum fdMCRecoOther( fTrueToRecoFDOther.Unoscillated() );
190 
191  Spectrum fdExtRecoPi( fReturnFDPi.Unoscillated() );
192  Spectrum fdExtRecoOther( fReturnFDOther.Unoscillated() );
193 
194  Spectrum decomppiresult( fDecomp->NuePiEstimate());
195  Spectrum decompotherresult( fDecomp->NueKaEstimate());
196 
197  double potND(decomppiresult.POT());
198 
200  fRecoNDPi, decomppiresult, potND,
201  "Data", kBlack, fLatex, "ND Reco Spectrum From Pi", "NDPi");
202 
204  fRecoNDOther, decompotherresult, potND,
205  "Data", kBlack, fLatex, "ND Reco Spectrum From Other", "NDOther");
206 
208  fdMCRecoPi, fdExtRecoPi, potFD,
209  "Extrapolated", kBlue, fLatex, "FD Reco Spectrum From Pi", "FDPi" );
210 
212  fdMCRecoOther, fdExtRecoOther, potFD,
213  "Extrapolated", kBlue, fLatex, "FD Reco Spectrum From Other", "FDOther" );
214 
215  TH1* hreturn = fdExtReco.ToTH1(potFD);
216  hreturn -> Draw("hist");
217  hreturn -> Write("FDReco");
218  tmp->cd();
219 
220  }
221 }
static std::unique_ptr< OscillatableSpectrum > LoadFrom(TDirectory *dir, const std::string &name)
const XML_Char * name
Definition: expat.h:151
tree Draw("slc.nhit")
void latex(double x, double y, std::string txt, double ang=0, int align=22, double size=0.042)
OscillatableSpectrum fTrueToRecoFDPi
Definition: FluxReweight.h:66
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
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
std::basic_string< char > fLatex
Definition: FluxReweight.h:73
const Var weight
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void SavePlots(TDirectory *dir, double potFD) const override
OscillatableSpectrum Eval() const override
Core extrapolation math.
const OscillatableSpectrum & Return() const
Interface so that result of Eval() is called only once and cached.
Float_t tmp
Definition: plot.C:36
void ReweightToRecoSpectrum(const Spectrum &target)
Recale bins so that Unweighted will return target.
void SaveTo(TDirectory *dir, const std::basic_string< char > &name) const override
static std::unique_ptr< FluxReweight > LoadFrom(TDirectory *dir, const std::basic_string< char > &name)
Spectrum fRecoNDOther
Definition: FluxReweight.h:65
Spectrum fRecoNDPi
Definition: FluxReweight.h:64
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
OscillatableSpectrum fReturnFDPi
Definition: FluxReweight.h:74
const Var kTrueE([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].E);})
Definition: Vars.cxx:85
Spectrum NueKaEstimate() const
Definition: FluxDecomp.h:104
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:546
FluxReweight(SpectrumLoaderBase &ndloader, const HistAxis &axis, const Cut &fdcut, const SystShifts &shiftMC, const Var &weight, std::string label, std::string latex, const Cut &ndcut, const FluxDecomp &decomposition, SpectrumLoaderBase &fdloader, const Cut &flavors, SpectrumLoaderBase &extrafdloaderswap=kNullLoader, SpectrumLoaderBase &extrafdloadertau=kNullLoader)
std::vector< std::string > flavors
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:517
OscillatableSpectrum fTrueToRecoFDOther
Definition: FluxReweight.h:67
std::basic_string< char > fLabel
Definition: FluxReweight.h:72
Represent the ratio between two spectra.
Definition: Ratio.h:8
virtual void AddReweightableSpectrum(ReweightableSpectrum &spect, const Var &xvar, const Var &yvar, const Cut &cut, const SystShifts &shift, const Var &wei)
For use by the constructors of ReweightableSpectrum subclasses.
Base class for the various types of spectrum loader.
OscillatableSpectrum fReturnFDOther
Definition: FluxReweight.h:75
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut kIsAllFromPions([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return( (sr->mc.nu[0].beam.ndecay==11|| sr->mc.nu[0].beam.ndecay==12|| sr->mc.nu[0].beam.ndecay==13|| sr->mc.nu[0].beam.ndecay==14)&& ((abs(sr->mc.nu[0].beam.tptype)==211)|| (abs(sr->mc.nu[0].beam.tptype)==3122)|| (abs(sr->mc.nu[0].beam.tptype)==3222)|| (abs(sr->mc.nu[0].beam.tptype)==3112)|| (abs(sr->mc.nu[0].beam.tptype)==3212)) );})
Definition: BeamCuts.h:35
TDirectory * dir
Definition: macro.C:5
const FluxDecomp * fDecomp
Definition: FluxReweight.h:70
assert(nhit_max >=nhit_nbins)
static Ratio FormSmartRatio(const Spectrum &num, const Spectrum &denom, std::string component, std::string location, const Spectrum &mult)
Form Ratio, but be aware of zero division.
Spectrum with true energy information, allowing it to be oscillated
Spectrum NuePiEstimate() const
Definition: FluxDecomp.h:103
Extrapolation for flux true space reweight.
Definition: FluxReweight.h:16
void SaveTo(TDirectory *dir, const std::string &name) const
static void ComparisonPlot(Spectrum mc, Spectrum notMC, double pot, std::string notMCLabel, int notMCColor, std::string latex, std::string title, std::string saveAs, bool restrictRange=false)
T GetVar1D() const
Definition: HistAxis.cxx:43
void SaveTo(TDirectory *dir, const std::string &name) const override
Definition: FluxDecomp.cxx:330
gm Write()