MadQEID.cxx

Go to the documentation of this file.

//
// MadQEID
//
// Mark's CC-QE analysis.
//
// Created:  M. Dorman -- July, 2005
//
// $Author: rhatcher $ 
//
// $Revision: 1.13 $
// 
// $Name:  $
//
// $Id: MadQEID.cxx,v 1.13 2006/05/22 19:16:38 rhatcher Exp $
//


#ifndef madqeid_cxx
#define madqeid_cxx

#include "Mad/MadQEID.h"
#include "Mad/NearbyEvents.h"
#include "Mad/HoughTrans.h"
#include "Mad/MadBase.h"
#include "CandNtupleSR/NtpSREvent.h"
#include <cmath>
#include "TFile.h"
#include "TDirectory.h"
#include "TH1.h"
#include "TH2.h"

//**********************************************************
MadQEID::MadQEID() 
  : ntrks(-1), nshows(-1), houghPeak(-1), recoWsqd(1000.0),
    rePHfrac(-1.0), PCAvar(-1.0), pdfs_normalised(kFALSE)
{

  for(Int_t i=0;i<18;i++){
    rePHfracCCQE[i] = 0;
    rePHfracCCR[i] = 0;
    rePHfracCCDIS[i] = 0;
    rePHfracNC[i] = 0;
    rePHfracBG[i] = 0;
    recoW2CCQE[i] = 0;
    recoW2CCR[i] = 0;
    recoW2CCDIS[i] = 0;
    recoW2NC[i] = 0;
    recoW2BG[i] = 0;
    ntrkCCQE[i] = 0;
    ntrkCCR[i] = 0;
    ntrkCCDIS[i] = 0;
    ntrkNC[i] = 0;
    ntrkBG[i] = 0;
    nshowCCQE[i] = 0;
    nshowCCR[i] = 0;    
    nshowCCDIS[i] = 0;
    nshowNC[i] = 0;
    nshowBG[i] = 0;
    hpeakCCQE[i] = 0;
    hpeakCCR[i] = 0;
    hpeakCCDIS[i] = 0;
    hpeakNC[i] = 0;
    hpeakBG[i] = 0;
    PCACCQE[i] = 0;
    PCACCR[i] = 0;
    PCACCDIS[i] = 0;
    PCANC[i] = 0;
    PCABG[i] = 0;
    trainingHists[i] = 0;
  }

  cutsHist = 0;

}

//**********************************************************
void MadQEID::InitPDFs()
{
  TH1::AddDirectory(kFALSE); // do not assign histograms to TDirectory

  for(Int_t i=0;i<18;i++){
    char name[30];
    sprintf(name,"rePHfracCCQE[%d]",i);
    rePHfracCCQE[i] = new TH1F(name,name,25,0.0,1.0);
    sprintf(name,"rePHfracCCR[%d]",i);
    rePHfracCCR[i] = new TH1F(name,name,25,0.0,1.0);
    sprintf(name,"rePHfracCCDIS[%d]",i);
    rePHfracCCDIS[i] = new TH1F(name,name,25,0.0,1.0);
    sprintf(name,"rePHfracNC[%d]",i);
    rePHfracNC[i] = new TH1F(name,name,25,0.0,1.0);
    sprintf(name,"rePHfracBG[%d]",i);
    rePHfracBG[i] = new TH1F(name,name,25,0.0,1.0);

    sprintf(name,"recoW2CCQE[%d]",i);
    recoW2CCQE[i] = new TH1F(name,name,75,-5,10);
    sprintf(name,"recoW2CCR[%d]",i);
    recoW2CCR[i] = new TH1F(name,name,75,-5,10);
    sprintf(name,"recoW2CCDIS[%d]",i);
    recoW2CCDIS[i] = new TH1F(name,name,75,-5,10);
    sprintf(name,"recoW2NC[%d]",i);
    recoW2NC[i] = new TH1F(name,name,75,-5,10);
    sprintf(name,"recoW2BG[%d]",i);
    recoW2BG[i] = new TH1F(name,name,75,-5,10);

    sprintf(name,"ntrkCCQE[%d]",i);
    ntrkCCQE[i] = new TH1F(name,name,5,0,5);
    sprintf(name,"ntrkCCR[%d]",i);
    ntrkCCR[i] = new TH1F(name,name,5,0,5);
    sprintf(name,"ntrkCCDIS[%d]",i);
    ntrkCCDIS[i] = new TH1F(name,name,5,0,5);
    sprintf(name,"ntrkNC[%d]",i);
    ntrkNC[i] = new TH1F(name,name,5,0,5);
    sprintf(name,"ntrkBG[%d]",i);
    ntrkBG[i] = new TH1F(name,name,5,0,5);

    sprintf(name,"nshowCCQE[%d]",i);
    nshowCCQE[i] = new TH1F(name,name,10,0,10);
    sprintf(name,"nshowCCR[%d]",i);
    nshowCCR[i] = new TH1F(name,name,10,0,10);
    sprintf(name,"nshowCCDIS[%d]",i);
    nshowCCDIS[i] = new TH1F(name,name,10,0,10);
    sprintf(name,"nshowNC[%d]",i);
    nshowNC[i] = new TH1F(name,name,10,0,10);
    sprintf(name,"nshowBG[%d]",i);
    nshowBG[i] = new TH1F(name,name,10,0,10);

    sprintf(name,"hpeakCCQE[%d]",i);
    hpeakCCQE[i] = new TH1F(name,name,25,0,50);
    sprintf(name,"hpeakCCR[%d]",i);
    hpeakCCR[i] = new TH1F(name,name,25,0,50);
    sprintf(name,"hpeakCCDIS[%d]",i);
    hpeakCCDIS[i] = new TH1F(name,name,25,0,50);
    sprintf(name,"hpeakNC[%d]",i);
    hpeakNC[i] = new TH1F(name,name,25,0,50);
    sprintf(name,"hpeakBG[%d]",i);
    hpeakBG[i] = new TH1F(name,name,25,0,50);

    sprintf(name,"PCACCQE[%d]",i);
    PCACCQE[i] = new TH1F(name,name,25,0,100000);
    sprintf(name,"PCACCR[%d]",i);
    PCACCR[i] = new TH1F(name,name,25,0,100000);
    sprintf(name,"PCACCDIS[%d]",i);
    PCACCDIS[i] = new TH1F(name,name,25,0,100000);
    sprintf(name,"PCANC[%d]",i);
    PCANC[i] = new TH1F(name,name,25,0,100000);
    sprintf(name,"PCABG[%d]",i);
    PCABG[i] = new TH1F(name,name,25,0,100000);

  }  

  TH1::AddDirectory(kTRUE); // revert above
}

