MadDpID Class Reference

#include <MadDpID.h>

List of all members.

Public Member Functions

 MadDpID ()
 ~MadDpID ()
Float_t CalcPID (MadQuantities *mb, Int_t event, Int_t method)
Float_t CalcPID (const NtpSRTrack *ntpTrack, const NtpSREvent *ntpEvent, const NtpSREventSummary *eventSummary, Detector::Detector_t det, Int_t method)
Float_t CalcPID (Float_t var1, Float_t var2, Float_t var3)
bool ReadPDFs (const char *name)
bool ChoosePDFs (const Detector::Detector_t &det, const BeamType::BeamType_t &beam, std::string recover="", std::string mcver="")
void InitPDFs (const Detector::Detector_t &det)
bool FillPDFs (MadQuantities *mb, Int_t event, Int_t method, Float_t weight=1)
void WritePDFs (const char *name)
void SetPHCorrection (double p=1.0)

Private Member Functions

void ResetPDFs ()
void CheckBinning (const Detector::Detector_t &det)
bool CalcVars (MadQuantities *mb, Int_t event, Int_t method, Float_t &var1, Float_t &var2, Float_t &var3)
bool 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)

Private Attributes

double ph_correction
Detector::Detector_t cache_det
BeamType::BeamType_t cache_beam
std::vector< TH1 * > fLikeliHist

Detailed Description

Definition at line 71 of file MadDpID.h.


Constructor & Destructor Documentation

MadDpID::MadDpID (  ) 

Definition at line 96 of file MadDpID.cxx.

References fLikeliHist.

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

MadDpID::~MadDpID (  ) 

Definition at line 101 of file MadDpID.cxx.

References fLikeliHist, and DirectoryHelpers::GetDirectory().

00101                   {
00102   for (unsigned int i=0; i<fLikeliHist.size(); i++){
00103     if(fLikeliHist[i]) {
00104       // check that histogram really isn't in a directory
00105       // before deleting
00106       if(fLikeliHist[i]->GetDirectory()!=0){
00107         delete fLikeliHist[i];      
00108         fLikeliHist[i]=0;
00109       }
00110     }
00111   }
00112 }


Member Function Documentation

Float_t MadDpID::CalcPID ( Float_t  var1,
Float_t  var2,
Float_t  var3 
)

Definition at line 199 of file MadDpID.cxx.

References fLikeliHist.

00200 {
00201   assert(fLikeliHist[0]); // essentially make sure file was read
00202 
00203   Float_t ProbNC = 1.;
00204   Float_t ProbCC = 1.;
00205   Float_t PidCC = 0.;
00206 
00207   Float_t prob1=fLikeliHist[0]->GetBinContent(fLikeliHist[0]->FindBin(var1));
00208   Float_t prob2=fLikeliHist[2]->GetBinContent(fLikeliHist[2]->FindBin(var2));
00209   Float_t prob3=fLikeliHist[4]->GetBinContent(fLikeliHist[4]->FindBin(var3));
00210 
00211   if (prob1==0) {prob1=0.0001;}
00212   if (prob2==0) {prob2=0.0001;}
00213   if (prob3==0) {prob3=0.0001;}
00214  
00215   ProbCC=prob1*prob2*prob3;
00216 
00217   prob1=fLikeliHist[1]->GetBinContent(fLikeliHist[1]->FindBin(var1));
00218   prob2=fLikeliHist[3]->GetBinContent(fLikeliHist[3]->FindBin(var2));
00219   prob3=fLikeliHist[5]->GetBinContent(fLikeliHist[5]->FindBin(var3));
00220  
00221   if (prob1==0) {prob1=0.0001;}
00222   if (prob2==0) {prob2=0.0001;}
00223   if (prob3==0) {prob3=0.0001;}
00224  
00225   ProbNC=prob1*prob2*prob3;
00226 
00227   PidCC=-(sqrt(-TMath::Log(ProbCC))-sqrt(-TMath::Log(ProbNC)));
00228 
00229   return PidCC;
00230 }

