NtpMCModule.cxx

Go to the documentation of this file.

//
// NtpMCModule.cxx
//
// A JobControl Module for filling an NtpMCRecord from data stored
// in the SimSnarlRecord.
//
#include <iostream>
#include <cassert>

using namespace std;

#include "TParticle.h"
#include "DataUtil/TruthHelper.h"

#include "MessageService/MsgService.h"
#include "JobControl/JobCModuleRegistry.h"
#include "JobControl/JobCommand.h"
#include "MinosObjectMap/MomNavigator.h"
#include "TClonesArray.h"
#include "MCNtuple/Module/NtpMCModule.h"
#include "MCNtuple/NtpMCRecord.h"
#include "MCNtuple/NtpMCTruth.h"
#include "MCNtuple/NtpMCStdHep.h"
#include "MCNtuple/NtpMCDigiScintHit.h"
#include "MCNtuple/NtpMCDetSimResult.h"
#include "MCNtuple/NtpMCPhotonResult.h"
#include "MCNtuple/NtpMCFluxInfo.h"
#include "MCNtuple/NtpMCFluxWgt.h"
#include "MCNtuple/NtpMCStdHepHit.h"
#include "StandardNtuple/NtpStRecord.h"
#include "Record/SimSnarlRecord.h"
#include "Record/SimSnarlHeader.h"
#include "Plex/PlexPlaneId.h"
#include "Plex/PlexStripEndId.h"
#include "REROOT_Classes/REROOT_NeuKin.h"
#include "REROOT_Classes/REROOT_FluxInfo.h"
#include "REROOT_Classes/REROOT_FluxWgt.h"
#include "Conventions/Munits.h"
#include "Conventions/Mphysical.h"
#include "Util/UtilMath.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "UgliGeometry/UgliStripHandle.h"
#include "Digitization/DigiScintHit.h"
#include "RawData/RawRecord.h"
#include "RawData/RawDigitDataBlock.h"
#include "RawData/RawDataBlock.h"
#include "RawData/RawDigit.h"
#include "DataUtil/Truthifier.h"
#include "PhotonTransport/PhotonEventResult.h"
#include "DetSim/SimEventResult.h"

ClassImp(NtpMCModule)

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

CVSID("$Id: NtpMCModule.cxx,v 1.41 2009/05/21 18:29:36 rhatcher Exp $");
JOBMODULE(NtpMCModule, "NtpMCModule",
         "A module for filling MC ntuple records.");


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

