MadDpID.cxx

Go to the documentation of this file.

//
// MadDpID
//
// A class to calculate Dave Petyt's CC PID
//
// Created:  D. Petyt -- some time ago
//
// Transcribed: M. Kordosky, 11 August 2005
//
// $Author: rhatcher $ 
//
// $Revision: 1.17 $
// 
// $Name:  $
//
// $Id: MadDpID.cxx,v 1.17 2008/11/19 17:16:21 rhatcher Exp $
//
// $Log: MadDpID.cxx,v $
// Revision 1.17  2008/11/19 17:16:21  rhatcher
// replace non-standard "uint" with "unsigned int".
//
// Revision 1.16  2007/05/22 14:33:40  vahle
// Update MadDpID to deal with non-standard names of far pdf files
//
// Revision 1.15  2007/02/06 16:44:42  vahle
// Changed MadDpID::ChoosePDF to choose the le010z185i pdf file if the beam type is kLE
//
// Revision 1.14  2007/02/06 13:34:13  vahle
// Changed MadDpID so that the pid for the horn off beam comes from the HE pdf file
//
// Revision 1.13  2007/02/06 13:25:15  vahle
// Fixed reco_y, reco_x initialization bug, changed dpid PHCorrection to 1 if the recoversion is cedar
//
// Revision 1.12  2007/02/05 18:43:30  vahle
// Changed ChoosePDFs in MadDpID to take reconstruction and MC versions as arguments.  This way it can figure out what pdf file to use based on the reco and mc versions.  Should be backwards compatible, if these two new parameters are not specified, the code will assume they are birch/carrot and load the pdf files that were used for the prl analysis
//
// Revision 1.11  2007/01/30 03:17:21  rhatcher
// No longer user BeamType "class" found in MCReweight, but rather the
// BeamType "namespace" found in Conventions package.
// Also covert to using "Detector::" rather than "DetectorType::", though
// for now they are synonyms.
//
// Revision 1.10  2007/01/29 23:05:13  boehm
// Adjusting how the data files are loaded into the two PIDs.  Now the code looks first in the SRT_PRIVATE_CONTEXT and then checks the SRT_PUBLIC_CONTEXT.  This further detachs the CC-pids from the Mad Framework.
//
// Revision 1.9  2006/05/08 21:27:44  boehm
// Primary change is the addition of functions to calculate the PID's from NtpSt objects.  Both PIDs may now be generated without relying on the rest of the Mad infrastructure.  Mad based code though should not be changed.
//
// Also removed the if(!ntpTrack.fit.pass) return;  from MadDpID.cxx as it does not reflect the current state of the CC analysis. This was removed in the branched version but remained in the development version.
//
// Revision 1.8  2006/05/04 21:32:41  boehm
// Minor adjustment so that there now exists a CalcPid(var1, var2, var3) function
// which generic non-Mad based programs can use to calculate a DP-PID value.
//
// Revision 1.7  2006/03/10 18:57:21  kordosky
// correct MadDpID and MadMKAnalysis to account for 1.8% shift in pulseheight related PID quantities
//
// Revision 1.6  2006/02/15 19:59:51  kordosky
// remove fiducial requirement to match jan-06 branch.
//
// Revision 1.5  2005/10/07 11:19:12  kordosky
// deal with BeamMonSpill::StatusBits.beam_type==0 without asserting. Clearly, it can be zero!
//
// Revision 1.4  2005/10/06 18:40:52  kordosky
// random additions
//
// Revision 1.3  2005/10/06 15:03:26  kordosky
// various improvements, pids finally integrated into pan construction routine, still not fully debugged though
//
// Revision 1.2  2005/09/14 13:32:45  kordosky
// DpID and NsID classes now fully working and validated with MadTestAnalysis by comparing nannpid vs. annpid and npid0 vs. pid0.  pdf files and weights stored in data subdirectory. FD pdfs needed for DpID.
//
//

#include "Mad/MadDpID.h"
#include "Mad/MadQuantities.h"

#include <cassert>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>