Float_t MadDpID::CalcPID ( const NtpSRTrack ntpTrack,
const NtpSREvent ntpEvent,
const NtpSREventSummary eventSummary,
Detector::Detector_t  det,
Int_t  method 
)

Definition at line 187 of file MadDpID.cxx.

References CalcPID(), and CalcVars().

00190 {
00191   Float_t var1, var2, var3;
00192   
00193   if(!CalcVars(ntpTrack,ntpEvent,eventSummary,det,method,var1,var2,var3) )
00194         return -10;
00195                                                                                 
00196   return CalcPID(var1, var2, var3);
00197 }

Float_t MadDpID::CalcPID ( MadQuantities mb,
Int_t  event,
Int_t  method 
)

Definition at line 176 of file MadDpID.cxx.

References CalcVars().

Referenced by MNtpModule::Ana(), ANtpAnalysisInfoAna::Analyze(), MuonRemovalInfoAna::Analyze(), AnalysisInfoAna::Analyze(), CalcPID(), NuAnalysis::ChargeSignCut(), MadTestAnalysis::CreatePAN(), MadMKAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), NuAnalysis::Efficiencies(), NuAnalysis::EnergySpect(), NuPIDInterface::GetDpID(), and NuAnalysis::N_1().

00177 {
00178   // calculate a PID or return -10 in case the calculation of variables
00179   // encounters an error. The -10 maintains historical convention
00180   
00181   Float_t var1,var2,var3;
00182   if(!CalcVars(mb,event,method,var1,var2,var3) ) return -10;
00183 
00184   return CalcPID(var1, var2, var3);
00185 }

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 
) [private]

Definition at line 114 of file MadDpID.cxx.

References NtpSRPlane::beg, NtpSRPlane::end, Detector::kNear, NtpSRPlane::n, NtpSRTrack::ph, NtpSREvent::ph, ph_correction, NtpSRTrack::plane, NtpSREvent::plane, NtpSREventSummary::plane, and NtpSRPulseHeight::sigcor.

00117 {
00118   if(method != 0) return false;
00119                                                                                 
00120   var1=var2=var3=0.0;
00121   if(ntpEvent==0 || ntpTrack==0) return false;
00122   if(eventSummary == 0) return false;                                                                              
00123   // basic checks on quality of event
00124 //  if(! ntpTrack->fit.pass) return false;
00125 
00126   //Calculate the Variables
00127 /*
00128   if(! MadDpAnalysis::InFidVol(ntpHeader->GetVldContext().GetDetector(),
00129                                ntpEvent->vtx.x, ntpEvent->vtx.y,
00130                                ntpEvent->vtx.z) ) return false;
00131   */
00132   // compute variables
00133 
00134   var1=eventSummary->plane.n;
00135   if(det == Detector::kNear){
00136     var1=ntpEvent->plane.end-ntpEvent->plane.beg+1;
00137   }
00138 
00139   Float_t phtrack=ntpTrack->ph.sigcor*ph_correction;
00140   Float_t phevent=ntpEvent->ph.sigcor*ph_correction;
00141                                                                                 
00142   if (phevent>0) {var2=phtrack/phevent;}
00143                                                                                 
00144   if (ntpTrack->plane.n!=0) var3 = float(ntpTrack->ph.sigcor*ph_correction)
00145                               /float(ntpTrack->plane.n);
00146   return true;
00147 }

bool MadDpID::CalcVars ( MadQuantities mb,
Int_t  event,
Int_t  method,
Float_t &  var1,
Float_t &  var2,
Float_t &  var3 
) [private]

Definition at line 149 of file MadDpID.cxx.

References VldContext::GetDetector(), MadBase::GetEvent(), MadBase::GetEventSummary(), MadBase::GetHeader(), MadBase::GetLargestTrackFromEvent(), and RecHeader::GetVldContext().

Referenced by CalcPID(), and FillPDFs().

