NuSCalculateCorr.C

Go to the documentation of this file.

#include "TCanvas.h"
#include "TF1.h"
#include "TFile.h"
#include "TGraphAsymmErrors.h"
#include "TH2.h"
#include "TTree.h"
#include "TTreeReader.h"
#include "TTreeReaderValue.h"
#include "TLegend.h"
#include "TProfile.h"

#include <iostream>
#include <string>

// HOW TO USE
// run CAFAna/nus/NuSTreeMaker.C over the cafe definition you want use for your correction
// this is done with: nohup sh -c "samweb list-definition-files <samDefName> | xargs -n 1 samweb2xrootd | xargs -n 1 cafe NuSTreeMaker.C" >& NuSTreeMaker.log &
// hadd the resultant root files. Name the new file with ND/FD and FHC/RHC in the name
// this ensures the titles are made correctly
// run this macro over the hadded file

// find the peak & FWHM of a histogram, ignoring one-off spikes
std::tuple<int, double, int, int, double> findPeakFWHM_ignoreSpikes(TH1D* hist)
{
  //std::cout << "......looking for peak..." << std::endl;
  int   nBins = hist->GetNbinsX();
  bool  foundPeak = false;
  bool  foundLow  = false;
  bool  foundHigh = false;
  int   maxBin, lBin, hBin;
  double peak, fwhm;

  // loop to find true peak
  for (int n = 0; n < nBins; n++)
    {
      maxBin = hist->GetMaximumBin();
      peak   = hist->GetBinContent(maxBin);
      if ( maxBin == 0                                 ||
           maxBin == nBins                             ||
           (hist->GetBinContent(maxBin-1) > 0.25 * peak ||
            hist->GetBinContent(maxBin+1) > 0.25 * peak ))
        {
          foundPeak = true;
          break;
        }
      hist->SetBinContent(maxBin,0);
    }
  //std::cout << "......found peak, looking for low edge of FWHM..." << std::endl;
  // loop to find true low bin
  for (int n = 0; n < nBins; n++)
    {
      lBin = hist->FindFirstBinAbove(peak/2);
      if ( lBin == 0                                   ||
           lBin == nBins                               ||
           (hist->GetBinContent(lBin-1) > 0.125 * peak ||
            hist->GetBinContent(lBin+1) > 0.125 * peak ))
        {
          foundLow = true;
          break;
        }
      hist->SetBinContent(lBin,0);
    }

  //std::cout << "......found looking for low edge of FWHM, looking for high edge of FWHM..." << std::endl;
  // loop to find true high bin
  for (int n = 0; n < nBins; n++)
    {
hBin = hist->FindLastBinAbove(peak/2);
      if ( hBin == 0                                   ||
           hBin == nBins                               ||
           (hist->GetBinContent(hBin-1) > 0.125 * peak &&
            hist->GetBinContent(hBin+1) > 0.125 * peak ))
        {
          foundHigh = true;
          break;
        }
      hist->SetBinContent(hBin,0);
    }

  //std::cout << "......found high edge of FWHM..." << std::endl;
  //std::cout << "......Done!" << std::endl;
  fwhm = hist->GetBinCenter(hBin) - hist->GetBinCenter(lBin);
  return std::tie(maxBin, peak, lBin, hBin, fwhm);
}