bool compareP4(Float_t* p4neu, TParticle* part)
{
  bool result = true;
  Float_t part_p4[4] = { part->Px(), part->Py(), part->Pz(), part->Energy() };
  for (int k=0; k<4; ++k) {
    // use float strictness scaled by 100 (vs default 10)
    // because for fava values are exchanged as ASCII messages
    // and more precision gets lost
    result &= UtilMath::sameValue(p4neu[k],part_p4[k],false,100.);
  }
  return result;
}

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

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

  MSG("NtpMC",Msg::kDebug) << 
    "NtpMCModule::DefaultConfig" << endl;

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

  r.UnLockValues();
  r.Set("WriteDigiHit",0);
  r.Set("UseStandard",0);
  r.Set("SimSnarlRecordName", "");
  r.Set("RecordName", "Primary");
  r.Set("RecordTitle", "Created by NtpMCModule");
  r.LockValues();

  return r;
}

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

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

  MSG("NtpMC",Msg::kDebug) << "NtpMCModule::Config" << endl;
  
  Int_t tmpi;
  //  Double_t tmpd;
  const Char_t* tmps;
  if ( r.Get("WriteDigiHit",tmpi) ) fWriteDigiHit = tmpi;
  if ( r.Get("UseStandard",tmpi) ) fUseStandard = tmpi;
  if ( r.Get("SimSnarlRecordName",tmps) ) fSimSnarlRecordName = tmps;
  if ( r.Get("RecordName", tmps) )  fRecordName = tmps;
  if ( r.Get("RecordTitle", tmps) ) fRecordTitle = tmps;

}

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

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

  JobCResult result(JobCResult::kPassed);  
  MSG("NtpMC",Msg::kDebug) << "NtpMCModule::Reco" << endl;

  // Check that mom exists.
  assert(mom);
  
  SimSnarlRecord* simrec = dynamic_cast<SimSnarlRecord*>
            (mom->GetFragment("SimSnarlRecord",fSimSnarlRecordName.c_str()));
  if ( !simrec ) {
    MSG("NtpMC",Msg::kWarning) << "No SimSnarlRecord of name "
                 << fSimSnarlRecordName.c_str() << " in Mom." << endl;
    result.SetWarning().SetFailed();
    return result;
  }

  const SimSnarlHeader* simhdr = simrec->GetSimSnarlHeader();

  if ( fStdHepRange ) delete [] fStdHepRange; fStdHepRange=0;
  
  TClonesArray* ntpmcarrptr = 0;
  TClonesArray* ntpstdheparrptr = 0;
  TClonesArray* ntpdigihitarrptr = 0;
  
  NtpStRecord* ntpstrec = 0;
  NtpMCRecord* ntpmcrec = 0;
  NtpMCSummary* ntpmcsummaryptr = 0;
  NtpMCPhotonResult* ntpmcphotonptr = 0;
  NtpMCDetSimResult* ntpmcdetsimptr = 0;
  
  if ( fUseStandard ) {
    ntpstrec = dynamic_cast<NtpStRecord*>(mom->GetFragment("NtpStRecord",
                                                fRecordName.c_str()));
    if ( !ntpstrec) {
      MSG("NtpMC",Msg::kWarning) << "No NtpStRecord in Mom of name "
                        << fRecordName.c_str() << " and UseStandard."
                        << "\nMust call NtpStModule::Get() first." << endl;
      result.SetWarning().SetFailed();
      return result;
    }
    ntpmcarrptr = ntpstrec->mc;
    ntpstdheparrptr = ntpstrec->stdhep;
    ntpdigihitarrptr = ntpstrec->digihit;
    ntpmcsummaryptr = &(ntpstrec->mchdr);
    ntpmcphotonptr = &(ntpstrec->photon);
    ntpmcdetsimptr = &(ntpstrec->detsim);
  }
  else {
    // Create a NtpMCRecord
    RecCandHeader ntphdr(simhdr->GetVldContext(),simhdr->GetRun(),
           simhdr->GetSubRun(),simhdr->GetRunType(),simhdr->GetErrorCode(),
           simhdr->GetSnarl(),simhdr->GetTrigSrc(),simhdr->GetTimeFrame(),
           simhdr->GetRemoteSpillType(),-1);
    ntpmcrec = new NtpMCRecord(ntphdr);
    ntpmcrec -> SetName(fRecordName.c_str());
    ntpmcrec -> SetTitle(fRecordTitle.c_str());

    RecJobHistory& jobhist 
          = const_cast<RecJobHistory&>(ntpmcrec->GetJobHistory());
    jobhist.Append(simrec->GetJobHistory());
    jobhist.CreateJobRecord(RecJobHistory::kNtpMC);
    
    ntpmcarrptr = ntpmcrec->mc;
    ntpstdheparrptr = ntpmcrec->stdhep;
    ntpdigihitarrptr = ntpmcrec->digihit;
    ntpmcsummaryptr = &(ntpmcrec->mchdr);
    ntpmcphotonptr = &(ntpmcrec->photon);
    ntpmcdetsimptr = &(ntpmcrec->detsim);
  }
  
  TClonesArray& ntpmctrutharray = *(ntpmcarrptr);
  TClonesArray& ntpstdheparray  = *(ntpstdheparrptr);
  NtpMCSummary& ntpmcsummary = *(ntpmcsummaryptr);
  NtpMCPhotonResult& ntpmcphoton = *(ntpmcphotonptr);
  NtpMCDetSimResult& ntpmcdetsim = *(ntpmcdetsimptr);
  
  ntpmcsummary.geninfo.time  = simhdr->GetMcGenTime();
  ntpmcsummary.geninfo.codename = simhdr->GetMcGenCodename();
  ntpmcsummary.geninfo.hostname = simhdr->GetMcGenHostname();

  this -> FillNtpMCTruth(ntpmctrutharray,ntpmcsummary,simrec,mom);

  this -> FillNtpMCStdHep(ntpstdheparray,ntpmcsummary,simrec);
  if ( fWriteDigiHit ) {
    TClonesArray& ntpdigihitarray = *(ntpdigihitarrptr);
    this -> FillNtpMCDigiScintHit(ntpdigihitarray,ntpmcsummary,simrec);
  }

  this -> FillNtpMCPhotonResult(ntpmcphoton,simrec);
  this -> FillNtpMCDetSimResult(ntpmcdetsim,simrec);

  // pass record to mom to own
  MSG("NtpMC",Msg::kDebug) << ntpmcsummary << endl;
  
  if ( !fUseStandard)  mom -> AdoptFragment(ntpmcrec);
  
  return result;

}

