4 #include "TGraphAsymmErrors.h" 7 #include "TTreeReader.h" 8 #include "TTreeReaderValue.h" 26 int nBins = hist->GetNbinsX();
27 bool foundPeak =
false;
28 bool foundLow =
false;
29 bool foundHigh =
false;
30 int maxBin, lBin, hBin;
36 maxBin = hist->GetMaximumBin();
37 peak = hist->GetBinContent(maxBin);
40 (hist->GetBinContent(maxBin-1) > 0.25 * peak ||
41 hist->GetBinContent(maxBin+1) > 0.25 * peak ))
46 hist->SetBinContent(maxBin,0);
52 lBin = hist->FindFirstBinAbove(peak/2);
55 (hist->GetBinContent(lBin-1) > 0.125 * peak ||
56 hist->GetBinContent(lBin+1) > 0.125 * peak ))
61 hist->SetBinContent(lBin,0);
68 hBin = hist->FindLastBinAbove(peak/2);
71 (hist->GetBinContent(hBin-1) > 0.125 * peak &&
72 hist->GetBinContent(hBin+1) > 0.125 * peak ))
77 hist->SetBinContent(hBin,0);
82 fwhm = hist->GetBinCenter(hBin) - hist->GetBinCenter(lBin);
83 return std::tie(maxBin, peak, lBin, hBin, fwhm);
89 TH1::SetDefaultSumw2(
true);
97 if (found != std::string::npos) det =
"ND";
99 if (found != std::string::npos) det =
"FD";
102 if (found != std::string::npos) cur =
"FHC";
104 if (found != std::string::npos) cur =
"RHC";
112 double binWidth = (double) 1.0 / invBinWidth;
113 int nLowBins = (
int) invBinWidth * highCalE;
114 if (highCalE >= 20.0) nLowBins = invBinWidth * 20 - 1;
115 double binEdges[nLowBins + 2];
116 for (
int binNum = 0; binNum < nLowBins + 1; binNum++)
118 binEdges[binNum] = binWidth * binNum;
120 binEdges[nLowBins + 1] = 20.0;
129 TTree* inTree = (TTree*) inFile->Get(
"outTree");
130 TTreeReader inReader(
"outTree", inFile);
133 TTreeReaderValue<bool> isNC (inReader,
"isNC");
134 TTreeReaderValue<bool> isInDet (inReader,
"isInDet");
135 TTreeReaderValue<bool> passedNDCuts(inReader,
"passedNDCuts");
136 TTreeReaderValue<bool> passedFDCuts(inReader,
"passedFDCuts");
137 TTreeReaderValue<double> eCal (inReader,
"eCal");
138 TTreeReaderValue<double> eDep (inReader,
"eDep");
139 TTreeReaderValue<double> eOrph (inReader,
"nonPngE");
140 TTreeReaderValue<double> eHad (inReader,
"eHad");
141 TTreeReaderValue<double> eEM (inReader,
"eEM");
142 TTreeReaderValue<double> ppfxWgt (inReader,
"ppfxWgt");
143 TTreeReaderValue<int>
mode (inReader,
"mode");
144 TTreeReaderValue<unsigned int> nPiZero (inReader,
"nPiZero");
145 TTreeReaderValue<unsigned int> nPiMinus (inReader,
"nPiMinus");
146 TTreeReaderValue<unsigned int> nPiPlus (inReader,
"nPiPlus");
147 TTreeReaderValue<unsigned int> nNeutron (inReader,
"nNeutron");
149 TTreeReaderValue<bool> passedCuts = (det ==
"ND") ? passedNDCuts : passedFDCuts;
153 TFile*
outFile =
new TFile(outFileName.c_str(),
"RECREATE");
157 double totEvents = 0;
160 std::string eCalSpecTitle = title +
";E_{Cal} (GeV);Events";
161 std::string eDepSpecTitle = title +
";y #times E_{#nu} (GeV);Events";
162 TH1D* eCalSpec =
new TH1D(
"eCalSpec", eCalSpecTitle.c_str(), nLowBins + 1, binEdges);
163 TH1D* eDepSpec =
new TH1D(
"eDepSpec", eDepSpecTitle.c_str(), nLowBins + 1, binEdges);
166 std::string emLikeCalTitle = title +
"#; At Least One #pi^{0}_{}#;No #pi^{-}_{}, #pi^{+}_{}, or Neutron;(E_{Cal} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu});Events";
167 std::string hadLikeCalTitle = title +
"#; At Least One #pi^{-}_{}, #pi^{+}_{}, or Neutron#;No #pi^{0}_{};(E_{Cal} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu});Events";
168 TH1D* emLikeCal =
new TH1D(
"emLikeCal" , emLikeCalTitle .c_str(), 200, -1, 1);
169 TH1D* hadLikeCal =
new TH1D(
"hadLikeCal", hadLikeCalTitle.c_str(), 200, -1, 1);
172 while (inReader.Next())
175 if (not (*isNC && *isInDet && *passedCuts && (*eCal>0.5)))
continue;
178 totEvents += *ppfxWgt;
181 eCalSpec->Fill(*eCal, *ppfxWgt);
182 eDepSpec->Fill(*eDep, *ppfxWgt);
183 if (*nPiZero > 0 && (*nPiMinus + *nPiPlus + *nNeutron) == 0) emLikeCal ->Fill((*eCal - *eDep) / *eDep, *ppfxWgt);
184 if (*nPiZero == 0 && (*nPiMinus + *nPiPlus + *nNeutron) > 0) hadLikeCal->Fill((*eCal - *eDep) / *eDep, *ppfxWgt);
194 std::string eCalRunFracTitle = title +
";E_{Cal} (GeV);Cumulative Fraction of Events";
195 std::string eDepRunFracTitle = title +
";y #times E_{#nu} (GeV);Cumulative Fraction of Events";
196 TH1D* eCalRunFrac =
new TH1D(
"eCalRunFrac", eCalRunFracTitle.c_str(), nLowBins + 1, binEdges);
197 TH1D* eDepRunFrac =
new TH1D(
"eDepRunFrac", eDepRunFracTitle.c_str(), nLowBins + 1, binEdges);
198 double runSumCal = 0;
199 double runSumDep = 0;
200 for (
int bin = 1;
bin <= nLowBins + 1;
bin++)
202 runSumCal += eCalSpec->GetBinContent(
bin);
203 runSumDep += eDepSpec->GetBinContent(
bin);
204 eCalRunFrac->SetBinContent(
bin, runSumCal);
205 eDepRunFrac->SetBinContent(
bin, runSumDep);
207 eCalRunFrac->Scale(1.0 / totEvents);
208 eDepRunFrac->Scale(1.0 / totEvents);
213 double emCalPeak = emLikeCal ->GetMaximum();
214 double hadCalPeak = hadLikeCal->GetMaximum();
215 double emLEdge = emLikeCal ->GetBinCenter(emLikeCal ->FindFirstBinAbove(0.75*emCalPeak ));
216 double hadLEdge = hadLikeCal->GetBinCenter(hadLikeCal->FindFirstBinAbove(0.75*hadCalPeak));
217 double emHEdge = emLikeCal ->GetBinCenter(emLikeCal ->FindLastBinAbove (0.75*emCalPeak ));
218 double hadHEdge = hadLikeCal->GetBinCenter(hadLikeCal->FindLastBinAbove (0.75*hadCalPeak));
219 TFitResultPtr emFit = emLikeCal ->Fit(
"gaus",
"SQO",
"",emLEdge ,emHEdge );
220 TFitResultPtr hadFit = hadLikeCal->Fit(
"gaus",
"SQO",
"",hadLEdge,hadHEdge);
221 double emScaleFactor = 1.0 / (1 + emFit ->Parameter(1));
222 double hadScaleFactor = 1.0 / (1 + hadFit->Parameter(1));
223 TF2* scaleFunc =
new TF2(
"scaleFunc",
"[0]*x + [1]*y");
224 scaleFunc->SetParName(0,
"hadScaleFactor");
225 scaleFunc->SetParName(1,
"emScaleFactor");
226 scaleFunc->SetParameters(hadScaleFactor, emScaleFactor);
227 scaleFunc->Write(
"scaleFunc");
230 std::string scaleCorVsCalTitle = title +
"#; Scaling Correction;E_{Cal} (GeV); (E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})";
231 std::string scaleCorVsDepTitle = title +
"#; Scaling Correction;y #times E_{#nu} (GeV); (E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})";
232 TH2D* scaleCorVsCal =
new TH2D(
"scaleCorVsCal" , scaleCorVsCalTitle.c_str(), nLowBins + 1, binEdges, 200, -1, 1);
233 TH2D* scaleCorVsDep =
new TH2D(
"scaleCorVsDep" , scaleCorVsDepTitle.c_str(), nLowBins + 1, binEdges, 200, -1, 1);
235 while (inReader.Next())
238 if (not (*isNC && *isInDet && *passedCuts && (*eCal>0.5)))
continue;
241 double eCorScale = hadScaleFactor * *eHad + emScaleFactor * *eEM;
242 double resScale = (eCorScale - *eDep) / *eDep;
245 scaleCorVsCal ->Fill(*eCal, resScale , *ppfxWgt);
246 scaleCorVsDep ->Fill(*eDep, resScale , *ppfxWgt);
248 TH1D* scaleCor = (TH1D*) scaleCorVsCal ->
ProjectionY(
"scaleCor");
249 scaleCor->GetYaxis()->SetTitle(
"Events");
255 TH1D*
temp =
new TH1D(
"temp",
"", 200, -1, 1);
256 TH1D* calEBinDists[nLowBins + 1];
257 TH1D* depEBinDists[nLowBins + 1];
258 std::string scaleCalBiasTitle = title +
"#; Scaling Correction;E_{Cal} (GeV); (E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})";
259 std::string scaleDepBiasTitle = title +
"#; Scaling Correction;y #times^{} E_{#nu} (GeV); (E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})";
260 TH1D* scaleCalBias =
new TH1D(
"scaleCalBias", scaleCalBiasTitle.c_str(), nLowBins + 1, binEdges);
261 TH1D* scaleDepBias =
new TH1D(
"scaleDepBias", scaleDepBiasTitle.c_str(), nLowBins + 1, binEdges);
267 for (
int eBin = 1; eBin <= nLowBins + 1; eBin++)
270 temp = (TH1D*) scaleCorVsCal->ProjectionY(calTempTitle.c_str(), eBin, eBin);
271 calEBinDists[eBin - 1] = (TH1D*) scaleCorVsCal->ProjectionY(calTempTitle.c_str(), eBin, eBin);
273 calEBinDists[eBin - 1]->SetTitle(calBinTitle.c_str());
274 calEBinDists[eBin - 1]->GetXaxis()->SetTitle(
"(E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})");
275 calEBinDists[eBin - 1]->GetYaxis()->SetTitle(
"Events");
276 double tempN = temp->Integral();
277 if (tempN < 1) tempN = 1;
279 double fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.5*peak));
280 double fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.5*peak));
288 fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.8*peak));
289 fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.8*peak));
292 TFitResultPtr fitResult = calEBinDists[eBin - 1]->Fit(
"gaus",
"SO",
"",fitL,fitH);
293 if (calEBinDists[eBin - 1]->
Integral() > 0 && (not fitResult->IsEmpty()))
295 double fitMean = fitResult->Parameter(1);
296 double fitMeanErr = fitResult->ParError (1);
297 scaleCalBias->SetBinContent(eBin, fitMean);
298 scaleCalBias->SetBinError (eBin, fitMeanErr);
301 temp = (TH1D*) scaleCorVsDep->ProjectionY(depTempTitle.c_str(), eBin, eBin);
302 depEBinDists[eBin - 1] = (TH1D*) scaleCorVsDep->ProjectionY(depTempTitle.c_str(), eBin, eBin);
304 depEBinDists[eBin - 1]->SetTitle(depBinTitle.c_str());
305 depEBinDists[eBin - 1]->GetXaxis()->SetTitle(
"(E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})");
306 depEBinDists[eBin - 1]->GetYaxis()->SetTitle(
"Events");
307 tempN = temp->Integral();
308 if (tempN < 1) tempN = 1;
310 fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.5*peak));
311 fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.5*peak));
319 fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.8*peak));
320 fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.8*peak));
322 fitResult = depEBinDists[eBin - 1]->Fit(
"gaus",
"SQO",
"",fitL,fitH);
323 if (depEBinDists[eBin - 1]->
Integral() > 0 && (not fitResult->IsEmpty()))
325 double fitMean = fitResult->Parameter(1);
326 double fitMeanErr = fitResult->ParError (1);
327 scaleDepBias->SetBinContent(eBin, fitMean);
328 scaleDepBias->SetBinError (eBin, fitMeanErr);
339 TF1* biasFit =
new TF1(
"biasFit",
"[0] * std::pow(x + [1],2) * exp(-[2]*x) + [3]");
340 biasFit->SetParameters(1,1,1,-1);
341 scaleCalBias->Fit(
"biasFit",
"",
"", 0.5, 20);
342 TF1* secondOrderFunc =
new TF1(
"secondOrderFunc",
"1 / (1 + [0] * std::pow(x + [1],2) * exp(-[2]*x) + [3])");
343 secondOrderFunc->SetParameters(biasFit->GetParameters());
344 biasFit->Write(
"biasFit");
345 secondOrderFunc->Write(
"secondOrderFunc");
348 std::string scaleCorVsCal2Title = title +
"#; Scaling with Second Order Bias Correction;E_{Cal} (GeV); (E_{Cor,2} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})";
349 std::string scaleCorVsDep2Title = title +
"#; Scaling with Second Order Bias Correction;y #times E_{#nu} (GeV); (E_{Cor,2} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})";
350 TH2D* scaleCorVsCal2 =
new TH2D(
"scaleCorVsCal2", scaleCorVsCal2Title.c_str(), nLowBins + 1, binEdges, 200, -1, 1);
351 TH2D* scaleCorVsDep2 =
new TH2D(
"scaleCorVsDep2", scaleCorVsDep2Title.c_str(), nLowBins + 1, binEdges, 200, -1, 1);
353 while (inReader.Next())
356 if (not (*isNC && *isInDet && *passedCuts && (*eCal>0.5)))
continue;
359 double eCorScale = hadScaleFactor * *eHad + emScaleFactor * *eEM;
360 double eCor2ndOrd = secondOrderFunc->Eval(*eCal) * eCorScale;
361 double resScale = (eCor2ndOrd - *eDep) / *eDep;
364 scaleCorVsCal2 ->Fill(*eCal, resScale , *ppfxWgt);
365 scaleCorVsDep2 ->Fill(*eDep, resScale , *ppfxWgt);
367 TH1* scaleCor2 = (TH1D*) scaleCorVsCal2 ->
ProjectionY(
"scaleCor2");
368 scaleCor2 ->GetYaxis()->SetTitle(
"Events");
374 TH1D* calEBinDists2[nLowBins + 1];
375 TH1D* depEBinDists2[nLowBins + 1];
376 std::string scaleCalBias2Title = title +
"#; Scaling with Second Order Bias Correction;E_{Cal} (GeV); (E_{Cor,2} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})";
377 std::string scaleDepBias2Title = title +
"#; Scaling with Second Order Bias Correction;y #times^{} E_{#nu} (GeV); (E_{Cor,2} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})";
378 TH1D* scaleCalBias2 =
new TH1D(
"scaleCalBias2", scaleCalBias2Title.c_str(), nLowBins + 1, binEdges);
379 TH1D* scaleDepBias2 =
new TH1D(
"scaleDepBias2", scaleDepBias2Title.c_str(), nLowBins + 1, binEdges);
380 for (
int eBin = 1; eBin <= nLowBins + 1; eBin++)
383 temp = (TH1D*) scaleCorVsCal2->ProjectionY(calTempTitle.c_str(), eBin, eBin);
384 calEBinDists2[eBin - 1] = (TH1D*) scaleCorVsCal2->ProjectionY(calTempTitle.c_str(), eBin, eBin);
386 calEBinDists2[eBin - 1]->SetTitle(calBinTitle.c_str());
387 calEBinDists2[eBin - 1]->GetXaxis()->SetTitle(
"(E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})");
388 calEBinDists2[eBin - 1]->GetYaxis()->SetTitle(
"Events");
389 double tempN = temp->Integral();
390 if (tempN < 1) tempN = 1;
392 double fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.5*peak));
393 double fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.5*peak));
401 fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.8*peak));
402 fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.8*peak));
404 std::cout <<
"-----Bias Fit For Bin " << eBin <<
" After Bias Corr-----" <<
std::endl;
405 TFitResultPtr fitResult = calEBinDists2[eBin - 1]->Fit(
"gaus",
"SO",
"",fitL,fitH);
406 if (calEBinDists2[eBin - 1]->
Integral() > 0 && (not fitResult->IsEmpty()))
408 double fitMean = fitResult->Parameter(1);
409 double fitMeanErr = fitResult->ParError (1);
410 scaleCalBias2->SetBinContent(eBin, fitMean);
411 scaleCalBias2->SetBinError (eBin, fitMeanErr);
414 temp = (TH1D*) scaleCorVsDep2->ProjectionY(depTempTitle.c_str(), eBin, eBin);
415 depEBinDists2[eBin - 1] = (TH1D*) scaleCorVsDep2->ProjectionY(depTempTitle.c_str(), eBin, eBin);
417 depEBinDists2[eBin - 1]->SetTitle(depBinTitle.c_str());
418 depEBinDists2[eBin - 1]->GetXaxis()->SetTitle(
"(E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})");
419 depEBinDists2[eBin - 1]->GetYaxis()->SetTitle(
"Events");
420 tempN = temp->Integral();
421 if (tempN < 1) tempN = 1;
423 fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.5*peak));
424 fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.5*peak));
432 fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.8*peak));
433 fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.8*peak));
435 fitResult = depEBinDists2[eBin - 1]->Fit(
"gaus",
"SQO",
"",fitL,fitH);
436 if (depEBinDists2[eBin - 1]->
Integral() > 0 && (not fitResult->IsEmpty()))
438 double fitMean = fitResult->Parameter(1);
439 double fitMeanErr = fitResult->ParError (1);
440 scaleDepBias2->SetBinContent(eBin, fitMean);
441 scaleDepBias2->SetBinError (eBin, fitMeanErr);
double Integral(const Spectrum &s, const double pot, int cvnregion)
std::tuple< int, double, int, int, double > findPeakFWHM_ignoreSpikes(TH1D *hist)
inFileName
if we get here, we're doing the base definitions:
Eigen::ArrayXd ProjectionY(const Eigen::MatrixXd &mat)
Helper for WeightingVariable.
void NuSCalculateCorr(const std::string inFileName="", const double highCalE=20.0)
std::string to_string(ModuleType mt)