void NuSCalculateCorr(const std::string inFileName = "", const double highCalE = 20.0)
{
  // Set Sumw2
  TH1::SetDefaultSumw2(true);

  // Generate title
  std::string det = "";
  std::string cur = "";
  std::size_t found;
  // Look for detector
  found = inFileName.find("ND");
  if (found != std::string::npos) det = "ND";
  found = inFileName.find("FD");
  if (found != std::string::npos) det = "FD";
  // Look for current
  found = inFileName.find("FHC");
  if (found != std::string::npos) cur = "FHC";
  found = inFileName.find("RHC");
  if (found != std::string::npos) cur = "RHC";  
  std::string title = det + " " + cur;

  // Set up binning
  // Go up to highCalE in steps of 0.5
  // if highCalE is not a integer multiple of 0.5 then round down
  // if highCalE > 20.0 just do 0.5 GeV bins all the way to 20
  int    invBinWidth = 2; // 1/invBinSize should be the size of the bins
  double binWidth    = (double) 1.0 / invBinWidth;
  int nLowBins = (int) invBinWidth * highCalE;
  if (highCalE >= 20.0) nLowBins = invBinWidth * 20 - 1;
  double binEdges[nLowBins + 2];
  for (int binNum = 0; binNum < nLowBins + 1; binNum++)
    {
      binEdges[binNum] = binWidth * binNum;
    }
  binEdges[nLowBins + 1] = 20.0;
  //std::cout << "There are " << nLowBins + 1 << " bins whose edges are..." << std::endl;
  //for (const double edge : binEdges)
  //  {
  //    std::cout << edge << std::endl;
  //  }

  // Load file and tree
  TFile* inFile = TFile::Open(inFileName.c_str(),  "READ");
  TTree* inTree = (TTree*) inFile->Get("outTree");
  TTreeReader inReader("outTree", inFile);

  // Set used variables
  TTreeReaderValue<bool>         isNC        (inReader,         "isNC");
  TTreeReaderValue<bool>         isInDet     (inReader,      "isInDet");
  TTreeReaderValue<bool>         passedNDCuts(inReader, "passedNDCuts");
  TTreeReaderValue<bool>         passedFDCuts(inReader, "passedFDCuts");
  TTreeReaderValue<double>       eCal        (inReader,         "eCal");
  TTreeReaderValue<double>       eDep        (inReader,         "eDep");
  TTreeReaderValue<double>       eOrph       (inReader,      "nonPngE");
  TTreeReaderValue<double>       eHad        (inReader,         "eHad");
  TTreeReaderValue<double>       eEM         (inReader,          "eEM");
  TTreeReaderValue<double>       ppfxWgt     (inReader,      "ppfxWgt");
  TTreeReaderValue<int>          mode        (inReader,         "mode");
  TTreeReaderValue<unsigned int> nPiZero     (inReader,      "nPiZero");
  TTreeReaderValue<unsigned int> nPiMinus    (inReader,     "nPiMinus");
  TTreeReaderValue<unsigned int> nPiPlus     (inReader,      "nPiPlus");
  TTreeReaderValue<unsigned int> nNeutron    (inReader,     "nNeutron");

  TTreeReaderValue<bool> passedCuts = (det == "ND") ? passedNDCuts : passedFDCuts;

  // Set outFile
  std::string outFileName = "CorrectionOut.root";
  TFile* outFile = new TFile(outFileName.c_str(),"RECREATE");
  outFile->cd();

  // Set total event counter. Technically the sum of event weights
  double totEvents = 0;

  // Set histograms for Ecal and Edep spectra
  std::string eCalSpecTitle = title + ";E_{Cal} (GeV);Events";
  std::string eDepSpecTitle = title + ";y #times E_{#nu} (GeV);Events";
  TH1D* eCalSpec = new TH1D("eCalSpec", eCalSpecTitle.c_str(), nLowBins + 1, binEdges);
  TH1D* eDepSpec = new TH1D("eDepSpec", eDepSpecTitle.c_str(), nLowBins + 1, binEdges);

  // Set histograms for getting corrections
  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";
  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";
  TH1D* emLikeCal  = new TH1D("emLikeCal" , emLikeCalTitle .c_str(), 200, -1, 1);
  TH1D* hadLikeCal = new TH1D("hadLikeCal", hadLikeCalTitle.c_str(), 200, -1, 1);

  // Loop to fill initial histograms
  while (inReader.Next())
    {
       // Take only selected events
       if (not (*isNC && *isInDet && *passedCuts && (*eCal>0.5))) continue;

       // Add event to total
       totEvents += *ppfxWgt;

       // Fill histograms
       eCalSpec->Fill(*eCal, *ppfxWgt);
       eDepSpec->Fill(*eDep, *ppfxWgt);
       if (*nPiZero > 0 && (*nPiMinus + *nPiPlus + *nNeutron) == 0) emLikeCal ->Fill((*eCal - *eDep) / *eDep, *ppfxWgt);
       if (*nPiZero == 0 && (*nPiMinus + *nPiPlus + *nNeutron) > 0) hadLikeCal->Fill((*eCal - *eDep) / *eDep, *ppfxWgt);
    }

  // Save histograms
  //eCalSpec->Write("eCalSpec");
  //eDepSpec->Write("eDepSpec");
  //emLikeCal ->Write("emLikeCal");
  //hadLikeCal->Write("hadLikeCal");

  // Make fraction of events per bin histograms
  std::string eCalRunFracTitle = title + ";E_{Cal} (GeV);Cumulative Fraction of Events";
  std::string eDepRunFracTitle = title + ";y #times E_{#nu} (GeV);Cumulative Fraction of Events";
  TH1D* eCalRunFrac = new TH1D("eCalRunFrac", eCalRunFracTitle.c_str(), nLowBins + 1, binEdges);
  TH1D* eDepRunFrac = new TH1D("eDepRunFrac", eDepRunFracTitle.c_str(), nLowBins + 1, binEdges);
  double runSumCal = 0;
  double runSumDep = 0;
  for (int bin = 1; bin <= nLowBins + 1; bin++)
    {
      runSumCal += eCalSpec->GetBinContent(bin);
      runSumDep += eDepSpec->GetBinContent(bin);
      eCalRunFrac->SetBinContent(bin, runSumCal);
      eDepRunFrac->SetBinContent(bin, runSumDep);
    }
  eCalRunFrac->Scale(1.0 / totEvents);
  eDepRunFrac->Scale(1.0 / totEvents);
  //eCalRunFrac->Write("eCalRunFrac");
  //eDepRunFrac->Write("eDepRunFrac");

  // Get corrections
  double emCalPeak  = emLikeCal ->GetMaximum();
  double hadCalPeak = hadLikeCal->GetMaximum();
  double emLEdge    = emLikeCal ->GetBinCenter(emLikeCal ->FindFirstBinAbove(0.75*emCalPeak ));
  double hadLEdge   = hadLikeCal->GetBinCenter(hadLikeCal->FindFirstBinAbove(0.75*hadCalPeak));
  double emHEdge    = emLikeCal ->GetBinCenter(emLikeCal ->FindLastBinAbove (0.75*emCalPeak ));
  double hadHEdge   = hadLikeCal->GetBinCenter(hadLikeCal->FindLastBinAbove (0.75*hadCalPeak));
  TFitResultPtr emFit  = emLikeCal ->Fit("gaus","SQO","",emLEdge ,emHEdge );
  TFitResultPtr hadFit = hadLikeCal->Fit("gaus","SQO","",hadLEdge,hadHEdge);
  double emScaleFactor  = 1.0 / (1 + emFit ->Parameter(1));
  double hadScaleFactor = 1.0 / (1 + hadFit->Parameter(1));
  TF2* scaleFunc = new TF2("scaleFunc", "[0]*x + [1]*y");
  scaleFunc->SetParName(0, "hadScaleFactor");
  scaleFunc->SetParName(1, "emScaleFactor");
  scaleFunc->SetParameters(hadScaleFactor, emScaleFactor);
  scaleFunc->Write("scaleFunc");

  // Make res vs eCal and vs eDep
  std::string scaleCorVsCalTitle = title + "#; Scaling Correction;E_{Cal} (GeV); (E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})";
  std::string scaleCorVsDepTitle = title + "#; Scaling Correction;y #times E_{#nu} (GeV); (E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})";
  TH2D* scaleCorVsCal  = new TH2D("scaleCorVsCal" , scaleCorVsCalTitle.c_str(), nLowBins + 1, binEdges, 200, -1, 1);
  TH2D* scaleCorVsDep  = new TH2D("scaleCorVsDep" , scaleCorVsDepTitle.c_str(), nLowBins + 1, binEdges, 200, -1, 1);
  inReader.Restart();
  while (inReader.Next())
    {
      // Take only selected events
      if (not (*isNC && *isInDet && *passedCuts && (*eCal>0.5))) continue;

      // calculate corrected energy for each method
      double eCorScale  = hadScaleFactor * *eHad + emScaleFactor * *eEM;
      double resScale   = (eCorScale  - *eDep) / *eDep;

      // Fill histograms
      scaleCorVsCal ->Fill(*eCal, resScale , *ppfxWgt);
      scaleCorVsDep ->Fill(*eDep, resScale , *ppfxWgt);
    }
  TH1D* scaleCor  = (TH1D*) scaleCorVsCal ->ProjectionY("scaleCor");
  scaleCor->GetYaxis()->SetTitle("Events");
  //scaleCorVsCal    ->Write("scaleCorVsCal");
  //scaleCorVsDep    ->Write("scaleCorVsDep");
  //scaleCor         ->Write("scaleCor");

  // Loop to find biases
  TH1D* temp = new TH1D("temp", "", 200, -1, 1);
  TH1D* calEBinDists[nLowBins + 1];
  TH1D* depEBinDists[nLowBins + 1];
  std::string scaleCalBiasTitle = title + "#; Scaling Correction;E_{Cal} (GeV); (E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})";
  std::string scaleDepBiasTitle = title + "#; Scaling Correction;y #times^{} E_{#nu} (GeV); (E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})";
  TH1D* scaleCalBias     = new TH1D("scaleCalBias", scaleCalBiasTitle.c_str(), nLowBins + 1, binEdges);
  TH1D* scaleDepBias     = new TH1D("scaleDepBias", scaleDepBiasTitle.c_str(), nLowBins + 1, binEdges);
  double peak;
  double FWHM;
  int    peakBin;
  int    FWHMBinL;
  int    FWHMBinH;
  for (int eBin = 1; eBin <= nLowBins + 1; eBin++)
    {
      std::string calTempTitle = "eCalBin" + std::to_string(eBin);
      temp = (TH1D*) scaleCorVsCal->ProjectionY(calTempTitle.c_str(), eBin, eBin);
      calEBinDists[eBin - 1] = (TH1D*) scaleCorVsCal->ProjectionY(calTempTitle.c_str(), eBin, eBin);
      std::string calBinTitle = title + "#; " + std::to_string(binEdges[eBin - 1]) + " GeV < E_{Cal} < " + std::to_string(binEdges[eBin]) + " GeV";
      calEBinDists[eBin - 1]->SetTitle(calBinTitle.c_str());
      calEBinDists[eBin - 1]->GetXaxis()->SetTitle("(E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})");
      calEBinDists[eBin - 1]->GetYaxis()->SetTitle("Events");
      double tempN = temp->Integral();
      if (tempN < 1) tempN = 1;
      std::tie(peakBin, peak, FWHMBinL, FWHMBinH, FWHM) = findPeakFWHM_ignoreSpikes(temp);
      double fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.5*peak));
      double fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.5*peak));
      if (tempN < 3200)
        {
          fitL = -1.0;
          fitH =  1.0;
        }
      if (tempN > 50000)
        {
          fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.8*peak));
          fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.8*peak));
        }
      std::cout << "-----Bias Fit For Bin " << eBin << "-----" << std::endl;
      TFitResultPtr fitResult = calEBinDists[eBin - 1]->Fit("gaus","SO","",fitL,fitH);
      if (calEBinDists[eBin - 1]->Integral() > 0 && (not fitResult->IsEmpty()))
        {
          double fitMean    = fitResult->Parameter(1);
          double fitMeanErr = fitResult->ParError (1);
          scaleCalBias->SetBinContent(eBin, fitMean);      
          scaleCalBias->SetBinError  (eBin, fitMeanErr);      
        }
      std::string depTempTitle = "eDepBin" + std::to_string(eBin);
      temp = (TH1D*) scaleCorVsDep->ProjectionY(depTempTitle.c_str(), eBin, eBin);
      depEBinDists[eBin - 1] = (TH1D*) scaleCorVsDep->ProjectionY(depTempTitle.c_str(), eBin, eBin);
      std::string depBinTitle = title + "#; " + std::to_string(binEdges[eBin - 1]) + " GeV < y #times E_{#nu} < " + std::to_string(binEdges[eBin]) + " GeV";
      depEBinDists[eBin - 1]->SetTitle(depBinTitle.c_str());
      depEBinDists[eBin - 1]->GetXaxis()->SetTitle("(E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})");
      depEBinDists[eBin - 1]->GetYaxis()->SetTitle("Events");
      tempN = temp->Integral();
      if (tempN < 1) tempN = 1;
      std::tie(peakBin, peak, FWHMBinL, FWHMBinH, FWHM) = findPeakFWHM_ignoreSpikes(temp);
      fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.5*peak));
      fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.5*peak));
      if (tempN < 3200)
        {
          fitL = -1.0;
          fitH =  1.0;
        }
      if (tempN > 50000)
        {
          fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.8*peak));
          fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.8*peak));
        }
      fitResult = depEBinDists[eBin - 1]->Fit("gaus","SQO","",fitL,fitH);
      if (depEBinDists[eBin - 1]->Integral() > 0 && (not fitResult->IsEmpty()))
        {
          double fitMean    = fitResult->Parameter(1);
          double fitMeanErr = fitResult->ParError (1);
          scaleDepBias->SetBinContent(eBin, fitMean);      
          scaleDepBias->SetBinError  (eBin, fitMeanErr);      
        }
    }

  //calEBinDists   ->Write("calEBinDists");
  //scaleCalBias    ->Write("scaleCalBias");
  //scaleCalBias_fit->Write("scaleCalBias_fit");
  //scaleDepBias    ->Write("scaleDepBias");
  //scalePeakFWHM   ->Write("scalePeakFWHM");

  // fit bias to obtain second order correction
  TF1* biasFit = new TF1("biasFit", "[0] * std::pow(x + [1],2) * exp(-[2]*x) + [3]");
  biasFit->SetParameters(1,1,1,-1);
  scaleCalBias->Fit("biasFit", "", "", 0.5, 20);
  TF1* secondOrderFunc = new TF1("secondOrderFunc", "1 / (1 + [0] * std::pow(x + [1],2) * exp(-[2]*x) + [3])");
  secondOrderFunc->SetParameters(biasFit->GetParameters());
  biasFit->Write("biasFit");
  secondOrderFunc->Write("secondOrderFunc");

  // Get second order corrected energy
  std::string scaleCorVsCal2Title = title + "#; Scaling with Second Order Bias Correction;E_{Cal} (GeV); (E_{Cor,2} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})";
  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})";
  TH2D* scaleCorVsCal2 = new TH2D("scaleCorVsCal2", scaleCorVsCal2Title.c_str(), nLowBins + 1, binEdges, 200, -1, 1);
  TH2D* scaleCorVsDep2 = new TH2D("scaleCorVsDep2", scaleCorVsDep2Title.c_str(), nLowBins + 1, binEdges, 200, -1, 1);
  inReader.Restart();
  while (inReader.Next())
    {
      // Take only selected events
      if (not (*isNC && *isInDet && *passedCuts && (*eCal>0.5))) continue;

      // calculate corrected energy for each method
      double eCorScale    = hadScaleFactor * *eHad + emScaleFactor * *eEM;
      double eCor2ndOrd   = secondOrderFunc->Eval(*eCal) * eCorScale;
      double resScale     = (eCor2ndOrd  - *eDep) / *eDep;

      // Fill histograms
      scaleCorVsCal2  ->Fill(*eCal, resScale  , *ppfxWgt);
      scaleCorVsDep2  ->Fill(*eDep, resScale  , *ppfxWgt);
    }
  TH1* scaleCor2   = (TH1D*) scaleCorVsCal2  ->ProjectionY("scaleCor2");
  scaleCor2 ->GetYaxis()->SetTitle("Events");
  //scaleCorVsCal2  ->Write("scaleCorVsCal2");
  //scaleCorVsDep2  ->Write("scaleCorVsDep2");
  //scaleCor2  ->Write("scaleCor2");

  // loop for biases
  TH1D* calEBinDists2[nLowBins + 1];
  TH1D* depEBinDists2[nLowBins + 1];
  std::string scaleCalBias2Title = title + "#; Scaling with Second Order Bias Correction;E_{Cal} (GeV); (E_{Cor,2} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})";
  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})";
  TH1D* scaleCalBias2 = new TH1D("scaleCalBias2", scaleCalBias2Title.c_str(), nLowBins + 1, binEdges);
  TH1D* scaleDepBias2 = new TH1D("scaleDepBias2", scaleDepBias2Title.c_str(), nLowBins + 1, binEdges);
  for (int eBin = 1; eBin <= nLowBins + 1; eBin++)
    {
      std::string calTempTitle = "eCalBin" + std::to_string(eBin) + "_2";
      temp = (TH1D*) scaleCorVsCal2->ProjectionY(calTempTitle.c_str(), eBin, eBin);
      calEBinDists2[eBin - 1] = (TH1D*) scaleCorVsCal2->ProjectionY(calTempTitle.c_str(), eBin, eBin);
      std::string calBinTitle = title + "#; " + std::to_string(binEdges[eBin - 1]) + " GeV < E_{Cal} < " + std::to_string(binEdges[eBin]) + " GeV#; Second Order Bias Correction";
      calEBinDists2[eBin - 1]->SetTitle(calBinTitle.c_str());
      calEBinDists2[eBin - 1]->GetXaxis()->SetTitle("(E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})");
      calEBinDists2[eBin - 1]->GetYaxis()->SetTitle("Events");
      double tempN = temp->Integral();
      if (tempN < 1) tempN = 1;
      std::tie(peakBin, peak, FWHMBinL, FWHMBinH, FWHM) = findPeakFWHM_ignoreSpikes(temp);
      double fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.5*peak));
      double fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.5*peak));
      if (tempN < 3200)
        {
          fitL = -1.0;
          fitH =  1.0;
        }
      if (tempN > 50000)
        {
          fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.8*peak));
          fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.8*peak));
        }
      std::cout << "-----Bias Fit For Bin " << eBin << " After Bias Corr-----" << std::endl;
      TFitResultPtr fitResult = calEBinDists2[eBin - 1]->Fit("gaus","SO","",fitL,fitH);
      if (calEBinDists2[eBin - 1]->Integral() > 0 && (not fitResult->IsEmpty()))
        {
          double fitMean    = fitResult->Parameter(1);
          double fitMeanErr = fitResult->ParError (1);
          scaleCalBias2->SetBinContent(eBin, fitMean);      
          scaleCalBias2->SetBinError  (eBin, fitMeanErr);      
        }
      std::string depTempTitle = "eDepBin" + std::to_string(eBin) + "_2";
      temp = (TH1D*) scaleCorVsDep2->ProjectionY(depTempTitle.c_str(), eBin, eBin);
      depEBinDists2[eBin - 1] = (TH1D*) scaleCorVsDep2->ProjectionY(depTempTitle.c_str(), eBin, eBin);
      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";
      depEBinDists2[eBin - 1]->SetTitle(depBinTitle.c_str());
      depEBinDists2[eBin - 1]->GetXaxis()->SetTitle("(E_{Cor} - y #times^{} E_{#nu}) / (y #times^{} E_{#nu})");
      depEBinDists2[eBin - 1]->GetYaxis()->SetTitle("Events");
      tempN = temp->Integral();
      if (tempN < 1) tempN = 1;
      std::tie(peakBin, peak, FWHMBinL, FWHMBinH, FWHM) = findPeakFWHM_ignoreSpikes(temp);
      fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.5*peak));
      fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.5*peak));
      if (tempN < 3200)
        {
          fitL = -1.0;
          fitH =  1.0;
        }
      if (tempN > 50000)
        {
          fitL = temp->GetBinCenter(temp->FindFirstBinAbove(0.8*peak));
          fitH = temp->GetBinCenter(temp->FindLastBinAbove (0.8*peak));
        }
      fitResult = depEBinDists2[eBin - 1]->Fit("gaus","SQO","",fitL,fitH);
      if (depEBinDists2[eBin - 1]->Integral() > 0 && (not fitResult->IsEmpty()))
        {
          double fitMean    = fitResult->Parameter(1);
          double fitMeanErr = fitResult->ParError (1);
          scaleDepBias2->SetBinContent(eBin, fitMean);      
          scaleDepBias2->SetBinError  (eBin, fitMeanErr);      
        }
    }
  //scaleCalBias2  ->Write("scaleCalBias2");
  //scaleDepBias2  ->Write("scaleDepBias2");

  //write file
  outFile->Write();

}