#include "TH1.h"
#include "TFile.h"
#include "TSystem.h"
#include "string"

#include "MadDpAnalysis.h"

MadDpID::MadDpID() : 
  ph_correction(1.0),cache_det(Detector::kUnknown), cache_beam(BeamType::kUnknown), fLikeliHist(6){
  for (unsigned int i=0; i<6; i++) fLikeliHist[i]=0;
}

MadDpID::~MadDpID() {
  for (unsigned int i=0; i<fLikeliHist.size(); i++){
    if(fLikeliHist[i]) {
      // check that histogram really isn't in a directory
      // before deleting
      if(fLikeliHist[i]->GetDirectory()!=0){
        delete fLikeliHist[i];      
        fLikeliHist[i]=0;
      }
    }
  }
}

bool MadDpID::CalcVars(const NtpSRTrack* ntpTrack, const NtpSREvent* ntpEvent, 
         const NtpSREventSummary *eventSummary, Detector::Detector_t det,
         Int_t method, Float_t& var1, Float_t& var2, Float_t& var3)
{
  if(method != 0) return false;
                                                                                
  var1=var2=var3=0.0;
  if(ntpEvent==0 || ntpTrack==0) return false;
  if(eventSummary == 0) return false;                                                                              
  // basic checks on quality of event
//  if(! ntpTrack->fit.pass) return false;

  //Calculate the Variables
/*
  if(! MadDpAnalysis::InFidVol(ntpHeader->GetVldContext().GetDetector(),
                               ntpEvent->vtx.x, ntpEvent->vtx.y,
                               ntpEvent->vtx.z) ) return false;
  */
  // compute variables

  var1=eventSummary->plane.n;
  if(det == Detector::kNear){
    var1=ntpEvent->plane.end-ntpEvent->plane.beg+1;
  }

  Float_t phtrack=ntpTrack->ph.sigcor*ph_correction;
  Float_t phevent=ntpEvent->ph.sigcor*ph_correction;
                                                                                
  if (phevent>0) {var2=phtrack/phevent;}
                                                                                
  if (ntpTrack->plane.n!=0) var3 = float(ntpTrack->ph.sigcor*ph_correction)
                              /float(ntpTrack->plane.n);
  return true;
}

bool MadDpID::CalcVars(MadQuantities* mb, Int_t event, Int_t method,
                       Float_t& var1, Float_t& var2, Float_t& var3)
{
  // calculate the variables used to construct the PID
  //
  // this routine should be used in PID calculation *and* to fill PDFs
  //
  // checks that the event passes certain criteria before doing the calc  
  // returns false in case of error, or if the event doesn't satisfy
  // 

  if(method!=0) return false;
  
  var1=var2=var3=0.0;

  const RecCandHeader* ntpHeader = mb->GetHeader();
  const NtpSREvent* ntpEvent=mb->GetEvent(event);
  Int_t track = -1;
  const NtpSRTrack* ntpTrack=mb->GetLargestTrackFromEvent(event,track);
  if(track==-1) return false;

  return CalcVars(ntpTrack, ntpEvent, mb->GetEventSummary(), 
                  ntpHeader->GetVldContext().GetDetector(), method, 
                  var1, var2, var3);
}


Float_t MadDpID::CalcPID(MadQuantities* mb, Int_t event, Int_t method)
{
  // calculate a PID or return -10 in case the calculation of variables
  // encounters an error. The -10 maintains historical convention
  
  Float_t var1,var2,var3;
  if(!CalcVars(mb,event,method,var1,var2,var3) ) return -10;

  return CalcPID(var1, var2, var3);
}

Float_t MadDpID::CalcPID(const NtpSRTrack* ntpTrack, const NtpSREvent* ntpEvent,
         const NtpSREventSummary *eventSummary, Detector::Detector_t det,
         Int_t method)
{
  Float_t var1, var2, var3;
  
  if(!CalcVars(ntpTrack,ntpEvent,eventSummary,det,method,var1,var2,var3) )
        return -10;
                                                                                
  return CalcPID(var1, var2, var3);
}

