Module/VtxModule.cxx

Go to the documentation of this file.

// $Id: VtxModule.cxx,v 1.2 2007/03/01 16:56:24 rhatcher Exp $
//
// FILL_IN: [Document your code!!]
//
// vahle@hep.ucl.ac.uk
#include <dirent.h>
#include "TChain.h"
#include "TClonesArray.h"
#include "TProfile2D.h"
#include "StandardNtuple/NtpStRecord.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
#include "Calibrator/CalMIPCalibration.h"
#include "MCNtuple/NtpMCRecord.h"
#include "MCNtuple/NtpMCTruth.h"
#include "TruthHelperNtuple/NtpTHRecord.h"
#include "CandNtupleSR/NtpSRRecord.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "NueAna/SntpHelpers.h"
#include "NueAna/NueModule.h"
#include "NueAna/NueRecord.h"
#include "NueAna/NueRecordAna.h"
#include "BeamData/ana/Summary/BeamSummary.h"
#include "NueAna/EventFilter.h"

#include "VertexFinder/Module/VtxModule.h"
#include "VertexFinder/Module/VtxRecordAna.h"


CVSID("$Id: VtxModule.cxx,v 1.2 2007/03/01 16:56:24 rhatcher Exp $");

#include "DatabaseInterface/DbiResultPtr.tpl"


JOBMODULE(VtxModule, "VtxModule",
          "");
//......................................................................

VtxModule::VtxModule():
   counter(0),
   passcounter(0),
   passcutcounter(0),
   failcounter(0),
   writecounter(0),
   foundmeu(false),
   SIGCORRMEU(1.),
   pot(0.),
   MEUPERGEV(25.66)
{}

//......................................................................

VtxModule::~VtxModule()
{
}

//......................................................................

