NuFluctuator.cxx

Go to the documentation of this file.

// $Id: NuFluctuator.cxx,v 1.10 2009/10/01 18:56:09 bckhouse Exp $
//
// Provides statistically fluctuated mock data sets for CPT and
// transition analyses.
//
// Justin Evans
// j.evans2@physics.ox.ac.uk
//
// $Log: NuFluctuator.cxx,v $
// Revision 1.10  2009/10/01 18:56:09  bckhouse
// Centralize calculation of oscillation functions into NuUtilities: function OscillationWeight, DecayWeight, DecoherenceWeight.
//
// Revision 1.9  2008/02/26 18:53:47  ahimmel
//
//
// Alternat flucutator implemented.  Seems to give the same results.
//
// Revision 1.8  2008/02/04 15:45:54  evans
// A function to make fake transition data.
//
// Revision 1.7  2008/01/23 21:20:55  evans
// The Oscillate() function was using the neutrino sn2 to oscillate the
// antineutrinos. It isn't now.
//
// Revision 1.6  2008/01/15 23:57:31  evans
// Adding the ability to make systematically shifted fake data. Also
// adding IsGoodMajorityCurvature to the list of standard cuts.
//
// Revision 1.5  2008/01/10 19:57:49  ahimmel
//
//
// Broke Fluctuate and Oscillate out into their own methods so they could be used externally.
//
// Revision 1.4  2008/01/04 04:30:08  evans
// Removing unnecessary commented out code.
//
// Revision 1.3  2007/12/25 21:12:56  evans
// In the transition analysis, using the true neutrino energy for the
// efficiency corrections instead of the reconstructed energy.
//
// Revision 1.2  2007/12/24 17:35:02  evans
// The Fluctuator can now produce fake data for a transition analysis. It
// needs to be provided with root files containing the cross sections and
// selection efficiencies.
//
// Revision 1.1  2007/12/18 12:36:02  evans
// Class to produce Poisson fluctuated sets of fake ND and FD data, from
// a high-statistics MC sample. Only produces FD CPT data so far:
// transition functionality will be added presently. Can produce ND
// nominal fake data at present. I'll add in functionality to produce
// systematic shifts.
//
//

#include "TString.h"
#include <vector>

#include "TGraph.h"
#include "TFile.h"
#include "TH1D.h"
#include "TMath.h"
#include "TRandom3.h"

#include "MessageService/MsgService.h"

#include "NtupleUtils/NuEvent.h"
#include "NtupleUtils/NuFluctuator.h"
#include "NtupleUtils/NuTreeWrapper.h"
#include "NtupleUtils/NuXMLConfig.h"
#include "NtupleUtils/NuSystematic.h"
#include "NtupleUtils/NuUtilities.h"

ClassImp(NuFluctuator)

CVSID("$Id: NuFluctuator.cxx,v 1.10 2009/10/01 18:56:09 bckhouse Exp $");

using std::vector;

//____________________________________________________________________72
NuFluctuator::NuFluctuator()
  : fPDFsFilled(false),
    fFakeDataPoT(0.0),
    fMCPoT(0.0),
    fNuMuDm2(0.0),
    fNuMuSn2(0.0),
    fNuBarDm2(0.0),
    fNuBarSn2(0.0),
    fTransitionProb(0.0),
    fSelection(NuCuts::kUnknown)
{
  fhNuMuCCRecoE = 0;
  fhNuBarCCRecoE = 0;
  fMonteCarlo = 0;
  fRandom = new TRandom3(0);
  fxsecfilename = "$SRT_PRIVATE_CONTEXT/NtupleUtils/data/xsec_minos_modbyrs4_v3_5_0_mk.root";
  fSelEffFileName = "/minos/scratch/ahimmel/SystematicsStudies/Transitions/CombinedHelpers.root";
  fNuBarCCCrossSectionGraph = 0;
  fNuMuCCCrossSectionGraph = 0;
  fhNuBarCCSelectionEfficiency = 0;
  fhNuMuCCSelectionEfficiency = 0;
  fhNuBarCCCrossSections = 0;
  fhNuMuCCCrossSections = 0;
}

//____________________________________________________________________72
NuFluctuator::~NuFluctuator()
{
  if (fhNuMuCCRecoE) {
    delete fhNuMuCCRecoE; fhNuMuCCRecoE = 0;
  }
  if (fhNuBarCCRecoE) {
    delete fhNuBarCCRecoE; fhNuBarCCRecoE = 0;
  }
  if (fMonteCarlo) {
    delete fMonteCarlo; fMonteCarlo = 0;
  }
}

