Vars.cxx
Go to the documentation of this file.
1 #include "CAFAna/Vars/Vars.h"
2 
4 
5 #include "Utilities/func/MathUtil.h"
6 
8 
9 #include "TFile.h"
10 #include "TH1.h"
11 
12 #include <iostream>
13 #include <vector>
14 
15 namespace ana
16 {
17 
18  const Cut kIsRHC([](const caf::SRProxy* sr) {return sr->spill.isRHC;});
19 
20  const Var kRun = SIMPLEVAR(hdr.run);
21  const Var kSubrun = SIMPLEVAR(hdr.subrun);
22  const Var kCycle = SIMPLEVAR(hdr.cycle);
23  const Var kEvt = SIMPLEVAR(hdr.evt);
24  const Var kSlc = SIMPLEVAR(hdr.subevt);
25 
26  const Var kLEM = SIMPLEVAR(sel.lem.pid);
27  const Var kLID = SIMPLEVAR(sel.lid.ann);
28  const Var kRVP = SIMPLEVAR(sel.rvp.pid);
29 
30  const Var kCVNSSe([](const caf::SRProxy* sr){
31  throw std::runtime_error("kCVNSSe is no longer available. Fix your macro so you don't use it.");
32  return -5.;
33  });
34 
35  const Var kCVNe = SIMPLEVAR(sel.cvn.nueid);
36  const Var kCVNe_looseptp = SIMPLEVAR(sel.cvnloosepreselptp.nueid);
37  const Var kCVNe_oldpresel = SIMPLEVAR(sel.cvnoldpresel.nueid);
38 
39  const Var kCVNm = SIMPLEVAR(sel.cvn.numuid);
40  const Var kCVNm_looseptp = SIMPLEVAR(sel.cvnloosepreselptp.numuid);
41  const Var kCVNm_oldpresel = SIMPLEVAR(sel.cvnoldpresel.numuid);
42 
43  const Var kCVNt = SIMPLEVAR(sel.cvn.nutauid);
44  const Var kCVNnc = SIMPLEVAR(sel.cvn.ncid);
45 
46  //const Var kCVNcos = SIMPLEVAR(sel.cvn.output[390]);
47  const Var kCVNcos_looseptp = SIMPLEVAR(sel.cvnloosepreselptp.cosmicid);
48  const Var kCVNcos_oldpresel = SIMPLEVAR(sel.cvnoldpresel .cosmicid);
49 
50  const Var kCosVeto = SIMPLEVAR(sel.veto.keep);
51 
52  const Var kCaloE = SIMPLEVAR(slc.calE);
53 
54  const Var kMeanTime = SIMPLEVAR(slc.meantime);
55 
56  const Var kNMichels([](const caf::SRProxy* sr){
57  int n_me = 0;
58  // Add MID selected SlcMEs
59  for (int i = 0; i < (int)sr->me.nslc; i++)
60  if (sr->me.slc[i].mid > 1.
61  && sr->me.slc[i].deltat > 1200.)
62  n_me++;
63  // Add MID selected TrkMEs
64  for (int i = 0; i < (int)sr->me.nkalman; i++)
65  if (sr->me.trkkalman[i].mid > 0.
66  && sr->me.trkkalman[i].deltat > 800.)
67  n_me++;
68  if (n_me > 2) n_me = 2; // For 0, 1, 2+ ME Binning
69  return n_me;});
70 
71  const Var kNHit = SIMPLEVAR(slc.nhit);
72 
73  const Var kLongestProng([](const caf::SRProxy* sr)
74  {
75  if(!sr->vtx.elastic.IsValid) return 0.f;
76  if(sr->vtx.elastic.fuzzyk.npng == 0) return 0.f;
77  auto idx = sr->vtx.elastic.fuzzyk.longestidx;
78  return float(sr->vtx.elastic.fuzzyk.png[idx].len);
79  });
80 
81  const Var kRemID = SIMPLEVAR(sel.remid.pid);
82 
83  const NuTruthVar kTrueE_NT([](const caf::SRNeutrinoProxy* sr)
84  {return float(sr->E);});
86 
87  const Var kMode ([](const caf::SRProxy* sr)
88  {return (sr->mc.nnu == 0) ? -1 : int(sr->mc.nu[0].mode);});
89 
90  const Var kTruetpT([](const caf::SRProxy* sr)
91  {return (sr->mc.nnu == 0) ? 0. : util::pythag(sr->mc.nu[0].beam.tp.X(),sr->mc.nu[0].beam.tp.Y());});
92 
93  const Var kTruetpz([](const caf::SRProxy* sr)
94  {return (sr->mc.nnu == 0) ? 0. : float(sr->mc.nu[0].beam.tp.z);});
95 
96 
97  const Var kPartonZ([](const caf::SRProxy * sr)
98  {
99  double hadEvis=0.;
100  double energy=-5.;
101  if(sr->mc.nu[0].prim.size() > 0){
102  hadEvis = sr->energy.numu.hadcalE;
103  for( unsigned int part = 0; part < sr->mc.nu[0].prim.size(); part++)
104  {
105  if(sr->mc.nu[0].prim[part].pdg == 2212){
106  //calculate KE if proton/neutrons
107  double KE = sr->mc.nu[0].prim[part].p.E - 938;
108  if( energy < KE)
109  energy = KE;
110  }
111  else if(sr->mc.nu[0].prim[part].pdg == 2112){
112  //calculate KE if proton/neutrons
113  double KE = sr->mc.nu[0].prim[part].p.E - 939.5;
114  if( energy < KE)
115  energy = KE;
116  }
117  else{
118  //calculate E
119  if( energy < sr->mc.nu[0].prim[part].p.E)
120  energy = sr->mc.nu[0].prim[part].p.E;
121  }
122  }//end for loop
123  }
124  if (hadEvis==0)
125  return -5.;
126  return energy/hadEvis;
127  });
128 
129 
130 
131 
132 
133 
134 
135 
136 
137 
138 
139 
140 
142  {
143  if(hdr.det != caf::kFARDET) return 1; // No problem in ND
144  //FD:
145  if(hdr.ismc){ // MC
146  if(hdr.run < 20753){ // low gain
147  return 0.9949/0.9844;
148  }
149  else{ // high gain
150  return 1;
151  }
152  }
153  else{ // DATA
154  return 0.9949;
155  }
156  }
157 
159  {
161  if(sr->hdr.det!=1 || sr->mc.nnu==0) return 1;
162 
163  double te = sr->mc.nu[0].E;
164  float weight = 1;
165  if(sr->mc.nu[0].pdgorig == 14) { // weights taken from Leo ratio plots
166  if(te < 10){weight = fRatiosm[0]->Interpolate(te);}
167  else {weight = 1;}
168  }
169  else if(sr->mc.nu[0].pdgorig == -14) {
170  if(te < 10) {weight = fRatiosm[1]->Interpolate(te);}
171  else {weight = 1;}
172  }
173  else if(sr->mc.nu[0].pdgorig == 12) {
174  if(te < 10) {weight = fRatiosm[2]->Interpolate(te);}
175  else {weight = 1;}
176  }
177  else if(sr->mc.nu[0].pdgorig == -12) {
178  if(te < 10) {weight = fRatiosm[3]->Interpolate(te);}
179  else {weight = 1;}
180  }
181  else weight = 1;
182  //std::cout<<"energy: "<<te<<" weight "<<weight<<std::endl;
183  return weight;
184  }
185 
187  if(fRatiosm[0]) return;
188  std::string fname = FindCAFAnaDir()+"/data/beam/MinervaFluxWtfromLeo.root";
189  TFile fin(fname.c_str(), "read");
190  if(fin.IsZombie()) {
191  std::cerr << "Could not open " << fname << std::endl;
192  return;
193  }
194  fRatiosm[0] = (TH1D*)fin.Get("NDnumu")->Clone();
195  fRatiosm[0]->SetDirectory(0);
196  fRatiosm[1] = (TH1D*)fin.Get("NDanumu")->Clone();
197  fRatiosm[1]->SetDirectory(0);
198  fRatiosm[2] = (TH1D*)fin.Get("NDnue")->Clone();
199  fRatiosm[2]->SetDirectory(0);
200  fRatiosm[3] = (TH1D*)fin.Get("NDanue")->Clone();
201  fRatiosm[3]->SetDirectory(0);
202  }
203 
204  //Mipp+NA49 from NOvA
205 
206 
208  {
209  InitializeMippNA49NOvAND();
210  if(sr->hdr.det!=1 || sr->mc.nnu==0) return 1;
211  double te = sr->mc.nu[0].E;
212  float weight = 1;
213  if(sr->mc.nu[0].pdgorig == 14) { // weights taken from Mipp+NA49 NOvA Xuebing & Kuldeep ratio plots
214  if(te < 20){weight = fRatiosn[0]->Interpolate(te);}
215  else {weight = 1;}
216  }
217  else if(sr->mc.nu[0].pdgorig == -14) {
218  if(te < 20) {weight = fRatiosn[1]->Interpolate(te);}
219  else {weight = 1;}
220  }
221  else if(sr->mc.nu[0].pdgorig == 12) {
222  if(te < 20) {weight = fRatiosn[2]->Interpolate(te);}
223  else {weight = 1;}
224  }
225 
226  else if(sr->mc.nu[0].pdgorig == -12) {
227  if(te < 20) {weight = fRatiosn[3]->Interpolate(te);}
228  else {weight = 1;}
229  }
230  else weight = 1;
231  return weight;
232  }
233  //
235  if(fRatiosn[0]) return;
236  std::string fname = FindCAFAnaDir()+"/data/beam/NOvAMippNA49NDFD.root";
237  TFile fin(fname.c_str(), "read");
238  if(fin.IsZombie()) {
239  std::cerr << "Could not open " << fname << std::endl;
240  return;
241  }
242  fRatiosn[0] = (TH1D*)fin.Get("NDnumu")->Clone();
243  fRatiosn[0]->SetDirectory(0);
244  fRatiosn[1] = (TH1D*)fin.Get("NDanumu")->Clone();
245  fRatiosn[1]->SetDirectory(0);
246  fRatiosn[2] = (TH1D*)fin.Get("NDnue")->Clone();
247  fRatiosn[2]->SetDirectory(0);
248  fRatiosn[3] = (TH1D*)fin.Get("NDanue")->Clone();
249  fRatiosn[3]->SetDirectory(0);
250  }
251 //FD Mipp+NA49
252 //
254  {
255  InitializeMippNA49NOvAFD();
256  if(sr->hdr.det!=2 || sr->mc.nnu==0) return 1;
257 
258  double te = sr->mc.nu[0].E;
259  float weight = 1;
260 
261  if(sr->mc.nu[0].pdgorig == 14) { // weights taken from Mipp+NA49 NOvA Xuebing & Kuldeep ratio plots
262  if(te < 20){weight = fRatiosf[0]->Interpolate(te);}
263  else {weight = 1;}
264  }
265  else if(sr->mc.nu[0].pdgorig == -14) {
266  if(te < 20) {weight = fRatiosf[1]->Interpolate(te);}
267  else {weight = 1;}
268  }
269  else if(sr->mc.nu[0].pdgorig == 12) {
270  if(te < 20) {weight = fRatiosf[2]->Interpolate(te);}
271  else {weight = 1;}
272  }
273  else if(sr->mc.nu[0].pdgorig == -12) {
274  if(te < 20) {weight = fRatiosf[3]->Interpolate(te);}
275  else {weight = 1;}
276  }
277  else weight = 1;
278  return weight;
279  }
281  if(fRatiosf[0]) return;
282  std::string fname = FindCAFAnaDir()+"/data/beam/NOvAMippNA49NDFD.root";
283  TFile fin(fname.c_str(), "read");
284  if(fin.IsZombie()) {
285  std::cerr << "Could not open " << fname << std::endl;
286  return;
287  }
288  fRatiosf[0] = (TH1D*)fin.Get("FDnumu")->Clone();
289  fRatiosf[0]->SetDirectory(0);
290  fRatiosf[1] = (TH1D*)fin.Get("FDanumu")->Clone();
291  fRatiosf[1]->SetDirectory(0);
292  fRatiosf[2] = (TH1D*)fin.Get("FDnue")->Clone();
293  fRatiosf[2]->SetDirectory(0);
294  fRatiosf[3] = (TH1D*)fin.Get("FDanue")->Clone();
295  fRatiosf[3]->SetDirectory(0);
296  }
297 }
caf::Proxy< size_t > npng
Definition: SRProxy.h:2020
caf::Proxy< caf::SRSpill > spill
Definition: SRProxy.h:2125
void InitializeMinervaND() const
Definition: Vars.cxx:186
TString fin
Definition: Style.C:24
const Var kMode([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?-1:int(sr->mc.nu[0].mode);})
Neutrino interaction mode.
Definition: Vars.h:99
const Var kLEM
PID
Definition: Vars.cxx:26
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2041
Far Detector at Ash River.
Definition: SREnums.h:11
caf::Proxy< size_t > nslc
Definition: SRProxy.h:695
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const NuTruthVar kTrueE_NT([](const caf::SRNeutrinoProxy *sr){return float(sr->E);})
True neutrino energy.
Definition: Vars.h:95
double operator()(const caf::SRProxy *sr) const
Definition: Vars.cxx:158
caf::Proxy< std::vector< caf::SRTrkME > > trkkalman
Definition: SRProxy.h:700
const Var kSlc
Definition: Vars.cxx:24
const Var kCVNnc
PID
Definition: Vars.cxx:44
const Var kCVNe
PID
Definition: Vars.cxx:35
caf::Proxy< std::vector< caf::SRSlcME > > slc
Definition: SRProxy.h:696
const Var weight
double CalibrationBugCorrectionFactor(const caf::SRHeaderProxy &hdr)
See docdb 23597.
Definition: Vars.cxx:141
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2119
caf::Proxy< caf::SRNumuEnergy > numu
Definition: SRProxy.h:213
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2108
void InitializeMippNA49NOvAND() const
Definition: Vars.cxx:234
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:597
caf::Proxy< caf::SRMichelE > me
Definition: SRProxy.h:2121
void InitializeMippNA49NOvAFD() const
Definition: Vars.cxx:280
OStream cerr
Definition: OStream.cxx:7
caf::Proxy< short int > nnu
Definition: SRProxy.h:596
const Var kSubrun
Definition: Vars.cxx:21
caf::Proxy< caf::SREnergyBranch > energy
Definition: SRProxy.h:2118
const Var kPartonZ([](const caf::SRProxy *sr){double hadEvis=0.;double energy=-5.;if(sr->mc.nu[0].prim.size() > 0){hadEvis=sr->energy.numu.hadcalE;for(unsigned int part=0;part< sr->mc.nu[0].prim.size();part++){if(sr->mc.nu[0].prim[part].pdg==2212){double KE=sr->mc.nu[0].prim[part].p.E-938;if(energy< KE) energy=KE;}else if(sr->mc.nu[0].prim[part].pdg==2112){double KE=sr->mc.nu[0].prim[part].p.E-939.5;if(energy< KE) energy=KE;}else{if(energy< sr->mc.nu[0].prim[part].p.E) energy=sr->mc.nu[0].prim[part].p.E;}}}if(hadEvis==0) return-5.;return energy/hadEvis;})
Definition: Vars.h:108
std::string FindCAFAnaDir()
Definition: Utilities.cxx:208
const Var kTruetpT([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:util::pythag(sr->mc.nu[0].beam.tp.X(), sr->mc.nu[0].beam.tp.Y());})
True neutrino ancestor (off-target) pT.
Definition: Vars.h:102
const Var kLID
PID
Definition: Vars.cxx:27
const Var kCVNt
PID
Definition: Vars.cxx:43
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2100
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2025
const Var kMeanTime
Definition: Vars.cxx:54
caf::Proxy< size_t > nkalman
Definition: SRProxy.h:694
TString part[npart]
Definition: Style.C:32
const Var kCosVeto
The result of CosVeto.
Definition: Vars.cxx:50
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
const Var kTrueE([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].E);})
Definition: Vars.cxx:85
const Var kLongestProng([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return 0.f;if(sr->vtx.elastic.fuzzyk.npng==0) return 0.f;auto idx=sr->vtx.elastic.fuzzyk.longestidx;return float(sr->vtx.elastic.fuzzyk.png[idx].len);})
Definition: Vars.h:89
const Var kEvt
Definition: Vars.cxx:23
const Var kRemID
PID
Definition: Vars.cxx:81
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
const Var kTruetpz([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].beam.tp.z);})
True neutrino ancestor (off-target) pz.
Definition: Vars.h:105
const Var kNHit
Definition: Vars.cxx:71
double energy
Definition: plottest35.C:25
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
caf::StandardRecord * sr
caf::Proxy< float > hadcalE
Definition: SRProxy.h:168
const Var kCycle
Definition: Vars.cxx:22
const Var kCVNSSe([](const caf::SRProxy *sr){throw std::runtime_error("kCVNSSe is no longer available. Fix your macro so you don't use it.");return-5.;})
2018 nue PID
Definition: Vars.h:52
const Var kNMichels([](const caf::SRProxy *sr){int n_me=0;for(int i=0;i< (int) sr->me.nslc;i++) if(sr->me.slc[i].mid > 1. &&sr->me.slc[i].deltat > 1200.) n_me++;for(int i=0;i< (int) sr->me.nkalman;i++) if(sr->me.trkkalman[i].mid > 0. &&sr->me.trkkalman[i].deltat > 800.) n_me++;if(n_me > 2) n_me=2;return n_me;})
Number of SlcME&#39;s assoc with slice.
Definition: Vars.h:85
const Var kCVNm_looseptp
Definition: Vars.cxx:40
const Var kCVNcos_looseptp
Definition: Vars.cxx:47
caf::Proxy< unsigned int > longestidx
Definition: SRProxy.h:2019
double operator()(const caf::SRProxy *sr) const
Definition: Vars.cxx:207
const Var kCVNm_oldpresel
Definition: Vars.cxx:41
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2120
const Cut kIsRHC([](const caf::SRProxy *sr){return sr->spill.isRHC;})
Definition: Vars.h:16
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:7
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2040
TH1D * fRatiosm[4]
Definition: Vars.h:117
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
const Var kRVP
PID
Definition: Vars.cxx:28
const Var kCVNcos_oldpresel
Definition: Vars.cxx:48
caf::Proxy< bool > isRHC
Definition: SRProxy.h:1355
const Var kRun
Definition: Vars.cxx:20
const Var kCVNe_oldpresel
Definition: Vars.cxx:37
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2128
const Var kCVNm
PID
Definition: Vars.cxx:39
const Var kCVNe_looseptp
Definition: Vars.cxx:36
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:231
double operator()(const caf::SRProxy *sr) const
Definition: Vars.cxx:253