ExtrapSterile.cxx
Go to the documentation of this file.
2 
3 #include <utility>
4 
7 #include "CAFAna/Core/Loaders.h"
9 
11 
12 #include "TDirectory.h"
13 #include "TH1D.h"
14 #include "TObject.h"
15 #include "TObjString.h"
16 
17 namespace ana
18 {
19  REGISTER_LOADFROM("ModularExtrapSterile", IExtrap, ModularExtrapSterile);
20 
21  //---------------------------------------------------------------------------
24  const IDecomp& NCSurvDecomp,
25  const IDecomp& NumuOscDecomp,
26  const HistAxis& axis,
27  const HistAxis& axisNumuND,
28  const Cut& fdcut,
29  const Cut& NCNDcut,
30  const Cut& NumuNDcut,
31  const SystShifts& shiftMC,
32  const Var& weight
33  ) {
35  loaders.GetLoader(
37  loaders.GetLoader(
39  loaders.GetLoader(
41  loaders.GetLoader(
43  NCSurvDecomp,
44  NumuOscDecomp,
45  axis,
46  axisNumuND,
47  fdcut,
48  NCNDcut,
49  NumuNDcut,
50  shiftMC,
51  weight
52  );
53  }
54 
55  //---------------------------------------------------------------------------
57  SpectrumLoaderBase& nearMC,
58  SpectrumLoaderBase& farMCswap,
59  SpectrumLoaderBase& farMCnonswap,
60  SpectrumLoaderBase& farMCtauswap,
61  const IDecomp& NCSurvDecomp,
62  const IDecomp& NumuOscDecomp,
63  const HistAxis& axis,
64  const HistAxis& axisNumuND,
65  const Cut& fdcut,
66  const Cut& NCNDcut,
67  const Cut& NumuNDcut,
68  const SystShifts& shiftMC,
69  const Var& weight
70  ) {
71  // Initialize all components with NoReweights,
72  // make sure other objects exist
73  ModularExtrapSterile extrap(
74  farMCswap,
75  farMCnonswap,
76  farMCtauswap,
77  axis,
78  fdcut,
79  shiftMC,
80  weight
81  );
82 
83  // NC -> NC ----
84  extrap.fNCextrap = std::unique_ptr<ModularExtrapComponent>(
85  new RecoReweight(
86  nearMC, axis, fdcut, shiftMC, weight,
87  "NC -> NC", "NC #rightarrow NC",
88  NCNDcut, NCSurvDecomp, // NC selection in ND
89  DecompResult::NCtot, kIsNC, // NC truth in ND
90  farMCswap, kIsNC, // NC->NC in FD
91  farMCnonswap, farMCtauswap // extra NC stats
92  )
93  );
94 
95  // mu -> mu ----
96  extrap.fMMextrap = std::unique_ptr<ModularExtrapComponent>(
97  new RecoReweight(
98  nearMC, axis, fdcut, shiftMC, weight,
99  "mu -> mu", "#mu #rightarrow #mu",
100  NCNDcut, NCSurvDecomp, // NC selection in ND
101  DecompResult::numu, kIsNumuCC && !kIsAntiNu, // numu truth in ND
102  farMCnonswap, kIsNumuCC && !kIsAntiNu // mu->mu in FD
103  )
104  );
105 
106  // mu -> e ----
107  extrap.fMEextrap = std::unique_ptr<ModularExtrapComponent>(
108  new TruthReweight(
109  nearMC, axis, axisNumuND, fdcut, shiftMC, weight,
110  "mu -> e", "#mu #rightarrow e",
111  NumuNDcut, NumuOscDecomp, // numu selection in ND
112  DecompResult::numu, kIsNumuCC && !kIsAntiNu, // numu truth in ND
113  farMCswap, kIsSig && !kIsAntiNu // mu->e in FD
114  )
115  );
116  extrap.fMEAntiextrap = std::unique_ptr<ModularExtrapComponent>(
117  new TruthReweight(
118  nearMC, axis, axisNumuND, fdcut, shiftMC, weight,
119  "mubar -> ebar", "#bar{#mu} #rightarrow #bar{e}",
120  NumuNDcut, NumuOscDecomp, // numu selection in ND
121  DecompResult::numubar, kIsNumuCC && kIsAntiNu, // numubar truth in ND
122  farMCswap, kIsSig && kIsAntiNu // mubar->ebar in FD
123  )
124  );
125 
126  // e -> e ----
127  extrap.fEEextrap = std::unique_ptr<ModularExtrapComponent>(
128  new RecoReweight(
129  nearMC, axis, fdcut, shiftMC, weight,
130  "e -> e", "e #rightarrow e",
131  NCNDcut, NCSurvDecomp, // NC selection in ND
132  DecompResult::nue, kIsBeamNue && !kIsAntiNu, // nue truth in ND
133  farMCnonswap, kIsBeamNue && !kIsAntiNu // e->e in FD
134  )
135  );
136 
137  return extrap;
138  }
139 
140  //---------------------------------------------------------------------------
142  Loaders& loaders,
143  const HistAxis& axis,
144  const Cut& fdcut,
145  const SystShifts& shiftMC,
146  const Var& weight
147  ) {
149  loaders.GetLoader(
151  loaders.GetLoader(
153  loaders.GetLoader(
155  axis,
156  fdcut,
157  shiftMC,
158  weight
159  );
160  }
161 
162  //---------------------------------------------------------------------------
164  SpectrumLoaderBase& farMCswap,
165  SpectrumLoaderBase& farMCnonswap,
166  SpectrumLoaderBase& farMCtauswap,
167  const HistAxis& axis,
168  const Cut& fdcut,
169  const SystShifts& shiftMC,
170  const Var& weight
171  ) {
172  // Initialize all components with NoReweights,
173  // make sure other objects exist
174  ModularExtrapSterile extrap(
175  farMCswap,
176  farMCnonswap,
177  farMCtauswap,
178  axis,
179  fdcut,
180  shiftMC,
181  weight
182  );
183 
184  return extrap;
185  }
186 
187  //---------------------------------------------------------------------------
189  SpectrumLoaderBase& farMCswap,
190  SpectrumLoaderBase& farMCnonswap,
191  SpectrumLoaderBase& farMCtauswap,
192  const HistAxis& axis,
193  const Cut& fdcut,
194  const SystShifts& shiftMC,
195  const Var& weight
196  ) : ModularExtrap(
197  farMCswap,
198  farMCnonswap,
199  farMCtauswap,
200  axis,
201  fdcut,
202  shiftMC,
203  weight
204  ), // Sets up the components with default NoReweights
205 
206  // Set up Spectrum objects for NC Proportion Ratios
207  fNCNueNumerator(std::make_unique<Spectrum>(farMCnonswap, axis,
208  fdcut && kIsNC && !kIsAntiNu && kNCBeamNue,
209  shiftMC, weight)),
210  fNCAntiNueNumerator(std::make_unique<Spectrum>(farMCnonswap, axis,
211  fdcut && kIsNC && kIsAntiNu && kNCBeamNue,
212  shiftMC, weight)),
213  fNCNumuNumerator(std::make_unique<Spectrum>(farMCnonswap, axis,
214  fdcut && kIsNC && !kIsAntiNu && kNCBeamNumu,
215  shiftMC, weight)),
216  fNCAntiNumuNumerator(std::make_unique<Spectrum>(farMCnonswap, axis,
217  fdcut && kIsNC && kIsAntiNu && kNCBeamNumu,
218  shiftMC, weight)),
219  fNCDenominator(std::make_unique<Spectrum>(farMCnonswap, axis,
220  fdcut && kIsNC,
221  shiftMC, weight)) {
222  }
223 
224  //---------------------------------------------------------------------------
226  : ModularExtrap(std::move(load)) { // Sets up the ModularExtrapComponents
227  }
228 
229  //---------------------------------------------------------------------------
231  return fNCextrap->Return();
232  }
233 
234  //---------------------------------------------------------------------------
236  return OscNCComponent().Unoscillated();
237  }
238 
239  //---------------------------------------------------------------------------
241  return Ratio(*fNCNueNumerator.get(), *fNCDenominator.get());
242  }
243 
244  //---------------------------------------------------------------------------
246  return Ratio(*fNCAntiNueNumerator.get(), *fNCDenominator.get());
247  }
248 
249  //---------------------------------------------------------------------------
251  return Ratio(*fNCNumuNumerator.get(), *fNCDenominator.get());
252  }
253 
254  //---------------------------------------------------------------------------
256  return Ratio(*fNCAntiNumuNumerator.get(), *fNCDenominator.get());
257  }
258 
259  //---------------------------------------------------------------------------
260  void ModularExtrapSterile::SaveTo(TDirectory* dir, const std::string& name) const {
261  TDirectory* tmp = gDirectory;
262 
263  dir = dir->mkdir(name.c_str()); // switch to subdir
264  dir->cd();
265 
266  ModularExtrap::SaveTo(dir, name);
267 
268  dir->cd();
269  TObjString("ModularExtrapSterile").Write("type", TObject::kOverwrite);
270 
271  fNCNueNumerator->SaveTo(dir, "NCNue");
272  fNCAntiNueNumerator->SaveTo(dir, "NCAntiNue");
273  fNCNumuNumerator->SaveTo(dir, "NCNumu");
274  fNCAntiNumuNumerator->SaveTo(dir, "NCAntiNumu");
275  fNCDenominator->SaveTo(dir, "AllNC");
276 
277  dir->Write();
278  delete dir;
279 
280  tmp->cd();
281  }
282 
283  //---------------------------------------------------------------------------
284  std::unique_ptr<ModularExtrapSterile> ModularExtrapSterile::LoadFrom(TDirectory* dir, const std::string& name) {
285  std::unique_ptr<ModularExtrapSterile> ret(new ModularExtrapSterile( std::move( *ModularExtrap::LoadFrom(dir, name).release() ) )
286  );
287 
288  dir = dir->GetDirectory(name.c_str()); // switch to subdir
289  assert(dir);
290 
291  ret->fNCNueNumerator = Spectrum::LoadFrom(dir, "NCNue");
292  ret->fNCAntiNueNumerator = Spectrum::LoadFrom(dir, "NCAntiNue");
293  ret->fNCNumuNumerator = Spectrum::LoadFrom(dir, "NCNumu");
294  ret->fNCAntiNumuNumerator = Spectrum::LoadFrom(dir, "NCAntiNumu");
295  ret->fNCDenominator = Spectrum::LoadFrom(dir, "AllNC");
296 
297  delete dir;
298 
299  return ret;
300  }
301 
302  //---------------------------------------------------------------------------
303  bool NCFlavSel::operator()(const caf::SRProxy* sr) const {
304  if(sr->mc.nnu == 0) return false;
305  assert(sr->mc.nnu == 1);
306  return (!sr->mc.nu[0].iscc &&
307  abs(sr->mc.nu[0].pdg) == fPdg &&
308  abs(sr->mc.nu[0].pdgorig) == fPdgOrig);
309  }
310 }
Near Detector underground.
Definition: SREnums.h:10
const XML_Char * name
Definition: expat.h:151
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
std::unique_ptr< Spectrum > fNCNueNumerator
Far Detector at Ash River.
Definition: SREnums.h:11
std::unique_ptr< ModularExtrapComponent > fMMextrap
Definition: ModularExtrap.h:92
void SaveTo(TDirectory *dir, const std::string &name) const override
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Ratio NCAntiNumuProportion() const
const Var weight
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
std::unique_ptr< Spectrum > fNCNumuNumerator
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2085
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
std::unique_ptr< ModularExtrapComponent > fMEextrap
Definition: ModularExtrap.h:94
std::unique_ptr< Spectrum > fNCAntiNumuNumerator
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:574
void load(std::string lib)
Definition: load_libs.C:3
const Cut kIsAntiNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return sr->mc.nu[0].pdg< 0;})
Is this truly an antineutrino?
Definition: TruthCuts.h:53
const Cut kIsBeamNue(CCFlavSel(12, 12))
Select CC .
Ratio NCNueProportion() const
Return the proportion of NCs that originate as a given neutrino flavor/sign.
void abs(TH1 *hist)
caf::Proxy< short int > nnu
Definition: SRProxy.h:573
OscillatableSpectrum OscNCComponent() const
Return the oscillatable NC component.
Float_t tmp
Definition: plot.C:36
std::unique_ptr< Spectrum > fNCDenominator
Ratio NCNumuProportion() const
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
std::unique_ptr< ModularExtrapComponent > fNCextrap
const Cut kNCBeamNue(NCFlavSel(12, 12))
Select NC .
ModularExtrapSterile(SpectrumLoaderBase &farMCswapLoader, SpectrumLoaderBase &farMCnonswapLoader, SpectrumLoaderBase &farMCtauswapLoader, const HistAxis &axis, const Cut &fdcut, const SystShifts &shiftMC, const Var &weight)
Sets up all components to use FD MC–internal use only.
Ratio NCAntiNueProportion() const
std::unique_ptr< ModularExtrapComponent > fMEAntiextrap
Definition: ModularExtrap.h:95
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:546
const Cut kIsNC([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return!sr->mc.nu[0].iscc;})
Is this a Neutral Current event?
Definition: TruthCuts.h:8
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
Definition: Loaders.cxx:129
bool operator()(const caf::SRProxy *sr) const
Extrapolates component using truth-over-truth method.
const Cut kIsSig(CCFlavSel(12, 14))
Select CC .
caf::StandardRecord * sr
const Cut kNCBeamNumu(NCFlavSel(14, 14))
Select NC .
Spectrum NCTotalComponent() override
Override the ModularExtrap method.
Represent the ratio between two spectra.
Definition: Ratio.h:8
static std::unique_ptr< ModularExtrap > LoadFrom(TDirectory *dir, const std::string &name)
Base class for the various types of spectrum loader.
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2097
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TDirectory * dir
Definition: macro.C:5
std::unique_ptr< Spectrum > fNCAntiNueNumerator
REGISTER_LOADFROM("BENDecomp", IDecomp, BENDecomp)
std::vector< Loaders * > loaders
Definition: syst_header.h:386
Standard interface to all decomposition techniques.
Definition: IDecomp.h:13
assert(nhit_max >=nhit_nbins)
static std::unique_ptr< ModularExtrapSterile > LoadFrom(TDirectory *dir, const std::string &name)
void SaveTo(TDirectory *dir, const std::string &name) const override
Spectrum with true energy information, allowing it to be oscillated
Extrapolate each component using a separate ModularExtrapComponent.
Definition: ModularExtrap.h:23
std::unique_ptr< ModularExtrapComponent > fEEextrap
Definition: ModularExtrap.h:90
static ModularExtrapSterile NCDisappearance(Loaders &loaders, const IDecomp &NCSurvDecomp, const IDecomp &NumuOscDecomp, const HistAxis &axis, const HistAxis &axisNumuND, const Cut &fdcut, const Cut &NCNDcut, const Cut &NumuNDcut, const SystShifts &shiftMC=kNoShift, const Var &weight=kUnweighted)
Creates a NC disappearance extrapolation.
static ModularExtrapSterile TrivialExtrapNC(Loaders &loaders, const HistAxis &axis, const Cut &fdcut, const SystShifts &shiftMC=kNoShift, const Var &weight=kUnweighted)
Creates a NC disappearance "extrapolation" by returning FD MC.
Extrapolates using reco-over-reco method.