//____________________________________________________________________72
TH1D NuFluctuator::Fluctuate(TH1D *hin) const
{
        TH1D hData(*hin);
  hData.Reset();
  //Set all bins, including overflow (not underflow)
  for (Int_t binCounter = 1;
       binCounter < hin->GetNbinsX() + 2;
       ++binCounter){
    hData.SetBinContent
                 (binCounter,
                 fRandom->Poisson(hin->GetBinContent(binCounter)));
  }
  return hData;
}

//____________________________________________________________________72
TH1D NuFluctuator::FluctuateSampling(TH1D *hin) const
{
    TH1D hData(*hin);
    hData.Reset();
    int entries = fRandom->Poisson(hin->Integral());
    //Set all bins, including overflow (not underflow)
    for (int i = 0; i < entries; i++) {
        hData.Fill(hin->GetRandom());
    }
    return hData;
}

//____________________________________________________________________72
const Bool_t NuFluctuator::CheckConfigForPDFMaking() const
{
  if (!fhNuMuCCRecoE){
    MSG("NuFluctuator",Msg::kError)
      << "NuMuCC binning not supplied"
      << endl;
    return false;
  }
  if (!fhNuBarCCRecoE){
    MSG("NuFluctuator",Msg::kError)
      << "NuBarCC binning not supplied"
      << endl;
    return false;
  }
  if (!fMonteCarlo){
    MSG("NuFluctuator",Msg::kError)
      << "MC not supplied"
      << endl;
    return false;
  }
  if (fMCPoT<0.001){
    MSG("NuFluctuator",Msg::kError)
      << "MC PoT are worryingly small"
      << endl;
    return false;
  }
  if (fFakeDataPoT<0.001){
    MSG("NuFluctuator",Msg::kWarning)
      << "Your requested fake data PoT are worryingly small"
      << endl;
    return true;
  }

  return true;
}

//____________________________________________________________________72
void NuFluctuator::ConfigHistosForTransitionAnalysis()
{
  //Get cross-sections
  TFile *xsecfile = new TFile(fxsecfilename.c_str(),"READ");

  if (fNuMuCCCrossSectionGraph){
    delete fNuMuCCCrossSectionGraph;
    fNuMuCCCrossSectionGraph = 0;
  }
  if (fhNuMuCCCrossSections){
    delete fhNuMuCCCrossSections;
    fhNuMuCCCrossSections = 0;
  }
  TH1F* XSec_CC = (TH1F*) xsecfile->Get("h_numu_cc_tot");
  fhNuMuCCCrossSections = new TH1F(*XSec_CC);
  fhNuMuCCCrossSections->SetDirectory(0);
  Float_t *x = new Float_t[XSec_CC->GetNbinsX()];
  Float_t *y = new Float_t[XSec_CC->GetNbinsX()];
  for(int i=0;i<XSec_CC->GetNbinsX();i++) {
    x[i] = XSec_CC->GetBinCenter(i+1);
    y[i] = XSec_CC->GetBinContent(i+1);
  }
  fNuMuCCCrossSectionGraph = new TGraph(XSec_CC->GetNbinsX(),x,y);
  if (x) {delete[] x; x = 0;}
  if (y) {delete[] y; y = 0;}

  if (fNuBarCCCrossSectionGraph){
    delete fNuBarCCCrossSectionGraph;
    fNuBarCCCrossSectionGraph = 0;
  }
  if (fhNuBarCCCrossSections){
    delete fhNuBarCCCrossSections;
    fhNuBarCCCrossSections = 0;
  }
  XSec_CC = (TH1F*) xsecfile->Get("h_numubar_cc_tot");
  fhNuBarCCCrossSections = new TH1F(*XSec_CC);
  fhNuBarCCCrossSections->SetDirectory(0);
  x = new Float_t[XSec_CC->GetNbinsX()];
  y = new Float_t[XSec_CC->GetNbinsX()];
  for(int i=0;i<XSec_CC->GetNbinsX();i++) {
    x[i] = XSec_CC->GetBinCenter(i+1);
    y[i] = XSec_CC->GetBinContent(i+1);
  }
  fNuBarCCCrossSectionGraph = new TGraph(XSec_CC->GetNbinsX(),x,y);
  if (x) {delete[] x; x = 0;}
  if (y) {delete[] y; y = 0;}
  
  xsecfile->Close();
  if (xsecfile){delete xsecfile; xsecfile = 0;}
  
  //Get selection efficiencies
  TFile* selEffFile = new TFile(fSelEffFileName.c_str(),"READ");

  if (fhNuBarCCSelectionEfficiency){
    delete fhNuBarCCSelectionEfficiency;
    fhNuBarCCSelectionEfficiency = 0;
  }
  fhNuBarCCSelectionEfficiency = (TH1D*) selEffFile->Get("EfficiencyPQ_FD");
  fhNuBarCCSelectionEfficiency->SetDirectory(0);

  if (fhNuMuCCSelectionEfficiency){
    delete fhNuMuCCSelectionEfficiency;
    fhNuMuCCSelectionEfficiency = 0;
  }
  fhNuMuCCSelectionEfficiency = (TH1D*) selEffFile->Get("Efficiency_FD");
  fhNuMuCCSelectionEfficiency->SetDirectory(0);

  selEffFile->Close();
  if (selEffFile){delete selEffFile; selEffFile = 0;}
}