00151 {
00152   // calculate the variables used to construct the PID
00153   //
00154   // this routine should be used in PID calculation *and* to fill PDFs
00155   //
00156   // checks that the event passes certain criteria before doing the calc  
00157   // returns false in case of error, or if the event doesn't satisfy
00158   // 
00159 
00160   if(method!=0) return false;
00161   
00162   var1=var2=var3=0.0;
00163 
00164   const RecCandHeader* ntpHeader = mb->GetHeader();
00165   const NtpSREvent* ntpEvent=mb->GetEvent(event);
00166   Int_t track = -1;
00167   const NtpSRTrack* ntpTrack=mb->GetLargestTrackFromEvent(event,track);
00168   if(track==-1) return false;
00169 
00170   return CalcVars(ntpTrack, ntpEvent, mb->GetEventSummary(), 
00171                   ntpHeader->GetVldContext().GetDetector(), method, 
00172                   var1, var2, var3);
00173 }

void MadDpID::CheckBinning ( const Detector::Detector_t det  )  [private]

Definition at line 377 of file MadDpID.cxx.

References fLikeliHist, Detector::kFar, and Detector::kNear.

Referenced by FillPDFs().

00377                                                        {
00378   // duplicate a funny little bit of code, which rebins
00379   // one of the pdfs in reaction to the detector type
00380   TH1* evlength_cc = fLikeliHist[0];
00381   TH1* evlength_nc = fLikeliHist[1];
00382   assert(evlength_cc);
00383   assert(evlength_nc);
00384 
00385   if(det==Detector::kFar){
00386     if (evlength_cc->GetNbinsX()==60) {
00387       evlength_cc->SetBins(50,0,500);
00388       evlength_nc->SetBins(50,0,500);
00389     }
00390   }
00391   else if(det==Detector::kNear){
00392     if (evlength_cc->GetBinLowEdge(50)>200) {
00393         evlength_cc->SetBins(60,0,300);
00394         evlength_nc->SetBins(60,0,300);
00395     }
00396   }
00397   
00398 }

bool MadDpID::ChoosePDFs ( const Detector::Detector_t det,
const BeamType::BeamType_t beam,
std::string  recover = "",
std::string  mcver = "" 
)

Definition at line 407 of file MadDpID.cxx.

References BeamType::AsString(), base, cache_beam, cache_det, BeamType::k010, Detector::kFar, BeamType::kHE, BeamType::kL000z200i, BeamType::kL010z000i, BeamType::kL010z170i, BeamType::kL010z185i, BeamType::kL010z185i_lowintensity, BeamType::kL010z200i, BeamType::kL250z200i, BeamType::kLE, BeamType::kM000z200i_nova, BeamType::kM000z200i_nova_rev, BeamType::kME, Detector::kNear, Detector::kUnknown, BeamType::kUnknown, ReadPDFs(), and stat.

Referenced by MNtpModule::Ana(), ANtpAnalysisInfoAna::Analyze(), MuonRemovalInfoAna::Analyze(), AnalysisInfoAna::Analyze(), NuAnalysis::ChargeSignCut(), MadMKAnalysis::CreatePAN(), MadTVAnalysis::CreatePAN(), NuAnalysis::Efficiencies(), NuAnalysis::EnergySpect(), NuPIDInterface::InitialiseDpID(), NuAnalysis::N_1(), and NuAnalysis::NuMuBarAppearance().

