CondensedNtpModuleNC.cxx

Go to the documentation of this file.

//$Id: CondensedNtpModuleNC.cxx,v 1.23 2007/09/05 15:28:24 brebel Exp $
//
//CondensedNtpModuleNC.cxx
//
//Module for making analysis trees from SR ntuples
//
//B. Rebel 10/2004

#include "AnalysisNtuples/Module/CondensedNtpModuleNC.h"

//User/LOON Includes
#include "AnalysisNtuples/ANtpEventInfoNC.h"
#include "AnalysisNtuples/ANtpShowerInfoNC.h"
#include "AnalysisNtuples/ANtpTrackInfoNC.h"
#include "AnalysisNtuples/ANtpTruthInfoBeam.h"
#include "AnalysisNtuples/ANtpBeamInfo.h"
#include "AnalysisNtuples/ANtpHeaderInfo.h"
#include "AnalysisNtuples/Module/ANtpInfoObjectFillerNC.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 "Validity/VldContext.h" 
#include "StandardNtuple/NtpStRecord.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 "Record/RecCandHeader.h"
#include "RecoBase/ArrivalTime.h"
#include "RecoBase/LinearFit.h"
#include "Conventions/PlaneView.h"
#include "Conventions/Detector.h"
#include "Conventions/SimFlag.h"
#include "BeamDataNtuple/NtpBDLiteRecord.h"

#include "TFolder.h"
#include "TDirectory.h"
#include "TObjArray.h"
#include "TCanvas.h"
#include "TPaveText.h"
#include "TPad.h"
#include "TGraph.h"
#include "TH2.h"
#include "TH1.h"
#include "Riostream.h"

#include <cassert>
#include <algorithm>
#include <vector>

using namespace std;

CVSID("$Id: CondensedNtpModuleNC.cxx,v 1.23 2007/09/05 15:28:24 brebel Exp $");

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

//......................................................................
CondensedNtpModuleNC::CondensedNtpModuleNC() :
  fFileName("analysisNtuple.root"),
  fTreeName("analysisNtuple"),
  fNtpFile(0),
  fNtuple(0),
  fEventInfo(0),
  fHeaderInfo(0),
  fShowerInfo(0),
  fTrackInfo(0),
  fTruthInfo(0),
  fVHSPlanes(20),
  fVHSStrips(20)
{
  MSG("JobC", Msg::kDebug) << "CondensedNtpModuleNC::Constructor" << endl;
    
  fHeaderInfo = new ANtpHeaderInfo();
  fEventInfo = new ANtpEventInfoNC();
  fTrackInfo = new ANtpTrackInfoNC();
  fShowerInfo = new ANtpShowerInfoNC();
  fTruthInfo = new ANtpTruthInfoBeam();
  fBeamInfo = new ANtpBeamInfo();
  
  fInfoFiller = new ANtpInfoObjectFillerNC();
  
  fFailedDeMux = 0;
}

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

}

