Fatras AlgTool to create a ISFParticle at a volume entry/exit. More...

#include <TransportTool.h>

Inheritance diagram for iFatras::TransportTool:

Collaboration diagram for iFatras::TransportTool:

Public Member Functions
	TransportTool (const std::string &, const std::string &, const IInterface *)
	Constructor.
virtual	~TransportTool ()
	Destructor.
virtual StatusCode	initialize () override
	AlgTool initialize method.
virtual StatusCode	finalize () override
	AlgTool finalize method.
virtual ISF::ISFParticle *	process (const ISF::ISFParticle &isp, CLHEP::HepRandomEngine *) const override
	Creates a new ParticleState from a given ParticleState, universal transport tool.

Private Member Functions
template<class T>
StatusCode	retrieveTool (ToolHandle< T > &thandle)
	templated Tool retrieval - gives unique handling & look and feel

Private Attributes
ServiceHandle< IAtRndmGenSvc >	m_rndGenSvc
	Random Generator service.
CLHEP::HepRandomEngine *	m_randomEngine
std::string	m_randomEngineName
	Name of the random number stream.
bool	m_validationOutput
	Validation output with histogram service.
ToolHandle< IPhysicsValidationTool >	m_validationTool
	the ntuple
ToolHandle< IParticleDecayHelper >	m_particleDecayHelper
	Particle Decay.
ToolHandle< ISimHitCreator >	m_simHitCreatorID
	The SimHit Creation.
ToolHandle< ISimHitCreator >	m_simHitCreatorMS
ToolHandle< Trk::ITimedExtrapolator >	m_extrapolator
	The Extrapolator setup.
ToolHandle< ISF::IParticleFilter >	m_trackFilter
	Filtering setup & other ISF stuff.
ToolHandle< ISF::IParticleFilter >	m_neutralHadronFilter
ToolHandle< ISF::IParticleFilter >	m_photonFilter
ToolHandle< ISF::IParticleHelper >	m_iparticleHelper
ToolHandle< IProcessSamplingTool >	m_samplingTool
Trk::PdgToParticleHypothesis	m_pdgToParticleHypothesis
bool	m_errorPropagation
bool	m_hitsOff

Detailed Description

Fatras AlgTool to create a ISFParticle at a volume entry/exit.

universal transport tool

Author: Sharka Todorova Sarka.Todorova -at- cern.ch; Andreas.Salzburger -at- cern.ch

Definition at line 54 of file ISF_Fatras/ISF_FatrasTools/src/TransportTool.h.

Constructor & Destructor Documentation

◆ TransportTool()

iFatras::TransportTool::TransportTool	(	const std::string &	t,
		const std::string &	n,
		const IInterface *	p )

Constructor.

Definition at line 43 of file ISF_Fatras/ISF_FatrasTools/src/TransportTool.cxx.

 : base_class(t,n,p),
   m_rndGenSvc("AtDSFMTGenSvc", n),
   m_randomEngine(nullptr),
   m_randomEngineName("FatrasRnd"),
   m_validationOutput(false),
   m_validationTool(""),
   m_particleDecayHelper(""),
   m_simHitCreatorID(""),
   m_simHitCreatorMS(""),
   m_extrapolator(""),
   m_trackFilter(""),
   m_neutralHadronFilter(""),
   m_photonFilter(""),
   m_iparticleHelper(""),
   m_samplingTool(""),
   m_errorPropagation(false),
   m_hitsOff(false)
{
  // validation output section
  declareProperty( "ValidationOutput",  m_validationOutput );
  declareProperty( "PhysicsValidationTool",     m_validationTool  );
  // tool handle for the particle decayer
  declareProperty( "ParticleDecayHelper",       m_particleDecayHelper      );
  // tool handle for the track creator
  declareProperty( "SimHitCreatorID",           m_simHitCreatorID          );
  declareProperty( "SimHitCreatorMS",           m_simHitCreatorMS          );
  // tool handles for extrapolator(s)
  declareProperty( "Extrapolator",              m_extrapolator        );
  // tool handles for the filtering setup
  declareProperty( "TrackFilter",               m_trackFilter              );
  declareProperty( "NeutralFilter",             m_neutralHadronFilter      );
  declareProperty( "PhotonFilter",              m_photonFilter             );
  // tool handles ISF Framework
  declareProperty( "ParticleHelper",            m_iparticleHelper          );
  declareProperty( "TurnOffHitCreation",        m_hitsOff                  );
  declareProperty( "ProcessSamplingTool",       m_samplingTool             );
  declareProperty("DoErrorPropagation",         m_errorPropagation         );
  // service handles
  declareProperty( "RandomNumberService",       m_rndGenSvc                );
}

