PTSimStack.cxx

Go to the documentation of this file.

//
// PTSimStack
//
// PTSimStack is a concrete implementation of MCAppStack to act as the
// particle stack during particle transport.
//
// S. Kasahara  05/04
//
// A PTSimStack extends a MCAppStack with methods related to the storage
// of primaries and secondaries in the output "stdhep" TClonesArray.

#include <TClonesArray.h>
#include <TVirtualMC.h>
#include <TPDGCode.h>
#include <TDatabasePDG.h>

#include "Util/UtilIstHEP.h"
#include "Conventions/Munits.h"
#include "Digitization/DigiScintHit.h"
#include "MessageService/MsgService.h"
#include "ParticleTransportSim/PTSim.h"
#include "ParticleTransportSim/PTSimStack.h"

ClassImp(PTSimStack)

CVSID("$Id: PTSimStack.cxx,v 1.33 2009/06/22 03:51:01 kasahara Exp $");

std::ostream& operator << (std::ostream& os, const PTSimStack& ps)
                          { return ps.Print(os); }
  
//_____________________________________________________________________________
PTSimStack::PTSimStack() : fStdHepSaveByRange(false),
                           fStdHepSaveSibling(1),
                           fStdHepSaveAncestor(1),
                           fStdHepSelectMask(PTSim::kMomentum),
                           fStdHepHitThr(0.001),
                           fStdHepThr(0.15) {
  // Default constructor

  MSG("PTSim",Msg::kVerbose) << "PTSimStack def ctor @ " << this << endl;
  
}

//_____________________________________________________________________________
PTSimStack::~PTSimStack() {
  // Destructor

  MSG("PTSim",Msg::kVerbose) << "PTSimStack dtor @ " << this << endl;
  Reset();
  
}

//_____________________________________________________________________________
void  PTSimStack::SetStdHepThrByType(double pthresh, TMCProcess process,
                                     Int_t pdgId, Int_t parentId) {
  // Sets the pthresh(GeV/c) threshold over which secondaries produced by the
  // particle type will be stored to the output stdhep array.
  // Particle type is defined as: pdgId*1000 + process.
  // Special cases:
  // 1)If pdgId=kRootino, the threshold is applied to all particles generated
  //   by the specified production mechanism process.  
  // 2)If the production mechanism is process=kPNoProcess, the threshold 
  //   is applied to all particles of the specified pdgId, regardless of 
  //   the production mechanism.
  // For the specified particle type, the specified pthresh overrides the
  // generic momentum threshold defined by fStdHepThr.

  // Particle Name
  std::string particlename = "all particle types";
  if ( pdgId != kRootino ) {
    const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
    particlename = "???";
    if ( dbpdg.GetParticle(pdgId) ) particlename 
                   = dbpdg.GetParticle(pdgId)->GetName();
  }
  
  // Production Process Name
  std::string processname = "all production processes";
  if ( process != kPNoProcess ) {
    processname = TMCProcessName[process];
  }
  
  if ( process == kPNoProcess ) {
    MSG("PTSim",Msg::kDebug) << "StdHepThrByType "
                             << pthresh << "(GeV/c)" 
                             << " for particle type " << pdgId << "("
                             << particlename.c_str() << ")" << endl;
  }
  else if ( pdgId == kRootino ) {
    MSG("PTSim",Msg::kDebug) << "StdHepThrByType "
                             << pthresh << "(GeV/c)" 
                             << " for production process "
                             << process << "(" << processname.c_str() << ")"
                             << endl;
  }
  else {                   
    MSG("PTSim",Msg::kDebug) << "StdHepThrByType "
                             << pthresh << "(GeV/c)" 
                             << " for particle type " << pdgId << "("
                             << particlename.c_str() << ")" 
                             << " and production process "
                             << process << "(" << processname.c_str() << ")"
                             << endl;
  }

  fStdHepThrByTypeMap[PTSimStdHepType(process,pdgId,parentId)] = pthresh;

}     

