getTimePeakPlots.C
Go to the documentation of this file.
1 #pragma once
2 
5 
6 
7 #include "Includes.h"
8 #include "VarsAndCuts.h"
9 #include "CustomFunctions.h"
10 #include "GetSpectra.h"
11 
12 using namespace ana;
13 
14 void getTimePeakPlots(std::string samdef, unsigned int timestamp, unsigned int nDays, std::string naming)
15 {
16 
17  //this should make us more economic and efficient
18  //only data spectra will be made every time
19  bool makePred = 0; // false because the predictions already exist.
20  if(makePred){
21  getPredictions(timestamp, nDays, naming);
22  }
23  bool makeData = 0;
24  if(makeData){
25  getData(timestamp, nDays, naming, samdef);
26  }
27 
28  std::cout << "Getting files ..." << std::endl;
29  std::string mcFile = "predictions_" + naming + ".root";
30  std::string dtFile = "dataspectra_" + naming + ".root";
31 
32 
33 
34  TFile* mcPred = new TFile(mcFile.c_str());
35  TFile* dtSpec = new TFile(dtFile.c_str());
36  std::cout << "Setting osc pars" << std::endl;
37  osc::OscCalcPMNSOpt calcMin;
38  osc::OscCalcPMNSOpt calcMax;
39  osc::OscCalcPMNSOpt calcDef;
43 
44  std::cout << "PASSED" << std::endl;
45  std::cout << " " << std::endl;
46 
47  calcMin.SetTh13(asin(sqrt(.084-0.008))/2); // (0.084) 0.076 to 0.092
48  calcMin.SetTh12(asin(sqrt(.846+0.021))/2); // (0.846) 0.825 to 0.867
49  calcMin.SetTh23(asin(sqrt(0.437))); // deltam2 > 0 -> sin2theta23 = 0.437 (PDG2016)
50  calcMin.SetDmsq21(7.71e-5); // (7.53) 7.35 to 7.71
51  calcMin.SetDmsq32(2.46e-3); // (2.40) 2.34 to 2.46
52 
53  calcMax.SetTh13(asin(sqrt(.084+0.008))/2); // (0.084) 0.076 to 0.092
54  calcMax.SetTh12(asin(sqrt(.846+0.021))/2); // (0.846) 0.825 to 0.867
55  calcMax.SetTh23(asin(sqrt(0.569))); // deltam2 < 0 -> sin2theta23 = 0.569 (PDG2016)
56  calcMax.SetDmsq21(7.35e-5); // (7.53) 7.35 to 7.71
57  calcMax.SetDmsq32(2.34e-3); // (2.40) 2.34 to 2.46
58 
59 
60  std::cout << "Defining plots ..." << std::endl;
62  {
63  // {"Slice time [#mus]", "tus", Binning::Simple(35,13,433), ana::kSliceTime},
64  {"Slice time [#mus]", "tus", Binning::Simple(45,1,541), ana::kSliceTime},
65  {"Run", "run", Binning::Simple(15000, 15000, 30000), SIMPLEVAR(hdr.run)},
66  {"Time", "time", Binning::Simple(nDays<3 ? 24*nDays : nDays,timestamp - 86400*nDays,timestamp), SIMPLEVAR(hdr.unixtime)},
67  {"Max y [m]", "maxy", kXYBins, kBoxMaxY},
68  {"Slice time [#mus]", "tighttus", Binning::Simple(360,205,241), ana::kSliceTime},
69  {"No. Kalman tracks", "nkal", Binning::Simple(10, 0, 10), SIMPLEVAR(trk.kalman.ntracks)},
70  {"No. Hits", "nhit", Binning::Simple(10, 0, 500), kNHit},
71  {"Kalman angle", "angkal", Binning::Simple(5, 0, 1), SIMPLEVAR(sel.cosrej.anglekal)},
72  {"Z direction", "dirz", Binning::Simple(5, 0.5, 1), ana::kDirZ},
73  {"BDT", "bdt", Binning::Simple(20, 0, 1), ana::kNumuContPID},
74  {"p_{T}/p", "ptp", Binning::Simple(20, 0, 1), ana::kPtP}
75  };
76 
77 
78 
79  std::cout << "Defining spectrums ..." << std::endl;
80  Spectrum* hData[kNumPlots];
81  Spectrum* hBkg[kNumPlots];
83 
84 
85  for(int i = 0; i < kNumPlots; ++i){
86  std::cout << "" << i << std::endl;
87  Plot p = plots[i];
88  std::string dataName = Form("data_%s", (p.fname).c_str());
89  std::string bkgName = Form("bkg_%s", (p.fname).c_str());
90  std::cout << "Data name :" << dataName << " and the Bkg Name :" << bkgName << std::endl;
91  hData[i] = LoadFrom<Spectrum>(dtSpec, dataName).release();
92  hBkg[i] = LoadFrom<Spectrum>(dtSpec, bkgName).release();
93  hPred[i] = LoadFrom<PredictionNoExtrap>(mcPred, (p.fname).c_str() ).release();
94  }
95 
96 
97  std::cout << "Getting miscellanius spectra ..." << std::endl;
98  Spectrum* potTime = LoadFrom<Spectrum>(dtSpec, "potTime").release();
99  Spectrum* massDenom = LoadFrom<Spectrum>(dtSpec, "massDenom").release();
100  Spectrum* massNum = LoadFrom<Spectrum>(dtSpec, "massNum").release();
101  Spectrum* nonSwapDenom = LoadFrom<Spectrum>(mcPred, "nonSwapDenom").release();
102  Spectrum* swapDenom = LoadFrom<Spectrum>(mcPred, "swapDenom").release();
103 
104  const double massFactor = massNum->ToTH1(1)->Integral(0, -1) / massDenom->ToTH1(1)->Integral(0, -1);
105  double signal = 0, mcsig = 0, mcsig_min = 0, mcsig_max = 0, bkg = 0;
106  double mc_numu = 0, mc_nue = 0, mc_nc = 0;
107 
108  std::cout << "Effective mass: " << massFactor*14 << " diblocks" << std::endl;
109 
110  for(int i = 0; i < kNumPlots; ++i){
111  new TCanvas;
112  TH1* hSig = hData[i]->ToTH1(hData[i]->POT());
113  hSig->SetName(("sig_"+plots[i].fname).c_str());
114 
115  TGraphAsymmErrors* grSig = GraphWithPoissonErrors(hSig, true, false);
116 
117  if(i==1 || i==2) // no errors in run or time plots
118  hSig->Draw("hist");
119  else
120  {
121  hSig->SetLineColor(0);
122  hSig->SetLineWidth(1);
123  hSig->SetMarkerColor(0);
124  hSig->GetXaxis()->SetTitle((plots[i].label).c_str());
125  hSig->GetYaxis()->SetRangeUser(0,hSig->GetMaximum()*1.12+4.);
126  hSig->Draw("hist");
127  grSig->SetMarkerStyle(20);
128  grSig->SetMarkerColor(kBlack);
129  grSig->SetLineColor(kBlack);
130  grSig->SetLineWidth(2);
131  grSig->SetName(("grSig_"+plots[i].fname).c_str());
132  grSig->GetXaxis()->SetTitle((plots[i].label).c_str());
133  grSig->Draw("AP");
134  }
135 
136  TH1* hMC = hPred[i]->Predict(&calcDef).ToTH1(hData[i]->POT()*massFactor);
137  TH1* hMC_min = hPred[i]->Predict(&calcMin).ToTH1(hData[i]->POT()*massFactor);
138  hMC_min->SetName(("MC_min_"+plots[i].fname).c_str());
139  hMC_min->SetLineColor(kRed-9);
140  hMC_min->SetLineWidth(1);
141  TH1* hMC_max = hPred[i]->Predict(&calcMax).ToTH1(hData[i]->POT()*massFactor);
142  hMC_max->SetName(("MC_max_"+plots[i].fname).c_str());
143  hMC_max->SetLineColor(kRed-9);
144  hMC_max->SetLineWidth(1);
145 
146  TH1* hMC_numu = hPred[i]->PredictComponent(&calcDef, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(hData[i]->POT()*massFactor);
147  hMC_numu->SetLineColor(kBlue);
148  TH1* hMC_nue = hPred[i]->PredictComponent(&calcDef, Flavors::kAllNuE, Current::kCC, Sign::kBoth).ToTH1(hData[i]->POT()*massFactor);
149  hMC_nue->SetLineColor(kGreen+2);
150  TH1* hMC_NC = hPred[i]->PredictComponent(&calcDef, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(hData[i]->POT()*massFactor);
151  hMC_NC->SetLineColor(kGray);
152 
153  TH1* hBKG = hBkg[i]->ToTH1(hData[i]->POT()*timeWindowFactor);
154  hBKG->SetName(("BKG_"+plots[i].fname).c_str());
155  hBKG->SetLineColor(kBlue);
156 
157  TH1* hExp = hPred[i]->Predict(&calcDef).ToTH1(hData[i]->POT()*massFactor);
158  TH1* hExp_min = hPred[i]->Predict(&calcMin).ToTH1(hData[i]->POT()*massFactor);
159  TH1* hExp_max = hPred[i]->Predict(&calcMax).ToTH1(hData[i]->POT()*massFactor);
160 
161  if(i == 0){
162  // But for which we should print s/b stats
163  signal = hSig->GetBinContent(19);
164  mcsig = hPred[i]->Predict(&calcDef).ToTH1(hData[i]->POT()*massFactor)->Integral(0, -1);
165  mc_numu = hMC_numu->Integral(0, -1);
166  mc_nue = hMC_nue->Integral(0, -1);
167  mc_nc = hMC_NC->Integral(0, -1);
168  mcsig_min = hPred[i]->Predict(&calcMin).ToTH1(hData[i]->POT()*massFactor)->Integral(0, -1);
169  mcsig_max = hPred[i]->Predict(&calcMax).ToTH1(hData[i]->POT()*massFactor)->Integral(0, -1);
170  bkg = (hSig->Integral(1, hSig->GetNbinsX())-signal)/(hSig->GetNbinsX()-1);
171  std::cout << "Signal: " << signal << std::endl;
172  std::cout << "Bkg estimate: " << bkg << std::endl;
173  std::cout << "MC signal min expectation: " << mcsig_min << std::endl;
174  std::cout << "MC signal expectation: " << mcsig << std::endl;
175  std::cout << "MC signal max expectation: " << mcsig_max << std::endl;
176  std::cout << "MC before cuts " << (*nonSwapDenom + *swapDenom).ToTH1(hData[i]->POT())->Integral(0, -1)*massFactor << std::endl;
177  std::cout << "MC expected numu : " << mc_numu << std::endl;
178  std::cout << "MC expected nue : " << mc_nue << std::endl;
179  std::cout << "MC expected NC : " << mc_nc << std::endl;
180 
181  hExp->SetBinContent(19,mcsig+bkg);
182  hExp_min->SetBinContent(19,mcsig_min+bkg);
183  hExp_max->SetBinContent(19,mcsig_max+bkg);
184 
185  TGraph* g = new TGraph;
186  g->SetPoint(0, -1000e3, bkg);
187  g->SetPoint(1, +1000e3, bkg);
188  g->SetLineColor(kBlue);
189  g->SetLineWidth(2);
190  g->SetLineStyle(7);
191  g->Draw("L");
192  }
193  //for time so we can reuse the predictions
194  //first clone BKG so the bins coincide
195  if(i==2){
196  TH1* hHelp = (TH1*)hBKG->Clone();
197  TH1* hHelpMin = (TH1*)hBKG->Clone();
198  TH1* hHelpMax = (TH1*)hBKG->Clone();
199 
200  for(int binId=0; binId<hExp->GetNbinsX(); binId++){
201  hHelp->SetBinContent(binId+1, hExp->GetBinContent(binId+1));
202  hHelpMin->SetBinContent(binId+1, hExp_min->GetBinContent(binId+1));
203  hHelpMax->SetBinContent(binId+1, hExp_max->GetBinContent(binId+1));
204  }
205 
206  hExp->Clear();
207  hExp = (TH1*)hHelp->Clone();
208  hExp->Add(hBKG);
209 
210  hExp_min->Clear();
211  hExp_min = (TH1*)hHelp->Clone();
212  hExp_min->Add(hBKG);
213 
214  hExp_max->Clear();
215  hExp_max = (TH1*)hHelp->Clone();
216  hExp_max->Add(hBKG);
217  }
218  else
219  {
220  hExp->Add(hBKG);
221  hExp_min->Add(hBKG);
222  hExp_max->Add(hBKG);
223  }
224 
225  hExp_max->SetName(("Exp_max_"+plots[i].fname).c_str());
226  hExp_max->SetLineColor(kRed-9);
227  hExp_min->SetName(("Exp_min_"+plots[i].fname).c_str());
228  hExp_max->SetLineColor(kRed-9);
229 
230  TGraph* grShade = ShadeBetweenHistograms(hExp_min, hExp_max);
231  grShade->SetName(("shade_"+plots[i].fname).c_str());
232  grShade->SetFillColor(kRed-9);
233  grShade->SetLineColor(kRed-9);
234  grShade->Draw("lf");
235 
236  double ymax = 1.3* std::max( hSig->GetMaximum(), hExp->GetMaximum());
237  hExp->GetYaxis()->SetRangeUser(0,ymax);
238  gPad->Update();
239  hExp->GetXaxis()->SetTitle((plots[i].label).c_str());
240  hExp->SetName(("Exp_"+plots[i].fname).c_str());
241  hExp->SetLineColor(kRed);
242  hExp->SetLineStyle(1);
243  hExp->Draw("hist same");
244  hMC_numu->Draw("hist same");
245  hMC_nue->Draw("hist same");
246  hMC_NC->Draw("hist same");
247  hExp->Draw("hist same");
248  if(i!=1 && i!=2){
249  hSig->Draw("hist same");
250  grSig->Draw("P same");
251  }
252 
253  TLegend* ll = new TLegend(0.15,0.60,0.45,0.89);
254  ll->AddEntry(grSig,"Data","l");
255  ll->AddEntry(grShade,"Osc. Uncty.","lf");
256  ll->AddEntry(hExp,"Expectation","l");
257  ll->AddEntry(hMC_numu,"#nu_{#mu}","l");
258  ll->AddEntry(hMC_nue,"#nu_{e}","l");
259  ll->AddEntry(hMC_NC,"NC","l");
260  ll->SetFillColor(0);
261  ll->Draw();
262 
263 
264  hSig->SaveAs(("/nova/ana/users/novadq/TimePeak/Plots/" + naming + "/Tools/sig_"+plots[i].fname+".root").c_str());
265  hMC->SaveAs(("/nova/ana/users/novadq/TimePeak/Plots/" + naming + "/Tools/MC_"+plots[i].fname+".root").c_str());
266  hBKG->SaveAs(("/nova/ana/users/novadq/TimePeak/Plots/" + naming + "/Tools/BKG_"+plots[i].fname+".root").c_str());
267 
268  lastUpdated(timestamp);
269  gPad->SaveAs(("/nova/ana/users/novadq/TimePeak/Plots/" + naming + "/Plots/"+plots[i].fname+".png").c_str());
270 
271  if(i == 1){
272  std::cout << "Bad runs:" << std::endl;
273  TH1* h = hData[i]->ToTH1(hData[i]->POT());
274  for(int j = 1; j <= h->GetNbinsX(); ++j){
275  if(h->GetBinContent(j) > 5)
276  std::cout << h->GetXaxis()->GetBinLowEdge(j) << std::endl;
277  }
278  }
279  } // End loop over i
280 
281 
282  TH2* hPOT = potTime->ToTH2(mcsig);
283  hPOT->SetName("hPOT");
284  hPOT->SaveAs(("/nova/ana/users/novadq/TimePeak/Plots/"+naming+"/Tools/hPOT.root").c_str());
285  if(naming == "All")
286  getSigmaPlotsAll(mcsig, timestamp, nDays, naming);
287  else
288  getSigmaPlots(mcsig, timestamp, nDays, naming);
289 
290 }
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
T max(const caf::Proxy< T > &a, T b)
enum BeamMode kRed
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetTh13(const T &th13) override
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
void getData(int sampleCut=9, std::string period="full", std::string anaCut="3A", bool energy3a=true, bool optBin=true)
Definition: getData.C:33
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
double Integral(const Spectrum &s, const double pot, int cvnregion)
const Var kNumuContPID
Definition: NumuVars.cxx:553
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
const Var kSliceTime([](const caf::SRProxy *sr){return sr->slc.meantime/1000;})
Definition: NumuVars.h:34
void getSigmaPlots(double MCSignal, unsigned int timestamp, unsigned int nDays, std::string naming)
TGraphAsymmErrors * GraphWithPoissonErrors(const TH1 *h, bool noErrorsXaxis, bool drawEmptyBins)
Calculate statistical errors appropriate for small Poisson numbers.
Definition: Plots.cxx:912
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:39
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
Double_t ymax
Definition: plot.C:25
void lastUpdated(unsigned int timestamp)
Track finder for cosmic rays.
const Var kPtP
Transverse momentum fraction in slice.
Definition: NueVars.cxx:90
Charged-current interactions.
Definition: IPrediction.h:39
const double timeWindowFactor
Definition: VarsAndCuts.h:55
void makePred(bool usesysts=true, bool useFakeNDData=false)
Definition: makePred.C:64
====================================================================== ///
Definition: CutFlow_Data.C:28
Spectrum Predict(osc::IOscCalc *calc) const override
const Var kBoxMaxY([](const caf::SRProxy *sr){return sr->slc.boxmax.Y()/100;})
const int kNumPlots
Definition: GetSpectra.h:2
const Var kNHit
Definition: Vars.cxx:71
const std::vector< Plot > plots
TGraph * ShadeBetweenHistograms(TH1 *hmin, TH1 *hmax)
Definition: Plots.cxx:1236
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
const Binning kXYBins
Definition: VarsAndCuts.h:102
const double j
Definition: BetheBloch.cxx:29
void SetTh23(const T &th23) override
void SetDmsq21(const T &dmsq21) override
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
void SetDmsq32(const T &dmsq32) override
std::string fname
Definition: CutFlow_Data.C:30
void SetTh12(const T &th12) override
void getSigmaPlotsAll(double MCSignal, unsigned int timestamp, unsigned int nDays, std::string naming)
Neutral-current interactions.
Definition: IPrediction.h:40
void getPredictions()
Definition: getTimePeak.C:88
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
std::string dataName
Prediction that just uses FD MC, with no extrapolation.
All neutrinos, any flavor.
Definition: IPrediction.h:26
Float_t e
Definition: plot.C:35
enum BeamMode kGreen
void getTimePeakPlots(std::string samdef, unsigned int timestamp, unsigned int nDays, std::string naming)
enum BeamMode kBlue
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string