//----------------------------------------------------------------------
void CondensedNtpModuleNC::BeginJob()
{
  MSG("CondensedNtpModuleNC", Msg::kDebug) 
    << "in BeginJob" << endl;
  
  // save the current working directory
  //TDirectory* savedir  = gDirectory;  
  
  //create the new TFile for holding the electronics tree
  // this changes the value of gDirectory
  fNtpFile = new TFile(fFileName,"RECREATE");  
//   if(!fNtpFile->IsOpen() ) MSG("CondensedNtpModuleNC", Msg::kFatal) << "file not open " 
//                                                                  << fFileName << endl;

  // create TTree, these will attach themselves to the current 
  //working directory   
  fNtuple = new TTree(fTreeName, "Analysis Tree");
  
  //-----------------------------------------------------------------------
  //here is where you define your tree
  fNtuple->Branch("header.", "ANtpHeaderInfo", &fHeaderInfo, 64000, 2);
  fNtuple->Branch("event.", "ANtpEventInfoNC", &fEventInfo, 64000, 2);
  fNtuple->Branch("shower.", "ANtpShowerInfoNC", &fShowerInfo, 64000, 2);
  fNtuple->Branch("track.", "ANtpTrackInfoNC", &fTrackInfo, 64000, 2);
  fNtuple->Branch("beam.", "ANtpBeamInfo", &fBeamInfo, 64000, 2);
  if(fDataType==1) 
    fNtuple->Branch("truth.", "ANtpTruthInfoBeam", &fTruthInfo, 64000, 2);
  
  MSG("CondensedNtpModuleNC", Msg::kDebug) << "got branches" << endl;
  
  //-----------------------------------------------------------------------

  // Setup Kregg's image recognition filling
  fInfoFiller->InitializeVHSTraining(fVHSPlanes, fVHSStrips);
    
  // 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 CondensedNtpModuleNC::Ana(const MomNavigator *mom)
{
  // bool mc = false;
  
  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
  
  NtpMCRecord *mcRecord = 0;
  NtpTHRecord *thRecord = 0;
  
  mcRecord = dynamic_cast<NtpMCRecord *>(mom->GetFragment("NtpMCRecord"));
  thRecord = dynamic_cast<NtpTHRecord *>(mom->GetFragment("NtpTHRecord"));

  //get the record for the beam tree
  NtpBDLiteRecord *bdRecord = 0;
  //bdRecord = dynamic_cast<NtpBDLiteRecord *>(mom->GetFragment("NtpBDLiteRecord"));
      
  //instantiate a NtpHelper 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);

  const VldContext *vldc = stRecord ? stRecord->GetVldContext() : record->GetVldContext();
  VldTimeStamp timeStamp = vldc->GetTimeStamp();

  //fill the header information
  Detector::Detector_t det = vldc->GetDetector();

  //get a sim flag
  SimFlag::SimFlag_t dataType = vldc->GetSimFlag();

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

  MSG("CondensedNtpModuleNC", Msg::kDebug) << "snarl " << fHeaderInfo->snarl << endl;
  
//   MAXMSG("CondensedNtpModuleNC",Msg::kDebug,20)
//     << "detector = " << fHeaderInfo->detector << " " << det << endl;

  if(fHeaderInfo->year > 2004){
    if(bdRecord) {
      fInfoFiller->FillBeamInformation(bdRecord, fBeamInfo);
      MAXMSG("CondensedNtpModuleNC",Msg::kDebug,1) 
        << "Using NtpBDLiteRecord to fill beam info." << endl;
      MAXMSG("CondensedNtpModuleNC", Msg::kDebug,100) << timeStamp << " " 
                                                     << " " << fHeaderInfo->utc 
                                                     << " " << fBeamInfo->tortgt
                                                     << endl;
    }
    else {
      fInfoFiller->FillBeamInformation(timeStamp, det, dataType, fBeamInfo);
      MAXMSG("CondensedNtpModuleNC",Msg::kDebug,1)
        << "Calling db method to fill beam info." << endl;
    }
  }

  if(fHeaderInfo->passedDeMux == 0)
    ++fFailedDeMux;
  
  //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

  //for now set the primary shower criteria both to 0, that will pick out the first
  //shower in the reco array, which has been determined to be the one to use.  some 
  //other cuts still apply, but can be put in at the analysis level, not the ntuple
  //making level.  this was settled at the reco meeting on 2/9/06 - bjr
  ntpManipulator->SetPrimaryShowerCriteria(0,0); // nplanes, total pulse height
  
  //--------------------------------------------------------------------
  //here is where you start looping over the events in the ntuple to fill 
  //your analysis tree
  
  //handle the special case where you have a spill but no events
  if(fHeaderInfo->events < 1){
    //reset all the event, track and shower info objects to their default values
    //the beam and header values remain
    fEventInfo->Reset();
    fTrackInfo->Reset();
    fShowerInfo->Reset();
    fNtuple->Fill();
    return result;
  }

  //initialize the kNN stuff as it should only be done once per snarl
  fInfoFiller->InitializekNN(ntpManipulator);

  //loop over the events in the snarl
  for(Int_t i = 0; i < fHeaderInfo->events; ++i){
    
    MSG("CondensedNtpModuleNC", Msg::kDebug) << "on event " << i << endl;
 
    //reset the variables for the tree
    ResetTreeVariables();

    //if this isnt the first event in the snarl, then set the fHeaderInfo->newSnarl
    //to 0
    if(i > 0) fHeaderInfo->newSnarl = 0;

    if( fInfoFiller->FillInformation(i, ntpManipulator, fEventInfo, fTrackInfo,
                                     fShowerInfo, fTruthInfo) ){
    
      fNtuple->Fill();
      
    }//end if filling succeeded
  }//end loop over events for this snarl
  //-------------------------------------------------------------------
  
  return result;
}

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

//----------------------------------------------------------------------
void CondensedNtpModuleNC::EndJob() 
{
  
  MSG("CondensedNtpModuleNC", Msg::kInfo) << "start end job method " << fFailedDeMux 
                                          << " failed demuxing" << endl;        
  
  //cd to the "directory" of the file
  fNtpFile->cd();
  
  // only write file  -- this is enough and prevents duplication
  fNtpFile->Write();
  fNtpFile->Close();
  
  return;
}

//......................................................................
const Registry& CondensedNtpModuleNC::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("DataType",               1);

  r.Set("VHSPlanes",20);
  r.Set("VHSStrips",20);

  r.LockValues();

  //quiet the compiler
  itrue = ifalse;
  
  return r;
}

//......................................................................
void CondensedNtpModuleNC::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("DataType",               tmpi)) fDataType             =  tmpi;

  if (r.Get("VHSPlanes",              tmpi)) fVHSPlanes            =  tmpi;
  if (r.Get("VHSStrips",              tmpi)) fVHSStrips            =  tmpi;

  tmpb = 0;
  tmpd = 1.;

  return;
}

//......................................................................
void CondensedNtpModuleNC::ResetTreeVariables()
{

  fEventInfo->Reset();
  fTrackInfo->Reset();
  fShowerInfo->Reset();
  fTruthInfo->Reset();

  return;
}

---