NuMMHelperPRL.cxx

Go to the documentation of this file.

#include <string>

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

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

ClassImp(NuMMHelperPRL)

CVSID("$Id: NuMMHelperPRL.cxx,v 1.5 2009/10/14 17:27:32 bckhouse Exp $");

//____________________________________________________________________72
NuMMHelperPRL::NuMMHelperPRL()
{
}

//____________________________________________________________________72
NuMMHelperPRL::NuMMHelperPRL(const std::string& helperFilename,
                             const std::string& xSecFilename)
{
  TFile helperFile(helperFilename.c_str(),"READ");

  fRecoVsTrueEnergy_ND = (TH2D*) helperFile.Get("RecoVsTrueEnergy_ND");
  fRecoVsTrueEnergy_ND->SetDirectory(0);
  fRecoVsTrueEnergy_FD = (TH2D*) helperFile.Get("RecoVsTrueEnergy_FD");
  fRecoVsTrueEnergy_FD->SetDirectory(0);

  fRecoVsTrueEnergy_ND_NC = (TH2D*) helperFile.Get("RecoVsTrueEnergy_ND_NC");
  if(fRecoVsTrueEnergy_ND_NC){
    fRecoVsTrueEnergy_ND_NC->SetDirectory(0);
    fRecoVsTrueEnergy_FD_NC = (TH2D*) helperFile.Get("RecoVsTrueEnergy_FD_NC");
    fRecoVsTrueEnergy_FD_NC->SetDirectory(0);

    for(int truth = 0; truth < ENCTruth::kNumTruths; ++truth){
      for(int det = 0; det <= 1; ++det){
        const TString detStr = det ? "FD" : "ND";
        fRecoVsTrueEnergy_NC_truly[truth][det] = (TH2D*)helperFile.Get("RecoVsTrueEnergy_truly"+ENCTruth::truthNames[truth]+"_"+detStr+"_NC");
        fRecoVsTrueEnergy_NC_truly[truth][det]->SetDirectory(0);
      }
    } // end for truth
  }
  else{ // We didn't get any NC plots, warn and null them out
    MSG("NuMMHelperPRL", Msg::kWarning) << "No NC plots found in helper file" << endl;
    fRecoVsTrueEnergy_FD_NC = 0;
    for(int truth = 0; truth < ENCTruth::kNumTruths; ++truth)
      for(int det = 0; det <= 1; ++det)
        fRecoVsTrueEnergy_NC_truly[truth][det] = 0;
  }
  
  fFDVsNDMatrixRW = (TH2D*) helperFile.Get("FDVsNDMatrixRW");
  if(fFDVsNDMatrixRW) fFDVsNDMatrixRW->SetDirectory(0);  
  fFDVsNDMatrixXSecRW = (TH2D*) helperFile.Get("FDVsNDMatrixXSecRW");
  if(fFDVsNDMatrixXSecRW) fFDVsNDMatrixXSecRW->SetDirectory(0);
  
  fEfficiency_ND = (TH1D*) helperFile.Get("Efficiency_ND");
  fEfficiency_ND->SetDirectory(0);
  fEfficiency_FD = (TH1D*) helperFile.Get("Efficiency_FD");
  fEfficiency_FD->SetDirectory(0);

  //The purity histograms in the helper file remove both the ND
  //background and the wrong sign CCs: the PRL analysis doesn't
  //remove the latter.
  fPurity_FD = (TH1D*) helperFile.Get("NCContamination_FD");
  fPurity_FD->SetDirectory(0);
  fPurity_ND = (TH1D*) helperFile.Get("NCContamination_ND");
  fPurity_ND->SetDirectory(0);
  for (int i=1; i<=fPurity_ND->GetNbinsX(); ++i){
    fPurity_ND->SetBinContent(i, 1.0 -
                              fPurity_ND->GetBinContent(i));
    fPurity_FD->SetBinContent(i, 1.0 -
                              fPurity_FD->GetBinContent(i));
  }

  //Tau helpers
  fEfficiencyTau_FD = (TH1D*) helperFile.Get("EfficiencyTau_FD");
  fEfficiencyTau_FD->SetDirectory(0);
  fRecoVsTrueEnergyTau_FD = (TH2D*)
    helperFile.Get("RecoVsTrueEnergyTau_FD");
  fRecoVsTrueEnergyTau_FD->SetDirectory(0);
  
  helperFile.Close();
  
  //Get cross-sections (numu)
  TFile *xsecfile = new TFile(xSecFilename.c_str(),"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 (tau)
  xsecfile = new TFile(xSecFilename.c_str(),"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;}
}

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

---