NuMMHelperCPT.cxx

Go to the documentation of this file.

#include <string>

#include "TFile.h"
#include "TGraph.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"

#include "MessageService/MsgService.h"
#include "NtupleUtils/NuMMHelperCPT.h"

ClassImp(NuMMHelperCPT)

CVSID("$Id: NuMMHelperCPT.cxx,v 1.8 2009/10/15 16:40:55 ahimmel Exp $");

//____________________________________________________________________72
NuMMHelperCPT::NuMMHelperCPT() : doTaus(false)
{
}

//____________________________________________________________________72
NuMMHelperCPT::NuMMHelperCPT(const char * helperFilename,
                             const char * xSecFilename)
{
  Construct(helperFilename, xSecFilename);
}

//____________________________________________________________________72
NuMMHelperCPT::NuMMHelperCPT(const std::string helperFilename,
                             const std::string xSecFilename)
{
  Construct(helperFilename.c_str(), xSecFilename.c_str());
}

//____________________________________________________________________72
void NuMMHelperCPT::Construct(const char * helperFilename,
                              const char * xSecFilename)
{
  TFile helperFile(helperFilename,"READ");

  fRecoVsTrueEnergy_ND = (TH2D*) helperFile.Get("RecoVsTrueEnergy_ND");
  fRecoVsTrueEnergy_ND->SetDirectory(0);
  fRecoVsTrueEnergy_FD = (TH2D*) helperFile.Get("RecoVsTrueEnergy_FD");
  fRecoVsTrueEnergy_FD->SetDirectory(0);
  
  fFDVsNDMatrixRW = (TH2D*) helperFile.Get("FDVsNDMatrixRW");
  if(fFDVsNDMatrixRW) fFDVsNDMatrixRW->SetDirectory(0);
  
  fEfficiency_ND = (TH1D*) helperFile.Get("Efficiency_ND");
  fEfficiency_ND->SetDirectory(0);
  fEfficiency_FD = (TH1D*) helperFile.Get("Efficiency_FD");
  fEfficiency_FD->SetDirectory(0);
  //Purity histograms remove both the NC and wrong-sign CC backgrounds
  fPurity_ND = (TH1D*) helperFile.Get("Purity_ND");
  fPurity_ND->SetDirectory(0);
  fPurity_FD = (TH1D*) helperFile.Get("Purity_FD");
  fPurity_FD->SetDirectory(0);
  
  fNCContamination_FD = (TH1D*) helperFile.Get("NCContamination_FD");
  fNCContamination_FD->SetDirectory(0);
  fRecoVsTrueCCContamination_FD = (TH2D*)
    helperFile.Get("CCContaminationRecoVsTrue_FD");
  fRecoVsTrueCCContamination_FD->SetDirectory(0);
  fOtherEfficiency_FD = (TH1D*)
    helperFile.Get("OtherEfficiency_FD");
  fOtherEfficiency_FD->SetDirectory(0);

  //Antineutrinos:
  fRecoVsTrueEnergyPQ_ND = (TH2D*) helperFile.Get("RecoVsTrueEnergyPQ_ND");
  fRecoVsTrueEnergyPQ_ND->SetDirectory(0);
  fRecoVsTrueEnergyPQ_FD = (TH2D*) helperFile.Get("RecoVsTrueEnergyPQ_FD");
  fRecoVsTrueEnergyPQ_FD->SetDirectory(0);
  
  fFDVsNDMatrixRWBar = (TH2D*) helperFile.Get("FDVsNDMatrixRWPQ");
  if(fFDVsNDMatrixRWBar) fFDVsNDMatrixRWBar->SetDirectory(0);
  
  fEfficiencyPQ_ND = (TH1D*) helperFile.Get("EfficiencyPQ_ND");
  fEfficiencyPQ_ND->SetDirectory(0);
  fEfficiencyPQ_FD = (TH1D*) helperFile.Get("EfficiencyPQ_FD");
  fEfficiencyPQ_FD->SetDirectory(0);
  //Purity histograms remove both the NC and wrong-sign CC backgrounds
  fPurityPQ_ND = (TH1D*) helperFile.Get("PurityPQ_ND");
  fPurityPQ_ND->SetDirectory(0);
  fPurityPQ_FD = (TH1D*) helperFile.Get("PurityPQ_FD");
  fPurityPQ_FD->SetDirectory(0);
  
  fNCContaminationPQ_FD = (TH1D*) helperFile.Get("NCContaminationPQ_FD");
  fNCContaminationPQ_FD->SetDirectory(0);
  fRecoVsTrueCCContaminationPQ_FD = (TH2D*)
    helperFile.Get("CCContaminationRecoVsTruePQ_FD");
  fRecoVsTrueCCContaminationPQ_FD->SetDirectory(0);
  fOtherEfficiencyPQ_FD = (TH1D*)
    helperFile.Get("OtherEfficiencyPQ_FD");
  fOtherEfficiencyPQ_FD->SetDirectory(0);

  //Tau helpers
  fEfficiencyTau_FD = (TH1D*) helperFile.Get("EfficiencyTau_FD");
  if (!fEfficiencyTau_FD)  doTaus = false;
  else                     doTaus = true;
  if (doTaus) {
    fEfficiencyTau_FD->SetDirectory(0);
    fRecoVsTrueEnergyTau_FD = (TH2D*)
    helperFile.Get("RecoVsTrueEnergyTau_FD");
    fRecoVsTrueEnergyTau_FD->SetDirectory(0);
    
    //TauBar helpers
    fEfficiencyTauPQ_FD = (TH1D*) helperFile.Get("EfficiencyTauPQ_FD");
    fEfficiencyTauPQ_FD->SetDirectory(0);
    fRecoVsTrueEnergyTauPQ_FD = (TH2D*)
    helperFile.Get("RecoVsTrueEnergyTauPQ_FD");
    fRecoVsTrueEnergyTauPQ_FD->SetDirectory(0);
  }
  
  helperFile.Close();
  
  //Get cross-sections (numu)
  TFile *xsecfile = new TFile(xSecFilename,"READ");
  TH1F* XSec_CC = (TH1F*) xsecfile->Get("h_numu_cc_tot");
  XSec_CC->SetDirectory(0);
  xsecfile->Close();
  if (xsecfile){delete xsecfile; xsecfile = 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);
  }
  fXSec_CC_Graph = new TGraph(XSec_CC->GetNbinsX(),x,y);
  if (x) {delete[] x; x = 0;}
  if (y) {delete[] y; y = 0;}
  
  //Get cross-sections (numubar)
  xsecfile = new TFile(xSecFilename,"READ");
  TH1F* XSec_CCBar = (TH1F*) xsecfile->Get("h_numubar_cc_tot");
  XSec_CCBar->SetDirectory(0);
  xsecfile->Close();
  if (xsecfile){delete xsecfile; xsecfile = 0;}
  x = new Float_t[XSec_CCBar->GetNbinsX()];
  y = new Float_t[XSec_CCBar->GetNbinsX()];
  for(int i=0;i<XSec_CCBar->GetNbinsX();i++) {
    x[i] = XSec_CCBar->GetBinCenter(i+1);
    y[i] = XSec_CCBar->GetBinContent(i+1);
  }
  fXSec_CCBar_Graph = new TGraph(XSec_CCBar->GetNbinsX(),x,y);
  if (x) {delete[] x; x = 0;}
  if (y) {delete[] y; y = 0;}
  
  if (doTaus) {
    //Get cross-sections (tau)
    xsecfile = new TFile(xSecFilename,"READ");
    TH1F* tau_XSec_CC = (TH1F*) xsecfile->Get("h_nutau_cc_tot");
    tau_XSec_CC->SetDirectory(0);
    xsecfile->Close();
    if (xsecfile){delete xsecfile; xsecfile = 0;}
    x = new Float_t[tau_XSec_CC->GetNbinsX()];
    y = new Float_t[tau_XSec_CC->GetNbinsX()];
    for(int i=0;i<tau_XSec_CC->GetNbinsX();i++) {
      x[i] = tau_XSec_CC->GetBinCenter(i+1);
      y[i] = tau_XSec_CC->GetBinContent(i+1);
    }
    fTau_XSec_CC_Graph = new TGraph(tau_XSec_CC->GetNbinsX(),x,y);
    if (x) {delete[] x; x = 0;}
    if (y) {delete[] y; y = 0;}
    
    //Get cross-sections (taubar)
    xsecfile = new TFile(xSecFilename,"READ");
    TH1F* taubar_XSec_CC = (TH1F*) xsecfile->Get("h_nutaubar_cc_tot");
    taubar_XSec_CC->SetDirectory(0);
    xsecfile->Close();
    if (xsecfile){delete xsecfile; xsecfile = 0;}
    x = new Float_t[taubar_XSec_CC->GetNbinsX()];
    y = new Float_t[taubar_XSec_CC->GetNbinsX()];
    for(int i=0;i<taubar_XSec_CC->GetNbinsX();i++) {
      x[i] = taubar_XSec_CC->GetBinCenter(i+1);
      y[i] = taubar_XSec_CC->GetBinContent(i+1);
    }
    fTauBar_XSec_CC_Graph = new TGraph(taubar_XSec_CC->GetNbinsX(),x,y);
    if (x) {delete[] x; x = 0;}
    if (y) {delete[] y; y = 0;}
  }
  
  if (fRecoVsTrueCCContamination_FD){ 
    //Normalise the FD CC contamination true to reco CC matrix
    //First loop over reco axis:
    for (Int_t i=1; i<=fRecoVsTrueCCContamination_FD->GetNbinsX()+1; ++i){
      //Find out the total number of true events in this column:
      Double_t trueEvents = 0.0;
      for (Int_t j=1; j<=fRecoVsTrueCCContamination_FD->GetNbinsY()+1; ++j){
        trueEvents += fRecoVsTrueCCContamination_FD->GetBinContent(i,j);
      }
      //Loop again, this time normalise.
      for (Int_t j=1; j<=fRecoVsTrueCCContamination_FD->GetNbinsY()+1; ++j){
        if (trueEvents>0){
          Double_t oldContent =
            fRecoVsTrueCCContamination_FD->GetBinContent(i,j);
          Double_t oldError =
            fRecoVsTrueCCContamination_FD->GetBinError(i,j);
          fRecoVsTrueCCContamination_FD->SetBinContent
            (i,j,oldContent/trueEvents);
          fRecoVsTrueCCContamination_FD->SetBinError
            (i,j,oldError/trueEvents);
        }
        else{
          fRecoVsTrueCCContamination_FD->SetBinContent(i,j,0);
          
          fRecoVsTrueCCContamination_FD->SetBinError(i,j,0);
        }
      }
    }
  }
  
  if (fRecoVsTrueCCContaminationPQ_FD){ 
    //Normalise the FD CC contamination true to reco CC matrix
    //First loop over reco axis:
    for (Int_t i=1; i<=fRecoVsTrueCCContaminationPQ_FD->GetNbinsX()+1; ++i){
      //Find out the total number of true events in this column:
      Double_t trueEvents = 0.0;
      for (Int_t j=1; j<=fRecoVsTrueCCContaminationPQ_FD->GetNbinsY()+1; ++j){
        trueEvents += fRecoVsTrueCCContaminationPQ_FD->GetBinContent(i,j);
      }
      //Loop again, this time normalise.
      for (Int_t j=1; j<=fRecoVsTrueCCContaminationPQ_FD->GetNbinsY()+1; ++j){
        if (trueEvents>0){
          Double_t oldContent =
            fRecoVsTrueCCContaminationPQ_FD->GetBinContent(i,j);
          Double_t oldError =
            fRecoVsTrueCCContaminationPQ_FD->GetBinError(i,j);
          fRecoVsTrueCCContaminationPQ_FD->SetBinContent
            (i,j,oldContent/trueEvents);
          fRecoVsTrueCCContaminationPQ_FD->SetBinError
            (i,j,oldError/trueEvents);
        }
        else{
          fRecoVsTrueCCContaminationPQ_FD->SetBinContent(i,j,0);
          
          fRecoVsTrueCCContaminationPQ_FD->SetBinError(i,j,0);
        }
      }
    }
  }
  
  TF1 *xsec_func = new TF1("xsec_func","[5]*(-[1]+([0]-1.+2.*(1.-[0])*x/25+([0]-1.)*x*x/(25*25))*(x<25)+([4]*([2]-x)^3+[3]*([2]-x)^2)*(x<[2]))",0,50);
  xsec_func->SetParameter(0, 0.895);
  xsec_func->SetParameter(1, 0.0617);
  
  xsec_func->SetParameter(2, 4.0);
  xsec_func->SetParameter(3, 7.59e-3);
  xsec_func->SetParameter(4, -8.05e-4);
  
  xsec_func->SetParameter(5, 0);
  fXSec_CCBar_Graph_orig = 0;
}

