CondensedNtpModule.cxx

Go to the documentation of this file.

//$Id: CondensedNtpModule.cxx,v 1.35 2007/02/12 21:25:37 asousa Exp $
//
//CondensedNtpModule.cxx
//
//Module for making analysis trees from SR ntuples
//
//B. Rebel 10/2004

#include "AnalysisNtuples/Module/CondensedNtpModule.h"
#include "AnalysisNtuples/Module/ANtpRecoNtpManipulator.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
#include "JobControl/JobCommand.h"
#include "Validity/VldTimeStamp.h" 
#include "CandNtupleSR/NtpSRRecord.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "CandNtupleSR/NtpSREventSummary.h"
#include "CandNtupleSR/NtpSRSlice.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "CandNtupleSR/NtpSRShower.h"
#include "CandNtupleSR/NtpSRShieldStrip.h"
#include "CandNtupleSR/NtpSRShieldSummary.h"
#include "CandNtupleSR/NtpSRCosmicRay.h"
#include "CandNtupleSR/NtpSRDmxStatus.h"
#include "MCNtuple/NtpMCRecord.h"
#include "MCNtuple/NtpMCTruth.h"
#include "MCNtuple/NtpMCStdHep.h"
#include "TruthHelperNtuple/NtpTHRecord.h"
#include "TruthHelperNtuple/NtpTHEvent.h"
#include "TruthHelperNtuple/NtpTHTrack.h"
#include "TruthHelperNtuple/NtpTHShower.h"
#include "TruthHelperNtuple/NtpTHStrip.h"
#include "Record/RecCandHeader.h"
#include "Conventions/Detector.h"
#include "Conventions/SimFlag.h"

#include "TFolder.h"
#include "TDirectory.h"
#include "TObjArray.h"

#include <cassert>
#include <algorithm>

ClassImp(CondensedNtpModule)

CVSID("$Id: CondensedNtpModule.cxx,v 1.35 2007/02/12 21:25:37 asousa Exp $");

// Declare this module to JobControll. Arguments are:
//  (1) The class name 
//  (2) The human-readable name 
//  (3) A short, human-readable description of what the module does
JOBMODULE(CondensedNtpModule, 
          "CondensedNtpModule",
          "A module used for making analysis tress from SR ntuples");

//......................................................................
CondensedNtpModule::CondensedNtpModule() :
  fFileName("analysisNtuple.root"),
  fTreeName("analysisNtuple"),
  fNtpFile(0),
  fNtuple(0),
  fDetector(1),
  fEventInfo(0),
  fHeaderInfo(0),
  fShowerInfo(0),
  fTrackInfo(0),
  fTruthInfo(0),
  fStripInfoPresent(1),
  fDataType(1),
  fReconstructed(1)
{

  MSG("JobC", Msg::kDebug) << "CondensedNtpModule::Constructor" << endl;
  
  fHeaderInfo = new ANtpHeaderInfo();
  fEventInfo = new ANtpEventInfo();
  fTrackInfo = new ANtpTrackInfo();
  fShowerInfo = new ANtpShowerInfo();
  fTruthInfo = new ANtpTruthInfo();
  
  //instantiate the ANtpInfoObjectFiller object
  fInfoFiller = new ANtpInfoObjectFiller();
        
}

//----------------------------------------------------------------------
CondensedNtpModule::~CondensedNtpModule()
{
  
  MSG("JobC", Msg::kDebug) << "CondensedNtpModule::Destructor" << endl;
  
  if(fHeaderInfo) delete fHeaderInfo;     
  if(fEventInfo) delete fEventInfo;       
  if(fShowerInfo) delete fShowerInfo;     
  if(fTrackInfo) delete fTrackInfo;       
  if(fTruthInfo) delete fTruthInfo;
    
}

