accum_nue_numu_equivalent.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Binning.h"
3 #include "CAFAna/Core/Spectrum.h"
5 #include "CAFAna/Cuts/Cuts.h"
11 #include "CAFAna/Vars/Vars.h"
16 
17 #include "TCanvas.h"
18 #include "TH2.h"
19 #include "TLegend.h"
20 #include "TLatex.h"
21 #include "TFile.h"
22 #include "TProfile.h"
23 #include "TString.h"
24 #include "TPaveText.h"
25 #include "TGaxis.h"
26 #include "TROOT.h"
27 #include "TStyle.h"
28 
29 #include "Utilities/rootlogon.C"
30 
31 using namespace ana;
32 
33 int countEvents(int iStart, int iRun, const long eventList[], int totalEvents);
34 int countFakeEvents(int iRun, long eventList[], int totalEvents);
35 double ksTest(TH1 *expCDF, TH1 *obsCDF, TH1 *expPDF, unsigned int nEvents, unsigned int nTrials);
36 
38  TLatex* CornLab = new TLatex(.1, .93, Str.c_str());
39  CornLab->SetTextColor(kGray+1);
40  CornLab->SetNDC();
41  CornLab->SetTextSize (2/30.);
42  CornLab->SetTextAlign(11);
43  CornLab->Draw();
44 }
45 
46 unsigned int nEvts;
47 unsigned int nEvtsNumu;
48 
49 std::vector<std::pair<int, int>> periodTime = {{1391000000, 1590000000}, // full
50  {1391000000, 1409875200}, // p1 - fhc
51  {1409875200, 1433721600}, // p2 - fhc
52  {1444608000, 1467229500}, // p3 - fhc
53  {1467158400, 1469750400}, // p4 - rhc
54  {1479081600, 1487548800}, // p5 - fhc
55  {1487548800, 1499385600}, // p6 - rhc
56  {1510012800, 1523318400}, // p7 - rhc
57  {1539993600, 1551139200}, // p8 - rhc
58  {1551139200, 1562371200}, // p9 - fhc
59  {1572307200, 1584662400}, // p10 - fhc
60 };
61 
62 // period 0 is the full sample.
64 {
65  nEvts = beam == "fhc" ? 82 : 33;
66  nEvtsNumu = beam == "fhc" ? 211 : 105;
67 
68  // Construct the vector of event times
69  long nueTime[nEvts];
70  long numuTime[nEvtsNumu];
71 
72  std::ifstream evtFileNue(("nue2020_"+beam+"_numi_vars.txt").c_str());
74  double junk;
75  double atime;
76  int i = 0;
77  while (getline(evtFileNue, line)) {
78  istringstream thisLine(line);
79  thisLine >> junk >> junk >> junk >> junk >>
80  atime >>
81  junk >> junk >> junk >> junk >> junk;
82  nueTime[i] = atime;
83  ++i;
84  }
85 
86  std::ifstream evtFileNumu(("numu2020_"+beam+"_numi_vars.txt").c_str());
87  i=0;
88  while (getline(evtFileNumu, line)) {
89  istringstream thisLine(line);
90  thisLine >> junk >> junk >> junk >> junk >>
91  atime >>
92  junk >> junk >> junk >> junk >> junk >> junk >> junk >> junk;
93  numuTime[i] = atime;
94  ++i;
95  }
96 
97  // Following Evan's guidelines more or less for time axis
98  gSystem->Setenv("TZ","UTC");
99  gStyle->SetTimeOffset(0);
100  std::string taxis_labels = "%b '%y";
101  // File with the exposure, livetime, POT plots
102  TFile *fAccounting = new TFile(("nue_ana2020_"+beam+"_pot.root").c_str(),"read");
103  TH1D* hPOTFile = (TH1D*)fAccounting->Get("hPOTMassT"); // switch to using kT*POT
104 
105  // Extract the pot for period
106  TH1D* hPOT;
107  int initTime = periodTime[period].first;
108  int finalTime = periodTime[period].second;
109  if(period == 0 && beam == "rhc"){
110  initTime = 1467000000;
111  finalTime = 1552000000;
112  }
113  int initBin = hPOTFile->GetXaxis()->FindBin(initTime);
114  int finalBin = hPOTFile->GetXaxis()->FindBin(finalTime);
115 
116  initTime = hPOTFile->GetXaxis()->GetBinLowEdge(initBin);
117  finalTime = hPOTFile->GetXaxis()->GetBinLowEdge(finalBin) + hPOTFile->GetXaxis()->GetBinWidth(finalBin);
118 
119  hPOT = new TH1D("hPOT", "", finalBin - initBin + 1, initTime, finalTime);
120 
121  for(int iBin = 1; iBin <= hPOT->GetNbinsX(); ++iBin)
122  hPOT->SetBinContent(iBin, hPOTFile->GetBinContent(iBin + initBin - 1));
123 
124  hPOT->Scale(1./10.43); // scale by 1./10.43kT to get equivalent POT (see docdb-23417)
125  int timeNBins = hPOT->GetXaxis()->GetNbins();
126  double loTime = hPOT->GetXaxis()->GetBinLowEdge(1);
127  double dxTime = hPOT->GetXaxis()->GetBinLowEdge(2) - hPOT->GetXaxis()->GetBinLowEdge(1);
128  double hiTime = hPOT->GetXaxis()->GetBinLowEdge(timeNBins) + dxTime;
129 
130  // Set up time axis
131  hPOT->SetTitle("");
132  hPOT->GetXaxis()->SetTitle("Time");
133  CenterTitles(hPOT);
134  hPOT->GetXaxis()->SetTimeDisplay(1);
135  hPOT->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
136  hPOT->GetXaxis()->SetTitle("Date");
137  hPOT->GetXaxis()->SetNdivisions(506,kFALSE);
138 
139  // Now make a plot of accumulated POT and accumulated event as function of time
140  // ---------------------------------------------------------------------------------
141  gStyle->SetPadTickY(0);
142  std::cout << "Making new histogram with min,max,nbins = " << loTime << "," << hiTime << "," << timeNBins << std::endl;
143  TH1D *hPOTAccumulated = new TH1D("hPOTAccumulated",";Time;Accumulated equivalent POT (x10^{20})",timeNBins,loTime,hiTime);
144  TH1D *hEvtAccumulated = new TH1D("hEvtAccumulated",";Time;Accumulated Events",timeNBins,loTime,hiTime);
145  TH1D *hEvtAccumulatedNumu = new TH1D("hEvtAccumulatedNumu",";Time;Accumulated Events",timeNBins,loTime,hiTime);
146 
147  for( int iRBin = 1; iRBin <= timeNBins; ++iRBin ){
148  if( iRBin==1 )
149  hPOTAccumulated->SetBinContent( iRBin, hPOT->GetBinContent(iRBin)/1.0e20 );
150  else
151  hPOTAccumulated->SetBinContent( iRBin, hPOTAccumulated->GetBinContent(iRBin-1)
152  +(hPOT->GetBinContent(iRBin)/1.0e20) );
153  hEvtAccumulated->SetBinContent( iRBin, countEvents(hEvtAccumulated->GetBinLowEdge(1),
154  hEvtAccumulated->GetBinLowEdge(iRBin)+dxTime,
155  nueTime,nEvts) );
156  hEvtAccumulatedNumu->SetBinContent( iRBin, countEvents(hEvtAccumulated->GetBinLowEdge(1),
157  hEvtAccumulatedNumu->GetBinLowEdge(iRBin)+dxTime,
158  numuTime,nEvtsNumu) );
159 
160  }
161 
162  nEvts = hEvtAccumulated->GetMaximum();
163  nEvtsNumu = hEvtAccumulatedNumu->GetMaximum();
164  double totPOT = hPOTAccumulated->GetMaximum();
165 
166  TCanvas *tcAccumulated = new TCanvas();
167 
168  // Blank plot to draw on
169  TH1D *hPOTAccumBlank = new TH1D("hPOTAccumBlank",";Time;Accumulated equivalent POT (x10^{20})",timeNBins,loTime,hiTime);
170  hPOTAccumBlank->GetYaxis()->SetRangeUser(0.,1.1*hPOTAccumulated->GetMaximum());
171  CenterTitles(hPOTAccumBlank);
172 
173  hPOTAccumBlank->GetXaxis()->SetTimeDisplay(1);
174  hPOTAccumBlank->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
175  hPOTAccumBlank->GetXaxis()->SetTitle("Date");
176  hPOTAccumBlank->GetXaxis()->SetNdivisions(506,kFALSE);
177 
178  hPOTAccumulated->GetXaxis()->SetTimeDisplay(1);
179  hPOTAccumulated->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
180  hPOTAccumulated->GetXaxis()->SetTitle("Date");
181  hPOTAccumulated->GetXaxis()->SetNdivisions(506,kFALSE);
182 
183  hPOTAccumBlank->Draw("hist");
184  hPOTAccumulated->Draw("hist same");
185 
186  // Set up the right axis and overlay the events
187  tcAccumulated->Update();
188  Float_t rightmax = 1.1*hEvtAccumulated->GetMaximum();
189  Float_t rightmaxNumu = 1.1*hEvtAccumulatedNumu->GetMaximum();
190  Float_t scale = gPad->GetUymax()/rightmax;
191  Float_t scaleNumu = gPad->GetUymax()/rightmaxNumu;
192  hEvtAccumulated->SetLineColor(kGreen+3);
193  hEvtAccumulated->SetMarkerColor(kGreen+3);
194  hEvtAccumulated->Scale(scale);
195 
196  hEvtAccumulated->GetXaxis()->SetTimeDisplay(1);
197  hEvtAccumulated->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
198  hEvtAccumulated->GetXaxis()->SetTitle("Date");
199  hEvtAccumulated->GetXaxis()->SetNdivisions(506,kFALSE);
200 
201  hEvtAccumulated->Draw("hist same");
202  hEvtAccumulatedNumu->SetLineColor(kRed);
203  hEvtAccumulatedNumu->SetMarkerColor(kRed);
204  hEvtAccumulatedNumu->Scale(scaleNumu);
205  hEvtAccumulatedNumu->Draw("hist same");
206  gStyle->SetLineWidth(3);
207  TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(), gPad->GetUxmax(),gPad->GetUymax(),0,1.1,506,"+L");
208  axis->SetLineColor(kBlack);
209  axis->SetLabelColor(kBlack);
210  // NOvA style axis options (copied from rootlogon file and set up to work)
211  // --Axis titles
212  axis->SetTitleSize(.055);
213  axis->SetTitleOffset(.8);
214  // --More space for y-axis to avoid clashing with big numbers
215  axis->SetTitleOffset(.9);
216  // --Axis labels (numbering)
217  axis->SetLabelSize(.04);
218  axis->SetLabelOffset(.005);
219  // --Font
220  axis->SetLabelFont(42);
221  axis->SetTitleFont(42);
222  // -----------------------------------------------------------------------
223  axis->SetTitleColor(kBlack);
224  axis->SetTitle("Fraction of Selected Events");
225  axis->CenterTitle();
226  axis->Draw();
227 
228  std::cout<< "Found " << nEvts << " nues," << std::endl;
229  std::cout<< "Found " << nEvtsNumu << " numus," << std::endl;
230  std::cout<< "In " << totPOT << " POT" << std::endl;
231 
232  // Get ks test values
233  double ksTestNue = ksTest(hPOTAccumulated,hEvtAccumulated,hPOT,nEvts,10000);
234  double ksTestNumu = ksTest(hPOTAccumulated,hEvtAccumulatedNumu,hPOT,nEvtsNumu,10000);
235 
236  // Legend
237  tcAccumulated->cd();
238  TLegend *lAccum = new TLegend(0.12,0.5,0.5,0.87);
239  lAccum->SetFillStyle(0);
240  lAccum->AddEntry("hPOTAccumBlank","Accumulated POT","l");
241  //lAccum->AddEntry(hEvtAccumulated->GetName(),TString::Format(beam == "fhc" ? "#nu_{e} candidates (KS prob = %.0f%%)" : "#bar{#nu}_{e} candidates (KS prob = %.0f%%)",ksTestNue*100.),"l");
242 //TString::Format(beam == "fhc" ? "#nu_{#mu} candidates (KS prob = %.0f%%)" : "#bar{#nu}_{#mu} candidates (KS prob = %.0f%%)",ksTestNumu*100.),"l");
243  lAccum->AddEntry(hEvtAccumulated->GetName(),
244  beam == "fhc" ? "#nu_{e} candidates" : "#bar{#nu}_{e} candidates","l");
245  lAccum->AddEntry((TObject*)0, TString::Format("KS prob = %.0f%%",ksTestNue*100), "");
246  lAccum->AddEntry(hEvtAccumulatedNumu->GetName(),
247  beam == "fhc" ? "#nu_{#mu} candidates" : "#bar{#nu}_{#mu} candidates", "l");
248  lAccum->AddEntry((TObject*)0, TString::Format("KS prob = %.0f%%",ksTestNumu*100), "");
249 
250  lAccum->Draw();
251 
252  tcAccumulated->cd();
253  // TPaveText *beam = new TPaveText(0.55,0.13,0.75,0.2,"NB NDC");
254  // beam->AddText("Neutrino Beam");
255  // beam->Draw();
256 
257  Preliminary();
258  MyCornerLabel(beam == "fhc" ? "#nu-beam" : "#bar{#nu}-beam");
259  tcAccumulated->Print(("AccumulatedEvents_Ana2020_"+beam+"_p"+std::to_string(period)+".pdf").c_str());
260 }
261 
262 int countEvents( int iStart, int iRun, const long eventList[], int totEvents ){
263  int accumulatedEvts = 0;
264  for( int iEvtCt=0; iEvtCt<totEvents; ++iEvtCt ){
265  if( eventList[iEvtCt] > iStart && eventList[iEvtCt] < iRun ) accumulatedEvts++; // for time version we pass bin high edge, so want explicitly lower than.
266  }
267  return accumulatedEvts;
268 }
269 
270 int countFakeEvents(int iRun, long eventList[], int totEvents ){
271  int accumulatedEvts = 0;
272  for( int iEvtCt=0; iEvtCt<totEvents; ++iEvtCt ){
273  if( eventList[iEvtCt] < iRun ) accumulatedEvts++; // for time version we pass bin high edge, so want explicitly lower than.
274  }
275  return accumulatedEvts;
276 }
277 
278 // Run pseudo-experiments with POT histograms (normalized here if needed) to test how likely the event distribution is
279 double ksTest( TH1 *expCDF, TH1 *obsCDF, TH1 *expPDF, unsigned int nEvents, unsigned int nTrials ){
280  if( nTrials==0 ) return -999.;
281 
282  // Set up histograms
283  int nCDFBins = expCDF->GetXaxis()->GetNbins();
284  double loBin = expCDF->GetXaxis()->GetBinLowEdge(1);
285  double dxBin = expCDF->GetXaxis()->GetBinLowEdge(2) - expCDF->GetXaxis()->GetBinLowEdge(1);
286  double hiBin = expCDF->GetXaxis()->GetBinLowEdge(nCDFBins);
287 
288  // Clone histograms to be acted upon
289  TH1D *expCDFcpy = (TH1D*)expCDF->Clone();
290  TH1D *expPDFcpy = (TH1D*)expPDF->Clone();
291  TH1D *obsCDFcpy = (TH1D*)obsCDF->Clone();
292  expCDFcpy->SetDirectory(0);
293  expPDFcpy->SetDirectory(0);
294  obsCDFcpy->SetDirectory(0);
295 
296  double expPDFcpyArea = expPDFcpy->Integral("width");
297  expPDFcpy->Scale(1./expPDFcpyArea); // area normalized
298  std::cout << "Check norm: " << fabs(expCDFcpy->GetMaximum()-1.0) << std::endl;
299  if( fabs(expCDFcpy->GetMaximum()-1.0)>1.0e-9 ){
300  // only normalize if not already normalized to have max at 1 (CDF)
301  expCDFcpy->Scale(1.0e20/expPDF->Integral()); // factor 1e20 is because it's factored out of POT histogram
302  }
303 
304  std::cout << "Check norm obsCDF: " << fabs(obsCDFcpy->GetMaximum()-1.0) << std::endl;
305  obsCDFcpy->Scale(1./obsCDFcpy->GetMaximum());
306  std::cout << " and again: " <<fabs(obsCDFcpy->GetMaximum()-1.0) << std::endl;
307 
308  std::cout << "Check PDF: " << expPDFcpy->Integral("width") << std::endl;
309 
310  // -----------------------------
311  // Make some debug plots
312  TCanvas *cExpPDF = new TCanvas();
313  expPDFcpy->Draw("hist");
314  if(nEvents==nEvts) cExpPDF->Print("debug_nue_expPDF_equivalent_time.pdf");
315  else if(nEvents==nEvtsNumu) cExpPDF->Print("debug_numu_expPDF_equivalent_time.pdf");
316 
317  TCanvas *cExpCDF = new TCanvas();
318  expCDFcpy->Draw("hist");
319  if(nEvents==nEvts) cExpCDF->Print("debug_nue_expCDF_equivalent_time.pdf");
320  else if(nEvents==nEvtsNumu) cExpCDF->Print("debug_numu_expCDF_equivalent_time.pdf");
321 
322  TCanvas *cObsCDF = new TCanvas();
323  obsCDFcpy->Draw("hist");
324  if(nEvents==nEvts) cObsCDF->Print("debug_nue_obsCDF_equivalent_time.pdf");
325  else if(nEvents==nEvtsNumu) cObsCDF->Print("debug_numu_obsCDF_equivalent_time.pdf");
326  // -----------------------------
327 
328  // Maximum difference between event CDF and POT CDF (value for KS test)
329  double obsMaxDiff = 0.;
330  double locMaxDiff = 0.;
331  int binMaxDiff = 0;
332  for( int iBin=1; iBin<=nCDFBins; ++iBin ){
333  if( fabs(expCDFcpy->GetBinContent(iBin)-obsCDFcpy->GetBinContent(iBin)) > obsMaxDiff){
334  obsMaxDiff = fabs(expCDFcpy->GetBinContent(iBin)-obsCDFcpy->GetBinContent(iBin));
335  locMaxDiff = expCDFcpy->GetBinCenter(iBin);
336  binMaxDiff = iBin;
337  }
338  }
339 
340  // -----------------------------
341  // Make some overlay debug plots
342  TCanvas *cDBOverlay = new TCanvas();
343  expCDFcpy->Draw("hist");
344  obsCDFcpy->Draw("hist same");
345  double minOfDiff = std::min(obsCDFcpy->GetBinContent(binMaxDiff),expCDFcpy->GetBinContent(binMaxDiff));
346  double maxOfDiff = std::max(obsCDFcpy->GetBinContent(binMaxDiff),expCDFcpy->GetBinContent(binMaxDiff));
347  TLine *overlayLine = new TLine(locMaxDiff,minOfDiff,locMaxDiff,maxOfDiff);
348  overlayLine->SetLineColor(kRed);
349  overlayLine->SetLineWidth(2);
350  overlayLine->SetLineStyle(7);
351  overlayLine->Draw();
352  TPaveText *txtOverlay = new TPaveText(0.15,0.6,0.4,0.87,"NB NDC");
353  txtOverlay->AddText( TString::Format("Max difference: %.3f",maxOfDiff-minOfDiff) );
354  txtOverlay->Draw();
355  if(nEvents==nEvts) cDBOverlay->Print("debug_nue_overlay_equivalent_time.pdf");
356  else if(nEvents==nEvtsNumu) cDBOverlay->Print("debug_numu_overlay_equivalent_time.pdf");
357  // -----------------------------
358 
359  std::cout << "Max Difference: " << obsMaxDiff << std::endl;
360 
361  unsigned int iPassTrial=0;
362 
363  TH1D *ksDistances;
364  if( nEvents==nEvts ) ksDistances = new TH1D("ksDistancesNue",";max separation;fake trials",50,0,1);
365  else if( nEvents==nEvtsNumu ) ksDistances = new TH1D("ksDistancesNumu",";max separation;fake trials",50,0,1);
366  else return -999.;
367  long fakeEvents[nEvents];
368  TH1D *trialCDF;
369  if( nEvents==nEvts ) trialCDF = new TH1D("trialCDFNue","",nCDFBins,loBin,hiBin+dxBin);
370  else if( nEvents==nEvtsNumu ) trialCDF = new TH1D("trialCDFNumu","",nCDFBins,loBin,hiBin+dxBin);
371  else return -999.;
372 
373  // Pseudo-experiments
374  for( unsigned int iTrial=0; iTrial<nTrials; ++iTrial ){
375  // Pick event times
376  for( unsigned int iEvt=0; iEvt<nEvents; ++iEvt ){
377  fakeEvents[iEvt] = expPDFcpy->GetRandom();
378  }
379  // Fill fake CDF
380  for( int iRBin = 1; iRBin <= nCDFBins; ++iRBin ){
381  trialCDF->SetBinContent( iRBin, countFakeEvents(trialCDF->GetBinCenter(iRBin),fakeEvents,nEvents) );
382  }
383  // Normalize CDF
384  if( nEvents>0 ) trialCDF->Scale(1./double(nEvents));
385  // Get the ks distance
386  double trialMaxDiff = 0.;
387  for( int iBin=1; iBin<=nCDFBins; ++iBin ){
388  if( fabs(expCDFcpy->GetBinContent(iBin)-trialCDF->GetBinContent(iBin)) > trialMaxDiff)
389  trialMaxDiff = fabs(expCDFcpy->GetBinContent(iBin)-trialCDF->GetBinContent(iBin));
390  }
391  ksDistances->Fill(trialMaxDiff);
392  // If greater then observed max difference, count
393  if(trialMaxDiff>=obsMaxDiff) iPassTrial+=1;
394  // Reset histogram
395  for( int iBin=1; iBin<=nCDFBins; ++iBin ){
396  trialCDF->SetBinContent(iBin,0);
397  }
398  }
399 
400  // -----------------------------
401  // Area normalize KS histogram and draw
402  double ksArea = ksDistances->Integral();
403  ksDistances->Scale(1./ksArea);
404 
405  TCanvas *cKSTest = new TCanvas();
406  ksDistances->Draw();
407  TLine *obsMaxDiffLine = new TLine(obsMaxDiff,0.,obsMaxDiff,0.95*ksDistances->GetMaximum());
408  obsMaxDiffLine->SetLineColor(kRed);
409  obsMaxDiffLine->Draw();
410  if( nEvents==nEvts )
411  cKSTest->Print("ksTestNue_equivalent_time.pdf");
412  if( nEvents==nEvtsNumu )
413  cKSTest->Print("ksTestNumu_equivalent_time.pdf");
414  // -----------------------------
415 
416  // Result of the test is the fraction of pseudo-
417  // experiments with KS distance >= observed KS distance
418  return double(iPassTrial)/double(nTrials);
419 }
void MyCornerLabel(std::string Str)
T max(const caf::Proxy< T > &a, T b)
void accum_nue_numu_equivalent(std::string beam="fhc", int period=0)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
unsigned int nEvts
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
double ksTest(TH1 *expCDF, TH1 *obsCDF, TH1 *expPDF, unsigned int nEvents, unsigned int nTrials)
unsigned int nEvtsNumu
Double_t scale
Definition: plot.C:25
int countEvents(int iStart, int iRun, const long eventList[], int totalEvents)
void Preliminary()
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::vector< std::pair< int, int > > periodTime
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
int countFakeEvents(int iRun, long eventList[], int totalEvents)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
T min(const caf::Proxy< T > &a, T b)