◆ ~TransportTool()

iFatras::TransportTool::~TransportTool ( )

virtual

Destructor.

Definition at line 91 of file ISF_Fatras/ISF_FatrasTools/src/TransportTool.cxx.

{
  delete m_randomEngine;
}

Member Function Documentation

◆ finalize()

StatusCode iFatras::TransportTool::finalize ( )

overridevirtual

AlgTool finalize method.

Definition at line 145 of file ISF_Fatras/ISF_FatrasTools/src/TransportTool.cxx.

{
  ATH_MSG_DEBUG( "finalize() successful" );
  return StatusCode::SUCCESS;
}

◆ initialize()

StatusCode iFatras::TransportTool::initialize ( )

overridevirtual

AlgTool initialize method.

Definition at line 101 of file ISF_Fatras/ISF_FatrasTools/src/TransportTool.cxx.

{
 
  ATH_MSG_DEBUG( "initialize()" );
 
  // Retrieve the tools one by one
  if (retrieveTool<iFatras::ISimHitCreator>(m_simHitCreatorID).isFailure())
      return StatusCode::FAILURE;
  if (retrieveTool<iFatras::ISimHitCreator>(m_simHitCreatorMS).isFailure())
      return StatusCode::FAILURE;
  if (retrieveTool<Trk::ITimedExtrapolator>(m_extrapolator).isFailure())
      return StatusCode::FAILURE;
  if (retrieveTool<iFatras::IParticleDecayHelper>(m_particleDecayHelper).isFailure())
      return StatusCode::FAILURE;
  if (retrieveTool<ISF::IParticleFilter>(m_trackFilter).isFailure())
      return StatusCode::FAILURE;
  if (retrieveTool<ISF::IParticleFilter>(m_photonFilter).isFailure())
      return StatusCode::FAILURE;
  if (retrieveTool<ISF::IParticleFilter>(m_neutralHadronFilter).isFailure())
      return StatusCode::FAILURE;
  if (retrieveTool<ISF::IParticleHelper>(m_iparticleHelper).isFailure())
      return StatusCode::FAILURE;
  if (retrieveTool<iFatras::IProcessSamplingTool>(m_samplingTool).isFailure())
      return StatusCode::FAILURE;
  ATH_CHECK( m_validationTool.retrieve( DisableTool{ m_validationTool.empty() || !m_validationOutput } ) );
 
  if ( m_rndGenSvc.retrieve().isFailure() ){
    ATH_MSG_FATAL( "Could not retrieve " << m_rndGenSvc );
    return StatusCode::FAILURE;
  }
  //Get own engine with own seeds:
  m_randomEngine = m_rndGenSvc->GetEngine(m_randomEngineName);
  if (!m_randomEngine) {
    ATH_MSG_FATAL( "Could not get random engine '" << m_randomEngineName << "'" );
    return StatusCode::FAILURE;
  }
  ATH_MSG_DEBUG( "finalize() successful" );
  return StatusCode::SUCCESS;
}

◆ process()

ISF::ISFParticle * iFatras::TransportTool::process	(	const ISF::ISFParticle &	isp,
		CLHEP::HepRandomEngine *	) const

overridevirtual

Creates a new ParticleState from a given ParticleState, universal transport tool.