//---------------------------------------------------------------------- 
void CondensedNtpModule::BeginJob()
{
  MSG("CondensedNtpModule", Msg::kDebug) 
    << "In BeginJob\n" << fFileName.c_str() 
    << " " << fTreeName.c_str() << endl;

  // Tell the InfoFillerObject which detector it is working on.
  // Module uses an int:     0-->NearDet & 1-->FarDet 
  // (Detector_t uses kNear = 1, kFar = 2) 
    switch (fDetector)
    {
    case 0:
      fInfoFiller->SetDetector(Detector::kNear);
      break;
    case 1:
      fInfoFiller->SetDetector(Detector::kFar);
      break;
    default:
      fInfoFiller->SetDetector(Detector::kUnknown); //Nothing is certain
      break;
    }
 
  // save the current working directory
  TDirectory* savedir  = gDirectory;  
  
  // create the new TFile for holding the electronics tree
  // this changes the value of gDirectoy
  fNtpFile = new TFile(fFileName.c_str(),"RECREATE");  
  
  // create TTree, these will attach themselves to the current 
  // working directory   
  fNtuple = new TTree(fTreeName.c_str(), "Analysis Tree");
  
  // the branches are all objects of type ANtp*Info*
  fNtuple->Branch("header.", "ANtpHeaderInfo", &fHeaderInfo, 64000, 2);
  fNtuple->Branch("event.", "ANtpEventInfo", &fEventInfo, 64000, 2);
  fNtuple->Branch("shower.", "ANtpShowerInfo", &fShowerInfo, 64000, 2);
  fNtuple->Branch("track.", "ANtpTrackInfo", &fTrackInfo, 64000, 2);
  if(fDataType>0){
    fNtuple->Branch("truth.", "ANtpTruthInfo", &fTruthInfo, 64000, 2);  
  } 
  
  // change back to current working directory before leaving constructor
  savedir->cd();   
  
  return;
}

