make_pi0_xcheck.h
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TH1.h"
3 #include "TF1.h"
4 #include "TFile.h"
5 #include "TLegend.h"
6 #include "TLegendEntry.h"
7 #include "TArrayD.h"
8 #include "TAxis.h"
9 #include "TLatex.h"
10 
11 #include "CAFAna/Core/Cut.h"
12 #include "CAFAna/Cuts/NueCutsFirstAna.h"
13 #include "CAFAna/Analysis/Plots.h"
14 #include "CAFAna/Analysis/Style.h"
16 #include "CAFAna/Core/Spectrum.h"
17 #include "CAFAna/Cuts/SpillCuts.h"
18 #include "CAFAna/Vars/Vars.h"
21 
23 
24 using namespace ana;
25 
26 const double deadscale = 0.8747;
27 
28 const Cut kNueVtxContain([](const caf::SRProxy* sr)
29  {
30  if( !sr->vtx.elastic.IsValid ) return false;
31 
32  if(fabs(sr->vtx.elastic.vtx.X()) > 180) return false;
33  if(fabs(sr->vtx.elastic.vtx.Y()) > 180) return false;
34  if(fabs(sr->vtx.elastic.vtx.Z()) < 50) return false;
35  if(fabs(sr->vtx.elastic.vtx.Z()) > 1000) return false;
36  return true;
37  });
38 
39 const Cut kNueVtxXWest([](const caf::SRProxy* sr)
40  {
41  if(sr->vtx.elastic.vtx.X() >= 0) return true;
42  return false;
43  });
44 
45 const Cut kNueVtxXEast([](const caf::SRProxy* sr)
46  {
47  if(sr->vtx.elastic.vtx.X() < 0) return true;
48  return false;
49  });
50 
51 const Cut kNueVtxYTop([](const caf::SRProxy* sr)
52  {
53  if(sr->vtx.elastic.vtx.Y() >= 0) return true;
54  return false;
55  });
56 
57 const Cut kNueVtxYBot([](const caf::SRProxy* sr)
58  {
59  if(sr->vtx.elastic.vtx.Y() < 0) return true;
60  return false;
61  });
62 
63 const Cut kNueShwContain([](const caf::SRProxy* sr)
64  {
65  if(sr->vtx.elastic.fuzzyk.nshwlid == 0) return false;
66 
67  if(fabs(sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.X()) > 180) return false;
68  if(fabs(sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.Y()) > 180) return false;
69  if(sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.Z() < 200) return false;
70  if(sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.Z() > 1200) return false;
71  return true;
72  });
73 
74 const Cut kRemidCut = SIMPLEVAR(sel.remid.pid) > 0 && SIMPLEVAR(sel.remid.pid) < 0.5;
75 
76 const Cut kdEdxCut = SIMPLEVAR(sand.nue.dedxpng1) < 0.003 && SIMPLEVAR(sand.nue.dedxpng2) < 0.003;
77 
78 const Cut kTwoProngs([](const caf::SRProxy* sr)
79  {
80  if( !sr->vtx.elastic.IsValid ) return false;
81  return sr->vtx.elastic.fuzzyk.npng == 2;
82  });
83 
84 const Cut kPhotonID([](const caf::SRProxy* sr)
85  {
86  if( !sr->vtx.elastic.IsValid ) return false;
87  if(sr->vtx.elastic.fuzzyk.npng != 2) return false;
88  return sr->vtx.elastic.fuzzyk.png[0].cvnpart.photonid > 0.75 &&
89  sr->vtx.elastic.fuzzyk.png[1].cvnpart.photonid > 0.75;
90  });
91 
92 const Cut kPhotonLow([](const caf::SRProxy* sr)
93  {
94  if( !sr->vtx.elastic.IsValid ) return false;
95  if(sr->vtx.elastic.fuzzyk.npng != 2) return false;
96  return sr->vtx.elastic.fuzzyk.png[0].cvnpart.photonid > 0.6 &&
97  sr->vtx.elastic.fuzzyk.png[1].cvnpart.photonid > 0.6;
98  });
99 
100 const Cut kPhotonHigh([](const caf::SRProxy* sr)
101  {
102  if( !sr->vtx.elastic.IsValid) return false;
103  if(sr->vtx.elastic.fuzzyk.npng != 2) return false;
104  return sr->vtx.elastic.fuzzyk.png[0].cvnpart.photonid > 0.9 &&
105  sr->vtx.elastic.fuzzyk.png[1].cvnpart.photonid > 0.9;
106  });
107 
108 const Cut kPhotonMax([](const caf::SRProxy* sr)
109  {
110  if( !sr->vtx.elastic.IsValid ) return false;
111  if(sr->vtx.elastic.fuzzyk.npng != 2) return false;
112  double pid1 = sr->vtx.elastic.fuzzyk.png[0].cvnpart.photonid;
113  double pid2 = sr->vtx.elastic.fuzzyk.png[1].cvnpart.photonid;
114  if(pid1 < sr->vtx.elastic.fuzzyk.png[0].cvnpart.muonid ||
115  pid1 < sr->vtx.elastic.fuzzyk.png[0].cvnpart.electronid ||
116  pid1 < sr->vtx.elastic.fuzzyk.png[0].cvnpart.neutronid ||
117  pid1 < sr->vtx.elastic.fuzzyk.png[0].cvnpart.pionid ||
118  pid1 < sr->vtx.elastic.fuzzyk.png[0].cvnpart.otherid)
119  return false;
120  if(pid2 < sr->vtx.elastic.fuzzyk.png[0].cvnpart.muonid ||
121  pid2 < sr->vtx.elastic.fuzzyk.png[0].cvnpart.electronid ||
122  pid2 < sr->vtx.elastic.fuzzyk.png[0].cvnpart.neutronid ||
123  pid2 < sr->vtx.elastic.fuzzyk.png[0].cvnpart.pionid ||
124  pid2 < sr->vtx.elastic.fuzzyk.png[0].cvnpart.otherid)
125  return false;
126  return true;
127  });
128 
129 const Cut kMissingPlanesCut([](const caf::SRProxy* sr)
130  {
131  if( !sr->vtx.elastic.IsValid ) return false;
132  for(auto& it: sr->vtx.elastic.fuzzyk.png)
133  if(it.maxplanegap <= 1) return false;
134  return true;
135  });
136 
137 const Cut kContigPlanesCut([](const caf::SRProxy* sr)
138  {
139  if( !sr->vtx.elastic.IsValid ) return false;
140  for(auto& it: sr->vtx.elastic.fuzzyk.png)
141  if(it.maxplanecont <= 4) return false;
142  return true;
143  });
144 
145 const Var kLongestProngVar = SIMPLEVAR(sel.rvp.longtr);
146 
147 const Var kNumHits = SIMPLEVAR(slc.nhit);
148 
149 const Var kNumPlanes([](const caf::SRProxy* sr)
150  {
151  return sr->slc.lastplane - sr->slc.firstplane + 1;
152  });
153 
154 const Var kVtxX([](const caf::SRProxy* sr)
155  {
156  return sr->vtx.elastic.vtx.x;
157  });
158 
159 const Var kVtxY([](const caf::SRProxy* sr)
160  {
161  return sr->vtx.elastic.vtx.y;
162  });
163 
164 const Var kVtxZ([](const caf::SRProxy* sr)
165  {
166  return sr->vtx.elastic.vtx.z;
167  });
168 
169 const Var kProng1E([](const caf::SRProxy* sr)
170  {
171  return deadscale*1000*sr->vtx.elastic.fuzzyk.png[0].calE;
172  });
173 
174 const Var kProng2E([](const caf::SRProxy* sr)
175  {
176  if(sr->vtx.elastic.fuzzyk.npng < 2) return -1.;
177  return (double)deadscale*1000*sr->vtx.elastic.fuzzyk.png[1].calE;
178  });
179 
180 const MultiVar kProngE([](const caf::SRProxy* sr){
181  std::vector<double> ret;
182  ret.push_back(kProng1E(sr));
183  ret.push_back(kProng2E(sr));
184  return ret;
185  });
186 
187 const Var kSumE([](const caf::SRProxy* sr)
188  {
189  if( !sr->vtx.elastic.IsValid ) return -1.;
190  if(sr->vtx.elastic.fuzzyk.npng != 2) return -1.;
191  const auto& p1 = sr->vtx.elastic.fuzzyk.png[0];
192  const auto& p2 = sr->vtx.elastic.fuzzyk.png[1];
193 
194  return (double)1000*deadscale*(p1.calE + p2.calE);
195  });
196 
197 const Var kPngAsymE([](const caf::SRProxy* sr)
198  {
199  if( !sr->vtx.elastic.IsValid ) return -1.f;
200  if(sr->vtx.elastic.fuzzyk.npng != 2) return -1.f;
201  const auto& p1 = sr->vtx.elastic.fuzzyk.png[0];
202  const auto& p2 = sr->vtx.elastic.fuzzyk.png[1];
203 
204  return abs(p1.calE - p2.calE)/(p1.calE + p2.calE);
205  });
206 
207 const Var kProng1L([](const caf::SRProxy* sr)
208  {
209  return sr->vtx.elastic.fuzzyk.png[0].len;
210  });
211 
212 const Var kProng2L([](const caf::SRProxy* sr)
213  {
214  if(sr->vtx.elastic.fuzzyk.npng < 2) return -1.;
215  return (double)sr->vtx.elastic.fuzzyk.png[1].len;
216  });
217 
218 const Var kProng1NHit([](const caf::SRProxy* sr)
219  {
220  return sr->vtx.elastic.fuzzyk.png[0].nhit;
221  });
222 
223 const Var kProng2NHit([](const caf::SRProxy* sr)
224  {
225  if(sr->vtx.elastic.fuzzyk.npng < 2) return -1.;
226  return (double)sr->vtx.elastic.fuzzyk.png[1].nhit;
227  });
228 
229 const MultiVar kProngNHit([](const caf::SRProxy* sr){
230  std::vector<double> ret;
231  ret.push_back(kProng1NHit(sr));
232  ret.push_back(kProng2NHit(sr));
233  return ret;
234  });
235 
236 const Var kProng1EHit([](const caf::SRProxy* sr)
237  {
238  return kProng1E(sr)/kProng1NHit(sr);
239  });
240 
241 const Var kProng2EHit([](const caf::SRProxy* sr)
242  {
243  return kProng2E(sr)/kProng2NHit(sr);
244  });
245 
246 const MultiVar kEHit([](const caf::SRProxy* sr){
247  std::vector<double> ret;
248  ret.push_back(kProng1EHit(sr));
249  ret.push_back(kProng2EHit(sr));
250  return ret;
251  });
252 
253 const Var kProng1NPlane([](const caf::SRProxy* sr)
254  {
255  return sr->vtx.elastic.fuzzyk.png[0].nplane;
256  });
257 
258 const Var kProng2NPlane([](const caf::SRProxy* sr)
259  {
260  if(sr->vtx.elastic.fuzzyk.npng < 2) return -1;
261  return (int)sr->vtx.elastic.fuzzyk.png[1].nplane;
262  });
263 
264 const Var kDotProd([](const caf::SRProxy* sr)
265  {
266  if( !sr->vtx.elastic.IsValid) return -5.;
267  if(sr->vtx.elastic.fuzzyk.npng != 2) return -5.;
268  return sr->vtx.elastic.fuzzyk.png[0].dir.Unit().Dot(sr->vtx.elastic.fuzzyk.png[1].dir.Unit());
269  });
270 
271 const Var kTotalP([](const caf::SRProxy* sr)
272  {
273  if( !sr->vtx.elastic.IsValid ) return -1.;
274  if(sr->vtx.elastic.fuzzyk.npng != 2) return -1.;
275 
276  double png1E = 1000*deadscale*sr->vtx.elastic.fuzzyk.png[0].calE;
277  double png2E = 1000*deadscale*sr->vtx.elastic.fuzzyk.png[1].calE;
278 
279  double px1 = png1E*sr->vtx.elastic.fuzzyk.png[0].dir.x;
280  double py1 = png1E*sr->vtx.elastic.fuzzyk.png[0].dir.y;
281  double pz1 = png1E*sr->vtx.elastic.fuzzyk.png[0].dir.z;
282  double px2 = png2E*sr->vtx.elastic.fuzzyk.png[1].dir.x;
283  double py2 = png2E*sr->vtx.elastic.fuzzyk.png[1].dir.y;
284  double pz2 = png2E*sr->vtx.elastic.fuzzyk.png[1].dir.z;
285 
286  return sqrt((px1+px2)*(px1+px2)+(py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2));
287  });
288 
289 const Var kCosTotal([](const caf::SRProxy* sr)
290  {
291  if( !sr->vtx.elastic.IsValid ) return -5.;
292  if(sr->vtx.elastic.fuzzyk.npng != 2) return -5.;
293 
294  double x1 = sr->vtx.elastic.fuzzyk.png[0].dir.x;
295  double y1 = sr->vtx.elastic.fuzzyk.png[0].dir.y;
296  double z1 = sr->vtx.elastic.fuzzyk.png[0].dir.z;
297  double x2 = sr->vtx.elastic.fuzzyk.png[1].dir.x;
298  double y2 = sr->vtx.elastic.fuzzyk.png[1].dir.y;
299  double z2 = sr->vtx.elastic.fuzzyk.png[1].dir.z;
300 
301  TVector3 tot(x1+x2,y1+y2,z1+z2);
302 
303  return tot.Unit().Dot(beamDirND);
304  });
305 
306 const Var kDotProdTruth([](const caf::SRProxy* sr)
307  {
308  if( !sr->vtx.elastic.IsValid ) return -1.;
309  if(sr->vtx.elastic.fuzzyk.npng != 2) return -1.;
310  TVector3 p1(sr->vtx.elastic.fuzzyk.png[0].truth.p.px,sr->vtx.elastic.fuzzyk.png[0].truth.p.py,sr->vtx.elastic.fuzzyk.png[0].truth.p.pz);
311  TVector3 p2(sr->vtx.elastic.fuzzyk.png[1].truth.p.px,sr->vtx.elastic.fuzzyk.png[1].truth.p.py,sr->vtx.elastic.fuzzyk.png[1].truth.p.pz);
312  return p1.Unit().Dot(p2.Unit());
313  });
314 
315 const Var thetaRes([](const caf::SRProxy* sr)
316  {
317  if( !sr->vtx.elastic.IsValid ) return -1.;
318  if(sr->vtx.elastic.fuzzyk.npng != 2) return -1.;
319  double a1 = sr->vtx.elastic.fuzzyk.png[0].dir.Unit().Dot(sr->vtx.elastic.fuzzyk.png[1].dir.Unit());
320 
321  TVector3 p1(sr->vtx.elastic.fuzzyk.png[0].truth.p.px,sr->vtx.elastic.fuzzyk.png[0].truth.p.py,sr->vtx.elastic.fuzzyk.png[0].truth.p.pz);
322  TVector3 p2(sr->vtx.elastic.fuzzyk.png[1].truth.p.px,sr->vtx.elastic.fuzzyk.png[1].truth.p.py,sr->vtx.elastic.fuzzyk.png[1].truth.p.pz);
323  double a2 = p1.Unit().Dot(p2.Unit());
324  return (a1-a2)/a2;
325  });
326 
327 const Var kNPng([](const caf::SRProxy* sr)
328  {
329  if( !sr->vtx.elastic.IsValid ) return -1;
330 
331  return (int)sr->vtx.elastic.fuzzyk.npng;
332  });
333 
334 const Var kEHitSlc([](const caf::SRProxy* sr)
335  {
336  return 1000*deadscale*sr->slc.calE/sr->slc.nhit;
337  });
338 
339 const Var kESlc([](const caf::SRProxy* sr)
340  {
341  return 1000*deadscale*sr->slc.calE;
342  });
343 
344 const Var kHitSlc([](const caf::SRProxy* sr)
345  {
346  return sr->slc.nhit;
347  });
348 
349 const Var kMass([](const caf::SRProxy* sr)
350  {
351  if( !sr->vtx.elastic.IsValid ) return -1.;
352  if(sr->vtx.elastic.fuzzyk.npng != 2) return -1.;
353  const auto& p1 = sr->vtx.elastic.fuzzyk.png[0];
354  const auto& p2 = sr->vtx.elastic.fuzzyk.png[1];
355 
356  const double dot = p1.dir.Unit().Dot(p2.dir.Unit());
357  return deadscale*1000*sqrt(2*p1.calE*p2.calE*(1-dot));
358  });
359 
360 const Var kMassNoScale([](const caf::SRProxy* sr) //Visible Energy
361  {
362  if(!sr->vtx.elastic.IsValid) return -1.;
363  if(sr->vtx.elastic.fuzzyk.npng != 2) return -1.;
364  const auto& p1 = sr->vtx.elastic.fuzzyk.png[0];
365  const auto& p2 = sr->vtx.elastic.fuzzyk.png[1];
366 
367  const double dot = p1.dir.Unit().Dot(p2.dir.Unit());
368  return 1000*sqrt(2*p1.calE*p2.calE*(1-dot))/1.78;
369  });
370 
371 const Cut kTruePi0 = SIMPLEVAR(sand.nue.npi0) > 0;
372 
373 const Var kPng1Pur([](const caf::SRProxy* sr)
374  {
375  return sr->vtx.elastic.fuzzyk.png[0].truth.pur;
376  });
377 
378 const Var kPng2Pur([](const caf::SRProxy* sr)
379  {
380  if(sr->vtx.elastic.fuzzyk.npng < 2) return -1.;
381  return (double)sr->vtx.elastic.fuzzyk.png[1].truth.pur;
382  });
383 
384 const MultiVar kProngPur([](const caf::SRProxy* sr){
385  std::vector<double> ret;
386  ret.push_back(kPng1Pur(sr));
387  ret.push_back(kPng2Pur(sr));
388  return ret;
389  });
390 
391 const Var kPng1Eff([](const caf::SRProxy* sr)
392  {
393  return sr->vtx.elastic.fuzzyk.png[0].truth.eff;
394  });
395 
396 const Var kPng2Eff([](const caf::SRProxy* sr)
397  {
398  if(sr->vtx.elastic.fuzzyk.npng < 2) return -1.0;
399  return (double)sr->vtx.elastic.fuzzyk.png[1].truth.eff;
400  });
401 
402 const MultiVar kProngEff([](const caf::SRProxy* sr){
403  std::vector<double> ret;
404  ret.push_back(kPng1Eff(sr));
405  ret.push_back(kPng2Eff(sr));
406  return ret;
407  });
408 
409 const Var kPng1CVN([](const caf::SRProxy* sr)
410  {
411  return sr->vtx.elastic.fuzzyk.png[0].cvnpart.photonid;
412  });
413 
414 const Var kPng2CVN([](const caf::SRProxy* sr)
415  {
416  if(sr->vtx.elastic.fuzzyk.npng < 2) return -1.0;
417  return (double)sr->vtx.elastic.fuzzyk.png[1].cvnpart.photonid;
418  });
419 
420 const MultiVar kProngCVN([](const caf::SRProxy* sr){
421  std::vector<double> ret;
422  ret.push_back(kPng1CVN(sr));
423  ret.push_back(kPng2CVN(sr));
424  return ret;
425  });
426 
427 const Cut kTwoPhoton([](const caf::SRProxy* sr)
428  {
429  if( !sr->vtx.elastic.IsValid ) return false;
430  int count=0;
431  for(auto& it: sr->vtx.elastic.fuzzyk.png)
432  if(it.cvnpart.photonid > 0.75) ++count;
433  if(count < 2) return false;
434  return true;
435  });
436 
437 const Var kOne([](const caf::SRProxy* sr)
438  {
439  return 1;
440  });
441 
442 struct PlotDef
443 {
444  std::string suffix = "";
446  Binning bins = Binning::Simple(50, 0, 500);
448  bool mcOnly = 0;
449 };
450 
452 {
453  std::string suffix = "";
457  bool mcOnly = 0;
458 };
459 
460 struct PlotDef2D
461 {
462  std::string suffix = "";
463  std::string label1 = "";
464  Binning bins1 = Binning::Simple(50,0,500);
465  Var var1 = kOne;
466  std::string label2 = "";
467  Binning bins2 = Binning::Simple(50,0,500);
468  Var var2 = kOne;
469 };
470 
471 struct SelDef
472 {
473  std::string suffix = "";
475  const Cut cut = kNoCut;
476 };
477 
480 const Cut SAsel = qual_base && kdEdxCut && kRemidCut;
481 const Cut CVN_two = qual_base && kPhotonID;
482 const Cut CVN_dedx = CVN_two && kdEdxCut;
483 const Cut CVN_flat = CVN_two && kCosTotal > 0.75 && kTotalP > 200 && kTotalP < 1000;
484 const Cut CVN_low = qual_base && kPhotonLow;
485 const Cut CVN_high = qual_base && kPhotonHigh;
486 const Cut CVN_max = qual_base && kPhotonMax;
487 const Cut CVNASA = CVN_two && SAsel;
488 const Cut CVNNSA = CVN_two && !SAsel;
489 const Cut NCVNSA = !CVN_two && SAsel;
490 
491 const unsigned int kNSels = 11;
492 const std::vector<SelDef> sels =
493  {{"qual_base","Base Selection Criteria",qual_base},
494  {"SA_sel", "SA Selection", SAsel},
495  {"CVN_two", "Two Prong CVN Selection", CVN_two},
496  {"CVN_dedx", "dedx Selection", CVN_dedx},
497  {"CVN_flat", "Flat Kinematics Selection", CVN_flat},
498  {"CVN_low", "Two Prong CVN Selection", CVN_low},
499  {"CVN_high", "Two Prong CVN Selection", CVN_high},
500  {"CVN_max", "Two Prong CVN Selection", CVN_max},
501  {"CVN_SA", "CVN and SA selection", CVNASA},
502  {"CVN_NSA", "CVN and not SA selection", CVNNSA},
503  {"NCVN_SA", "Not CVN and SA selection", NCVNSA}
504  };
505 
506 const unsigned int kNPlots = 21;
507 const std::vector<PlotDef> plots =
508  {{"counting", "", Binning::Simple(3, 0, 3), kOne,0},
509  {"mass", "M_{#gamma#gamma} (MeV)", Binning::Simple(50, 0, 500), kMass,0},
510  {"asymE", "|E1-E2|/E1+E2", Binning::Simple(50, 0, 1), kPngAsymE,0},
511  {"Ptot", "P_{#gamma#gamma} (MeV)", Binning::Simple(50, 0, 4000), kTotalP,0},
512  {"Etot", "E_{#gamma#gamma} (MeV)", Binning::Simple(50, 0, 4000), kSumE,0},
513  {"Costot", "cos(#theta)", Binning::Simple(50, 0, 1), kCosTotal,0},
514  {"Open", "Cos(#theta_{#gamma#gamma})", Binning::Simple(50, 0, 1), kDotProd,0},
515  {"PngE1", "E_{p1} (MeV)", Binning::Simple(50, 0, 2000), kProng1E,0},
516  {"PngE2", "E_{p2} (MeV)", Binning::Simple(50, 0, 2000), kProng2E,0},
517  {"PngHit1", "hit_{p1}", Binning::Simple(50, 0, 100), kProng1NHit,0},
518  {"PngHit2", "hit_{p2}", Binning::Simple(50, 0, 100), kProng2NHit,0},
519  {"SlcEHit", "E_{slc}/hit_{slc} (MeV)", Binning::Simple(50, 0, 50), kEHitSlc,0},
520  {"SlcE", "E_{slc} (MeV)", Binning::Simple(50, 0, 4000), kESlc,0},
521  {"SlcHit", "hit_{slc}", Binning::Simple(50, 0, 200), kHitSlc,0},
522  {"PngEHit1", "E_{p1}/hit_{p1} (MeV)", Binning::Simple(50, 0, 50), kProng1EHit,0},
523  {"PngEHit2", "E_{p2}/hit_{p2} (MeV)", Binning::Simple(50, 0, 50), kProng2EHit,0},
524  {"PngPur1", "Prong 1 Purity", Binning::Simple(50, 0, 1), kPng1Pur,1},
525  {"PngPur2", "Prong 2 Purity", Binning::Simple(50, 0, 1), kPng2Pur,1},
526  {"vtxx", "x (cm)", Binning::Simple(50, -200, 200), kVtxX,0},
527  {"vtxy", "y (cm)", Binning::Simple(50, -200, 200), kVtxY,0},
528  {"vtxz", "z (cm)", Binning::Simple(50, 0, 1600), kVtxZ,0}
529  };
530 
531 const unsigned int kN2D = 11;
532 std::vector<PlotDef2D> plots2d =
533  {{"vertices_XY", "x (cm)", Binning::Simple(50,-200,200), kVtxX, "y (cm)", Binning::Simple(50,-200,200), kVtxY},
534  {"vertices_ZX", "z (cm)", Binning::Simple(50,0,1600), kVtxZ, "x (cm)", Binning::Simple(50,-200,200), kVtxX},
535  {"vertices_ZY", "z (cm)", Binning::Simple(50,0,1600), kVtxZ, "y (cm)", Binning::Simple(50,-200,200), kVtxY},
536  {"peak_Ptot", "P_{#gamma#gamma} (MeV)", Binning::Simple(50,0,4000), kTotalP, "M_{#gamma#gamma} (MeV)",Binning::Simple(50,0,500), kMass},
537  {"peak_Etot", "E_{#gamma#gamma} (MeV)", Binning::Simple(50,0,4000), kSumE, "M_{#gamma#gamma} (MeV)", Binning::Simple(50,0,500), kMass},
538  {"peak_asymE", "|E1-E2|/E1+E2", Binning::Simple(50,0,1), kPngAsymE, "M_{#gamma#gamma} (MeV)", Binning::Simple(50,0,500), kMass},
539  {"peak_Costot", "cos(#theta)", Binning::Simple(50,0,1), kCosTotal, "M_{#gamma#gamma} (MeV)", Binning::Simple(50,0,500), kMass},
540  {"peak_open", "cos(#theta_{#gamma#gamma})", Binning::Simple(50,0,1), kDotProd, "M_{#gamma#gamma} (MeV)", Binning::Simple(50,0,500), kMass},
541  {"peak_x", "x (cm)", Binning::Simple(50,-200,200), kVtxX, "M_{#gamma#gamma} (MeV)", Binning::Simple(50,0,500), kMass},
542  {"peak_y", "y (cm)", Binning::Simple(50,-200,200), kVtxY, "M_{#gamma#gamma} (MeV)", Binning::Simple(50,0,500), kMass},
543  {"peak_z", "z (cm)", Binning::Simple(50,0,1600), kVtxZ, "M_{#gamma#gamma} (MeV)", Binning::Simple(50,0,500), kMass}
544  };
545 
546 const unsigned int kNMulti = 6;
547 const std::vector<PlotDefMulti> plotsMulti =
548  {{"PngPur", "Prong Purity", Binning::Simple(50, 0, 1), kProngPur,1},
549  {"PngEff", "Prong Efficiency", Binning::Simple(50, 0, 1), kProngEff,1},
550  {"PngE", "E_{png} (MeV)", Binning::Simple(50, 0, 2000), kProngE,0},
551  {"PngHits", "hit_{png}", Binning::Simple(50, 0, 100), kProngNHit,0},
552  {"PngEHit", "E_{png}/hit_{png} (MeV)", Binning::Simple(50, 0, 50), kEHit,0},
553  {"PngCVN", "CVN_{png}", Binning::Simple(50, 0.5, 1), kProngCVN,0}
554  };
const Cut kMissingPlanesCut([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;for(auto &it:sr->vtx.elastic.fuzzyk.png) if(it.maxplanegap<=1) return false;return true;})
caf::Proxy< size_t > npng
Definition: SRProxy.h:2038
caf::Proxy< unsigned int > nshwlid
Definition: SRProxy.h:2040
std::vector< PlotDef2D > plots2d
const Var kProng2EHit([](const caf::SRProxy *sr){return kProng2E(sr)/kProng2NHit(sr);})
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
const Cut CVN_low
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
const Cut CVN_two
TH1F * a2
Definition: f2_nu.C:545
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const Cut kPhotonHigh([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.npng!=2) return false;return sr->vtx.elastic.fuzzyk.png[0].cvnpart.photonid > 0.9 &&sr->vtx.elastic.fuzzyk.png[1].cvnpart.photonid > 0.9;})
Float_t y1[n_points_granero]
Definition: compare.C:5
const Cut NCVNSA
const unsigned int kN2D
Float_t x1[n_points_granero]
Definition: compare.C:5
const Cut kNueShwContain([](const caf::SRProxy *sr){if(sr->vtx.elastic.fuzzyk.nshwlid==0) return false;if(fabs(sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.X()) > 180) return false;if(fabs(sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.Y()) > 180) return false;if(sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.Z()< 200) return false;if(sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.Z() > 1200) return false;return true;})
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
T sqrt(T number)
Definition: d0nt_math.hpp:156
const Var kPngAsymE([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.f;if(sr->vtx.elastic.fuzzyk.npng!=2) return-1.f;const auto &p1=sr->vtx.elastic.fuzzyk.png[0];const auto &p2=sr->vtx.elastic.fuzzyk.png[1];return abs(p1.calE-p2.calE)/(p1.calE+p2.calE);})
const Cut kContigPlanesCut([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;for(auto &it:sr->vtx.elastic.fuzzyk.png) if(it.maxplanecont<=4) return false;return true;})
const Var kPng1Pur([](const caf::SRProxy *sr){return sr->vtx.elastic.fuzzyk.png[0].truth.pur;})
const Var kPng2Eff([](const caf::SRProxy *sr){if(sr->vtx.elastic.fuzzyk.npng< 2) return-1.0;return(double) sr->vtx.elastic.fuzzyk.png[1].truth.eff;})
const Cut kPhotonMax([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.npng!=2) return false;double pid1=sr->vtx.elastic.fuzzyk.png[0].cvnpart.photonid;double pid2=sr->vtx.elastic.fuzzyk.png[1].cvnpart.photonid;if(pid1< sr->vtx.elastic.fuzzyk.png[0].cvnpart.muonid||pid1< sr->vtx.elastic.fuzzyk.png[0].cvnpart.electronid||pid1< sr->vtx.elastic.fuzzyk.png[0].cvnpart.neutronid||pid1< sr->vtx.elastic.fuzzyk.png[0].cvnpart.pionid||pid1< sr->vtx.elastic.fuzzyk.png[0].cvnpart.otherid) return false;if(pid2< sr->vtx.elastic.fuzzyk.png[0].cvnpart.muonid||pid2< sr->vtx.elastic.fuzzyk.png[0].cvnpart.electronid||pid2< sr->vtx.elastic.fuzzyk.png[0].cvnpart.neutronid||pid2< sr->vtx.elastic.fuzzyk.png[0].cvnpart.pionid||pid2< sr->vtx.elastic.fuzzyk.png[0].cvnpart.otherid) return false;return true;})
const unsigned int kNSels
const Var kPng1Eff([](const caf::SRProxy *sr){return sr->vtx.elastic.fuzzyk.png[0].truth.eff;})
const MultiVar kProngE([](const caf::SRProxy *sr){std::vector< double > ret;ret.push_back(kProng1E(sr));ret.push_back(kProng2E(sr));return ret;})
void abs(TH1 *hist)
const unsigned int kNPlots
const Var thetaRes([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.;if(sr->vtx.elastic.fuzzyk.npng!=2) return-1.;double a1=sr->vtx.elastic.fuzzyk.png[0].dir.Unit().Dot(sr->vtx.elastic.fuzzyk.png[1].dir.Unit());TVector3 p1(sr->vtx.elastic.fuzzyk.png[0].truth.p.px, sr->vtx.elastic.fuzzyk.png[0].truth.p.py, sr->vtx.elastic.fuzzyk.png[0].truth.p.pz);TVector3 p2(sr->vtx.elastic.fuzzyk.png[1].truth.p.px, sr->vtx.elastic.fuzzyk.png[1].truth.p.py, sr->vtx.elastic.fuzzyk.png[1].truth.p.pz);double a2=p1.Unit().Dot(p2.Unit());return(a1-a2)/a2;})
const MultiVar kProngPur([](const caf::SRProxy *sr){std::vector< double > ret;ret.push_back(kPng1Pur(sr));ret.push_back(kPng2Pur(sr));return ret;})
const Cut kPhotonID([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.npng!=2) return false;return sr->vtx.elastic.fuzzyk.png[0].cvnpart.photonid > 0.75 &&sr->vtx.elastic.fuzzyk.png[1].cvnpart.photonid > 0.75;})
const Var kVtxX
const Var kProng2E([](const caf::SRProxy *sr){if(sr->vtx.elastic.fuzzyk.npng< 2) return-1.;return(double) deadscale *1000 *sr->vtx.elastic.fuzzyk.png[1].calE;})
const MultiVar kProngNHit([](const caf::SRProxy *sr){std::vector< double > ret;ret.push_back(kProng1NHit(sr));ret.push_back(kProng2NHit(sr));return ret;})
const Var kLongestProngVar
const Cut kTwoProngs([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;return sr->vtx.elastic.fuzzyk.npng==2;})
const Cut qual_base
const char * label
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2043
caf::Proxy< unsigned int > nhit
Definition: SRProxy.h:1315
const Var kDotProdTruth([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.;if(sr->vtx.elastic.fuzzyk.npng!=2) return-1.;TVector3 p1(sr->vtx.elastic.fuzzyk.png[0].truth.p.px, sr->vtx.elastic.fuzzyk.png[0].truth.p.py, sr->vtx.elastic.fuzzyk.png[0].truth.p.pz);TVector3 p2(sr->vtx.elastic.fuzzyk.png[1].truth.p.px, sr->vtx.elastic.fuzzyk.png[1].truth.p.py, sr->vtx.elastic.fuzzyk.png[1].truth.p.pz);return p1.Unit().Dot(p2.Unit());})
const Var kPng2Pur([](const caf::SRProxy *sr){if(sr->vtx.elastic.fuzzyk.npng< 2) return-1.;return(double) sr->vtx.elastic.fuzzyk.png[1].truth.pur;})
caf::Proxy< float > z
Definition: SRProxy.h:108
const MultiVar kProngEff([](const caf::SRProxy *sr){std::vector< double > ret;ret.push_back(kPng1Eff(sr));ret.push_back(kPng2Eff(sr));return ret;})
const Var kProng1L([](const caf::SRProxy *sr){return sr->vtx.elastic.fuzzyk.png[0].len;})
const Cut sels[kNumSels]
Definition: vars.h:44
caf::Proxy< float > x
Definition: SRProxy.h:106
TH1F * a1
Definition: f2_nu.C:476
const Var kProng2NHit([](const caf::SRProxy *sr){if(sr->vtx.elastic.fuzzyk.npng< 2) return-1.;return(double) sr->vtx.elastic.fuzzyk.png[1].nhit;})
const Var kProng1EHit([](const caf::SRProxy *sr){return kProng1E(sr)/kProng1NHit(sr);})
const Var kProng1NPlane([](const caf::SRProxy *sr){return sr->vtx.elastic.fuzzyk.png[0].nplane;})
const Var kEHitSlc([](const caf::SRProxy *sr){return 1000 *deadscale *sr->slc.calE/sr->slc.nhit;})
const Cut kTwoPhoton([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;int count=0;for(auto &it:sr->vtx.elastic.fuzzyk.png) if(it.cvnpart.photonid > 0.75)++count;if(count< 2) return false;return true;})
const Cut kNueVtxYTop([](const caf::SRProxy *sr){if(sr->vtx.elastic.vtx.Y() >=0) return true;return false;})
caf::StandardRecord * sr
const Var kProng2L([](const caf::SRProxy *sr){if(sr->vtx.elastic.fuzzyk.npng< 2) return-1.;return(double) sr->vtx.elastic.fuzzyk.png[1].len;})
caf::Proxy< unsigned int > lastplane
Definition: SRProxy.h:1309
const Var kOne([](const caf::SRProxy *sr){return 1;})
const Cut SAsel
const Cut CVN_flat
const Var kPng2CVN([](const caf::SRProxy *sr){if(sr->vtx.elastic.fuzzyk.npng< 2) return-1.0;return(double) sr->vtx.elastic.fuzzyk.png[1].cvnpart.photonid;})
const Cut kTruePi0([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;if(sr->vtx.elastic.fuzzyk.nshwlid==0) return false;if(sr->vtx.elastic.fuzzyk.npng==0) return false;if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);if(sr->mc.nu[0].prim.size()==0) return false;float Pi0Mass=0.134977;float KE=-10.0;float TE=-99.0;int nprim=sr->mc.nu[0].prim.size();int npi0=0;for(int i=0;i< nprim;i++){if(sr->mc.nu[0].prim[i].pdg==111){TE=sr->mc.nu[0].prim[i].p.E;KE=pow(pow(TE, 2.0)-pow(Pi0Mass, 2.0), 0.5);if(KE >=0.0) npi0++;}}if(npi0){return true;}return false;})
const Cut kNueVtxXWest([](const caf::SRProxy *sr){if(sr->vtx.elastic.vtx.X() >=0) return true;return false;})
const Cut qual
const Var kHitSlc([](const caf::SRProxy *sr){return sr->slc.nhit;})
const Var kTotalP([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.;if(sr->vtx.elastic.fuzzyk.npng!=2) return-1.;double png1E=1000 *deadscale *sr->vtx.elastic.fuzzyk.png[0].calE;double png2E=1000 *deadscale *sr->vtx.elastic.fuzzyk.png[1].calE;double px1=png1E *sr->vtx.elastic.fuzzyk.png[0].dir.x;double py1=png1E *sr->vtx.elastic.fuzzyk.png[0].dir.y;double pz1=png1E *sr->vtx.elastic.fuzzyk.png[0].dir.z;double px2=png2E *sr->vtx.elastic.fuzzyk.png[1].dir.x;double py2=png2E *sr->vtx.elastic.fuzzyk.png[1].dir.y;double pz2=png2E *sr->vtx.elastic.fuzzyk.png[1].dir.z;return sqrt((px1+px2)*(px1+px2)+(py1+py2)*(py1+py2)+(pz1+pz2)*(pz1+pz2));})
const std::vector< PlotDefMulti > plotsMulti
const TVector3 beamDirND
const std::vector< PlotDef > plots
const Cut CVN_dedx
const Binning bins
Definition: NumuCC_CPiBin.h:8
const Var kProng2NPlane([](const caf::SRProxy *sr){if(sr->vtx.elastic.fuzzyk.npng< 2) return-1;return(int) sr->vtx.elastic.fuzzyk.png[1].nplane;})
double dot(const std::vector< double > &x, const std::vector< double > &y)
Definition: dot.hpp:10
const Cut CVNNSA
const Cut cut
Definition: exporter_fd.C:30
const Var kProng1NHit([](const caf::SRProxy *sr){return sr->vtx.elastic.fuzzyk.png[0].nhit;})
const Var kESlc([](const caf::SRProxy *sr){return 1000 *deadscale *sr->slc.calE;})
const Var kNPng([](const caf::SRProxy *sr){if(!(sr->vtx.elastic.IsValid)) return 0.0f;return float(sr->vtx.elastic.fuzzyk.npng);})
const Var kPng1CVN([](const caf::SRProxy *sr){return sr->vtx.elastic.fuzzyk.png[0].cvnpart.photonid;})
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
const Var kVtxY
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2142
const Cut kRemidCut
const Var kSumE([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.;if(sr->vtx.elastic.fuzzyk.npng!=2) return-1.;const auto &p1=sr->vtx.elastic.fuzzyk.png[0];const auto &p2=sr->vtx.elastic.fuzzyk.png[1];return(double) 1000 *deadscale *(p1.calE+p2.calE);})
caf::Proxy< float > y
Definition: SRProxy.h:107
caf::Proxy< float > calE
Definition: SRProxy.h:1292
const unsigned int kNMulti
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
const Var kVtxZ
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2073
const Var kDotProd([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-5.;if(sr->vtx.elastic.fuzzyk.npng!=2) return-5.;return sr->vtx.elastic.fuzzyk.png[0].dir.Unit().Dot(sr->vtx.elastic.fuzzyk.png[1].dir.Unit());})
const Cut kNueVtxYBot([](const caf::SRProxy *sr){if(sr->vtx.elastic.vtx.Y()< 0) return true;return false;})
const Var kProng1E([](const caf::SRProxy *sr){return deadscale *1000 *sr->vtx.elastic.fuzzyk.png[0].calE;})
caf::Proxy< unsigned int > firstplane
Definition: SRProxy.h:1305
const double deadscale
const Cut CVNASA
const Cut kPhotonLow([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.npng!=2) return false;return sr->vtx.elastic.fuzzyk.png[0].cvnpart.photonid > 0.6 &&sr->vtx.elastic.fuzzyk.png[1].cvnpart.photonid > 0.6;})
const MultiVar kEHit([](const caf::SRProxy *sr){std::vector< double > ret;ret.push_back(kProng1EHit(sr));ret.push_back(kProng2EHit(sr));return ret;})
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
const Cut kdEdxCut
const Var kCosTotal([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-5.;if(sr->vtx.elastic.fuzzyk.npng!=2) return-5.;double x1=sr->vtx.elastic.fuzzyk.png[0].dir.x;double y1=sr->vtx.elastic.fuzzyk.png[0].dir.y;double z1=sr->vtx.elastic.fuzzyk.png[0].dir.z;double x2=sr->vtx.elastic.fuzzyk.png[1].dir.x;double y2=sr->vtx.elastic.fuzzyk.png[1].dir.y;double z2=sr->vtx.elastic.fuzzyk.png[1].dir.z;TVector3 tot(x1+x2, y1+y2, z1+z2);return tot.Unit().Dot(beamDirND);})
const Cut kNueVtxXEast([](const caf::SRProxy *sr){if(sr->vtx.elastic.vtx.X()< 0) return true;return false;})
const Var kNumPlanes([](const caf::SRProxy *sr){return sr->slc.lastplane-sr->slc.firstplane+1;})
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
const Var kMass([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.;if(sr->vtx.elastic.fuzzyk.npng!=2) return-1.;const auto &p1=sr->vtx.elastic.fuzzyk.png[0];const auto &p2=sr->vtx.elastic.fuzzyk.png[1];const double dot=p1.dir.Unit().Dot(p2.dir.Unit());return deadscale *1000 *sqrt(2 *p1.calE *p2.calE *(1-dot));})
const MultiVar kProngCVN([](const caf::SRProxy *sr){std::vector< double > ret;ret.push_back(kPng1CVN(sr));ret.push_back(kPng2CVN(sr));return ret;})
const Var kMassNoScale([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.;if(sr->vtx.elastic.fuzzyk.npng!=2) return-1.;const auto &p1=sr->vtx.elastic.fuzzyk.png[0];const auto &p2=sr->vtx.elastic.fuzzyk.png[1];const double dot=p1.dir.Unit().Dot(p2.dir.Unit());return 1000 *sqrt(2 *p1.calE *p2.calE *(1-dot))/1.78;})
const Cut CVN_high
const Cut kNueVtxContain([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(fabs(sr->vtx.elastic.vtx.X()) > 180) return false;if(fabs(sr->vtx.elastic.vtx.Y()) > 180) return false;if(fabs(sr->vtx.elastic.vtx.Z())< 50) return false;if(fabs(sr->vtx.elastic.vtx.Z()) > 1000) return false;return true;})
const Cut CVN_max
enum BeamMode string
const Var kNumHits