2 #include "CAFAna/Core/Binning.h" 24 #include "TPaveText.h" 29 #include "/nova/app/users/howard/Dev_Nue2018_Commit/CAFAna_includes/fhc_results.h" 35 int countEvents(
int iRun,
const long eventList[],
int totalEvents);
37 double ksTest(TH1 *expCDF, TH1 *obsCDF, TH1 *expPDF,
unsigned int nEvents,
unsigned int nTrials);
42 gSystem->Setenv(
"TZ",
"UTC");
43 gStyle->SetTimeOffset(0);
47 TFile *fAccounting =
new TFile(
"/nova/ana/users/howard/Ana2018-pot/nue_ana2018_pot-fhc.root",
"read");
48 TH1D* hPOT = (TH1D*)fAccounting->Get(
"hPOTMassT");
49 hPOT->Scale(1./10.43);
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;
57 hPOT->GetXaxis()->SetTitle(
"Time");
59 hPOT->GetXaxis()->SetTimeDisplay(1);
60 hPOT->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
61 hPOT->GetXaxis()->SetTitle(
"Date");
62 hPOT->GetXaxis()->SetNdivisions(506,kFALSE);
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);
72 for(
int iRBin = 1; iRBin <= timeNBins; ++iRBin ){
74 hPOTAccumulated->SetBinContent( iRBin, hPOT->GetBinContent(iRBin)/1.0e20 );
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) );
81 TCanvas *tcAccumulated =
new TCanvas();
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());
89 hPOTAccumBlank->GetXaxis()->SetTimeDisplay(1);
90 hPOTAccumBlank->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
91 hPOTAccumBlank->GetXaxis()->SetTitle(
"Date");
92 hPOTAccumBlank->GetXaxis()->SetNdivisions(506,kFALSE);
94 hPOTAccumulated->GetXaxis()->SetTimeDisplay(1);
95 hPOTAccumulated->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
96 hPOTAccumulated->GetXaxis()->SetTitle(
"Date");
97 hPOTAccumulated->GetXaxis()->SetNdivisions(506,kFALSE);
99 hPOTAccumBlank->Draw(
"hist");
100 hPOTAccumulated->Draw(
"hist same");
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);
112 hEvtAccumulated->GetXaxis()->SetTimeDisplay(1);
113 hEvtAccumulated->GetXaxis()->SetTimeFormat(taxis_labels.c_str());
114 hEvtAccumulated->GetXaxis()->SetTitle(
"Date");
115 hEvtAccumulated->GetXaxis()->SetNdivisions(506,kFALSE);
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);
128 axis->SetTitleSize(.055);
129 axis->SetTitleOffset(.8);
131 axis->SetTitleOffset(.9);
133 axis->SetLabelSize(.04);
134 axis->SetLabelOffset(.005);
136 axis->SetLabelFont(42);
137 axis->SetTitleFont(42);
139 axis->SetTitleColor(kBlack);
140 axis->SetTitle(
"Fraction of Selected Events");
145 double ksTestNue =
ksTest(hPOTAccumulated,hEvtAccumulated,hPOT,
nEvts,10000);
146 double ksTestNumu =
ksTest(hPOTAccumulated,hEvtAccumulatedNumu,hPOT,
nEvtsNumu,10000);
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");
158 TPaveText *
beam =
new TPaveText(0.55,0.13,0.75,0.2,
"NB NDC");
159 beam->AddText(
"Neutrino Beam");
163 tcAccumulated->Print(
"AccumulatedEvents_Ana2018_fhc_equivalent_time.pdf");
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++;
171 return accumulatedEvts;
175 int accumulatedEvts = 0;
176 for(
int iEvtCt=0; iEvtCt<totEvents; ++iEvtCt ){
177 if( eventList[iEvtCt] < iRun ) accumulatedEvts++;
179 return accumulatedEvts;
183 double ksTest( TH1 *expCDF, TH1 *obsCDF, TH1 *expPDF,
unsigned int nEvents,
unsigned int nTrials ){
184 if( nTrials==0 )
return -999.;
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);
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);
200 double expPDFcpyArea = expPDFcpy->Integral(
"width");
201 expPDFcpy->Scale(1./expPDFcpyArea);
203 if( fabs(expCDFcpy->GetMaximum()-1.0)>1.0e-9 ){
205 expCDFcpy->Scale(1.0e20/expPDF->Integral());
207 obsCDFcpy->Scale(1./obsCDFcpy->GetMaximum());
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");
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");
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");
230 double obsMaxDiff = 0.;
231 double locMaxDiff = 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);
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);
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) );
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");
262 unsigned int iPassTrial=0;
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);
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);
275 for(
unsigned int iTrial=0; iTrial<nTrials; ++iTrial ){
277 for(
unsigned int iEvt=0; iEvt<
nEvents; ++iEvt ){
278 fakeEvents[iEvt] = expPDFcpy->GetRandom();
281 for(
int iRBin = 1; iRBin <= nCDFBins; ++iRBin ){
285 if( nEvents>0 ) trialCDF->Scale(1./
double(nEvents));
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));
292 ksDistances->Fill(trialMaxDiff);
294 if(trialMaxDiff>=obsMaxDiff) iPassTrial+=1;
296 for(
int iBin=1; iBin<=nCDFBins; ++iBin ){
297 trialCDF->SetBinContent(iBin,0);
303 double ksArea = ksDistances->Integral();
304 ksDistances->Scale(1./ksArea);
306 TCanvas *cKSTest =
new TCanvas();
308 TLine *obsMaxDiffLine =
new TLine(obsMaxDiff,0.,obsMaxDiff,0.95*ksDistances->GetMaximum());
309 obsMaxDiffLine->SetLineColor(
kRed);
310 obsMaxDiffLine->Draw();
312 cKSTest->Print(
"ksTestNue_equivalent_time_fhc.pdf");
314 cKSTest->Print(
"ksTestNumu_equivalent_time_fhc.pdf");
319 return double(iPassTrial)/double(nTrials);
T max(const caf::Proxy< T > &a, T b)
Cuts and Vars for the 2020 FD DiF Study.
double ksTest(TH1 *expCDF, TH1 *obsCDF, TH1 *expPDF, unsigned int nEvents, unsigned int nTrials)
void CenterTitles(TH1 *histo)
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)
int countFakeEvents(int iRun, long eventList[], int totalEvents)
T min(const caf::Proxy< T > &a, T b)
void accum_nue_numu_equivalent_fhc()