//----------------------------------------------------------------------
//here is where you would loop over the events in the snarl to fill the 
//analysis tree
JobCResult CondensedNtpModule::Ana(const MomNavigator *mom)
{
  MSG("CondensedNtpModule", Msg::kDebug) << "start ana" << endl;
  
  JobCResult result(JobCResult::kPassed);
  
  //get the records from MOM
  assert(mom);
  
  NtpSRRecord *record = dynamic_cast<NtpSRRecord *>(mom->GetFragment("NtpSRRecord"));
  NtpStRecord *stRecord = dynamic_cast<NtpStRecord *>(mom->GetFragment("NtpStRecord"));
  
  if(!record && !stRecord){
    MSG("CondensedNtpModuleAtm", Msg::kWarning) << "Could not get NtpSR or NtpSt Record from MOM" << endl;
    result.SetWarning().SetFailed();
    return result;
  }//end if no NtpSRRecord
  
  //if there is no NtpMCRecord, then mom returns a null pointer
  NtpMCRecord *mcRecord = 0;
  mcRecord = dynamic_cast<NtpMCRecord *>(mom->GetFragment("NtpMCRecord"));
  if(!mcRecord&&fDataType&&!stRecord){
    MSG("CondensedNtpModule", Msg::kWarning) << "Could not get NtpMCRecord from MOM " 
                                             << " - hopefully you werent expecting one " << endl;
    result.SetWarning().SetFailed();
    return result;
  }
  
  NtpTHRecord *thRecord = 0;
  thRecord = dynamic_cast<NtpTHRecord *>(mom->GetFragment("NtpTHRecord"));
  
  //fill the header information
  Detector::Detector_t det = Detector::kFar;
  if(fDetector<1) det = Detector::kNear;
  //get a sim flag
  SimFlag::SimFlag_t dataType = SimFlag::kData;
  if(fDataType) dataType = SimFlag::kMC;
  
  //instantiate a ANtpRecoNtpManipulator object to help you get the info you want
  ANtpRecoNtpManipulator *ntpManipulator = 0;
  if(record) ntpManipulator = new ANtpRecoNtpManipulator(record, mcRecord, thRecord);
  else if(stRecord) ntpManipulator = new ANtpRecoNtpManipulator(stRecord);

  fInfoFiller->SetStripArray(ntpManipulator->GetStripArray());
  fInfoFiller->FillHeaderInformation(ntpManipulator, det, dataType, fHeaderInfo);  

  //set up which flags you want to use to determine the primary shower or track
  //a value of 0 for a flag means it will not be used
  ntpManipulator->SetPrimaryTrackCriteria(0,1,0); // nplanes, length, total pulse height
  ntpManipulator->SetPrimaryShowerCriteria(0,1); // nplanes, total pulse height
  
  //mc data because fDataType = 1, but no reconstructed event because 
  //record doesnt exist.  still want to have the mc info in the tree though 
  //to understand events that werent reconstructed
  if((mcRecord&&!record&&fDataType) 
     || (stRecord&&fDataType&&ntpManipulator->GetSnarlEventSummary().nevent<1)){
    
    MSG("CondensedNtpModule", Msg::kDebug) << "why am i here?" << endl;
    
    fReconstructed = 0;
    
    //flag the run, subRun values to indicated that 
    //this snarl was not reconstructed.  keep track of 
    //how many times that happens
    FillUnRecoedMCInformation(ntpManipulator->GetMCArray(), ntpManipulator->GetStdHepArray());
    fNtuple->Fill();
    ResetTreeVariables();
    
    return result;              
  }
  
  fReconstructed = 1;
  
  //-----------------------------------------------------------------------------------
  //here is where you start looping over the events in the ntuple to fill your analysis tree
  
  //set up some pointers to typical NtpSR objects
  NtpSRTrack *ntpTrack = 0;
  NtpSRStrip *ntpStrip = 0;
  NtpSREvent *ntpEvent = 0;
  NtpSRShower *ntpShower = 0;
  NtpMCTruth *ntpMCTruth = 0;
  NtpMCStdHep *ntpMCStdHep = 0;
  NtpTHEvent *ntpTHEvent = 0;
  //NtpTHStrip *ntpTHStrip = 0;
  //NtpTHTrack *ntpTHTrack = 0;
  //NtpTHShower *ntpTHShower = 0;
  
  //loop over the events in the snarl - it is also possible to just loop over 
  //showers, tracks or slices in the record using a loop like
  //for(Int_t i = 0; i < record->evthdr.ntrack; ++i){
  //  ntpTrack = ntpManipulator->GetSnarlTrack(i);
  //}
  //
  //you can get the primary track/shower from the snarl using 
  //ntpManipulator->GetSnarlPrimaryTrack()/GetSnarlPrimaryShower()
  
  for(Int_t i = 0; i < ntpManipulator->GetSnarlEventSummary().nevent; ++i){
    
    MSG("CondensedNtpModule", Msg::kDebug) << "on event " << i << endl;
    
    fEventInfo->index = i;
    
    //get event i.  the call to NtpHelper::GetEventInSnarl sets the 
    //NtpSREvent data member in that object to the current event so that
    //you can later call for the primary track, shower, etc.
    ntpManipulator->GetEventManipulator()->SetEventInSnarl(i);
    ntpEvent = ntpManipulator->GetEventManipulator()->GetEvent();               
    if(ntpEvent) fInfoFiller->FillEventInformation(ntpManipulator, 
                                                   ntpEvent, 
                                                   fEventInfo);
    
    //get the primary track for the event - if no track is present it 
    //returns 0
    ntpTrack = ntpManipulator->GetEventManipulator()->GetPrimaryTrack();                
    if(ntpTrack) fInfoFiller->FillTrackInformation(ntpTrack, fTrackInfo);
    
    //get the primary shower for the event - if no shower is present it 
    //returns 0
    ntpShower = ntpManipulator->GetEventManipulator()->GetPrimaryShower();              
    if(ntpShower) fInfoFiller->FillShowerInformation(ntpShower, fShowerInfo);
    
    //why not get a strip as well - if there is a strip in the event, the 
    //following line will get it
    if(fStripInfoPresent) ntpStrip = ntpManipulator->GetEventManipulator()->GetStrip(0);
    
    MSG("CondensedNtpModule", Msg::kDebug) << "CondensedNtpModule::Ana mcRecord="<<mcRecord << endl;
    MSG("CondensedNtpModule", Msg::kDebug) << "CondensedNtpModule::Ana thRecord="<<thRecord << endl;
    //check to see if there is MC information around and get the NtpMCTruth and
    //NtpTHEvent objects

    //the GetNtpTHEvent method takes as its argument the same index as the 
    //GetEventInSnarl method, ie the index of the current event.  you can 
    //then use the NtpTHEvent object returned to figure out which index to 
    //give the GetNtpMCTruth and GetNtpMCStdHep methods to return the appropriate 
    //truth information corresponding to this event.  what is going on behind the 
    //scenes is that there are arrays of NtpMCTruth and NtpMCStdHep objects in 
    //the NtpMCRecord but to figure out which entry is the appropriate one you need 
    //to have the NtpTHRecord array of thevt which is filled in the same order as 
    //the evt array in the NtpSRRecord.
    ntpTHEvent = ntpManipulator->GetMCManipulator()->GetNtpTHEvent(i);
    if(ntpTHEvent){
      ntpMCTruth = ntpManipulator->GetMCManipulator()->GetNtpMCTruth(ntpTHEvent->neumc);        
      ntpMCStdHep = ntpManipulator->GetMCManipulator()->GetNtpMCStdHep(ntpTHEvent->neustdhep);  
    }
    
    //can also do the above using NtpTHStrips, NtpTHShowers or NtpTHTracks.  get the index 
    //of the strip, track or shower in the NtpSRRecord stp, trk or shw arrays first, then 
    //get the NtpTHTrack/Shower object using the NtpHelper.  have to uncomment the variable 
    //declarations in lines 223-225 to use this
    //uncomment out the disired lines below
    /*
    
      if(ntpStrip) ntpTHStrip = ntpManipulator->GetMCManipulator()->GetNtpTHStrip(ntpStrip->index);
      if(ntpTHStrip){
      ntpMCTruth = ntpManipulator->GetMCManipulator()->GetNtpMCTruth(ntpTHStrip->neumc);        
      ntpMCStdHep = ntpManipulator->GetMCManipulator()->GetNtpMCStdHep(ntpTHStrip->neustdhep);  
      }
    
      if(ntpTrack) ntpTHTrack = ntpManipulator->GetMCManipulator()->GetNtpTHTrack((Int_t)ntpTrack->index);
      if(ntpTHTrack){
      ntpMCTruth = ntpManipulator->GetMCManipulator()->GetNtpMCTruth(ntpTHTrack->neumc);        
      ntpMCStdHep = ntpManipulator->GetMCManipulator()->GetNtpMCStdHep(ntpTHTrack->neustdhep);  
      }
      if(ntpShower) ntpTHShower = ntpManipulator->GetMCManipulator()->GetNtpTHShower((Int_t)ntpShower->index);   
      if(ntpTHShower){
      ntpMCTruth = ntpManipulator->GetMCManipulator()->GetNtpMCTruth(ntpTHShower->neumc);       
      ntpMCStdHep = ntpManipulator->GetMCManipulator()->GetNtpMCStdHep(ntpTHShower->neustdhep); 
      }
    */
      
    //got the NtpMCTruth object now, so fill the desired info
    MSG("CondensedNtpModule", Msg::kDebug) << "CondensedNtpModule::Ana ntpMCTruth="<<ntpMCTruth<< endl;
    if(ntpMCTruth){
      fInfoFiller->FillMCTruthInformation(ntpMCTruth, fTruthInfo);
    }                           
  
    
    fNtuple->Fill();
    ResetTreeVariables();
  }//end loop over events for this snarl
  //-----------------------------------------------------------------------------------

  delete ntpManipulator;

  return result;
}