//_____________________________________________________________________________
void PTSimStack::PushTrack(Int_t toBeDone, Int_t parentId, Int_t pdg,
                           Double_t px, Double_t py, Double_t pz, Double_t e,
                           Double_t vx, Double_t vy, Double_t vz, Double_t tof,
                           Double_t polx, Double_t poly, Double_t polz,
                           TMCProcess mech, Int_t& ntr, Double_t weight,
                           Int_t status) {
  // Specialized MCAppStack<T>::PushTrack(...) to handle case of storing
  // decay parent at point of decay.  Needs to be pushed before decay
  // children to keep things in the right order.
  // See MCAppStack<T>::PushTrack(...) for full list of arguments.
  // 
  // ntr - track number of this particle in fParticles, output to caller
  // 

  Int_t mechstatus = status;
  if ( mech != kPPrimary) mechstatus = (Int_t)UtilIstHEP::GetProdMethod(mech)
                                     + UtilIstHEP::kProdMethodOffset;
  
  if ( mech == kPDecay ) {
    // Check last daughter position of parent.  
    const PTSimParticle* parent = GetParticle(parentId);
    const PTSimParticle* lastchild = dynamic_cast<const PTSimParticle*>
                            (parent->GetChild(parent->GetNChildren()-1));
    Int_t decayparentId = -1;
    Int_t decayparentcode = UtilIstHEP::kMDecay 
                          + UtilIstHEP::kProdMethodOffset + 1000;
    if ( lastchild != 0 ) {
      // Check to see if parent at decay point has already been pushed
      if ( lastchild -> GetStatusCode() == decayparentcode )
                                      decayparentId = lastchild->GetID();
    }
    if ( decayparentId < 0 ) {
      // Push parent at point of decay now
      TLorentzVector pos;
      gMC->TrackPosition(pos);
      TLorentzVector mom;
      gMC->TrackMomentum(mom);
      Int_t parentpdg = gMC->TrackPid();
      // First argument 0 means already tracked (don't push to stack)
      MCAppStack<PTSimParticle>::PushTrack(0,parentId,parentpdg,
          mom.X(),mom.Y(),mom.Z(),mom.E(),pos.X(),pos.Y(),pos.Z(),pos.T(),
          0.,0.,0.,kPTransportation,decayparentId,1.,decayparentcode);
    }
    // Now push the decay product with the decayparentId as parent
    MCAppStack<PTSimParticle>::PushTrack(toBeDone,decayparentId,pdg,
                                        px,py,pz,e,vx,vy,vz,tof,
                                        polx,poly,polz,mech,ntr,weight,
                                        mechstatus);
  }
  else {
    // All other cases (not decay)
    MCAppStack<PTSimParticle>::PushTrack(toBeDone,parentId,pdg,
                                        px,py,pz,e,vx,vy,vz,tof,
                                        polx,poly,polz,mech,ntr,weight,
                                        mechstatus);
  }
  
}

//_____________________________________________________________________________
Double_t PTSimStack::GetStdHepThrByType(TMCProcess process, Int_t pdgId,
                                        Int_t parentId) const {
  // Return momentum threshold for particle produced by production
  // mechanism process, pdg code pdgId, and parentId.   
  // The momentum threshold is determined by looking at the 
  // user configured momentum thresholds and taking the minimum
  // of the thresholds of all those that match the user's type.
  // A "match" to a user-configured threshold by type is made if:
  // 1)The process is an exact match or kPNoProcess has been specified in
  //   user configured momentum threshold by type, and:
  // 2)The pdgId is an exact match or kRootino has been specified in
  //   user configured momentum threshold by type, and:
  // 3)The parentId is an exact match (can be kRootino)
  // If no match, -1 is returned.
  //

  Double_t pthreshbytype = -1.;
  
  PTSimStdHepType thrtype(process,pdgId,parentId);
  std::map<PTSimStdHepType,double>::const_iterator threshMapItr;
  threshMapItr = fStdHepThrByTypeMap.find(thrtype);
  if ( threshMapItr != fStdHepThrByTypeMap.end() ) {
    // exact match
    pthreshbytype = threshMapItr -> second;
  }
  
  // Try again with particle undefined (kRootino) & process
  thrtype = PTSimStdHepType(process,kRootino,parentId);
  threshMapItr = fStdHepThrByTypeMap.find(thrtype);
  if ( threshMapItr != fStdHepThrByTypeMap.end() ) {
    if ( pthreshbytype < 0. ) pthreshbytype = threshMapItr -> second;
    else pthreshbytype = TMath::Min(pthreshbytype,threshMapItr->second);
  }

  // And again with particle & process undefined (kPNoProcess)
  thrtype = PTSimStdHepType(kPNoProcess,pdgId,parentId);
  threshMapItr = fStdHepThrByTypeMap.find(thrtype);
  if ( threshMapItr != fStdHepThrByTypeMap.end() ) {
    if ( pthreshbytype < 0. ) pthreshbytype = threshMapItr -> second;
    // If particle passed two different particle type definitions
    // using generic particle type for first and a generic production 
    // mechanism for second, take minimum of two thresholds.
    else pthreshbytype = TMath::Min(pthreshbytype,threshMapItr->second);
  }

  return pthreshbytype;
  
}

