pi0_xcheck.C
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"
20 
21 using namespace ana;
22 
24 
25 const double deadscale = 0.855;
26 
27 // Put a "NOvA Preliminary" tag in the corner
29 {
30  TLatex* prelim = new TLatex(.9, .95, "NOvA Preliminary");
31  prelim->SetTextColor(kBlue);
32  prelim->SetNDC();
33  prelim->SetTextSize(2/30.);
34  prelim->SetTextAlign(32);
35  prelim->Draw();
36 }
37 
38 
39 const Cut kNueVtxContain([](const caf::SRProxy* sr)
40  {
41  if( !sr->vtx.elastic.IsValid ) return false;
42 
43  if(fabs(sr->vtx.elastic.vtx.X()) > 180) return false;
44  if(fabs(sr->vtx.elastic.vtx.Y()) > 180) return false;
45  if(fabs(sr->vtx.elastic.vtx.Z()) < 50) return false;
46  if(fabs(sr->vtx.elastic.vtx.Z()) > 1000) return false;
47  return true;
48  });
49 
50 const Cut kNueShwContain([](const caf::SRProxy* sr)
51  {
52  if(sr->vtx.elastic.fuzzyk.nshwlid == 0) return false;
53 
54  if(fabs(sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.X()) > 180) return false;
55  if(fabs(sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.Y()) > 180) return false;
56  if(sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.Z() < 200) return false;
57  if(sr->vtx.elastic.fuzzyk.png[0].shwlid.stop.Z() > 1200) return false;
58  return true;
59  });
60 
61 const Cut kRemidCut = SIMPLEVAR(sel.remid.pid) > 0 && SIMPLEVAR(sel.remid.pid) < 0.5;
62 const Cut kdEdxCut = SIMPLEVAR(sand.nue.dedxpng1) < 0.003 && SIMPLEVAR(sand.nue.dedxpng2) < 0.003;
63 
64 const Cut kTwoProngs([](const caf::SRProxy* sr)
65  {
66  if( !sr->vtx.elastic.IsValid ) return false;
67  return sr->vtx.elastic.fuzzyk.npng == 2;
68  });
69 
70 
71 const Cut kMissingPlanesCut([](const caf::SRProxy* sr)
72  {
73  if( !sr->vtx.elastic.IsValid ) return false;
74  for(auto& it: sr->vtx.elastic.fuzzyk.png)
75  if(it.maxplanegap <= 1) return false;
76  return true;
77  });
78 
79 const Cut kContigPlanesCut([](const caf::SRProxy* sr)
80  {
81  if( !sr->vtx.elastic.IsValid ) return false;
82  for(auto& it: sr->vtx.elastic.fuzzyk.png)
83  if(it.maxplanecont <= 4) return false;
84  return true;
85  });
86 
87 const Var kLongestProngVar = SIMPLEVAR(sel.rvp.longtr);
88 const Var kNumHits = SIMPLEVAR(slc.nhit);
89 const Var kNumPlanes([](const caf::SRProxy* sr)
90  {
91  return sr->slc.lastplane - sr->slc.firstplane + 1;
92  });
93 
94 const Var kProng1E([](const caf::SRProxy* sr)
95  {
96  return sr->vtx.elastic.fuzzyk.png[0].calE*deadscale;
97  });
98 
99 const Var kProng2E([](const caf::SRProxy* sr)
100  {
101  return sr->vtx.elastic.fuzzyk.png[1].calE*deadscale;
102  });
103 
104 const Var kProng1L([](const caf::SRProxy* sr)
105  {
106  return sr->vtx.elastic.fuzzyk.png[0].len;
107  });
108 
109 const Var kProng2L([](const caf::SRProxy* sr)
110  {
111  return sr->vtx.elastic.fuzzyk.png[1].len;
112  });
113 
114 const Var kProng1NHit([](const caf::SRProxy* sr)
115  {
116  return sr->vtx.elastic.fuzzyk.png[0].nhit;
117  });
118 
119 const Var kProng2NHit([](const caf::SRProxy* sr)
120  {
121  return sr->vtx.elastic.fuzzyk.png[1].nhit;
122  });
123 
124 const Var kProng1NPlane([](const caf::SRProxy* sr)
125  {
126  return sr->vtx.elastic.fuzzyk.png[0].nplane;
127  });
128 
129 const Var kProng2NPlane([](const caf::SRProxy* sr)
130  {
131  return sr->vtx.elastic.fuzzyk.png[1].nplane;
132  });
133 
134 const Var kDotProd([](const caf::SRProxy* sr)
135  {
136  if( !sr->vtx.elastic.IsValid ) return -1.;
137  if(sr->vtx.elastic.fuzzyk.npng != 2) return -1.;
138  return sr->vtx.elastic.fuzzyk.png[0].dir.Unit().Dot(sr->vtx.elastic.fuzzyk.png[1].dir.Unit());
139  });
140 
141 const Cut kCalE500([](const caf::SRProxy* sr)
142  {
143  if( !sr->vtx.elastic.IsValid ) return false;
144  if(sr->vtx.elastic.fuzzyk.npng < 2 ) return false;
145  if(sr->slc.calE*deadscale >= 0.5) return false;
146  return true;
147  });
148 
149 const Cut kCalE500_1k([](const caf::SRProxy* sr)
150  {
151  if( !sr->vtx.elastic.IsValid ) return false;
152  if(sr->vtx.elastic.fuzzyk.npng < 2) return false;
153  if(sr->slc.calE*deadscale < 0.5 || sr->slc.calE*deadscale >= 1.0) return false;
154  return true;
155  });
156 
157 const Cut kCalE1k_1500([](const caf::SRProxy* sr)
158  {
159  if( !sr->vtx.elastic.IsValid ) return false;
160  if(sr->vtx.elastic.fuzzyk.npng < 2) return false;
161  if(sr->slc.calE*deadscale < 1.0 || sr->slc.calE*deadscale >= 1.5) return false;
162  return true;
163  });
164 
165 const Cut kCalE1500_2k([](const caf::SRProxy* sr)
166  {
167  if( !sr->vtx.elastic.IsValid ) return false;
168  if(sr->vtx.elastic.fuzzyk.npng < 2) return false;
169  if(sr->slc.calE < 1.5*deadscale || sr->slc.calE*deadscale >= 2.0) return false;
170  return true;
171  });
172 
173 const Cut kCalE2k_2500([](const caf::SRProxy* sr)
174  {
175  if( !sr->vtx.elastic.IsValid ) return false;
176  if(sr->vtx.elastic.fuzzyk.npng < 2) return false;
177  if(sr->slc.calE*deadscale < 2.0 || sr->slc.calE*deadscale >= 2.5) return false;
178  return true;
179  });
180 
181 const Cut kCalE2500_3k([](const caf::SRProxy* sr)
182  {
183  if( !sr->vtx.elastic.IsValid ) return false;
184  if(sr->vtx.elastic.fuzzyk.npng < 2) return false;
185  if(sr->slc.calE*deadscale < 2.5 || sr->slc.calE*deadscale >= 3.0) return false;
186  return true;
187  });
188 
189 const Var kMass([](const caf::SRProxy* sr)
190  {
191  if( !sr->vtx.elastic.IsValid ) return -1.;
192  if(sr->vtx.elastic.fuzzyk.npng != 2) return -1.;
193  const auto& p1 = sr->vtx.elastic.fuzzyk.png[0];
194  const auto& p2 = sr->vtx.elastic.fuzzyk.png[1];
195 
196  const double dot = p1.dir.Unit().Dot(p2.dir.Unit());
197  return deadscale*1000*sqrt(2*p1.calE*p2.calE*(1-dot));
198  });
199 
200 const Var kSumE([](const caf::SRProxy* sr)
201  {
202  if( !sr->vtx.elastic.IsValid) return -1.f;
203  if(sr->vtx.elastic.fuzzyk.npng != 2) return -1.f;
204  const auto& p1 = sr->vtx.elastic.fuzzyk.png[0];
205  const auto& p2 = sr->vtx.elastic.fuzzyk.png[1];
206 
207  return p1.calE + p2.calE;
208  });
209 
210 const Cut kTruePi0 = SIMPLEVAR(sand.nue.npi0) > 0;
211 
212 void legend(TH1* data, TH1* mc, TH1* bkgd){
213  TLegend *leg = new TLegend(.60,.64,.84,.84);
214 
215  TLegendEntry* dl = leg->AddEntry(data,"Data","lpe");
216  dl->SetTextColor(data->GetLineColor());
217  TLegendEntry* ml = leg->AddEntry(mc,"MC #pi^{0} signal","l");
218  ml->SetTextColor(mc->GetLineColor());
219  TLegendEntry* bl = leg->AddEntry(bkgd,"MC bkgd","l");
220  bl->SetTextColor(bkgd->GetLineColor());
221  leg->SetBorderSize(0); //no border for legend
222  leg->SetFillColor(0); //fill colour is white
223  leg->Draw();
224 }
225 
226 void FitValues(TF1* data, TF1* mc, TH1* hdata, TH1* hmc) {
227 
228  TLegend *leg_dm = new TLegend(.50, .30, .85, 0.55);
229 
230  char dmean[50], dsigma[50], mmean[50], msigma[50];
231  sprintf(dmean,"Data #mu: #bf{%2.1f} #pm %2.1f MeV", data->GetParameter(1), data->GetParError(1));
232  sprintf(dsigma,"Data #sigma: %2.1f #pm %2.1f MeV", data->GetParameter(2), data->GetParError(2));
233  sprintf(mmean,"MC #mu: #bf{%2.1f} #pm %2.1f MeV", mc->GetParameter(1), mc->GetParError(1));
234  sprintf(msigma,"MC #sigma %2.1f #pm %2.1f MeV", mc->GetParameter(2), mc->GetParError(2));
235 
236 
237  TLegendEntry* dm = leg_dm->AddEntry((TObject*)0, dmean, "");
238  dm->SetTextColor(hdata->GetLineColor());
239  TLegendEntry* ds = leg_dm->AddEntry((TObject*)0, dsigma, "");
240  ds->SetTextColor(hdata->GetLineColor());
241  TLegendEntry* mm = leg_dm->AddEntry((TObject*)0, mmean, "");
242  mm->SetTextColor(hmc->GetLineColor());
243  TLegendEntry* ms = leg_dm->AddEntry((TObject*)0, msigma, "");
244  ms->SetTextColor(hmc->GetLineColor());
245 
246  leg_dm->SetBorderSize(0); //no border for legend
247  leg_dm->SetFillColor(0); //fill colour is white
248  //leg_dm->SetTextSize(2/40.);
249  leg_dm->Draw();
250 
251 }
252 
253 void pi0_xcheck(int period, bool mec = false, bool mrcc = false)
254 {
256  SpectrumLoaderBase* loaderMC;
257 
258  std::string pStr = "";
259  std::string fmc = "";
260  std::string fdata = "";
261 
262  //Datasets from Gavin's original prod2 blessed plots
263 
264  //const std::string fAllMC = "defname: prod_decaf_R16-03-03-prod2reco.a_nd_genie_nonswap_genierw_fhc_nova_v08_full_nue_contain_v1";
265  const std::string fAllMC = "/pnfs/nova/persistent/production/concat/R16-03-03-prod2reco.a/prod_decaf_R16-03-03-prod2reco.a_nd_genie_nonswap_genierw_fhc_nova_v08_full_nue_contain_v1.root";
266  const std::string fAllData = "defname: gsdavies_decaf_R16-03-03-prod2reco.a_nd_numi_p1-3c_nue_contain_goodruns";
267 
268  const std::string fe3MC = "defname: gsdavies_decaf_R16-03-03-prod2reco.a_nd_genie_nonswap_genierw_fhc_nova_v08_e3b3c_nue_contain_v1";
269  const std::string fe3Data = "defname: gsdavies_decaf_R16-03-03-prod2reco.a_nd_numi_e3b3c_nue_contain_goodruns";
270 
271  const std::string fp1MC = "defname: prod_decaf_R16-03-03-prod2reco.a_nd_genie_nonswap_genierw_fhc_nova_v08_period1_nue_contain_v1";
272  const std::string fp1Data = "defname: prod_decaf_R16-03-03-prod2reco.a_nd_numi_period1_nue_contain_goodruns";
273 
274  const std::string fp2MC = "defname: prod_decaf_R16-03-03-prod2reco.a_nd_genie_nonswap_genierw_fhc_nova_v08_period2_nue_contain_v1";
275 
276  const std::string fp2Data = "defname: prod_decaf_R16-03-03-prod2reco.a_nd_numi_period2_nue_contain_goodruns";
277 
278  //example prod3 data files used for testing
279  //cafe pi0_xcheck.C 5
280 
281  const std::string testData = "/pnfs/nova/persistent/production/concat/R17-03-01-prod3reco.d/period5/prod_decaf_R17-03-01-prod3reco.d_nd_numi_fhc_period5_v1_goodruns_numu2017/prod_decaf_R17-03-01-prod3reco.d_nd_numi_fhc_period5_v1_goodruns_numu2017_1_of_20.root";
282 
283  const std::string testMC = "/pnfs/nova/persistent/production/concat/R17-03-01-prod3reco.d/period5/prod_decaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_period5_v1_numu2017/prod_decaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_period5_v1_numu2017_2_of_70.root";
284 
285 
286  switch(period){
287  case 1:
288  fmc = fp1MC;
289  fdata = fp1Data;
290  pStr = "p1";
291  break;
292  case 2:
293  fmc = fp2MC;
294  fdata = fp2Data;
295  pStr = "p2";
296  break;
297  case 3:
298  fmc = fe3MC;
299  fdata = fe3Data;
300  pStr = "e3b3c";
301  break;
302  case 4:
303  fmc = fAllMC;
304  fdata = fAllData;
305  pStr = "p1-3c";
306  break;
307  case 5:
308  fmc = testMC;
309  fdata = testData;
310  pStr = "testing";
311  break;
312  default:
313  fmc = fAllMC;
314  fdata = fAllData;
315  pStr = "p1-3c";
316  break;
317  }
318 
319  const std::string fMC = fmc;
320  const std::string fData = fdata;
321 
322  std::cout << "MC: " << fMC << std::endl;
323  std::cout << "Data: " << fData << std::endl;
324  std::cout << "PStr: " << pStr << std::endl;
325 
326  loader = new SpectrumLoader(fData);
327  loaderMC = new SpectrumLoader(fMC);
328 
330  loaderMC->SetSpillCut(kStandardDQCuts);
331 
332 
333  const Cut sel = kTwoProngs && kNueVtxContain && kNueShwContain && kRemidCut
335 
336  const unsigned int kNumCuts = 7;
337  const Cut sels[kNumCuts] =
338  {sel,
339  sel && kCalE500,
340  sel && kCalE500_1k,
341  sel && kCalE1k_1500,
342  sel && kCalE1500_2k,
343  sel && kCalE2k_2500,
344  sel
345  };
346 
347 
348  struct PlotDef
349  {
350  std::string suffix;
352  Binning bins;
353  Var var;
354  };
355 
356  const unsigned int kNumPlots = 19;
358  {{"mass1", "M_{#gamma#gamma} (MeV)", Binning::Simple(50, 0, 500), kMass},
359  {"mass2", "M_{#gamma#gamma} (MeV)", Binning::Simple(50, 0, 500), kMass},
360  {"mass3", "M_{#gamma#gamma} (MeV)", Binning::Simple(50, 0, 500), kMass},
361  {"mass4", "M_{#gamma#gamma} (MeV)", Binning::Simple(50, 0, 500), kMass},
362  {"mass5", "M_{#gamma#gamma} (MeV)", Binning::Simple(50, 0, 500), kMass},
363  {"mass6", "M_{#gamma#gamma} (MeV)", Binning::Simple(50, 0, 500), kMass},
364  {"mass_eta", "M_{#gamma#gamma} (MeV)", Binning::Simple(100, 0, 1000), kMass},
365  {"calE", "Prong calorimetric energy (GeV)", Binning::Simple(20, 0, 2), kProng1E},
366  {"junk", "junk", Binning::Simple(20, 0, 2), kProng2E},
367  {"length", "Prong length (cm)", Binning::Simple(25, 0, 500), kProng1L},
368  {"junk", "junk", Binning::Simple(25, 0, 500), kProng2L},
369  {"nhits", "Number of hits in prong", Binning::Simple(50, 0, 100), kProng1NHit},
370  {"junk", "junk", Binning::Simple(50, 0, 100), kProng2NHit},
371  {"nplanes", "Number of prong planes", Binning::Simple(50, 0, 100), kProng1NPlane},
372  {"junk", "junk", Binning::Simple(50, 0, 100), kProng2NPlane},
373  {"dot", "cos#theta", Binning::Simple(20, 0, 1), kDotProd},
374  {"Epi0", "E_{#gamma#gamma} (GeV)", Binning::Simple(50, 0, 2), kSumE},
375  {"lid", "LID", Binning::Simple(24,-0.1,1.1), kElecID},
376  {"lem", "LEM", Binning::Simple(24,-0.1,1.1), kLEM}
377  };
378 
379  Spectrum* dataSpects[kNumPlots];
380  Spectrum* mcSpects[kNumPlots];
381  Spectrum* bkgSpects[kNumPlots];
382 
383  for(unsigned int i = 0; i < kNumPlots; ++i){
384  PlotDef def = plots[i];
385 
386  if(i < 7){
387  dataSpects[i] = new Spectrum(def.label, def.bins, *loader, def.var, sels[i]);
388  if(mec){
389  mcSpects[i] = new Spectrum(def.label, def.bins, *loaderMC, def.var, sels[i], kNoShift, kTuftsWeightCC);
390  bkgSpects[i] = new Spectrum(def.label, def.bins, *loaderMC, def.var, sels[i] && !kTruePi0, kNoShift, kTuftsWeightCC);
391  }
392  else{
393  mcSpects[i] = new Spectrum(def.label, def.bins, *loaderMC, def.var, sels[i]);
394  bkgSpects[i] = new Spectrum(def.label, def.bins, *loaderMC, def.var, sels[i] && !kTruePi0);
395 
396  }
397  }
398  else{
399  dataSpects[i] = new Spectrum(def.label, def.bins, *loader, def.var, sel);
400  if(mec){
401  std::cout << "MEC EVENTS. APPLY MEC REWEIGHTING" << std::endl;
402  mcSpects[i] = new Spectrum(def.label, def.bins, *loaderMC, def.var, sel, kNoShift, kTuftsWeightCC);
403  bkgSpects[i] = new Spectrum(def.label, def.bins, *loaderMC, def.var, sel && !kTruePi0, kNoShift, kTuftsWeightCC);
404  }
405  else{
406  std::cout << "IGNORE MEC EVENTS. NO MEC REWEIGHTING" << std::endl;
407  mcSpects[i] = new Spectrum(def.label, def.bins, *loaderMC, def.var, sel);
408  bkgSpects[i] = new Spectrum(def.label, def.bins, *loaderMC, def.var, sel && !kTruePi0);
409  }
410  }
411  }
412 
413  loader->Go();
414  loaderMC->Go();
415 
416  const std::string mrccStr = mrcc ? "_mrcc" : "";
417  const std::string mecStr = mec ? "MEC" : "noMEC";
418  TFile fout(("pi0_xcheck_output_"+pStr+mrccStr+"_"+mecStr+".root").c_str(), "RECREATE");
419 
420  const double deadscale = 1.522;
421 
422  for(unsigned int i = 0; i < kNumPlots; ++i){
423 
424  if(plots[i].suffix == "junk") continue;
425 
426  TCanvas *c = new TCanvas((plots[i].suffix+mrccStr).c_str());
427  //DataMCComparison(*dataSpects[i], *mcSpects[i]);
428  TH1* hMC = mcSpects[i]->ToTH1(dataSpects[i]->POT());
429  hMC->SetLineColor(kTotalMCColor);
430 
431  TH1* hData = dataSpects[i]->ToTH1(dataSpects[i]->POT());
432  //hData->Sumw2();
433  hData->SetMarkerStyle(kFullCircle);
434 
435  if(i < 7){
436  TF1 *fData = new TF1("fData","gaus",100,200);
437  fData->SetLineColor(kBlack);
438  hData->Fit("fData","R0","");
439 
440  TF1 *fMC = new TF1("fMC","gaus",100,200);
441  fMC->SetLineColor(kRed);
442  hMC->Fit("fMC","R0+","");
443 
444  if(hMC->GetMaximum() > hData->GetMaximum() + sqrt(hData->GetMaximum())){
445  hMC->Draw("hist");
446  hData->Draw("ep same");
447  }
448  else{
449  hData->Draw("ep");
450  hMC->Draw("hist same");
451  }
452 
453  TH1* h = bkgSpects[i]->ToTH1(dataSpects[i]->POT());
454  h->SetLineColor(kBlue);
455  h->Draw("hist same");
456 
457  FitValues(fData, fMC, hData, hMC);
458  legend(hData,hMC,h);
459  Preliminary();
460  double mData = fData->GetParameter(1);
461  double meData = fData->GetParError(1);
462  double sData = fData->GetParameter(2);
463  double seData = fData->GetParError(2);
464 
465  std::cout << "Iteration: " << i << "\n";
466  std::cout << "Data Mean: " << mData << " +- " << meData << "\n";
467  std::cout << "Data Sigma: " << sData << " +- " << seData << "\n";
468 
469  double mMC = fMC->GetParameter(1);
470  double meMC = fMC->GetParError(1);
471  double sMC = fMC->GetParameter(2);
472  double seMC = fMC->GetParError(2);
473 
474  std::cout << "MC Mean: " << mMC << " +- " << meMC << "\n";
475  std::cout << "MC Sigma: " << sMC << " +- " << seMC << "\n";
476  }
477  else{
478  if(hMC->GetMaximum() > hData->GetMaximum() + sqrt(hData->GetMaximum())){
479  hMC->Draw("hist");
480  hData->Draw("ep same");
481  }
482  else{
483  hData->Draw("ep");
484  hMC->Draw("hist same");
485  }
486 
487  TH1* h = bkgSpects[i]->ToTH1(dataSpects[i]->POT());
488  h->SetLineColor(kBlue);
489  h->Draw("hist same");
490 
491  legend(hData,hMC,h);
492  Preliminary();
493  }
494 
495  fout.cd();
496  c->Write();
497  c->Print((pStr+"/data_mc_"+plots[i].suffix+mrccStr+"_"+pStr+"_"+mecStr+".png").c_str());
498  c->Print((pStr+"/data_mc_"+plots[i].suffix+mrccStr+"_"+pStr+"_"+mecStr+".root").c_str());
499  }
500 }
caf::Proxy< size_t > npng
Definition: SRProxy.h:2038
caf::Proxy< unsigned int > nshwlid
Definition: SRProxy.h:2040
static constexpr Double_t ms
Definition: Munits.h:192
const Var kLEM
PID
Definition: Vars.cxx:26
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
enum BeamMode kRed
const double deadscale
Definition: pi0_xcheck.C:25
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;})
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
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const Var kDotProd([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.;if(sr->vtx.elastic.fuzzyk.npng!=2) return-1.;return sr->vtx.elastic.fuzzyk.png[0].dir.Unit().Dot(sr->vtx.elastic.fuzzyk.png[1].dir.Unit());})
const Var kProng2NPlane([](const caf::SRProxy *sr){return sr->vtx.elastic.fuzzyk.png[1].nplane;})
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
const Var kSumE([](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 p1.calE+p2.calE;})
const Var kProng1NHit([](const caf::SRProxy *sr){return sr->vtx.elastic.fuzzyk.png[0].nhit;})
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
T sqrt(T number)
Definition: d0nt_math.hpp:156
const Var kProng2E([](const caf::SRProxy *sr){return sr->vtx.elastic.fuzzyk.png[1].calE *deadscale;})
void Preliminary()
Definition: pi0_xcheck.C:28
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 Var kProng1NPlane([](const caf::SRProxy *sr){return sr->vtx.elastic.fuzzyk.png[0].nplane;})
TFile * fmc
Definition: bdt_com.C:14
void SetSpillCut(const SpillCut &cut)
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;})
const Var kLongestProngVar
Definition: pi0_xcheck.C:87
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
const XML_Char const XML_Char * data
Definition: expat.h:268
const SpillCut kTightBeamQualityCuts([](const caf::SRSpillProxy *s){if(s->ismc) return true; if(s->trigger==2) return true;if(s->spilltimesec==0 &&s->deltaspilltimensec==0 &&s->widthx==0) return bool(s->isgoodspill);if(std::abs(s->deltaspilltimensec) > 0.5e9) return false;if(s->spillpot< 2e12) return false;if(s->hornI< -202|| s->hornI >-198) return false;if(s->posx< -2.00|| s->posx >+2.00) return false;if(s->posy< -2.00|| s->posy >+2.00) return false;return kBeamWidthCut(s);})
Definition: SpillCuts.h:10
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
void FitValues(TF1 *data, TF1 *mc, TH1 *hdata, TH1 *hmc)
Definition: pi0_xcheck.C:226
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2043
virtual void Go()=0
Load all the registered spectra.
const Cut kCalE2k_2500([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.npng< 2) return false;if(sr->slc.calE *deadscale< 2.0||sr->slc.calE *deadscale >=2.5) return false;return true;})
const Cut kTwoProngs([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;return sr->vtx.elastic.fuzzyk.npng==2;})
const Cut kdEdxCut
Definition: pi0_xcheck.C:62
const Cut sels[kNumSels]
Definition: vars.h:44
const int kNumPlots
Definition: GetSpectra.h:1
const Var kProng2NHit([](const caf::SRProxy *sr){return sr->vtx.elastic.fuzzyk.png[1].nhit;})
const std::vector< Plot > plots
const SpillCut kStandardDQCuts([](const caf::SRSpillProxy *spill){if(spill->dcmedgematchfrac==0 &&spill->fracdcm3hits==0 &&spill->nmissingdcmslg==0) return bool(spill->isgoodspill); if(spill->det==caf::kNEARDET && (spill->fracdcm3hits > 0.45|| spill->nmissingdcms > 0)) return false; if(spill->eventincomplete) return false; if(spill->det==caf::kFARDET && spill->nmissingdcmslg > 0) return false; if(spill->det==caf::kFARDET && !spill->ismc && spill->dcmedgematchfrac<=0.2) return false;return true;})
Cut out events with a noisy detector or with parts missing.
Definition: SpillCuts.h:16
const Cut kCalE1500_2k([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.npng< 2) return false;if(sr->slc.calE< 1.5 *deadscale||sr->slc.calE *deadscale >=2.0) return false;return true;})
void pi0_xcheck(int period, bool mec=false, bool mrcc=false)
Definition: pi0_xcheck.C:253
caf::StandardRecord * sr
caf::Proxy< unsigned int > lastplane
Definition: SRProxy.h:1309
const Cut kCalE1k_1500([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.npng< 2) return false;if(sr->slc.calE *deadscale< 1.0||sr->slc.calE *deadscale >=1.5) return false;return true;})
const Cut kCalE2500_3k([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.npng< 2) return false;if(sr->slc.calE *deadscale< 2.5||sr->slc.calE *deadscale >=3.0) return false;return true;})
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:610
TLatex * prelim
Definition: Xsec_final.C:133
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;})
A very simple service to remember what detector we&#39;re working in.
const Var kNumHits
Definition: pi0_xcheck.C:88
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
const Var kProng1E([](const caf::SRProxy *sr){return sr->vtx.elastic.fuzzyk.png[0].calE *deadscale;})
Base class for the various types of spectrum loader.
static constexpr Double_t mm
Definition: Munits.h:136
std::string label
const Cut kRemidCut
Definition: pi0_xcheck.C:61
const Binning bins
Definition: NumuCC_CPiBin.h:8
double dot(const std::vector< double > &x, const std::vector< double > &y)
Definition: dot.hpp:10
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
const Var kProng1L([](const caf::SRProxy *sr){return sr->vtx.elastic.fuzzyk.png[0].len;})
const Var kProng2L([](const caf::SRProxy *sr){return sr->vtx.elastic.fuzzyk.png[1].len;})
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2142
const Var kTuftsWeightCC
Definition: XsecTunes.h:31
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
caf::Proxy< float > calE
Definition: SRProxy.h:1292
const Cut kCalE500([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.npng< 2) return false;if(sr->slc.calE *deadscale >=0.5) return false;return true;})
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
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;})
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2073
void legend(TH1 *data, TH1 *mc, TH1 *bkgd)
Definition: pi0_xcheck.C:212
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 kNumPlanes([](const caf::SRProxy *sr){return sr->slc.lastplane-sr->slc.firstplane+1;})
caf::Proxy< unsigned int > firstplane
Definition: SRProxy.h:1305
const Color_t kTotalMCColor
Definition: Style.h:16
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
enum BeamMode kBlue
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
Binning bins
const Cut kCalE500_1k([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.npng< 2) return false;if(sr->slc.calE *deadscale< 0.5||sr->slc.calE *deadscale >=1.0) return false;return true;})
enum BeamMode string