Definition at line 155 of file ISF_Fatras/ISF_FatrasTools/src/TransportTool.cxx.

{
  // copy the current particle onto the particle clipboard
  ISF::ParticleClipboard::getInstance().setParticle( isp);
 
  // if process/generation not set, assume primary
  //if (m_validationOutput && !isp.getUserInformation()) {
  //  ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
  //  validInfo->setProcess(0);
  //  validInfo->setGeneration(0);
  //  isp.setUserInformation(validInfo);
  //}
 
  ATH_MSG_VERBOSE ("[ fatras transport ] processing particle "<<isp.pdgCode() );
 
  std::vector<Trk::HitInfo>*     hitVector = nullptr;    // hit or no hit collection
 
  int  absPdg = abs(isp.pdgCode());
  bool photon         = (absPdg == 22);
  bool geantino       = (absPdg == 999);
  bool neutrino       = (absPdg == 12) || (absPdg == 14) || (absPdg == 16);
  bool charged        = photon || geantino ? false : (isp.charge()*isp.charge() > 0) ;
 
  // particle Hypothesis for the extrapolation
 
  Trk::ParticleHypothesis pHypothesis = m_pdgToParticleHypothesis.convert(isp.pdgCode(),isp.charge());
 
  // geantinos not handled by PdgToParticleHypothesis - fix there
  if ( absPdg == 999 ) pHypothesis = Trk::geantino;
 
  // chose the filter & extrapolator
  const Trk::ITimedExtrapolator* processor = &(*m_extrapolator);
  const ISF::IParticleFilter* filter  = nullptr;
 
  if ( photon ) {
    ATH_MSG_VERBOSE( "[ fatras transport ] Photon - use transport tool for photons");
    filter        =  m_photonFilter.empty() ? nullptr :  &(*m_photonFilter);
  } else if  ( neutrino ) {
      ATH_MSG_VERBOSE( "[ fatras transport ] Particle is neutrino -> particle ignored.");
      return nullptr;
  } else if ( !charged ){
    ATH_MSG_VERBOSE( "[ fatras transport ] Particle not charged - use transport tool for neutrals");
    filter        =  m_neutralHadronFilter.empty() ? nullptr :  &(*m_neutralHadronFilter);
  } else {
    ATH_MSG_VERBOSE( "[ fatras transport ] Particle is a lepton or charged hadron");
    filter    =  m_trackFilter.empty() ? nullptr :  &(*m_trackFilter);
  }
 
  if (filter)
      ATH_MSG_DEBUG( "[ fatras transport ] Determined processor and filter.");
  else
      ATH_MSG_DEBUG( "[ fatras transport ] Determined processor w/o filter.");
 
  // now transport the particle - if it passes the filter
  if ( !processor || ( filter && !filter->passFilter(isp)) ) {
    ATH_MSG_VERBOSE( "[ fatras transport ] Filter not passed, ignore particle.");
    return nullptr;
  }
 
  ATH_MSG_VERBOSE( "[ fatras transport ] The StackParticle passed filter - starting transport.");
  // input parameters : curvilinear parameters
  Trk::CurvilinearParameters inputPar(isp.position(),isp.momentum(),isp.charge());
 
  // stable vs. unstable check
  double freepath = ( !m_particleDecayHelper.empty()) ? m_particleDecayHelper->freePath(isp) : - 1.;
  ATH_MSG_VERBOSE( "[ fatras transport ] Particle free path : " << freepath);
  // path limit -> time limit  ( TODO : extract life-time directly from decay helper )
  double tDec = freepath > 0. ? freepath : -1.;
  int decayProc = 0;
 
  // beta calculated here for further use in validation
  double mass = Trk::ParticleMasses::mass[pHypothesis];
  double mom = isp.momentum().mag();
  double beta = mom/sqrt(mom*mom+mass*mass);
 
  if ( tDec>0.) {
    tDec = tDec/beta/CLHEP::c_light + isp.timeStamp();
    decayProc = 201;
  }
 
  Trk::TimeLimit timeLim(tDec,isp.timeStamp(),decayProc);
 
  // prompt decay
  if ( freepath>0. && freepath<0.01 ) {
    if (!m_particleDecayHelper.empty()) {
      ATH_MSG_VERBOSE( "[ fatras transport ] Decay is triggered for input particle.");
      m_particleDecayHelper->decay(isp,isp.position(),isp.momentum(),isp.timeStamp()); // Registers TruthIncident internally
    }
 
    // validation mode - for all particle registered into stack
    if ( m_validationOutput && m_validationTool.isEnabled() ) {
      int endProcess = decayProc;
      m_validationTool->saveISFParticleInfo(isp,endProcess,&inputPar,timeLim.time,0.);
    }
 
    return nullptr;
  }
 
  // presample interactions if not done already
  Trk::PathLimit pathLim(-1.,0);
  const ISF::MaterialPathInfo* matLimit = isp.getUserInformation() ? isp.getUserInformation()->materialLimit() : nullptr;
  if (matLimit) {
    pathLim=Trk::PathLimit( matLimit->dMax,matLimit->process);
    pathLim.updateMat(matLimit->dCollected,13.,0.);          // arbitrary Z choice : update MaterialPathInfo
  } else if (absPdg!=999 && pHypothesis<99) { // need to resample
    pathLim = m_samplingTool->sampleProcess(m_randomEngine, isp.momentum().mag(),isp.charge(),pHypothesis);
  }
 
  // use extrapolation with path limit - automatic exit at subdetector boundary
  // additional exercise due to the current mismatch in geoID
  Trk::GeometrySignature nextGeoID=static_cast<Trk::GeometrySignature>(isp.nextGeoID());
 
  std::unique_ptr<const Trk::TrackParameters> eParameters = nullptr;
 
  // hit creation/energy deposit
  hitVector = (!m_hitsOff) ? new std::vector<Trk::HitInfo> : nullptr;
  if (hitVector) ATH_MSG_VERBOSE( "[ fatras transport ] hit creation activated");
 
  if ( !charged ) {
 
    eParameters = processor->transportNeutralsWithPathLimit(inputPar,
                                                            pathLim,
                                                            timeLim,
                                                            Trk::alongMomentum,
                                                            pHypothesis,
                                                            hitVector,
                                                            nextGeoID);
 
  } else {
 
    if (m_validationOutput && m_errorPropagation ) {
      // input covariance matrix
      AmgSymMatrix(5) inputCov;
      inputCov.setZero();
      std::unique_ptr<Trk::TrackParameters> measuredInputPar =
        inputPar.associatedSurface()
          .createUniqueTrackParameters(inputPar.parameters()[0],
                                       inputPar.parameters()[1],
                                       inputPar.parameters()[2],
                                       inputPar.parameters()[3],
                                       inputPar.parameters()[4],
                                       std::move(inputCov));
 
      eParameters = processor->extrapolateWithPathLimit(*measuredInputPar,
                                                        pathLim,
                                                        timeLim,
                                                        Trk::alongMomentum,
                                                        pHypothesis,
                                                        hitVector,
                                                        nextGeoID);
 
    } else {
 
      eParameters = processor->extrapolateWithPathLimit(inputPar,
                                                        pathLim,
                                                        timeLim,
                                                        Trk::alongMomentum,
                                                        pHypothesis,
                                                        hitVector,
                                                        nextGeoID);
    }
  }
 
  // create hits
  if (hitVector && !hitVector->empty()){
    ATH_MSG_VERBOSE( "[ fatras transport ] processing .");
    // create the hits from the parameters vector
    if ( isp.nextGeoID()==AtlasDetDescr::fAtlasID ) {
      if (!m_simHitCreatorID.empty()) m_simHitCreatorID->createHits(isp, *hitVector);
      ATH_MSG_VERBOSE( "[ fatras transport ] ID hits processed.");
    } else if ( isp.nextGeoID()==AtlasDetDescr::fUndefinedAtlasRegion ) {
      ATH_MSG_WARNING( "[ fatras transport ] hit vector without Geo info");
    } else if ( isp.nextGeoID()==AtlasDetDescr::fAtlasMS ) {
      if (!m_simHitCreatorMS.empty()) m_simHitCreatorMS->createHits(isp, *hitVector);
      ATH_MSG_VERBOSE( "[ fatras transport ] MS hits processed.");
    }
    // memory cleanup
    delete hitVector; hitVector = nullptr;
  }
 
  // additional exercise due to the current mismatch in geoID
  AtlasDetDescr::AtlasRegion geoID=AtlasDetDescr::AtlasRegion(5);
  if (nextGeoID<99) {
    geoID = static_cast<AtlasDetDescr::AtlasRegion>(nextGeoID);
  }
 
  // validation mode - for all particle registered into stack
  if ( m_validationOutput && m_validationTool.isEnabled() ) {
 
    int dProc = ( timeLim.tMax>0. && timeLim.tMax<=timeLim.time ) ? timeLim.process : 0;
    int mProc = ( pathLim.x0Max>0. && ( pathLim.x0Max <= pathLim.x0Collected ||
                 (pathLim.process>100 && pathLim.x0Max <= pathLim.l0Collected))) ? pathLim.process : 0;
    int endProcess = eParameters ? 0 : ( dProc > mProc ? dProc : mProc );
 
    m_validationTool->saveISFParticleInfo(isp,endProcess,eParameters.get(),timeLim.time,pathLim.x0Collected);
 
  }
 
  // create the updated particle at the end of processing step
  ISF::ISFParticle* uisp = eParameters ? m_iparticleHelper->updatedParticle(isp,
                                        eParameters->position(),
                                        eParameters->momentum(),
                                        timeLim.time-isp.timeStamp()) : nullptr;     // update expects time difference
  // free memory
  if ( hitVector ) {
    delete hitVector;
  }
 
  if (uisp && m_validationOutput) {
    // save validation info
    ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
    validInfo->setProcess(0);
    if (isp.getUserInformation()) validInfo->setGeneration(isp.getUserInformation()->generation());
    else validInfo->setGeneration(0);     // assume primary parent
    uisp->setUserInformation(validInfo);
  }
 
  // TODO save material ? if (uisp && m_validationOutput) uisp->setX0(pathLim.x0Collected,m_eloss);
 
  // trigger decay or return to stack
 
  if (uisp && timeLim.tMax>0.  && timeLim.time >=timeLim.tMax ) {
    if (!m_particleDecayHelper.empty()) {
      ATH_MSG_VERBOSE( "[ fatras transport ] Decay is triggered for input particle.");
      m_particleDecayHelper->decay(*uisp,uisp->position(),uisp->momentum(),uisp->timeStamp()); // Registers TruthIncident internally
    }
    delete uisp;
    return nullptr;
  } else if (uisp) {
    ATH_MSG_VERBOSE( "[ fatras transport ] Particle transported to detector boundary, return to stack, timing:" << uisp->timeStamp());
    uisp->setNextGeoID(geoID);
    return uisp;
  }
 
  // nothing to be returned to the stack (and we should not get here at all)
  return nullptr;
}

◆ retrieveTool()

template<class T>

StatusCode iFatras::TransportTool::retrieveTool ( ToolHandle< T > & thandle )

inlineprivate

templated Tool retrieval - gives unique handling & look and feel

Definition at line 76 of file ISF_Fatras/ISF_FatrasTools/src/TransportTool.h.

                                                                       {
       if (!thandle.empty() && thandle.retrieve().isFailure()){
         ATH_MSG_FATAL( "[ fatras setup ] Cannot retrieve " << thandle << ". Abort.");
         return StatusCode::FAILURE;
       }
       else {
         ATH_MSG_DEBUG("[ fatras setup ] Successfully retrieved " << thandle);
       }
       return StatusCode::SUCCESS;
     }