//____________________________________________________________________72
void NuMMHelperCPT::OverrideNuBeamMatrix(const TH2D* beamMatrix)
{
  if (fFDVsNDMatrixRW){delete fFDVsNDMatrixRW; fFDVsNDMatrixRW=0;}
  fFDVsNDMatrixRW = new TH2D(*beamMatrix);
}

//____________________________________________________________________72
void NuMMHelperCPT::OverrideBarBeamMatrix(const TH2D* beamMatrix)
{
  if (fFDVsNDMatrixRWBar)
    {delete fFDVsNDMatrixRWBar; fFDVsNDMatrixRWBar=0;}
  fFDVsNDMatrixRWBar = new TH2D(*beamMatrix);
}


//____________________________________________________________________

void NuMMHelperCPT::ShiftBarXSecGraph(double scale)
{
  if (fXSec_CCBar_Graph_orig) {
    MAXMSG("NuMMHelperCPT",Msg::kWarning, 1) << "You did not reset the Bar xsec after modifying it.  Doing so now." << endl;
    ResetBarXSecGraph();
  }
  
  MAXMSG("NuMMHelperCPT",Msg::kInfo, 1) << "Making a Bar xsec backup copy." << endl;
  
  // Backup the original object, create a new one for modifying
  fXSec_CCBar_Graph_orig = fXSec_CCBar_Graph;
  fXSec_CCBar_Graph = new TGraph(*fXSec_CCBar_Graph);
  
  // Set the scale of the xsec shift
  xsec_func->SetParameter(5, scale);
  
  Double_t x, y;
  for (int i = 0; i < fXSec_CCBar_Graph->GetN(); ++i) {
    fXSec_CCBar_Graph->GetPoint(i, x, y);
    y *= xsec_func->Eval(x);
    fXSec_CCBar_Graph->SetPoint(i, x, y);    
  }
}

//____________________________________________________________________

void NuMMHelperCPT::ResetBarXSecGraph()
{
  if (fXSec_CCBar_Graph_orig) {
    MAXMSG("NuMMHelperCPT",Msg::kInfo, 1) << "Reverting Bar to backup copy." << endl;
    
    delete fXSec_CCBar_Graph;
    fXSec_CCBar_Graph = fXSec_CCBar_Graph_orig;
    fXSec_CCBar_Graph_orig = 0;
  }
  else {
    MAXMSG("NuMMHelperCPT",Msg::kWarning, 1) << "Resetting unmodified Bars data.  Not harmful, just superfluous." << endl;
  }
}

---