//____________________________________________________________________72
void NuFluctuator::MakeCPTPDFs() const
{
  if (!this->CheckConfigForPDFMaking()){return;}
        
  fhNuMuCCRecoE->Reset();
  fhNuBarCCRecoE->Reset();
        
  Int_t numMCEntries = fMonteCarlo->GetEntries();
  for (Int_t eventCounter = 0;
       eventCounter < numMCEntries;
       ++eventCounter){
    if (!(eventCounter%10000)){
      MSG("NuFluctuator",Msg::kInfo)
                        << eventCounter << " of " << numMCEntries
                        << endl
                        << "NuMu PDF has " << fhNuMuCCRecoE->GetEntries()
                        << " entries; NuBar PDF has "
                        << fhNuBarCCRecoE->GetEntries() << " entries."
                        << endl;
    }
    NuEvent event = fMonteCarlo->GetInfoObject(eventCounter);
    NuFluctuator::NuSelectionResult_t recoType = this->SelectedAs(event);
                
//     Double_t eventWeight = event.rw; 
    Double_t eventWeight = event.beamWeightRunI;
                
    //Oscillate
                eventWeight *= Oscillate(event);
    
    if (recoType==NuFluctuator::kCCNuMu){
      fhNuMuCCRecoE->Fill(event.energy,eventWeight);
    }
    if (recoType==NuFluctuator::kCCNuMuBar){
      fhNuBarCCRecoE->Fill(event.energy,eventWeight);
    }
  }
  fhNuMuCCRecoE->Scale(fFakeDataPoT/fMCPoT);
  fhNuBarCCRecoE->Scale(fFakeDataPoT/fMCPoT);
}

//____________________________________________________________________72
void NuFluctuator::MakeShiftedCPTPDFs(const TString& xmlConfigFile) const
{
  NuXMLConfig xmlConfig(xmlConfigFile);

  if (!this->CheckConfigForPDFMaking()){return;}
        
  fhNuMuCCRecoE->Reset();
  fhNuBarCCRecoE->Reset();

  NuSystematic nuSyst(xmlConfig);
        
  Int_t numMCEntries = fMonteCarlo->GetEntries();
  for (Int_t eventCounter = 0;
       eventCounter < numMCEntries;
       ++eventCounter){
    if (!(eventCounter%10000)){
      MSG("NuFluctuator",Msg::kInfo)
                        << eventCounter << " of " << numMCEntries
                        << endl
                        << "NuMu PDF has " << fhNuMuCCRecoE->GetEntries()
                        << " entries; NuBar PDF has "
                        << fhNuBarCCRecoE->GetEntries() << " entries."
                        << endl;
    }
    NuEvent event = fMonteCarlo->GetInfoObject(eventCounter);
    nuSyst.Shift(event);
    NuFluctuator::NuSelectionResult_t recoType = this->SelectedAs(event);
                
    Double_t eventWeight = event.rw;
                
    //Oscillate
                eventWeight *= Oscillate(event);
    
    if (recoType==NuFluctuator::kCCNuMu){
      fhNuMuCCRecoE->Fill(event.energy,eventWeight);
    }
    if (recoType==NuFluctuator::kCCNuMuBar){
      fhNuBarCCRecoE->Fill(event.energy,eventWeight);
    }
  }
  fhNuMuCCRecoE->Scale(fFakeDataPoT/fMCPoT);
  fhNuBarCCRecoE->Scale(fFakeDataPoT/fMCPoT);
}

