accum_nue_numu_equivalent_fhc.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 "/nova/app/users/howard/Dev_Nue2018_Commit/CAFAna_includes/fhc_results.h"
30 
31 #include "Utilities/rootlogon.C"
32 
33 using namespace ana;
34 
35 int countEvents(int iRun, const long eventList[], int totalEvents);
36 int countFakeEvents(int iRun, long eventList[], int totalEvents);
37 double ksTest(TH1 *expCDF, TH1 *obsCDF, TH1 *expPDF, unsigned int nEvents, unsigned int nTrials);
38 
40 {
41  // Following Evan's guidelines more or less for time axis
42  gSystem->Setenv("TZ","UTC");
43  gStyle->SetTimeOffset(0);
44  std::string taxis_labels = "%b '%y";
45 
46  // File with the exposure, livetime, POT plots
47  TFile *fAccounting = new TFile("/nova/ana/users/howard/Ana2018-pot/nue_ana2018_pot-fhc.root","read");
48  TH1D* hPOT = (TH1D*)fAccounting->Get("hPOTMassT"); // switch to using kT*POT
49  hPOT->Scale(1./10.43); // scale by 1./10.43kT to get equivalent POT (see docdb-23417)
50  int timeNBins = hPOT->GetXaxis()->GetNbins();
51  double loTime = hPOT->GetXaxis()->GetBinLowEdge(1);
52  double dxTime = hPOT->GetXaxis()->GetBinLowEdge(2) - hPOT->GetXaxis()->GetBinLowEdge(1);
53  double hiTime = hPOT->GetXaxis()->GetBinLowEdge(timeNBins) + dxTime;
54 
55  // Set up time axis
56  hPOT->SetTitle("");
57  hPOT->GetXaxis()->SetTitle("Time");
58  CenterTitles(hPOT);
59  hPOT->GetXaxis()->SetTimeDisplay(1);
60  hPOT->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
61  hPOT->GetXaxis()->SetTitle("Date");
62  hPOT->GetXaxis()->SetNdivisions(506,kFALSE);
63 
64  // Now make a plot of accumulated POT and accumulated event as function of time
65  // ---------------------------------------------------------------------------------
66  gStyle->SetPadTickY(0);
67  std::cout << "Making new histogram with min,max,nbins = " << loTime << "," << hiTime << "," << timeNBins << std::endl;
68  TH1D *hPOTAccumulated = new TH1D("hPOTAccumulated",";Time;Accumulated equivalent POT (x10^{20})",timeNBins,loTime,hiTime);
69  TH1D *hEvtAccumulated = new TH1D("hEvtAccumulated",";Time;Accumulated Events",timeNBins,loTime,hiTime);
70  TH1D *hEvtAccumulatedNumu = new TH1D("hEvtAccumulatedNumu",";Time;Accumulated Events",timeNBins,loTime,hiTime);
71 
72  for( int iRBin = 1; iRBin <= timeNBins; ++iRBin ){
73  if( iRBin==1 )
74  hPOTAccumulated->SetBinContent( iRBin, hPOT->GetBinContent(iRBin)/1.0e20 );
75  else
76  hPOTAccumulated->SetBinContent( iRBin, hPOTAccumulated->GetBinContent(iRBin-1)+(hPOT->GetBinContent(iRBin)/1.0e20) );
77  hEvtAccumulated->SetBinContent( iRBin, countEvents(hEvtAccumulated->GetBinLowEdge(iRBin)+dxTime,nueTime,nEvts) );
78  hEvtAccumulatedNumu->SetBinContent( iRBin, countEvents(hEvtAccumulatedNumu->GetBinLowEdge(iRBin)+dxTime,numuTime,nEvtsNumu) );
79  }
80 
81  TCanvas *tcAccumulated = new TCanvas();
82 
83  // Blank plot to draw on
84  TH1D *hPOTAccumBlank = new TH1D("hPOTAccumBlank",";Time;Accumulated equivalent POT (x10^{20})",timeNBins,loTime,hiTime);
85  hPOTAccumBlank->GetXaxis()->SetRangeUser(loTime,1490000000);
86  hPOTAccumBlank->GetYaxis()->SetRangeUser(0.,1.1*hPOTAccumulated->GetMaximum());
87  CenterTitles(hPOTAccumBlank);
88 
89  hPOTAccumBlank->GetXaxis()->SetTimeDisplay(1);
90  hPOTAccumBlank->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
91  hPOTAccumBlank->GetXaxis()->SetTitle("Date");
92  hPOTAccumBlank->GetXaxis()->SetNdivisions(506,kFALSE);
93 
94  hPOTAccumulated->GetXaxis()->SetTimeDisplay(1);
95  hPOTAccumulated->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
96  hPOTAccumulated->GetXaxis()->SetTitle("Date");
97  hPOTAccumulated->GetXaxis()->SetNdivisions(506,kFALSE);
98 
99  hPOTAccumBlank->Draw("hist");
100  hPOTAccumulated->Draw("hist same");
101 
102  // Set up the right axis and overlay the events
103  tcAccumulated->Update();
104  Float_t rightmax = 1.1*hEvtAccumulated->GetMaximum();
105  Float_t rightmaxNumu = 1.1*hEvtAccumulatedNumu->GetMaximum();
106  Float_t scale = gPad->GetUymax()/rightmax;
107  Float_t scaleNumu = gPad->GetUymax()/rightmaxNumu;
108  hEvtAccumulated->SetLineColor(kGreen+3);
109  hEvtAccumulated->SetMarkerColor(kGreen+3);
110  hEvtAccumulated->Scale(scale);
111 
112  hEvtAccumulated->GetXaxis()->SetTimeDisplay(1);
113  hEvtAccumulated->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
114  hEvtAccumulated->GetXaxis()->SetTitle("Date");
115  hEvtAccumulated->GetXaxis()->SetNdivisions(506,kFALSE);
116 
117  hEvtAccumulated->Draw("hist same");
118  hEvtAccumulatedNumu->SetLineColor(kRed);
119  hEvtAccumulatedNumu->SetMarkerColor(kRed);
120  hEvtAccumulatedNumu->Scale(scaleNumu);
121  hEvtAccumulatedNumu->Draw("hist same");
122  gStyle->SetLineWidth(3);
123  TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(), gPad->GetUxmax(),gPad->GetUymax(),0,1.1,506,"+L");
124  axis->SetLineColor(kBlack);
125  axis->SetLabelColor(kBlack);
126  // NOvA style axis options (copied from rootlogon file and set up to work)
127  // --Axis titles
128  axis->SetTitleSize(.055);
129  axis->SetTitleOffset(.8);
130  // --More space for y-axis to avoid clashing with big numbers
131  axis->SetTitleOffset(.9);
132  // --Axis labels (numbering)
133  axis->SetLabelSize(.04);
134  axis->SetLabelOffset(.005);
135  // --Font
136  axis->SetLabelFont(42);
137  axis->SetTitleFont(42);
138  // -----------------------------------------------------------------------
139  axis->SetTitleColor(kBlack);
140  axis->SetTitle("Fraction of Selected Events");
141  axis->CenterTitle();
142  axis->Draw();
143 
144  // Get ks test values
145  double ksTestNue = ksTest(hPOTAccumulated,hEvtAccumulated,hPOT,nEvts,10000);
146  double ksTestNumu = ksTest(hPOTAccumulated,hEvtAccumulatedNumu,hPOT,nEvtsNumu,10000);
147 
148  // Legend
149  tcAccumulated->cd();
150  TLegend *lAccum = new TLegend(0.12,0.55,0.55,0.87);
151  lAccum->SetFillStyle(0);
152  lAccum->AddEntry("hPOTAccumBlank","Accumulated POT","l");
153  lAccum->AddEntry(hEvtAccumulated->GetName(),TString::Format("#nu_{e} candidates (KS prob = %.0f%%)",ksTestNue*100.),"l");
154  lAccum->AddEntry(hEvtAccumulatedNumu->GetName(),TString::Format("#nu_{#mu} candidates (KS prob = %.0f%%)",ksTestNumu*100.),"l");
155  lAccum->Draw();
156 
157  tcAccumulated->cd();
158  TPaveText *beam = new TPaveText(0.55,0.13,0.75,0.2,"NB NDC");
159  beam->AddText("Neutrino Beam");
160  beam->Draw();
161 
162  Preliminary();
163  tcAccumulated->Print("AccumulatedEvents_Ana2018_fhc_equivalent_time.pdf");
164 }
165 
166 int countEvents( int iRun, const long eventList[], int totEvents ){
167  int accumulatedEvts = 0;
168  for( int iEvtCt=0; iEvtCt<totEvents; ++iEvtCt ){
169  if( eventList[iEvtCt] < iRun ) accumulatedEvts++; // for time version we pass bin high edge, so want explicitly lower than.
170  }
171  return accumulatedEvts;
172 }
173 
174 int countFakeEvents( int iRun, long eventList[], int totEvents ){
175  int accumulatedEvts = 0;
176  for( int iEvtCt=0; iEvtCt<totEvents; ++iEvtCt ){
177  if( eventList[iEvtCt] < iRun ) accumulatedEvts++; // for time version we pass bin high edge, so want explicitly lower than.
178  }
179  return accumulatedEvts;
180 }
181 
182 // Run pseudo-experiments with POT histograms (normalized here if needed) to test how likely the event distribution is
183 double ksTest( TH1 *expCDF, TH1 *obsCDF, TH1 *expPDF, unsigned int nEvents, unsigned int nTrials ){
184  if( nTrials==0 ) return -999.;
185 
186  // Set up histograms
187  int nCDFBins = expCDF->GetXaxis()->GetNbins();
188  double loBin = expCDF->GetXaxis()->GetBinLowEdge(1);
189  double dxBin = expCDF->GetXaxis()->GetBinLowEdge(2) - expCDF->GetXaxis()->GetBinLowEdge(1);
190  double hiBin = expCDF->GetXaxis()->GetBinLowEdge(nCDFBins);
191 
192  // Clone histograms to be acted upon
193  TH1D *expCDFcpy = (TH1D*)expCDF->Clone();
194  TH1D *expPDFcpy = (TH1D*)expPDF->Clone();
195  TH1D *obsCDFcpy = (TH1D*)obsCDF->Clone();
196  expCDFcpy->SetDirectory(0);
197  expPDFcpy->SetDirectory(0);
198  obsCDFcpy->SetDirectory(0);
199 
200  double expPDFcpyArea = expPDFcpy->Integral("width");
201  expPDFcpy->Scale(1./expPDFcpyArea); // area normalized
202  std::cout << "Check norm: " << fabs(expCDFcpy->GetMaximum()-1.0) << std::endl;
203  if( fabs(expCDFcpy->GetMaximum()-1.0)>1.0e-9 ){
204  // only normalize if not already normalized to have max at 1 (CDF)
205  expCDFcpy->Scale(1.0e20/expPDF->Integral()); // factor 1e20 is because it's factored out of POT histogram
206  }
207  obsCDFcpy->Scale(1./obsCDFcpy->GetMaximum());
208 
209  std::cout << "Check PDF: " << expPDFcpy->Integral("width") << std::endl;
210 
211  // -----------------------------
212  // Make some debug plots
213  TCanvas *cExpPDF = new TCanvas();
214  expPDFcpy->Draw("hist");
215  if(nEvents==nEvts) cExpPDF->Print("debug_nue_expPDF_equivalent_time_fhc.pdf");
216  else if(nEvents==nEvtsNumu) cExpPDF->Print("debug_numu_expPDF_equivalent_time_fhc.pdf");
217 
218  TCanvas *cExpCDF = new TCanvas();
219  expCDFcpy->Draw("hist");
220  if(nEvents==nEvts) cExpCDF->Print("debug_nue_expCDF_equivalent_time_fhc.pdf");
221  else if(nEvents==nEvtsNumu) cExpCDF->Print("debug_numu_expCDF_equivalent_time_fhc.pdf");
222 
223  TCanvas *cObsCDF = new TCanvas();
224  obsCDFcpy->Draw("hist");
225  if(nEvents==nEvts) cObsCDF->Print("debug_nue_obsCDF_equivalent_time_fhc.pdf");
226  else if(nEvents==nEvtsNumu) cObsCDF->Print("debug_numu_obsCDF_equivalent_time_fhc.pdf");
227  // -----------------------------
228 
229  // Maximum difference between event CDF and POT CDF (value for KS test)
230  double obsMaxDiff = 0.;
231  double locMaxDiff = 0.;
232  int binMaxDiff = 0;
233  for( int iBin=1; iBin<=nCDFBins; ++iBin ){
234  if( fabs(expCDFcpy->GetBinContent(iBin)-obsCDFcpy->GetBinContent(iBin)) > obsMaxDiff){
235  obsMaxDiff = fabs(expCDFcpy->GetBinContent(iBin)-obsCDFcpy->GetBinContent(iBin));
236  locMaxDiff = expCDFcpy->GetBinCenter(iBin);
237  binMaxDiff = iBin;
238  }
239  }
240 
241  // -----------------------------
242  // Make some overlay debug plots
243  TCanvas *cDBOverlay = new TCanvas();
244  expCDFcpy->Draw("hist");
245  obsCDFcpy->Draw("hist same");
246  double minOfDiff = std::min(obsCDFcpy->GetBinContent(binMaxDiff),expCDFcpy->GetBinContent(binMaxDiff));
247  double maxOfDiff = std::max(obsCDFcpy->GetBinContent(binMaxDiff),expCDFcpy->GetBinContent(binMaxDiff));
248  TLine *overlayLine = new TLine(locMaxDiff,minOfDiff,locMaxDiff,maxOfDiff);
249  overlayLine->SetLineColor(kRed);
250  overlayLine->SetLineWidth(2);
251  overlayLine->SetLineStyle(7);
252  overlayLine->Draw();
253  TPaveText *txtOverlay = new TPaveText(0.15,0.6,0.4,0.87,"NB NDC");
254  txtOverlay->AddText( TString::Format("Max difference: %.3f",maxOfDiff-minOfDiff) );
255  txtOverlay->Draw();
256  if(nEvents==nEvts) cDBOverlay->Print("debug_nue_overlay_equivalent_time_fhc.pdf");
257  else if(nEvents==nEvtsNumu) cDBOverlay->Print("debug_numu_overlay_equivalent_time_fhc.pdf");
258  // -----------------------------
259 
260  std::cout << "Max Difference: " << obsMaxDiff << std::endl;
261 
262  unsigned int iPassTrial=0;
263 
264  TH1D *ksDistances;
265  if( nEvents==nEvts ) ksDistances = new TH1D("ksDistancesNue",";max separation;fake trials",50,0,1);
266  else if( nEvents==nEvtsNumu ) ksDistances = new TH1D("ksDistancesNumu",";max separation;fake trials",50,0,1);
267  else return -999.;
268  long fakeEvents[nEvtsNumu];
269  TH1D *trialCDF;
270  if( nEvents==nEvts ) trialCDF = new TH1D("trialCDFNue","",nCDFBins,loBin,hiBin+dxBin);
271  else if( nEvents==nEvtsNumu ) trialCDF = new TH1D("trialCDFNumu","",nCDFBins,loBin,hiBin+dxBin);
272  else return -999.;
273 
274  // Pseudo-experiments
275  for( unsigned int iTrial=0; iTrial<nTrials; ++iTrial ){
276  // Pick event times
277  for( unsigned int iEvt=0; iEvt<nEvents; ++iEvt ){
278  fakeEvents[iEvt] = expPDFcpy->GetRandom();
279  }
280  // Fill fake CDF
281  for( int iRBin = 1; iRBin <= nCDFBins; ++iRBin ){
282  trialCDF->SetBinContent( iRBin, countFakeEvents(trialCDF->GetBinCenter(iRBin),fakeEvents,nEvents) );
283  }
284  // Normalize CDF
285  if( nEvents>0 ) trialCDF->Scale(1./double(nEvents));
286  // Get the ks distance
287  double trialMaxDiff = 0.;
288  for( int iBin=1; iBin<=nCDFBins; ++iBin ){
289  if( fabs(expCDFcpy->GetBinContent(iBin)-trialCDF->GetBinContent(iBin)) > trialMaxDiff)
290  trialMaxDiff = fabs(expCDFcpy->GetBinContent(iBin)-trialCDF->GetBinContent(iBin));
291  }
292  ksDistances->Fill(trialMaxDiff);
293  // If greater then observed max difference, count
294  if(trialMaxDiff>=obsMaxDiff) iPassTrial+=1;
295  // Reset histogram
296  for( int iBin=1; iBin<=nCDFBins; ++iBin ){
297  trialCDF->SetBinContent(iBin,0);
298  }
299  }
300 
301  // -----------------------------
302  // Area normalize KS histogram and draw
303  double ksArea = ksDistances->Integral();
304  ksDistances->Scale(1./ksArea);
305 
306  TCanvas *cKSTest = new TCanvas();
307  ksDistances->Draw();
308  TLine *obsMaxDiffLine = new TLine(obsMaxDiff,0.,obsMaxDiff,0.95*ksDistances->GetMaximum());
309  obsMaxDiffLine->SetLineColor(kRed);
310  obsMaxDiffLine->Draw();
311  if( nEvents==nEvts )
312  cKSTest->Print("ksTestNue_equivalent_time_fhc.pdf");
313  if( nEvents==nEvtsNumu )
314  cKSTest->Print("ksTestNumu_equivalent_time_fhc.pdf");
315  // -----------------------------
316 
317  // Result of the test is the fraction of pseudo-
318  // experiments with KS distance >= observed KS distance
319  return double(iPassTrial)/double(nTrials);
320 }
T max(const caf::Proxy< T > &a, T b)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double ksTest(TH1 *expCDF, TH1 *obsCDF, TH1 *expPDF, unsigned int nEvents, unsigned int nTrials)
unsigned int nEvts
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
unsigned int nEvtsNumu
Double_t scale
Definition: plot.C:25
void Preliminary()
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
int countEvents(int iRun, const long eventList[], int totalEvents)
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)
T min(const caf::Proxy< T > &a, T b)
void accum_nue_numu_equivalent_fhc()