//_____________________________________________________________________________
void PTSimStack::FillStdHepArray(TClonesArray* stdheparray, Int_t snarl,
                                 TClonesArray* hitarray, Int_t firsthit) {
  // Fill stdheparray with all primaries and selected secondaries.
  // The snarl number is used to determine if secondaries were requested
  // for storage for this snarl.
  // Adjust hitarray hit track indices starting with hit "firsthit" to 
  // match those in newly filled section of stdheparray.
  //
  // Note: Neither hitarray or stdheparray are adopted by this method.
  //       Ownership belongs to caller.

  MSG("PTSim",Msg::kDebug) << "PTSimStack::FillStdHepArray for snarl "
                             << snarl << endl;

  if ( !stdheparray ) {
    MSG("PTSim",Msg::kFatal) << "PTSimStack::FillStdHepArray: Null "
                             << "stdheparray TClonesArray ptr. Abort." << endl;
    abort();
  }
  if ( !hitarray ) {
    MSG("PTSim",Msg::kFatal) << "PTSimStack::FillStdHepArray: Null "
                             << "hitarray TClonesArray ptr. Abort." << endl;
    abort();
  }
  
  Int_t firsttrk = stdheparray -> GetEntriesFast();
  Int_t nstdhep = firsttrk;
  Int_t npart = GetNtrack();

  MSG("PTSim",Msg::kDebug) << "PTSimStack::FillStdHepArray.  Start with "
                      << nstdhep << " particles in input stdhep array and " 
                      << npart << " particles in stack." << endl;
  
  // Determine if secondaries are to be stored for this snarl
  bool savesecond = SaveSecondariesForSnarl(snarl);

  Int_t decayparentcode = UtilIstHEP::kMDecay 
                        + UtilIstHEP::kProdMethodOffset + 1000;

  std::stack<const PTSimParticle*> ministack;
  // First go through and tag all particles that are to be saved to
  // secondary array
  for ( Int_t ip = 0; ip < npart; ip++ ) {
    // Determine if particle is primary or passes user-defined selection cuts
    PTSimParticle* particle = const_cast<PTSimParticle*>(GetParticle(ip));
    
    if ( particle -> GetProcess() == kPPrimary ) {
      // primary

      fPartToStdHep.insert(make_pair(particle->GetID(),nstdhep));
      TParticle* stdhep = 
       new((*stdheparray)[nstdhep++])TParticle(*(particle->GetTParticle()));
      // Convert vertex to meters
      stdhep -> SetProductionVertex(stdhep->Vx()*Munits::cm,
                                    stdhep->Vy()*Munits::cm,
                                    stdhep->Vz()*Munits::cm,stdhep->T());

      // Adjust indices of parents/children
      if ( particle->GetParentID(0) != -1 ) 
        stdhep->SetFirstMother(particle->GetParentID(0)+firsttrk);
      if ( particle->GetParentID(1) != -1 ) 
        stdhep->SetLastMother(particle->GetParentID(1)+firsttrk);
      
      if ( particle -> GetNChildren() > 0 ) {
        if ( particle->GetStatusCode() == 1 ) {
          ministack.push(particle); // prime stack for pushing secondaries
        }
        else {
          Int_t nchild = particle->GetNChildren();
          if ( particle->GetChild(nchild-1)->GetStatusCode() == 999 ) {
            // last child is rootino
            nchild += -1;
          }
          if (  nchild > 0 ) {
            stdhep->SetFirstDaughter(particle->GetChildID(0)+firsttrk);
            stdhep->SetLastDaughter(particle->GetChildID(nchild-1)+firsttrk);
          }
        }
      }      
    }
    else if ( savesecond ) {
      // Never save decay parents directly unless selected by decay products
      if ( decayparentcode != particle->GetStatusCode() ) {
        if ( ParticleIsSelected(particle) ) {
          // Save secondary.  Mark ancestors to save, and optionally
          // siblings and siblings of direct ancestors to save as well.
          particle -> SetToSave(true,fStdHepSaveAncestor,fStdHepSaveSibling);
        }
      }
    }
  }
    
  while ( !ministack.empty() ) {
    const PTSimParticle* particle = ministack.top();
    ministack.pop();  
    
    // Loop over children, push to stdhep array so that they are in sequence
    // because VMC w/Geant4 doesn't guarantee this.
    // Also push children with children to ministack for further processing
    
    // Secondary set to save should have mother already saved
    Int_t motherId = particle->GetID();
    std::map<int,int>::iterator idItr = fPartToStdHep.find(motherId);
    if (idItr == fPartToStdHep.end() ) {
      if (idItr == fPartToStdHep.end() ) {
        MSG("PTSim",Msg::kError) << "PTSimStack::FillStdHepArray motherId " 
                 << motherId << " not in map. abort." << endl;
        abort();
      }
    }
    Int_t motherstdindex = idItr->second;

    Int_t nchildren = particle->GetNChildren();
    for ( Int_t ic = 0; ic < nchildren; ic++ ) {
      const PTSimParticle* child = dynamic_cast<const PTSimParticle*>
                                             (particle->GetChild(ic));

      if ( child -> IsSetToSave() ) {
        // Particle is to be filled in stdhep array
        TParticle* stdhep 
          = new((*stdheparray)[nstdhep])TParticle(*(child->GetTParticle()));
          
        fPartToStdHep.insert(make_pair(child->GetID(),nstdhep));
        stdhep -> SetFirstMother(TMath::Abs(motherstdindex)); 
        stdhep -> SetLastMother(TMath::Abs(motherstdindex));
        // Convert vertex to meters
        stdhep -> SetProductionVertex(stdhep->Vx()*Munits::cm,
                                      stdhep->Vy()*Munits::cm,
                                      stdhep->Vz()*Munits::cm,stdhep->T());
        TParticle* ancestor = dynamic_cast<TParticle*>
                              (stdheparray->At(TMath::Abs(motherstdindex)));
        if(ancestor->GetFirstDaughter()< 0)
                                  ancestor->SetFirstDaughter(nstdhep);
        ancestor -> SetLastDaughter(nstdhep);
        nstdhep++;
      }
      else {
        // Required because digihit may not have direct parent stored
        // in output stdhep array
        fPartToStdHep.insert(make_pair(child->GetID(),
                                      -TMath::Abs(motherstdindex)));
      }

      // push particle to stack if it has children
      if ( child->GetNChildren() > 0 ) {
        ministack.push(child);
      }

    }
  }
  
  // Now adjust hit track indices
  AdjustHitTrackId(hitarray,firsthit);
  
  Int_t nfinal = stdheparray->GetEntriesFast();
  MSG("PTSim",Msg::kDebug) << "FillStdHepArray.  Finish with "
                           << nfinal << " particles in stdheparray." << endl;

 
}