Float_t MadDpID::CalcPID(Float_t var1, Float_t var2, Float_t var3)
{
  assert(fLikeliHist[0]); // essentially make sure file was read

  Float_t ProbNC = 1.;
  Float_t ProbCC = 1.;
  Float_t PidCC = 0.;

  Float_t prob1=fLikeliHist[0]->GetBinContent(fLikeliHist[0]->FindBin(var1));
  Float_t prob2=fLikeliHist[2]->GetBinContent(fLikeliHist[2]->FindBin(var2));
  Float_t prob3=fLikeliHist[4]->GetBinContent(fLikeliHist[4]->FindBin(var3));

  if (prob1==0) {prob1=0.0001;}
  if (prob2==0) {prob2=0.0001;}
  if (prob3==0) {prob3=0.0001;}
 
  ProbCC=prob1*prob2*prob3;

  prob1=fLikeliHist[1]->GetBinContent(fLikeliHist[1]->FindBin(var1));
  prob2=fLikeliHist[3]->GetBinContent(fLikeliHist[3]->FindBin(var2));
  prob3=fLikeliHist[5]->GetBinContent(fLikeliHist[5]->FindBin(var3));
 
  if (prob1==0) {prob1=0.0001;}
  if (prob2==0) {prob2=0.0001;}
  if (prob3==0) {prob3=0.0001;}
 
  ProbNC=prob1*prob2*prob3;

  PidCC=-(sqrt(-TMath::Log(ProbCC))-sqrt(-TMath::Log(ProbNC)));

  return PidCC;
}

bool MadDpID::ReadPDFs(const char* name){
  

  // check if this file or directory exists
  Long_t id,size,flags,modtime;
  if(gSystem->GetPathInfo(name,&id,&size,&flags,&modtime) == 1) {
    return false;
  }
  // check that it is a file and not a directory
  if(flags>1) return false;
  

  TDirectory* sd = gDirectory;
  TFile fLikeliFile(name,"READ");
  bool found=true;
  if(fLikeliFile.IsOpen() && !fLikeliFile.IsZombie()) { //if file is found
    
    fLikeliFile.cd();
    
    fLikeliHist[0]=dynamic_cast<TH1*>( fLikeliFile.Get("Event length - cc"));
    fLikeliHist[1]=dynamic_cast<TH1*>( fLikeliFile.Get("Event length - nc"));
    fLikeliHist[2]=dynamic_cast<TH1*>( fLikeliFile.Get("Track ph frac - cc"));
    fLikeliHist[3]=dynamic_cast<TH1*>( fLikeliFile.Get("Track ph frac - nc"));
    fLikeliHist[4]=dynamic_cast<TH1*>( fLikeliFile.Get("ph per plane (cc)"));
    fLikeliHist[5]=dynamic_cast<TH1*>( fLikeliFile.Get("ph per plane (nc)"));
    
    for (unsigned int i=0;i<fLikeliHist.size();i++){
      assert(fLikeliHist[i]);
      fLikeliHist[i]->SetDirectory(0); // histogram owned by this object
      float integ = fLikeliHist[i]->Integral(); 
      fLikeliHist[i]->Scale(1./integ);
    }
    std::cout<<"Read MadDpID pdfs from : "<<name<<std::endl;
  }
  else found=false;
  if(sd) sd->cd();
  return found;
}

void MadDpID::WritePDFs(const char* name){
  // write histograms to a file but retain ownership with this object
  TDirectory* sd = gDirectory;
  TFile f(name,"recreate");
  for (unsigned int i=0; i<fLikeliHist.size(); i++){
    fLikeliHist[i]->Write();
  }
  if(sd) sd->cd();
}


