SideBandLoad.C
Go to the documentation of this file.
1 // This macro creates spectra for OBSERVED
2 // and predicted event rates in sideband regions
3 // This uses UNBLIND DATA. Edit with extreme caution.
4 
9 #include "NuXAna/Cuts/NusCuts.h"
10 #include "CAFAna/Cuts/SpillCuts.h"
11 #include "CAFAna/Cuts/TimingCuts.h"
13 #include "CAFAna/Extrap/ExtrapSterile.h"
17 #include "CAFAna/Vars/HistAxes.h"
18 #include "NuXAna/Vars/HistAxes.h"
21 
22 #include <string>
23 #include <utility>
24 
25 using namespace ana;
26 
27 struct Sample {
28  HistAxis* axnc;
29  HistAxis* axnumu;
30  Cut* nd;
31  Cut* fd;
32 };
33 
34 Sample MakeSample(const HistAxis& axnc, const HistAxis& axnm,
35  const Cut& nd, const Cut& fd);
36 
37 std::vector<std::string> MakeUnblindList();
38 
40 {
41  TH1::AddDirectory(0);
42 
43  // ND DATA + MC
44  SpectrumLoader loaderNDMC(fnamenear_concat);
46  loaderNDMC.SetSpillCut(kStandardSpillCuts);
47  loaderNDdata.SetSpillCut(kStandardSpillCuts);
48 
49  // PERIOD 1, 2 + EPOCH 3B
50  SpectrumLoader loaderFDMC_non(fFDMC_non);
51  SpectrumLoader loaderFDMC_swp(fFDMC_swp);
52  SpectrumLoader loaderFDMC_tau(fFDMC_tau);
53  loaderFDMC_non.SetSpillCut(kStandardSpillCuts);
54  loaderFDMC_swp.SetSpillCut(kStandardSpillCuts);
55  loaderFDMC_tau.SetSpillCut(kStandardSpillCuts);
56 
57  // PERIOD 3C (and implied 3D)
58  SpectrumLoader loaderFDMC_non_3c(fFDMC_non_3c);
59  SpectrumLoader loaderFDMC_swp_3c(fFDMC_swp_3c);
60  SpectrumLoader loaderFDMC_tau_3c(fFDMC_tau_3c);
61  loaderFDMC_non_3c.SetSpillCut(kStandardSpillCuts);
62  loaderFDMC_swp_3c.SetSpillCut(kStandardSpillCuts);
63  loaderFDMC_tau_3c.SetSpillCut(kStandardSpillCuts);
64 
65  // FD DATA
66  const std::vector<std::string> fnameunblind(MakeUnblindList());
67  SpectrumLoader loaderFDdata(fnameunblind);
68  loaderFDdata.SetSpillCut(kStandardSpillCuts);
69 
70  const Binning kNCHighEBins = Binning::Simple(8, 4, 6);
71  const HistAxis kNCHighEAxis( "Calorimetric Energy (GeV)",
72  kNCHighEBins, kCaloE );
73  const HistAxis kHighENumuCCAxis( "Reconstructed Neutrino Energy (GeV)",
74  kNCHighEBins, kCCE );
75 
76  std::string labelRecoE = "Calorimetric Energy (GeV)";
77  const Cut kLowCVN = kNHit < 200 && kCVNnc < 0.2;
78  const Cut kMidBDT({
79  "sel.cosrej.numucontpid", "sel.nuecosrej.partptp",
80  "vtx.elastic.fuzzyk.nshwlid", "vtx.elastic.fuzzyk.png.shwlid.lid.ismuon",
81  "slc.calE", "slc.nhit","sel.nuecosrej.starttop",
82  "sel.nuecosrej.stoptop"
83  },
84  [](const caf::SRProxy* sr)
85  {
86  double numucontpid = sr->sel.cosrej.numucontpid;
87  double partptp = sr->sel.nuecosrej.partptp;
88  if(sr->vtx.elastic.fuzzyk.nshwlid == 0) return false;
89  if(numucontpid > 0.5) return false;
90  if(numucontpid <= 0.42) return false;
91  if(partptp >= 0.8) return false;
92  if(sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.ismuon == 1) return false;
93  double nhit = (double)sr->slc.nhit;
94  if(nhit <= 0.) return false;
95  if(sr->slc.calE/nhit <= 0.018) return false;
96  if(std::min(sr->sel.nuecosrej.starttop, sr->sel.nuecosrej.stoptop) < 480) return false;
97  return true;
98  });
99 
100  std::map<std::string, Sample> sidebands;
101  sidebands["HighE"] = MakeSample(
102  kNCHighEAxis, kHighENumuCCAxis,
103  kNusNDAllE && kCaloE > 4.,
104  kNusFDAllE && kCaloE > 4.
105  );
106  sidebands["LowCVN"] = MakeSample(
108  kNusNDPresel && kNusCosRejMod && kLowCVN && kNCCalERange,
109  kNusFDPresel && kNusCosRej && kLowCVN && kNCCalERange
110  );
111  sidebands["MidBDT"] = MakeSample(
113  kNusND,
114  kNusFDPresel && kNusNCSel && kMidBDT && kNCCalERange
115  );
116 
117  std::map<std::string, IDecomp*> ncdecomps;
118  std::map<std::string, IDecomp*> numudecomps;
119  std::map<std::string, ModularExtrapSterile*> extrap123bs;
120  std::map<std::string, ModularExtrapSterile*> extrap3c3ds;
121  std::map<std::string, PredictionSterile*> predst123bs;
122  std::map<std::string, PredictionSterile*> predst3c3ds;
123  std::map<std::string, PredictionCombinePeriods*> preds;
124  std::map<std::string, Spectrum*> cosmics;
125  std::map<std::string, Spectrum*> dataspecs;
126 
127  for(const auto& sideband : sidebands) {
128  std::string samplelabel = sideband.first;
129  Sample sample = sideband.second;
130 
131  ncdecomps[samplelabel] = new ProportionalDecomp(
132  loaderNDMC, loaderNDdata,
133  *sample.axnc, *sample.nd,
135  );
136 
137  numudecomps[samplelabel] = new ProportionalDecomp(
138  loaderNDMC, loaderNDdata,
139  *sample.axnumu, kNumuND,
141  );
142 
143  extrap123bs[samplelabel] = new ModularExtrapSterile(
145  loaderNDMC, loaderFDMC_swp, loaderFDMC_non, loaderFDMC_tau,
146  *ncdecomps[samplelabel], *numudecomps[samplelabel],
147  *sample.axnc, *sample.axnumu,
148  *sample.fd, *sample.nd, kNumuND,
150  )
151  );
152 
153  extrap3c3ds[samplelabel] = new ModularExtrapSterile(
155  loaderNDMC, loaderFDMC_swp_3c, loaderFDMC_non_3c, loaderFDMC_tau_3c,
156  *ncdecomps[samplelabel], *numudecomps[samplelabel],
157  *sample.axnc, *sample.axnumu,
158  *sample.fd, *sample.nd, kNumuND,
160  )
161  );
162 
163  predst123bs[samplelabel] = new PredictionSterile(extrap123bs[samplelabel]);
164  predst3c3ds[samplelabel] = new PredictionSterile(extrap3c3ds[samplelabel]);
165 
166  cosmics[samplelabel] = new Spectrum(loaderFDdata, *sample.axnc,
167  kInTimingSideband && *sample.fd,
168  kNoShift, kTimingSidebandWeight);
169 
170  dataspecs[samplelabel] = new Spectrum(loaderFDdata, *sample.axnc, kInBeamSpill && *sample.fd);
171  }
172 
173  // Fill the spectra!
174  loaderNDMC.Go();
175  loaderNDdata.Go();
176  loaderFDMC_non.Go();
177  loaderFDMC_swp.Go();
178  loaderFDMC_tau.Go();
179  loaderFDMC_non_3c.Go();
180  loaderFDMC_swp_3c.Go();
181  loaderFDMC_tau_3c.Go();
182  loaderFDdata.Go();
183 
184  // Two scale factors for scaling spectra, 1 data yr and 3 data yr
186  const double pot_e3ce3d = kSecondAnaEpoch3cPOT+kSecondAnaEpoch3dPOT;
187 
188  for(const auto& sideband : sidebands) {
189  std::string samplelabel = sideband.first;
190  preds[samplelabel] = new PredictionCombinePeriods({
191  std::make_pair(predst123bs[samplelabel], pot_p1p2e3b),
192  std::make_pair(predst3c3ds[samplelabel], pot_e3ce3d)
193  });
194  }
195 
196  // Set up output filename
197  TFile* rootF = new TFile(outfile.c_str(), "RECREATE");
198 
199  // Save all of the objects
200  TDirectory* tmp = gDirectory;
201  TDirectory* saveDir = gDirectory;
203  std::string sep = "__";
204 
205  for(const auto& cafanaobj : ncdecomps) {
206  dir = "ncdecomp" + sep + cafanaobj.first; // sample name
207  saveDir = rootF->mkdir(dir.c_str());
208  cafanaobj.second->SaveTo(saveDir);
209  }
210 
211  for(const auto& cafanaobj : numudecomps) {
212  dir = "numudecomp" + sep + cafanaobj.first; // sample name
213  saveDir = rootF->mkdir(dir.c_str());
214  cafanaobj.second->SaveTo(saveDir);
215  }
216 
217  for(const auto& cafanaobj : extrap123bs) {
218  dir = "extrap123b" + sep + cafanaobj.first; // sample name
219  saveDir = rootF->mkdir(dir.c_str());
220  cafanaobj.second->SaveTo(saveDir);
221  }
222 
223  for(const auto& cafanaobj : extrap3c3ds) {
224  dir = "extrap3c3d" + sep + cafanaobj.first; // sample name
225  saveDir = rootF->mkdir(dir.c_str());
226  cafanaobj.second->SaveTo(saveDir);
227  }
228 
229  for(const auto& cafanaobj : predst123bs) {
230  dir = "predst123b" + sep + cafanaobj.first; // sample name
231  saveDir = rootF->mkdir(dir.c_str());
232  cafanaobj.second->SaveTo(saveDir);
233  }
234 
235  for(const auto& cafanaobj : predst3c3ds) {
236  dir = "predst3c3d" + sep + cafanaobj.first; // sample name
237  saveDir = rootF->mkdir(dir.c_str());
238  cafanaobj.second->SaveTo(saveDir);
239  }
240 
241  for(const auto& cafanaobj : preds) {
242  dir = "pred" + sep + cafanaobj.first; // sample name
243  saveDir = rootF->mkdir(dir.c_str());
244  cafanaobj.second->SaveTo(saveDir);
245  }
246 
247  for(const auto& cafanaobj : cosmics) {
248  dir = "cosmic" + sep + cafanaobj.first; // sample name
249  saveDir = rootF->mkdir(dir.c_str());
250  cafanaobj.second->SaveTo(saveDir);
251  }
252 
253  for(const auto& cafanaobj : dataspecs) {
254  dir = "dataspec" + sep + cafanaobj.first; // sample name
255  saveDir = rootF->mkdir(dir.c_str());
256  cafanaobj.second->SaveTo(saveDir);
257  }
258 
259  tmp->cd();
260  rootF->Close(); // Close the file
261 }
262 
263 Sample MakeSample(const HistAxis& axnc, const HistAxis& axnm,
264  const Cut& nd, const Cut& fd)
265 {
266  Sample ret;
267  ret.axnc = new HistAxis(axnc.label, axnc.bins, axnc.var);
268  ret.axnumu = new HistAxis(axnm.label, axnm.bins, axnm.var);
269  ret.nd = new Cut(nd);
270  ret.fd = new Cut(fd);
271 
272  return ret;
273 }
274 
275 std::vector<std::string> MakeUnblindList()
276 {
277  std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
278  std::cout << " WARNING: YOU ARE USING UNBLIND DATA." << std::endl;
279  std::cout << " Be very cautious about what distributions you are making," << std::endl;
280  std::cout << " where the results are going, and who will see them." << std::endl;
281  std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
282 
283  std::vector<std::string> fname1 = Wildcard(
284  diskdirFD + "by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_period1_nus_contain_v1_goodruns_prod2-snapshot/"
285  + "prod_restricteddecaf_S16-05-20_fd_numi_fhc_period1_nus_contain_v1_goodruns_prod2-snapshot*.root"
286  );
287  std::vector<std::string> fname2 = Wildcard(
288  diskdirFD + "by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_period2_nus_contain_v1_goodruns_prod2-snapshot/"
289  + "prod_restricteddecaf_S16-05-20_fd_numi_fhc_period2_nus_contain_v1_goodruns_prod2-snapshot*.root"
290  );
291  std::vector<std::string> fname3b = Wildcard(
292  diskdirFD + "by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3*_nus_contain_v1_goodruns_prod2-snapshot.root"
293  );
294  std::vector<std::string> fname3c = Wildcard(
295  diskdirFD + "by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3c_nus_contain_v1_goodruns_prod2-snapshot/"
296  + "prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3c_nus_contain_v1_goodruns_prod2-snapshot*.root"
297  );
298  std::vector<std::string> fname3d = Wildcard(
299  diskdirFD + "by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3d_nus_contain_v1_goodruns_prod2-snapshot/"
300  + "prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3d_nus_contain_v1_goodruns_prod2-snapshot*.root"
301  );
302  std::vector<std::string> ret;
303  for(const auto& file : fname1) { ret.push_back(file); }
304  for(const auto& file : fname2) { ret.push_back(file); }
305  for(const auto& file : fname3b) { ret.push_back(file); }
306  for(const auto& file : fname3c) { ret.push_back(file); }
307  for(const auto& file : fname3d) { ret.push_back(file); }
308  return ret;
309 }
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
float partptp
Definition: NusVarsTemp.cxx:50
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const double kSecondAnaPeriod2POT
Definition: Exposures.h:74
Sample(std::string shortName, std::string longName, Cut cut)
Definition: numu_cut_flow.C:98
const Var kCVNnc
PID
Definition: Vars.cxx:44
SRCosRej cosrej
Output from CosRej (Cosmic Rejection)
Definition: SRIDBranch.h:47
std::vector< std::string > MakeUnblindList()
Definition: Ana01Data.h:7
std::vector< SRFuzzyKProng > png
Vector of 3D prong objects.
Definition: SRFuzzyK.h:19
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2085
const double kSecondAnaEpoch3dPOT
Definition: Exposures.h:77
const double kSecondAnaEpoch3bPOT
Definition: Exposures.h:75
const Cut kNusFDAllE
Definition: NusCuts.h:45
const Cut kNusNDPresel
Definition: NusCuts.h:69
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 Cut kNusNCSel([](const caf::SRProxy *sr){if(sr->slc.nhit >=200) return false;if(sr->slc.nhit< 20) return false;if(sr->sel.cvn.ncid< 0.2) return false;return true;})
Cut that is more CC rejection than NC selection from docdb 15285.
Definition: NusCuts.h:27
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
const std::string fnamenear_concat
const Cut kNusNDAllE
Definition: NusCuts.h:70
std::vector< std::string > Wildcard(const std::string &wildcardString)
Find files matching a UNIX glob, plus expand environment variables.
Definition: UtilsExt.cxx:268
unsigned int nshwlid
number of shwlid showers - either 0 or number of 3d prongs
Definition: SRFuzzyK.h:24
const Cut kNCCalERange
Definition: NusCuts.h:35
const std::string fnameneardata_concat
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
const std::vector< std::string > fFDMC_swp
const Var kNHit
Definition: Vars.cxx:71
SRNueCosRej nuecosrej
Output from NueCosRej (Nue Cosmic Rejection)
Definition: SRIDBranch.h:48
const std::string fFDMC_non_3c
float calE
Calorimetric energy of the cluster [GeV].
Definition: SRSlice.h:38
const Cut kNumuND
Definition: NumuCuts.h:55
const Var kCCE
Definition: NumuVars.h:21
Sample MakeSample(const HistAxis &axnc, const HistAxis &axnm, const Cut &nd, const Cut &fd)
Definition: SideBandLoad.C:263
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
const std::string diskdirFD
caf::StandardRecord * sr
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
unsigned int nhit
number of hits
Definition: SRSlice.h:22
const char sep
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
::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
Splits Data proportionally according to MC.
const Var kTuftsWeightCC
Definition: XsecTunes.h:31
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
TDirectory * dir
Definition: macro.C:5
void SideBandLoad(std::string outfile)
Definition: SideBandLoad.C:39
const std::vector< std::string > fFDMC_tau
SRIDBranch sel
Selector (PID) branch.
SRElastic elastic
Single vertex found by Elastic Arms.
const std::string fFDMC_swp_3c
TFile * file
Definition: cellShifts.C:17
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
const Cut kNusFDPresel
Definition: NusCuts.h:38
const Cut kNusCosRejMod([](const caf::SRProxy *sr){double partptp=sr->sel.nuecosrej.partptp;if(sr->vtx.elastic.fuzzyk.nshwlid==0) return false;if(partptp >=0.8) return false;if(sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.ismuon==1) return false;double nhit=(double) sr->slc.nhit;if(nhit<=0.) return false;if(sr->slc.calE/nhit<=0.009) return false;return true;})
Definition: NusCuts.h:66
SRSlice slc
Slice branch: nhit, extents, time, etc.
A prediction object compatible with sterile oscillations.
const double kSecondAnaEpoch3cPOT
Definition: Exposures.h:76
SRFuzzyK fuzzyk
Primary 3D prong object.
Definition: SRElastic.h:44
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.
T min(const caf::Proxy< T > &a, T b)
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 Cut kNusCosRej([](const caf::SRProxy *sr){double numucontpid2019=sr->sel.cosrej.numucontpid2019;double partptp=sr->sel.nuecosrej.partptp;if(sr->vtx.elastic.fuzzyk.nshwlid==0) return false;if(numucontpid2019<=0.5) return false;if(partptp >=0.8) return false;if(sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.ismuon==1) return false;double nhit=(double) sr->slc.nhit;if(nhit<=0.) return false;if(sr->slc.calE/nhit<=0.018) return false;if(std::min(sr->sel.nuecosrej.starttop, sr->sel.nuecosrej.stoptop)< 480) return false;return true;})
Cosmic rejection for the NC sample from docdb 15241.
Definition: NusCuts.h:33
SRVertexBranch vtx
Vertex branch: location, time, etc.