NuSCalculateCorr.C
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TF1.h"
3 #include "TFile.h"
4 #include "TGraphAsymmErrors.h"
5 #include "TH2.h"
6 #include "TTree.h"
7 #include "TTreeReader.h"
8 #include "TTreeReaderValue.h"
9 #include "TLegend.h"
10 #include "TProfile.h"
11 
12 #include <iostream>
13 #include <string>
14 
15 // HOW TO USE
16 // run CAFAna/nus/NuSTreeMaker.C over the cafe definition you want use for your correction
17 // this is done with: nohup sh -c "samweb list-definition-files <samDefName> | xargs -n 1 samweb2xrootd | xargs -n 1 cafe NuSTreeMaker.C" >& NuSTreeMaker.log &
18 // hadd the resultant root files. Name the new file with ND/FD and FHC/RHC in the name
19 // this ensures the titles are made correctly
20 // run this macro over the hadded file
21 
22 // find the peak & FWHM of a histogram, ignoring one-off spikes
23 std::tuple<int, double, int, int, double> findPeakFWHM_ignoreSpikes(TH1D* hist)
24 {
25  //std::cout << "......looking for peak..." << std::endl;
26  int nBins = hist->GetNbinsX();
27  bool foundPeak = false;
28  bool foundLow = false;
29  bool foundHigh = false;
30  int maxBin, lBin, hBin;
31  double peak, fwhm;
32 
33  // loop to find true peak
34  for (int n = 0; n < nBins; n++)
35  {
36  maxBin = hist->GetMaximumBin();
37  peak = hist->GetBinContent(maxBin);
38  if ( maxBin == 0 ||
39  maxBin == nBins ||
40  (hist->GetBinContent(maxBin-1) > 0.25 * peak ||
41  hist->GetBinContent(maxBin+1) > 0.25 * peak ))
42  {
43  foundPeak = true;
44  break;
45  }
46  hist->SetBinContent(maxBin,0);
47  }
48  //std::cout << "......found peak, looking for low edge of FWHM..." << std::endl;
49  // loop to find true low bin
50  for (int n = 0; n < nBins; n++)
51  {
52  lBin = hist->FindFirstBinAbove(peak/2);
53  if ( lBin == 0 ||
54  lBin == nBins ||
55  (hist->GetBinContent(lBin-1) > 0.125 * peak ||
56  hist->GetBinContent(lBin+1) > 0.125 * peak ))
57  {
58  foundLow = true;
59  break;
60  }
61  hist->SetBinContent(lBin,0);
62  }
63 
64  //std::cout << "......found looking for low edge of FWHM, looking for high edge of FWHM..." << std::endl;
65  // loop to find true high bin
66  for (int n = 0; n < nBins; n++)
67  {
68 hBin = hist->FindLastBinAbove(peak/2);
69  if ( hBin == 0 ||
70  hBin == nBins ||
71  (hist->GetBinContent(hBin-1) > 0.125 * peak &&
72  hist->GetBinContent(hBin+1) > 0.125 * peak ))
73  {
74  foundHigh = true;
75  break;
76  }
77  hist->SetBinContent(hBin,0);
78  }
79 
80  //std::cout << "......found high edge of FWHM..." << std::endl;
81  //std::cout << "......Done!" << std::endl;
82  fwhm = hist->GetBinCenter(hBin) - hist->GetBinCenter(lBin);
83  return std::tie(maxBin, peak, lBin, hBin, fwhm);
84 }
85 
86 void NuSCalculateCorr(const std::string inFileName = "", const double highCalE = 20.0)
87 {
88  // Set Sumw2
89  TH1::SetDefaultSumw2(true);
90 
91  // Generate title
92  std::string det = "";
93  std::string cur = "";
94  std::size_t found;
95  // Look for detector
96  found = inFileName.find("ND");
97  if (found != std::string::npos) det = "ND";
98  found = inFileName.find("FD");
99  if (found != std::string::npos) det = "FD";
100  // Look for current
101  found = inFileName.find("FHC");
102  if (found != std::string::npos) cur = "FHC";
103  found = inFileName.find("RHC");
104  if (found != std::string::npos) cur = "RHC";
105  std::string title = det + " " + cur;
106 
107  // Set up binning
108  // Go up to highCalE in steps of 0.5
109  // if highCalE is not a integer multiple of 0.5 then round down
110  // if highCalE > 20.0 just do 0.5 GeV bins all the way to 20
111  int invBinWidth = 2; // 1/invBinSize should be the size of the bins
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++)
117  {
118  binEdges[binNum] = binWidth * binNum;
119  }
120  binEdges[nLowBins + 1] = 20.0;
121  //std::cout << "There are " << nLowBins + 1 << " bins whose edges are..." << std::endl;
122  //for (const double edge : binEdges)
123  // {
124  // std::cout << edge << std::endl;
125  // }
126 
127  // Load file and tree
128  TFile* inFile = TFile::Open(inFileName.c_str(), "READ");
129  TTree* inTree = (TTree*) inFile->Get("outTree");
130  TTreeReader inReader("outTree", inFile);
131 
132  // Set used variables
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");
148 
149  TTreeReaderValue<bool> passedCuts = (det == "ND") ? passedNDCuts : passedFDCuts;
150 
151  // Set outFile
152  std::string outFileName = "CorrectionOut.root";
153  TFile* outFile = new TFile(outFileName.c_str(),"RECREATE");
154  outFile->cd();
155 
156  // Set total event counter. Technically the sum of event weights
157  double totEvents = 0;
158 
159  // Set histograms for Ecal and Edep spectra
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);
164 
165  // Set histograms for getting corrections
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);
170 
171  // Loop to fill initial histograms
172  while (inReader.Next())
173  {
174  // Take only selected events
175  if (not (*isNC && *isInDet && *passedCuts && (*eCal>0.5))) continue;
176 
177  // Add event to total
178  totEvents += *ppfxWgt;
179 
180  // Fill histograms
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);
185  }
186 
187  // Save histograms
188  //eCalSpec->Write("eCalSpec");
189  //eDepSpec->Write("eDepSpec");
190  //emLikeCal ->Write("emLikeCal");
191  //hadLikeCal->Write("hadLikeCal");
192 
193  // Make fraction of events per bin histograms
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++)
201  {
202  runSumCal += eCalSpec->GetBinContent(bin);
203  runSumDep += eDepSpec->GetBinContent(bin);
204  eCalRunFrac->SetBinContent(bin, runSumCal);
205  eDepRunFrac->SetBinContent(bin, runSumDep);
206  }
207  eCalRunFrac->Scale(1.0 / totEvents);
208  eDepRunFrac->Scale(1.0 / totEvents);
209  //eCalRunFrac->Write("eCalRunFrac");
210  //eDepRunFrac->Write("eDepRunFrac");
211 
212  // Get corrections
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");
228 
229  // Make res vs eCal and vs eDep
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);
234  inReader.Restart();
235  while (inReader.Next())
236  {
237  // Take only selected events
238  if (not (*isNC && *isInDet && *passedCuts && (*eCal>0.5))) continue;
239 
240  // calculate corrected energy for each method
241  double eCorScale = hadScaleFactor * *eHad + emScaleFactor * *eEM;
242  double resScale = (eCorScale - *eDep) / *eDep;
243 
244  // Fill histograms
245  scaleCorVsCal ->Fill(*eCal, resScale , *ppfxWgt);
246  scaleCorVsDep ->Fill(*eDep, resScale , *ppfxWgt);
247  }
248  TH1D* scaleCor = (TH1D*) scaleCorVsCal ->ProjectionY("scaleCor");
249  scaleCor->GetYaxis()->SetTitle("Events");
250  //scaleCorVsCal ->Write("scaleCorVsCal");
251  //scaleCorVsDep ->Write("scaleCorVsDep");
252  //scaleCor ->Write("scaleCor");
253 
254  // Loop to find biases
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);
262  double peak;
263  double FWHM;
264  int peakBin;
265  int FWHMBinL;
266  int FWHMBinH;
267  for (int eBin = 1; eBin <= nLowBins + 1; eBin++)
268  {
269  std::string calTempTitle = "eCalBin" + std::to_string(eBin);
270  temp = (TH1D*) scaleCorVsCal->ProjectionY(calTempTitle.c_str(), eBin, eBin);
271  calEBinDists[eBin - 1] = (TH1D*) scaleCorVsCal->ProjectionY(calTempTitle.c_str(), eBin, eBin);
272  std::string calBinTitle = title + "#; " + std::to_string(binEdges[eBin - 1]) + " GeV < E_{Cal} < " + std::to_string(binEdges[eBin]) + " GeV";
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;
278  std::tie(peakBin, peak, FWHMBinL, FWHMBinH, FWHM) = findPeakFWHM_ignoreSpikes(temp);
279  double fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.5*peak));
280  double fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.5*peak));
281  if (tempN < 3200)
282  {
283  fitL = -1.0;
284  fitH = 1.0;
285  }
286  if (tempN > 50000)
287  {
288  fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.8*peak));
289  fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.8*peak));
290  }
291  std::cout << "-----Bias Fit For Bin " << eBin << "-----" << std::endl;
292  TFitResultPtr fitResult = calEBinDists[eBin - 1]->Fit("gaus","SO","",fitL,fitH);
293  if (calEBinDists[eBin - 1]->Integral() > 0 && (not fitResult->IsEmpty()))
294  {
295  double fitMean = fitResult->Parameter(1);
296  double fitMeanErr = fitResult->ParError (1);
297  scaleCalBias->SetBinContent(eBin, fitMean);
298  scaleCalBias->SetBinError (eBin, fitMeanErr);
299  }
300  std::string depTempTitle = "eDepBin" + std::to_string(eBin);
301  temp = (TH1D*) scaleCorVsDep->ProjectionY(depTempTitle.c_str(), eBin, eBin);
302  depEBinDists[eBin - 1] = (TH1D*) scaleCorVsDep->ProjectionY(depTempTitle.c_str(), eBin, eBin);
303  std::string depBinTitle = title + "#; " + std::to_string(binEdges[eBin - 1]) + " GeV < y #times E_{#nu} < " + std::to_string(binEdges[eBin]) + " GeV";
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;
309  std::tie(peakBin, peak, FWHMBinL, FWHMBinH, FWHM) = findPeakFWHM_ignoreSpikes(temp);
310  fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.5*peak));
311  fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.5*peak));
312  if (tempN < 3200)
313  {
314  fitL = -1.0;
315  fitH = 1.0;
316  }
317  if (tempN > 50000)
318  {
319  fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.8*peak));
320  fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.8*peak));
321  }
322  fitResult = depEBinDists[eBin - 1]->Fit("gaus","SQO","",fitL,fitH);
323  if (depEBinDists[eBin - 1]->Integral() > 0 && (not fitResult->IsEmpty()))
324  {
325  double fitMean = fitResult->Parameter(1);
326  double fitMeanErr = fitResult->ParError (1);
327  scaleDepBias->SetBinContent(eBin, fitMean);
328  scaleDepBias->SetBinError (eBin, fitMeanErr);
329  }
330  }
331 
332  //calEBinDists ->Write("calEBinDists");
333  //scaleCalBias ->Write("scaleCalBias");
334  //scaleCalBias_fit->Write("scaleCalBias_fit");
335  //scaleDepBias ->Write("scaleDepBias");
336  //scalePeakFWHM ->Write("scalePeakFWHM");
337 
338  // fit bias to obtain second order correction
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");
346 
347  // Get second order corrected energy
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);
352  inReader.Restart();
353  while (inReader.Next())
354  {
355  // Take only selected events
356  if (not (*isNC && *isInDet && *passedCuts && (*eCal>0.5))) continue;
357 
358  // calculate corrected energy for each method
359  double eCorScale = hadScaleFactor * *eHad + emScaleFactor * *eEM;
360  double eCor2ndOrd = secondOrderFunc->Eval(*eCal) * eCorScale;
361  double resScale = (eCor2ndOrd - *eDep) / *eDep;
362 
363  // Fill histograms
364  scaleCorVsCal2 ->Fill(*eCal, resScale , *ppfxWgt);
365  scaleCorVsDep2 ->Fill(*eDep, resScale , *ppfxWgt);
366  }
367  TH1* scaleCor2 = (TH1D*) scaleCorVsCal2 ->ProjectionY("scaleCor2");
368  scaleCor2 ->GetYaxis()->SetTitle("Events");
369  //scaleCorVsCal2 ->Write("scaleCorVsCal2");
370  //scaleCorVsDep2 ->Write("scaleCorVsDep2");
371  //scaleCor2 ->Write("scaleCor2");
372 
373  // loop for biases
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++)
381  {
382  std::string calTempTitle = "eCalBin" + std::to_string(eBin) + "_2";
383  temp = (TH1D*) scaleCorVsCal2->ProjectionY(calTempTitle.c_str(), eBin, eBin);
384  calEBinDists2[eBin - 1] = (TH1D*) scaleCorVsCal2->ProjectionY(calTempTitle.c_str(), eBin, eBin);
385  std::string calBinTitle = title + "#; " + std::to_string(binEdges[eBin - 1]) + " GeV < E_{Cal} < " + std::to_string(binEdges[eBin]) + " GeV#; Second Order Bias Correction";
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;
391  std::tie(peakBin, peak, FWHMBinL, FWHMBinH, FWHM) = findPeakFWHM_ignoreSpikes(temp);
392  double fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.5*peak));
393  double fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.5*peak));
394  if (tempN < 3200)
395  {
396  fitL = -1.0;
397  fitH = 1.0;
398  }
399  if (tempN > 50000)
400  {
401  fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.8*peak));
402  fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.8*peak));
403  }
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()))
407  {
408  double fitMean = fitResult->Parameter(1);
409  double fitMeanErr = fitResult->ParError (1);
410  scaleCalBias2->SetBinContent(eBin, fitMean);
411  scaleCalBias2->SetBinError (eBin, fitMeanErr);
412  }
413  std::string depTempTitle = "eDepBin" + std::to_string(eBin) + "_2";
414  temp = (TH1D*) scaleCorVsDep2->ProjectionY(depTempTitle.c_str(), eBin, eBin);
415  depEBinDists2[eBin - 1] = (TH1D*) scaleCorVsDep2->ProjectionY(depTempTitle.c_str(), eBin, eBin);
416  std::string depBinTitle = title + "#; " + std::to_string(binEdges[eBin - 1]) + " GeV < y #times E_{#nu} < " + std::to_string(binEdges[eBin]) + " GeV#; Second Order Bias Correction";
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;
422  std::tie(peakBin, peak, FWHMBinL, FWHMBinH, FWHM) = findPeakFWHM_ignoreSpikes(temp);
423  fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.5*peak));
424  fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.5*peak));
425  if (tempN < 3200)
426  {
427  fitL = -1.0;
428  fitH = 1.0;
429  }
430  if (tempN > 50000)
431  {
432  fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.8*peak));
433  fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.8*peak));
434  }
435  fitResult = depEBinDists2[eBin - 1]->Fit("gaus","SQO","",fitL,fitH);
436  if (depEBinDists2[eBin - 1]->Integral() > 0 && (not fitResult->IsEmpty()))
437  {
438  double fitMean = fitResult->Parameter(1);
439  double fitMeanErr = fitResult->ParError (1);
440  scaleDepBias2->SetBinContent(eBin, fitMean);
441  scaleDepBias2->SetBinError (eBin, fitMeanErr);
442  }
443  }
444  //scaleCalBias2 ->Write("scaleCalBias2");
445  //scaleDepBias2 ->Write("scaleDepBias2");
446 
447  //write file
448  outFile->Write();
449 
450 }
int nBins
Definition: plotROC.py:16
double Integral(const Spectrum &s, const double pot, int cvnregion)
std::tuple< int, double, int, int, double > findPeakFWHM_ignoreSpikes(TH1D *hist)
ifstream inFile
Definition: AnaPlotMaker.h:34
TFile * outFile
Definition: PlotXSec.C:135
inFileName
if we get here, we&#39;re doing the base definitions:
Definition: mkDefs.py:168
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
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)
Definition: ModuleType.h:32
enum BeamMode string