void MadDpID::InitPDFs(const Detector::Detector_t& det){
  // create pdf histograms and assign to fLikeliHist array

  int nbins=60;
  float low=0; 
  float high=300;
  if(det==Detector::kFar){ nbins=50; low=0; high=500;}
  
  TH1F *evlength_nc = new TH1F("Event length - nc",
                               "Event length - nc",
                               nbins,low,high);
  evlength_nc->SetXTitle("Event length (planes)");
  TH1F *evlength_cc = new TH1F("Event length - cc",
                               "Event length - cc",
                               nbins,low,high);
  evlength_cc->SetXTitle("Event length (planes)");


  TH1F *phfrac_nc = new TH1F("Track ph frac - nc",
                               "Track ph frac - nc",
                               50,0,1);
  phfrac_nc->SetXTitle("pulse height fraction");


  TH1F *phfrac_cc = new TH1F("Track ph frac - cc",
                               "Track ph frac - cc",
                               50,0,1);
  phfrac_cc->SetXTitle("pulse height fraction");


  TH1F *phplane_nc = new TH1F("ph per plane (nc)",
                              "ph per plane (nc)",
                               50,0,5000);
  phplane_nc->SetXTitle("pulse height per plane");
  TH1F *phplane_cc = new TH1F("ph per plane (cc)",
                              "ph per plane (cc)",
                               50,0,5000);
  phplane_cc->SetXTitle("pulse height per plane");


  fLikeliHist[0] = evlength_cc;
  fLikeliHist[1] = evlength_nc;
  
  fLikeliHist[2] = phfrac_cc;
  fLikeliHist[3] = phfrac_nc;
  
  fLikeliHist[4] = phplane_cc;
  fLikeliHist[5] = phplane_nc;

  for (unsigned int i=0; i<fLikeliHist.size(); i++) fLikeliHist[i]->SetDirectory(0);
  
}

bool MadDpID::FillPDFs(MadQuantities* mb, Int_t event, Int_t method, 
                       Float_t weight){
  
  // fill PDFs for this event
  // the code checks truth information to determine
  // if the event is numu-CC or NC
  

  assert(fLikeliHist[0]);// assure histograms are there

  Float_t var1,var2,var3;
  if(!CalcVars(mb,event,method,var1,var2,var3) ) return false;
  
  Int_t mcevent=-1;
  const NtpMCTruth* ntpTruth = mb->GetTruthForReco(event,mcevent);
  if(mcevent == -1  || ntpTruth==0) return false;

  // check binning of histograms, possibly rebin 
  // not in love with this bit of the code!
  const RecCandHeader* header = mb->GetHeader();
  assert(header); // something deeply wrong
  CheckBinning(header->GetVldContext().GetDetector());
  
  //
  Int_t cc_nc = mb->IAction(mcevent);
  Int_t flavour = mb->Flavour(mcevent);
  
  if(flavour==2 && cc_nc==1){
    fLikeliHist[0]->Fill(var1,weight);
    fLikeliHist[2]->Fill(var2,weight);
    fLikeliHist[4]->Fill(var3,weight);
  }
  else if(cc_nc==0) {
    fLikeliHist[1]->Fill(var1,weight);
    fLikeliHist[3]->Fill(var2,weight);
    fLikeliHist[5]->Fill(var3,weight);
  }
  
  return true;
}


void MadDpID::CheckBinning(const Detector::Detector_t& det){
  // duplicate a funny little bit of code, which rebins
  // one of the pdfs in reaction to the detector type
  TH1* evlength_cc = fLikeliHist[0];
  TH1* evlength_nc = fLikeliHist[1];
  assert(evlength_cc);
  assert(evlength_nc);

  if(det==Detector::kFar){
    if (evlength_cc->GetNbinsX()==60) {
      evlength_cc->SetBins(50,0,500);
      evlength_nc->SetBins(50,0,500);
    }
  }
  else if(det==Detector::kNear){
    if (evlength_cc->GetBinLowEdge(50)>200) {
        evlength_cc->SetBins(60,0,300);
        evlength_nc->SetBins(60,0,300);
    }
  }
  
}

void MadDpID::ResetPDFs(){
  for (unsigned int i=0; i<fLikeliHist.size(); i++){
    TH1* h = fLikeliHist[i];
    delete h; fLikeliHist[i]=0;
  }
}