JobCResult VtxModule::Reco(MomNavigator* mom)
{
//======================================================================
// Reads in sue-style ntuples from mom, calculates a bunch of useful
// nue quantites
//======================================================================
   MSG("VtxModule",Msg::kDebug)<<"In VtxModule::Reco"<<endl;

   if(counter%1000==0){      
      MSG("VtxModule",Msg::kInfo)<<"On entry "<<counter<<endl;
   }
   counter++;
   bool foundST=false;
      
   VldContext vc;
   //first ask mom for a NtpStRecord
   NtpStRecord *str=dynamic_cast<NtpStRecord *>(mom->GetFragment("NtpStRecord"));
   if(str){
     foundST=true;
     vc=str->GetHeader().GetVldContext();
   }

   //if we didn't find a NtpStRecord, it could be because we're reading an 
   //old tree, check for the other records
                                                
   if(!foundST){
      //got nothing useful from mom
      MSG("VtxModule",Msg::kWarning)<<"Got Nothing from mom"<<endl;
      failcounter++;
      return JobCResult::kFailed;
   }
   
   //UGLY!!!!
   //need to be able to compare single strip 
   //pulseheight between ND and FD, but we can't
   //in the ntuple.  Will ask the db for the
   //conversion from sigcor to MEU 
   //on the first event, and use throughout.  
   //Since sigcorrs are already corrected for all strip variations, 
   //it's ok to ask for 1 number.  
   //Should be ok for MC, maybe ok for 1 run of data at a time, but will
   //definitely be wrong if many weeks worth of data are strung together
   //we should revisit this and make it better in the future!!!!!!!
   if(!foundmeu){
     DbiResultPtr<CalMIPCalibration> dbp(vc);
     if(dbp.GetNumRows()>0){
       const CalMIPCalibration *m = dbp.GetRow(0);
       float mip=m->GetMIP(1.); 
       if(mip>0){
         SIGCORRMEU=1./mip;
         foundmeu=true;
       }
     }
   }

   if(foundST){
     //find number of events in the snarl
      int evtn=str->evthdr.nevent;
      MSG("VtxModule",Msg::kDebug)<<"Events in this snarl: "<<evtn<<endl;
      
      //find vldcontext to make the header
      VldContext vc = str->GetHeader().GetVldContext();
      
      if(evtn==0){
        //we found an SR, but nothing got reco'ed in it
        
        //make a header
        NueHeader h(vc);
      
      //fill header info
        h.SetSnarl(str->GetHeader().GetSnarl());
        h.SetRun(str->GetHeader().GetRun());
        h.SetSubRun(str->GetHeader().GetSubRun());
        h.SetEventNo(-1);
        h.SetEvents(0);
        h.SetTrackLength(0);
        h.SetFoundBits(false,true,false);
      
        //make an ana_nue object
        VtxRecord *vtxRec = new VtxRecord(h);
        vtxRec->SetName("VtxRecord");

        //make a nuerecordana object
        VtxRecordAna ana_vtx(*vtxRec);
        
        //fill the truth info
        //if it's data, we can also do beam mon analyze
        ana_vtx.FillTrue(0,str);

        //give the ana_nue object to mom, she'll write it to the file
         mom->AdoptFragment(vtxRec);
         writecounter++;
         failcounter++;
        
         return JobCResult::kPassed;
      }      
      
      //if events were reco'd, loop over events in a snarl
      for(int i=0;i<evtn;i++){
              
        //make a header
        NueHeader h(vc);
        
        //fill header info
        h.SetSnarl(str->GetHeader().GetSnarl());
        h.SetRun(str->GetHeader().GetRun());
        h.SetSubRun(str->GetHeader().GetSubRun());
        h.SetEventNo(i);
        h.SetEvents(evtn);
        h.SetFoundBits(true,true,true);
        
        MSG("VtxModule",Msg::kDebug)<<"Getting event"<<endl;
        NtpSREvent *event = SntpHelpers::GetEvent(i,str);
      
        //loop over tracks in this event, find longest
        int longtrack=0;
        for(int j=0;j<event->ntrack;j++){
          int tindex = SntpHelpers::GetTrackIndex(j,event);
          NtpSRTrack *track = SntpHelpers::GetTrack(tindex,str);
          if(longtrack<track->plane.n){
            longtrack = track->plane.n;
          }
        }
        h.SetTrackLength(longtrack);
        
        VtxRecord *vtxRec = new VtxRecord(h);
        vtxRec->SetName("VtxRecord");

        //make a nuerecordana object
        VtxRecordAna ana_vtx(*vtxRec);

        //configure BeamMonAna
         passcutcounter++;
          //analyze calculates everybodies
          //variables, then calls the Truth and Reco Object Fillers
          //configure MSTCalcAna routine
         ana_vtx.ansia.SetParams(SIGCORRMEU, MEUPERGEV);
         ana_vtx.antia.SetParams(SIGCORRMEU, MEUPERGEV);          
         ana_vtx.aneia.SetParams(SIGCORRMEU, MEUPERGEV);
          //analysize
         ana_vtx.FillTrue(i,str);
         ana_vtx.FillReco(i,str);                        
         ana_vtx.Analyze(i,str);
        
        //give the ana_nue object to mom, she'll write it to the file
        MSG("VtxModule",Msg::kDebug)<<"Giving Fragment to mom"<<endl;
        writecounter++;
        mom->AdoptFragment(vtxRec);
        MSG("VtxModule",Msg::kDebug)<<"Mom took it"<<endl;
      }
      MSG("VtxModule",Msg::kDebug)<<"Done with snarl"<<endl;
      passcounter++;
      return JobCResult::kPassed; // kNoDecision, kFailed, etc.       
   }

   return JobCResult::kFailed;
}

//......................................................................

const Registry& VtxModule::DefaultConfig() const
{
//======================================================================
// Supply the default configuration for the module
//======================================================================
   MSG("VtxModule",Msg::kDebug)<<"In VtxModule::DefaultConfig"<<endl;

  static Registry r; // Default configuration for module

  // Set name of config
  std::string name = this->GetName();
  name += ".config.default";
  r.SetName(name.c_str());

  // Set values in configuration
  r.UnLockValues();
  r.LockValues();

  return r;
}

//......................................................................

void VtxModule::Config(const Registry& r)
{
//======================================================================
// Configure the module given the Registry r
//======================================================================
  MSG("VtxModule",Msg::kDebug)<<"In VtxModule::Config"<<endl;

}

void VtxModule::BeginJob()
{
  MSG("VtxModule",Msg::kDebug)<<"In VtxModule::BeginJob"<<endl;
}

void VtxModule::EndJob()
{

   MSG("VtxModule",Msg::kInfo)<<"Counter "<<counter
                              <<" passcutcounter "<<passcutcounter
                              <<" passcounter "<<passcounter
                              <<" failcounter "<<failcounter
                              <<" writecounter "<<writecounter<<endl;
   MSG("VtxModule",Msg::kInfo)<<"Number of POT in this run: "<<pot<<endl;

}

---