void NtpMCModule::FillNtpMCTruth(TClonesArray& ntpmctrutharray,
                                 NtpMCSummary& ntpmcsummary,
                                 const SimSnarlRecord* simrec,
                                 const MomNavigator* mom) {
  //
  //  Purpose:  Private method used to fill mc truth portion of ntuple record.
  //
  //  Arguments: NtpMCTruth, SimSnarlRecord and Mom for Truthifier (for filling
  //             ndigu(v)/tphu(v) per true event
  //  
  //  Return: status.
  // 

  MSG("NtpMC",Msg::kDebug) << "NtpMCModule::FillNtpMCTruth" << endl;

  if ( !simrec ) {
    MSG("NtpMC",Msg::kWarning) 
    << "NtpMCModule::FillNtpMCTruth called with null argument pointer." 
    << endl;
    return;
  }

  // SimSnarlRecord does *NOT* have REROOT_NeuVtx objects
  // (these rely on gMINFast to give correct z offset)
  // so we must retrieve the vertex info from the StdHep array

  const TClonesArray* neukinarray = dynamic_cast<const TClonesArray*>
                          (simrec->FindComponent("TClonesArray","NeuKinList"));
  const TClonesArray* stdheparray = dynamic_cast<const TClonesArray*>
                          (simrec->FindComponent("TClonesArray","StdHep"));
  const TClonesArray* finfoarray = dynamic_cast<const TClonesArray *>
    (simrec->FindComponent("TClonesArray","FluxInfoList"));
  const TClonesArray* fwgtarray = dynamic_cast<const TClonesArray *>
    (simrec->FindComponent("TClonesArray","FluxWgtList"));
  
  if ( neukinarray && neukinarray->GetEntriesFast() != 0 ) {
    REROOT_NeuKin* rneukin 
        = dynamic_cast<REROOT_NeuKin*>(neukinarray->At(0));
    if ( rneukin -> INu() != 0 ) {
      // Only complain if not a dummy NeuKin entry, such is generated
      // when user uses KINE card to set particle initial state
      if ( !finfoarray || finfoarray->GetEntriesFast() == 0 ) {
        MAXMSG("NtpMC",Msg::kWarning,1) << "FluxInfoList is empty.\n"
          << "mc.flux will not be properly filled!" << endl;
      }
      if ( !fwgtarray || fwgtarray->GetEntriesFast() == 0 ) {
        MAXMSG("NtpMC",Msg::kWarning,1) << "FluxWgtList is empty.\n"
          << "mc.fluxwgt will not be properly filled!" << endl;
      }
    }
  }


  if ( neukinarray && stdheparray ) {
    Int_t nneukin = neukinarray->GetLast()+1;
    ntpmcsummary.nmc = nneukin;
    // Set error bit if build failed
    if ( !BuildMCToStdHep(neukinarray,stdheparray) ) ntpmcsummary.error = 1;
    
    for ( Int_t ineu = 0; ineu < nneukin; ineu++ ) {
      NtpMCTruth* ntpmctruth = new(ntpmctrutharray[ineu])NtpMCTruth();
      ntpmctruth->index = ineu;
      if ( fStdHepRange ) {
        ntpmctruth->stdhep[0] = fStdHepRange[ineu*2+0];
        ntpmctruth->stdhep[1] = fStdHepRange[ineu*2+1];
      }
      
      REROOT_NeuKin* rneukin 
        = dynamic_cast<REROOT_NeuKin*>(neukinarray->At(ineu));
      if ( rneukin ) {
        ntpmctruth->inu        = rneukin->INu();
        ntpmctruth->inunoosc   = rneukin->INuNoOsc();
        ntpmctruth->itg        = rneukin->ITg();
        ntpmctruth->iboson     = rneukin->IBoson();
        ntpmctruth->iresonance = rneukin->IResonance();
        ntpmctruth->iaction    = rneukin->IAction();
        ntpmctruth->istruckq   = rneukin->IStruckQ();
        ntpmctruth->iflags     = rneukin->IFlags();
        ntpmctruth->a          = rneukin->A();
        ntpmctruth->z          = rneukin->Z();
        ntpmctruth->sigma      = rneukin->Sigma();
        ntpmctruth->sigmadiff  = rneukin->SigmaDiff();
        for (int i = 0; i < 4; i++) {
          ntpmctruth->p4neu[i] = rneukin->P4Neu()[i];
          ntpmctruth->p4neunoosc[i] = rneukin->P4NeuNoOsc()[i];
          ntpmctruth->p4tgt[i] = rneukin->P4Tgt()[i];
          ntpmctruth->p4shw[i] = rneukin->P4Shw()[i];
          ntpmctruth->p4mu1[i] = rneukin->P4Mu1()[i];
          ntpmctruth->p4mu2[i] = rneukin->P4Mu2()[i];
          ntpmctruth->p4el1[i] = rneukin->P4El1()[i];
          ntpmctruth->p4el2[i] = rneukin->P4El2()[i];
          ntpmctruth->p4tau[i] = rneukin->P4Tau()[i];
        }
        ntpmctruth->x  = rneukin->X();
        ntpmctruth->y  = rneukin->Y();
        ntpmctruth->q2 = rneukin->Q2();
        ntpmctruth->w2 = rneukin->W2();
        ntpmctruth->emfrac = rneukin->EMFrac();
      }
    
      RawRecord* rawrec = dynamic_cast<RawRecord*>
       (mom->GetFragment("RawRecord","","DaqSnarl"));

      if ( rawrec && ntpmctruth->stdhep[0] >= 0 ) {
        const Truthifier& truth = Truthifier::Instance(mom);
        TIter itr = rawrec -> GetRawBlockIter();
        while (RawDataBlock* rawblk=dynamic_cast<RawDataBlock*>(itr.Next())){
          RawDigitDataBlock* rddb = dynamic_cast<RawDigitDataBlock*>(rawblk);
          if ( rddb ) {
            TIter diter = rddb -> GetDatumIter();
            while (RawDigit* rawdigit=dynamic_cast<RawDigit*>(diter.Next())) {
              RawChannelId digitRCId = rawdigit->GetChannel();
              Double_t adc_offset = 0;
              if ( digitRCId.GetElecType() == ElecType::kQIE ) adc_offset=50;
              Bool_t digfound = false;
              Bool_t inu = false;
              Bool_t inv = false;
              if ( const DigiSignal* signal = truth.GetSignal(rawdigit) ) {
                if ( signal ) {
                  for(UInt_t ihit=0; ihit<signal->GetNumberOfHits(); ihit++) {
                    const DigiScintHit* hit = signal->GetHit(ihit);
                    if ( hit ) {
                      Int_t trkid = abs(hit->TrackId());
                      if ( trkid >= ntpmctruth->stdhep[0] && 
                           trkid <= ntpmctruth->stdhep[1] ) {
                        digfound = true;
                        PlexStripEndId seid = hit->StripEndId();
                        if ( seid.GetPlaneView()==PlaneView::kU )inu = true;
                        if ( seid.GetPlaneView()==PlaneView::kV )inv = true;
                      }
                    }
                  }
                }
              }
              if ( digfound ) {
                if ( inu ) {
                  ntpmctruth->ndigu += 1;
                  ntpmctruth->tphu += (rawdigit->GetADC()-adc_offset);
                }
                else if ( inv ) {
                  ntpmctruth->ndigv += 1;
                  ntpmctruth->tphv += (rawdigit->GetADC()-adc_offset);
                }
              }
            }
          }
        }         
      }

      TruthHelper th;  // used to access methos IsNeutrino & IsDocStatus
      // look for the neutrino vector in StdHep that matches NeuKin
      for (int ihep=0; ihep<= stdheparray->GetLast(); ++ihep) {
        TParticle* part = 
          dynamic_cast<TParticle*>(stdheparray->At(ihep));
        if ( !part || !th.IsNeutrino(part) || !th.IsDocStatus(part) ) continue;
        if ( !compareP4(ntpmctruth->p4neu,part) ) continue;
        // this should be the matching entry for neutrino
        
        // take the vertex from the primary lepton origin 
        // should be first daughter of neutrino
        int ihep_lepton_out = part->GetFirstDaughter();
        TParticle* outlep = 
          dynamic_cast<TParticle*>(stdheparray->At(ihep_lepton_out));
        ntpmctruth->vtxx = outlep->Vx();
        ntpmctruth->vtxy = outlep->Vy();
        ntpmctruth->vtxz = outlep->Vz();
      }

      //Fill gnumi flux info.  Added by TV, 10-31-2005
      if(finfoarray && finfoarray->GetEntriesFast() != 0){
        if(finfoarray->GetEntriesFast()<=ineu) {
          MAXMSG("NtpMC",Msg::kWarning,5)
            <<"Can not get "<<ineu<<" entry from FluxInfo Array "
            <<" which only has "<<finfoarray->GetEntriesFast()<<endl;
        }
        else {
          const REROOT_FluxInfo *rrfi = dynamic_cast<const REROOT_FluxInfo*>
            (finfoarray->At(ineu));
          FillNtpFluxInfo(ntpmctruth->flux,rrfi);
        }
      }
      if(fwgtarray && fwgtarray->GetEntriesFast() != 0){
        if(fwgtarray->GetEntriesFast()<=ineu){
          MAXMSG("NtpMC",Msg::kWarning,5)
            <<"Can not get "<<ineu<<" entry from FluxWgt Array "
            <<" which only has "<<fwgtarray->GetEntriesFast()<<endl;
        }
        else{
          const REROOT_FluxWgt *rrfw = dynamic_cast<const REROOT_FluxWgt*>
            (fwgtarray->At(ineu));
          FillNtpFluxWgt(ntpmctruth->fluxwgt,rrfw);
        }
      }
      
      MSG("NtpMC",Msg::kDebug) << *ntpmctruth << endl;
    }
    

  }
  else {
    MSG("NtpMC",Msg::kWarning) 
     << "Either NeuKin list " << neukinarray 
     << " or StdHep list " << stdheparray
     << " is missing in SimSnarlRecord.\nNtpMCTruth array will not be filled."
     << endl;
    return;
  }

  return;

}