//**********************************************************
void MadQEID::InitTrainingHists()
{
  for(Int_t i=0;i<18;i++){
    trainingHists[i] = 0;
  }
  TH1::AddDirectory(kFALSE); // do not assign histograms to TDirectory
  for(Int_t i=0;i<18;i++){
    char name[30];
    sprintf(name,"trainingHists[%d]",i);
    trainingHists[i] = new TH1F(name,name,1200,-6.0,6.0);
  }
  Float_t axarray[19] = {0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,
                         5.0,6.0,7.0,8.0,10.0,12.0,14.0,17.0,20.0};
  cutsHist = new TH1F("cuthist","cuthist",18,axarray);
  TH1::AddDirectory(kTRUE); // revert above
}



//**********************************************************
void MadQEID::ResetPDFs()
{
  for(Int_t i=0;i<18;i++){
    rePHfracCCQE[i]->Reset();
    rePHfracCCR[i]->Reset();
    rePHfracCCDIS[i]->Reset();
    rePHfracNC[i]->Reset();
    rePHfracBG[i]->Reset();
    recoW2CCQE[i]->Reset();
    recoW2CCR[i]->Reset();
    recoW2CCDIS[i]->Reset();
    recoW2NC[i]->Reset();
    recoW2BG[i]->Reset();
    ntrkCCQE[i]->Reset();
    ntrkCCR[i]->Reset();
    ntrkCCDIS[i]->Reset();
    ntrkNC[i]->Reset();
    ntrkBG[i]->Reset();
    nshowCCQE[i]->Reset();
    nshowCCR[i]->Reset();
    nshowCCDIS[i]->Reset();
    nshowNC[i]->Reset();
    nshowBG[i]->Reset();
    hpeakCCQE[i]->Reset();
    hpeakCCR[i]->Reset();
    hpeakCCDIS[i]->Reset();
    hpeakNC[i]->Reset();
    hpeakBG[i]->Reset();
    PCACCQE[i]->Reset();
    PCACCR[i]->Reset();
    PCACCDIS[i]->Reset();
    PCANC[i]->Reset();
    PCABG[i]->Reset();
  }
}
        
//**********************************************************
Bool_t MadQEID::CalcQEVars(MadBase* mb, Int_t evtindex)
{
  ntrks = -1;
  nshows = -1;
  houghPeak = -1;
  rePHfrac = -1.0;
  PCAvar = -1.0;

  if(!mb) return false;
  const NtpSREvent* evt = mb->GetEvent(evtindex);
  if(!evt) return false;
  ntrks = evt->ntrack;
  nshows = evt->nshower;
  Int_t track_index = -1;
  const NtpSRTrack* trk = mb->GetLargestTrackFromEvent(evtindex,track_index);
  if(track_index==-1) return false;
  Float_t EvtPH=0.0;
  Float_t PHremaining=0.0;
  Int_t numbighits=0;

  // MAK 02/08/05: no reason to create HoughTrans objects on the heap with new
  // make these objects static so we don't have to keep deleting them
  static HoughTrans HTUpln(80,-4.0,4.0);
  static HoughTrans HTVpln(80,-4.0,4.0);
  HTUpln.ResetHough();
  HTVpln.ResetHough();

  HTUpln.SetVtxz(evt->vtx.z);
  HTVpln.SetVtxz(evt->vtx.z);
  for(Int_t istp=0;istp<evt->nstrip;istp++){
    int stripindex = evt->stp[istp];
    const NtpSRStrip* stp = mb->GetStrip(stripindex);
    if(!stp) continue;
    Float_t zPos = stp->z;
    Float_t tPos = stp->tpos;
    Float_t pes = (stp->ph0.pe)+(stp->ph1.pe);
    Float_t scor = (stp->ph0.sigcor)+(stp->ph1.sigcor);
    Float_t vtxz = evt->vtx.z;
    Int_t pln = stp->plane;
    EvtPH+=scor;
    if(pes<1.5) continue;
    if((zPos-vtxz)>2) continue;
    if(!IsTrkHit(zPos,tPos,trk,mb)){
      PHremaining+=scor;
      if(pes>20) numbighits++;
      if(IsVPln(pln)) HTVpln.FillHough(zPos,tPos);
      if(!IsVPln(pln)) HTUpln.FillHough(zPos,tPos);
    }
    for(Int_t ishw=0;ishw<nshows;ishw++){
      int shwindex = evt->shw[ishw];    
      const NtpSRShower* shwr = mb->GetShower(shwindex);
      if(!shwr) continue;
      if(IsShowHit(zPos,tPos,shwr,mb) && IsTrkHit(zPos,tPos,trk,mb)){
        PHremaining+=(scor-545);  //how much is a MIP in sigcor??? ~545(ND)
        if(pes>26) numbighits++;  //and in pes??? ~6 (gain ~100 in ND)
      }
    }
  }

  houghPeak = (HTUpln.GetPeakHeight())+(HTVpln.GetPeakHeight());

  if(EvtPH!=0) rePHfrac = (PHremaining/EvtPH);
  if(numbighits<=0 || PHremaining<=0) PCAvar = 0.0;
  else PCAvar = (((2500*numbighits)/PHremaining)*(2500*numbighits+PHremaining));

  return true;
}

