CondensedNtpModuleNC.cxx

Go to the documentation of this file.

//$Id: CondensedNtpModuleNC.cxx,v 1.29 2009/11/24 16:50:39 rodriges Exp $
//
//CondensedNtpModuleNC.cxx
//
//Module for making analysis trees from SR ntuples
//
//B. Rebel 10/2004

#include "AnalysisNtuples/Module/CondensedNtpModuleNC.h"

#include "AnalysisNtuples/ANtpNueInfo.h"
#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/ANtpDefaultValue.h"
#include "AnalysisNtuples/ANtpEventInfoMRCC.h"

#include "AnalysisNtuples/Module/ANtpInfoObjectFillerNC.h"
#include "AnalysisNtuples/Module/ANtpInfoObjectFillerMRCC.h"
#include "AnalysisNtuples/Module/ANtpInfoObjectFillerNue.h"
#include "AnalysisNtuples/Module/ANtpRecoNtpManipulator.h"

#include "Conventions/Detector.h"
#include "Conventions/SimFlag.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro

#include "TFile.h"
#include "TTree.h"

#include <cassert>

class NtpSRRecord;
class NtpMCRecord;
class NtpMRRecord;
class NtpStRecord;
class NtpTHRecord;


CVSID("$Id: CondensedNtpModuleNC.cxx,v 1.29 2009/11/24 16:50:39 rodriges 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),
  fVHSPlanes(20),
  fVHSStrips(20)
{
  MSG("JobC", Msg::kDebug) << "CondensedNtpModuleNC::Constructor" << endl;

  fHeaderInfo = new ANtpHeaderInfo;
  fEventInfo = new ANtpEventInfoNC;
  fNueInfo = new ANtpNueInfo;
  fTrackInfo = new ANtpTrackInfoNC;
  fShowerInfo = new ANtpShowerInfoNC;
  fTruthInfo = new ANtpTruthInfoBeam;
  fBeamInfo = new ANtpBeamInfo;
  fEventInfoMRCC = new ANtpEventInfoMRCC;

  fInfoFiller = new ANtpInfoObjectFillerNC;
  fInfoFillerNue = new ANtpInfoObjectFillerNue;
  fInfoFillerMRCC = new ANtpInfoObjectFillerMRCC;

  fFailedDeMux = 0;
}

//----------------------------------------------------------------------
CondensedNtpModuleNC::~CondensedNtpModuleNC()
{

  MSG("JobC", Msg::kDebug) << "CondensedNtpModuleNC::Destructor" << endl;

  if(fHeaderInfo) delete fHeaderInfo;
  if(fEventInfo) delete fEventInfo;
  if(fEventInfoMRCC) delete fEventInfoMRCC;
  if(fShowerInfo) delete fShowerInfo;
  if(fTrackInfo) delete fTrackInfo;
  if(fTruthInfo) delete fTruthInfo;
  if(fBeamInfo) delete fBeamInfo;
  if(fNueInfo) delete fNueInfo;
}

//----------------------------------------------------------------------
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(fMakeNueBranch)fNtuple->Branch("nue.", "ANtpNueInfo", &fNueInfo, 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)
{
  JobCResult result(JobCResult::kPassed);

  //get the records from MOM
  assert(mom);

  NtpSRRecord* record     = (NtpSRRecord*)(mom->GetFragment("NtpSRRecord"));
  NtpMRRecord* mrRecord   = (NtpMRRecord*)(mom->GetFragment("NtpMRRecord"));

  NtpStRecord* stRecord=0;
  NtpStRecord* stRecordOrig=0;
  // Only want to create the MRCC output branch if the input file had MRCC in
  // This information isn't available to us in BeginJob, I can't think of a
  // better way than this.
  if(mrRecord){

    // This is a little fiddly: for MRCC processing, we need two
    // NtpStRecords; the one that got produced when we processed with
    // MRCC, and the one for the same snarl as processed without
    // MRCC. The former has name "MuonRemoved", while the latter has
    // name "Primary"
    stRecordOrig = (NtpStRecord*)(mom->GetFragment("NtpStRecord",
                                                   "Primary"));
    stRecord     = (NtpStRecord*)(mom->GetFragment("NtpStRecord",
                                                   "MuonRemoved"));
    
    if(!stRecordOrig){
      // Would really like to exit if this happens, to avoid the case
      // where only one file is given (although this check is done in
      // the macro). However, concatenation is not our friend, and
      // sometimes the concatenated sntp file gets chopped before the
      // mrnt file, so it doesn't have all the same subruns in it. In
      // this case, just ignore this snarl.
      MAXMSG("CondensedNtpModuleNC", Msg::kWarning, 10)
        << "No NtpStRecord found. Probably the sntp file is shorter than the mrnt."
        << "Ignoring this snarl" << endl;
      return JobCResult::kFailed;
    }
    static bool once = true;
    if(once){
      fNtuple->Branch("mrcc.", "ANtpEventInfoMRCC", &fEventInfoMRCC, 64000, 2);
      once = false;
    }
  }
  else{
    // If we're not running MRCC, the Primary NtpStRecord is the one to use
    stRecord = (NtpStRecord*)(mom->GetFragment("NtpStRecord",
                                               "Primary"));
  }

  if(!record && !stRecord){
    MSG("CondensedNtpModuleNC", Msg::kWarning) << "Could not get NtpSR or NtpSt "
                                                << "Record from MOM" << endl;
    result.SetWarning().SetFailed();
    return result;
  }//end if no NtpSRRecord

  NtpMCRecord* mcRecord = (NtpMCRecord*)(mom->GetFragment("NtpMCRecord"));
  NtpTHRecord* thRecord = (NtpTHRecord*)(mom->GetFragment("NtpTHRecord"));

  //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);

  if(fHeaderInfo->snarl%500==0)
    MSG("CondensedNtpModuleNC", Msg::kInfo) << "snarl " << fHeaderInfo->snarl << endl;

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

  if(fHeaderInfo->year > 2004){
    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();
    fEventInfoMRCC->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(fMakeNueBranch)fInfoFillerNue->Analyze(i, ntpManipulator->GetNtpStRecord(), fNueInfo);

    // HACK HACK HACK: 
    //
    // With MRCC files, this line will fill fTruthInfo
    // with the truth found from the mrnt file. Sadly, due to the
    // processing snafu documented in docdb 6599, the truth
    // information in mrnts isn't filled, so we'll later on fill the
    // truth from fInfoFillerMRCC::FillEventInformation
    if( fInfoFiller->FillInformation(i, ntpManipulator, fEventInfo, fTrackInfo,
                                     fShowerInfo, fTruthInfo) ){
      // The index field should actually be the index in this snarl. If this
      // isn't true then I don't know how the indexing works.
      assert(fEventInfo->index == i);
      // Fill the MRCC branch if this appears to be an MRCC file. Will
      // also fill the truth information, because of the problem
      // explained a couple of lines above
      if(mrRecord){
        fInfoFillerMRCC->FillEventInformation(i, mrRecord, stRecordOrig, 
                                              fEventInfoMRCC, 
                                              fDataType==1 ? fTruthInfo : 0);
        // Make sure it did the right thing
        assert(fEventInfoMRCC->bestEvent == fEventInfo->index || 
               fEventInfoMRCC->bestEvent == ANtpDefaultValue::kInt);
      } // end if mrRecord

      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.Set("MakeNueBranch",          false);

  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.
  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;

  if (r.Get("MakeNueBranch",          tmpb)) fMakeNueBranch        =  tmpb;

}

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

  fEventInfo->Reset();
  fTrackInfo->Reset();
  fShowerInfo->Reset();
  fTruthInfo->Reset();
  if(fMakeNueBranch)fNueInfo->Reset();
  fEventInfoMRCC->Reset();

  return;
}

---