void NtpMCModule::FillNtpMCStdHep(TClonesArray& ntpstdheparray,
                                  NtpMCSummary& ntpmcsummary,
                                  const SimSnarlRecord* simrec) {
  //
  //  Purpose:  Private method used to fill shield strip portion of ntuple 
  //            record.
  //
  //  Arguments: TClonesArray of NtpMCStdHep and SimSnarlRecord
  //  
  //  Return: status.
  // 

  MSG("NtpMC",Msg::kDebug) << "NtpMCModule::FillNtpMCStdHep" << endl;

  if ( !simrec ) {
    MSG("NtpMC",Msg::kWarning) 
    << "NtpMCModule::FillNtpMCStdHep called with null argument pointer." 
    << endl;
    return;
  }

  const TClonesArray* simstdheparray = dynamic_cast<const TClonesArray*>
                           (simrec -> FindComponent("TClonesArray","StdHep"));
  if ( !simstdheparray ) {
    MSG("NtpMC",Msg::kWarning) << "No StdHep array in SimSnarlRecord." << endl;
    return;
  }
  const TClonesArray* simdigihitarray = dynamic_cast<const TClonesArray*>
                    (simrec -> FindComponent("TClonesArray","DigiScintHits"));
  if ( !simdigihitarray ) {
    MSG("NtpMC",Msg::kWarning) << "No DigiScintHits array in SimSnarlRecord." 
                               << endl;
    return;
  }

  // Extract vldcontext and create an UgliGeomHandle
  const VldContext& vldc = *(simrec -> GetVldContext()); 
  UgliGeomHandle geom(vldc);

  // Fill fStdHepToDigiHitMap with particle associated digihits
  BuildStdHepToDigiHit(simstdheparray,simdigihitarray);  

  Int_t nstdhep = simstdheparray->GetLast()+1;
  ntpmcsummary.nstdhep = nstdhep;
  Int_t mcindex = 0;
  for ( Int_t istd = 0; istd < nstdhep; istd++ ) {
    NtpMCStdHep* ntpstdhep = new(ntpstdheparray[istd])NtpMCStdHep();
    ntpstdhep->index = istd;
    TParticle* simstdhep = dynamic_cast<TParticle*>(simstdheparray->At(istd));
    if ( simstdhep ) {
      ntpstdhep->parent[0] = simstdhep->GetMother(0);
      ntpstdhep->parent[1] = simstdhep->GetMother(1);
      ntpstdhep->child[0] = simstdhep->GetDaughter(0);
      ntpstdhep->child[1] = simstdhep->GetDaughter(1);
      ntpstdhep->IdHEP  = simstdhep->GetPdgCode(); 
      ntpstdhep->IstHEP = simstdhep->GetStatusCode();
      ntpstdhep->mass   = simstdhep->GetCalcMass();
      ntpstdhep->p4[0]  = simstdhep->Px(); // GeV
      ntpstdhep->p4[1]  = simstdhep->Py(); // GeV
      ntpstdhep->p4[2]  = simstdhep->Pz(); // GeV
      ntpstdhep->p4[3]  = simstdhep->Energy();  // GeV
      ntpstdhep->vtx[0] = simstdhep->Vx(); // Meter 
      ntpstdhep->vtx[1] = simstdhep->Vy(); // Meter 
      ntpstdhep->vtx[2] = simstdhep->Vz(); // Meter 
      ntpstdhep->vtx[3] = simstdhep->T(); // Sec (previously mm/c)
      if ( fStdHepRange ) {
        if ( istd > fStdHepRange[mcindex*2+1] ) mcindex++;
        ntpstdhep->mc = mcindex;
      }

      // Add digiscinthit information for first/last hit in the detector
      std::map<Int_t,std::vector<DigiScintHit*> >::iterator mapItr 
                                         = fStdHepToDigiHitMap.find(istd);
      if ( mapItr != fStdHepToDigiHitMap.end() ) {
        std::vector<DigiScintHit*>& hitList = mapItr->second;
        ntpstdhep->ndethit = hitList.size();
        for ( int iend=0; iend < 2; iend++ ) {
          const DigiScintHit* simdigihit = 0;
          if ( iend == 0 ) simdigihit = hitList[iend];
          else simdigihit = hitList[hitList.size()-1];
          
          NtpMCStdHepHit& ntpstdhephit = ntpstdhep->dethit[iend];

          ntpstdhephit.planeview = simdigihit->StripEndId().GetPlaneView();
          ntpstdhephit.plane = simdigihit->Plane();
          ntpstdhephit.strip = simdigihit->Strip();

          const UgliStripHandle& uglistp
                     = geom.GetStripHandle(simdigihit->StripEndId());
          // DigiScintHit stores X1/X2 instead of X0/X1
          TVector3 vlocal0(simdigihit->X1(),simdigihit->Y1(),simdigihit->Z1());

          const TVector3& vglobal0 = uglistp.LocalToGlobal(vlocal0);
          ntpstdhephit.vtx[0] = vglobal0.X();
          ntpstdhephit.vtx[1] = vglobal0.Y();
          ntpstdhephit.vtx[2] = vglobal0.Z();
          ntpstdhephit.vtx[3] = simdigihit->T1(); // DigiScintHit stores T1/T2
          ntpstdhephit.mom[0] = simdigihit->ParticlePX1();
          ntpstdhephit.mom[1] = simdigihit->ParticlePY1();
          ntpstdhephit.mom[2] = simdigihit->ParticlePZ1();
          ntpstdhephit.mom[3] = simdigihit->ParticleEnergy();
        }
      }
    }
  }

  Int_t logLevel = MsgService::Instance()->GetStream("NtpMC")->GetLogLevel();
  if ( logLevel <= Msg::kDebug ) {
    std::string entrystr = " entries:";
    if ( nstdhep == 1 ) entrystr = " entry:";
    MSG("NtpMC",Msg::kDebug) << "\nPrint stdhep array of NtpMCStdHep w/" 
                             << nstdhep << entrystr.c_str() << endl;
    for ( Int_t istd = 0; istd < nstdhep; istd++ ) {
      const NtpMCStdHep* ntpstdhep 
      = dynamic_cast<const NtpMCStdHep*>(ntpstdheparray.At(istd));
      if ( ntpstdhep->parent[0] == -1 ) 
         ntpstdhep->Print(std::cout,"","",&ntpstdheparray);
      
    }
  }

  fStdHepToDigiHitMap.clear();
  
  return;

}

   
bool NtpMCModule::BuildMCToStdHep(const TClonesArray* neukinarray,
                                  const TClonesArray* stdheparray) {
//
//  Purpose:  Private method used to build mc to stdhep map
//
//  Return: false if build failed (an error message is also printed), else true
//

   MSG("NtpMC",Msg::kDebug) << "NtpMCModule::BuildMCToStdHep()"
                            << endl;
 
   if ( fStdHepRange ) delete [] fStdHepRange; fStdHepRange = 0;
   
   Int_t nneukin = neukinarray->GetLast()+1;
   if ( nneukin <= 0 ) return true; // nothing to build
   
   Int_t nstdhep = stdheparray->GetLast()+1;
   fStdHepRange = new Int_t[nneukin*2];
   for ( int ineu = 0; ineu < nneukin; ineu++ ) {
     fStdHepRange[ineu*2+0] = -1;
     fStdHepRange[ineu*2+1] = -1;
   }
   
   Int_t mcindex = -1;
   Int_t nprim = 0;
   bool  prevchild = true;
   for ( Int_t istd = 0; istd < nstdhep; istd++ ) {
     TParticle* simstdhep = dynamic_cast<TParticle*>(stdheparray->At(istd));
     if ( simstdhep->GetMother(0) == -1 && simstdhep->GetMother(1) == -1 ) {
       // primary, check to make sure its new
       nprim++;
       if ( nprim > 2 || prevchild ) { 
         mcindex++; 
         nprim = 1; 
         if ( mcindex < nneukin ) {
           fStdHepRange[mcindex*2+0] = istd; // range min
         }
         else {
           MSG("NtpMC",Msg::kError) 
           << "FillNtpMCTruth:\nBreakdown in procedure to match stdhep to "
           << "associated mc entry.\nmc to stdhep indexing will not be done!"
           << endl;
           if ( fStdHepRange ) delete [] fStdHepRange; fStdHepRange = 0;
           return false;
         }
       }
       prevchild = false; 
     }
     else {
       prevchild = true;
       nprim = 0;
     }
      
     fStdHepRange[mcindex*2+1] = istd; // range max
     
   }   

   // At end, check grand total
   if ( mcindex != nneukin - 1  ) {
     MSG("NtpMC",Msg::kError) 
      << "FillNtpMCTruth:\nBreakdown in procedure to match stdhep to " 
      << "associated mc entry.\nFound " << mcindex + 1 << " primaries "
      << "in stdhep array, but nneukin = " << nneukin << ".\n"
      << "mc to stdhep indexing will not be done!" << endl;
     if ( fStdHepRange ) delete [] fStdHepRange; fStdHepRange = 0;
     return false;
   }

   return true;
   
}

