AnaResultsLoad.C
Go to the documentation of this file.
1 // This macro creates spectra for OBSERVED
2 // and predicted event rates in analysis region
3 // This uses UNBLIND DATA. Proceed with extreme caution!
4 
8 #include "NuXAna/Cuts/NusCuts.h"
10 #include "CAFAna/Cuts/TimingCuts.h"
12 #include "CAFAna/Extrap/ExtrapSterile.h"
16 #include "CAFAna/Vars/HistAxes.h"
19 #include "NuXAna/Vars/HistAxes.h"
20 
21 #include <string>
22 #include <utility>
23 
24 using namespace ana;
25 
27 {
28  TH1::AddDirectory(0);
29 
30  // ND DATA + MC
31  SpectrumLoader loaderNDMC(fnamenear_concat);
33  loaderNDMC.SetSpillCut(kStandardSpillCuts);
34  loaderNDdata.SetSpillCut(kStandardSpillCuts);
35 
36  // PERIOD 1, 2 + EPOCH 3B
37  SpectrumLoader loaderFDMC_non(fFDMC_non);
38  SpectrumLoader loaderFDMC_swp(fFDMC_swp);
39  SpectrumLoader loaderFDMC_tau(fFDMC_tau);
40  loaderFDMC_non.SetSpillCut(kStandardSpillCuts);
41  loaderFDMC_swp.SetSpillCut(kStandardSpillCuts);
42  loaderFDMC_tau.SetSpillCut(kStandardSpillCuts);
43 
44  // PERIOD 3C (and implied 3D)
45  SpectrumLoader loaderFDMC_non_3c(fFDMC_non_3c);
46  SpectrumLoader loaderFDMC_swp_3c(fFDMC_swp_3c);
47  SpectrumLoader loaderFDMC_tau_3c(fFDMC_tau_3c);
48  loaderFDMC_non_3c.SetSpillCut(kStandardSpillCuts);
49  loaderFDMC_swp_3c.SetSpillCut(kStandardSpillCuts);
50  loaderFDMC_tau_3c.SetSpillCut(kStandardSpillCuts);
51 
52  // FD DATA
54  loaderFDdata.SetSpillCut(kStandardSpillCuts);
55 
56  std::string labelVtxX = "X Vertex Position (cm)";
57  std::string labelVtxY = "Y Vertex Position (cm)";
58  std::string labelVtxZ = "Z Vertex Position (cm)";
59  std::string labelHits = "Number of Hits";
60  std::string labelCVN = "CVN NCID";
61  std::string labelNCP = "NumucontPID";
62  std::string labelPTP = "PT/P";
63  std::string labelEPerHit = "Average E/Hit (GeV)";
64  std::string labelTop = "Dist. to Det. Top (cm)";
65 
66  Binning kXYBinsFD = Binning::Simple(64, -800., 800.);
67  Binning kZBinsFD = Binning::Simple(60, 0., 6000.);
68  Binning kNearEdgeBins = Binning::Simple(64, 0., 1600.);
69  Binning kFatPIDBins = Binning::Simple(25, 0., 1.);
70  Binning kHitBins = Binning::Simple(50, 0., 200.);
71  Binning kEDepBins = Binning::Simple(25, 0., .1);
72 
73  const Var kRecoVtxX(
74  [](const caf::SRProxy* sr)
75  {
76  if(!sr->vtx.elastic.IsValid) { return 9000.; }
77  return sr->vtx.elastic.vtx.x;
78  });
79  const Var kRecoVtxY(
80  [](const caf::SRProxy* sr)
81  {
82  if(!sr->vtx.elastic.IsValid) { return 9000.; }
83  return sr->vtx.elastic.vtx.y;
84  });
85  const Var kRecoVtxZ(
86  [](const caf::SRProxy* sr)
87  {
88  if(!sr->vtx.elastic.IsValid) { return 9000.; }
89  return sr->vtx.elastic.vtx.z;
90  });
91 
92  const Var kNCP = SIMPLEVAR(sel.cosrej.numucontpid);
93  const Var kPartPTP(
94  [](const caf::SRProxy* sr)
95  {
96  if(std::isnan(sr->sel.nuecosrej.partptp)) { return -5.f; }
97  return sr->sel.nuecosrej.partptp;
98  });
99  const Var kEPerHit(
100  [](const caf::SRProxy* sr)
101  {
102  double nhit = (double)sr->slc.nhit;
103  if(nhit <= 0.) { return 0.; }
104  return sr->slc.calE/nhit;
105  });
106  const Var kDistTop(
107  [](const caf::SRProxy* sr)
108  {
110  });
111 
112  std::map<std::string, HistAxis*> axes;
113  axes["CalE"] = new HistAxis(kNCAxis);
114  axes["VtxX"] = new HistAxis(labelVtxX, kXYBinsFD, kRecoVtxX);
115  axes["VtxY"] = new HistAxis(labelVtxY, kXYBinsFD, kRecoVtxY);
116  axes["VtxZ"] = new HistAxis(labelVtxZ, kZBinsFD, kRecoVtxZ);
117  axes["CVN"] = new HistAxis(labelCVN, kFatPIDBins, kCVNnc);
118  axes["NHit"] = new HistAxis(labelHits, kHitBins, kNHit);
119  axes["CosBDT"] = new HistAxis(labelNCP, kFatPIDBins, kNCP);
120  axes["PTP"] = new HistAxis(labelPTP, kFatPIDBins, kPartPTP);
121  axes["EperHit"] = new HistAxis(labelEPerHit, kEDepBins, kEPerHit);
122  axes["DistTop"] = new HistAxis(labelTop, kNearEdgeBins, kDistTop);
123 
124  std::map<std::string, IDecomp*> ncdecomps;
125  std::map<std::string, IDecomp*> numudecomps;
126  std::map<std::string, ModularExtrapSterile*> extrap123bs;
127  std::map<std::string, ModularExtrapSterile*> extrap3c3ds;
128  std::map<std::string, PredictionSterile*> predst123bs;
129  std::map<std::string, PredictionSterile*> predst3c3ds;
130  std::map<std::string, PredictionCombinePeriods*> preds;
131  std::map<std::string, Spectrum*> cosmics;
132  std::map<std::string, Spectrum*> dataspecs;
133 
134  for(const auto& ax : axes) {
135  std::string axlabel = ax.first;
136 
137  ncdecomps[axlabel] = new ProportionalDecomp(
138  loaderNDMC, loaderNDdata,
139  *ax.second, kNusND,
141  );
142 
143  numudecomps[axlabel] = new ProportionalDecomp(
144  loaderNDMC, loaderNDdata,
147  );
148 
149  extrap123bs[axlabel] = new ModularExtrapSterile(
151  loaderNDMC, loaderFDMC_swp, loaderFDMC_non, loaderFDMC_tau,
152  *ncdecomps[axlabel], *numudecomps[axlabel],
153  *ax.second, kNCBinsNumuCCAxis,
156  )
157  );
158 
159  extrap3c3ds[axlabel] = new ModularExtrapSterile(
161  loaderNDMC, loaderFDMC_swp_3c, loaderFDMC_non_3c, loaderFDMC_tau_3c,
162  *ncdecomps[axlabel], *numudecomps[axlabel],
163  *ax.second, kNCBinsNumuCCAxis,
166  )
167  );
168 
169  predst123bs[axlabel] = new PredictionSterile(extrap123bs[axlabel]);
170  predst3c3ds[axlabel] = new PredictionSterile(extrap3c3ds[axlabel]);
171 
172  cosmics[axlabel] = new Spectrum(loaderFDdata, *ax.second,
174  kNoShift, kTimingSidebandWeight);
175 
176  dataspecs[axlabel] = new Spectrum(loaderFDdata, *ax.second, kInBeamSpill && kNusFD);
177  }
178 
179  // Fill the spectra!
180  loaderNDMC.Go();
181  loaderNDdata.Go();
182  loaderFDMC_non.Go();
183  loaderFDMC_swp.Go();
184  loaderFDMC_tau.Go();
185  loaderFDMC_non_3c.Go();
186  loaderFDMC_swp_3c.Go();
187  loaderFDMC_tau_3c.Go();
188  loaderFDdata.Go();
189 
190  // Two scale factors for scaling spectra, 1 data yr and 3 data yr
192  const double pot_e3ce3d = kSecondAnaEpoch3cPOT+kSecondAnaEpoch3dPOT;
193 
194  for(const auto& ax : axes) {
195  std::string axlabel = ax.first;
196 
197  preds[axlabel] = new PredictionCombinePeriods({
198  std::make_pair(predst123bs[axlabel], pot_p1p2e3b),
199  std::make_pair(predst3c3ds[axlabel], pot_e3ce3d)
200  });
201  }
202 
203  // Set up output filename
204  TFile* rootF = new TFile(outfile.c_str(), "RECREATE");
205 
206  // Save all of the objects
207  TDirectory* tmp = gDirectory;
208  TDirectory* saveDir = gDirectory;
210  std::string sep = "__";
211 
212  for(const auto& ax : axes) {
213  std::string axlabel = ax.first;
214 
215  dir = "ncdecomp" + sep + axlabel;
216  saveDir = rootF->mkdir(dir.c_str());
217  ncdecomps[axlabel]->SaveTo(saveDir);
218 
219  dir = "numudecomp" + sep + axlabel;
220  saveDir = rootF->mkdir(dir.c_str());
221  numudecomps[axlabel]->SaveTo(saveDir);
222 
223  dir = "extrap123b" + sep + axlabel;
224  saveDir = rootF->mkdir(dir.c_str());
225  extrap123bs[axlabel]->SaveTo(saveDir);
226 
227  dir = "extrap3c3d" + sep + axlabel;
228  saveDir = rootF->mkdir(dir.c_str());
229  extrap3c3ds[axlabel]->SaveTo(saveDir);
230 
231  dir = "predst123b" + sep + axlabel;
232  saveDir = rootF->mkdir(dir.c_str());
233  predst123bs[axlabel]->SaveTo(saveDir);
234 
235  dir = "predst3c3d" + sep + axlabel;
236  saveDir = rootF->mkdir(dir.c_str());
237  predst3c3ds[axlabel]->SaveTo(saveDir);
238 
239  dir = "pred" + sep + axlabel;
240  saveDir = rootF->mkdir(dir.c_str());
241  preds[axlabel]->SaveTo(saveDir);
242 
243  dir = "cosmic" + sep + axlabel;
244  saveDir = rootF->mkdir(dir.c_str());
245  cosmics[axlabel]->SaveTo(saveDir);
246 
247  dir = "dataspec" + sep + axlabel;
248  saveDir = rootF->mkdir(dir.c_str());
249  dataspecs[axlabel]->SaveTo(saveDir);
250  }
251 
252  tmp->cd();
253  rootF->Close(); // Close the file
254 }
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
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
caf::Proxy< float > stoptop
Definition: SRProxy.h:1036
const double kSecondAnaPeriod2POT
Definition: Exposures.h:74
const Var kCVNnc
PID
Definition: Vars.cxx:44
const Var kDistTop([](const caf::SRProxy *sr){return std::min(sr->sel.nuecosrej.starttop, sr->sel.nuecosrej.stoptop);})
Definition: NusVars.h:76
const Cut kNusFD
Definition: NusCuts.h:46
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2085
const double kSecondAnaEpoch3dPOT
Definition: Exposures.h:77
const double kSecondAnaEpoch3bPOT
Definition: Exposures.h:75
nhit
Definition: demo1.py:25
const std::vector< std::string > fFDMC_non
Float_t tmp
Definition: plot.C:36
void SetSpillCut(const SpillCut &cut)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
const double kSecondAnaPeriod1POT
Definition: Exposures.h:73
const Var kEPerHit([](const caf::SRProxy *sr){if(sr->slc.nhit >0) return 1000.0 *(sr->slc.calE/sr->slc.nhit);else return-5.;})
Definition: NueVarsExtra.h:14
caf::Proxy< float > starttop
Definition: SRProxy.h:1030
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
const Cut kInTimingSideband([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return(kInTimingSideband_before(sr)|| kInTimingSideband_after(sr));else return(kInTimingSideband_before(sr)|| kInTimingSideband_afterA(sr)|| kInTimingSideband_afterB(sr));}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return(kInTimingSideband_before.Livetime(spill)+ kInTimingSideband_after.Livetime(spill));else return(kInTimingSideband_before.Livetime(spill)+ kInTimingSideband_afterA.Livetime(spill)+ kInTimingSideband_afterB.Livetime(spill));}, [](const caf::SRSpillProxy *spill){return 0;})
Definition: TimingCuts.h:12
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2077
caf::Proxy< caf::SRNueCosRej > nuecosrej
Definition: SRProxy.h:1221
const std::string fnamenear_concat
const Var kNCP
caf::Proxy< unsigned int > nhit
Definition: SRProxy.h:1271
const Binning kXYBinsFD
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
caf::Proxy< float > z
Definition: SRProxy.h:107
const std::string fnameneardata_concat
const Binning kZBinsFD
const std::vector< std::string > fFDMC_swp
caf::Proxy< float > x
Definition: SRProxy.h:105
const Var kNHit
Definition: Vars.cxx:71
const std::string fFDMC_non_3c
const Cut kNumuND
Definition: NumuCuts.h:55
void AnaResultsLoad(std::string outfile)
caf::StandardRecord * sr
caf::Proxy< float > partptp
Definition: SRProxy.h:1013
virtual void Go() override
Load all the registered spectra.
const HistAxis kNCAxis("Calorimetric Energy (GeV)", kNCDisappearanceEnergyBinning, kCaloE)
Axes used in Ana01 analysis by nus group.
Definition: HistAxes.h:8
std::vector< float > Spectrum
Definition: Constants.h:570
const char sep
const SystShifts kNoShift
Definition: SystShifts.cxx:21
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const std::string fFDMC_tau_3c
const HistAxis kNCBinsNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNCDisappearanceEnergyBinning, kCCE)
Definition: HistAxes.h:9
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2017
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2101
Splits Data proportionally according to MC.
const Var kTuftsWeightCC
Definition: XsecTunes.h:31
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
caf::Proxy< float > y
Definition: SRProxy.h:106
TDirectory * dir
Definition: macro.C:5
caf::Proxy< float > calE
Definition: SRProxy.h:1248
const std::vector< std::string > fFDMC_tau
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2032
const std::string fFDMC_swp_3c
const Cut kInBeamSpill([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return kInBeamSpill_main(sr);else return kInBeamSpill_main(sr)||kInBeamSpill_shifted(sr);}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return kInBeamSpill_main.Livetime(spill);else return kInBeamSpill_main.Livetime(spill)+kInBeamSpill_shifted.Livetime(spill);}, [](const caf::SRSpillProxy *spill) -> double{return spill->spillpot;})
Does the event fall inside the window we call the beam spill?
Definition: TimingCuts.h:8
A prediction object compatible with sterile oscillations.
const double kSecondAnaEpoch3cPOT
Definition: Exposures.h:76
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Sum MC predictions from different periods scaled according to data POT targets.
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2100
T min(const caf::Proxy< T > &a, T b)
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2105
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 Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
const Cut kNusND
Definition: NusCuts.h:71
FILE * outfile
Definition: dump_event.C:13
const std::vector< std::string > fnamefardata_unblind(MakeUnblindList())