fa_fd_data_mc_systs.C
Go to the documentation of this file.
1 #include "CAFAna/Cuts/Cuts.h"
4 #include "CAFAna/Cuts/NumuCuts.h"
5 #include "CAFAna/Core/Loaders.h"
7 #include "CAFAna/Core/Spectrum.h"
9 #include "CAFAna/Vars/Vars.h"
10 #include "CAFAna/Vars/NumuVars.h"
12 #include "CAFAna/Analysis/Calcs.h"
13 #include "CAFAna/Analysis/Plots.h"
14 
21 #include "CAFAna/Systs/NumuSysts.h"
22 #include "CAFAna/Vars/FitVars.h"
23 
24 #include "TStyle.h"
25 #include "TCanvas.h"
26 #include "TGraph.h"
27 #include "TGraphAsymmErrors.h"
28 #include "TH1.h"
29 #include "TH2.h"
30 #include "TLegend.h"
31 #include "TLatex.h"
32 #include "TArrow.h"
33 
34 using namespace ana;
35 
37  const std::vector<const ISyst*>& systs,
39  double pot,
40  int col = -1, int errCol = -1,
41  float ymax = 60.0);
42 
44  std::vector<Spectrum> upShifts,
45  std::vector<Spectrum> downShifts,
46  double pot,
47  int col = -1, int errCol = -1,
48  float ymax = 60.0);
49 
50 void Preliminary();
51 void RemoveBadPoints(TGraphAsymmErrors* gr, TH1* pred);
52 void draw_ratio(TCanvas* c1, TH1* MC, TH1* Unosc, TH1* data, TGraphAsymmErrors* grData);
53 
54 // And some binning definitions
55 const Binning kXYBins = Binning::Simple(12,-9,9);
56 const Binning kZBins = Binning::Simple(12,0,60);
58 
60 {
61  std::cout << "Drawing plots..." << std::endl;
62 
63  gStyle->SetStatFormat("6.3g"); // Significant figures. Fix if this is not satisfactory
64  gStyle->SetPaintTextFormat("6.3g");
65  gStyle->SetFitFormat("6.3g");
66 
67  // Matthew's all-time favourite hadded files
68  std::string directory = "/nova/ana/users/tamsett/data/sam_definitions/";
69  // Data
70  std::string nd_data = directory+"hadd_prod_decaf_S15-05-22a_nd_numi_numu_contain_goodruns.root";
71  std::string fd_data = directory+"hadd_prod_restricteddecaf_S15-05-22a_fd_numi_full_numu_contain_goodruns_all.root";
72 
73  // Baseline MC
74  std::string nd_nonswap = directory+"hadd_rocco_decaf_S15-05-22_nd_genie_fhc_nonswap_genierw_numu_contain.root";
75  // std::string fd_nonswap = directory+"hadd_prod_decaf_S15-05-22a_fd_genie_fhc_nonswap_numu_contain_matchedruns.root";
76  // std::string fd_fluxswap = directory+"hadd_prod_decaf_S15-05-22a_fd_genie_fhc_fluxswap_numu_contain_matchedruns.root";
77  // std::string fd_tau = directory+"hadd_prod_decaf_S15-05-22a_fd_genie_fhc_tau_numu_contain_matchedruns.root";
78 
79  // New Chris' hadded FD MC
80  std::string directory_new = "/nova/ana/nu_e_ana/concat/";
81  std::string fd_nonswap_new = directory_new+"prod_decaf_S15-05-22a_fd_genie_fhc_nonswap_fdfluxv08_numu_contain.root";
82  std::string fd_fluxswap_new = directory_new+"prod_decaf_S15-05-22a_fd_genie_fhc_fluxswap_fdfluxv08_numu_contain.root";
83  std::string fd_tau_new = directory_new+"prod_decaf_S15-05-22a_fd_genie_fhc_tau_fdfluxv08_numu_contain.root";
84 
92 
93  SpectrumLoader loader(fd_data);
94 
97 
100 
101  calc.SetTh23(TMath::ASin(sqrt(0.428474))); // Updated 2015-01-13
102  calc.SetDmsq32(2.52314e-3); // Same
103 
104  const Cut kCCEshiftBelow5(
105  [](const caf::SRProxy* sr)
106  {
107  float muE = sr->energy.numu.recomuonE;
108  float hadE = sr->energy.numu.recotrkcchadE;
109  return ( (sr->hdr.ismc && sr->energy.numu.trkccE < 5.0)
110  || (!sr->hdr.ismc && muE+1.21f*hadE < 5.0) );
111  });
112 
113  const Cut FDCut = kNumuFD && kCCEshiftBelow5;
114 
115  struct Plot
116  {
119  Binning bins;
120  Var var;
121  };
122 
123  const int kNumPlots = 32; // Doing all would take ages!! Consider removing some
124 
125  Plot plots[kNumPlots] = {
126  {";Reconstructed Neutrino Energy (GeV);Events / (0.25 GeV)", "ccEshift", kEnergyBinning, kCCEshift},
127  {";#Deltat from t_{0} (#mus);Events / (12 #mus)", "tus", Binning::Simple(35,26.125,446.125), kSliceTime},
128  {";Reconstructed Neutrino Energy (GeV);Events / (0.25 GeV)", "ccE", kEnergyBinning, kCCE},
129  {";Calorimetric Energy (GeV);Events / (0.25 GeV)", "calE", kEnergyBinning, SIMPLEVAR(slc.calE)},
130  {";Hadronic Energy (GeV);Events / (0.25 GeV)", "hadEshift", kEnergyBinning, kHadEshift},
131  {";Calorimetric Hadronic Energy (GeV);Events / (0.25 GeV)", "hadcalE", kEnergyBinning, kNumuHadCalE},
132  {";Maximum y (m);Events / (1.5 m)", "maxy", kXYBins, kSlcMaxY},
133  {";Number of Tracks in Slice;Events / bin", "nkal", Binning::Simple(12, 1, 13), kNKalman},
134  {";Total Number of Hits;Events / bin", "nhit", Binning::Simple(8, 0, 400), kNHit},
135  {";Cosmic Rejection BDT Output;Events / bin", "bdt", Binning::Simple(10, 0.53, 0.73), kNumuContPID},
136  {";Track Start X Position (m);Events / (1.5 m)", "trkStartX", kXYBins, kTrkStartX},
137  {";Track Start Y Position (m);Events / (1.5 m)", "trkStartY", kXYBins, kTrkStartY},
138  {";Track Start Z Position (m);Events / (5 m)", "trkStartZ", kZBins, kTrkStartZ},
139  {";Track End X Position (m);Events / (1.5 m)", "trkEndX", kXYBins, kTrkEndX},
140  {";Track End Y Position (m);Events / (1.5 m)", "trkEndY", kXYBins, kTrkEndY},
141  {";Track End Z Position (m);Events / (5 m)", "trkEndZ", kZBins, kTrkEndZ},
142  {";Number of Planes to Front;Events / bin", "nPlanesFront", Binning::Simple(15,0,900), kNplanesToFront},
143  {";Number of Planes to Back;Events / bin", "nPlanesBack", Binning::Simple(15,0,900), kNplanesToBack},
144  {";Number of Cells to Edge;Events / bin", "nCellsFromEdge", Binning::Simple(18,0,180), kNcellsFromEdge},
145  {";Basic Track Forward Distance to Edge (m);Events / (5 m)", "cosFwdDist", kZBins, kCosFwdDist},
146  {";Basic Track Backward Distance to Edge (m);Events / (5 m)", "cosBackDist", kZBins, kCosBackDist},
147  {";Kalman Track Forward Distance to Edge (m);Events / (5 m)", "kalFwdDist", kZBins, kKalmanFwdDist},
148  {";Kalman Track Backward Distance to Edge (m);Events / (5 m)", "kalBackDist", kZBins, kKalmanBackDist},
149  {";Slice Duration (ns);Events / (200 ns)", "sliceDuration", Binning::Simple(15,0,3000), kSliceDuration},
150  {";Number of Hits in Main Track;Events / bin", "trkNhits", Binning::Simple(20,0,500), kTrkNhits},
151  {";Number of Hits not in Main Track;Events / bin", "hadNhit", Binning::Simple(8, 0, 160), kHadNHit},
152  {";Length of Main Track (m);Events / (2 m)", "trkLength", Binning::Simple(15,0,30), kTrkLength},
153  {";Scattering Log-Likelihood;Events / bin", "scatLL", Binning::Simple(11,-0.2,0.3), kReMIdScatLLH},
154  {";dE/dx Log-Likelihood;Events / bin", "dedxLL", Binning::Simple(12,-0.1,0.5), kReMIdDEDxLLH},
155  {";cos#theta_{beam};Events / bin", "angkal", Binning::Simple(12, 0.3, 1), kNumuCosRejAngleKal},
156  {";cos#theta_{Z};Events / bin", "dirZ", Binning::Simple(12, 0.3, 1), kDirZ},
157  {";Muon ID;Events / bin", "remid", Binning::Simple(14,0.75-1./160,1.+1./160), kRemID}
158  };
159 
160  Spectrum* hData[kNumPlots];
161  Spectrum* hCosmic[kNumPlots];
162 
163  PredictionInterp* prediction[kNumPlots];
164  std::vector<const ISyst*> systs = getAllNumuSystsSA();
165 
166  for(int i = 0; i < kNumPlots; ++i){
167 
168  Plot p = plots[i];
169 
170  hData[i] = new Spectrum(p.label, p.bins, loader, p.var, i == 1? FDCut : FDCut && kInBeamSpillShifted);
171  hCosmic[i] = new Spectrum(p.label, p.bins, loader, p.var, FDCut && kInTimingSideband, kNoShift, kTimingSidebandWeight);
172 
173  HistAxis axis( p.label, p.bins, p.var );
174  NumuExtrapGenerator NumuCCExtrap( axis, FDCut , kNumuND, kNoShift, kNDPosRwtVar );
175  prediction[i] = new PredictionInterp(systs, &calc, NumuCCExtrap, loaders);
176  }
177 
178  // GO!
179  loaders.Go();
180  loader.Go();
181 
182  double ymax;
183  bool drawUnosc;
184 
185  for(int i = 0; i < kNumPlots; ++i)
186  {
187  std::string fullName = plots[i].fname+"_unblind";
188 
189  TCanvas* c = new TCanvas("c",(fullName).c_str());
190 
191  c->SetBottomMargin(0.15);
192  c->SetLeftMargin(0.15);
193 
194  TH1* hSig = hData[i]->ToTH1(hData[i]->POT());
195  hSig->SetName(("Sig_"+fullName).c_str());
196  hSig->GetXaxis()->CenterTitle();
197  hSig->GetYaxis()->CenterTitle();
198 
199  TGraphAsymmErrors* grSig = GraphWithPoissonErrors(hSig, true, true);
200 
201  hSig->SetLineColor(0);
202  hSig->SetLineWidth(1);
203  hSig->SetMarkerColor(0);
204  hSig->GetYaxis()->SetRangeUser(0,hSig->GetMaximum()*1.12+4.);
205 
206  grSig->SetMarkerColor(kBlack);
207  grSig->SetLineColor(kBlack);
208  grSig->SetLineWidth(2);
209  grSig->SetName(("grSig_"+fullName).c_str());
210 
211  TH1* hBKG = hCosmic[i]->ToTH1(hData[i]->POT());
212  hBKG->SetLineWidth(1);
213  hBKG->SetLineColor(kGray+1);
214  hBKG->SetFillColor(kGray+1);
215  hBKG->SetName(("BKG_"+fullName).c_str());
216 
217  if(i == 1)
218  {
219  hBKG->Reset(); // not for the timing peak
220  hBKG->SetLineColor(10);
221  hBKG->SetFillColor(10);
222  hBKG->SetLineWidth(0);
223  }
224 
225  TH1* hExp = prediction[i]->Predict(&calc).ToTH1(hData[i]->POT());
226  TH1* hUnosc = prediction[i]->PredictUnoscillated().ToTH1(hData[i]->POT());
227  TH1* hNue = prediction[i]->PredictComponent(&calc, ana::Flavors::kAllNuE, ana::Current::kCC, ana::Sign::kBoth).ToTH1(hData[i]->POT());
228  TH1* hNutau = prediction[i]->PredictComponent(&calc, ana::Flavors::kAllNuTau, ana::Current::kCC, ana::Sign::kBoth).ToTH1(hData[i]->POT());
229  TH1* hNC = prediction[i]->PredictComponent(&calc, ana::Flavors::kAll, ana::Current::kNC, ana::Sign::kBoth).ToTH1(hData[i]->POT());
230 
231  hNC->Add(hNue); // All background (but cosmics) drawn together
232  hNC->Add(hNutau);
233  hNC->SetLineColor(kGreen+2);
234  hNC->SetLineWidth(2);
235  hNC->SetName(("NC_"+fullName).c_str());
236 
237  hExp->SetName(("Exp_"+fullName).c_str());
238  hExp->SetLineColor(kRed);
239  hExp->SetLineStyle(1);
240 
241  hUnosc->SetName(("Unosc_"+fullName).c_str());
242  hUnosc->SetLineColor(kBlue);
243  hUnosc->SetLineStyle(2);
244 
245  if(i<=1)
246  drawUnosc = true;
247  else
248  drawUnosc = false;
249 
250  ymax = 2.25 * TMath::Max(hExp->GetBinContent(hExp->GetMaximumBin()),
251  hSig->GetBinContent(hSig->GetMaximumBin()));
252 
253  if(i==1 || i>=29) ymax = ymax * 0.75;
254 
255  myPlotWithSystErrorBand(prediction[i], hCosmic[i], systs,
256  &calc, hData[i]->POT(), kRed, kRed-9, ymax);
257 
258  hBKG->Add(hNC); //It's important not to add the other currents to the expectation
259 
260  if(i != 1) // important not to do this trick for the peak
261  RemoveBadPoints(grSig, hExp);
262 
263  hBKG->Draw("hist same ][");
264  grSig->SetMarkerStyle(20);
265  c->RedrawAxis();
266  grSig->Draw("P");
267 
268  TLegend* ll = new TLegend(0.18,0.68,0.48,0.88);
269  ll->AddEntry(grSig,"FD Data","lp");
270  ll->AddEntry(hExp,"Best-fit prediction","lf");
271  if(i != 1)
272  ll->AddEntry(hBKG,"Background","f");
273  ll->SetFillColor(0);
274  ll->Draw();
275 
276  hExp->GetXaxis()->SetRangeUser(hSig->GetXaxis()->GetXmin(),
277  hSig->GetXaxis()->GetXmax()); // Fix weird mean bug
278 
279  // chi2 = LogLikelihood(hExp, hSig);
280  TLatex* ExposureLabel = new TLatex(0.55,0.80,"2.74#times10^{20} POT-equiv." );
281 
282  ExposureLabel->SetNDC();
283  ExposureLabel->SetTextSize(0.05);
284  ExposureLabel->Draw();
285 
286  Preliminary();
287 
288  c->SaveAs(("plots/"+fullName+".pdf").c_str());
289  c->SaveAs(("plots/"+fullName+".root").c_str());
290 
291  if(drawUnosc)
292  {
293  ymax = 1.2 * hUnosc->GetBinContent(hUnosc->GetMaximumBin());
294  myPlotWithSystErrorBand(prediction[i], hCosmic[i], systs,
295  &calc, hData[i]->POT(), kRed, kRed-9, ymax);
296 
297  hBKG->Draw("hist same ][");
298  grSig->SetMarkerStyle(20);
299  c->RedrawAxis();
300  grSig->Draw("P");
301 
302  ll->SetX1NDC(0.55);
303  ll->SetY1NDC(0.45);
304  ll->SetX2NDC(0.85);
305  ll->SetY2NDC(0.65);
306  ll->AddEntry(hUnosc,"Unoscillated","l");
307  ll->Draw();
308 
309  TLatex* ExposureLabel = new TLatex(0.55,0.80,"2.74#times10^{20} POT-equiv." );
310 
311  ExposureLabel->SetNDC();
312  ExposureLabel->SetTextSize(0.05);
313  ExposureLabel->Draw();
314 
315  Preliminary();
316 
317  c->SaveAs(("plots/"+fullName+"_wUnosc.pdf").c_str());
318  c->SaveAs(("plots/"+fullName+"_wUnosc.root").c_str());
319 
320  c->Clear();
321 
322  // Remove the background "everywhere"
323  // That way we're essentially plotting the disappearance probability
324  for(int k = 1; k <= hExp->GetNbinsX(); k++)
325  {
326  hUnosc->SetBinContent(k, hUnosc->GetBinContent(k) - hBKG->GetBinContent(k) );
327  hExp->SetBinContent(k, hExp->GetBinContent(k) - hBKG->GetBinContent(k) );
328  }
329  for(int k = 0; k < grSig->GetN(); k++)
330  {
331  double x, y;
332  grSig->GetPoint(k, x, y);
333  y = y - hBKG->GetBinContent(hBKG->FindBin(x));
334  grSig->SetPoint(k, x, y);
335  }
336 
337  draw_ratio(c, hExp, hUnosc, hSig, grSig);
338 
339  Preliminary();
340 
341  c->SaveAs(("plots/"+fullName+"_wUnosc_subtractBKG.pdf").c_str());
342  c->SaveAs(("plots/"+fullName+"_wUnosc_subtractBKG.root").c_str());
343  }
344 
345  c->Close();
346 
347  std::cout << "All plots finished!" << std::endl;
348 
349  } // End loop over i
350 
351 } // End of function
352 
353 
354 //----------------------------------------------------------------------
355 
356 void Preliminary() // For some reason it doesn't include it automatically
357 {
358  TLatex* prelim = new TLatex(.9, .95, "NOvA Preliminary");
359  prelim->SetTextColor(kBlue);
360  prelim->SetNDC();
361  prelim->SetTextSize(2/30.);
362  prelim->SetTextAlign(32);
363  prelim->Draw();
364 }
365 
366 //----------------------------------------------------------------------
367 
368 void draw_ratio(TCanvas* c1, TH1* MC, TH1* Unosc, TH1* data, TGraphAsymmErrors* grData)
369 {
370  grData->SetMarkerStyle(20);
371 
372  TH1* Ratio2 = (TH1*)data->Clone("Ratio2");
373  Ratio2->Divide(Unosc);
374  TGraphAsymmErrors* grRatio = (TGraphAsymmErrors*)grData->Clone("grRatio");
375 
376  Ratio2->GetYaxis()->SetRangeUser(-0.2,3.0);
377  Ratio2->GetYaxis()->SetTitle("Data / Unoscillated prediction");
378 
379  Ratio2->SetMarkerColor(10);
380  Ratio2->SetLineColor(10);
381  Ratio2->SetMarkerStyle(20);
382  Ratio2->GetYaxis()->SetNdivisions(506);
383  Ratio2->GetXaxis()->SetLabelSize(0.06);
384  Ratio2->GetYaxis()->SetLabelSize(0.06);
385  Ratio2->GetXaxis()->SetTitleOffset(1.00);
386  Ratio2->GetYaxis()->SetTitleOffset(1.00);
387  Ratio2->Draw("hist");
388 
389  for(int i = 0; i < grRatio->GetN(); ++i )
390  {
391  double x,y,dy_low,dy_high;
392  grRatio->GetPoint(i,x,y);
393  dy_low = grRatio->GetErrorYlow(i);
394  dy_high = grRatio->GetErrorYhigh(i);
395 
396  double fUnosc = Unosc->GetBinContent(Unosc->FindBin(x));
397 
398  if(fUnosc > 0){
399  grRatio->SetPoint(i, x, y / fUnosc);
400  grRatio->SetPointEYhigh(i, dy_high / fUnosc);
401  grRatio->SetPointEYlow(i, dy_low / fUnosc);
402  if((y+dy_high)/fUnosc > 3)
403  {
404  TArrow* arr = new TArrow(x,y/fUnosc,x,3.,0.03);
405  arr->SetLineColor(kBlack);
406  arr->SetLineWidth(2);
407  arr->Draw();
408  }
409  }
410  else if(y == 0.)
411  {
412  grRatio->SetPoint(i, x, 1.);
413  grRatio->SetPointEYhigh(i, 0);
414  grRatio->SetPointEYlow(i, 0);
415  }
416  }
417 
418  TH1* hProb = (TH1*)MC->Clone("hProb");
419  hProb->Divide(Unosc);
420  hProb->SetLineColor(kBlue);
421  hProb->SetFillColor(10);
422  hProb->Draw("hist same");
423 
424  TGraph* g1 = new TGraph;
425  g1->SetPoint(0,-1e9,1);
426  g1->SetPoint(1,1e9,1);
427  g1->SetLineStyle(2);
428  g1->SetLineColor(kBlack);
429  g1->Draw("L");
430 
431  grRatio->Draw("P");
432 
433  TLatex* ExposureLabel = new TLatex(0.55,0.80,"2.74#times10^{20} POT-equiv." );
434 
435  ExposureLabel->SetNDC();
436  ExposureLabel->SetTextSize(0.05);
437  ExposureLabel->Draw();
438 
439 }
440 
441 //----------------------------------------------------------------------
442 
444  const std::vector<const ISyst*>& systs,
446  double pot,
447  int col, int errCol, float ymax)
448 {
449  Spectrum nom = pred->Predict(calc);
450 
451  std::vector<Spectrum> ups, dns;
452 
453  for(const ISyst* syst: systs){
454  SystShifts shifts;
455  shifts.SetShift(syst, +1);
456  ups.push_back(pred->PredictSyst(calc, shifts));
457  shifts.SetShift(syst, -1);
458  dns.push_back(pred->PredictSyst(calc, shifts));
459  }
460 
461  myPlotWithSystErrorBand(nom, bkg, ups, dns, pot, col, errCol, ymax);
462 }
463 
464 //----------------------------------------------------------------------
465 
467  std::vector<Spectrum> upShifts,
468  std::vector<Spectrum> downShifts,
469  double pot,
470  int col, int errCol, float ymax)
471 {
472 
473  TH1* nom = nominal.ToTH1(pot);
474  TH1* hBKG = bkg->ToTH1(pot);
475 
476  std::vector<TH1*> ups, dns;
477 
478  for(const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(pot));
479  for(const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot));
480 
481  TGraphAsymmErrors* g = new TGraphAsymmErrors;
482 
483  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
484  const double y = nom->GetBinContent(binIdx);
485  const double y_bkg = hBKG->GetBinContent(binIdx);
486  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y + y_bkg);
487 
488  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
489 
490  double errUp = 0;
491  double errDn = 0;
492 
493  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
494  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
495  double lo = y-dns[systIdx]->GetBinContent(binIdx);
496 
497  if(lo <= 0 && hi <= 0) std::swap(lo, hi);
498 
499  errUp += hi*hi;
500  errDn += lo*lo;
501 
502  // TODO: what happens if they're both high or both low?
503  } // end for systIdx
504  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
505  } // end for i
506 
507  nom->Add(hBKG);
508  nom->SetLineColor(col);
509  nom->GetYaxis()->SetTitle(hBKG->GetYaxis()->GetTitle());
510  nom->GetYaxis()->SetTitleSize(0.06);
511  nom->GetYaxis()->SetTitleOffset(1.00);
512  nom->GetYaxis()->CenterTitle();
513  nom->GetYaxis()->SetLabelSize(0.06);
514  nom->GetXaxis()->SetTitle(hBKG->GetXaxis()->GetTitle());
515  nom->GetXaxis()->SetTitleOffset(1.00);
516  nom->GetXaxis()->CenterTitle();
517  nom->GetXaxis()->SetLabelSize(0.06);
518  nom->Draw("hist");
519 
520  g->SetFillColor(errCol);
521  g->Draw("e2 same");
522  g->GetYaxis()->SetRangeUser(0, ymax);
523  nom->GetYaxis()->SetRangeUser(0, ymax);
524 
525  nom->Draw("hist same");
526 
527  for(TH1* up: ups) delete up;
528  for(TH1* dn: dns) delete dn;
529 }
530 
531 //----------------------------------------------------------------------
532 
533 void RemoveBadPoints(TGraphAsymmErrors* gr, TH1* pred)
534 {
535  int n = gr->GetN();
536 
537  double x, y;
538 
539  for(int i = n-1; i>=0; i--)
540  {
541  gr->GetPoint(i,x,y);
542  if(pred->GetBinContent(pred->FindBin(x)) < 1.e-20) // zero is sometimes problematic
543  gr->RemovePoint(i);
544  }
545 }
Near Detector underground.
Definition: SREnums.h:10
const Var kHadNHit([](const caf::SRProxy *sr){unsigned int nought=0;if(sr->trk.kalman.ntracks< 1) return nought;return sr->slc.nhit-sr->trk.kalman.tracks[0].nhit;})
Definition: NumuVars.h:41
SREnergyBranchProxy energy
Definition: SRProxy.h:2251
TSpline3 lo("lo", xlo, ylo, 12,"0")
const Var kNKalman
Definition: NumuVars.cxx:517
Implements systematic errors by interpolation between shifted templates.
Far Detector at Ash River.
Definition: SREnums.h:11
osc::OscCalculatorDumb calc
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
const Var kNcellsFromEdge
Definition: NumuVars.cxx:537
Oscillation analysis framework, runs over CAF files outside of ART.
const Var kReMIdScatLLH
Definition: NumuVars.cxx:532
std::vector< SystGroupDef > systs
Definition: syst_header.h:384
void RemoveBadPoints(TGraphAsymmErrors *gr, TH1 *pred)
General interface to oscillation calculators.
Definition: FwdDeclare.h:15
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:553
std::string label
Definition: CutFlow_Data.C:28
Optimized version of OscCalculatorPMNS.
Definition: StanTypedefs.h:28
const Var kSlcMaxY([](const caf::SRProxy *sr){return sr->slc.boxmax.Y()/100.;})
Definition: NumuVars.h:74
const Var kNumuContPID
Definition: NumuVars.cxx:530
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Var kNumuCosRejAngleKal
Definition: NumuVars.cxx:524
const char * p
Definition: xmltok.h:285
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:17
T sqrt(T number)
Definition: d0nt_math.hpp:156
const Var kTrkStartY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Y()/100;})
Definition: NumuVars.h:32
virtual Spectrum PredictUnoscillated() const
Definition: IPrediction.cxx:82
const Var kSliceTime([](const caf::SRProxy *sr){return sr->slc.meantime/1000;})
Definition: NumuVars.h:18
Proxy< float > recotrkcchadE
Definition: SRProxy.h:1736
const Var kCosBackDist([](const caf::SRProxy *sr){return sr->sel.contain.cosbakdist/100.;})
Definition: NumuVars.h:51
string directory
projection from multiple chains
var
const Var kSliceDuration([](const caf::SRProxy *sr){return(sr->slc.endtime-sr->slc.starttime);})
Definition: NumuVars.h:19
const Var kCosFwdDist([](const caf::SRProxy *sr){return sr->sel.contain.cosfwddist/100.;})
Definition: NumuVars.h:52
void ResetOscCalcToDefault(osc::IOscCalculatorAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
Proxy for StandardRecord.
Definition: SRProxy.h:2237
const Cut kNumuFD
Definition: NumuCuts.h:52
void SetSpillCut(const SpillCut &cut)
TGraphAsymmErrors * GraphWithPoissonErrors(const TH1 *h, bool noErrorsXaxis, bool drawEmptyBins)
Calculate statistical errors appropriate for small Poisson numbers.
Definition: Plots.cxx:842
void SetTh23(const T &th23) override
const Var kDirZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].dir.Z();})
Definition: NumuVars.h:23
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
Definition: NumuVars.h:45
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
Generates extrapolated Numu predictions.
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
void SetDmsq32(const T &dmsq32) override
const Var kNumuHadCalE
Definition: NumuVars.cxx:515
virtual Spectrum PredictSyst(osc::IOscCalculator *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:98
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:170
const char * label
const Binning kZBins
const XML_Char const XML_Char * data
Definition: expat.h:268
void draw_ratio(TCanvas *c1, TH1 *MC, TH1 *Unosc, TH1 *data, TGraphAsymmErrors *grData)
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
Double_t ymax
Definition: plot.C:25
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
Charged-current interactions.
Definition: IPrediction.h:39
const Var kTimingSidebandWeight([](const caf::SRProxy *sr){if(sr->spill.run<=util::kLastBadTimingRun) return util::kBeamWindowBadPeriodMicroSec/util::kTimingSidebandBadPeriodMicroSec;return util::kBeamWindowMicroSec/util::kTimingSidebandMicroSec;})
Definition: TimingCuts.h:29
TSpline3 hi("hi", xhi, yhi, 18,"0")
const Cut kCCEshiftBelow5({"energy.numusimp.trkccE","sel.remid.bestidx","energy.mutrkE.E","hdr.ismc"}, [](const caf::StandardRecord *sr){float muE=sr->energy.mutrkE[sr->sel.remid.bestidx].E;float hadE=sr->energy.numusimp.trkccE-muE;return((sr->hdr.ismc &&sr->energy.numusimp.trkccE< 5.0) ||(!sr->hdr.ismc &&muE+1.14f *hadE< 5.0));})
====================================================================== ///
Definition: CutFlow_Data.C:27
const Var kKalmanFwdDist([](const caf::SRProxy *sr){return sr->sel.contain.kalfwddist/100.;})
Definition: NumuVars.h:54
const Var kRemID
PID
Definition: Vars.cxx:79
SRHeaderProxy hdr
Definition: SRProxy.h:2245
Int_t col[ntarg]
Definition: Style.C:29
const int kNumPlots
Definition: GetSpectra.h:1
const Var kNHit
Definition: Vars.cxx:69
const Var kTrkStartZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Z()/100;})
Definition: NumuVars.h:33
const Cut kNumuND
Definition: NumuCuts.h:54
const Var kCCE
Definition: Vars.h:90
const Binning kEnergyBinning
#define pot
virtual void Go() override
Load all the registered spectra.
void fa_fd_data_mc_systs()
loader
Definition: demo0.py:10
const Var kKalmanBackDist([](const caf::SRProxy *sr){return sr->sel.contain.kalbakdist/100.;})
Definition: NumuVars.h:55
TLatex * prelim
Definition: Xsec_final.C:133
void Preliminary()
const Var kTrkStartX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.X()/100;})
Definition: NumuVars.h:31
const SystShifts kNoShift
Definition: SystShifts.h:112
const Var kTrkEndZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.Z()/100;})
Definition: NumuVars.h:37
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
SRNumuEnergyProxy numu
Definition: SRProxy.h:1916
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:94
const std::vector< PlotDef > plots
const Binning bins
Definition: NumuCC_CPiBin.h:8
std::string fname
Definition: CutFlow_Data.C:29
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::vector< const ISyst * > getAllNumuSystsSA(bool useGENIENC)
Get a vector of all the numu group systs.
Definition: NumuSysts.cxx:216
Proxy< float > trkccE
Definition: SRProxy.h:1729
const Var kTrkEndY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.Y()/100;})
Definition: NumuVars.h:36
const Var kTrkEndX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.X()/100;})
Definition: NumuVars.h:35
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Spectrum Predict(osc::IOscCalculator *calc) const override
const Var kNplanesToBack
Definition: NumuVars.cxx:536
std::vector< Loaders * > loaders
Definition: syst_header.h:385
string syst
Definition: plotSysts.py:176
Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
const Cut kInTimingSideband([](const caf::SRProxy *sr){const double t=sr->slc.meantime;if(t< 1000 *util::kMinTimingSidebandBeforeMicroSec|| t > 1000 *util::kMaxTimingSidebandAfterMicroSec) return false;if(t > 1000 *util::kMaxTimingSidebandBeforeMicroSec && t< 1000 *util::kMinTimingSidebandAfterMicroSec) return false;if(sr->spill.run<=util::kLastBadTimingRun && t > 1000 *util::kMaxTimingSidebandBeforeShiftedWindowMicroSec && t< 1000 *util::kMinTimingSidebandAfterShiftedWindowMicroSec) return false; return true;}, [](const caf::SRSpillProxy *spill){if(spill->run<=util::kLastBadTimingRun) return 1e-6 *util::kTimingSidebandBadPeriodMicroSec;else return 1e-6 *util::kTimingSidebandMicroSec;}, [](const caf::SRSpillProxy *spill){return 0;})
Definition: TimingCuts.h:13
Neutral-current interactions.
Definition: IPrediction.h:40
Proxy< bool > ismc
Definition: SRProxy.h:1007
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
c1
Definition: demo5.py:24
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const Var kReMIdDEDxLLH
Definition: NumuVars.cxx:533
const Var kNplanesToFront
Definition: NumuVars.cxx:535
Proxy< float > recomuonE
Definition: SRProxy.h:1733
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:28
Binning bins
Definition: CutFlow_Data.C:30
All neutrinos, any flavor.
Definition: IPrediction.h:26
Var var
Definition: CutFlow_Data.C:31
Float_t e
Definition: plot.C:35
void myPlotWithSystErrorBand(IPrediction *pred, Spectrum *bkg, const std::vector< const ISyst * > &systs, osc::IOscCalculator *calc, double pot, int col=-1, int errCol=-1, float ymax=60.0)
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
Definition: Loaders.cxx:25
Float_t w
Definition: plot.C:20
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:70
const Binning kXYBins
static constexpr Double_t sr
Definition: Munits.h:164
const Var kTrkNhits([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 65535;return int(sr->trk.kalman.tracks[0].nhit);})
Definition: NumuVars.h:39