//**********************************************************
Float_t MadQEID::CalcQEID(MadBase* mb, Int_t evtindex, Float_t reco_enu, Float_t recoW2)
{
  if(pdfs_normalised==kFALSE) NormalisePDFs();

  if(reco_enu>=20.0) return -1000.0;
  if(!CalcQEVars(mb,evtindex)) return -1000.0;

  recoWsqd = recoW2;

  if(recoW2>=10.0) return -5.0;
  if(PCAvar>=50000.0) return -5.0;

  Int_t currentBin = GetEnergyBin(reco_enu);

  //QE
  Float_t QErePHfrac=rePHfracCCQE[currentBin]->GetBinContent(rePHfracCCQE[currentBin]->FindBin(rePHfrac));
  Float_t QErecoW2=recoW2CCQE[currentBin]->GetBinContent(recoW2CCQE[currentBin]->FindBin(recoWsqd));
  Float_t QEnshow=nshowCCQE[currentBin]->GetBinContent(nshowCCQE[currentBin]->FindBin(nshows));
  //Float_t QEntrk=ntrkCCQE[currentBin]->GetBinContent(ntrkCCQE[currentBin]->FindBin(ntrks));
  Float_t QEhpk=hpeakCCQE[currentBin]->GetBinContent(hpeakCCQE[currentBin]->FindBin(houghPeak));
  Float_t QEpca=PCACCQE[currentBin]->GetBinContent(PCACCQE[currentBin]->FindBin(PCAvar));
  
  //Res
  Float_t RrePHfrac=rePHfracCCR[currentBin]->GetBinContent(rePHfracCCR[currentBin]->FindBin(rePHfrac));
  Float_t RrecoW2=recoW2CCR[currentBin]->GetBinContent(recoW2CCR[currentBin]->FindBin(recoWsqd));
  Float_t Rnshow=nshowCCR[currentBin]->GetBinContent(nshowCCR[currentBin]->FindBin(nshows));
  //Float_t Rntrk=ntrkCCR[currentBin]->GetBinContent(ntrkCCR[currentBin]->FindBin(ntrks));
  Float_t Rhpk=hpeakCCR[currentBin]->GetBinContent(hpeakCCR[currentBin]->FindBin(houghPeak));
  Float_t Rpca=PCACCR[currentBin]->GetBinContent(PCACCR[currentBin]->FindBin(PCAvar));
  
  //DIS
  Float_t DISrePHfrac=rePHfracCCDIS[currentBin]->GetBinContent(rePHfracCCDIS[currentBin]->FindBin(rePHfrac));
  Float_t DISrecoW2=recoW2CCDIS[currentBin]->GetBinContent(recoW2CCDIS[currentBin]->FindBin(recoWsqd));
  Float_t DISnshow=nshowCCDIS[currentBin]->GetBinContent(nshowCCDIS[currentBin]->FindBin(nshows));
  //Float_t DISntrk=ntrkCCDIS[currentBin]->GetBinContent(ntrkCCDIS[currentBin]->FindBin(ntrks));
  Float_t DIShpk=hpeakCCDIS[currentBin]->GetBinContent(hpeakCCDIS[currentBin]->FindBin(houghPeak));
  Float_t DISpca=PCACCDIS[currentBin]->GetBinContent(PCACCDIS[currentBin]->FindBin(PCAvar));
  
  //NC
  Float_t NCrePHfrac=rePHfracNC[currentBin]->GetBinContent(rePHfracNC[currentBin]->FindBin(rePHfrac));
  Float_t NCrecoW2=recoW2NC[currentBin]->GetBinContent(recoW2NC[currentBin]->FindBin(recoWsqd));
  Float_t NCnshow=nshowNC[currentBin]->GetBinContent(nshowNC[currentBin]->FindBin(nshows));
  //Float_t NCntrk=ntrkNC[currentBin]->GetBinContent(ntrkNC[currentBin]->FindBin(ntrks));
  Float_t NChpk=hpeakNC[currentBin]->GetBinContent(hpeakNC[currentBin]->FindBin(houghPeak));
  Float_t NCpca=PCANC[currentBin]->GetBinContent(PCANC[currentBin]->FindBin(PCAvar));

  Float_t likliQE = QErePHfrac*QErecoW2*QEnshow*QEhpk*QEpca;
  Float_t liklinotQE = (RrePHfrac*RrecoW2*Rnshow*Rhpk*Rpca)
    +(DISrePHfrac*DISrecoW2*DISnshow*DIShpk*DISpca)
    +(NCrePHfrac*NCrecoW2*NCnshow*NChpk*NCpca);

  if(likliQE==0 && liklinotQE==0) return -1000.0;
  if(likliQE==0 && liklinotQE>0) return -5.0;
  if(likliQE>0 && liklinotQE==0) return 5.0;
  Float_t PID = -(sqrt(-log(likliQE)))+(sqrt(-log(liklinotQE)));

  return PID;  

}

