NtpTHModule.cxx

Go to the documentation of this file.

//
// NtpTHModule.cxx
//
// A JobControl Module for filling an NtpTHRecord
//
#include <iostream>
using namespace std;
#include <cassert>

#include "TClonesArray.h"

#include "MessageService/MsgService.h"
#include "JobControl/JobCModuleRegistry.h"
#include "JobControl/JobCommand.h"
#include "MinosObjectMap/MomNavigator.h"
#include "Record/SimSnarlRecord.h"
#include "Record/SimSnarlHeader.h"
#include "CandData/CandRecord.h"
#include "RecoBase/CandStripListHandle.h"
#include "RecoBase/CandStripHandle.h"
#include "RecoBase/CandSliceListHandle.h"
#include "RecoBase/CandSliceHandle.h"
#include "RecoBase/CandShowerListHandle.h"
#include "RecoBase/CandShowerHandle.h"
#include "RecoBase/CandTrackListHandle.h"
#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandFitTrackListHandle.h"
#include "RecoBase/CandFitTrackHandle.h"
#include "RecoBase/CandEventListHandle.h"
#include "RecoBase/CandEventHandle.h"
#include "DataUtil/TruthHelper.h"
#include "DataUtil/Truthifier.h"
#include "TruthHelperNtuple/Module/NtpTHModule.h"
#include "TruthHelperNtuple/NtpTHRecord.h"
#include "TruthHelperNtuple/NtpTHStrip.h"
#include "TruthHelperNtuple/NtpTHSlice.h"
#include "TruthHelperNtuple/NtpTHShower.h"
#include "TruthHelperNtuple/NtpTHTrack.h"
#include "TruthHelperNtuple/NtpTHEvent.h"
#include "StandardNtuple/NtpStRecord.h"
#include "Digitization/DigiSignal.h"

ClassImp(NtpTHModule)

//   Definition of static data members
//   *********************************

CVSID("$Id: NtpTHModule.cxx,v 1.17 2008/11/05 17:45:13 rodriges Exp $");
JOBMODULE(NtpTHModule, "NtpTHModule",
         "A module for filling truth helper ntuple records.");

//-------------------------------------------------------------------------------------
namespace NtpTH
{
   struct Fraction
   {
      Fraction() : stdhep(0), adc(0.0), frac(0.0) {}
      int stdhep;
      double adc;
      double frac;
   };

   typedef std::vector<Fraction> FVec;
   typedef FVec::iterator FIter;

   struct ByFrac : public std::binary_function<Fraction, Fraction, bool>
   {
      bool operator()(const Fraction &lhs, const Fraction &rhs)
         {
            return (lhs.frac < rhs.frac);
         }
   };

   struct ById : public std::binary_function<Fraction, Fraction, bool>
   {
      bool operator()(const Fraction &lhs, const Fraction &rhs)
         {
            return (lhs.stdhep < rhs.stdhep);
         }
   };
}

// Definition of methods (alphabetical order)
// ***************************************************


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

const Registry& NtpTHModule::DefaultConfig() const {
  //
  // Purpose: Method to return default configuration.
  // 
  // Arguments: none.
  //
  // Return: Registry containing default configuration
  //

  MSG("NtpTH",Msg::kDebug) << "NtpTHModule::DefaultConfig" << endl;

  static Registry r; 
  std::string name = this->JobCModule::GetName();
  name += ".config.default";
  r.SetName(name.c_str());

  r.UnLockValues();
  r.Set("WriteStrip",1);
  r.Set("UseStandard",0);
  r.Set("CandRecordName","PrimaryCandidateRecord");
  r.Set("SimSnarlRecordName","");
  r.Set("RecordName","Primary");
  r.Set("RecordTitle","Created by NtpTHModule");
  r.LockValues();

  return r;
}

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