//_____________________________________________________________________________
void PTSimStack::Reset() {
  MSG("PTSim",Msg::kVerbose) << "PTSimStack::Reset " << endl;
  
  MCAppStack<PTSimParticle>::Reset();
  fPartToStdHep.clear();
  
}

  
//_____________________________________________________________________________
bool PTSimStack::SaveSecondariesForSnarl(Int_t snarl) const {
  // Protected method.  
  // Used to determine if secondaries should be saved for input argument 
  // snarl given secondary save configuration parameters.
  
  bool savesecond = false;
  if ( !fStdHepSaveByRange ) {
    for ( unsigned int is = 0; is < fStdHepSave.size(); is++ ) {
      if ( snarl == fStdHepSave[is] ) { savesecond = true; break; }
    }
  }
  else {
    if ( snarl <= fStdHepSave[1] && snarl >= fStdHepSave[0] ) {
      savesecond = true;
    }
  }

  return savesecond;
  
}

  
//_____________________________________________________________________________
bool PTSimStack::ParticleIsSelected(const PTSimParticle* particle) const {
  // Protected method.
  // According to StdHepXXX configuration parameters, pass or fail
  // input particle.
  // The criteria for passing a particle are:
  // 1) Primaries are always passed.
  // 2) Secondaries are passed according to the bits in the fStdHepSelectMask.
  //    a) If (fStdHepSelectMask & kMomentum), all secondaries with 
  //       momentum above threshold are stored. The momentum threshold is 
  //       determined using:
  //       i) A momentum threshold appropriate for the secondary particle type/
  //          production process (TMCProcess) that resulted in the generation 
  //          of the secondary, if specified by the user, OR
  //      ii) A default momentum threshold (fStdHepThr) applied to all
  //          secondaries generated by processes not covered by i).
  //    b) If (fStdHepSelectMask & kHit), all secondaries 
  //       resulting in energy deposition hit above threshold are stored.
  //       The hit threshold is determined using:
  //       i) An energy threshold (fStdHepHitThr).
  //    The user may specify any combination of the 2 options and an OR
  //    will be taken of the results of the individual options to determine
  //    the final pass/fail.

  // Primaries always pass
  if ( particle -> GetProcess() == kPPrimary ) return true;
  
  // Secondaries
  if ( fStdHepSelectMask & PTSim::kHit ) {
    if ( particle -> HasHitAboveThresh() ) return true;
  }
  
  if ( fStdHepSelectMask & PTSim::kMomentum ) {

    // Pass if above momentum threshold
    // Look for a match to StdHepThrByType specification. If multiple matches,
    // take minimum threshold.  If no match, take default momentum threshold
    // specified by StdHepThr.
    Int_t pdgId = particle -> GetPdgCode();
    TMCProcess process = particle -> GetProcess();
    // try once with no defined parent
    Double_t pthresh = GetStdHepThrByType(process,pdgId,kRootino);
    // Loop over parents 
    for ( UInt_t ip = 0; ip < particle->GetNParents(); ip++ ) {
      Int_t parentId = particle->GetParent(ip)->GetPdgCode();
      Double_t pthreshbytype = GetStdHepThrByType(process,pdgId,parentId);
      if ( pthreshbytype >= 0 ) {
        if ( pthresh >= 0 ) pthresh = TMath::Min(pthresh,pthreshbytype);
        else pthresh = pthreshbytype;
      }
    }
    if ( pthresh < 0 ) pthresh = GetStdHepThr(); // take def thresh if no match
    if ( particle -> GetTParticle() -> P() >= pthresh ) return true; 
  }
    
  // Default is no secondaries stored
  return false;

}

  
//_____________________________________________________________________________
void PTSimStack::AdjustHitTrackId(TClonesArray* hitarray, 
                                  Int_t firsthit) const {
  // Protected method.  Adjust digiscinthits' track id to that of the
  // tracks stored in the output stdhep array.

  for ( Int_t ihit = firsthit; ihit < hitarray->GetEntriesFast(); ihit++ ) {
    DigiScintHit* hit = dynamic_cast<DigiScintHit*>(hitarray->At(ihit));
    Int_t trackid = hit->TrackId();
    Int_t stdhepid = GetStdHepIndex(trackid);
    // transition from fParticles Id to stdhep index
    hit -> SetTrackId(stdhepid); 
  }

}

//_____________________________________________________________________________
Int_t PTSimStack::GetStdHepIndex(UInt_t id) const {
  // Protected method.
  // Return stdhep array index of particle corresponding to 
  // fParticles id.  The map fPartToStdHep is created during FillStdHepArray

  std::map<int,int>::const_iterator idItr = fPartToStdHep.find(id);
  if (idItr == fPartToStdHep.end() ) {
    MSG("PTSim",Msg::kError) << "PTSimStack::GetStdHepIndex for trk Id " << id
                             << " not in map. abort." << endl;
    abort();
  }
  return idItr->second;

}

---