//**********************************************************
Float_t MadQEID::AltCalcQEID(MadBase* mb, Int_t evtindex, Float_t reco_enu, Float_t recoW2)
{
  if(pdfs_normalised==kFALSE) NormalisePDFs();

  if(reco_enu>=20.0) return -1000.0;
  if(!CalcQEVars(mb,evtindex)) return -1000.0;

  recoWsqd = recoW2;

  if(recoW2>=10.0) return -5.0;
  if(PCAvar>=50000.0) return -5.0;

  Int_t currentBin = GetEnergyBin(reco_enu);

  //QE
  Float_t QErePHfrac=rePHfracCCQE[currentBin]->GetBinContent(rePHfracCCQE[currentBin]->FindBin(rePHfrac));
  Float_t QErecoW2=recoW2CCQE[currentBin]->GetBinContent(recoW2CCQE[currentBin]->FindBin(recoWsqd));
  Float_t QEnshow=nshowCCQE[currentBin]->GetBinContent(nshowCCQE[currentBin]->FindBin(nshows));
  //Float_t QEntrk=ntrkCCQE[currentBin]->GetBinContent(ntrkCCQE[currentBin]->FindBin(ntrks));
  Float_t QEhpk=hpeakCCQE[currentBin]->GetBinContent(hpeakCCQE[currentBin]->FindBin(houghPeak));
  Float_t QEpca=PCACCQE[currentBin]->GetBinContent(PCACCQE[currentBin]->FindBin(PCAvar));
  
  //all backgrounds together 
  Float_t BGrePHfrac=rePHfracBG[currentBin]->GetBinContent(rePHfracBG[currentBin]->FindBin(rePHfrac));
  Float_t BGrecoW2=recoW2BG[currentBin]->GetBinContent(recoW2BG[currentBin]->FindBin(recoWsqd));
  Float_t BGnshow=nshowBG[currentBin]->GetBinContent(nshowBG[currentBin]->FindBin(nshows));
  //Float_t BGntrk=ntrkBG[currentBin]->GetBinContent(ntrkBG[currentBin]->FindBin(ntrks));
  Float_t BGhpk=hpeakBG[currentBin]->GetBinContent(hpeakBG[currentBin]->FindBin(houghPeak));
  Float_t BGpca=PCABG[currentBin]->GetBinContent(PCABG[currentBin]->FindBin(PCAvar));
  
  Float_t likliQE = QErePHfrac*QErecoW2*QEnshow*QEhpk*QEpca;
  Float_t liklinotQE = BGrePHfrac*BGrecoW2*BGnshow*BGhpk*BGpca;

  if(likliQE==0 && liklinotQE==0) return -1000.0;
  if(likliQE==0 && liklinotQE>0) return -5.0;
  if(likliQE>0 && liklinotQE==0) return 5.0;
  Float_t PID = -(sqrt(-log(likliQE)))+(sqrt(-log(liklinotQE)));

  return PID;  

}

//**********************************************************
Int_t MadQEID::PassMEDQECut(Float_t reco_enu, Float_t QEIDval)
{
  if(reco_enu>=20.0) return -1;
  if(QEIDval==-1000.0) return -1;
  if(QEIDval==-5.0) return 0;
  Int_t currentBin = GetEnergyBin(reco_enu);
  Float_t currentBinRange = GetBinRange(currentBin);
  Float_t currentBinStart = GetBinStart(currentBin);
  Float_t lookup = (currentBinStart+(currentBinRange/2.0));
  Float_t cutVal = cutsHist->GetBinContent(cutsHist->FindBin(lookup));
  if(QEIDval>cutVal) return 1;
  return 0;
}

