nue_michelDataMC.C
Go to the documentation of this file.
1 #include "OscLib/OscCalcDumb.h"
2 
3 #include "CAFAna/Core/Spectrum.h"
5 #include "CAFAna/Systs/Systs.h"
7 #include "CAFAna/Core/MultiVar.h"
8 
9 #include "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/Cuts/SpillCuts.h"
11 #include "CAFAna/Cuts/NueCutsSecondAna.h"
14 
20 
22 
23 #include "CAFAna/Vars/Vars.h"
28 #include "CAFAna/Vars/HistAxes.h"
29 
30 
33 
34 #include "Utilities/func/MathUtil.h"
35 
36 #include "TCanvas.h"
37 #include "TH2.h"
38 #include "TH1.h"
39 #include "TFile.h"
40 #include "TLegend.h"
41 
42 #include <iostream>
43 #include <cmath>
44 
45 using namespace ana;
46 
47 
48 const MultiVar kMECalE([](const caf::SRProxy *sr){
49  std::vector<double> ret;
50  for (unsigned int i = 0; i < sr->me.nslc; i++){
51  if (sr->me.slc[i].deltat > 800)
52  ret.push_back(1000*sr->me.slc[i].calE);
53  }
54  for (unsigned int i = 0; i < sr->me.nkalman; i++){
55  if (sr->me.trkkalman[i].deltat > 800)
56  ret.push_back(1000*sr->me.trkkalman[i].calE);
57  }
58  return ret;
59  });
60 const MultiVar kMECalETrueME([](const caf::SRProxy *sr){
61  std::vector<double> ret;
62  for (unsigned int i = 0; i < sr->me.nslc; i++){
63  if (sr->me.slc[i].deltat > 800 &&
64  std::abs(sr->me.slc[i].truth.pdg) == 11)
65  ret.push_back(1000*sr->me.slc[i].calE);
66  }
67  for (unsigned int i = 0; i < sr->me.nkalman; i++){
68  if (sr->me.trkkalman[i].deltat > 800 &&
69  std::abs(sr->me.trkkalman[i].truth.pdg) == 11)
70  ret.push_back(1000*sr->me.trkkalman[i].calE);
71  }
72  return ret;
73  });
74 
75 
76 const MultiVar kMENHit([](const caf::SRProxy *sr){
77  std::vector<double> ret;
78  for (unsigned int i = 0; i < sr->me.nslc; i++){
79  if (sr->me.slc[i].deltat > 800)
80  ret.push_back(sr->me.slc[i].nhitx+sr->me.slc[i].nhity);
81  }
82  for (unsigned int i = 0; i < sr->me.nkalman; i++){
83  if (sr->me.trkkalman[i].deltat > 800)
84  ret.push_back(sr->me.trkkalman[i].nhitx+sr->me.trkkalman[i].nhity);
85  }
86  return ret;
87  });
88 const MultiVar kMENHitTrueME([](const caf::SRProxy *sr){
89  std::vector<double> ret;
90  for (unsigned int i = 0; i < sr->me.nslc; i++){
91  if (sr->me.slc[i].deltat > 800 &&
92  std::abs(sr->me.slc[i].truth.pdg) == 11)
93  ret.push_back(sr->me.slc[i].nhitx+sr->me.slc[i].nhity);
94  }
95  for (unsigned int i = 0; i < sr->me.nkalman; i++){
96  if (sr->me.trkkalman[i].deltat > 800 &&
97  std::abs(sr->me.trkkalman[i].truth.pdg) == 11)
98  ret.push_back(sr->me.trkkalman[i].nhitx+sr->me.trkkalman[i].nhity);
99  }
100  return ret;
101  });
102 
103 
104 const MultiVar kMEDeltaT([](const caf::SRProxy *sr){
105  std::vector<double> ret;
106  for (unsigned int i = 0; i < sr->me.nslc; i++){
107  if (sr->me.slc[i].deltat > 800)
108  ret.push_back(0.001*sr->me.slc[i].deltat);
109  }
110  for (unsigned int i = 0; i < sr->me.nkalman; i++){
111  if (sr->me.trkkalman[i].deltat > 800)
112  ret.push_back(0.001*sr->me.trkkalman[i].deltat);
113  }
114  return ret;
115  });
116 const MultiVar kMEDeltaTTrueME([](const caf::SRProxy *sr){
117  std::vector<double> ret;
118  for (unsigned int i = 0; i < sr->me.nslc; i++){
119  if (sr->me.slc[i].deltat > 800 &&
120  std::abs(sr->me.slc[i].truth.pdg) == 11)
121  ret.push_back(0.001*sr->me.slc[i].deltat);
122  }
123  for (unsigned int i = 0; i < sr->me.nkalman; i++){
124  if (sr->me.trkkalman[i].deltat > 800 &&
125  std::abs(sr->me.trkkalman[i].truth.pdg) == 11)
126  ret.push_back(0.001*sr->me.trkkalman[i].deltat);
127  }
128  return ret;
129  });
130 
131 
132 const MultiVar kMEDist([](const caf::SRProxy *sr){
133  std::vector<double> ret;
134  for (unsigned int i = 0; i < sr->me.nslc; i++){
135  if (sr->me.slc[i].deltat > 800)
136  ret.push_back(sr->me.slc[i].disttoslc);
137  }
138  for (unsigned int i = 0; i < sr->me.nkalman; i++){
139  if (sr->me.trkkalman[i].deltat > 800)
140  ret.push_back(sr->me.trkkalman[i].disttotrack);
141  }
142  return ret;
143  });
144 const MultiVar kMEDistTrueME([](const caf::SRProxy *sr){
145  std::vector<double> ret;
146  for (unsigned int i = 0; i < sr->me.nslc; i++){
147  if (sr->me.slc[i].deltat > 800 &&
148  std::abs(sr->me.slc[i].truth.pdg) == 11)
149  ret.push_back(sr->me.slc[i].disttoslc);
150  }
151  for (unsigned int i = 0; i < sr->me.nkalman; i++){
152  if (sr->me.trkkalman[i].deltat > 800 &&
153  std::abs(sr->me.trkkalman[i].truth.pdg) == 11)
154  ret.push_back(sr->me.trkkalman[i].disttotrack);
155  }
156  return ret;
157  });
158 
159 
160 const MultiVar kMEMID([](const caf::SRProxy *sr){
161  std::vector<double> ret;
162  for (unsigned int i = 0; i < sr->me.nslc; i++){
163  if (sr->me.slc[i].deltat > 800)
164  ret.push_back(sr->me.slc[i].mid);
165  }
166  for (unsigned int i = 0; i < sr->me.nkalman; i++){
167  if (sr->me.trkkalman[i].deltat > 800)
168  ret.push_back(sr->me.trkkalman[i].mid);
169  }
170  return ret;
171  });
172 const MultiVar kMEMIDTrueME([](const caf::SRProxy *sr){
173  std::vector<double> ret;
174  for (unsigned int i = 0; i < sr->me.nslc; i++){
175  if (sr->me.slc[i].deltat > 800 &&
176  std::abs(sr->me.slc[i].truth.pdg) == 11)
177  ret.push_back(sr->me.slc[i].mid);
178  }
179  for (unsigned int i = 0; i < sr->me.nkalman; i++){
180  if (sr->me.trkkalman[i].deltat > 800 &&
181  std::abs(sr->me.trkkalman[i].truth.pdg) == 11)
182  ret.push_back(sr->me.trkkalman[i].mid);
183  }
184  return ret;
185  });
186 
187 
188 void StandardSavePlots(TDirectory *dir,
189  std::vector<Spectrum*> spects, std::string name)
190 {
191  for (int sIdx = 0; sIdx < 3; sIdx++){
192  dir->cd();
193  std::string idxstr = "DA";
194  if (sIdx == 1) idxstr = "MC";
195  if (sIdx == 2) idxstr = "MCTrueME";
196 
197  spects[sIdx]->SaveTo(dir, (name+"_"+idxstr).c_str());
198  }
199 }
200 
201 
202 std::vector<Spectrum*> MakeSpectra(SpectrumLoader *lMC, SpectrumLoader *lDa,
204  MultiVar var, MultiVar varME,
205  Cut sel, Var wei)
206 {
207  std::vector<Spectrum*> ret;
208  ret.push_back(new Spectrum(label, bins, *lDa, var, sel, kNoShift, wei));
209  ret.push_back(new Spectrum(label, bins, *lMC, var, sel, kNoShift, wei));
210  ret.push_back(new Spectrum(label, bins, *lMC, varME, sel, kNoShift, wei));
211 
212  return ret;
213 }
214 
215 
216 
217 
218 void nue_michelDataMC(bool doNumu)
219 {
220  double pot = 9e20; // Pick something
221 
222  SpectrumLoader lNDMC("prod_decaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1");
223  SpectrumLoader lNDDa("prod_caf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns");
226 
227  // Configure weights / cuts to 2017 analysis
229  Cut sel = doNumu ? kNumuND : kNue2017NDPresel;
230 
231  std::vector<Spectrum*> calEs =
232  MakeSpectra(&lNDMC, &lNDDa,"Cal E [MeV]", Binning::Simple(100,0,100),
233  kMECalE, kMECalETrueME, sel, wei);
234  std::vector<Spectrum*> NHits =
235  MakeSpectra(&lNDMC, &lNDDa,"NHit", Binning::Simple(10,0,10),
236  kMENHit, kMENHitTrueME, sel, wei);
237  std::vector<Spectrum*> DeltaTs =
238  MakeSpectra(&lNDMC, &lNDDa,"DeltaT [#mus]", Binning::Simple(200,0,10),
239  kMEDeltaT, kMEDeltaTTrueME, sel, wei);
240  std::vector<Spectrum*> Dists =
241  MakeSpectra(&lNDMC, &lNDDa,"DistToSlc [cm]", Binning::Simple(20,0,20),
242  kMEDist, kMEDistTrueME, sel, wei);
243  std::vector<Spectrum*> MIDs =
244  MakeSpectra(&lNDMC, &lNDDa,"MID", Binning::Simple(200,-8,8),
245  kMEMID, kMEMIDTrueME, sel, wei);
246 
247 
248 
249  lNDMC.Go();
250  lNDDa.Go();
251 
252  std::string nStr = doNumu ? "_numu" : "";
253  TFile *outFile = new TFile(("out_MEDataMC"+nStr+".root").c_str(),"recreate");
254 
255  StandardSavePlots(outFile->mkdir("calE"),calEs,"calE");
256  StandardSavePlots(outFile->mkdir("NHit"),NHits,"NHit");
257  StandardSavePlots(outFile->mkdir("DeltaT"),DeltaTs,"DeltaT");
258  StandardSavePlots(outFile->mkdir("Dist"),Dists,"Dist");
259  StandardSavePlots(outFile->mkdir("MID"),MIDs,"MID");
260 
261  outFile->Close();
262 }
void StandardSavePlots(TDirectory *dir, std::vector< Spectrum * > spects, std::string name)
const XML_Char * name
Definition: expat.h:151
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
caf::Proxy< size_t > nslc
Definition: SRProxy.h:716
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const MultiVar kMEDeltaTTrueME([](const caf::SRProxy *sr){std::vector< double > ret;for(unsigned int i=0;i< sr->me.nslc;i++){if(sr->me.slc[i].deltat > 800 && std::abs(sr->me.slc[i].truth.pdg)==11) ret.push_back(0.001 *sr->me.slc[i].deltat);}for(unsigned int i=0;i< sr->me.nkalman;i++){if(sr->me.trkkalman[i].deltat > 800 && std::abs(sr->me.trkkalman[i].truth.pdg)==11) ret.push_back(0.001 *sr->me.trkkalman[i].deltat);}return ret;})
caf::Proxy< std::vector< caf::SRTrkME > > trkkalman
Definition: SRProxy.h:721
const MultiVar kMEDeltaT([](const caf::SRProxy *sr){std::vector< double > ret;for(unsigned int i=0;i< sr->me.nslc;i++){if(sr->me.slc[i].deltat > 800) ret.push_back(0.001 *sr->me.slc[i].deltat);}for(unsigned int i=0;i< sr->me.nkalman;i++){if(sr->me.trkkalman[i].deltat > 800) ret.push_back(0.001 *sr->me.trkkalman[i].deltat);}return ret;})
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
caf::Proxy< std::vector< caf::SRSlcME > > slc
Definition: SRProxy.h:717
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< caf::SRMichelE > me
Definition: SRProxy.h:2139
void SetSpillCut(const SpillCut &cut)
const MultiVar kMEMIDTrueME([](const caf::SRProxy *sr){std::vector< double > ret;for(unsigned int i=0;i< sr->me.nslc;i++){if(sr->me.slc[i].deltat > 800 && std::abs(sr->me.slc[i].truth.pdg)==11) ret.push_back(sr->me.slc[i].mid);}for(unsigned int i=0;i< sr->me.nkalman;i++){if(sr->me.trkkalman[i].deltat > 800 && std::abs(sr->me.trkkalman[i].truth.pdg)==11) ret.push_back(sr->me.trkkalman[i].mid);}return ret;})
float abs(float number)
Definition: d0nt_math.hpp:39
const char * label
caf::Proxy< size_t > nkalman
Definition: SRProxy.h:715
TFile * outFile
Definition: PlotXSec.C:135
const MultiVar kMENHit([](const caf::SRProxy *sr){std::vector< double > ret;for(unsigned int i=0;i< sr->me.nslc;i++){if(sr->me.slc[i].deltat > 800) ret.push_back(sr->me.slc[i].nhitx+sr->me.slc[i].nhity);}for(unsigned int i=0;i< sr->me.nkalman;i++){if(sr->me.trkkalman[i].deltat > 800) ret.push_back(sr->me.trkkalman[i].nhitx+sr->me.trkkalman[i].nhity);}return ret;})
const Cut kNumuND
Definition: NumuCuts.h:55
#define pot
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
std::vector< float > Spectrum
Definition: Constants.h:759
void nue_michelDataMC(bool doNumu)
const SystShifts kNoShift
Definition: SystShifts.cxx:22
const Cut kNue2017NDPresel
Definition: NueCuts2017.h:285
const Binning bins
Definition: NumuCC_CPiBin.h:8
const MultiVar kMENHitTrueME([](const caf::SRProxy *sr){std::vector< double > ret;for(unsigned int i=0;i< sr->me.nslc;i++){if(sr->me.slc[i].deltat > 800 && std::abs(sr->me.slc[i].truth.pdg)==11) ret.push_back(sr->me.slc[i].nhitx+sr->me.slc[i].nhity);}for(unsigned int i=0;i< sr->me.nkalman;i++){if(sr->me.trkkalman[i].deltat > 800 && std::abs(sr->me.trkkalman[i].truth.pdg)==11) ret.push_back(sr->me.trkkalman[i].nhitx+sr->me.trkkalman[i].nhity);}return ret;})
const MultiVar kMEMID([](const caf::SRProxy *sr){std::vector< double > ret;for(unsigned int i=0;i< sr->me.nslc;i++){if(sr->me.slc[i].deltat > 800) ret.push_back(sr->me.slc[i].mid);}for(unsigned int i=0;i< sr->me.nkalman;i++){if(sr->me.trkkalman[i].deltat > 800) ret.push_back(sr->me.trkkalman[i].mid);}return ret;})
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
TDirectory * dir
Definition: macro.C:5
const MultiVar kMEDist([](const caf::SRProxy *sr){std::vector< double > ret;for(unsigned int i=0;i< sr->me.nslc;i++){if(sr->me.slc[i].deltat > 800) ret.push_back(sr->me.slc[i].disttoslc);}for(unsigned int i=0;i< sr->me.nkalman;i++){if(sr->me.trkkalman[i].deltat > 800) ret.push_back(sr->me.trkkalman[i].disttotrack);}return ret;})
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const MultiVar kMEDistTrueME([](const caf::SRProxy *sr){std::vector< double > ret;for(unsigned int i=0;i< sr->me.nslc;i++){if(sr->me.slc[i].deltat > 800 && std::abs(sr->me.slc[i].truth.pdg)==11) ret.push_back(sr->me.slc[i].disttoslc);}for(unsigned int i=0;i< sr->me.nkalman;i++){if(sr->me.trkkalman[i].deltat > 800 && std::abs(sr->me.trkkalman[i].truth.pdg)==11) ret.push_back(sr->me.trkkalman[i].disttotrack);}return ret;})
const Var kXSecCVWgt2017
Definition: XsecTunes.h:36
std::vector< Spectrum * > MakeSpectra(SpectrumLoader *lMC, SpectrumLoader *lDa, std::string label, Binning bins, MultiVar var, MultiVar varME, Cut sel, Var wei)
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const MultiVar kMECalETrueME([](const caf::SRProxy *sr){std::vector< double > ret;for(unsigned int i=0;i< sr->me.nslc;i++){if(sr->me.slc[i].deltat > 800 && std::abs(sr->me.slc[i].truth.pdg)==11) ret.push_back(1000 *sr->me.slc[i].calE);}for(unsigned int i=0;i< sr->me.nkalman;i++){if(sr->me.trkkalman[i].deltat > 800 && std::abs(sr->me.trkkalman[i].truth.pdg)==11) ret.push_back(1000 *sr->me.trkkalman[i].calE);}return ret;})
const MultiVar kMECalE([](const caf::SRProxy *sr){std::vector< double > ret;for(unsigned int i=0;i< sr->me.nslc;i++){if(sr->me.slc[i].deltat > 800) ret.push_back(1000 *sr->me.slc[i].calE);}for(unsigned int i=0;i< sr->me.nkalman;i++){if(sr->me.trkkalman[i].deltat > 800) ret.push_back(1000 *sr->me.trkkalman[i].calE);}return ret;})
enum BeamMode string