//......................................................................
void CondensedNtpModule::Help() 
{
  MSG("JobC", Msg::kInfo) 
    << "NearElectronicsCheck::Help\n"
    <<"CondensedNtpModule is a module which demultiplexes events "
    <<"in the far detector."
    << endl;
}

//----------------------------------------------------------------------
void CondensedNtpModule::EndJob() 
{
  
  MSG("CondensedNtpModule", Msg::kInfo) << "start end job method" << endl;      
  
  //save current working directory
  TDirectory *savedir = gDirectory;
  
  //cd to the "directory" of the file
  fNtpFile->cd();

  fNtuple->Write();
  
  //go back to original working directory
  savedir->cd();  
  
  fNtpFile->Write();
  fNtpFile->Close();
  
  return;
}

//......................................................................
const Registry& CondensedNtpModule::DefaultConfig() const
{
  
  int itrue = 1;  // work around for lack of bool in registry
  int ifalse = 0; // work around for lack of bool in registry
  
  static Registry r;
  
  r.UnLockValues();
  
  r.Set("FileName",           "analysisTree.root");
  r.Set("TreeName",           "analysisTree");
  r.Set("EndPlane",           485);
  r.Set("Detector",           1);
  r.Set("StripInfoPresent",   1);
  r.Set("DataType",           1);
  
  r.LockValues();
  
  //quiet compiler warnings
  itrue = ifalse;
  
  return r;
}