//**********************************************************
void MadQEID::FillPDFs(Int_t cc_nc, Int_t process, MadBase* mb, Int_t evtindex, Float_t reco_enu, Float_t recoW2)
{
  static bool first=true;
  if(first){ InitPDFs(); first=false;}
  if(!CalcQEVars(mb,evtindex)) return;
  if(reco_enu>=20.0) return;
  Int_t currentBin = GetEnergyBin(reco_enu);
  recoWsqd = recoW2;

  if(cc_nc==1 && process==1001){
    rePHfracCCQE[currentBin]->Fill(rePHfrac);
    ntrkCCQE[currentBin]->Fill(ntrks);
    recoW2CCQE[currentBin]->Fill(recoWsqd);
    nshowCCQE[currentBin]->Fill(nshows);
    hpeakCCQE[currentBin]->Fill(houghPeak);
    PCACCQE[currentBin]->Fill(PCAvar);
  }
  if(cc_nc==1 && process==1002){ 
    rePHfracCCR[currentBin]->Fill(rePHfrac);
    ntrkCCR[currentBin]->Fill(ntrks);
    recoW2CCR[currentBin]->Fill(recoWsqd);
    nshowCCR[currentBin]->Fill(nshows);
    hpeakCCR[currentBin]->Fill(houghPeak);
    PCACCR[currentBin]->Fill(PCAvar);
    rePHfracBG[currentBin]->Fill(rePHfrac);
    ntrkBG[currentBin]->Fill(ntrks);
    recoW2BG[currentBin]->Fill(recoWsqd);
    nshowBG[currentBin]->Fill(nshows);
    hpeakBG[currentBin]->Fill(houghPeak);
    PCABG[currentBin]->Fill(PCAvar);
  }   
  if(cc_nc==1 && process==1003){
    rePHfracCCDIS[currentBin]->Fill(rePHfrac);
    ntrkCCDIS[currentBin]->Fill(ntrks);
    recoW2CCDIS[currentBin]->Fill(recoWsqd);
    nshowCCDIS[currentBin]->Fill(nshows);
    hpeakCCDIS[currentBin]->Fill(houghPeak);
    PCACCDIS[currentBin]->Fill(PCAvar);
    rePHfracBG[currentBin]->Fill(rePHfrac);
    ntrkBG[currentBin]->Fill(ntrks);
    recoW2BG[currentBin]->Fill(recoWsqd);
    nshowBG[currentBin]->Fill(nshows);
    hpeakBG[currentBin]->Fill(houghPeak);
    PCABG[currentBin]->Fill(PCAvar);
  }
  if(cc_nc==0){
    rePHfracNC[currentBin]->Fill(rePHfrac);
    ntrkNC[currentBin]->Fill(ntrks);
    recoW2NC[currentBin]->Fill(recoWsqd);
    nshowNC[currentBin]->Fill(nshows);
    hpeakNC[currentBin]->Fill(houghPeak);
    PCANC[currentBin]->Fill(PCAvar);
    rePHfracBG[currentBin]->Fill(rePHfrac);
    ntrkBG[currentBin]->Fill(ntrks);
    recoW2BG[currentBin]->Fill(recoWsqd);
    nshowBG[currentBin]->Fill(nshows);
    hpeakBG[currentBin]->Fill(houghPeak);
    PCABG[currentBin]->Fill(PCAvar);
  }
}

//********************************************************
void MadQEID::NormalisePDFs()
{
  Float_t integ;
  for(Int_t i=0;i<18;i++){

    integ=rePHfracCCQE[i]->Integral(0,((rePHfracCCQE[i]->GetNbinsX())+1));
    if(integ!=0) rePHfracCCQE[i]->Scale(1./integ);
    integ=rePHfracCCR[i]->Integral(0,((rePHfracCCR[i]->GetNbinsX())+1));
    if(integ!=0) rePHfracCCR[i]->Scale(1./integ);
    integ=rePHfracCCDIS[i]->Integral(0,((rePHfracCCDIS[i]->GetNbinsX())+1));
    if(integ!=0) rePHfracCCDIS[i]->Scale(1./integ);
    integ=rePHfracNC[i]->Integral(0,((rePHfracNC[i]->GetNbinsX())+1));
    if(integ!=0) rePHfracNC[i]->Scale(1./integ);
    integ=rePHfracBG[i]->Integral(0,((rePHfracBG[i]->GetNbinsX())+1));
    if(integ!=0) rePHfracBG[i]->Scale(1./integ);

    integ=recoW2CCQE[i]->Integral(0,((recoW2CCQE[i]->GetNbinsX())+1));
    if(integ!=0) recoW2CCQE[i]->Scale(1./integ);
    integ=recoW2CCR[i]->Integral(0,((recoW2CCR[i]->GetNbinsX())+1));
    if(integ!=0) recoW2CCR[i]->Scale(1./integ);
    integ=recoW2CCDIS[i]->Integral(0,((recoW2CCDIS[i]->GetNbinsX())+1));
    if(integ!=0) recoW2CCDIS[i]->Scale(1./integ);
    integ=recoW2NC[i]->Integral(0,((recoW2NC[i]->GetNbinsX())+1));
    if(integ!=0) recoW2NC[i]->Scale(1./integ);
    integ=recoW2BG[i]->Integral(0,((recoW2BG[i]->GetNbinsX())+1));
    if(integ!=0) recoW2BG[i]->Scale(1./integ);
    
    integ=ntrkCCQE[i]->Integral(0,((ntrkCCQE[i]->GetNbinsX())+1));
    if(integ!=0) ntrkCCQE[i]->Scale(1./integ);
    integ=ntrkCCR[i]->Integral(0,((ntrkCCR[i]->GetNbinsX())+1));
    if(integ!=0) ntrkCCR[i]->Scale(1./integ);
    integ=ntrkCCDIS[i]->Integral(0,((ntrkCCDIS[i]->GetNbinsX())+1));
    if(integ!=0) ntrkCCDIS[i]->Scale(1./integ);
    integ=ntrkNC[i]->Integral(0,((ntrkNC[i]->GetNbinsX())+1));
    if(integ!=0) ntrkNC[i]->Scale(1./integ);
    integ=ntrkBG[i]->Integral(0,((ntrkBG[i]->GetNbinsX())+1));
    if(integ!=0) ntrkBG[i]->Scale(1./integ);
    
    integ=nshowCCQE[i]->Integral(0,((nshowCCQE[i]->GetNbinsX())+1));
    if(integ!=0) nshowCCQE[i]->Scale(1./integ);
    integ=nshowCCR[i]->Integral(0,((nshowCCR[i]->GetNbinsX())+1));
    if(integ!=0) nshowCCR[i]->Scale(1./integ);
    integ=nshowCCDIS[i]->Integral(0,((nshowCCDIS[i]->GetNbinsX())+1));
    if(integ!=0) nshowCCDIS[i]->Scale(1./integ);
    integ=nshowNC[i]->Integral(0,((nshowNC[i]->GetNbinsX())+1));
    if(integ!=0) nshowNC[i]->Scale(1./integ);
    integ=nshowBG[i]->Integral(0,((nshowBG[i]->GetNbinsX())+1));
    if(integ!=0) nshowBG[i]->Scale(1./integ);
    
    integ=hpeakCCQE[i]->Integral(0,((hpeakCCQE[i]->GetNbinsX())+1));
    if(integ!=0) hpeakCCQE[i]->Scale(1./integ);
    integ=hpeakCCR[i]->Integral(0,((hpeakCCR[i]->GetNbinsX())+1));
    if(integ!=0) hpeakCCR[i]->Scale(1./integ);
    integ=hpeakCCDIS[i]->Integral(0,((hpeakCCDIS[i]->GetNbinsX())+1));
    if(integ!=0) hpeakCCDIS[i]->Scale(1./integ);
    integ=hpeakNC[i]->Integral(0,((hpeakNC[i]->GetNbinsX())+1));
    if(integ!=0) hpeakNC[i]->Scale(1./integ);
    integ=hpeakBG[i]->Integral(0,((hpeakBG[i]->GetNbinsX())+1));
    if(integ!=0) hpeakBG[i]->Scale(1./integ);

    integ=PCACCQE[i]->Integral(0,((PCACCQE[i]->GetNbinsX())+1));
    if(integ!=0) PCACCQE[i]->Scale(1./integ);
    integ=PCACCR[i]->Integral(0,((PCACCR[i]->GetNbinsX())+1));
    if(integ!=0) PCACCR[i]->Scale(1./integ);
    integ=PCACCDIS[i]->Integral(0,((PCACCDIS[i]->GetNbinsX())+1));
    if(integ!=0) PCACCDIS[i]->Scale(1./integ);
    integ=PCANC[i]->Integral(0,((PCANC[i]->GetNbinsX())+1));
    if(integ!=0) PCANC[i]->Scale(1./integ);
    integ=PCABG[i]->Integral(0,((PCABG[i]->GetNbinsX())+1));
    if(integ!=0) PCABG[i]->Scale(1./integ);
  }
  pdfs_normalised=kTRUE;
}