00410 {
00411   // check to see if we already have this file
00412   if(det==cache_det && beam==cache_beam){
00413     // no action is needed 
00414     return true;
00415   }
00416   if(det==Detector::kUnknown || beam==BeamType::kUnknown){
00417     return false;
00418   }
00419   std::string priv="";
00420   std::string pub="";
00421 
00422   std::string base="";
00423 
00424   priv=getenv("SRT_PRIVATE_CONTEXT");
00425   pub = getenv("SRT_PUBLIC_CONTEXT");
00426  
00427   if(priv == "" && pub == ""){
00428     std::cerr<<"No SRT_PUBLIC_CONTEXT set"<<std::endl;
00429     assert(false);
00430   }
00431 
00432   
00433   priv += "/Mad/data/";
00434   pub += "/Mad/data/";
00435 
00436   struct stat ss;
00437   if(stat(priv.c_str(), &ss) == -1) {
00438    std::cout<<"Data not present in SRT_PRIVATE_CONTEXT trying SRT_PUBLIC"<<std::endl;
00439    if(stat(pub.c_str(), &ss) == -1) {
00440      std::cout<<"Data not present in SRT_PUBLIC - doesn't seem to exist"<<std::endl;
00441       assert(false);
00442    }else{ base = pub; }
00443   }else { base = priv; }
00444 
00445   base+="dp_pdf_";
00446 
00447   if((recover==""||recover.find("birch")<recover.size()||
00448      recover.find("Birch")<recover.size()||
00449      recover.find("BIRCH")<recover.size()||
00450      recover.find("R1.18")<recover.size())&&
00451      (mcver==""||mcver.find("carrot")<mcver.size()||
00452       mcver.find("Carrot")<mcver.size()||
00453       mcver.find("CARROT")<mcver.size())){
00454     if(det==Detector::kNear){
00455       base+="near";
00456     }
00457     else if(det==Detector::kFar){
00458       base+="far";
00459     }
00460     //  else return false;
00461     
00462     if(beam==BeamType::kLE || beam==BeamType::k010){
00463       base+="_le";
00464     }
00465     else if(beam==BeamType::kME){
00466       base+="_me";
00467     }
00468     else if(beam==BeamType::kHE){
00469       base+="_he";
00470     }
00471   }
00472   else if((mcver.find("daikon")<mcver.size()||
00473           mcver.find("Daikon")<mcver.size()||
00474           mcver.find("DAIKON")<mcver.size())&&
00475           (recover.find("cedar")<recover.size()||
00476            recover.find("Cedar")<recover.size()||
00477            recover.find("CEDAR")<recover.size()||
00478            recover.find("1.24")<recover.size())){
00479     BeamType::BeamType_t tmpbeam = beam;
00480     switch(beam){
00481     case BeamType::kL010z170i: tmpbeam = BeamType::kL010z185i; break;
00482     case BeamType::kL010z200i: tmpbeam = BeamType::kL010z185i; break;
00483       //forcing le to go to le010
00484     case BeamType::kL000z200i: tmpbeam = BeamType::kL010z185i; break;
00485     case BeamType::kL010z185i_lowintensity: 
00486       tmpbeam = BeamType::kL010z185i; break;
00487       //for birch, we used le010 pdfs for the horn off running,
00488       //but, I think we should use the he since the
00489       //mean neutrino energy in the horn off running is 
00490       //higher than le, more like he.
00491     case BeamType::kL010z000i: tmpbeam = BeamType::kL250z200i; break;
00492     case BeamType::kM000z200i_nova: tmpbeam = BeamType::kL250z200i; break;
00493     case BeamType::kM000z200i_nova_rev: tmpbeam = BeamType::kL250z200i; break;
00494     default: tmpbeam = beam; break;
00495     }
00496     if(det==Detector::kNear){
00497       base+="near";
00498       base+=("_"+(string)BeamType::AsString(tmpbeam));
00499     }
00500     else if(det==Detector::kFar){
00501       base+="far";
00502       if(tmpbeam==BeamType::kL250z200i){
00503         base+="_phe";
00504       }
00505       else{
00506         base+="_le";
00507       }
00508     }
00509 
00510     base+="_cedar_daikon";
00511   }
00512 
00513   //  else return false;
00514   
00515   base+=".root";
00516   if(stat(base.c_str(),&ss) == -1) {
00517     std::cerr<<"The file '"<<base<<"' doesn't seem to exist"<<std::endl;
00518     //    assert(false);
00519     return false;
00520   }
00521   cout<<"Reading pdfs from "<<base<<endl;
00522   //  assert(ReadPDFs(base.c_str()));
00523   if(ReadPDFs(base.c_str())){
00524   // set these so we don't have to reread each event
00525     cache_det=det;
00526     cache_beam=beam;    
00527     return true;
00528   }
00529   else {
00530     return false;
00531   }
00532 }