void NtpTHModule::Config(const Registry& r) {
  //
  // Purpose: Configure the module given a registry.
  //
  // Arguments: Registry to use to configure the module.
  //
  // Return: none.
  //

  MSG("NtpTH",Msg::kDebug) << "NtpTHModule::Config" << endl;
  
  Int_t tmpi;
  //  Double_t tmpd;
  const Char_t* tmps;
  if (r.Get("WriteStrip",tmpi)) fWriteStrip = tmpi;
  if (r.Get("UseStandard",tmpi)) fUseStandard = tmpi;
  if ( r.Get("CandRecordName",tmps) ) fCandRecordName = tmps;
  if ( r.Get("SimSnarlRecordName",tmps) ) fSimSnarlRecordName = tmps;
  if ( r.Get("RecordName", tmps) )  fRecordName = tmps;
  if ( r.Get("RecordTitle", tmps) ) fRecordTitle = tmps;

}

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

JobCResult NtpTHModule::Reco(MomNavigator *mom) {
  //
  //  Purpose:  Create and fill ntuple record.
  //
  //  Arguments: mom.
  //  
  //  Return: status.
  // 

  JobCResult result(JobCResult::kPassed);  
  MSG("NtpTH",Msg::kDebug) << "NtpTHModule::Reco" << endl;

  // Check that mom exists.
  assert(mom);

  SimSnarlRecord* simrec = dynamic_cast<SimSnarlRecord*>
              (mom->GetFragment("SimSnarlRecord",fSimSnarlRecordName.c_str()));
  CandRecord* cndrec     = dynamic_cast<CandRecord*>
              (mom->GetFragment("CandRecord",fCandRecordName.c_str()));
  if ( !simrec ) {
    MSG("NtpTH",Msg::kWarning) << " No SimSnarlRecord of name " 
                               << fSimSnarlRecordName.c_str() << endl;
    result.SetWarning().SetFailed();
    return result; 
  }
 
  const SimSnarlHeader* simhdr = simrec->GetSimSnarlHeader();

  TClonesArray* ntpthstp_arrptr = 0;
  TClonesArray* ntpthslc_arrptr = 0;
  TClonesArray* ntpthshw_arrptr = 0;
  TClonesArray* ntpthtrk_arrptr = 0;
  TClonesArray* ntpthevt_arrptr = 0;
  
  NtpStRecord* ntpstrec = 0;
  NtpTHRecord* ntpthrec = 0;
  if ( fUseStandard ) {
    ntpstrec = dynamic_cast<NtpStRecord*>(mom->GetFragment("NtpStRecord",
                                                  fRecordName.c_str()));
    if (!ntpstrec) {
      MSG("NtpSt",Msg::kWarning) << "No NtpStRecord in Mom and UseStandard."
                      << "\nMust call NtpStModule::Get() first." << endl;
      result.SetWarning().SetFailed();
      return result;
    }
    ntpthstp_arrptr = ntpstrec->thstp;
    ntpthslc_arrptr = ntpstrec->thslc;
    ntpthshw_arrptr = ntpstrec->thshw;
    ntpthtrk_arrptr = ntpstrec->thtrk;
    ntpthevt_arrptr = ntpstrec->thevt;
  }
  else {
    // Create a NtpTHRecord 
    RecCandHeader ntphdr(simhdr->GetVldContext(),simhdr->GetRun(),
           simhdr->GetSubRun(),simhdr->GetRunType(),simhdr->GetErrorCode(),
           simhdr->GetSnarl(),simhdr->GetTrigSrc(),simhdr->GetTimeFrame(),
           simhdr->GetRemoteSpillType(),-1);
    ntpthrec = new NtpTHRecord(ntphdr);

    RecJobHistory& jobhist
      = const_cast<RecJobHistory&>(ntpthrec->GetJobHistory());
    jobhist.Append(simrec->GetJobHistory());
    if ( cndrec ) jobhist.Append(cndrec->GetJobHistory());
    jobhist.CreateJobRecord(RecJobHistory::kNtpTH);
    
    ntpthrec->SetName(fRecordName.c_str());
    ntpthrec->SetTitle(fRecordTitle.c_str());
    ntpthstp_arrptr = ntpthrec->thstp;
    ntpthslc_arrptr = ntpthrec->thslc;
    ntpthshw_arrptr = ntpthrec->thshw;
    ntpthtrk_arrptr = ntpthrec->thtrk;
    ntpthevt_arrptr = ntpthrec->thevt;
  }
  
  if ( cndrec ) {
    // Only if CandRecord (event triggered)
    TClonesArray& ntpthstp_array = *(ntpthstp_arrptr);
    TClonesArray& ntpthslc_array = *(ntpthslc_arrptr);
    TClonesArray& ntpthshw_array = *(ntpthshw_arrptr);
    TClonesArray& ntpthtrk_array = *(ntpthtrk_arrptr);
    TClonesArray& ntpthevt_array = *(ntpthevt_arrptr);

    TruthHelper* helper = new TruthHelper(mom);
    if ( fWriteStrip ) {
      const Truthifier& truth = Truthifier::Instance(mom);
      this -> FillNtpTHStrip(ntpthstp_array,cndrec,truth,helper);
    }
    this -> FillNtpTHSlice(ntpthslc_array,cndrec,helper);
    this -> FillNtpTHShower(ntpthshw_array,cndrec,helper);
    this -> FillNtpTHTrack(ntpthtrk_array,cndrec,helper);
    this -> FillNtpTHEvent(ntpthevt_array,cndrec,helper);
    if ( helper ) delete helper; helper = 0;
  }
  
  // pass record to mom to own
  if ( !fUseStandard ) mom -> AdoptFragment(ntpthrec);  

  return result;

}

