NueUtilities.cxx

Go to the documentation of this file.

// $Id: NueUtilities.cxx,v 1.2 2009/09/29 11:42:16 pawloski Exp $
//
// The NueUtilities namespace is intended to contain common functions
//  and utilities that analyzers frequently find themselves implementing 
//  For example:
//       Chaining files together
//       Getting the total POT for the chained files
//       Looping over events in a chain
//
// Author: Gregory Pawloski 
// Created: August 28, 2009

#include "NueAna/NueUtilities.h"

#include "MessageService/MsgService.h"

CVSID("");

 TChain* NueUtilities::CreateChain(std::string treeName,std::string fileListName)
 {
  TChain *chain= new TChain(treeName.c_str());
  ifstream fileList(fileListName.c_str());
  if( fileList )
  {
   while( fileList.good() && !fileList.eof() )
   {
     std::string line;
     getline(fileList,line,'\n');
     if ( line.size()==0 ) continue;
     chain->Add(line.c_str());
   }
  }
  return chain;
 }


 Double_t NueUtilities::GetChainPOT(TChain *chain)
 {
  NuePOT *nuepot = NULL;
  chain->SetBranchAddress("NuePOT",&nuepot);
 
  Double_t pot = 0;
  for( int i = 0; i<chain->GetEntries(); ++i )
  {
    chain->GetEntry(i);
    pot += nuepot->pot;
  }
  chain->ResetBranchAddresses();

  return(pot);
 }


 NueUtilities::AnaNueProcessor::AnaNueProcessor()
 {
   filelist = "";
   ana_nue_chain = NULL;
   pottree_chain = NULL;

   totalPOTNotFilled=true;

   nueRecord = NULL;
   numEntries = 0;
   currentEntry = -1;
   lastfraction = -1;

   ProcessedMC = false;
   ProcessedData = false;
   ProcessedFar = false;
   ProcessedNear = false;
   ProcessedNC = false;
   ProcessedCC = false;
   ProcessedBeamNue = false;
   ProcessedTau = false;
   ProcessedSignal = false;
   ProcessedRun1 = false;
   ProcessedRun2 = false;
   ProcessedRun3 = false;
 }

 NueUtilities::AnaNueProcessor::AnaNueProcessor(std::string fileListName)
 {
   filelist = fileListName;
   ana_nue_chain= NueUtilities::CreateChain("ana_nue",fileListName);
   pottree_chain= NueUtilities::CreateChain("pottree",fileListName);

   totalPOTNotFilled=true;

   nueRecord = NULL;
   ana_nue_chain->SetBranchAddress("NueRecord",&nueRecord);
   
   numEntries = ana_nue_chain->GetEntries();
   currentEntry = -1;
   lastfraction = -1;

   ProcessedMC = false;
   ProcessedData = false;
   ProcessedFar = false;
   ProcessedNear = false;
   ProcessedNC = false;
   ProcessedCC = false;
   ProcessedBeamNue = false;
   ProcessedTau = false;
   ProcessedSignal = false;
   ProcessedRun1 = false;
   ProcessedRun2 = false;
   ProcessedRun3 = false;
 }


 //Copy constructor
 NueUtilities::AnaNueProcessor::AnaNueProcessor(const NueUtilities::AnaNueProcessor &original)
 {
   *this = original;
 }

 NueUtilities::AnaNueProcessor & NueUtilities::AnaNueProcessor::operator=(const NueUtilities::AnaNueProcessor &original)   
 {  
  // check for "self assignment" and do nothing in that case
  if (this != &original)
  {
   filelist = original.GetFileListName();

   ana_nue_chain= NueUtilities::CreateChain("ana_nue",filelist);
   pottree_chain= NueUtilities::CreateChain("pottree",filelist);

   totalPOTNotFilled=true;

   nueRecord = NULL;
   ana_nue_chain->SetBranchAddress("NueRecord",&nueRecord);

   ProcessedMC = original.didMC();
   ProcessedData = original.didData();
   ProcessedFar = original.didFar();
   ProcessedNear = original.didNear();
   ProcessedNC = original.didNC();
   ProcessedCC = original.didCC();
   ProcessedBeamNue = original.didBeamNue();
   ProcessedTau = original.didTau();
   ProcessedSignal = original.didSignal();
   ProcessedRun1 = original.didRun1();
   ProcessedRun2 = original.didRun2();
   ProcessedRun3 = original.didRun3();
   
   numEntries = ana_nue_chain->GetEntries();
   lastfraction = original.GetLastFraction();

   this->SetEntry(original.CurrentEntry());
  }// if not itself
  return *this;

 }
 


 NueUtilities::AnaNueProcessor::~AnaNueProcessor()
 {
   delete ana_nue_chain;
   delete pottree_chain;

 }

 Double_t NueUtilities::AnaNueProcessor::GetTotalPOT()
 {
  if(totalPOTNotFilled)
  {
   totalPOT = NueUtilities::GetChainPOT(pottree_chain);
   totalPOTNotFilled = false;
  }
  return (totalPOT);
  
 }

 Bool_t NueUtilities::AnaNueProcessor::NextEntry()
 {
  ++currentEntry;
  return (SetEntry(currentEntry));
 }

 Bool_t NueUtilities::AnaNueProcessor::SetEntry(Long64_t entry)
 {
  currentEntry = entry;
  if(currentEntry <0 || currentEntry >= numEntries)
  {
    nueRecord= NULL;
    return(false);
  }
  ana_nue_chain->GetEntry(currentEntry);
   ProcessedMC |= isMC();
   ProcessedData |= isData();
   ProcessedFar |= isFar();
   ProcessedNear |= isNear();
   ProcessedNC |= isNC();
   ProcessedCC |= isCC();
   ProcessedBeamNue |= isBeamNue();
   ProcessedTau |= isTau();
   ProcessedSignal |= isSignal();
   ProcessedRun1 |= isRun1();
   ProcessedRun2 |= isRun2();
   ProcessedRun3 |= isRun3();

  //Could of had this function return the result of GetEntry(currentEntry),
  //but that would just cause loops to end and people would probably think that
  //they just hit the end of the entries. However, at this point in the code if
  //GetEntry(currentEntry) returns a nonpositve value, then something wrong happend
  //and I want the program to crash

  return(true);
  
 }

 void NueUtilities::AnaNueProcessor::PrintProgress(unsigned int percentInterval)
 {
   Int_t fraction = (100*currentEntry)/numEntries;
   
   if ( (fraction-lastfraction) >= (Int_t)percentInterval )
   {
    std::cout<<"\n"<<fraction<<"% done"<<std::endl;
    lastfraction = fraction;
   }

   if( currentEntry == (numEntries-1) ) 
   {
    std::cout<<"\nOn last entry"<<std::endl;
    lastfraction = fraction;
   }
 }

 void NueUtilities::AnaNueProcessor::PrintFileNames()
 {
  Long64_t num = 0;
  std::cout << "\nProcessing the following files:";
  ifstream fileList(filelist.c_str());
  if( fileList )
  {
   while( fileList.good() && !fileList.eof() )
   {
     std::string line;
     getline(fileList,line,'\n');
     if ( line.size()==0 ) continue;
     std::cout << "\n " << line;
     ++num;
   }
  }
  std::cout << "\nTotal number of files = " << num << endl;
 }


 Bool_t NueUtilities::AnaNueProcessor::isMC()
 {
  return(nueRecord->GetHeader().GetVldContext().GetSimFlag() == SimFlag::kMC);
 }

 Bool_t NueUtilities::AnaNueProcessor::isData()
 {
  return(nueRecord->GetHeader().GetVldContext().GetSimFlag() == SimFlag::kData);
 }

 Bool_t NueUtilities::AnaNueProcessor::isFar()
 {
  return(nueRecord->GetHeader().GetVldContext().GetDetector() == Detector::kFar);
 }

 Bool_t NueUtilities::AnaNueProcessor::isNear()
 {
  return(nueRecord->GetHeader().GetVldContext().GetDetector() == Detector::kNear);
 }

 Bool_t NueUtilities::AnaNueProcessor::isNC()
 {
  return(isMC() && nueRecord->mctrue.interactionType==0);
 }

 Bool_t NueUtilities::AnaNueProcessor::isCC()
 {
  return(isMC() && !isNC() && abs(nueRecord->mctrue.nuFlavor)==14);
 }

 Bool_t NueUtilities::AnaNueProcessor::isBeamNue()
 {
  return(isMC() && !isNC() && abs(nueRecord->mctrue.nuFlavor)==12 && abs(nueRecord->mctrue.nonOscNuFlavor)==12);
 }

 Bool_t NueUtilities::AnaNueProcessor::isTau()
 {
  return(isMC() && !isNC() && abs(nueRecord->mctrue.nuFlavor)==16);
 }

 Bool_t NueUtilities::AnaNueProcessor::isSignal()
 {
  return( isMC() && !isNC() && abs(nueRecord->mctrue.nuFlavor)==12 && abs(nueRecord->mctrue.nonOscNuFlavor)==14 );
 }

 Bool_t NueUtilities::AnaNueProcessor::isRun1()
 {
  return( NueStandard::IsRun1(nueRecord) );
 }

 Bool_t NueUtilities::AnaNueProcessor::isRun2()
 {
  return( NueStandard::IsRun2(nueRecord) );
 }

 Bool_t NueUtilities::AnaNueProcessor::isRun3()
 {
  return( NueStandard::IsRun3(nueRecord) );
 }

 Double_t NueUtilities::AnaNueProcessor::Normalization_RunSeparatedSample()
 {
  //normalize POT assuming sample only contains events from one run period
   Double_t eventWeight =1;
   if(isNear())
   {
     //Normalize Near Data and MC
      eventWeight = NueStandard::kNormalizedNearPOT/GetTotalPOT();
   }
   else if(isMC())
   {
     //Only Normalize Far MC
     //This normalization is only appropriate if Run 1, 2, and 3 are in separate samples
     if( NueStandard::IsRun1(nueRecord) )      eventWeight = NueStandard::kNormalizedFarPOT_Run1/GetTotalPOT();
     else if( NueStandard::IsRun2(nueRecord) ) eventWeight = NueStandard::kNormalizedFarPOT_Run2/GetTotalPOT();
     else if( NueStandard::IsRun3(nueRecord) ) eventWeight = NueStandard::kNormalizedFarPOT_Run3/GetTotalPOT();
   }
    
   return(eventWeight);
 }
 
 Double_t NueUtilities::AnaNueProcessor::Normalization_MixedRunSample()
 {
  //normalize POT assuming sample only contains mixture of events from different run periods
   Double_t eventWeight = 1;
   if(isNear())
   {
     //Normalize Near Data and MC
      eventWeight = NueStandard::kNormalizedNearPOT/GetTotalPOT();
   }
   else if(isMC())
   {
     //Only Normalize Far MC
     //This normalization is only appropriate if Run 1, 2, and 3 are in separate samples
     eventWeight = NueStandard::kNormalizedFarPOT/GetTotalPOT();
   }
    
   return(eventWeight);
 }

 Double_t NueUtilities::AnaNueProcessor::GetOscWeight_f210f213f214Separate()
 {
  //This method of applying oscillations depends on how you do add your samples together
  //This method works if you keep the f210*, f213*, and f214* files separate
   Double_t eventWeight =1;

   if( isMC() && isFar() && !isNC() )
   {
    eventWeight = NueStandard::GetOscWeight(nueRecord->mctrue.nuFlavor, nueRecord->mctrue.nonOscNuFlavor, nueRecord->mctrue.nuEnergy);
   }

   return(eventWeight);
 }


 Double_t NueUtilities::AnaNueProcessor::GetOldANN11()
 {
   return( NueStandard::GetPIDValue(nueRecord, Selection::kANN2PE) );
 }

 Double_t NueUtilities::AnaNueProcessor::GetANN11()
 {
   return( NueStandard::GetPIDValue(nueRecord, Selection::kANN2PE_DAIKON04) );
 }

 Double_t NueUtilities::AnaNueProcessor::GetANN14()
 {
   return( NueStandard::GetPIDValue(nueRecord, Selection::kANN14_DAIKON04) );
 }

 Double_t NueUtilities::AnaNueProcessor::GetSSPID()
 {
   return( NueStandard::GetPIDValue(nueRecord, Selection::kSSPID) );
 }

 Double_t NueUtilities::AnaNueProcessor::GetParticlePID()
 {
   return( NueStandard::GetPIDValue(nueRecord, Selection::kParticlePID) );
 }

 Double_t NueUtilities::AnaNueProcessor::GetRecoE()
 {
   return(nueRecord->srevent.phNueGeV);
 }


 Double_t NueUtilities::AnaNueProcessor::GetTruncatedOldANN11()
 {
   Double_t truncatedPID = GetOldANN11();
   if( truncatedPID < -0.049)
   {
    truncatedPID = -0.049; 
   }
   else if(truncatedPID > 0.99)
   {
    truncatedPID = 0.99; 
   }
   return(truncatedPID);
 }

 Double_t NueUtilities::AnaNueProcessor::GetTruncatedANN11()
 {
   Double_t truncatedPID = GetANN11();
   if( truncatedPID < -0.049)
   {
    truncatedPID = -0.049; 
   }
   else if(truncatedPID > 0.99)
   {
    truncatedPID = 0.99; 
   }
   return(truncatedPID);
 }

 Double_t NueUtilities::AnaNueProcessor::GetTruncatedANN14()
 {
   Double_t truncatedPID = GetANN14();
   if( truncatedPID < -0.049)
   {
    truncatedPID = -0.049; 
   }
   else if(truncatedPID > 0.99)
   {
    truncatedPID = 0.99; 
   }
   return(truncatedPID);
 }

 Double_t NueUtilities::AnaNueProcessor::GetTruncatedSSPID()
 {
   Double_t truncatedPID = GetSSPID();
   if( truncatedPID < -0.049)
   {
    truncatedPID = -0.049; 
   }
   else if(truncatedPID > 0.99)
   {
    truncatedPID = 0.99; 
   }
   return(truncatedPID);
 }

 Double_t NueUtilities::AnaNueProcessor::GetTruncatedParticlePID()
 {
   Double_t truncatedPID = GetParticlePID();
   if( truncatedPID < -0.049)
   {
    truncatedPID = -0.049; 
   }
   else if(truncatedPID > 0.99)
   {
    truncatedPID = 0.99; 
   }
   return(truncatedPID);
 }

---