bool NtpMCModule::BuildStdHepToDigiHit(const TClonesArray* simstdheparray,
                                       const TClonesArray* simdigihitarray) {
//
//  Purpose:  Private method to fill fStdHepToDigiHitMap with
//            <stdhep particle index, vector of digiscinthits> pairs.
//            Only those digihits directly generated by the stdhep particle
//            are stored in the list of digihits associated with the
//            particle.  Veto shield hits are excluded.
//
//  Return: false if build failed (an error message is also printed), else true
//

   MSG("NtpMC",Msg::kDebug) << "NtpMCModule::BuildStdHepToDigiHit()"
                            << endl;
 
   fStdHepToDigiHitMap.clear(); 
   
   if ( !simstdheparray ) {
     MSG("NtpMC",Msg::kWarning)
     << "NtpMCModule::BuildStdHepToDigiHit received null stdheparray." 
     << endl;
     return false;
   }
   if ( !simdigihitarray ) {
     MSG("NtpMC",Msg::kWarning)
     << "NtpMCModule::BuildStdHepToDigiHit received null digihitarray." 
     << endl;
     return false;
   }

   Int_t ndigihit = simdigihitarray->GetLast()+1;
   for ( Int_t ihit = 0; ihit < ndigihit; ihit++ ) {
     DigiScintHit* simdigihit
     = dynamic_cast<DigiScintHit*>(simdigihitarray->At(ihit));
     Int_t trkId = simdigihit->TrackId(); //id of trk that produced hit
     if ( trkId < 0 ) continue; // track not present in stdhep array
     if ( (simdigihit->StripEndId()).IsVetoShield() ) continue; // skip VS hits

     // Because cosmic reroot files place the primary muon in the 0th slot 
     // in stdhep array, such that hits off the secondaries produced by
     // the muon will have trkId's 0 instead of a negative trkId, check
     // to make sure particleId in this slot is that of primary.
     const TParticle* stdhep = dynamic_cast<const TParticle*>
                               (simstdheparray->At(trkId));
     if ( stdhep->GetPdgCode() != simdigihit->ParticleId() ) continue;

     // The following will create an entry in map if not present already
     std::vector<DigiScintHit*>& hitlist = fStdHepToDigiHitMap[trkId];
     hitlist.push_back(simdigihit);
   }

   return true;
   
}