//____________________________________________________________________72
Double_t NuFluctuator::Oscillate(NuEvent& event) const
{
  if (1==event.iaction){
    if (14==event.inu){
      return NuUtilities::OscillationWeight(event.neuEnMC,
                                            fNuMuDm2, fNuMuSn2);
    }
    if (-14==event.inu){
      return NuUtilities::OscillationWeight(event.neuEnMC,
                                            fNuBarDm2, fNuBarSn2);
    }
  }
  return 1; // Doesn't oscillate
}

//____________________________________________________________________72
void NuFluctuator::MakeTransitionPDFs()
{
  if (!this->CheckConfigForPDFMaking()){return;}
  this->ConfigHistosForTransitionAnalysis();

  fhNuMuCCRecoE->Reset();
  fhNuBarCCRecoE->Reset();

  Int_t numMCEntries = fMonteCarlo->GetEntries();
  for (Int_t eventCounter = 0;
       eventCounter < numMCEntries;
       ++eventCounter){
    if (!(eventCounter%10000)){
      MSG("NuFluctuator",Msg::kInfo)
        << eventCounter << " of " << numMCEntries
        << endl;
    }
    NuEvent event = fMonteCarlo->GetInfoObject(eventCounter);
    NuFluctuator::NuSelectionResult_t recoType = this->SelectedAs(event);

    Double_t eventWeight = event.rw;
    Double_t transitionWeight = event.rw;

    //Oscillate
                eventWeight *= Oscillate(event);
    
    Double_t nuBarSelEff = this->NuBarCCSelectionEfficiency(event.neuEnMC);
    Double_t nuBarXSec = this->NuBarCCCrossSection(event.neuEnMC);
    Double_t nuMuSelEff = this->NuMuCCSelectionEfficiency(event.neuEnMC);
    Double_t nuMuXSec = this->NuMuCCCrossSection(event.neuEnMC);
    
    if (recoType==NuFluctuator::kCCNuMu){
      fhNuMuCCRecoE->Fill(event.energy,eventWeight);
      if (1==event.iaction && 14==event.inu){
        transitionWeight *= fTransitionProb*(1-NuUtilities::OscillationWeight(event.neuEnMC, fNuMuDm2, fNuMuSn2));
        if (!nuMuXSec || !nuMuSelEff){
          MSG("NuFluctuator.cxx",Msg::kWarning)
            << "Trying to apply transitions to a NuMu event "
            << "that shouldn't exist."
            << endl;
        }
        else{
          transitionWeight /= nuMuXSec;
          transitionWeight /= nuMuSelEff;
          transitionWeight *= nuBarSelEff;
          transitionWeight *= nuBarXSec;
          fhNuBarCCRecoE->Fill(event.energy,transitionWeight);
        }
      }
    }
    if (recoType==NuFluctuator::kCCNuMuBar){
      fhNuBarCCRecoE->Fill(event.energy,eventWeight);
      if (1==event.iaction && -14==event.inu){
        transitionWeight *= fTransitionProb*(1-NuUtilities::OscillationWeight(event.neuEnMC, fNuMuDm2, fNuMuSn2));
        if (!nuMuXSec || !nuMuSelEff){
          MSG("NuFluctuator.cxx",Msg::kWarning)
            << "Trying to apply transitions to a NuMu event "
            << "that shouldn't exist."
            << endl;
        }
        else{
          transitionWeight *= nuMuXSec;
          transitionWeight *= nuMuSelEff;
          transitionWeight /= nuBarSelEff;
          transitionWeight /= nuBarXSec;
          fhNuMuCCRecoE->Fill(event.energy,transitionWeight);
        }
      }
    }
  }
  fhNuMuCCRecoE->Scale(fFakeDataPoT/fMCPoT);
  fhNuBarCCRecoE->Scale(fFakeDataPoT/fMCPoT);
}

