9 #include "CAFAna/Vars/NumuVars.h" 17 #include "TGraphErrors.h" 18 #include "TGraphAsymmErrors.h" 25 #include "TSpectrum.h" 36 gStyle->SetStatFormat(
"6.3g");
37 gStyle->SetPaintTextFormat(
"6.3g");
38 gStyle->SetFitFormat(
"6.3g");
41 TFile *f_Out =
new TFile((TString)(outDir+
"BeamMonitoring_LongTerm.root"),
"RECREATE");
42 TDirectory *dir_Profile = f_Out->mkdir(
"Profiles");
43 TDirectory *dir_1D = f_Out->mkdir(
"1Ds");
44 TDirectory *dir_Ratio = f_Out->mkdir(
"Ratios");
47 unsigned int startTime = endTime - 86400*nDays;
48 const std::string fFiles =
"dataset_def_name_newest_snapshot " + samDef + Form(
" and online.subrunendtime >= '%u' and online.subrunendtime <= '%u'", startTime, endTime);
53 const Cut kNumuNDCut = kNumuCut2020ND;
57 std::vector<TH1D*> vec_hPOT;
61 TH1D *hPOT =
new TH1D(Form(
"hPOT_Period%i",
i),
";Central Time;POT (10^{15})", nDays, startTime, endTime);
62 vec_hPOT.push_back(hPOT);
69 {
"Central Time",
"Reconstructed Neutrino Energy, (GeV)",
"ccE",
Binning::Simple(10,0,5),
kCCE, kNumuNDCut, 1.5, 2.2 },
70 {
"Central Time",
"Reconstructed Hadronic Energy, (GeV)",
"hadE",
Binning::Simple(12,0,3),
kHadE, kNumuNDCut, 0.2, 0.6 },
71 {
"Central Time",
"Reconstructed Muon Energy, (GeV)",
"muonE",
Binning::Simple(15,0,4.5),
kMuE, kNumuNDCut, 0.9, 1.7 },
72 {
"Central Time",
"Hadronic Energy Fraction",
"hadEfrac",
Binning::Simple(10,0,1.0),
kHadEFrac, kNumuNDCut, 0.05, 0.35},
73 {
"Central Time",
"No. of Hits",
"nhit",
Binning::Simple(10,0,250),
kNHit, kNumuNDCut, 85, 105 },
74 {
"Central Time",
"Vertex X, (m)",
"startX",
Binning::Simple(10,-2,2),
kTrkStartX, kNumuNDCut, -0.40, 0.0 },
75 {
"Central Time",
"Vertex Y, (m)",
"startY",
Binning::Simple(10,-2,2),
kTrkStartY, kNumuNDCut, -0.25, 0.15},
76 {
"Central Time",
"Vertex Z, (m)",
"startZ",
Binning::Simple(12,0,12),
kTrkStartZ, kNumuNDCut, 4.5, 5.8 }
79 std::map<std::pair<unsigned int,unsigned int>,
Spectrum*> map_PeriodAndVarToSpectrum;
80 std::map<std::string,TCanvas*> map_VarToCanvas;
81 std::map<std::string,TH1D*> map_VarToMeanHist;
83 for(
unsigned int i = 0;
i < vec_hPOT.size();
i++)
85 for(
unsigned int j = 0;
j < vec_Plots.size();
j++)
87 Spectrum *h_Spectrum =
new Spectrum(vec_Plots.at(
j).fLabelX, vec_Plots.at(
j).fLabelY,
loader, kTimeBinning,
kTime, vec_Plots.at(
j).fBins, vec_Plots.at(
j).fVar,
89 map_PeriodAndVarToSpectrum[{
i,
j}] = h_Spectrum;
92 map_VarToCanvas[vec_Plots.at(
j).fName] =
new TCanvas((TString)(
"c_"+vec_Plots.at(
j).fName+
"_Profile"),(TString)(
"c_"+vec_Plots.at(
j).fName+
"_Profile"));
93 map_VarToCanvas[vec_Plots.at(
j).fName]->SetBottomMargin(0.15);
94 map_VarToCanvas[vec_Plots.at(
j).fName]->SetLeftMargin (0.15);
95 map_VarToCanvas[vec_Plots.at(
j).fName]->SetRightMargin (0.15);
96 map_VarToCanvas[vec_Plots.at(
j).fName+
"_1D"] =
new TCanvas((TString)(
"c_"+vec_Plots.at(
j).fName+
"_1D"),(TString)(
"c_"+vec_Plots.at(
j).fName+
"_1D"),2*map_VarToCanvas[vec_Plots.at(
j).fName]->GetWw(),
97 map_VarToCanvas[vec_Plots.at(
j).fName]->GetWh());
98 map_VarToCanvas[vec_Plots.at(
j).fName+
"_1D"]->SetBottomMargin(0.15);
99 map_VarToCanvas[vec_Plots.at(
j).fName+
"_1D"]->SetLeftMargin (0.15);
100 map_VarToCanvas[vec_Plots.at(
j).fName+
"_1D"]->SetRightMargin (0.15);
101 map_VarToCanvas[vec_Plots.at(
j).fName+
"_Ratio"] =
new TCanvas((TString)(
"c_"+vec_Plots.at(
j).fName+
"_Ratio"),(TString)(
"c_"+vec_Plots.at(
j).fName+
"_Ratio"),2*map_VarToCanvas[vec_Plots.at(
j).fName]->GetWw(),
102 map_VarToCanvas[vec_Plots.at(
j).fName]->GetWh());
103 map_VarToCanvas[vec_Plots.at(
j).fName+
"_Ratio"]->SetBottomMargin(0.15);
104 map_VarToCanvas[vec_Plots.at(
j).fName+
"_Ratio"]->SetLeftMargin (0.15);
105 map_VarToCanvas[vec_Plots.at(
j).fName+
"_Ratio"]->SetRightMargin (0.15);
107 TH1D *h_Mean =
new TH1D();
108 h_Mean->SetName((TString)vec_Plots.at(
j).fName);
109 h_Mean->SetLineWidth(1);
110 h_Mean->SetLineColor(kBlack);
111 map_VarToMeanHist[vec_Plots.at(
j).fName] = h_Mean;
119 TLegend*
leg =
new TLegend(0.17,0.17,0.83,0.4);
120 leg->SetFillColor(0);
122 TLegend* leg_1D_1 =
new TLegend(0.6,0.6,0.85,0.85);
123 leg_1D_1->SetFillColor(0);
124 leg_1D_1->SetNColumns(3);
125 TLegend* leg_1D_2 =
new TLegend(0.35,0.20,0.62,0.48);
126 leg_1D_2->SetFillColor(0);
127 leg_1D_2->SetNColumns(3);
129 for(
auto spectrum : map_PeriodAndVarToSpectrum)
131 TH2 *h_2D = spectrum.second->ToTH2(spectrum.second->POT());
132 TProfile *h_Prof = h_2D->ProfileX((TString)(
"prof_"+vec_Plots.at(spectrum.first.second).fName+
"_"+Form(
"%i_%i",
datasets[spectrum.first.first].fYear,
datasets[spectrum.first.first].fMonth)));
133 h_Prof->SetMarkerStyle(
datasets[spectrum.first.first].fStyle);
134 h_Prof->SetMarkerColor(
datasets[spectrum.first.first].fColor);
135 h_Prof->SetLineColor (
datasets[spectrum.first.first].fColor);
136 h_Prof->SetMarkerSize (1.4);
137 if(spectrum.first.first==0)
139 map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName]->cd();
140 h_Prof->SetTitle(
"");
141 h_Prof->GetXaxis()->SetTitle((TString)spectrum.second->GetLabels().at(0));
142 h_Prof->GetYaxis()->SetTitle((TString)spectrum.second->GetLabels().at(1));
143 h_Prof->GetXaxis()->CenterTitle();
144 h_Prof->GetYaxis()->CenterTitle();
145 h_Prof->GetXaxis()->SetTimeDisplay(1);
146 h_Prof->GetXaxis()->SetTimeFormat(
"%y-%m-%d%F1970-01-01 00:00:00");
149 unsigned int minBin = h_Prof->GetXaxis()->FindBin(startTime);
150 unsigned int maxBin = h_Prof->GetXaxis()->FindBin(endTime);
151 double minZoom = h_Prof->GetXaxis()->GetBinLowEdge(minBin-1);
152 double maxZoom = h_Prof->GetXaxis()->GetBinLowEdge(maxBin+2);
153 h_Prof->GetXaxis()->SetRangeUser(minZoom, maxZoom);
154 h_Prof->GetYaxis()->SetRangeUser(vec_Plots.at(spectrum.first.second).fYRangeLow, vec_Plots.at(spectrum.first.second).fYRangeHigh);
158 map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName]->cd();
159 h_Prof->Draw(
"SAME");
162 if(spectrum.first.second==0)
164 leg->AddEntry(h_Prof, (TString)
datasets[spectrum.first.first].fLabel,
"P");
165 leg_1D_1->AddEntry(h_Prof, (TString)datasets[spectrum.first.first].fLabel,
"P");
166 leg_1D_2->AddEntry(h_Prof, (TString)datasets[spectrum.first.first].fLabel,
"P");
168 if(spectrum.first.first==
datasets.size()-1)
174 TH1D *h_Dummy =
new TH1D(
"h_Dummy",
"h_Dummy",1,0,1);
175 h_Dummy->SetLineColor(kBlack);
176 leg_1D_1->AddEntry(h_Dummy,
"Average (Mean)",
"L");
177 leg_1D_2->AddEntry(h_Dummy,
"Average (Mean)",
"L");
180 for(
auto spectrum : map_PeriodAndVarToSpectrum)
182 TH2 *h_2D = spectrum.second->ToTH2(spectrum.second->POT());
183 TH1D *h_1D = h_2D->ProjectionY((TString)(
"1D_"+vec_Plots.at(spectrum.first.second).fName+
"_"+Form(
"%i_%i",
datasets[spectrum.first.first].fYear,
datasets[spectrum.first.first].fMonth)));
184 if(spectrum.first.first==0)
186 map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName] = (TH1D*)h_1D->Clone((TString)(
"h_Mean"+vec_Plots.at(spectrum.first.second).fName));
190 map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->Add(h_1D);
193 if(vec_hPOT.at(spectrum.first.first)->Integral()!=0.){h_1D->Scale(1.0/vec_hPOT.at(spectrum.first.first)->Integral());}
194 TGraphErrors *gr =
new TGraphErrors(h_1D);
195 gr->SetMarkerStyle(
datasets[spectrum.first.first].fStyle);
196 gr->SetMarkerColor(
datasets[spectrum.first.first].fColor);
197 gr->SetLineColor (
datasets[spectrum.first.first].fColor);
199 for(
int j = 0;
j < gr->GetN();
j++)
202 double shift = (h_1D->GetBinWidth(1)/(double)
datasets.size())*((spectrum.first.first+1)*1. -
datasets.size()*0.5);
203 gr->GetPoint(
j, x, y);
204 dy = gr->GetErrorY(
j);
205 gr->SetPoint(
j, x + shift, y);
206 gr->SetPointError(
j, 0, dy);
208 if(spectrum.first.first==0)
211 gr->GetXaxis()->SetTitle((TString)spectrum.second->GetLabels().at(1));
212 gr->GetYaxis()->SetTitle(
"Slices / 10^{15} POT");
213 gr->GetXaxis()->CenterTitle();
214 gr->GetYaxis()->CenterTitle();
215 map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+
"_1D"]->cd();
220 map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+
"_1D"]->cd();
224 if(spectrum.first.first==
datasets.size()-1)
226 map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->Scale(1.0e15/map_PeriodAndVarToSpectrum[{0, spectrum.first.second}]->POT());
227 map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->SetTitle(
"");
228 map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->Draw(
"HIST SAME");
229 if(vec_Plots.at(spectrum.first.second).fName==
"ccE" ||
230 vec_Plots.at(spectrum.first.second).fName==
"muonE" ||
231 vec_Plots.at(spectrum.first.second).fName==
"hadE" ||
232 vec_Plots.at(spectrum.first.second).fName==
"hadEfrac" ||
233 vec_Plots.at(spectrum.first.second).fName==
"nhit") {
244 for(
auto spectrum : map_PeriodAndVarToSpectrum)
246 TH2 *h_2D = spectrum.second->ToTH2(spectrum.second->POT());
247 TH1D *h_1D = h_2D->ProjectionY((TString)(
"1D_"+vec_Plots.at(spectrum.first.second).fName+
"_"+Form(
"%i_%i",
datasets[spectrum.first.first].fYear,
datasets[spectrum.first.first].fMonth)));
248 if(vec_hPOT.at(spectrum.first.first)->Integral()!=0.){h_1D->Scale(1.0/vec_hPOT.at(spectrum.first.first)->Integral());}
249 TGraphErrors *gr =
new TGraphErrors(h_1D);
250 gr->SetMarkerStyle(
datasets[spectrum.first.first].fStyle);
251 gr->SetMarkerColor(
datasets[spectrum.first.first].fColor);
252 gr->SetLineColor (
datasets[spectrum.first.first].fColor);
254 for(
int j = 0;
j < gr->GetN();
j++)
257 double shift = (h_1D->GetBinWidth(1)/(double)
datasets.size())*((spectrum.first.first+1)*1. -
datasets.size()*0.5);
258 my = map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName]->GetBinContent(
j+1);
259 gr->GetPoint(
j, x, y);
263 dy = gr->GetErrorY(
j)/my;
264 gr->SetPoint(
j, x + shift, y);
265 gr->SetPointError(
j, 0, dy);
269 if(spectrum.first.first==0)
271 map_VarToMeanHist[vec_Plots.at(spectrum.first.second).fName] = (TH1D*)h_1D->Clone((TString)(
"h_Mean"+vec_Plots.at(spectrum.first.second).fName));
273 gr->GetXaxis()->SetTitle((TString)spectrum.second->GetLabels().at(1));
274 gr->GetYaxis()->SetTitle(
"Slices per 10^{15} POT / Average");
275 gr->GetXaxis()->CenterTitle();
276 gr->GetYaxis()->CenterTitle();
277 map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+
"_Ratio"]->cd();
279 gr->GetYaxis()->SetRangeUser(0.8,1.2);
281 double binsArrary[h_1D->GetNbinsX()];
282 h_1D->GetLowEdge(binsArrary);
283 for(
int j = 0;
j < h_1D->GetNbinsX();
j++)
285 TLine *
l =
new TLine(binsArrary[
j], 0.8, binsArrary[j], 1.2);
286 l->SetLineStyle(kDashed);
292 map_VarToCanvas[vec_Plots.at(spectrum.first.second).fName+
"_Ratio"]->cd();
296 if(spectrum.first.first==
datasets.size()-1)
299 TGraph *xAxis =
new TGraph;
300 xAxis->SetPoint(0, -1.e99, 1);
301 xAxis->SetPoint(1, 1.e99, 1);
306 for(
auto it : map_VarToCanvas)
308 it.second->SaveAs((TString)(outDir+
it.first+
"_LongTerm.png"));
309 if (((
std::string)(
it.second->GetName())).find(
"Profile")!=std::string::npos){dir_Profile->cd();}
310 else if(((
std::string)(
it.second->GetName())).find(
"1D") !=std::string::npos){dir_1D ->cd();}
311 else if(((
std::string)(
it.second->GetName())).find(
"Ratio") !=std::string::npos){dir_Ratio ->cd();}
Represent the binning of a Spectrum's x-axis.
Cuts and Vars for the 2020 FD DiF Study.
const SpillVar kSpillTime([](const caf::SRSpillProxy *spill){return spill->spilltimesec;})
std::vector< Plot > vec_Plots
std::vector< Dataset > datasets
const Var kTrkStartY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Y()/100;})
void SetSpillCut(const SpillCut &cut)
Representation of a spectrum in any variable, with associated POT.
void lastUpdated(unsigned int timestamp)
virtual void AddSpillHistogram(TH1 *h, const SpillVar &var, const SpillCut &cut, const SpillVar &wei=kSpillUnweighted)
Uses include counting the total POT or spills in a run.
const Var kTrkStartZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Z()/100;})
SpillCut getYearAndMonth(unsigned int iYear, unsigned int iMonth)
virtual void Go() override
Load all the registered spectra.
std::vector< float > Spectrum
const Var kTrkStartX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.X()/100;})
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const SpillVar kPOTe15([](const caf::SRSpillProxy *spill){return spill->spillpot/1e15;})
const SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
void drawLongTerm_2019(std::string samDef, unsigned int endTime, unsigned int nDays, std::string outDir)