void NtpMCModule::FillNtpMCDigiScintHit(TClonesArray& ntpdigihitarray,
                                        NtpMCSummary& ntpmcsummary,
                                        const SimSnarlRecord* simrec) {
//
//  Purpose:  Private method used to fill DigiScintHit portion of ntuple
//            record.
//
//  Arguments: TClonesArray of digihits and SimSnarlRecord
//
//  Return: status.
//

   MSG("NtpMC",Msg::kDebug) << "NtpMCModule::FillNtpMCDigiScintHit"
                                                             << endl;

   if ( !simrec ) {
     MSG("NtpMC",Msg::kWarning)
     << "NtpMCModule::FillNtpMCDigiScintHit called with null argument pointer."
     << endl;
     return;
   }

   // Extract vldcontext and create an UgliGeomHandle
   const VldContext& vldc = *(simrec -> GetVldContext());
   
   UgliGeomHandle geom(vldc);

   const TClonesArray* simdigihitarray = dynamic_cast<const TClonesArray*>
                    (simrec -> FindComponent("TClonesArray","DigiScintHits"));
   if ( !simdigihitarray ) {
     MSG("NtpMC",Msg::kWarning)
     << "No DigiScintHit array in SimSnarlRecord." << endl;
     return;
   }

   Int_t ndigihit = simdigihitarray->GetLast()+1;
   ntpmcsummary.ndigihit = ndigihit;
   for ( Int_t ihit = 0; ihit < ndigihit; ihit++ ) {
     NtpMCDigiScintHit* ntpdigihit
     = new(ntpdigihitarray[ihit])NtpMCDigiScintHit();
     ntpdigihit->index = ihit;
     DigiScintHit* simdigihit
     = dynamic_cast<DigiScintHit*>(simdigihitarray->At(ihit));
     if ( simdigihit ) {
       ntpdigihit->planeview = simdigihit->StripEndId().GetPlaneView();
       ntpdigihit->strip = simdigihit->Strip();
       ntpdigihit->plane = simdigihit->Plane();
       ntpdigihit->t0 = simdigihit->T1(); // time at start of segment
       ntpdigihit->t1 = simdigihit->T2(); // time at end of segment
       ntpdigihit->dS = simdigihit->DS(); // total path length
       ntpdigihit->dE = simdigihit->DE(); // total energy deposited
       ntpdigihit->pE = simdigihit->ParticleEnergy(); // parent particle energy
       ntpdigihit->pId = simdigihit->ParticleId(); // parent particle Id
       ntpdigihit->trkId = simdigihit->TrackId(); //id of trk that produced hit
       ntpdigihit->failbits = simdigihit->FailBits();

       const UgliStripHandle& uglistp
        = geom.GetStripHandle(simdigihit->StripEndId());
       TVector3 vlocal0(simdigihit->X1(),simdigihit->Y1(),simdigihit->Z1());
       const TVector3& vglobal0 = uglistp.LocalToGlobal(vlocal0);
       ntpdigihit->x0 = vglobal0.X();
       ntpdigihit->y0 = vglobal0.Y();
       ntpdigihit->z0 = vglobal0.Z();
       TVector3 vlocal1(simdigihit->X2(),simdigihit->Y2(),simdigihit->Z2());
       const TVector3& vglobal1 = uglistp.LocalToGlobal(vlocal1);
       ntpdigihit->x1 = vglobal1.X();
       ntpdigihit->y1 = vglobal1.Y();
       ntpdigihit->z1 = vglobal1.Z();
     }
     MSG("NtpMC",Msg::kVerbose) << *ntpdigihit << endl;
   }

   return;

}