//********************************************************
void MadQEID::RebinPDFs(Int_t nbins)
{
  for(Int_t i=0;i<18;i++){
    rePHfracCCQE[i]->Rebin(nbins);
    rePHfracCCR[i]->Rebin(nbins);
    rePHfracCCDIS[i]->Rebin(nbins);
    rePHfracNC[i]->Rebin(nbins);
    rePHfracBG[i]->Rebin(nbins);
    recoW2CCQE[i]->Rebin(nbins);
    recoW2CCR[i]->Rebin(nbins);
    recoW2CCDIS[i]->Rebin(nbins);
    recoW2NC[i]->Rebin(nbins);
    recoW2BG[i]->Rebin(nbins);
    PCACCQE[i]->Rebin(nbins);
    PCACCR[i]->Rebin(nbins);
    PCACCDIS[i]->Rebin(nbins);
    PCANC[i]->Rebin(nbins);
    PCABG[i]->Rebin(nbins);
  }
}

//********************************************************
void MadQEID::WritePDFs(const char* name)
{
  TDirectory* dir = gDirectory;
  TFile* f = new TFile(name,"recreate");
  dir->cd();
  for(Int_t i=0;i<18;i++){
    rePHfracCCQE[i]->SetDirectory(f);
    rePHfracCCR[i]->SetDirectory(f);
    rePHfracCCDIS[i]->SetDirectory(f);
    rePHfracNC[i]->SetDirectory(f);
    rePHfracBG[i]->SetDirectory(f);
    recoW2CCQE[i]->SetDirectory(f);
    recoW2CCR[i]->SetDirectory(f);
    recoW2CCDIS[i]->SetDirectory(f);
    recoW2NC[i]->SetDirectory(f);
    recoW2BG[i]->SetDirectory(f);
    ntrkCCQE[i]->SetDirectory(f);
    ntrkCCR[i]->SetDirectory(f);
    ntrkCCDIS[i]->SetDirectory(f);
    ntrkNC[i]->SetDirectory(f);
    ntrkBG[i]->SetDirectory(f);
    nshowCCQE[i]->SetDirectory(f);
    nshowCCR[i]->SetDirectory(f);
    nshowCCDIS[i]->SetDirectory(f);
    nshowNC[i]->SetDirectory(f);
    nshowBG[i]->SetDirectory(f);
    hpeakCCQE[i]->SetDirectory(f);
    hpeakCCR[i]->SetDirectory(f);
    hpeakCCDIS[i]->SetDirectory(f);
    hpeakNC[i]->SetDirectory(f);
    hpeakBG[i]->SetDirectory(f);
    PCACCQE[i]->SetDirectory(f);
    PCACCR[i]->SetDirectory(f);
    PCACCDIS[i]->SetDirectory(f);
    PCANC[i]->SetDirectory(f);
    PCABG[i]->SetDirectory(f);
    if(trainingHists[i]) trainingHists[i]->SetDirectory(f);
  } 
  if(cutsHist) cutsHist->SetDirectory(f); 
  f->Write();
}