//____________________________________________________________________72
void NuFluctuator::MakeShiftedTransitionPDFs(const TString& xmlConfigFile)
{
  NuXMLConfig xmlConfig(xmlConfigFile);

  if (!this->CheckConfigForPDFMaking()){return;}
  this->ConfigHistosForTransitionAnalysis();

  fhNuMuCCRecoE->Reset();
  fhNuBarCCRecoE->Reset();

  NuSystematic nuSyst(xmlConfig);

  Int_t numMCEntries = fMonteCarlo->GetEntries();
  for (Int_t eventCounter = 0;
       eventCounter < numMCEntries;
       ++eventCounter){
    if (!(eventCounter%10000)){
      MSG("NuFluctuator",Msg::kInfo)
        << eventCounter << " of " << numMCEntries
        << endl;
    }
    NuEvent event = fMonteCarlo->GetInfoObject(eventCounter);
    nuSyst.Shift(event);
    NuFluctuator::NuSelectionResult_t recoType = this->SelectedAs(event);

    Double_t eventWeight = event.rw;
    Double_t transitionWeight = event.rw;

    //Oscillate
                eventWeight *= Oscillate(event);
    
    Double_t nuBarSelEff = this->NuBarCCSelectionEfficiency(event.neuEnMC);
    Double_t nuBarXSec = this->NuBarCCCrossSection(event.neuEnMC);
    Double_t nuMuSelEff = this->NuMuCCSelectionEfficiency(event.neuEnMC);
    Double_t nuMuXSec = this->NuMuCCCrossSection(event.neuEnMC);
    
    if (recoType==NuFluctuator::kCCNuMu){
      fhNuMuCCRecoE->Fill(event.energy,eventWeight);
      if (1==event.iaction && 14==event.inu){
        transitionWeight *= fTransitionProb*(1-NuUtilities::OscillationWeight(event.neuEnMC, fNuMuDm2, fNuMuSn2));
        if (!nuMuXSec || !nuMuSelEff){
          MSG("NuFluctuator.cxx",Msg::kWarning)
            << "Trying to apply transitions to a NuMu event "
            << "that shouldn't exist."
            << endl;
        }
        else{
          transitionWeight /= nuMuXSec;
          transitionWeight /= nuMuSelEff;
          transitionWeight *= nuBarSelEff;
          transitionWeight *= nuBarXSec;
          fhNuBarCCRecoE->Fill(event.energy,transitionWeight);
        }
      }
    }
    if (recoType==NuFluctuator::kCCNuMuBar){
      fhNuBarCCRecoE->Fill(event.energy,eventWeight);
      if (1==event.iaction && -14==event.inu){
        transitionWeight *= fTransitionProb*(1-NuUtilities::OscillationWeight(event.neuEnMC, fNuMuDm2, fNuMuSn2));
        if (!nuMuXSec || !nuMuSelEff){
          MSG("NuFluctuator.cxx",Msg::kWarning)
            << "Trying to apply transitions to a NuMu event "
            << "that shouldn't exist."
            << endl;
        }
        else{
          transitionWeight *= nuMuXSec;
          transitionWeight *= nuMuSelEff;
          transitionWeight /= nuBarSelEff;
          transitionWeight /= nuBarXSec;
          fhNuMuCCRecoE->Fill(event.energy,transitionWeight);
        }
      }
    }
  }
  fhNuMuCCRecoE->Scale(fFakeDataPoT/fMCPoT);
  fhNuBarCCRecoE->Scale(fFakeDataPoT/fMCPoT);
}

//____________________________________________________________________72
const Double_t NuFluctuator
::NuBarCCCrossSection(const Double_t trueNuE) const
{
  //  return fNuBarCCCrossSectionGraph->Eval(trueNuE,0,"");
  return fhNuBarCCCrossSections->
    GetBinContent(fhNuBarCCCrossSections->FindBin(trueNuE));
}

//____________________________________________________________________72
const Double_t NuFluctuator
::NuBarCCSelectionEfficiency(const Double_t trueNuE) const
{
  return fhNuBarCCSelectionEfficiency->
    GetBinContent(fhNuBarCCSelectionEfficiency->FindBin(trueNuE));
}

//____________________________________________________________________72
const Double_t NuFluctuator
::NuMuCCCrossSection(const Double_t trueNuE) const
{
  //  return fNuMuCCCrossSectionGraph->Eval(trueNuE,0,"");
  return fhNuMuCCCrossSections->
    GetBinContent(fhNuMuCCCrossSections->FindBin(trueNuE));
}

//____________________________________________________________________72
const Double_t NuFluctuator
::NuMuCCSelectionEfficiency(const Double_t trueNuE) const
{
  return fhNuMuCCSelectionEfficiency->
    GetBinContent(fhNuMuCCSelectionEfficiency->FindBin(trueNuE));
}