void NtpTHModule::FillNtpTHStrip(TClonesArray& ntpstriparray, 
                      const CandRecord* cndrec,
                      const Truthifier& truth, TruthHelper* helper) {
  //
  //  Purpose:  Private method used to fill strip portion of ntuple record.
  //
  //  Arguments: TClonesArray of NtpTHStrip, CandRecord, Truthifier 
  //             & TruthHelper
  //  
  //  Return: status.
  // 

  MSG("NtpTH",Msg::kDebug) << "NtpTHModule::FillNtpTHStrip" << endl;

  const CandStripListHandle* striplisthandle 
    = dynamic_cast < const CandStripListHandle* >
      (cndrec -> FindCandHandle("CandStripListHandle"));
  if ( !striplisthandle ) return; // all done

  Int_t nstrip = 0;
  TIter stripItr(striplisthandle->GetDaughterIterator());
  while(CandStripHandle* strip = dynamic_cast<CandStripHandle*>(stripItr())){
    NtpTHStrip* ntpstrip = new(ntpstriparray[nstrip++])NtpTHStrip();
    Int_t neuId = helper -> GetStripNeuIndex(strip);
    ntpstrip -> index = nstrip - 1;
    ntpstrip -> neumc = helper -> GetNeuKinIndex(neuId);
    ntpstrip -> nneu = (Short_t)helper -> StripPurity(strip);

    double strip_adc = 0.0;
    NtpTH::FVec svec;

    TIter digitItr(strip -> GetDaughterIterator());
    while (CandDigitHandle* digit=dynamic_cast<CandDigitHandle*>(digitItr()))
    {
       Int_t truth_flags = (Int_t)truth.GetSignalTruthFlags(*digit);
       ntpstrip -> sigflg |= truth_flags; // OR of all digits on strip
       
       
       const double adc = digit -> GetCharge(CalDigitType::kNone);
       
       strip_adc += adc;
       
       // get DigiSignal that matches this digit
       const DigiSignal *signal = truth.GetSignal(*digit);
       if(!signal)
       {
          continue;
       }
       
       // vector of hit fractions for this digit
       NtpTH::FVec dvec;
       
       // sum of all DigiScintHit weights
       double weight_sum = 0;
       
       for(unsigned int ihit = 0; ihit < signal -> GetNumberOfHits(); ++ihit)
       {
          const DigiScintHit* dhit = signal -> GetHit(ihit);
          if(!dhit)
          {
             continue;
          }
          
          const double weight = signal -> GetHitWeight(ihit);
          if(!(weight > 0.0))
          {
             continue;
          }
          
          NtpTH::Fraction frac;
          frac.stdhep = dhit -> TrackId();
          frac.adc    = adc;
          frac.frac   = weight;
          
          dvec.push_back(frac);
          
          weight_sum += weight;
       }
      
       // no particle weights present - this strip is probably pure noise
       if(!(weight_sum > 0.0))
       {
          continue;
       }

       // get total weight of all DigiScintHits in DigiSignal
       // prior to 2006.09.21 candidates had incorrect total weight
       // for such candidate files use weight_sum as total weight
       // this slightly underestimates total weight
       const double total_weight = signal -> GetTotalHitWeight();       
       if(total_weight > weight_sum)
       {
          weight_sum = total_weight;
       }
       
       //normalize fractions to total weight
       for(NtpTH::FIter fit = dvec.begin(); fit != dvec.end(); ++fit)
       {
          fit -> frac /= weight_sum;
       }
       
       // add digit fractions to all strip digit fractions
       svec.insert(svec.begin(), dvec.begin(), dvec.end());
       
    } // end of CandDigit loop

    // strange case - no charge in this CandStrip
    if(!(strip_adc > 0.0))
    {
       continue;
    }

    // sort vector by stdhep index
    std::sort(svec.begin(), svec.end(), NtpTH::ById());

    // vector for total reconstructed CandStrip pulse height fractions
    NtpTH::FVec frac_vec;
    
    // merge together fractions with the same stdhep index
    // that are attributed to different CandDigits
    NtpTH::FIter frac_it = svec.begin();
    while(frac_it != svec.end())
    {
       // find vector range with equal values of stdhep index
       const pair<NtpTH::FIter, NtpTH::FIter> prange =
          std::equal_range(frac_it, svec.end(), *frac_it, NtpTH::ById());
       
       // should never happen but check anyway just in case
       if(prange.first == prange.second)
       {
          break;
       }
       
       NtpTH::Fraction frac_all;
       for(NtpTH::FIter fit = prange.first; fit != prange.second; ++fit)
       {
          if(fit == prange.first)
          {
             frac_all.stdhep = fit -> stdhep;
          }
          else if(frac_all.stdhep != fit -> stdhep)
          {
             MSG("NtpTH",Msg::kError) << "Mismatched stdhep index..." << endl;
             continue;
          }
          
          frac_all.adc += (fit -> frac) * (fit -> adc);
       }
       
       frac_all.frac = frac_all.adc/strip_adc;
       
       frac_vec.push_back(frac_all);
       
       frac_it = prange.second;
    }

    // sort vector by pulse height fraction
    std::sort(frac_vec.begin(), frac_vec.end(), NtpTH::ByFrac());
    
    // record stdhep particles with highest fractions
    for(vector<NtpTH::Fraction>::reverse_iterator rit = frac_vec.rbegin();
        rit != frac_vec.rend(); ++rit)
    {
       const unsigned int idist = std::distance(frac_vec.rbegin(), rit);
       
       if(idist > 2)
       {
          break;
       }
       
       ntpstrip -> stdhep[idist] = rit -> stdhep;
       ntpstrip -> phfrac[idist] = static_cast<float> (rit -> frac);
       
       MSG("NtpTH",Msg::kVerbose) << idist
                                  << ": stdhep = " << rit -> stdhep
                                    << ", frac = " << rit -> frac << endl;
    }

    MSG("NtpTH",Msg::kVerbose) << *(ntpstrip) << endl;
    
  } // of CandStripList loop

  return;
  
}