//********************************************************
void MadQEID::ReadPDFs(const char* name)
{
  TFile* QEfile = new TFile(name,"READ");
  if(QEfile->IsOpen() && !QEfile->IsZombie()){
    std::cout<<"Read MadQEID PDFs from : "<<name<<std::endl;
    QEfile->cd();
    for(Int_t i=0;i<18;i++){
      char name[30];

      sprintf(name,"rePHfracCCQE[%d]",i);
      rePHfracCCQE[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"rePHfracCCR[%d]",i);
      rePHfracCCR[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"rePHfracCCDIS[%d]",i);
      rePHfracCCDIS[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"rePHfracNC[%d]",i);
      rePHfracNC[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"rePHfracBG[%d]",i);
      rePHfracBG[i] = (TH1F*) QEfile->Get(name);

      sprintf(name,"recoW2CCQE[%d]",i);
      recoW2CCQE[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"recoW2CCR[%d]",i);
      recoW2CCR[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"recoW2CCDIS[%d]",i);
      recoW2CCDIS[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"recoW2NC[%d]",i);
      recoW2NC[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"recoW2BG[%d]",i);
      recoW2BG[i] = (TH1F*) QEfile->Get(name);

      sprintf(name,"hpeakCCQE[%d]",i);
      hpeakCCQE[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"hpeakCCR[%d]",i);
      hpeakCCR[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"hpeakCCDIS[%d]",i);
      hpeakCCDIS[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"hpeakNC[%d]",i);
      hpeakNC[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"hpeakBG[%d]",i);
      hpeakBG[i] = (TH1F*) QEfile->Get(name);

      sprintf(name,"ntrkCCQE[%d]",i);
      ntrkCCQE[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"ntrkCCR[%d]",i);
      ntrkCCR[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"ntrkCCDIS[%d]",i);
      ntrkCCDIS[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"ntrkNC[%d]",i);
      ntrkNC[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"ntrkBG[%d]",i);
      ntrkBG[i] = (TH1F*) QEfile->Get(name);

      sprintf(name,"nshowCCQE[%d]",i);
      nshowCCQE[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"nshowCCR[%d]",i);
      nshowCCR[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"nshowCCDIS[%d]",i);
      nshowCCDIS[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"nshowNC[%d]",i);
      nshowNC[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"nshowBG[%d]",i);
      nshowBG[i] = (TH1F*) QEfile->Get(name);

      sprintf(name,"PCACCQE[%d]",i);
      PCACCQE[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"PCACCR[%d]",i);
      PCACCR[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"PCACCDIS[%d]",i);
      PCACCDIS[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"PCANC[%d]",i);
      PCANC[i] = (TH1F*) QEfile->Get(name);
      sprintf(name,"PCABG[%d]",i);
      PCABG[i] = (TH1F*) QEfile->Get(name);

    }
  }
  else {
    std::cout<<"Failed to read MadQEID PDFs from : "<<name<<std::endl;
    assert(false);
  }
  NormalisePDFs();
  cutsHist = (TH1F*) QEfile->Get("cuthist");
}

//********************************************************
void MadQEID::FillCutTrainingHists(MadBase* mb, Int_t evtindex, Float_t reco_enu, Float_t recoW2)
{
  if(reco_enu>=20.0) return;
  static bool first=true;
  if(first){ InitTrainingHists(); first=false;}
  Float_t currentID = CalcQEID(mb,evtindex,reco_enu,recoW2);
  Int_t currentBin = GetEnergyBin(reco_enu);
  trainingHists[currentBin]->Fill(currentID);
}

//********************************************************
void MadQEID::AltFillCutTrainingHists(MadBase* mb, Int_t evtindex, Float_t reco_enu, Float_t recoW2)
{
  if(reco_enu>=20.0) return;
  static bool first=true;
  if(first){ InitTrainingHists(); first=false;}
  Float_t currentID = AltCalcQEID(mb,evtindex,reco_enu,recoW2);
  Int_t currentBin = GetEnergyBin(reco_enu);
  trainingHists[currentBin]->Fill(currentID);
}