bool MadDpID::FillPDFs ( MadQuantities mb,
Int_t  event,
Int_t  method,
Float_t  weight = 1 
)

Definition at line 335 of file MadDpID.cxx.

References CalcVars(), CheckBinning(), MadQuantities::Flavour(), fLikeliHist, VldContext::GetDetector(), MadBase::GetHeader(), MadBase::GetTruthForReco(), RecHeader::GetVldContext(), header, and MadQuantities::IAction().

00336                                       {
00337   
00338   // fill PDFs for this event
00339   // the code checks truth information to determine
00340   // if the event is numu-CC or NC
00341   
00342 
00343   assert(fLikeliHist[0]);// assure histograms are there
00344 
00345   Float_t var1,var2,var3;
00346   if(!CalcVars(mb,event,method,var1,var2,var3) ) return false;
00347   
00348   Int_t mcevent=-1;
00349   const NtpMCTruth* ntpTruth = mb->GetTruthForReco(event,mcevent);
00350   if(mcevent == -1  || ntpTruth==0) return false;
00351 
00352   // check binning of histograms, possibly rebin 
00353   // not in love with this bit of the code!
00354   const RecCandHeader* header = mb->GetHeader();
00355   assert(header); // something deeply wrong
00356   CheckBinning(header->GetVldContext().GetDetector());
00357   
00358   //
00359   Int_t cc_nc = mb->IAction(mcevent);
00360   Int_t flavour = mb->Flavour(mcevent);
00361   
00362   if(flavour==2 && cc_nc==1){
00363     fLikeliHist[0]->Fill(var1,weight);
00364     fLikeliHist[2]->Fill(var2,weight);
00365     fLikeliHist[4]->Fill(var3,weight);
00366   }
00367   else if(cc_nc==0) {
00368     fLikeliHist[1]->Fill(var1,weight);
00369     fLikeliHist[3]->Fill(var2,weight);
00370     fLikeliHist[5]->Fill(var3,weight);
00371   }
00372   
00373   return true;
00374 }

void MadDpID::InitPDFs ( const Detector::Detector_t det  ) 

Definition at line 282 of file MadDpID.cxx.

References fLikeliHist, and Detector::kFar.

00282                                                    {
00283   // create pdf histograms and assign to fLikeliHist array
00284 
00285   int nbins=60;
00286   float low=0; 
00287   float high=300;
00288   if(det==Detector::kFar){ nbins=50; low=0; high=500;}
00289   
00290   TH1F *evlength_nc = new TH1F("Event length - nc",
00291                                "Event length - nc",
00292                                nbins,low,high);
00293   evlength_nc->SetXTitle("Event length (planes)");
00294   TH1F *evlength_cc = new TH1F("Event length - cc",
00295                                "Event length - cc",
00296                                nbins,low,high);
00297   evlength_cc->SetXTitle("Event length (planes)");
00298 
00299 
00300   TH1F *phfrac_nc = new TH1F("Track ph frac - nc",
00301                                "Track ph frac - nc",
00302                                50,0,1);
00303   phfrac_nc->SetXTitle("pulse height fraction");
00304 
00305 
00306   TH1F *phfrac_cc = new TH1F("Track ph frac - cc",
00307                                "Track ph frac - cc",
00308                                50,0,1);
00309   phfrac_cc->SetXTitle("pulse height fraction");
00310 
00311 
00312   TH1F *phplane_nc = new TH1F("ph per plane (nc)",
00313                               "ph per plane (nc)",
00314                                50,0,5000);
00315   phplane_nc->SetXTitle("pulse height per plane");
00316   TH1F *phplane_cc = new TH1F("ph per plane (cc)",
00317                               "ph per plane (cc)",
00318                                50,0,5000);
00319   phplane_cc->SetXTitle("pulse height per plane");
00320 
00321 
00322   fLikeliHist[0] = evlength_cc;
00323   fLikeliHist[1] = evlength_nc;
00324   
00325   fLikeliHist[2] = phfrac_cc;
00326   fLikeliHist[3] = phfrac_nc;
00327   
00328   fLikeliHist[4] = phplane_cc;
00329   fLikeliHist[5] = phplane_nc;
00330 
00331   for (unsigned int i=0; i<fLikeliHist.size(); i++) fLikeliHist[i]->SetDirectory(0);
00332   
00333 }