void NtpTHModule::FillNtpTHSlice(TClonesArray& ntpslicearray,
                              const CandRecord* cndrec, TruthHelper* helper ) {
  //
  //  Purpose:  Private method used to fill slice portion of ntuple record.
  //
  //  Arguments: TClonesArray of NtpTHSlice, CandRecord & TruthHelper
  //  
  //  Return: status.
  // 

  MSG("NtpTH",Msg::kDebug) << "NtpTHModule::FillNtpTHSlice" << endl;

  const CandSliceListHandle* slicelisthandle 
    = dynamic_cast < const CandSliceListHandle* >
      (cndrec -> FindCandHandle("CandSliceListHandle"));
  if ( !slicelisthandle ) return; // all done

  Int_t nslice = 0;
  TIter sliceItr(slicelisthandle->GetDaughterIterator());
  while(CandSliceHandle* slc = dynamic_cast<CandSliceHandle*>(sliceItr())){
    NtpTHSlice* ntpslice = new(ntpslicearray[nslice++])NtpTHSlice();
    Int_t neuId = helper -> GetBestSliceNeuMatch(*slc);
    ntpslice -> index = nslice - 1;
    ntpslice -> neumc = helper -> GetNeuKinIndex(neuId);
    ntpslice -> neustdhep = neuId;
    ntpslice -> nneu = helper -> NumNeu(*slc);
    ntpslice -> purity = helper -> SlicePurity(*slc);
    ntpslice -> secondpurity = helper -> secondNEU(*slc);
    ntpslice -> complete = helper -> SliceCompleteness(*slc);
    MSG("NtpTH",Msg::kDebug) << *(ntpslice) << endl;
  }

  return;

}