void NtpMCModule::FillNtpMCPhotonResult(NtpMCPhotonResult& ntpmcphoton,
                                        const SimSnarlRecord* simrec) {
//
//  Purpose:  Private method used to fill photon result portion of ntuple
//            record.
//
//  Arguments: reference to photon result and SimSnarlRecord
//
//  Return: none.
//

   MSG("NtpMC",Msg::kDebug) << "NtpMCModule::FillNtpMCPhotonResult"
                            << endl;

   if ( !simrec ) {
     MSG("NtpMC",Msg::kWarning)
     << "NtpMCModule::FillNtpMCPhotonResult called with null argument pointer."
     << endl;
     return;
   }

   const PhotonEventResult* per = dynamic_cast<const PhotonEventResult*>
                              (simrec->FindComponent("PhotonEventResult"));
   if ( !per ) {
     MAXMSG("NtpMC",Msg::kWarning,1) 
      << "No PhotonEventResult component in SimSnarlRecord.\n"
      << "NtpMC photon data member will not be filled." << endl;
      return;
   }
 
   ntpmcphoton.hitsDiscardedGeom   = per->hitsDiscardedGeom;
   ntpmcphoton.hitsDiscardedBad    = per->hitsDiscardedBad;
   ntpmcphoton.totalHits           = per->totalHits;
   ntpmcphoton.totalStripsHit      = per->totalStripsHit;
   ntpmcphoton.bluePhotons         = per->bluePhotons;
   ntpmcphoton.greenPhotons        = per->greenPhotons;
   ntpmcphoton.bluePhotons_nonprescaled = per->bluePhotons_nonprescaled;
   ntpmcphoton.greenPhotons_nonprescaled = per->greenPhotons_nonprescaled;
   ntpmcphoton.totalPE             = per->totalPE;
   ntpmcphoton.totalPixels         = per->totalPixels;
   ntpmcphoton.totalHitEnergy      = per->totalHitEnergy;
   ntpmcphoton.energyDiscardedGeom = per->energyDiscardedGeom;
   ntpmcphoton.energyDiscardedBad  = per->energyDiscardedBad;

   MSG("NtpMC",Msg::kVerbose) << ntpmcphoton << endl;

   return;

}

void NtpMCModule::FillNtpMCDetSimResult(NtpMCDetSimResult& ntpmcdetsim,
                                        const SimSnarlRecord* simrec) {
//
//  Purpose:  Private method used to fill detsim result portion of ntuple
//            record.
//
//  Arguments: reference to detsim result and SimSnarlRecord
//
//  Return: none.
//

   MSG("NtpMC",Msg::kDebug) << "NtpMCModule::FillNtpMCDetSimResult"
                            << endl;

   if ( !simrec ) {
     MSG("NtpMC",Msg::kWarning)
     << "NtpMCModule::FillNtpMCDetSimResult called with null argument pointer."
     << endl;
     return;
   }

   const SimEventResult* det = dynamic_cast<const SimEventResult*>
                            (simrec->FindComponent("SimEventResult"));
   if ( !det ) {
     MAXMSG("NtpMC",Msg::kWarning,1) 
      << "No SimEventResult component in SimSnarlRecord.\n"
      << "NtpMC detsim data member will not be filled." << endl;
      return;
   }
 
   ntpmcdetsim.timeShift             = det->timeShift;
   ntpmcdetsim.nPE                   = det->npe;
   ntpmcdetsim.hitPixels             = det->hitPixels;
   ntpmcdetsim.hitPixelsWithXtalk    = det->hitPixelsWithXtalk;
   ntpmcdetsim.digitsAfterFETrigger  = det->digitsAfterFETrigger;
   ntpmcdetsim.digitsAfterSpars      = det->digitsAfterSpars;
   ntpmcdetsim.digitsAfterDaqTrigger = det->digitsAfterDaqTrigger;
   ntpmcdetsim.totalPE               = det->totalPe;
   ntpmcdetsim.totalCharge           = det->totalCharge;
   ntpmcdetsim.adcsAfterFETrigger    = det->adcsAfterFETrigger;
   ntpmcdetsim.adcsAfterSpars        = det->adcsAfterSpars;
   ntpmcdetsim.adcsAfterDaqTrigger   = det->adcsAfterDaqTrigger;
   ntpmcdetsim.bigSnarl              = det->bigSnarl;
   ntpmcdetsim.snarls                = det->snarls;
   ntpmcdetsim.snarlDigits           = new UInt_t[ntpmcdetsim.snarls];
   ntpmcdetsim.snarlTrigger          = new UInt_t[ntpmcdetsim.snarls];
   ntpmcdetsim.snarlAdcs             = new Float_t[ntpmcdetsim.snarls];
   
   for (int is = 0; is < ntpmcdetsim.snarls; is++ ) {
     ntpmcdetsim.snarlDigits[is]  = det->snarl_digits[is];
     ntpmcdetsim.snarlTrigger[is] = det->snarl_trigger[is];
     ntpmcdetsim.snarlAdcs[is]    = det->snarl_adcs[is];
   }

   MSG("NtpMC",Msg::kVerbose) << ntpmcdetsim << endl;

   return;

}