//____________________________________________________________________72
const NuFluctuator::NuSelectionResult_t
NuFluctuator::SelectedAs(NuEvent& nuEvent) const
{
  NuCuts nuCuts;
  if (NuCuts::kUnknown != this->SelectionScheme()){
    nuEvent.anaVersion = this->SelectionScheme();
  }
  
  //cut on LI
  if (nuCuts.IsLI(nuEvent)) return NuFluctuator::kUnknown;
  //cut on the data quality
  if (!nuCuts.IsGoodDataQuality(nuEvent)) return NuFluctuator::kUnknown;
  //cut on the spill time
  if (!nuCuts.IsGoodTimeToNearestSpill(nuEvent)) return NuFluctuator::kUnknown;
  //cut on the beam
  if (!nuCuts.IsGoodBeam(nuEvent)) return NuFluctuator::kUnknown;

  //cut on whether event vtx is in fid vol (does nothing for CC ana)
  //    if (!nuCuts.IsInFidVol(evt,nu)) continue;
  //ensure good number of tracks in event (note preselection above)
  if (!nuCuts.IsGoodNumberOfTracks(nuEvent)) return NuFluctuator::kUnknown;
  //check trk is in fiducial volume (note preselection above)
  if (!nuCuts.IsInFidVolTrk(nuEvent)) return NuFluctuator::kUnknown;

  //require a good trk fit
  if (!nuCuts.IsGoodTrackFitPass(nuEvent)) return NuFluctuator::kUnknown;
  //require a forward going neutrino about beam direction
  if (!nuCuts.IsGoodDirCos(nuEvent)) return NuFluctuator::kUnknown;

  //cut on the PID
  if (!nuCuts.IsGoodPID(nuEvent)) {
    return NuFluctuator::kUnknown;
  }
  //cut on the fractional track momentum and sign error      
  if (!nuCuts.IsGoodSigmaQP_QP(nuEvent)) {
    return NuFluctuator::kUnknown;
  }
  //cut on the track fit probability      
  if (!nuCuts.IsGoodFitProb(nuEvent)) {
    return NuFluctuator::kUnknown;
  }
  if (!nuCuts.IsGoodMajorityCurvature(nuEvent)){
    return NuFluctuator::kUnknown;
  }
  
  if (-1==nuEvent.charge){return NuFluctuator::kCCNuMu;}
  if (1==nuEvent.charge){return NuFluctuator::kCCNuMuBar;}
  return NuFluctuator::kUnknown;
}

//____________________________________________________________________72
void NuFluctuator::NuMuCCRecoEBinning(const vector<Double_t>& numuBins)
{
  this->PDFsFilled(false);

  if (fhNuMuCCRecoE) {
    delete fhNuMuCCRecoE; fhNuMuCCRecoE = 0;
  }

  Int_t size = numuBins.size();
  if (!size){
    MSG("NuCrossSectionFitter",Msg::kError)
      << "NuMu reconstructed energy binning not supplied."
      << endl;
    return;
  }

  Double_t* recoEBinEdges = new Double_t[size];

  Int_t i=0;
  for (vector<Double_t>::const_iterator it =
         numuBins.begin();
       it != numuBins.end();
       ++it, ++i){
    recoEBinEdges[i] = *it;
  }

  fhNuMuCCRecoE = new TH1D("fhNuMuCCRecoE",
                           "",
                           size-1,
                           recoEBinEdges);
  return;
}

//____________________________________________________________________72
void NuFluctuator::NuBarCCRecoEBinning(const vector<Double_t>& nubarBins)
{
  this->PDFsFilled(false);

  if (fhNuBarCCRecoE) {
    delete fhNuBarCCRecoE; fhNuBarCCRecoE = 0;
  }

  Int_t size = nubarBins.size();
  if (!size){
    MSG("NuCrossSectionFitter",Msg::kError)
      << "NuMu reconstructed energy binning not supplied."
      << endl;
    return;
  }

  Double_t* recoEBinEdges = new Double_t[size];

  Int_t i=0;
  for (vector<Double_t>::const_iterator it =
         nubarBins.begin();
       it != nubarBins.end();
       ++it, ++i){
    recoEBinEdges[i] = *it;
  }

  fhNuBarCCRecoE = new TH1D("fhNuBarCCRecoE",
                               "",
                               size-1,
                               recoEBinEdges);
  return;
}

//____________________________________________________________________72
void NuFluctuator::MonteCarlo(const string& mcFileName)
{
  this->PDFsFilled(false);

  if (fMonteCarlo){
    delete fMonteCarlo; fMonteCarlo = 0;
  }
  fMonteCarlo = new NuTreeWrapper(mcFileName);
  fMCPoT = fMonteCarlo->GetPoT();
}

---