void NtpTHModule::FillNtpTHShower(TClonesArray& ntpshowerarray,
  const CandRecord* cndrec, TruthHelper* helper ) {
  //
  //  Purpose:  Private method used to fill shower portion of ntuple record.
  //
  //  Arguments: NtpTHRecord, CandRecord and TruthHelper
  //  
  //  Return: status.
  // 

  MSG("NtpTH",Msg::kDebug) << "NtpTHModule::FillNtpTHShower" << endl;

  const CandShowerListHandle* showerlisthandle 
    = dynamic_cast < const CandShowerListHandle* >
      (cndrec -> FindCandHandle("CandShowerListHandle"));
  if ( !showerlisthandle ) return; // all done

  Int_t nshower = 0;
  TIter showerItr(showerlisthandle->GetDaughterIterator());
  while(CandShowerHandle* shw = dynamic_cast<CandShowerHandle*>(showerItr())){
    NtpTHShower* ntpshower = new(ntpshowerarray[nshower++])NtpTHShower();
    Int_t neuId = helper -> GetBestShowerNeuMatch(*shw);
    ntpshower -> index    = nshower - 1;
    ntpshower -> neumc  = helper -> GetNeuKinIndex(neuId);
    ntpshower -> neustdhep    = neuId;
    ntpshower -> purity   = helper -> ShowerPurity(*shw,neuId);
    ntpshower -> completeall = helper -> ShowerCompleteness(*shw);
    ntpshower -> completeallnopecut = helper -> ShowerCompletenessAllStrips(*shw);

    CandSliceHandle* slch = const_cast<CandSliceHandle*>(shw -> GetCandSlice());
    if ( slch ) {
      ntpshower -> completeslc = helper -> ShowerCompleteness(*shw,*slch);
      ntpshower -> completeslcnopecut = helper -> ShowerCompletenessAllStrips(*shw,*slch);
    }

    MSG("NtpTH",Msg::kDebug) << *(ntpshower) << endl;
  }

  return;

}