//********************************************************
void MadQEID::MakeCutHist(Float_t eff)
{

  for(Int_t i=0;i<18;i++){
    Int_t totQE = static_cast<Int_t>(trainingHists[i]->GetEntries());
    Float_t currentBinStart = GetBinStart(i);
    Float_t currentBinRange = GetBinRange(i);
    Float_t lookup = (currentBinStart + (currentBinRange/2.0));
    if(totQE==0){
      cutsHist->Fill(lookup,0.3);
      //cout << "Training energy bin number: " << i << endl;
      //cout << "Cut value:                   0.3" << endl;
      //cout << endl;
      continue;
    }
    Int_t currentNumQE = 0;
    Float_t currentEff = 0.0;
    Bool_t foundEff = false;
    Int_t trainingHistBins = trainingHists[i]->GetNbinsX();

    //cout << "-------------------------------" << endl;
    //cout << "Training energy bin number: " << i << endl;
    //cout << "Lookup energy:              " << lookup << endl;
    //cout << endl; 

    for(Int_t bin=1;bin<trainingHistBins;bin++){
      if(foundEff) continue;
      Int_t binNumFromEnd = trainingHistBins-bin;
      currentNumQE+=static_cast<Int_t>(trainingHists[i]->GetBinContent(binNumFromEnd));
      if(currentNumQE==0){ 
        currentEff=0.0;
        continue;
      }
      else currentEff = (static_cast<Float_t>(currentNumQE)/static_cast<Float_t>(totQE));

      //cout << "Current QEID value:   " << trainingHists[i]->GetBinCenter(binNumFromEnd) << endl;
      //cout << "Total # QE:           " << totQE << endl;
      //cout << "Current # QE:         " << currentNumQE << endl;
      //cout << "Current efficiency:   " << currentEff << endl;
      //cout << endl;

      if(currentEff>eff && !foundEff){
        Float_t cutVal = static_cast<Float_t>((trainingHists[i]->GetBinCenter(binNumFromEnd)));
        cutsHist->Fill(lookup,cutVal);
        //cout << "Training energy bin number: " << i << endl;
        //cout << "Cut value:                  " << cutVal << endl;
        //cout << endl;
        foundEff = true;
      }
      if(bin==(trainingHistBins-1) && !foundEff){ 
        cutsHist->Fill(lookup,4.5); //still include definite QE in sample
        //cout << "Training energy bin number: " << i << endl;
        //cout << "Cut value:                   4.5" << endl;
        //cout << endl;
      }
    }
  }
}

//********************************************************
Float_t MadQEID::GetRePHfrac()
{
  return rePHfrac;
}

//********************************************************
Bool_t MadQEID::IsTrkHit(Float_t zPos,Float_t tPos, const NtpSRTrack* ntpTrk, MadBase* madb)
{
  for(Int_t i=0;i<ntpTrk->nstrip;i++){
    const NtpSRStrip* strp = madb->GetStrip(ntpTrk->stp[i]);
    if(!strp) continue;
    if(strp->tpos==tPos && strp->z==zPos){
      return true;
    }
  }
  return false;
}

//********************************************************
Bool_t MadQEID::IsShowHit(Float_t zPos,Float_t tPos, const NtpSRShower* ntpShow, MadBase* madb)
{
  for(Int_t i=0;i<ntpShow->nstrip;i++){
    const NtpSRStrip* strp = madb->GetStrip(ntpShow->stp[i]);
    if(!strp) continue;
    if(strp->tpos==tPos && strp->z==zPos){
      return true;
    }
  }
  return false;
}

//********************************************************
Bool_t MadQEID::IsVPln(Int_t plane)
{
  for(int i=0;i<246;i++){
    if(plane==2*i) return true;
  }
  return false;
}

//********************************************************
Int_t MadQEID::GetEnergyBin(Float_t enu)
{
  Int_t bin = 0;
  if(enu>=0.0 && enu<0.5) bin=0;
  if(enu>=0.5 && enu<1.0) bin=1;
  if(enu>=1.0 && enu<1.5) bin=2;
  if(enu>=1.5 && enu<2.0) bin=3;
  if(enu>=2.0 && enu<2.5) bin=4;
  if(enu>=2.5 && enu<3.0) bin=5;
  if(enu>=3.0 && enu<3.5) bin=6;
  if(enu>=3.5 && enu<4.0) bin=7;
  if(enu>=4.0 && enu<4.5) bin=8;
  if(enu>=4.5 && enu<5.0) bin=9;
  if(enu>=5.0 && enu<6.0) bin=10;
  if(enu>=6.0 && enu<7.0) bin=11;
  if(enu>=7.0 && enu<8.0) bin=12;
  if(enu>=8.0 && enu<10.0) bin=13;
  if(enu>=10.0 && enu<12.0) bin=14;
  if(enu>=12.0 && enu<14.0) bin=15;
  if(enu>=14.0 && enu<17.0) bin=16;
  if(enu>=17.0 && enu<20.0) bin=17;
  return bin;
}

//********************************************************
Float_t MadQEID::GetBinStart(Int_t bin)
{
  Float_t binstart=0.;
  if(bin==0) binstart=0.0;
  if(bin==1) binstart=0.5;
  if(bin==2) binstart=1.0;
  if(bin==3) binstart=1.5;
  if(bin==4) binstart=2.0;
  if(bin==5) binstart=2.5;
  if(bin==6) binstart=3.0;
  if(bin==7) binstart=3.5;
  if(bin==8) binstart=4.0;
  if(bin==9) binstart=4.5;
  if(bin==10) binstart=5.0;
  if(bin==11) binstart=6.0;
  if(bin==12) binstart=7.0;
  if(bin==13) binstart=8.0;
  if(bin==14) binstart=10.0;
  if(bin==15) binstart=12.0;
  if(bin==16) binstart=14.0;
  if(bin==17) binstart=17.0;
  return binstart;
}

//********************************************************
Float_t MadQEID::GetBinRange(Int_t bin)
{
  Float_t binrange = 0.;
  if(bin<=9){
    binrange=0.5;
  } 
  if(bin>9 && bin<=12){
    binrange=1.0;
  }
  if(bin>=13 && bin<16){
    binrange=2.0;
  }
  if(bin>=16){
    binrange=3.0;
  }
  return binrange;
}


//********************************************************

#endif

---