bool MadDpID::ReadPDFs ( const char *  name  ) 

Definition at line 232 of file MadDpID.cxx.

References fLikeliHist, gSystem(), id, and size.

Referenced by ChoosePDFs().

00232                                       {
00233   
00234 
00235   // check if this file or directory exists
00236   Long_t id,size,flags,modtime;
00237   if(gSystem->GetPathInfo(name,&id,&size,&flags,&modtime) == 1) {
00238     return false;
00239   }
00240   // check that it is a file and not a directory
00241   if(flags>1) return false;
00242   
00243 
00244   TDirectory* sd = gDirectory;
00245   TFile fLikeliFile(name,"READ");
00246   bool found=true;
00247   if(fLikeliFile.IsOpen() && !fLikeliFile.IsZombie()) { //if file is found
00248     
00249     fLikeliFile.cd();
00250     
00251     fLikeliHist[0]=dynamic_cast<TH1*>( fLikeliFile.Get("Event length - cc"));
00252     fLikeliHist[1]=dynamic_cast<TH1*>( fLikeliFile.Get("Event length - nc"));
00253     fLikeliHist[2]=dynamic_cast<TH1*>( fLikeliFile.Get("Track ph frac - cc"));
00254     fLikeliHist[3]=dynamic_cast<TH1*>( fLikeliFile.Get("Track ph frac - nc"));
00255     fLikeliHist[4]=dynamic_cast<TH1*>( fLikeliFile.Get("ph per plane (cc)"));
00256     fLikeliHist[5]=dynamic_cast<TH1*>( fLikeliFile.Get("ph per plane (nc)"));
00257     
00258     for (unsigned int i=0;i<fLikeliHist.size();i++){
00259       assert(fLikeliHist[i]);
00260       fLikeliHist[i]->SetDirectory(0); // histogram owned by this object
00261       float integ = fLikeliHist[i]->Integral(); 
00262       fLikeliHist[i]->Scale(1./integ);
00263     }
00264     std::cout<<"Read MadDpID pdfs from : "<<name<<std::endl;
00265   }
00266   else found=false;
00267   if(sd) sd->cd();
00268   return found;
00269 }

void MadDpID::ResetPDFs (  )  [private]

Definition at line 400 of file MadDpID.cxx.

References fLikeliHist.

00400                        {
00401   for (unsigned int i=0; i<fLikeliHist.size(); i++){
00402     TH1* h = fLikeliHist[i];
00403     delete h; fLikeliHist[i]=0;
00404   }
00405 }

void MadDpID::SetPHCorrection ( double  p = 1.0  )  [inline]
void MadDpID::WritePDFs ( const char *  name  ) 

Definition at line 271 of file MadDpID.cxx.

References fLikeliHist.

00271                                        {
00272   // write histograms to a file but retain ownership with this object
00273   TDirectory* sd = gDirectory;
00274   TFile f(name,"recreate");
00275   for (unsigned int i=0; i<fLikeliHist.size(); i++){
00276     fLikeliHist[i]->Write();
00277   }
00278   if(sd) sd->cd();
00279 }


Member Data Documentation

Definition at line 107 of file MadDpID.h.

Referenced by ChoosePDFs().

Definition at line 106 of file MadDpID.h.

Referenced by ChoosePDFs().

std::vector<TH1*> MadDpID::fLikeliHist [private]
double MadDpID::ph_correction [private]

Definition at line 104 of file MadDpID.h.

Referenced by CalcVars(), and SetPHCorrection().


The documentation for this class was generated from the following files:

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1