void NtpTHModule::FillNtpTHTrack(TClonesArray& ntptrackarray,
  const CandRecord* cndrec, TruthHelper* helper ) {
  //
  //  Purpose:  Private method used to fill track portion of ntuple record.
  //
  //  Arguments: NtpTHRecord, CandRecord and TruthHelper
  //  
  //  Return: status.
  // 

  MSG("NtpTH",Msg::kDebug) << "NtpTHModule::FillNtpTHTrack" << endl;

  const CandTrackListHandle* tracklisthandle 
    = dynamic_cast < const CandTrackListHandle* >
      (cndrec -> FindCandHandle("CandTrackListHandle"));
  if ( !tracklisthandle ) return; // all done

  Int_t ntrack = 0;
  TIter trackItr(tracklisthandle->GetDaughterIterator());
  while(CandTrackHandle* trk = dynamic_cast<CandTrackHandle*>(trackItr())){
    if(trk->GetNDaughters()==0 &&
      trk->InheritsFrom("CandFitTrackHandle") &&
       dynamic_cast<CandFitTrackHandle*>(trk)->GetFinderTrack())
      trk=dynamic_cast<CandFitTrackHandle*>(trk)->GetFinderTrack();

    NtpTHTrack* ntptrack = new(ntptrackarray[ntrack++])NtpTHTrack();
    Int_t bestTrack = helper -> GetBestTrackIdMatch(*trk);
    ntptrack -> index = ntrack - 1;
    ntptrack -> neumc  = helper -> GetNeuKinIndex(bestTrack);
    ntptrack -> trkstdhep = bestTrack;
    ntptrack -> neustdhep = helper -> GetNeuId(bestTrack);
    ntptrack -> completeall = helper -> TrackCompleteness(*trk);
    ntptrack -> purity = helper -> TrackPurity(*trk,bestTrack);
    CandSliceHandle* slch = const_cast<CandSliceHandle*>(trk -> GetCandSlice());
    if ( slch ) {
      ntptrack -> completeslc = helper -> TrackCompleteness(*trk,*slch);
    }  
    MSG("NtpTH",Msg::kDebug) << *(ntptrack) << endl;
  }

  return;

}

void NtpTHModule::FillNtpTHEvent(TClonesArray& ntpeventarray,
  const CandRecord* cndrec, TruthHelper* helper ) {
  //
  //  Purpose:  Private method used to fill event portion of ntuple record.
  //
  //  Arguments: TClonesArray, CandRecord and TruthHelper
  //  
  //  Return: status.
  // 

  MSG("NtpTH",Msg::kDebug) << "NtpTHModule::FillNtpTHEvent" << endl;

  const CandEventListHandle* eventlisthandle 
    = dynamic_cast < const CandEventListHandle* >
      (cndrec -> FindCandHandle("CandEventListHandle"));
  if ( !eventlisthandle ) return; // all done

  Int_t nevent = 0;
  TIter eventItr(eventlisthandle->GetDaughterIterator());

  while(CandEventHandle* evt = dynamic_cast<CandEventHandle*>(eventItr())){
    NtpTHEvent* ntpevent = new(ntpeventarray[nevent++])NtpTHEvent();
    Int_t bestEvent = helper -> GetBestEventNeuMatch(*evt);
    ntpevent -> index = nevent - 1;
    ntpevent -> neumc = helper -> GetNeuKinIndex(bestEvent);
    ntpevent -> neustdhep = bestEvent;
    ntpevent -> purity = helper -> EventPurity(*evt,bestEvent);
    ntpevent -> completeall = helper -> EventCompleteness(*evt);
    ntpevent -> completeallnopecut = helper -> EventCompletenessAllStrips(*evt);
    CandSliceHandle* slch = const_cast<CandSliceHandle*>(evt -> GetCandSlice());
    if ( slch ) {
      ntpevent -> completeslc = helper -> EventCompleteness(*evt,*slch);
      ntpevent -> completeslcnopecut = helper -> EventCompletenessAllStrips(*evt,*slch);
    }  
    MSG("NtpTH",Msg::kDebug) << *(ntpevent) << endl;
  }

  return;

}

---