bool MadDpID::ChoosePDFs(const Detector::Detector_t& det, 
                         const BeamType::BeamType_t& beam,
                         std::string recover, std::string mcver)
{
  // check to see if we already have this file
  if(det==cache_det && beam==cache_beam){
    // no action is needed 
    return true;
  }
  if(det==Detector::kUnknown || beam==BeamType::kUnknown){
    return false;
  }
  std::string priv="";
  std::string pub="";

  std::string base="";

  priv=getenv("SRT_PRIVATE_CONTEXT");
  pub = getenv("SRT_PUBLIC_CONTEXT");
 
  if(priv == "" && pub == ""){
    std::cerr<<"No SRT_PUBLIC_CONTEXT set"<<std::endl;
    assert(false);
  }

  
  priv += "/Mad/data/";
  pub += "/Mad/data/";

  struct stat ss;
  if(stat(priv.c_str(), &ss) == -1) {
   std::cout<<"Data not present in SRT_PRIVATE_CONTEXT trying SRT_PUBLIC"<<std::endl;
   if(stat(pub.c_str(), &ss) == -1) {
     std::cout<<"Data not present in SRT_PUBLIC - doesn't seem to exist"<<std::endl;
      assert(false);
   }else{ base = pub; }
  }else { base = priv; }

  base+="dp_pdf_";

  if((recover==""||recover.find("birch")<recover.size()||
     recover.find("Birch")<recover.size()||
     recover.find("BIRCH")<recover.size()||
     recover.find("R1.18")<recover.size())&&
     (mcver==""||mcver.find("carrot")<mcver.size()||
      mcver.find("Carrot")<mcver.size()||
      mcver.find("CARROT")<mcver.size())){
    if(det==Detector::kNear){
      base+="near";
    }
    else if(det==Detector::kFar){
      base+="far";
    }
    //  else return false;
    
    if(beam==BeamType::kLE || beam==BeamType::k010){
      base+="_le";
    }
    else if(beam==BeamType::kME){
      base+="_me";
    }
    else if(beam==BeamType::kHE){
      base+="_he";
    }
  }
  else if((mcver.find("daikon")<mcver.size()||
          mcver.find("Daikon")<mcver.size()||
          mcver.find("DAIKON")<mcver.size())&&
          (recover.find("cedar")<recover.size()||
           recover.find("Cedar")<recover.size()||
           recover.find("CEDAR")<recover.size()||
           recover.find("1.24")<recover.size())){
    BeamType::BeamType_t tmpbeam = beam;
    switch(beam){
    case BeamType::kL010z170i: tmpbeam = BeamType::kL010z185i; break;
    case BeamType::kL010z200i: tmpbeam = BeamType::kL010z185i; break;
      //forcing le to go to le010
    case BeamType::kL000z200i: tmpbeam = BeamType::kL010z185i; break;
    case BeamType::kL010z185i_lowintensity: 
      tmpbeam = BeamType::kL010z185i; break;
      //for birch, we used le010 pdfs for the horn off running,
      //but, I think we should use the he since the
      //mean neutrino energy in the horn off running is 
      //higher than le, more like he.
    case BeamType::kL010z000i: tmpbeam = BeamType::kL250z200i; break;
    default: tmpbeam = beam; break;
    }
    if(det==Detector::kNear){
      base+="near";
      base+=("_"+(string)BeamType::AsString(tmpbeam));
    }
    else if(det==Detector::kFar){
      base+="far";
      if(tmpbeam==BeamType::kL250z200i){
        base+="_phe";
      }
      else{
        base+="_le";
      }
    }

    base+="_cedar_daikon";
  }

  //  else return false;
  
  base+=".root";
  if(stat(base.c_str(),&ss) == -1) {
    std::cerr<<"The file '"<<base<<"' doesn't seem to exist"<<std::endl;
    //    assert(false);
    return false;
  }
  cout<<"Reading pdfs from "<<base<<endl;
  //  assert(ReadPDFs(base.c_str()));
  if(ReadPDFs(base.c_str())){
  // set these so we don't have to reread each event
    cache_det=det;
    cache_beam=beam;    
    return true;
  }
  else {
    return false;
  }
}

---