void NtpMCModule::FillNtpFluxInfo(NtpMCFluxInfo &ntpfinfo,
                                  const REROOT_FluxInfo *rrfinfo)
{
  MSG("NtpMC",Msg::kDebug) << "NtpMCModule::FillNtpFlux" << endl;

  if ( !rrfinfo ) {
    MAXMSG("NtpMC",Msg::kWarning,5) 
    << "NtpMCModule::FillNtpFlux called with null argument pointer." 
    << endl;
    return;
  }

  ntpfinfo.index = rrfinfo->ID();
  ntpfinfo.fluxrun = rrfinfo->FluxRun();
  ntpfinfo.fluxevtno = rrfinfo->FluxEvtNo();
  ntpfinfo.ndxdz = rrfinfo->Ndxdz();
  ntpfinfo.ndydz = rrfinfo->Ndydz();
  ntpfinfo.npz = rrfinfo->Npz();
  ntpfinfo.nenergy = rrfinfo->Nenergy();
  ntpfinfo.ndxdznear = rrfinfo->NdxdzNear();
  ntpfinfo.ndydznear = rrfinfo->NdydzNear();
  ntpfinfo.nenergynear = rrfinfo->NenergyN();
  ntpfinfo.nwtnear = rrfinfo->NWtNear();
  ntpfinfo.ndxdzfar = rrfinfo->NdxdzFar();
  ntpfinfo.ndydzfar = rrfinfo->NdydzFar();
  ntpfinfo.nenergyfar = rrfinfo->NenergyF();
  ntpfinfo.nwtfar = rrfinfo->NWtFar();
  ntpfinfo.norig = rrfinfo->Norig();
  ntpfinfo.ndecay = rrfinfo->Ndecay();
  ntpfinfo.ntype = rrfinfo->Ntype();
  ntpfinfo.vx = rrfinfo->Vx();
  ntpfinfo.vy = rrfinfo->Vy();
  ntpfinfo.vz = rrfinfo->Vz();
  ntpfinfo.pdpx = rrfinfo->pdPx();
  ntpfinfo.pdpy = rrfinfo->pdPy();
  ntpfinfo.pdpz = rrfinfo->pdPz();
  ntpfinfo.ppdxdz = rrfinfo->ppdxdz();
  ntpfinfo.ppdydz = rrfinfo->ppdydz();
  ntpfinfo.pppz = rrfinfo->pppz();
  ntpfinfo.ppenergy = rrfinfo->ppenergy();
  ntpfinfo.ppmedium = rrfinfo->ppmedium();
  ntpfinfo.ptype = rrfinfo->ptype();
  ntpfinfo.ppvx = rrfinfo->ppvx();
  ntpfinfo.ppvy = rrfinfo->ppvy();
  ntpfinfo.ppvz = rrfinfo->ppvz();
  ntpfinfo.muparpx = rrfinfo->muparpx();
  ntpfinfo.muparpy = rrfinfo->muparpy();
  ntpfinfo.muparpz = rrfinfo->muparpz();
  ntpfinfo.mupare = rrfinfo->mupare();
  ntpfinfo.necm = rrfinfo->Necm();
  ntpfinfo.nimpwt = rrfinfo->Nimpwt();
  ntpfinfo.xpoint = rrfinfo->xpoint();
  ntpfinfo.ypoint = rrfinfo->ypoint();
  ntpfinfo.zpoint = rrfinfo->zpoint();
  ntpfinfo.tvx = rrfinfo->tvx();
  ntpfinfo.tvy = rrfinfo->tvy();
  ntpfinfo.tvz = rrfinfo->tvz();
  ntpfinfo.tpx = rrfinfo->tpx();
  ntpfinfo.tpy = rrfinfo->tpy();
  ntpfinfo.tpz = rrfinfo->tpz();
  ntpfinfo.tptype = rrfinfo->tptype();
  ntpfinfo.tgen = rrfinfo->tgen();
  // NUGPINFO data
  ntpfinfo.tgptype = rrfinfo->tgptype();
  ntpfinfo.tgppx = rrfinfo->tgppx();
  ntpfinfo.tgppy = rrfinfo->tgppy();
  ntpfinfo.tgppz = rrfinfo->tgppz();
  ntpfinfo.tprivx = rrfinfo->tprivx();
  ntpfinfo.tprivy = rrfinfo->tprivy();
  ntpfinfo.tprivz = rrfinfo->tprivz();
  ntpfinfo.beamx = rrfinfo->beamx();
  ntpfinfo.beamy = rrfinfo->beamy();
  ntpfinfo.beamz = rrfinfo->beamz();
  ntpfinfo.beampx = rrfinfo->beampx();
  ntpfinfo.beampy = rrfinfo->beampy();
  ntpfinfo.beampz = rrfinfo->beampz();

  return;
}

void NtpMCModule::FillNtpFluxWgt(NtpMCFluxWgt &ntpfwgt,
                                 const REROOT_FluxWgt *rrfwgt)
{
  MSG("NtpMC",Msg::kDebug) << "NtpMCModule::FillNtpFluxWgt" << endl;

  if ( !rrfwgt ) {
    MAXMSG("NtpMC",Msg::kWarning,5) 
    << "NtpMCModule::FillNtpFluxWgt called with null argument pointer." 
    << endl;
    return;
  }

  ntpfwgt.index = rrfwgt->ID();
  ntpfwgt.version = rrfwgt->Version();
  ntpfwgt.weight = rrfwgt->Weight();
  ntpfwgt.weighterr = rrfwgt->WeightErr();
  strcpy(ntpfwgt.beam,rrfwgt->Beam());  

  return;
}

---