selection_story_plots.C
Go to the documentation of this file.
1 
2 #include <fstream>
8 #include "CAFAna/Core/Ratio.h"
11 #include "CAFAna/Core/Spectrum.h"
12 #include "CAFAna/Vars/Vars.h"
13 #include "CAFAna/Analysis/Style.h"
14 #include "CAFAna/Analysis/Plots.h"
15 #include "CAFAna/Analysis/Calcs.h"
20 //#include "./vars.h"
21 
22 using namespace ana;
23 
24 #include "OscLib/IOscCalc.h"
25 #include "OscLib/OscCalcDumb.h"
26 #include "OscLib/OscCalcPMNSOpt.h"
27 
28 #include "TLine.h"
29 #include "TArrow.h"
30 #include "TCanvas.h"
31 #include "TColor.h"
32 #include "TFile.h"
33 #include "TH1.h"
34 #include "TGraphAsymmErrors.h"
35 #include "TLegend.h"
36 #include "TPad.h"
37 #include "TLatex.h"
38 
39 #include <iostream>
40 #include <iomanip>
41 //std::ofstream tables;
42 
43 TLegend* Legend(double x0, double y0, double x1, double y1, bool useData, bool smartLeg=false, bool mc=false)
44 { TLegend* leg;
45  if (!smartLeg) {leg = new TLegend(x0, y0, x1, y1);}
46  if (smartLeg) {leg = AutoPlaceLegend(x1-x0,y1-y0,0.6); }
47  leg->SetFillStyle(0);
48 
49  TH1* dummy = new TH1F("", "", 1, 0, 1);
50  TH1* dummyFill = new TH1F("", "", 1, 0, 1);
51  if(!mc){
52  dummyFill->SetFillColor(kTotalMCErrorBandColor);
53  dummy->SetMarkerStyle(kFullCircle);
54  if( useData) leg->AddEntry(dummy->Clone(), "FD data", "lep");
55  if(!useData) leg->AddEntry(dummy->Clone(), "FD mock data", "lep");
56  dummy->SetLineColor(kViolet-1);
57  leg->AddEntry(dummy->Clone(), "Total Pred. (Best fit)", "l");
58  dummy->SetLineColor(kGray+1);
59  FillWithDimColor(dummy);
60  leg->AddEntry(dummy->Clone(), "Total Background", "bf");
61  dummy->SetLineColor(kAzure+1);
62  FillWithDimColor(dummy);
63  leg->AddEntry(dummy->Clone(), "Cosmic Background", "bf");}
64 
65  if(mc){
66  dummy->SetLineColor(kPink-3);
67  FillWithDimColor(dummy);
68  leg->AddEntry(dummy->Clone(), "signal #nu_{e} CC", "bf");
69  dummy->SetLineColor(kMagenta+3);
70  FillWithDimColor(dummy);
71  leg->AddEntry(dummy->Clone(), "beam #nu_{e} CC", "bf");
72  dummy->SetLineColor(kYellow-9);
73  FillWithDimColor(dummy);
74  leg->AddEntry(dummy->Clone(), "NC", "bf");
75  TColor *color4 = new TColor(232, 0.85, 0.75, 0.58);
76  dummy->SetLineColor(kRed+1);
77  FillWithDimColor(dummy);
78  leg->AddEntry(dummy->Clone(), "#nu_{#mu} CC", "bf");
79  dummy->SetLineColor(kGreen+3);
80  FillWithDimColor(dummy);
81  leg->AddEntry(dummy->Clone(), "#nu_{#tau} CC", "bf");
82  dummy->SetLineColor(kMagenta);
83  FillWithDimColor(dummy);
84  leg->AddEntry(dummy->Clone(), "#bar{#nu_{e}} CC", "bf");
85  dummy->SetLineColor(kAzure+1);
86  FillWithDimColor(dummy);
87  leg->AddEntry(dummy->Clone(), "cosmic", "bf");
88  }
89 
90  return leg;
91 }
92 
94 {
95  TLatex* prelim = new TLatex(.935, .95, "NOvA Preliminary");
96  prelim->SetTextColor(kBlue);
97  prelim->SetNDC();
98  prelim->SetTextSize(2/30.);
99  prelim->SetTextAlign(32);
100  prelim->Draw();
101 }
102 
103 
104 void selection_story_plots(std::string fnameMC, std::string fnameCo, std::string fnameDa, std::string fnameRo, bool useData=false)
105 {
106 
107  //tables.open(((TString)"plots/tables_"+(useData?"data":"fake") +".txt").Data(), std::ofstream::out);
108 
112  PredictionExtrap* extrapCombo[kNumSels][kNumVars];
114 
115  //osc::OscCalcDumb calc; //don't forget about it
116  auto calc = DefaultOscCalc();
117  calc->SetdCP(1.21*M_PI);
118  calc->SetTh23(asin(sqrt(0.556)));
119  calc->SetDmsq32(2.44e-3);
120  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
121  const char* selName = selNames[selIdx].c_str();
122  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
123  const char* varName = defs[varIdx].name.c_str();
124  spectMC[selIdx][varIdx] = LoadFromFile<PredictionNoExtrap>(fnameMC, TString::Format("%s/%s", selName, varName).Data()).release();
125  cosmic[selIdx][varIdx] = LoadFromFile<Spectrum>(fnameCo, TString::Format("%s/cosm_%s", selName, varName).Data()).release();
126  data[selIdx][varIdx] = LoadFromFile<Spectrum>(fnameDa, TString::Format("%s/spect_%s", selName, varName).Data()).release();
127  rock[selIdx][varIdx] = LoadFromFile<PredictionNoExtrap>(fnameRo, TString::Format("%s/rock_%s", selName, varName).Data()).release();
128  extrapCombo[selIdx][varIdx]=LoadFromFile<PredictionExtrap>(fnameMC, TString::Format("%s/%s_extrap_combo", selName, varName).Data()).release();
129  }
130  }
131  TH1* hNue; //signal
132  TH1* hMC; //total MC
133  TH1* hTotbkg; //total sum of all bkg
134  TH1* hbkg; // all bkg except cosmics
135 
136  TH1* hNC; //1
137  TH1* hCC; //2
138  TH1* hBeam; //3
139  TH1* hTau; //4
140  TH1* hCos; //5
141  TH1* hRock; //6 rock events
142  TH1* hNue_bg; //7 wrong sign nue
143 
144  TH1* hData; //data
145  TH1* hFakeData;
146 
147  int start = 0; int end = kNumSels;
148  bool mc = 0;
149  bool sideband = 0;
150  if(sideband) {end = 6;}
151  if(!sideband) {start = 6;}
152  bool useExtrap = 1;
153 
154  for(int selIdx = 0; selIdx < end; ++selIdx){
155  const std::string selName = selNames[selIdx];
156  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
157  const char* varName = defs[varIdx].name.c_str();
158  double nonsScale = 1./10.45;
159  double swapScale = 1./12.91;
160  IPrediction* spectMC2 = spectMC[selIdx][varIdx]; // just FD MC
161  auto data1 = data[selIdx][varIdx];
162  auto cosmic1 = cosmic[selIdx][varIdx];
163 
164  if(useExtrap) {
165  /*if(varIdx==0){
166  auto noextrap = spectMC[selIdx][varIdx];
167  spectMC2 = new PredictionExtendToPeripheral(extrapCombo[selIdx][varIdx],noextrap);
168  data1 = data[selIdx][21];
169  cosmic1 = cosmic[selIdx][21];
170  }*/
171  auto spectMC2 = extrapCombo[selIdx][varIdx];
172  }
173  auto spectMC1 = new PredictionAddRock(spectMC2, rock[selIdx][varIdx], nonsScale, swapScale); //with rock events
174  const double POT=data[selIdx][varIdx]->POT();
175  hNue = spectMC1->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu).ToTH1(POT); //signal
176  hNue_bg = spectMC1->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu).ToTH1(POT); //wrong sign bkg (1)
177  hTau = spectMC1->PredictComponent(calc, Flavors::kAllNuTau, Current::kCC, Sign::kBoth).ToTH1(POT); // tau CC bkg (2)
178  hBeam = spectMC1->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).ToTH1(POT); // beam bkg (3)
179  hNC = spectMC1->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(POT); // NC bkg (4)
180  hCC = spectMC1->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(POT); // numu CC bkg (5)
181  hMC = spectMC1->Predict(calc).ToTH1(POT);
182 
183  bool useScales = !useExtrap && true; //flat scaling to mimic the extrap
184  if(useScales){
185  std::cout << "Rescaling MC components to Extrap/No Extrap TODO";
186  hMC->Scale(57.5/57.1);
187  hNue->Scale(42.2/43.1);
188  hBeam->Scale(7.09/7.03);
189  hNC->Scale(6.62/5.52);
190  hCC->Scale(1.1/1.0);
191  hTau->Scale(0.446/0.446);
192  hNue_bg->Scale(1.0/1.0);
193  }
194  hTotbkg = (TH1*)hNue_bg->Clone(UniqueName().c_str());
195  hTotbkg->Add(hTau);
196  hTotbkg->Add(hBeam);
197  hTotbkg->Add(hNC);
198  hTotbkg->Add(hCC);
199 
200  hTotbkg->SetLineColor(kGray+1);
201  FillWithDimColor(hTotbkg);
202  hMC->SetLineColor(kViolet-1);
203  hRock = spectMC1->Predict(calc).ToTH1(POT);
204  TH1* hJustMC = spectMC2->Predict(calc).ToTH1(POT);
205  hRock->Add(hJustMC, -1); // just rock events (bkg?)
206  double Livetime = data[selIdx][varIdx]->Livetime();
207  cout<<Livetime<<endl;
208  hCos = cosmic1->ToTH1(Livetime, kGray+1, kSolid, ana::EExposureType::kLivetime); //cosmic bkg (6)
209  cout<<hCos->Integral()<<endl;
210  hbkg = (TH1D*) hTotbkg->Clone(UniqueName().c_str());
211  hTotbkg->Add(hCos);
212  hMC->Add(hCos);
213  hCos->SetLineColor(kAzure+1);
214  FillWithDimColor(hCos);
215  if(!mc) {hData = data1->ToTH1(Livetime, kBlack, kSolid, ana::EExposureType::kLivetime);
216  if(!useData) {
217  std::cout << "\n\nWARNING useData=false; plotting mock data\n\n";
218  hData = spectMC1->Predict(calc).ToTH1(POT,kBlack,kSolid);
219  hData->Add(hCos);
220  Spectrum temp (hData,POT,0);
221  hData = temp.MockData(POT).ToTH1(POT);
222  }}
223 
224 /* for (auto &h:{hData, hMC, hNue, hbkg, hCos}){
225  std::cout<<" 1 bin " << h->Integral(1,9) << std::endl;
226  std::cout<<" 2 bin " << h->Integral(10,18) << std::endl;
227  std::cout<<" 3 bin " << h->Integral(19,27) << std::endl;
228  std::cout<<" periph " << h->Integral(28, 36) << std::endl;
229  std::cout<<" all " << h->Integral() << std::endl;
230  }
231 */
232  //HACK HACK HACK
233 /* if((varIdx==2 || varIdx==20 )&& selIdx==4){
234  for(auto &h:{hNue,hbkg,hTotbkg,hCos,hMC}){
235  if(varIdx==2) h->Rebin(10);
236  if(varIdx==20) h->Rebin(15);
237  }
238  if(!mc) {if(varIdx==2) hData->Rebin(10);
239  if(varIdx==20) hData->Rebin(15);}
240  }
241 */
242  if(!mc) {hData->SetMarkerStyle(kFullCircle);}
243 
244  int rebin = 1;
245  if (varIdx == 1 && selIdx == 0 || varIdx == 1 && selIdx == 1 || varIdx == 1 && selIdx == 2) {rebin = 3;}
246  if (varIdx == 1 && selIdx == 3 || varIdx == 1 && selIdx == 4 || varIdx == 1 && selIdx == 5) {rebin = 5;}
247  if (varIdx==20) {rebin = 3;}
248  if (varIdx==19 || varIdx==6 || varIdx==7 || varIdx==8) {rebin = 4;}
249  if (varIdx==11 || varIdx==12 || varIdx==13 || varIdx==14 || varIdx==15 || varIdx==16 || varIdx==17) {rebin = 10;}
250  if (!mc) {hData->Rebin(rebin);}
251  for (auto &h:{hMC, hNue, hbkg, hTotbkg, hNue_bg, hNC, hCC, hBeam, hTau, hRock, hCos}){
252  h->Rebin(rebin);
253  }
254  // if(varIdx==0) {cout<<"MC "<<hMC->Integral()<< " signal "<<hNue->Integral()<<" total beam "<<hbkg->Integral()<<" cosmic "<<hCos->Integral()<<endl;}
255  cout<<"MC "<<hMC->Integral()<< " signal "<<hNue->Integral()<<" total beam "<<hbkg->Integral()<<" cosmic "<<hCos->Integral()<<endl;
256  //Print results
257 /*
258  std::cout <<"============== Information =============="<<std::endl;
259  std::cout << "\nSelection " << selName << "; Variable " << varName
260  << (useExtrap?"; Extrap.":";No Extrap.\n")
261  << "Compare data to prediction - dumb oscillations"
262  << std::endl;
263  std::cout << "Components including rock\n"
264  << "Nue Antinue Beam NC Numu NuTau Total Rock-only\n";
265  for (auto &h:{hNue,hNue_bg,hBeam,hNC,hCC,hTau,hbkg,hRock}){
266  std::cout << h->Integral() << "\t";
267  }
268  TString ltxcommand = selName+varName; ltxcommand.ReplaceAll("_","");
269  TString shname = selName+" " + varName; shname.ReplaceAll("_"," ");
270  shname.ReplaceAll (" custom","");
271  if(varIdx==2 || varIdx==3 || varIdx==4 || varIdx==20 || varIdx==0 || varIdx==1 || varIdx==21 || varIdx==22){
272  tables << "\n\\newcommand{\\" << ltxcommand << (useData?"data":"fake")
273  << "}\n{\\centerline{";
274  tables << "\n\\begin{tabular}{p{7em}|rrr|rr}\n";
275  tables << "& \t Signal \t & \t " << "Bkg. \t & \t" << "Cosmic \t & \t"
276  << "Total \t & \t" << (useData?"Data":"FakeData") << "\t \\\\ \\hline \n";
277  tables <<shname;
278  for (auto &h:{hNue,hbkg,hCos,hMC,hData}){
279  tables << std::setprecision(3);
280  tables << "\t & \t" <<h->Integral();
281  }
282  tables << "\\end{tabular}\n}\n}\n"<<std::endl;
283  }
284  std::cout <<"============== And that's it =============="<<std::endl;
285 */
286  hMC->SetMaximum(1.1*hMC->GetMaximum());
287  if(!mc) {int n = hData->GetNbinsX();
288  for (int i=1; i<= n; i++) {if(hData->GetBinContent(i) <0){
289  for (auto &h:{hMC,hData,hTotbkg,hCos}) h->SetBinContent(i,0);}}
290  double mult = 1.4;
291  if (varIdx==0) {mult = 1.11;}
292  if ( varIdx==7 || varIdx==8) {mult = 1.6;}
293  if (varIdx==6) {mult = 1.8;}
294  if (varIdx == 1 && selIdx == 2) {mult = 1.0;}
295  hMC->SetMaximum(mult*max(hMC->GetMaximum(),(hData->GetMaximum()+sqrt(hData->GetMaximum()))));}
296  hMC->GetXaxis()->CenterTitle();
297  if (varIdx == 1 && selIdx == 0 || varIdx == 1 && selIdx == 1 || varIdx == 1 && selIdx == 2) {hMC->GetXaxis()->SetRange(50.0/rebin+2,100.0/rebin);}
298  if (varIdx == 1 && selIdx == 3 || varIdx == 1 && selIdx == 4 || varIdx == 1 && selIdx == 5) {hMC->GetXaxis()->SetRange(50.0/rebin+1,100.0/rebin);}
299  if (varIdx==5) {hMC->GetXaxis()->SetRange(1,20);}
300  if (selIdx==6 && varIdx==2) {hMC->GetXaxis()->SetRange(50.0/rebin+1,100.0/rebin);}
301  hMC->GetYaxis()->CenterTitle();
302  hMC->GetYaxis()->SetTitleOffset(1.15);
303  hMC->GetYaxis()->SetTitle("Events / 8.85 #times 10^{20} POT-equiv");
304  if(varIdx == 19){hMC->GetXaxis()->SetTitle("p_{T}/p");}
305 
306  TCanvas *c =new TCanvas("c","c",1000,800);
307  c->SetLeftMargin(0.15);
308  c->SetRightMargin(0.05);
309 
310  gPad->SetFillStyle(0);
311  if(mc){
312  hMC->SetLineColor(kPink-3);
313  FillWithDimColor(hMC);
314  hMC->Draw("hist");
315  hTotbkg->SetLineColor(kMagenta+3);
316  FillWithDimColor(hTotbkg);
317  hTotbkg->Draw("hist same");//beam
318  TH1* wobeam = (TH1*)hTotbkg->Clone(UniqueName().c_str());
319  wobeam->SetLineColor(kYellow-9);
320  FillWithDimColor(wobeam);
321  wobeam->Add(hBeam,-1);
322  wobeam->Draw("hist same");//nc
323  TH1* wonc = (TH1*)wobeam->Clone(UniqueName().c_str());
324  wonc->SetLineColor(kRed+1);
325  FillWithDimColor(wonc);
326  wonc->Add(hNC,-1);
327  wonc->Draw("hist same");//cc
328  TH1* wocc = (TH1*)wonc->Clone(UniqueName().c_str());
329  wocc->SetLineColor(kGreen+3);
330  FillWithDimColor(wocc);
331  wocc->Add(hCC,-1);
332  wocc->Draw("hist same");//tau
333  TH1* wotau = (TH1*)wocc->Clone(UniqueName().c_str());
334  wotau->SetLineColor(kMagenta);
335  FillWithDimColor(wotau);
336  wotau->Add(hTau,-1);
337  wotau->Draw("hist same");//nue wrong sign
338  TH1* wown = (TH1*)wotau->Clone(UniqueName().c_str());
339  wown->SetLineColor(kGray);
340  FillWithDimColor(wown);
341  wown->Add(hNue_bg,-1);
342  wown->Draw("hist same");//rock
343  TH1* worock = (TH1*)wown->Clone(UniqueName().c_str());
344  worock->SetLineColor(kAzure+1);
345  FillWithDimColor(worock);
346  wown->Add(hRock,-1);
347  worock->Draw("hist same");//cosmic
348  }
349  if(!mc) {
350  hMC->Draw("hist ][");
351  hTotbkg->Draw("hist same");
352  hCos->Draw("hist same");
353  hMC->Draw("hist ][ same");
354  //TGraphAsymmErrors* gr = GraphWithPoissonErrors(hData, false, false);
355  // gr->SetMarkerStyle(kFullCircle);
356  // gr->SetLineWidth(2);
357  //gr->Draw("ep same");
358  }
359 
360  Preliminary();
361 
362  double legendxlow=0.2;
363  double legendxhigh=0.5;
364  double legendylow=0.5;
365  double legendyhigh=0.7;
366  bool marker = true;
367  if (varIdx==1) {marker = false; legendxlow=0.2; legendxhigh=0.5; legendylow=0.7; legendyhigh=0.88;}
368  if (varIdx==0) { marker = false; legendxlow=0.2; legendxhigh=0.5; legendylow=0.6; legendyhigh=0.8; }
369  if (selIdx == 5 && varIdx == 4 || selIdx == 4 && varIdx == 20 || selIdx == 4 && varIdx == 2 || selIdx == 0 && varIdx == 3 || selIdx == 0 && varIdx == 4 ||
370  selIdx == 2 && varIdx == 3) { marker = false; legendxlow=0.17; legendxhigh=0.465; legendylow=0.6; legendyhigh=0.8; }
371  if (varIdx==21) { marker = false; legendxlow=0.2; legendxhigh=0.5; legendylow=0.5; legendyhigh=0.7; }
372  if (mc) {legendylow=0.5; legendyhigh=0.8;}
373  auto leg = Legend(legendxlow, legendylow, legendxhigh, legendyhigh, useData, marker, mc);
374  leg->SetTextSize(hMC->GetXaxis()->GetLabelSize());
375  leg->Draw();
376 
377  if(varIdx == 0){
378  Nue2017FourBinDivisions(kBlack, 3);
379  Nue2017FourBinLabels(0.98, 0.031, kGray+3, true);
380  Nue2017FourBinAxis(hMC,true, true);
381  }
382 
383  gPad->Update();
384 
385  if (selName=="core_basic") {auto *tt = new TText(gPad->GetUxmin(), gPad->GetUymax() * 1.05,"Basic cuts");
386  int n = hData->GetNbinsX();
387  for (int i=1; i<= n; i++) {cout<<hData->GetBinContent(i); if(i>25) {hData->SetBinContent(i,-1); cout<<"in loop"<<endl;} cout<<i<<" hi "<<hData->GetBinContent(i)<<endl;}
388  tt->Draw();
389  TLine *l=new TLine(0.87,0,0.87,gPad->GetUymax());
390  l->SetLineStyle(3);
391  l->SetLineColor(kGray+2);
392  l->SetLineWidth(2);
393  l->Draw();
394  l=new TLine(0.95,0,0.95,gPad->GetUymax());
395  l->SetLineStyle(3);
396  l->SetLineColor(kGray+2);
397  l->SetLineWidth(2);
398  l->Draw();
399  l=new TLine(0.75,0,0.75,gPad->GetUymax());
400  // l->SetLineStyle(2);
401  l->SetLineWidth(2);
402  l->SetLineColor(kGray+2);
403  l->Draw();}
404  if (selName=="core_presel") {auto *tt = new TText(gPad->GetUxmin(), gPad->GetUymax() * 1.05,"Preselection cut, Core");
405  int n = hData->GetNbinsX();
406  for (int i=1; i<= n; i++) {cout<<hData->GetBinContent(i); if(i>25) {hData->SetBinContent(i,-1); cout<<"in loop"<<endl;} cout<<i<<" hi "<<hData->GetBinContent(i)<<endl;}
407  tt->Draw();
408  TLine *l=new TLine(0.87,0,0.87,gPad->GetUymax());
409  l->SetLineStyle(3);
410  l->SetLineColor(kGray+2);
411  l->SetLineWidth(2);
412  l->Draw();
413  l=new TLine(0.95,0,0.95,gPad->GetUymax());
414  l->SetLineStyle(3);
415  l->SetLineColor(kGray+2);
416  l->SetLineWidth(2);
417  l->Draw();
418  l=new TLine(0.75,0,0.75,gPad->GetUymax());
419  l->SetLineWidth(2);
420  l->SetLineColor(kGray+2);
421  l->Draw();}
422  if (selName=="core_full_presel") {auto *tt = new TText(gPad->GetUxmin(), gPad->GetUymax() * 1.05,"Full Preselection, Core");
423  int n = hData->GetNbinsX();
424  for (int i=1; i<= n; i++) {cout<<hData->GetBinContent(i); if(i>25) {hData->SetBinContent(i,-1); cout<<"in loop"<<endl;} cout<<i<<" hi "<<hData->GetBinContent(i)<<endl;}
425  tt->Draw();
426  TLine *l=new TLine(0.87,0,0.87,gPad->GetUymax());
427  l->SetLineStyle(3);
428  l->SetLineColor(kGray+2);
429  l->SetLineWidth(2);
430  l->Draw();
431  l=new TLine(0.95,0,0.95,gPad->GetUymax());
432  l->SetLineStyle(3);
433  l->SetLineColor(kGray+2);
434  l->SetLineWidth(2);
435  l->Draw();
436  l=new TLine(0.75,0,0.75,gPad->GetUymax());
437  l->SetLineWidth(2);
438  l->SetLineColor(kGray+2);
439  l->Draw();
440  TArrow *ar = new TArrow(0.75,25,0.8,25, 0.04);
441  ar->SetLineWidth(2);
442  ar->SetLineColor(kGray+2);
443  ar->Draw();}
444  if (selName=="perip_presel") {auto *tt = new TText(gPad->GetUxmin(), gPad->GetUymax() * 1.05,"Preselection, Peripheral");
445  int n = hData->GetNbinsX();
446  for (int i=1; i<= n; i++) {cout<<hData->GetBinContent(i); if(i>19) {hData->SetBinContent(i,-1); cout<<"in loop "<<endl;} cout<<" "<<i<<" hi "<<hData->GetBinContent(i)<<endl;}
447  tt->Draw();
448  TLine *l=new TLine(0.95,0,0.95,gPad->GetUymax());
449  l->SetLineWidth(2);
450  l->SetLineColor(kGray+2);
451  l->Draw();}
452  if (selName=="perip_presel_bdt_test_cvn") {auto *tt = new TText(gPad->GetUxmin(), gPad->GetUymax() * 1.05,"BDT cut, Peripheral");
453  int n = hData->GetNbinsX();
454  for (int i=1; i<= n; i++) {cout<<hData->GetBinContent(i); if(i>19) {hData->SetBinContent(i,-1); cout<<"in loop"<<endl;} cout<<i<<" hi "<<hData->GetBinContent(i)<<endl;}
455  tt->Draw();
456  TLine *l=new TLine(0.95,0,0.95,gPad->GetUymax());
457  l->SetLineWidth(2);
458  l->SetLineColor(kGray+2);
459  TArrow *ar = new TArrow(0.95,18,0.98,18, 0.03);
460  ar->SetLineWidth(2);
461  ar->SetLineColor(kGray+2);
462  ar->Draw();
463  l->Draw();
464  cout<<hCos->GetBinContent(20)<<" "<<hNue->GetBinContent(20)<<" "<<hMC->GetBinContent(20)<<" "<<hTotbkg->GetBinContent(20)<<endl;}
465  if (varIdx==0) cout<<"4 bin "<<hCos->Integral(28,36)<<" "<<hNue->Integral(28,36)<<" "<<hMC->Integral(28,36)<<" "<<hTotbkg->Integral(28,36)<<endl;
466  TGraphAsymmErrors* gr = GraphWithPoissonErrors(hData, false, false);
467  gr->SetMarkerStyle(kFullCircle);
468  gr->SetLineWidth(2);
469  gr->Draw("ep same");
470 
471  if (selName=="Peripheral") {auto *tt = new TText(gPad->GetUxmin(), gPad->GetUymax() * 1.05,"Peripheral Sample");
472  tt->Draw();}
473  if (selName=="Core") {auto *tt = new TText(gPad->GetUxmin(), gPad->GetUymax() * 1.05 ,"Core Sample");
474  tt->Draw();}
475  if (selName=="Presel") {auto *tt = new TText(gPad->GetUxmin(), gPad->GetUymax() * 1.05 ,"Preselection cut tier");
476  tt->Draw();}
477  c->Update();
478  gPad->RedrawAxis();
479  /*if (selIdx == 5) { auto *tt = new TText(gPad->GetUxmin(), gPad->GetUymax() * 1.02 ,"High Energy Sideband (Core)"); tt->SetTextSize(0.03); tt->Draw();}
480  if (selIdx == 4) { auto *tt = new TText(gPad->GetUxmin(), gPad->GetUymax() * 1.02 ,"Low PID Sideband (Core)"); tt->SetTextSize(0.03); tt->Draw();}
481  if (selIdx == 0) { auto *tt = new TText(gPad->GetUxmin(), gPad->GetUymax() * 1.01 ,"High Energy, low PID Sideband (Peripheral)"); tt->SetTextSize(0.03); tt->Draw(); }
482  if (selIdx == 2) { auto *tt = new TText(gPad->GetUxmin(), gPad->GetUymax() * 1.02 ,"Low PID Sideband (Peripheral)"); tt->SetTextSize(0.03); tt->Draw();}
483  */
484 
485  if(sideband){
486  ofstream caption(("plots/sideband/" + selName + "_" + varName + ".txt").c_str());
487  caption<<selName<<" cut applied, sideband data/prediction study, with extrapolation from the ND and combo decomposition, 2017 Best Fit oscillation parameters.";
488  c->Print(("plots/sideband/" + selName + "_" + varName + ".pdf").c_str());}
489  if(mc) {
490  ofstream caption(("plots/fdmc/fdmc_" + selName + "_" + varName + ".txt").c_str());
491  caption<<selName<<" cut applied, FDMC with extrapolation from the ND and combo decomposition, 2017 Best Fit oscillation parameters.";
492  c->Print(("plots/fdmc/fdmc_" + selName + "_" + varName + ".pdf").c_str());}
493  if(useExtrap & !mc & !sideband) {
494  ofstream caption(("plots/extrap/decomp_" + selName + "_" + varName + (useData?"_data":"_fake") + ".txt").c_str());
495  caption<<selName<<" cut applied, data/prediction comparisons, with extrapolation from the ND and combo decomposition, 2017 Best Fit oscillation parameters.";
496  c->Print(("plots/extrap/decomp_" + selName + "_" + varName + (useData?"_data":"_fake") + "_blind.pdf").c_str());}
497  // else c->Print(("plots/" + selName + "_" + varName + (useData?"_data":"_fake") + ".pdf").c_str());
498 
499  }//end var
500  }//end selection
501 }
502 
void Nue2017FourBinLabels(const double yNDC, const double textsize, const int color, const bool merged)
void Nue2017FourBinAxis(TH1 *axes, bool drawLabels, bool merged)
const int kNumVars
Definition: vars.h:14
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Antineutrinos-only.
Definition: IPrediction.h:50
Float_t y1[n_points_granero]
Definition: compare.C:5
(&#39; appearance&#39;)
Definition: IPrediction.h:18
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
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
Float_t x1[n_points_granero]
Definition: compare.C:5
(&#39;beam &#39;)
Definition: IPrediction.h:15
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::string name
Definition: NuePlotLists.h:12
TGraphAsymmErrors * GraphWithPoissonErrors(const TH1 *h, bool noErrorsXaxis, bool drawEmptyBins)
Calculate statistical errors appropriate for small Poisson numbers.
Definition: Plots.cxx:910
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
#define M_PI
Definition: SbMath.h:34
Spectrum MockData(double pot, int seed=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:300
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const XML_Char const XML_Char * data
Definition: expat.h:268
Definition: type_traits.h:56
Charged-current interactions.
Definition: IPrediction.h:39
Sum up livetimes from individual cosmic triggers.
const HistDef defs[kNumVars]
Definition: vars.h:15
void selection_story_plots(std::string fnameMC, std::string fnameCo, std::string fnameDa, std::string fnameRo, bool useData=false)
const int kNumSels
Definition: vars.h:43
TLatex * prelim
Definition: Xsec_final.C:133
double POT() const
Definition: Spectrum.h:227
TLegend * Legend(double x0, double y0, double x1, double y1, bool useData, bool smartLeg=false, bool mc=false)
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
void Preliminary()
enum BeamMode kViolet
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:715
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void Nue2017FourBinDivisions(const int color, const int style)
Prediction that just uses FD MC, with no extrapolation.
Take the output of an extrapolation and oscillate it as required.
All neutrinos, any flavor.
Definition: IPrediction.h:26
const std::string selNames[kNumSels]
Definition: vars.h:46
Float_t e
Definition: plot.C:35
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
enum BeamMode kGreen
enum BeamMode kBlue
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
void FillWithDimColor(TH1 *h, bool usealpha, float dim)
void rock(std::string suffix="full")
Definition: rock.C:28
T asin(T number)
Definition: d0nt_math.hpp:60
TH1F * hbkg
Definition: Xdiff_gwt.C:58
enum BeamMode string