//......................................................................
void CondensedNtpModule::Config(const Registry& r)
{
  int         tmpb;  // a temp bool. See comment under DefaultConfig...
  int         tmpi;  // a temp int.
  double      tmpd;  // a temp double.
  const char* tmps;  // a temp string
  
  if (r.Get("FileName",           tmps)) fFileName              =  tmps;
  if (r.Get("TreeName",           tmps)) fTreeName              =  tmps;
  if (r.Get("EndPlane",           tmpi)) fEndPlane              =  tmpi;
  if (r.Get("Detector",           tmpi)) fDetector              =  tmpi;
  if (r.Get("StripInfoPresent",   tmpi)) fStripInfoPresent      =  tmpi;
  if (r.Get("DataType",           tmpi)) fDataType              =  tmpi;
  
  //quiet compiler warnings
  tmpb = 0;
  tmpd = 0.;
  
  return;
}

//----------------------------------------------------------------------
void CondensedNtpModule::FillUnRecoedMCInformation(TClonesArray *mcArray, TClonesArray *stdArray)
{
  NtpMCTruth *ntpMCTruth = 0;
  NtpMCStdHep *ntpMCStdHep = 0;
  
  //get the number of entries in the arrays
  Int_t numMCEntries = mcArray->GetLast()+1;
  Int_t stdHepLow = 0;
  Int_t stdHepHigh = 0;
  
  //loop over the individual events - only one in the far detector, 
  //many in the near
  for(Int_t i = 0; i < numMCEntries; ++i){
    
    fEventInfo->index = i;
    
    ntpMCTruth = dynamic_cast<NtpMCTruth *>(mcArray->At(i));
    if(ntpMCTruth) fInfoFiller->FillMCTruthInformation(ntpMCTruth, fTruthInfo);
    
    stdHepLow = ntpMCTruth->stdhep[0];
    stdHepHigh = ntpMCTruth->stdhep[1];
    
    //loop over the stdhep entries for the current mc event
    for(Int_t j = stdHepLow; j < stdHepHigh+1; ++j){
      ntpMCStdHep = dynamic_cast<NtpMCStdHep *>(stdArray->At(j));
      
      //add the code for what you want to fill here 
    }                   
    
  }//end loop over mc events
  
  
  return;
}

//----------------------------------------------------------------------
void CondensedNtpModule::ResetTreeVariables()
{
    
  fEventInfo->Reset();
  fShowerInfo->Reset();
  fTrackInfo->Reset();
  fTruthInfo->Reset();
  
  return;
}

---