diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/HadIntProcessorParametric.h b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/HadIntProcessorParametric.h
index 763bd4eeae23eac05bc0e4a14b273b476ecbbdc2..c00aba52724335d6f93a6ab6ae71bd5f6c2ada0b 100644
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/HadIntProcessorParametric.h
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/HadIntProcessorParametric.h
@@ -75,24 +75,26 @@ namespace iFatras {
StatusCode finalize();
/** interface for processing of the nuclear interactions */
- bool hadronicInteraction(const Trk::TrackParameters& parm, double p, double E,
+ bool hadronicInteraction(const Amg::Vector3D& position, const Amg::Vector3D& momentum,
+ double p, double E, double charge,
const Trk::MaterialProperties& mprop, double pathCorrection,
Trk::ParticleHypothesis particle=Trk::pion) const;
/** interface for processing of the presampled nuclear interaction */
bool recordHadState(double time, double p,
- const Amg::Vector3D& vertex,
- const Amg::Vector3D& particleDir,
- Trk::ParticleHypothesis particle ) const;
+ const Amg::Vector3D& vertex,
+ const Amg::Vector3D& particleDir,
+ Trk::ParticleHypothesis particle ) const;
- bool doHadronicInteraction(double time,const Trk::TrackParameters& parm,
- const Trk::Material* ematprop,
+ bool doHadronicInteraction(double time, const Amg::Vector3D& position, const Amg::Vector3D& momentum,
+ const Trk::Material* emat,
Trk::ParticleHypothesis particle,
bool processSecondaries) const;
- ISF::ISFParticleVector doHadIntOnLayer(const ISF::ISFParticle* parent, double time,const Trk::TrackParameters& parm,
- const Trk::MaterialProperties* ematprop,
- Trk::ParticleHypothesis particle) const;
+ ISF::ISFParticleVector doHadIntOnLayer(const ISF::ISFParticle* parent, double time,
+ const Amg::Vector3D& position, const Amg::Vector3D& momentum,
+ const Trk::Material* emat,
+ Trk::ParticleHypothesis particle) const;
private:
/** interface for calculation of absorption length */
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McBetheHeitlerEnergyLossUpdator.h b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McBetheHeitlerEnergyLossUpdator.h
deleted file mode 100755
index 3168e672184b1a2b8b3a9a22765b20347cd913fd..0000000000000000000000000000000000000000
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McBetheHeitlerEnergyLossUpdator.h
+++ /dev/null
@@ -1,95 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-/*********************************************************************************
- McBetheHeitlerEnergyLossUpdator.h - description
- -------------------------------------------------
-begin : Wednesday 28th June 2006
-author : atkinson
-email : Tom.Atkinson@cern.ch
-decription : Class for updating energy loss based on the Bethe-Heitler
- distribution
-*********************************************************************************/
-
-#ifndef ISF_Fatras_McBetheHeitlerEnergyLossUpdator_H
-#define ISF_Fatras_McBetheHeitlerEnergyLossUpdator_H
-
-// Gaudi & Athena
-#include "AthenaBaseComps/AthAlgTool.h"
-#include "GaudiKernel/ServiceHandle.h"
-#include "AthenaKernel/IAtRndmGenSvc.h"
-// Trk
-#include "TrkExInterfaces/IEnergyLossUpdator.h"
-#include "TrkEventPrimitives/ParticleHypothesis.h"
-#include "TrkEventPrimitives/PropDirection.h"
-
-namespace Trk{
- class MaterialProperties;
- class EnergyLoss;
-}
-
-namespace iFatras{
-
-class McBetheHeitlerEnergyLossUpdator : public AthAlgTool, virtual public Trk::IEnergyLossUpdator {
-
- public:
-
- /** Constructor with AlgTool parameters */
- McBetheHeitlerEnergyLossUpdator( const std::string&, const std::string&, const IInterface* );
-
- /** Destructor */
- ~McBetheHeitlerEnergyLossUpdator();
-
- /** AlgTool initialise method */
- StatusCode initialize();
-
- /** AlgTool finalize method */
- StatusCode finalize();
-
- /** dEdx method - not used */
- double dEdX( const Trk::MaterialProperties& materialProperties,
- double momentum,
- Trk::ParticleHypothesis particleHypothesis = Trk::electron ) const;
-
- /** Method to determine the change in q/p and the variance */
- Trk::EnergyLoss* energyLoss( const Trk::MaterialProperties& materialProperties,
- double momentum,
- double pathCorrection,
- Trk::PropDirection direction = Trk::alongMomentum,
- Trk::ParticleHypothesis particleHypothesis = Trk::electron,
- bool mpv = false) const;
-
- /** Method to recalculate Eloss values for the fit setting an elossFlag using as an input
- the detailed Eloss information Calorimeter energy, error momentum and momentum error */
-
- Trk::EnergyLoss* updateEnergyLoss( Trk::EnergyLoss*, double, double, double, double, int&) const override { return 0; }
-
- /** Routine to calculate X0 and Eloss scale factors for the Calorimeter and Muon System */
- void getX0ElossScales(int, double, double, double&, double& ) const override {}
-
- private:
-
- /** Private method to compute the Bethe-Heitler PDF */
- std::vector<double> betheHeitlerPDF( double pathLength ) const;
-
- ServiceHandle<IAtRndmGenSvc> m_rndmGenSvc; //!< Random number generator
- CLHEP::HepRandomEngine* m_randomEngine;
- std::string m_randomEngineName; //!< Name of the random number stream
-
- mutable std::vector<double> m_pdf;
- mutable double m_cashedPathLength;
- int m_numberPointsInPDF;
-
- double m_scaleFactor; //!< the one free parameter to scale
-
- static Trk::ParticleMasses s_particleMasses; //!< struct of Particle masses
-
-};
-
-} // end iFatras namespace
-
-inline double iFatras::McBetheHeitlerEnergyLossUpdator::dEdX( const Trk::MaterialProperties&, double, Trk::ParticleHypothesis ) const
-{ return 0.; }
-
-#endif
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McEnergyLossUpdator.h b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McEnergyLossUpdator.h
index 446f7d2a539570d4854bd1258cc92755d5c02dac..c590e605e1803b0244395c10f08e5334064cd136 100755
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McEnergyLossUpdator.h
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McEnergyLossUpdator.h
@@ -67,22 +67,20 @@ namespace iFatras{
Trk::ParticleHypothesis particleHypothesis = Trk::pion,
bool mpv=true) const override;
- /** Method to recalculate Eloss values for the fit setting an elossFlag using as an input
- the detailed Eloss information Calorimeter energy, error momentum and momentum error */
-
+ /** Dummy methodes imposed by public interface - cleanup */
+ /** Method to recalculate Eloss values for the fit setting an elossFlag using as an input
+ the detailed Eloss information Calorimeter energy, error momentum and momentum error */
Trk::EnergyLoss* updateEnergyLoss( Trk::EnergyLoss*, double, double, double, double, int&) const override { return 0; }
- /** Routine to calculate X0 and Eloss scale factors for the Calorimeter and Muon System */
- void getX0ElossScales(int, double, double, double&, double& ) const override {}
-
+ /** Routine to calculate X0 and Eloss scale factors for the Calorimeter and Muon System */
+ void getX0ElossScales(int, double, double, double&, double& ) const override {}
+ /** Dummy methods end here */
+
private:
ToolHandle<IEnergyLossUpdator> m_energyLossUpdator; //!< Pointer to the energy loss updator
int m_energyLossDistribution; //!< include energy loss straggling or not ( 0 == none, 1 == gauss, 2 == landau)
- bool m_dedicatedElectronUpdator; //!< boolean switch for use of a dedicated eloss updator
- ToolHandle<IEnergyLossUpdator> m_elEnergyLossUpdator; //!< Pointer to the energy loss updator - electrons
-
/** Random Generator service */
ServiceHandle<IAtRndmGenSvc> m_rndGenSvc;
/** Random engine */
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsEngine.h b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsEngine.h
index 3dd8158fdef5a91aef69515f9224239966d93d74..1c73288f95ce74a98c2d0223b3b17aa49a0f9532 100644
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsEngine.h
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsEngine.h
@@ -2,6 +2,7 @@
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
+
///////////////////////////////////////////////////////////////////
// McMaterialEffectsEngine.h, (c) ATLAS Detector software
///////////////////////////////////////////////////////////////////
@@ -17,11 +18,12 @@
// Trk
#include "TrkExInterfaces/IMaterialEffectsEngine.h"
+#include "TrkExInterfaces/IEnergyLossUpdator.h"
#include "TrkExUtils/ExtrapolationCell.h"
-#include "TrkExUtils/MaterialInteraction.h"
#include "TrkEventPrimitives/PropDirection.h"
#include "TrkParameters/TrackParameters.h"
#include "TrkExUtils/MaterialUpdateMode.h"
+#include "TrkEventPrimitives/PdgToParticleHypothesis.h"
// ISF
#include "ISF_Event/ITruthIncident.h"
@@ -32,129 +34,190 @@
class StoreGateSvc;
-namespace Trk {
+namespace Trk{
class EnergyLoss;
class CylinderVolumeBounds;
- class IEnergyLossUpdator;
- class IMultipleScatteringUpdator;
class TrackingGeometry;
class ITrackingGeometrySvc;
}
namespace ISF {
- class IParticleBroker;
- class ISFParticle;
- class ITruthSvc;
+ class IParticleBroker;
+ class ISFParticle;
+ class ITruthSvc;
}
namespace iFatras {
+
+ class IPhysicsValidationTool;
+ class IEnergyLossSampler;
+ class IMultipleScatteringSampler;
+ class IHadronicInteractionProcessor;
+ class IPhotonConversionSampler;
+ class IParticleDecayHelper;
+ class IProcessSamplingTool;
/** @class McMaterialEffectsEngine
-
- Updator for a track on a Trk::Layer, it extends the IMaterialEffectsEngine to
- be used inside the Extrapolator with an update based on a Random number.
-
- The McMaterialEffectsEngine uses both an extended EnergyLossUpdator
- and the standard ATLAS MultipleScatteringUpdator configured for the Gaussian-Mixture-Model.
-
- @author Andreas.Salzburger@cern.ch, Noemi.Calace@cern.ch
- */
-
- class McMaterialEffectsEngine : public AthAlgTool, virtual public Trk::IMaterialEffectsEngine {
- public:
- /**AlgTool constructor for McMaterialEffectsEngine*/
- McMaterialEffectsEngine(const std::string&,const std::string&,const IInterface*);
- /**Destructor*/
- virtual ~McMaterialEffectsEngine();
-
- /** AlgTool initailize method.*/
- StatusCode initialize();
- /** AlgTool finalize method */
- StatusCode finalize();
-
- /** Updator interface (full update for a layer):
- given track parameters are deleted internally,
- if no update has to be done,
- the pointer is returned
- */
- virtual Trk::ExtrapolationCode handleMaterial(Trk::ExCellCharged& ecCharged,
- Trk::PropDirection dir=Trk::alongMomentum,
- Trk::MaterialUpdateStage matupstage=Trk::fullUpdate) const;
-
- virtual Trk::ExtrapolationCode handleMaterial(Trk::ExCellNeutral&,
- Trk::PropDirection,
- Trk::MaterialUpdateStage) const {return Trk::ExtrapolationCode::InProgress; }
-
- private:
-
- const Trk::TrackParameters* updateTrackParameters(const Trk::TrackParameters& parameters,
- Trk::ExCellCharged& eCell,
- Trk::PropDirection dir,
- Trk::MaterialUpdateStage matupstage) const;
-
- double simulate_theta_space(double beta, double p,double dOverX0,double Z, double scale) const;
- double * get_gaussian(double beta, double p,double dOverX0, double scale) const;
- double * get_gaussmix(double beta, double p,double dOverX0,double Z, double scale) const;
- double * get_semigauss(double beta,double p,double dOverX0,double Z, double scale) const;
- double sim_gaussmix(double * scattering_params) const;
- double sim_semigauss(double * scattering_params) const;
-
- void ionize(const Trk::TrackParameters& parm, AmgVector(5)& updatedPar, double dInX0,
- Trk::PropDirection pdir, Trk::ParticleHypothesis particle ) const ;
-
- /** the private multiple scattering sigma calculation*/
- double msSigma(double dInX0, double p,Trk::ParticleHypothesis particle) const;
-
- /** the private multiple Scattering update method, thetaMs is the projected random number*/
- void multipleScatteringUpdate( const Trk::TrackParameters& parm,
- AmgVector(5)& parameters,
- double sigmaMSproj, double pathCorrection) const;
-
- /** IEnergyLossUpdator */
- bool m_eLoss;
- ToolHandle<Trk::IEnergyLossUpdator> m_eLossUpdator;
- /** IMultipleScatteringUpdator */
- bool m_ms;
- ToolHandle<Trk::IMultipleScatteringUpdator> m_msUpdator;
+ Mc Material effects engine for charged and neutral (fast track simulation) used inside the ExtrapolatorEngine. It extends the IMaterialEffectsEngine.
- /** Minimum momentum cut */
- double m_minimumMomentum;
+ The McMaterialEffectsEngine uses both an extended EnergyLossSampler
+ and the standard ATLAS MultipleScatteringSampler configured for the Gaussian-Mixture-Model.
- /** switch between MSUpdator and local parametrization */
- bool m_use_msUpdator;
-
- /** Switch to use reference material */
- bool m_referenceMaterial;
-
- /** Switch to use bending correction */
- bool m_bendingCorrection;
-
- /** Random Generator service */
- ServiceHandle<IAtRndmGenSvc> m_rndGenSvc;
- /** Random engine */
- CLHEP::HepRandomEngine* m_randomEngine;
- std::string m_randomEngineName; //!< Name of the random number stream
-
- mutable const Trk::TrackingGeometry* m_trackingGeometry; //!< the tracking geometry owned by the navigator
- ServiceHandle<Trk::ITrackingGeometrySvc> m_trackingGeometrySvc; //!< ToolHandle to the TrackingGeometrySvc
- std::string m_trackingGeometryName; //!< default name of the TrackingGeometry
-
- /** struct of Particle Masses */
- Trk::ParticleMasses m_particleMasses;
-
- ServiceHandle<ISF::IParticleBroker> m_particleBroker;
- ServiceHandle<ISF::ITruthSvc> m_truthRecordSvc;
-
- /** cache incoming particle */
- mutable const ISF::ISFParticle* m_isp;
-
- /** cache layer properties */
- mutable const Trk::Layer* m_layer;
- mutable const Trk::MaterialProperties* m_matProp;
-
- };
+ @author Andreas.Salzburger@cern.ch, Noemi.Calace@cern.ch
+ */
+ class McMaterialEffectsEngine : public AthAlgTool, virtual public Trk::IMaterialEffectsEngine {
+ public:
+ /**AlgTool constructor for McMaterialEffectsEngine*/
+ McMaterialEffectsEngine(const std::string&,const std::string&,const IInterface*);
+ /**Destructor*/
+ virtual ~McMaterialEffectsEngine();
+
+ /** AlgTool initailize method.*/
+ StatusCode initialize();
+ /** AlgTool finalize method */
+ StatusCode finalize();
+
+ /** charged extrapolation */
+ virtual Trk::ExtrapolationCode handleMaterial(Trk::ExCellCharged& ecCharged,
+ Trk::PropDirection dir=Trk::alongMomentum,
+ Trk::MaterialUpdateStage matupstage=Trk::fullUpdate) const;
+
+ /** charged extrapolation */
+ virtual Trk::ExtrapolationCode handleMaterial(Trk::ExCellNeutral& ecNeutral,
+ Trk::PropDirection dir= Trk::alongMomentum,
+ Trk::MaterialUpdateStage matupstage=Trk::fullUpdate) const;
+
+ private:
+
+ template <class T> Trk::ExtrapolationCode handleMaterialT( Trk::ExtrapolationCell<T>& eCell,
+ Trk::PropDirection dir=Trk::alongMomentum,
+ Trk::MaterialUpdateStage matupstage=Trk::fullUpdate) const;
+
+ Trk::ExtrapolationCode processMaterialOnLayer(const ISF::ISFParticle* isp,
+ Trk::ExCellCharged& eCell,
+ Trk::PropDirection dir,
+ float& mFraction) const;
+
+ Trk::ExtrapolationCode processMaterialOnLayer(const ISF::ISFParticle* isp,
+ Trk::ExCellNeutral& eCell,
+ Trk::PropDirection dir,
+ float& mFraction) const;
+
+ void multipleScatteringUpdate(const Trk::TrackParameters& pars,
+ AmgVector(5)& parameters,
+ double simTheta,
+ double num_deltaPhi) const;
+
+ void recordBremPhoton(const ISF::ISFParticle* parent,
+ double time,
+ double pElectron,
+ double gammaE,
+ const Amg::Vector3D& vertex,
+ Amg::Vector3D& particleDir,Trk::ParticleHypothesis ) const;
+
+ void recordBremPhotonLay(const ISF::ISFParticle* parent,
+ double time,
+ double pElectron,
+ double gammaE,
+ const Amg::Vector3D& vertex,
+ Amg::Vector3D& particleDir,Trk::ParticleHypothesis,
+ Trk::PropDirection dir,float mFraction ) const;
+
+ ISF::ISFParticle* bremPhoton(const ISF::ISFParticle* parent,
+ double time,
+ double pElectron,
+ double gammaE,
+ const Amg::Vector3D& vertex,
+ Amg::Vector3D& particleDir, Trk::ParticleHypothesis ) const;
+
+ void radiate(const ISF::ISFParticle* parent, AmgVector(5)& parm ,
+ Trk::ExCellCharged& eCell, float pathlim, float mFr,Trk::PropDirection dir,float refX) const;
+
+ /** IProcessSamplig */
+ ToolHandle<IProcessSamplingTool> m_samplingTool;
+
+ /** IEnergyLossSampler */
+ bool m_eLoss;
+ //ToolHandle<iFatras::IEnergyLossSampler> m_eLossSampler;
+ ToolHandle<Trk::IEnergyLossUpdator> m_eLossSampler;
+ /** Boolean switch for use of a dedicated eloss updator */
+ bool m_dedicatedElectronSampler;
+ /** Pointer to the energy loss sampler - electrons */
+ ToolHandle<iFatras::IEnergyLossSampler> m_elEnergyLossSampler;
+
+ /** IMultipleScatteringSampler */
+ bool m_ms;
+ ToolHandle<iFatras::IMultipleScatteringSampler> m_msSampler;
+
+ /** Particle Decay */
+ ToolHandle<IParticleDecayHelper> m_particleDecayHelper;
+
+ /** describe deflection parametric/do real deflection */
+ bool m_parametricScattering;
+
+ /** use the relativistic hertz dipole for brem photon radiation */
+ bool m_uniformHertzDipoleAngle;
+
+ /** MCTruth process code for TruthIncidents created by this tool */
+ Barcode::PhysicsProcessCode m_processCode;
+
+ /** Minimum momentum cut */
+ double m_minimumMomentum;
+
+ /** Minimum momentum for brem photons */
+ double m_minimumBremPhotonMomentum;
+
+ /** Create brem photons flag */
+ bool m_createBremPhoton;
+
+ /** Switch to use bending correction */
+ bool m_bendingCorrection;
+
+ /** Random Generator service */
+ ServiceHandle<IAtRndmGenSvc> m_rndGenSvc;
+ /** Random engine */
+ CLHEP::HepRandomEngine* m_randomEngine;
+ /** Name of the random number stream */
+ std::string m_randomEngineName;
+
+ /** the tracking geometry owned by the navigator */
+ mutable const Trk::TrackingGeometry* m_trackingGeometry;
+ /** ToolHandle to the TrackingGeometrySvc */
+ ServiceHandle<Trk::ITrackingGeometrySvc> m_trackingGeometrySvc;
+ /** default name of the TrackingGeometry */
+ std::string m_trackingGeometryName;
+
+ /** struct of Particle Masses */
+ Trk::ParticleMasses m_particleMasses;
+ Trk::PdgToParticleHypothesis m_pdgToParticleHypothesis;
+
+ /** Validation section */
+ bool m_validationMode;
+ ToolHandle<IPhysicsValidationTool> m_validationTool;
+
+ ServiceHandle<ISF::IParticleBroker> m_particleBroker;
+ ServiceHandle<ISF::ITruthSvc> m_truthRecordSvc;
+
+ /** useful for the angle calculation of the brem photon */
+ double m_oneOverThree;
+
+ /** cache incoming particle */
+ mutable const ISF::ISFParticle* m_isp;
+
+ /** cache layer properties */
+ mutable const Trk::Layer* m_layer;
+ mutable const Trk::MaterialProperties* m_matProp;
+ mutable float m_thicknessInX0;
+ mutable float m_thicknessInL0;
+
+ double m_projectionFactor;
+ };
+
}
+//!< define the templated function
+#include "McMaterialEffectsEngine.icc"
+
#endif // ISF_FATRASTOOLS_MCMATERIALEFFECTENIGINE_H
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsEngine.icc b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsEngine.icc
new file mode 100644
index 0000000000000000000000000000000000000000..2eea3be465e631fb54bede7f43334aabbd33723e
--- /dev/null
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsEngine.icc
@@ -0,0 +1,50 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// McMaterialEffectsEngine.icc, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#include "TrkExInterfaces/IMaterialEffectsEngine.h"
+#include "ISF_Event/ParticleClipboard.h"
+#include "TrkGeometry/Layer.h"
+#include <iostream>
+
+template <class T> Trk::ExtrapolationCode iFatras::McMaterialEffectsEngine::handleMaterialT( Trk::ExtrapolationCell<T>& eCell,
+ Trk::PropDirection dir,
+ Trk::MaterialUpdateStage matupstage) const
+{
+ // get parent particle
+ const ISF::ISFParticle *isp = ISF::ParticleClipboard::getInstance().getParticle();
+ // something is seriously wrong if there is no parent particle
+ assert(isp);
+ m_isp = isp;
+
+ // the Extrapolator made sure that the layer is the lead layer && the parameters are the lead parameters
+ m_layer = eCell.leadLayer;
+ if (m_layer && eCell.leadParameters) {
+ // path correction
+ double pathCorrection = m_layer->surfaceRepresentation().pathCorrection(eCell.leadParameters->position(),dir*(eCell.leadParameters->momentum()));
+ // the relative direction wrt with the layer
+ Trk::PropDirection rlDir = (pathCorrection > 0. ? Trk::alongMomentum : Trk::oppositeMomentum);
+ // multiply by the pre-and post-update factor
+ double mFactor = m_layer->layerMaterialProperties()->factor(rlDir, matupstage);
+ //EX_MSG_DEBUG(eCell.navigationStep, "Update", "char", "Update on layer with index " << m_layer->layerIndex().value() << " - corr factor = " << pathCorrection*mFactor);
+
+ const Trk::MaterialProperties* materialProperties = m_layer->layerMaterialProperties()->fullMaterial(eCell.leadParameters->position());
+ if ( materialProperties ) {
+ m_matProp = materialProperties;
+ m_thicknessInX0 = mFactor*m_matProp->thicknessInX0();
+ m_thicknessInL0 = mFactor*m_matProp->thicknessInL0();
+
+ float mFraction=0.;
+ return processMaterialOnLayer(isp,eCell,dir,mFraction);
+ }
+ }
+ return Trk::ExtrapolationCode::InProgress;
+}
+
+
+
+
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsUpdator.h b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsUpdator.h
index aa89a366e20e8a506e284e2c88d71ea2166daba9..0e3ac799713abeee0bd45e21744239c5507e5a9f 100755
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsUpdator.h
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/McMaterialEffectsUpdator.h
@@ -61,7 +61,7 @@ namespace iFatras {
class IPhysicsValidationTool;
class IHadronicInteractionProcessor;
class IPhotonConversionTool;
- class IParticleDecayer;
+ class IParticleDecayHelper;
/** @class McMaterialEffectsUpdator
@@ -217,7 +217,7 @@ namespace iFatras {
ToolHandle<IHadronicInteractionProcessor> m_hadIntProcessor;
/** Particle Decay */
- ToolHandle<IParticleDecayer> m_particleDecayer;
+ ToolHandle<IParticleDecayHelper> m_particleDecayer;
/** Minimum momentum cut */
double m_minimumMomentum;
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/MultipleScatteringSamplerGaussianMixture.h b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/MultipleScatteringSamplerGaussianMixture.h
new file mode 100644
index 0000000000000000000000000000000000000000..a7032a4953d03164a3372bc51887f6b936c24614
--- /dev/null
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/MultipleScatteringSamplerGaussianMixture.h
@@ -0,0 +1,105 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// MultipleScatteringSamplerGaussianMixture.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#ifndef ISF_FATRASTOOLS_MULTIPLESCATTERINGSAMPLERGAUSSIANMIXTURE_H
+#define ISF_FATRASTOOLS_MULTIPLESCATTERINGSAMPLERGAUSSIANMIXTURE_H
+
+// Gaudi
+#include "AthenaBaseComps/AthAlgTool.h"
+// Trk
+#include "ISF_FatrasInterfaces/IMultipleScatteringSampler.h"
+#include "TrkEventPrimitives/PropDirection.h"
+#include "TrkEventPrimitives/ParticleHypothesis.h"
+#include "GaudiKernel/ServiceHandle.h"
+
+#include "AthenaKernel/IAtRndmGenSvc.h"
+
+#include "CLHEP/Random/RandFlat.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+#include "EventPrimitives/EventPrimitives.h"
+
+namespace Trk {
+ class MaterialProperties;
+}
+
+namespace iFatras {
+
+ /**@class MultipleScatteringSamplerGaussianMixture
+
+ @author Noemi.Calace@cern.ch , Andreas.Salzburger@cern.ch
+ */
+
+ class MultipleScatteringSamplerGaussianMixture : public AthAlgTool,
+ virtual public IMultipleScatteringSampler {
+
+ public:
+ /** AlgTool like constructor */
+ MultipleScatteringSamplerGaussianMixture(const std::string&,const std::string&,const IInterface*);
+
+ /**Virtual destructor*/
+ virtual ~MultipleScatteringSamplerGaussianMixture();
+
+ /** AlgTool initailize method.*/
+ StatusCode initialize();
+
+ /** AlgTool finalize method */
+ StatusCode finalize();
+
+ /** Calculate the sigma on theta introduced by multiple scattering,
+ according to the RutherFord-Scott Formula
+ */
+ double simTheta(const Trk::MaterialProperties& mat,
+ double p,
+ double pathcorrection,
+ Trk::ParticleHypothesis particle=Trk::pion,
+ double deltaE=0.) const;
+
+ private:
+ bool m_log_include; //!< boolean switch to include log term
+
+ bool m_optGaussianMixtureG4; //!< modifies the Fruehwirth/Regler model to fit with G4
+
+ //========== used for Gaussian mixture model =================================================
+
+ /** Random Generator service */
+ ServiceHandle<IAtRndmGenSvc> m_rndGenSvc;
+ /** Random engine */
+ CLHEP::HepRandomEngine* m_randomEngine;
+ std::string m_randomEngineName; //!< Name of the random number stream
+
+ static Trk::ParticleMasses s_particleMasses; //!< struct of Particle masses
+
+ static double s_main_RutherfordScott; //!< main factor of Rutherford-Scott formula
+ static double s_log_RutherfordScott; //!< log factor of Rutherford-Scott formula
+
+ static double s_main_RossiGreisen; //!< main factor for Rossi-Greisen formula
+ static double s_log_RossiGreisen; //!< main factor for Rossi-Greisen formula
+
+
+ // ========= Gaussian mixture model Fruehwirth, Regler Nucl. Inst. Methods A 456(2001) =========
+ static double s_gausMixSigma1_a0; //!< Gaussian mixture model: Sigma parameter a0
+ static double s_gausMixSigma1_a1; //!< Gaussian mixture model: Sigma parameter a1
+ static double s_gausMixSigma1_a2; //!< Gaussian mixture model: Sigma parameter a2
+
+ static double s_gausMixEpsilon_a0; //!< Gaussian mixture model: Epsilon parameter a0
+ static double s_gausMixEpsilon_a1; //!< Gaussian mixture model: Epsilon parameter a1
+ static double s_gausMixEpsilon_a2; //!< Gaussian mixture model: Epsilon parameter a2
+
+ static double s_gausMixEpsilon_b0; //!< Gaussian mixture model: Epsilon parameter b0
+ static double s_gausMixEpsilon_b1; //!< Gaussian mixture model: Epsilon parameter b1
+ static double s_gausMixEpsilon_b2; //!< Gaussian mixture model: Epsilon parameter b2
+
+ static double s_projectionFactor; //!< projection factor to scale the projected angle out of the plane
+
+ };
+
+
+} // end of namespace
+
+
+#endif // ISF_FATRASTOOLS_MULTIPLESCATTERINGSAMPLERGAUSSIANMIXTURE_H
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/MultipleScatteringSamplerGeneralMixture.h b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/MultipleScatteringSamplerGeneralMixture.h
new file mode 100644
index 0000000000000000000000000000000000000000..113ee363c12e9e18634fc3ebce6cfc51b0ca6c78
--- /dev/null
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/MultipleScatteringSamplerGeneralMixture.h
@@ -0,0 +1,99 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// MultipleScatteringSamplerGeneralMixture.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#ifndef ISF_FATRASTOOLS_MULTIPLESCATTERINGSAMPLERGENERALMIXTURE_H
+#define ISF_FATRASTOOLS_MULTIPLESCATTERINGSAMPLERGENERALMIXTURE_H
+
+// Gaudi
+#include "AthenaBaseComps/AthAlgTool.h"
+// Trk
+#include "ISF_FatrasInterfaces/IMultipleScatteringSampler.h"
+#include "TrkEventPrimitives/PropDirection.h"
+#include "TrkEventPrimitives/ParticleHypothesis.h"
+#include "GaudiKernel/ServiceHandle.h"
+
+#include "AthenaKernel/IAtRndmGenSvc.h"
+
+#include "CLHEP/Random/RandFlat.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+#include "EventPrimitives/EventPrimitives.h"
+
+namespace Trk {
+ class MaterialProperties;
+}
+
+namespace iFatras {
+
+ /**@class MultipleScatteringSamplerGeneralMixture
+
+ @author Noemi.Calace@cern.ch , Andreas.Salzburger@cern.ch, Artem.Basalaev@cern.ch
+ */
+
+ class MultipleScatteringSamplerGeneralMixture : public AthAlgTool,
+ virtual public IMultipleScatteringSampler {
+
+ public:
+ /** AlgTool like constructor */
+ MultipleScatteringSamplerGeneralMixture(const std::string&,const std::string&,const IInterface*);
+
+ /**Virtual destructor*/
+ virtual ~MultipleScatteringSamplerGeneralMixture();
+
+ /** AlgTool initailize method.*/
+ StatusCode initialize();
+
+ /** AlgTool finalize method */
+ StatusCode finalize();
+
+ /** Calculate theta introduced by multiple scattering,
+ according to the RutherFord-Scott Formula
+ */
+ double simTheta(const Trk::MaterialProperties& mat,
+ double p,
+ double pathcorrection,
+ Trk::ParticleHypothesis particle=Trk::pion,
+ double deltaE=0.) const;
+
+ private:
+ bool m_log_include; //!< boolean switch to include log term
+ static double s_projectionFactor; //!< projection factor to scale the projected angle out of the plane
+
+ //========== used for General mixture model =================================================
+
+ /** Random Generator service */
+ ServiceHandle<IAtRndmGenSvc> m_rndGenSvc;
+ /** Random engine */
+ CLHEP::HepRandomEngine* m_randomEngine;
+ std::string m_randomEngineName; //!< Name of the random number stream
+
+ static Trk::ParticleMasses s_particleMasses; //!< struct of Particle masses
+
+ static double s_main_RossiGreisen; //!< main factor for Rossi-Greisen formula
+ static double s_log_RossiGreisen; //!< main factor for Rossi-Greisen formula
+
+
+ // ========= General mixture model R.Frühwirth, M. Liendl. Comp. Phys. Comm. 141 (2001) 230–246 =========
+ static double s_genMixScale; //!< General mixture model: Scaling factor
+
+ //!< General mixture model: get parameters for single gaussian simulation
+ double *getGaussian(double beta, double p,double dOverX0, double scale) const;
+ //!< General mixture model: get parameters for gaussian mixture
+ double *getGaussmix(double beta, double p,double dOverX0,double Z, double scale) const;
+ //!< General mixture model: get parameters for semi-gaussian mixture
+ double *getSemigauss(double beta,double p,double dOverX0,double Z, double scale) const;
+ //!< General mixture model: simulate semi-gaussian mixture
+ double simGaussmix(double * scattering_params) const;
+ //!< General mixture model: simulate gaussian mixture
+ double simSemigauss(double * scattering_params) const;
+ };
+
+
+} // end of namespace
+
+
+#endif // ISF_FATRASTOOLS_MULTIPLESCATTERINGSAMPLERGENERALMIXTURE_H
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/MultipleScatteringSamplerHighland.h b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/MultipleScatteringSamplerHighland.h
new file mode 100644
index 0000000000000000000000000000000000000000..1ede9ea15b95700ef1b135e3285bcea28999a912
--- /dev/null
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/MultipleScatteringSamplerHighland.h
@@ -0,0 +1,95 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// MultipleScatteringSamplerHighland.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#ifndef ISF_FATRASTOOLS_MULTIPLESCATTERINGSAMPLERHIGHLAND_H
+#define ISF_FATRASTOOLS_MULTIPLESCATTERINGSAMPLERHIGHLAND_H
+
+// Gaudi
+#include "AthenaBaseComps/AthAlgTool.h"
+// Trk
+#include "ISF_FatrasInterfaces/IMultipleScatteringSampler.h"
+#include "TrkEventPrimitives/PropDirection.h"
+#include "TrkEventPrimitives/ParticleHypothesis.h"
+#include "GaudiKernel/ServiceHandle.h"
+
+#include "AthenaKernel/IAtRndmGenSvc.h"
+
+#include "CLHEP/Random/RandFlat.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+#include "EventPrimitives/EventPrimitives.h"
+
+namespace Trk {
+ class MaterialProperties;
+}
+
+namespace iFatras {
+
+ /**@class MultipleScatteringSamplerHighland
+
+ The Formula used is the highland formula for the projected scattering angle :
+
+ @f$ \theta_{ms} = \frac{13.6MeV}{p}\cdot\sqrt{t/X_{0}}[1 + 0.038\ln(t/X_{0})] @f$
+
+ What is returned is the square of the expectation value of the deflection
+ @f$ < (\theta_ms)^2 > = \sigma_ms^2 @f$
+
+ @author Noemi.Calace@cern.ch , Andreas.Salzburger@cern.ch
+ */
+
+ class MultipleScatteringSamplerHighland : public AthAlgTool,
+ virtual public IMultipleScatteringSampler {
+
+ public:
+ /** AlgTool like constructor */
+ MultipleScatteringSamplerHighland(const std::string&,const std::string&,const IInterface*);
+
+ /**Virtual destructor*/
+ virtual ~MultipleScatteringSamplerHighland();
+
+ /** AlgTool initailize method.*/
+ StatusCode initialize();
+
+ /** AlgTool finalize method */
+ StatusCode finalize();
+
+ /** Calculate the theta introduced by multiple scattering,
+ according to the RutherFord-Scott Formula
+ */
+ double simTheta(const Trk::MaterialProperties& mat,
+ double p,
+ double pathcorrection,
+ Trk::ParticleHypothesis particle=Trk::pion,
+ double deltaE=0.) const;
+
+ private:
+ bool m_log_include; //!< boolean switch to include log term
+
+ /** Random Generator service */
+ ServiceHandle<IAtRndmGenSvc> m_rndGenSvc;
+ /** Random engine */
+ CLHEP::HepRandomEngine* m_randomEngine;
+ std::string m_randomEngineName; //!< Name of the random number stream
+
+ static Trk::ParticleMasses s_particleMasses; //!< struct of Particle masses
+
+ static double s_main_RutherfordScott; //!< main factor of Rutherford-Scott formula
+ static double s_log_RutherfordScott; //!< log factor of Rutherford-Scott formula
+
+ static double s_main_RossiGreisen; //!< main factor for Rossi-Greisen formula
+ static double s_log_RossiGreisen; //!< main factor for Rossi-Greisen formula
+
+ static double s_projectionFactor; //!< projection factor to scale the projected angle out of the plane
+
+
+ };
+
+
+} // end of namespace
+
+
+#endif // ISF_FATRASTOOLS_MULTIPLESCATTERINGSAMPLERHIGHLAND_H
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/PhotonConversionTool.h b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/PhotonConversionTool.h
index a7807c02065eb4e0001a9684b7d1036db296f21a..88ce25922621200c699147e0003fdf80c5070628 100755
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/PhotonConversionTool.h
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/PhotonConversionTool.h
@@ -17,7 +17,6 @@
// Trk
#include "TrkEventPrimitives/PropDirection.h"
#include "TrkEventPrimitives/ParticleHypothesis.h"
-#include "TrkParameters/TrackParameters.h"
#include "TrkExUtils/MaterialUpdateMode.h"
#include "TrkDetDescrUtils/GeometrySignature.h"
// ISF
@@ -75,12 +74,12 @@ namespace iFatras {
double p) const;
/** interface for processing of the presampled pair production */
- bool doConversion(double time, const Trk::TrackParameters& parm,
+ bool doConversion(double time, const Trk::NeutralParameters& parm,
const Trk::ExtendedMaterialProperties* extMatProp=0) const;
/** interface for processing of the presampled nuclear interactions on layer*/
ISF::ISFParticleVector doConversionOnLayer(const ISF::ISFParticle* parent,
- double time, const Trk::TrackParameters& parm,
+ double time, const Trk::NeutralParameters& parm,
const Trk::ExtendedMaterialProperties *ematprop=0) const;
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/PhysicsValidationTool.h b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/PhysicsValidationTool.h
index b3d1194bb6cae414a79520eaee73a61f119e2af0..032616b056b541c31ef790655555d7db8b66be08 100755
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/PhysicsValidationTool.h
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/PhysicsValidationTool.h
@@ -67,9 +67,18 @@ namespace iFatras
/** AlgTool finalize method */
virtual StatusCode finalize();
-
- /** Fill ntuple */
+
+ /** ISFParticle info: old transport tool */
void saveISFParticleInfo(const ISF::ISFParticle& isp, int endProcess, const Trk::TrackParameters* ePar, double time, double dX0 );
+
+ /** ISFParticle info: new transport tool */
+ void saveISFParticleInfo(const ISF::ISFParticle& isp, const Trk::ExtrapolationCell<Trk::TrackParameters>& ec,
+ Trk::ExtrapolationCode ecode );
+
+ /** ISFParticle info: new transport tool */
+ void saveISFParticleInfo(const ISF::ISFParticle& isp, const Trk::ExtrapolationCell<Trk::NeutralParameters>& ec,
+ Trk::ExtrapolationCode ecode );
+
void saveISFVertexInfo(int process,Amg::Vector3D vertex,const ISF::ISFParticle& isp,Amg::Vector3D primIn,
Amg::Vector3D* primOut, const ISF::ISFParticleVector children);
@@ -85,6 +94,8 @@ namespace iFatras
ATH_MSG_DEBUG("[ fatras setup ] Successfully retrieved " << thandle);
return StatusCode::SUCCESS;
}
+
+ void saveInfo(const ISF::ISFParticle& isp);
/*---------------------------------------------------------------------
* Private members
@@ -105,6 +116,10 @@ namespace iFatras
mutable float m_pph;
mutable float m_p;
mutable float m_eloss;
+ mutable float m_ionloss;
+ mutable float m_radloss;
+ mutable float m_zOaTr;
+ mutable float m_wZ;
mutable float m_thIn;
mutable float m_phIn;
mutable float m_dIn;
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/ProcessSamplingTool.h b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/ProcessSamplingTool.h
index 30c003917dab69d7dea10d6a7a60d160a2ae462e..a0210c6b2343d9870bcc0240f5b3d9357c09a0a7 100755
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/ProcessSamplingTool.h
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/ProcessSamplingTool.h
@@ -23,8 +23,16 @@
#include "TrkParameters/TrackParameters.h"
#include "TrkEventPrimitives/ParticleHypothesis.h"
+namespace ISF {
+ class ITruthSvc;
+}
+
namespace iFatras
{
+ class IHadronicInteractionProcessor;
+ class IPhotonConversionTool;
+ class IParticleDecayHelper;
+ class IPhysicsValidationTool;
/** @class ProcessSamplingTool
@@ -51,8 +59,15 @@ namespace iFatras
virtual StatusCode finalize();
/** Process pre-sampling : to be moved into material updators eventually */
- Trk::PathLimit sampleProcess(double mom, double charge, Trk::ParticleHypothesis pHypothesis);
+ Trk::PathLimit sampleProcess(double mom, double charge, Trk::ParticleHypothesis pHypothesis) const;
+ /** Process simulation */
+ ISF::ISFParticleVector interact(const ISF::ISFParticle* parent,
+ Trk::ExCellCharged& eCell,
+ const Trk::Material* mat) const;
+ ISF::ISFParticleVector interact(const ISF::ISFParticle* parent,
+ Trk::ExCellNeutral& eCell,
+ const Trk::Material* mat) const;
private:
/** templated Tool retrieval - gives unique handling & look and feel */
template <class T> StatusCode retrieveTool(ToolHandle<T>& thandle){
@@ -72,8 +87,25 @@ namespace iFatras
ServiceHandle<IAtRndmGenSvc> m_rndGenSvc;
CLHEP::HepRandomEngine* m_randomEngine;
std::string m_randomEngineName; //!< Name of the random number stream
+
+ /** Truth record */
+ ServiceHandle<ISF::ITruthSvc> m_truthRecordSvc;
+
+ /** hadronic interaction */
+ bool m_hadInt;
+ ToolHandle<IHadronicInteractionProcessor> m_hadIntProcessor;
+
+ /** decay */
+ ToolHandle<IParticleDecayHelper> m_particleDecayer;
+
+ /** IPhotonConversionTool */
+ ToolHandle<iFatras::IPhotonConversionTool> m_conversionTool;
+
+ Trk::ParticleMasses m_particleMasses; //!< Struct of Particle masses
- Trk::ParticleMasses m_particleMasses; //!< Struct of Particle masses
+ /** validation */
+ bool m_validationMode;
+ ToolHandle<IPhysicsValidationTool> m_validationTool;
};
}
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/TransportEngine.h b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/TransportEngine.h
index 583fbc635b1fe3c194d1406025f9d0a2dcdb5b7f..8ae9d784196e70049d821d151e00462653781757 100644
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/TransportEngine.h
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/TransportEngine.h
@@ -16,12 +16,13 @@
#include "AthenaKernel/IAtRndmGenSvc.h"
#include "TrkExInterfaces/IExtrapolationEngine.h"
#include "TrkExUtils/ExtrapolationCell.h"
-// iFatras
-#include "ISF_FatrasInterfaces/ITransportTool.h"
+// ISF
+#include "ISF_Interfaces/IParticleProcessor.h"
// Tracking
#include "TrkEventPrimitives/PdgToParticleHypothesis.h"
#include "TrkEventPrimitives/ParticleHypothesis.h"
+#include "TrkDetDescrUtils/GeometrySignature.h"
class IIncidentSvc;
@@ -39,6 +40,7 @@ namespace iFatras
class ISimHitCreator;
class IParticleDecayHelper;
class IProcessSamplingTool;
+ class IPhysicsValidationTool;
/** @class TransportEngine
@@ -48,7 +50,7 @@ namespace iFatras
@author Noemi Calace -at- cern.ch, Andreas Salzburger -at cern.ch
*/
- class TransportEngine : virtual public ITransportTool,
+ class TransportEngine : virtual public ISF::IParticleProcessor,
public AthAlgTool
{
public:
@@ -74,8 +76,8 @@ namespace iFatras
Trk::ExtrapolationCode eCode,
const Amg::Vector3D& position,
const Amg::Vector3D& momentum,
- double charge,
- double stime);
+ double stime,
+ Trk::GeometrySignature geoID);
/** templated Tool retrieval - gives unique handling & look and feel */
template <class T> StatusCode retrieveTool(ToolHandle<T>& thandle)
@@ -122,6 +124,9 @@ namespace iFatras
Trk::PdgToParticleHypothesis m_pdgToParticleHypothesis;
Trk::ParticleMasses m_particleMasses; //!< Struct of Particle masses
+
+ bool m_validationMode;
+ ToolHandle<IPhysicsValidationTool> m_validationTool;
};
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/TransportTool.h b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/TransportTool.h
index cb5eb58f8158c323cbb4f6d7b55584e78cf320ff..08debd043beefe45892c13edb3d60dd2a9eaf58d 100755
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/TransportTool.h
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/ISF_FatrasTools/TransportTool.h
@@ -16,8 +16,8 @@
#include "AthenaKernel/IAtRndmGenSvc.h"
#include "TrkExInterfaces/ITimedExtrapolator.h"
-// iFatras
-#include "ISF_FatrasInterfaces/ITransportTool.h"
+// IFS
+#include "ISF_Interfaces/IParticleProcessor.h"
// Tracking
#include "TrkEventPrimitives/PdgToParticleHypothesis.h"
@@ -53,7 +53,7 @@ namespace iFatras
@author Andreas.Salzburger -at- cern.ch
*/
- class TransportTool : virtual public ITransportTool,
+ class TransportTool : virtual public ISF::IParticleProcessor,
public AthAlgTool
{
public:
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/cmt/requirements b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/cmt/requirements
index 6454dd51256d41c9da3c5c1ae4cea0299739d0a9..81bf736a0b4899bc45b165820f2b320d7dbe5328 100644
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/cmt/requirements
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/cmt/requirements
@@ -7,46 +7,50 @@ manager Wolfgang Lukas <Wolfgang.Lukas@cern.ch>
public
########## Control #############################################################
use AtlasPolicy AtlasPolicy-*
-use AthenaBaseComps AthenaBaseComps-* Control
-use AthenaKernel AthenaKernel-* Control
-use GaudiInterface GaudiInterface-* External
+use AthenaBaseComps AthenaBaseComps-* Control
+use AthenaKernel AthenaKernel-* Control
+use GaudiInterface GaudiInterface-* External
-########## ISF #################################################################
-use ISF_Event ISF_Event-* Simulation/ISF/ISF_Core
-use ISF_FatrasInterfaces ISF_FatrasInterfaces-* Simulation/ISF/ISF_Fatras
+########## External ############################################################
+use AtlasCLHEP AtlasCLHEP-* External
+########## ISF #################################################################
+use ISF_Event ISF_Event-* Simulation/ISF/ISF_Core
+use ISF_FatrasInterfaces ISF_FatrasInterfaces-* Simulation/ISF/ISF_Fatras
+use ISF_Interfaces ISF_Interfaces-* Simulation/ISF/ISF_Core
+
########## Barcode #############################################################
-use BarcodeInterfaces BarcodeInterfaces-* Simulation/Barcode
+use BarcodeInterfaces BarcodeInterfaces-* Simulation/Barcode
########## Tracking ############################################################
-use TrkDetDescrUtils TrkDetDescrUtils-* Tracking/TrkDetDescr
-use TrkEventPrimitives TrkEventPrimitives-* Tracking/TrkEvent
-use TrkExInterfaces TrkExInterfaces-* Tracking/TrkExtrapolation
-use TrkExUtils TrkExUtils-* Tracking/TrkExtrapolation
-use TrkParameters TrkParameters-* Tracking/TrkEvent
+use TrkDetDescrUtils TrkDetDescrUtils-* Tracking/TrkDetDescr
+use TrkEventPrimitives TrkEventPrimitives-* Tracking/TrkEvent
+use TrkExInterfaces TrkExInterfaces-* Tracking/TrkExtrapolation
+use TrkExUtils TrkExUtils-* Tracking/TrkExtrapolation
+use TrkParameters TrkParameters-* Tracking/TrkEvent
+use EventPrimitives EventPrimitives-* Event
+use GeoPrimitives GeoPrimitives-* DetectorDescription
+use TrkGeometry TrkGeometry-* Tracking/TrkDetDescr
+
private
########## Control #############################################################
use StoreGate StoreGate-* Control
########## External ############################################################
-use AtlasCLHEP AtlasCLHEP-* External
use AtlasROOT AtlasROOT-* External
########## RandomNumbers #######################################################
use AtlasCLHEP_RandomGenerators AtlasCLHEP_RandomGenerators-* Simulation/Tools
########## Tracking ############################################################
-use TrkGeometry TrkGeometry-* Tracking/TrkDetDescr
-use TrkSurfaces TrkSurfaces-* Tracking/TrkDetDescr
-use TrkVolumes TrkVolumes-* Tracking/TrkDetDescr
-use TrkDetDescrInterfaces TrkDetDescrInterfaces-* Tracking/TrkDetDescr
-use TrkTrack TrkTrack-* Tracking/TrkEvent
-use TrkMaterialOnTrack TrkMaterialOnTrack-* Tracking/TrkEvent
-use TrkNeutralParameters TrkNeutralParameters-* Tracking/TrkEvent
-
-########## ISF #################################################################
-use ISF_Interfaces ISF_Interfaces-* Simulation/ISF/ISF_Core
+use TrkVolumes TrkVolumes-* Tracking/TrkDetDescr
+use TrkDetDescrInterfaces TrkDetDescrInterfaces-* Tracking/TrkDetDescr
+use TrkTrack TrkTrack-* Tracking/TrkEvent
+use TrkMaterialOnTrack TrkMaterialOnTrack-* Tracking/TrkEvent
+use TrkNeutralParameters TrkNeutralParameters-* Tracking/TrkEvent
+use TrkSurfaces TrkSurfaces-* Tracking/TrkDetDescr
+use TrkExEngine TrkExEngine-* Tracking/TrkExtrapolation
########## ISF #################################################################
@@ -67,6 +71,6 @@ library ISF_FatrasTools *.cxx components/*.cxx
apply_pattern component_library
# use the following to compile with debug information
-#private
+private
#macro cppdebugflags '$(cppdebugflags_s)'
#macro_remove componentshr_linkopts "-Wl,-s"
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/HadIntProcessorParametric.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/HadIntProcessorParametric.cxx
index efae1ce5241196fa1515e5148b9c9348674d11a0..7bcdb02f57b89f85d87dc5606f92d3d3a1cb0277 100644
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/HadIntProcessorParametric.cxx
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/HadIntProcessorParametric.cxx
@@ -199,9 +199,10 @@ StatusCode iFatras::HadIntProcessorParametric::finalize()
}
-bool iFatras::HadIntProcessorParametric::hadronicInteraction(const Trk::TrackParameters& parm, double p, double /*E*/,
- const Trk::MaterialProperties& mprop, double pathCorrection,
- Trk::ParticleHypothesis particle) const
+bool iFatras::HadIntProcessorParametric::hadronicInteraction(const Amg::Vector3D& position, const Amg::Vector3D& momentum,
+ double p, double /*E*/, double charge,
+ const Trk::MaterialProperties& mprop, double pathCorrection,
+ Trk::ParticleHypothesis particle) const
{
const Trk::MaterialProperties* extMprop = dynamic_cast<const Trk::MaterialProperties*>(&mprop);
double prob = 0.;
@@ -209,7 +210,7 @@ bool iFatras::HadIntProcessorParametric::hadronicInteraction(const Trk::TrackPar
// m_hadIntProbScale is used later, not here
if (extMprop && !m_hadIntFromX0) {
- double al = absorptionLength(extMprop, p, parm.charge(), particle); // in mm
+ double al = absorptionLength(extMprop, p, charge, particle); // in mm
if (al>0.) prob = exp( (-1.)*pathCorrection*extMprop->thickness()/al );
else prob = exp( (-1.)*pathCorrection*extMprop->thicknessInL0());
@@ -228,8 +229,8 @@ bool iFatras::HadIntProcessorParametric::hadronicInteraction(const Trk::TrackPar
// (1. - prob) is generally O(0.01), so this is the right way to scale it
// TODO fix time info (if needed)
if (CLHEP::RandFlat::shoot(m_randomEngine) < (1. - prob) * m_hadIntProbScale ) return recordHadState(0.,p,
- parm.position(),
- parm.momentum().unit(),
+ position,
+ momentum.unit(),
particle);
// no hadronic interactions were computed
@@ -239,7 +240,7 @@ bool iFatras::HadIntProcessorParametric::hadronicInteraction(const Trk::TrackPar
/** absorption length */
double iFatras::HadIntProcessorParametric::absorptionLength(const Trk::MaterialProperties* mat,
- double p, double, Trk::ParticleHypothesis particle) const {
+ double p, double, Trk::ParticleHypothesis particle) const {
double al = mat->l0();
/* // these parametrization comes from comparison with G4 sampler, but give too many interactions
@@ -676,7 +677,7 @@ ISF::ISFParticleVector iFatras::HadIntProcessorParametric::getHadState(const ISF
}
-bool iFatras::HadIntProcessorParametric::doHadronicInteraction(double time,const Trk::TrackParameters& parm,
+bool iFatras::HadIntProcessorParametric::doHadronicInteraction(double time, const Amg::Vector3D& position, const Amg::Vector3D& momentum,
const Trk::Material* /*ematprop*/,
Trk::ParticleHypothesis particle, bool processSecondaries) const {
// get parent particle
@@ -684,7 +685,7 @@ bool iFatras::HadIntProcessorParametric::doHadronicInteraction(double time,const
// something is seriously wrong if there is no parent particle
assert(parent);
- ISF::ISFParticleVector ispVec=getHadState(parent,time, parm.momentum().mag(), parm.position(), parm.momentum().unit(),particle);
+ ISF::ISFParticleVector ispVec=getHadState(parent, time, momentum.mag(), position, momentum.unit(), particle);
// having no secondaries does not necessarily mean the interaction did not take place : TODO : add flag into ::getHadState
// if (!ispVec.size()) return false;
@@ -699,12 +700,13 @@ bool iFatras::HadIntProcessorParametric::doHadronicInteraction(double time,const
}
-ISF::ISFParticleVector iFatras::HadIntProcessorParametric::doHadIntOnLayer(const ISF::ISFParticle* parent, double time,const Trk::TrackParameters& parm,
- const Trk::MaterialProperties* /*ematprop*/,
- Trk::ParticleHypothesis particle) const {
-
- return getHadState(parent, time, parm.momentum().mag(), parm.position(), parm.momentum().unit(),particle);
+ISF::ISFParticleVector iFatras::HadIntProcessorParametric::doHadIntOnLayer(const ISF::ISFParticle* parent, double time,
+ const Amg::Vector3D& position, const Amg::Vector3D& momentum,
+ const Trk::Material* /*emat*/,
+ Trk::ParticleHypothesis particle) const {
+ return getHadState(parent, time, momentum.mag(), position, momentum.unit(), particle);
+
}
/** interface for processing of the presampled nuclear interaction */
@@ -716,17 +718,17 @@ bool iFatras::HadIntProcessorParametric::recordHadState(double time, double p,
const ISF::ISFParticle *parent = ISF::ParticleClipboard::getInstance().getParticle();
// something is seriously wrong if there is no parent particle
assert(parent);
-
- ISF::ISFParticleVector ispVec=getHadState(parent,time, p, vertex, particleDir, particle);
-
+
+ ISF::ISFParticleVector ispVec=getHadState(parent, time, p, vertex, particleDir, particle);
+
// having no secondaries does not necessarily mean the interaction did not take place : TODO : add flag into ::getHadState
// if (!ispVec.size()) return false;
-
+
// push onto ParticleStack
if (ispVec.size() ) {
for (unsigned int ic=0; ic<ispVec.size(); ic++) m_particleBroker->push(ispVec[ic], parent);
}
-
+
return true;
}
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McBetheHeitlerEnergyLossUpdator.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McBetheHeitlerEnergyLossUpdator.cxx
deleted file mode 100755
index 8eac108424c4aaa0e2b57bd1eb755e581d8498d0..0000000000000000000000000000000000000000
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McBetheHeitlerEnergyLossUpdator.cxx
+++ /dev/null
@@ -1,138 +0,0 @@
-/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
-
-/***************************************************************************************
- McBetheHeitlerEnergyLossUpdator.cxx - description
- ---------------------------------------------------
-begin : Wednesday 28th June 2006
-author : atkinson
-email : Tom.Atkinson@cern.ch
-decription : Implementation code for McBetheHeitlerEnergyLossUpdator class
-***************************************************************************************/
-
-// class header include
-#include "ISF_FatrasTools/McBetheHeitlerEnergyLossUpdator.h"
-
-#include "TrkGeometry/MaterialProperties.h"
-#include "TrkMaterialOnTrack/EnergyLoss.h"
-#include "CLHEP/Random/RandGamma.h"
-#include <cmath>
-
-
-// static partilce masses
-Trk::ParticleMasses iFatras::McBetheHeitlerEnergyLossUpdator::s_particleMasses;
-
-iFatras::McBetheHeitlerEnergyLossUpdator::McBetheHeitlerEnergyLossUpdator( const std::string& type, const std::string& name, const IInterface* parent )
- :
- AthAlgTool( type, name, parent ),
- m_rndmGenSvc("AtDSFMTGenSvc", name),
- m_randomEngine(0),
- m_randomEngineName("FatrasRnd"),
- m_cashedPathLength(0.),
- m_numberPointsInPDF(100000),
- m_scaleFactor(1.0)
-{
-
- declareInterface<IEnergyLossUpdator>(this);
- // get the property from outside
- declareProperty("ScaleFactor", m_scaleFactor);
- declareProperty("RandomNumberService" , m_rndmGenSvc , "Random number generator");
- declareProperty("RandomStreamName" , m_randomEngineName , "Name of the random number stream");
-
-}
-
-iFatras::McBetheHeitlerEnergyLossUpdator::~McBetheHeitlerEnergyLossUpdator(){}
-
-StatusCode iFatras::McBetheHeitlerEnergyLossUpdator::initialize()
-{
-
- ATH_MSG_INFO( "initialize()" );
-
- if ( m_rndmGenSvc.retrieve().isFailure() ){
- ATH_MSG_FATAL( "Could not retrieve " << m_rndmGenSvc );
- return StatusCode::FAILURE;
- }
- //Get own engine with own seeds:
- m_randomEngine = m_rndmGenSvc->GetEngine(m_randomEngineName);
- if (!m_randomEngine) {
- ATH_MSG_FATAL( "Could not get random engine '" << m_randomEngineName << "'" );
- return StatusCode::FAILURE;
- }
- return StatusCode::SUCCESS;
-}
-
-StatusCode iFatras::McBetheHeitlerEnergyLossUpdator::finalize()
-{
-
- ATH_MSG_INFO( "finalize() successful" );
- return StatusCode::SUCCESS;
-}
-
-Trk::EnergyLoss* iFatras::McBetheHeitlerEnergyLossUpdator::energyLoss(
- const Trk::MaterialProperties& materialProperties,
- double pInitial,
- double pathCorrection,
- Trk::PropDirection direction,
- Trk::ParticleHypothesis,
- bool) const
-{
-
- double pathLength = pathCorrection * materialProperties.thickness() / materialProperties.x0();
-
- double u = CLHEP::RandGamma::shoot(m_randomEngine, pathLength / log(2.), 1.);
-
- double z = exp( -1. * u );
-
- // initial Energy
- double E = sqrt(pInitial*pInitial+s_particleMasses.mass[Trk::electron]*s_particleMasses.mass[Trk::electron]);
-
-
- /*
- if ( fabs( pathLength - m_cashedPathLength ) > 1.e-6 || !m_pdf.size() ){
-
- m_pdf = betheHeitlerPDF( pathLength );
-
- m_cashedPathLength = pathLength;
-
- }
-
- Rndm::Numbers betheHeitler( m_rndmGenSvc, Rndm::DefinedPdf( m_pdf, 0 ) );
-
- double z = betheHeitler.shoot();
-
- */
-
- double deltaE(0.);
-
- if ( direction == Trk::alongMomentum )
- deltaE = m_scaleFactor*E * ( z - 1. );
-
- else
- deltaE = m_scaleFactor*E * ( 1. / z - 1. );
-
- return new Trk::EnergyLoss(deltaE, z);
-}
-
-std::vector<double> iFatras::McBetheHeitlerEnergyLossUpdator::betheHeitlerPDF( double pathLength ) const
-{
-
- std::vector<double> pdf;
-
- double c = pathLength / log(2);
-
- double sizeOfIncrement = 1. / (double) m_numberPointsInPDF;
-
- double currentZ(0.);
-
- int index(0);
-
- for( ; index < m_numberPointsInPDF ; ++index ){
- double pdfEntry = pow( -1. * log( currentZ ), c - 1. ) / CLHEP::RandGamma::shoot(m_randomEngine, c, 1.);
- pdf.push_back( pdfEntry );
- currentZ += sizeOfIncrement;
- }
-
- return pdf;
-
-}
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McEnergyLossUpdator.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McEnergyLossUpdator.cxx
index 9c85c5cde374c37049b13c2ec0b5cb3f96069029..b3313722438dda4d40601076226d8535ffdcecd5 100755
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McEnergyLossUpdator.cxx
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McEnergyLossUpdator.cxx
@@ -24,8 +24,6 @@ iFatras::McEnergyLossUpdator::McEnergyLossUpdator( const std::string& type, cons
AthAlgTool( type, name, parent ),
m_energyLossUpdator("Trk::EnergyLossUpdator/AtlasEnergyLossUpdator" ),
m_energyLossDistribution(2),
- m_dedicatedElectronUpdator(true),
- m_elEnergyLossUpdator("iFatras::McBetheHeitlerEnergyLossUpdator/FatrasBetheHeitlerEnergyLossUpdator" ),
m_rndGenSvc("AtDSFMTGenSvc", name),
m_randomEngine(0),
m_randomEngineName("FatrasRnd")
@@ -33,9 +31,6 @@ iFatras::McEnergyLossUpdator::McEnergyLossUpdator( const std::string& type, cons
declareInterface<IEnergyLossUpdator>(this);
// heavy particles energy loss
declareProperty( "EnergyLossUpdator", m_energyLossUpdator);
- // dedicated electron updates
- declareProperty( "UseElectronUpdator", m_dedicatedElectronUpdator );
- declareProperty( "ElectronEnergyLossUpdator", m_elEnergyLossUpdator);
// random number service
declareProperty( "RandomNumberService", m_rndGenSvc);
declareProperty( "RandomStreamName", m_randomEngineName, "Name of the random number stream");
@@ -55,16 +50,6 @@ StatusCode iFatras::McEnergyLossUpdator::initialize()
return StatusCode::FAILURE;
}
-
- // Retrieve the dedicated electron energz loss updator
- if (m_dedicatedElectronUpdator){
- // Retrieve the energy loss updator tool
- if ( m_elEnergyLossUpdator.retrieve().isFailure() ){
- ATH_MSG_FATAL( "Could not retrieve " << m_elEnergyLossUpdator );
- return StatusCode::FAILURE;
- }
- }
-
// get the random generator serice
if (m_rndGenSvc.retrieve().isFailure()){
ATH_MSG_FATAL( "Could not retrieve " << m_rndGenSvc );
@@ -108,7 +93,7 @@ Trk::EnergyLoss* iFatras::McEnergyLossUpdator::energyLoss(
bool) const
{
- double simulatedDeltaE = 0.;
+ bool mpvSwitch = m_energyLossDistribution<2 ? false : true;
// get the number of the material effects distribution
Trk::EnergyLoss* sampledEloss = m_energyLossUpdator->energyLoss(
@@ -117,64 +102,37 @@ Trk::EnergyLoss* iFatras::McEnergyLossUpdator::energyLoss(
pathCorrection,
direction,
particleHypothesis,
- true);
+ mpvSwitch);
-
- double energyLoss = 0.;
- double energyLossSigma = 0.;
- double m = s_particleMasses.mass[particleHypothesis];
- double E = sqrt(momentum*momentum+m*m);
-
- // get the sampled energy loss
- if(sampledEloss){
-
- // get the energy loss and its sigma
- // for gauss: mean, sigma
- // for landau: mpv, sigma
- energyLoss = fabs(sampledEloss->deltaE());
- energyLossSigma = sampledEloss->sigmaDeltaE();
- // memory cleanup
- delete sampledEloss;
- }
+ if (!sampledEloss) return 0;
+ // smear according to preselected distribution
switch (m_energyLossDistribution) {
// no straggling
- case 0 : { simulatedDeltaE += energyLoss; } break;
+ case 0 : { } break;
// gaussian smearing
- case 1 : { simulatedDeltaE += energyLoss + energyLossSigma * CLHEP::RandGaussZiggurat::shoot(m_randomEngine); } break;
+ case 1 : {
+ float deIoni = sampledEloss->sigmaIoni() * CLHEP::RandGaussZiggurat::shoot(m_randomEngine);
+ float deRad = sampledEloss->sigmaRad() * CLHEP::RandGaussZiggurat::shoot(m_randomEngine);
+ sampledEloss->update(deIoni,0.,deRad,0.,false);
+ } break;
// landau smearing
default : {
- simulatedDeltaE += (energyLoss+energyLossSigma*CLHEP::RandLandau::shoot(m_randomEngine));
-
+ float deIoni = -sampledEloss->sigmaIoni() * CLHEP::RandLandau::shoot(m_randomEngine); // TODO :check sign
+ float deRad = -sampledEloss->sigmaRad() * CLHEP::RandLandau::shoot(m_randomEngine); // TODO :check sign
+ sampledEloss->update(deIoni,0.,deRad,0.,false);
} break;
}
- // transform into negative nergyloss
- simulatedDeltaE *= direction;
-
- /* ST radiative corrections treated separately
- // it is an electron and has radiative corrections
- if (m_dedicatedElectronUpdator && particleHypothesis==Trk::electron) {
- // !< @ todo introduce cut for the calculation of bremsstrahlung
- Trk::EnergyLoss* radiationLoss = m_elEnergyLossUpdator->energyLoss(
- materialProperties,
- momentum,
- pathCorrection,
- direction,
- particleHypothesis);
- // add the simulated bethe heitler energy loss
- double radiationDeltaE = radiationLoss ? radiationLoss->deltaE() : 0.;
-
- simulatedDeltaE += (-1*radiationDeltaE);
- delete radiationLoss;
- }
- */
+ // protection due to straggling - maximum energy loss is E-m
+ double m = s_particleMasses.mass[particleHypothesis];
+ double E = sqrt(momentum*momentum+m*m);
- // protection due to straggling - maximum energy loss is E
- simulatedDeltaE = simulatedDeltaE*simulatedDeltaE > E*E ? E : simulatedDeltaE;
-
- ATH_MSG_VERBOSE( "[+] FatrasEnergyUpdator created random deltaP as : " << -1.*simulatedDeltaE );
+ if (sampledEloss->deltaE()+E<m ) { // particle stopping - rest energy
+ float dRad_rest = m-E-sampledEloss->deltaE();
+ sampledEloss->update(0.,0.,dRad_rest,0.,false);
+ }
- return new Trk::EnergyLoss( -1.*simulatedDeltaE, 0. );
+ return sampledEloss;
}
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsEngine.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsEngine.cxx
index 83b5e76867635ed2bbe7dbd4848c338093f2de7d..3f9888685cad482a6d22fa77802f7a897e551264 100644
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsEngine.cxx
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsEngine.cxx
@@ -10,64 +10,62 @@
#include "ISF_FatrasTools/McMaterialEffectsEngine.h"
// Gaudi Kernel
-#include "GaudiKernel/IRndmGenSvc.h"
-#include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/DataSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
+
// ISF includes
#include "ISF_Interfaces/IParticleBroker.h"
#include "ISF_Interfaces/ITruthSvc.h"
#include "ISF_Event/ISFParticle.h"
#include "ISF_Event/ISFTruthIncident.h"
-#include "ISF_Event/ParticleClipboard.h"
#include "ISF_Event/ParticleUserInformation.h"
// iFatras
-// #include "ISF_FatrasInterfaces/IHadronicInteractionProcessor.h"
-// #include "ISF_FatrasInterfaces/IProcessSamplingTool.h"
-// #include "ISF_FatrasInterfaces/IPhysicsValidationTool.h"
-// #include "ISF_FatrasInterfaces/IPhotonConversionTool.h"
-// #include "ISF_FatrasInterfaces/IParticleDecayer.h"
-// FastSimulation
-// #include "FastSimulationEvent/GenParticleEnergyDepositMap.h"
+#include "ISF_FatrasInterfaces/IHadronicInteractionProcessor.h"
+#include "ISF_FatrasInterfaces/IProcessSamplingTool.h"
+#include "ISF_FatrasInterfaces/IPhysicsValidationTool.h"
+#include "ISF_FatrasInterfaces/IParticleDecayHelper.h"
+#include "ISF_FatrasInterfaces/IEnergyLossSampler.h"
+#include "ISF_FatrasInterfaces/IMultipleScatteringSampler.h"
// Trk inlcude
-#include "TrkExInterfaces/IEnergyLossUpdator.h"
-#include "TrkExInterfaces/IMultipleScatteringUpdator.h"
#include "TrkParameters/TrackParameters.h"
#include "TrkEventPrimitives/ParamDefs.h"
#include "TrkSurfaces/Surface.h"
#include "TrkDetDescrInterfaces/ITrackingGeometrySvc.h"
-#include "TrkGeometry/Layer.h"
#include "TrkGeometry/TrackingGeometry.h"
#include "TrkGeometry/TrackingVolume.h"
#include "TrkGeometry/MaterialProperties.h"
#include "TrkVolumes/CylinderVolumeBounds.h"
#include "TrkMaterialOnTrack/EnergyLoss.h"
+#include "TrkExEngine/ExtrapolationMacros.h"
+#include "TrkExInterfaces/ITimedExtrapolator.h"
+
// CLHEP
#include "CLHEP/Units/SystemOfUnits.h"
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Random/RandFlat.h"
-#include "CLHEP/Random/RandPoisson.h"
#include "CLHEP/Random/RandLandau.h"
-#include "AtlasCLHEP_RandomGenerators/RandGaussZiggurat.h"
+
// Validation mode - TTree includes
#include "TTree.h"
#include "GaudiKernel/ITHistSvc.h"
// STD
#include <math.h>
-// temporary
-#include "TrkGeometry/TrackingVolume.h"
-
// constructor
iFatras::McMaterialEffectsEngine::McMaterialEffectsEngine(const std::string& t, const std::string& n, const IInterface* p) :
AthAlgTool(t,n,p),
+ m_samplingTool(""),
m_eLoss(true),
- m_eLossUpdator("iFatras::McEnergyLossUpdator/FatrasEnergyLossUpdator"),
+ m_eLossSampler(""),
m_ms(true),
- m_msUpdator("Trk::MultipleScatteringUpdator/AtlasMultipleScatteringUpdator"),
+ m_msSampler(""),
+ m_particleDecayHelper(""),
+ m_parametricScattering(true),
+ m_uniformHertzDipoleAngle(false),
+ m_processCode(3),
m_minimumMomentum(50.*CLHEP::MeV),
- m_use_msUpdator(false),
- m_referenceMaterial(true),
+ m_minimumBremPhotonMomentum(50.*CLHEP::MeV),
+ m_createBremPhoton(true),
m_bendingCorrection(false),
m_rndGenSvc("AtDSFMTGenSvc", n),
m_randomEngine(0),
@@ -75,24 +73,50 @@ iFatras::McMaterialEffectsEngine::McMaterialEffectsEngine(const std::string& t,
m_trackingGeometry(0),
m_trackingGeometrySvc("AtlasTrackingGeometrySvc", n),
m_trackingGeometryName("AtlasTrackingGeometry"),
+ m_validationMode(false),
+ m_validationTool(""),
m_particleBroker("ISF_ParticleParticleBroker", n),
m_truthRecordSvc("ISF_TruthRecordSvc", n),
- m_layer(0)
+ m_oneOverThree(1./3.),
+ m_projectionFactor(sqrt(2.)/2.)
{
declareInterface<Trk::IMaterialEffectsEngine>(this);
- // the tool parameters -----------------------------------------------------
+
+ // steering of the screen outoput (SOP)
+ declareProperty("OutputPrefix" , m_sopPrefix);
+ declareProperty("OutputPostfix" , m_sopPostfix);
+
+ // Tool Parameters
+ // Energy Loss Samplers
declareProperty("EnergyLoss" , m_eLoss);
- declareProperty("EnergyLossUpdator" , m_eLossUpdator);
+ declareProperty("EnergyLossSampler" , m_eLossSampler);
+ // Multiple Scattering Sampler
declareProperty("MultipleScattering" , m_ms);
- declareProperty("MultipleScatteringUpdator" , m_msUpdator);
- // the steering --------------------------------------------------------------
- declareProperty("UseMultipleScatteringUpdator" , m_use_msUpdator);
+ declareProperty("MultipleScatteringSampler" , m_msSampler);
+ // tool handle for the particle decayer
+ declareProperty( "ParticleDecayHelper", m_particleDecayHelper );
+ // Process sampling
+ declareProperty( "ProcessSamplingTool", m_samplingTool );
+
+ // MC Truth Properties
+ declareProperty("BremProcessCode" , m_processCode, "MCTruth Physics Process Code");
+
+ // The Steering
+ declareProperty("ParametericScattering" , m_parametricScattering);
declareProperty("MomentumCut" , m_minimumMomentum);
+ declareProperty("MinimumBremPhotonMomentum" , m_minimumBremPhotonMomentum);
+ declareProperty("CreateBremPhotons" , m_createBremPhoton);
+ declareProperty("BendingCorrection" , m_bendingCorrection);
+
+ // Validation mode
+ declareProperty("ValidationMode" , m_validationMode);
+ declareProperty("PhysicsValidationTool" , m_validationTool);
+
// TrackingGeometry Service
- declareProperty("TrackingGeometrySvc", m_trackingGeometrySvc);
+ declareProperty("TrackingGeometrySvc" , m_trackingGeometrySvc);
declareProperty("RandomNumberService" , m_rndGenSvc , "Random number generator");
declareProperty("RandomStreamName" , m_randomEngineName , "Name of the random number stream");
-
+
// ISF Services and Tools
declareProperty("ParticleBroker" , m_particleBroker , "ISF Particle Broker Svc");
declareProperty("TruthRecordSvc" , m_truthRecordSvc , "ISF Particle Truth Svc");
@@ -107,414 +131,811 @@ iFatras::McMaterialEffectsEngine::~McMaterialEffectsEngine()
StatusCode iFatras::McMaterialEffectsEngine::initialize()
{
- ATH_MSG_INFO( "initialize()" );
+ ATH_MSG_INFO( ""<<"init"<<""<<"starting initialize()" );
- // retrieve the energy loss updator
+ // retrieve the energy loss sampler
if (m_eLoss){
- if (m_eLossUpdator.retrieve().isFailure()){
- ATH_MSG_FATAL( "Could not retrieve " << m_eLossUpdator );
+ ATH_MSG_VERBOSE(""<<"init"<<""<<"Running WITH EnergyLossSampler ");
+ if (m_eLossSampler.retrieve().isFailure()){
+ ATH_MSG_FATAL("Could not retrieve " << m_eLossSampler );
return StatusCode::FAILURE;
- } else
- ATH_MSG_VERBOSE( "Successfully retrieved " << m_eLossUpdator );
- }
-
- // retrieve the multiple scattering updator tool
- if (m_ms && m_use_msUpdator){
- if (m_msUpdator.retrieve().isFailure()){
- ATH_MSG_FATAL( "Could not retrieve " << m_msUpdator );
+ } else ATH_MSG_VERBOSE(""<<"init"<<""<<"Successfully retrieved " << m_eLossSampler );
+ } else ATH_MSG_VERBOSE(""<<"init"<<""<<"Running WITHOUT EnergyLossSampler ");
+
+ // retrieve the multiple scattering sampler tool
+ if (m_ms){
+ ATH_MSG_VERBOSE(""<<"init"<<""<<"Running WITH MultipleScatteringSamper ");
+ if (m_msSampler.retrieve().isFailure()){
+ ATH_MSG_FATAL( "Could not retrieve " << m_msSampler );
return StatusCode::FAILURE;
- } else
- ATH_MSG_VERBOSE( "Successfully retrieved " << m_msUpdator );
- }
-
+ } else ATH_MSG_VERBOSE(""<<"init"<<""<<"Successfully retrieved " << m_msSampler );
+ } else ATH_MSG_VERBOSE(""<<"init"<<""<<"Running WITHOUT MultipleScatteringSamper ");
+
+ // retrieve the process sampling tool
+ if (m_samplingTool.retrieve().isFailure()){
+ ATH_MSG_FATAL( "Could not retrieve " << m_samplingTool );
+ return StatusCode::FAILURE;
+ } else ATH_MSG_VERBOSE(""<<"init"<<""<<"Successfully retrieved " << m_samplingTool );
+
+ // retrieve the decayHelper
+ if (m_particleDecayHelper.retrieve().isFailure()){
+ ATH_MSG_FATAL( "Could not retrieve " << m_particleDecayHelper );
+ return StatusCode::FAILURE;
+ } else ATH_MSG_VERBOSE(""<<"init"<<""<<"Successfully retrieved " << m_particleDecayHelper );
+
// get the random generator service
if (m_rndGenSvc.retrieve().isFailure()){
ATH_MSG_FATAL( "Could not retrieve " << m_rndGenSvc );
return StatusCode::FAILURE;
} else
- ATH_MSG_VERBOSE( "Successfully retrieved " << m_rndGenSvc );
-
- //Get own engine with own seeds:
+ ATH_MSG_VERBOSE(""<<"init"<<""<<"Successfully retrieved " << m_rndGenSvc );
+
+ // 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;
- }
-
+ } else
+ ATH_MSG_VERBOSE(""<<"init"<<""<<"Successfully got random engine '" << m_randomEngineName << "'" );
+
// get the tracking geometry for layer lookup
// get the TrackingGeometrySvc
if (m_trackingGeometrySvc.retrieve().isSuccess()){
- ATH_MSG_INFO( "Successfully retrieved " << m_trackingGeometrySvc );
+ ATH_MSG_VERBOSE(""<<"init"<<""<<"Successfully retrieved " << m_trackingGeometrySvc );
m_trackingGeometryName = m_trackingGeometrySvc->trackingGeometryName();
} else {
ATH_MSG_WARNING( "Couldn't retrieve " << m_trackingGeometrySvc << ". " );
ATH_MSG_WARNING( " -> Trying to retrieve default '" << m_trackingGeometryName << "' from DetectorStore." );
- }
-
+ }
+
// ISF Services
if (m_particleBroker.retrieve().isFailure()){
ATH_MSG_FATAL( "Could not retrieve " << m_particleBroker );
return StatusCode::FAILURE;
- }
+ } else ATH_MSG_VERBOSE(""<<"init"<<""<<"Successfully retrieved " << m_particleBroker );
+
if (m_truthRecordSvc.retrieve().isFailure()){
ATH_MSG_FATAL( "Could not retrieve " << m_truthRecordSvc );
return StatusCode::FAILURE;
+ } else ATH_MSG_VERBOSE(""<<"init"<<""<<"Successfully retrieved " << m_truthRecordSvc );
+
+ if (m_validationMode){
+ // retrieve the physics validation tool
+ if (m_validationTool.retrieve().isFailure()){
+ ATH_MSG_FATAL( "Could not retrieve " << m_validationTool );
+ return StatusCode::FAILURE;
+ } else
+ ATH_MSG_VERBOSE( "Successfully retrieved " << m_validationTool );
}
- return StatusCode::SUCCESS;
+ return StatusCode::SUCCESS;
}
// finalize
StatusCode iFatras::McMaterialEffectsEngine::finalize()
{
- ATH_MSG_INFO( "finalize() successful" );
+ ATH_MSG_INFO(""<<"fini "<<""<<"finalize() successful" );
return StatusCode::SUCCESS;
}
-Trk::ExtrapolationCode iFatras::McMaterialEffectsEngine::handleMaterial(Trk::ExCellCharged& eCell,
+Trk::ExtrapolationCode iFatras::McMaterialEffectsEngine::handleMaterial(Trk::ExCellCharged& ecCharged,
Trk::PropDirection dir,
Trk::MaterialUpdateStage matupstage) const
-{
- m_layer = eCell.leadLayer;
-
- // get the material properties
- m_matProp = 0;
+{
+ // ATH_MSG_DEBUG(++ecCharged.navigationStep, "handleMaterial"<<"char"<<"handleMaterial for charge particle called.");
+ return handleMaterialT<Trk::TrackParameters> (ecCharged, dir, matupstage);
- if (m_referenceMaterial){
- if (m_layer && m_layer->layerMaterialProperties()) {
- // get the LayerMaterialProperties
- const Trk::LayerMaterialProperties* layerMaterial = m_layer->layerMaterialProperties();
- m_matProp = layerMaterial ? layerMaterial->fullMaterial( eCell.leadParameters->position() ) : 0;
- }
- }
+}
- // the Extrapolator made sure that the layer is the lead layer && the parameters are the lead parameters
- if (eCell.leadLayer && eCell.leadLayer->layerMaterialProperties()){
- ATH_MSG_DEBUG( "handleMaterial for TrackParameters.");
- // update the track parameters
-
- eCell.leadParameters = updateTrackParameters(*eCell.leadParameters,eCell,dir,matupstage);
- }
+Trk::ExtrapolationCode iFatras::McMaterialEffectsEngine::handleMaterial(Trk::ExCellNeutral& ecNeutral,
+ Trk::PropDirection dir,
+ Trk::MaterialUpdateStage matupstage) const
+{
+ //ATH_MSG_DEBUG(++ecNeutral.navigationStep, "handleMaterial"<<"neut"<<"handleMaterial for neutral particle called.");
+ return handleMaterialT<Trk::NeutralParameters> (ecNeutral, dir, matupstage);
+}
+
+void iFatras::McMaterialEffectsEngine::multipleScatteringUpdate(const Trk::TrackParameters& pars,
+ AmgVector(5)& parameters,
+ double simTheta, double num_deltaPhi) const
+{
+ // parametric scattering - independent in x/y
+ if (m_parametricScattering){
+ ATH_MSG_VERBOSE("[msupdate]"<<"MultipleScatteringUpdate"<<""<<"Using parametric scattering." );
+ // the initial values
+ double theta = parameters[Trk::theta];
+ double phi = parameters[Trk::phi];
+ double sinTheta = (theta*theta > 10e-10) ? sin(theta) : 1.;
+
+ // sample them in an independent way
+ double deltaTheta = m_projectionFactor*simTheta;
+ double deltaPhi = m_projectionFactor*num_deltaPhi/sinTheta;
+
+ phi += deltaPhi;
+ if (phi >= M_PI) phi -= M_PI;
+ else if (phi < -M_PI) phi += M_PI;
+ if (theta > M_PI) theta -= M_PI;
+
+ ATH_MSG_VERBOSE("[msupdate]"<<"MultipleScatteringUpdate"<<""<<"deltaPhi / deltaTheta = " << deltaPhi << " / " << deltaTheta );
+
+ // assign the new values
+ parameters[Trk::phi] = phi;
+ parameters[Trk::theta] = fabs(theta + deltaTheta);
- return Trk::ExtrapolationCode::InProgress;
-}
+ } else {
+ double thetaMs = simTheta;
+ double psi = 2.*M_PI*CLHEP::RandFlat::shoot(m_randomEngine);
+ // more complex but "more true" -
+ CLHEP::Hep3Vector parsMomHep( pars.momentum().x(), pars.momentum().y(), pars.momentum().z() );
+ CLHEP::Hep3Vector newDirectionHep( parsMomHep.unit().x(), parsMomHep.unit().y(), parsMomHep.unit().z());
+ double x = -newDirectionHep.y();
+ double y = newDirectionHep.x();
+ double z = 0.;
+ // if it runs along the z axis - no good ==> take the x axis
+ if (newDirectionHep.z()*newDirectionHep.z() > 0.999999)
+ x = 1.; y=0.;
+ // deflector direction
+ CLHEP::Hep3Vector deflector(x,y,z);
+ // rotate the new direction for scattering
+ newDirectionHep.rotate(thetaMs, deflector);
+ // and arbitrarily in psi
+ newDirectionHep.rotate(psi, parsMomHep);
+
+ ATH_MSG_VERBOSE("[msupdate]"<<"MultipleScatteringUpdate"<<""<<"deltaPsi / deltaTheta = " << psi << " / " << thetaMs );
+
+ // assign the new values
+ parameters[Trk::phi] = newDirectionHep.phi();
+ parameters[Trk::theta] = newDirectionHep.theta();
+ }
+}
-const Trk::TrackParameters* iFatras::McMaterialEffectsEngine::updateTrackParameters( const Trk::TrackParameters& parm,
- Trk::ExCellCharged& eCell,
- // const Trk::Layer& lay,
- // Trk::TimeLimit& timeLim, Trk::PathLimit& pathLim,
- // Trk::GeometrySignature /*geoID*/,
- Trk::PropDirection dir,
- Trk::MaterialUpdateStage matupstage) const
- // Trk::ParticleHypothesis particle) const
-{
+void iFatras::McMaterialEffectsEngine::recordBremPhoton(const ISF::ISFParticle* parent,
+ double time,
+ double pElectron,
+ double gammaE,
+ const Amg::Vector3D& vertex,
+ Amg::Vector3D& particleDir,
+ Trk::ParticleHypothesis particle) const
+{
+ ISF::ISFParticle *bremPhot = bremPhoton(parent,time,pElectron,gammaE,vertex,particleDir,particle);
+ m_particleBroker->push(bremPhot, parent);
- double matFraction = 0.;
+}
- const Trk::TrackParameters* tParameters = &(parm);
+void iFatras::McMaterialEffectsEngine::recordBremPhotonLay(const ISF::ISFParticle* parent,
+ double time,
+ double pElectron,
+ double gammaE,
+ const Amg::Vector3D& vertex,
+ Amg::Vector3D& particleDir,
+ Trk::ParticleHypothesis particle,
+ Trk::PropDirection dir,
+ float matFraction ) const {
+
+ ISF::ISFParticle *bremPhot = bremPhoton(parent,time,pElectron,gammaE,vertex,particleDir,particle);
+
+ // process the photon on layer
+ if (gammaE<m_minimumMomentum) return;
+ Trk::PathLimit pLim = m_samplingTool->sampleProcess(gammaE,0,Trk::photon);
+
+ // save presampled information
+ ISF::ParticleUserInformation* childInfo = bremPhot->getUserInformation();
+ if (!childInfo) {
+ childInfo = new ISF::ParticleUserInformation();
+ bremPhot->setUserInformation(childInfo);
+ }
+ childInfo->setMaterialLimit(pLim.process,pLim.x0Max,0.);
+
+ // material fraction : flip if direction of propagation changed
+ float ci = m_layer->surfaceRepresentation().normal().dot( particleDir );
+ float co = m_layer->surfaceRepresentation().normal().dot( bremPhot->momentum().unit() );
+ float remMat = ci*co <0 ? (1-matFraction) : matFraction;
+
+ Trk::NeutralCurvilinearParameters nParm(bremPhot->position(),bremPhot->momentum(),bremPhot->charge());
+ Trk::ExtrapolationCell< Trk::NeutralParameters > enc(nParm);
+ enc.leadLayer = m_layer;
+
+ if (pLim.x0Max>0) {
+ enc.addConfigurationMode(Trk::ExtrapolationMode::StopWithMaterialLimitX0);
+ enc.materialLimitX0 = pLim.x0Max;
+ }
+ enc.materialProcess = pLim.process;
+ enc.materialX0 = 0.;
+
+ // loop back
+ processMaterialOnLayer(bremPhot,enc,dir,remMat);
- // get parent particle
- const ISF::ISFParticle *isp = ISF::ParticleClipboard::getInstance().getParticle();
- // something is seriously wrong if there is no parent particle
- assert(isp);
+}
- m_isp = isp;
+ISF::ISFParticle* iFatras::McMaterialEffectsEngine::bremPhoton(const ISF::ISFParticle* parent,
+ double time,
+ double pElectron,
+ double gammaE,
+ const Amg::Vector3D& vertex,
+ Amg::Vector3D& particleDir,Trk::ParticleHypothesis particle) const
+{
+ // ------------------------------------------------------
+ // simple approach
+ // (a) simulate theta uniform within the opening angle of the relativistic Hertz dipole
+ // theta_max = 1/gamma
+ // (b)Following the Geant4 approximation from L. Urban -> encapsulate that later
+ // the azimutal angle
+
+ double psi = 2.*M_PI*CLHEP::RandFlat::shoot(m_randomEngine);
+
+ // the start of the equation
+ double theta = 0.;
+ double m = m_particleMasses.mass[particle];
+ double E = sqrt(pElectron*pElectron + m*m);
+
+ if (m_uniformHertzDipoleAngle) {
+ // the simplest simulation
+ theta = m/E * CLHEP::RandFlat::shoot(m_randomEngine);
+ } else {
+ // ----->
+ theta = m/E;
+ // follow
+ double a = 0.625; // 5/8
+
+ double r1 = CLHEP::RandFlat::shoot(m_randomEngine);
+ double r2 = CLHEP::RandFlat::shoot(m_randomEngine);
+ double r3 = CLHEP::RandFlat::shoot(m_randomEngine);
+
+ double u = -log(r2*r3)/a;
+
+ theta *= (r1 < 0.25 ) ? u : u*m_oneOverThree; // 9./(9.+27) = 0.25
+ }
- //std::cout << "updateTrackParameters : processing isp:position within layer: "<<isp <<","<< isp->generation() <<","<< particle <<std::endl;
- //std::cout << "position: global, within layer:"<<parm->position()<<","<< matFraction<< std::endl;
+ ATH_MSG_VERBOSE("[ brem ]"<<"BremPhoton"<<""<<"Simulated angle to electron = " << theta << "." );
- // recalculate if missing
- m_matProp = m_matProp ? m_matProp : m_layer->fullUpdateMaterialProperties(*tParameters);
- double pathCorrection = m_layer->surfaceRepresentation().pathCorrection(tParameters->position(),dir*(tParameters->momentum()));
- //pathCorrection = pathCorrection > 0. ? pathCorrection : m_layer->pathCorrection(parm);
+ double th = particleDir.theta()-theta;
+ double ph = particleDir.phi();
+ if ( th<0.) { th *=-1; ph += M_PI; }
+ CLHEP::Hep3Vector newDirectionHep( sin(th)*cos(ph),sin(th)*sin(ph),cos(th) );
+ CLHEP::Hep3Vector particleDirHep( particleDir.x(), particleDir.y(), particleDir.z() );
+ newDirectionHep.rotate(psi,particleDirHep);
+ Amg::Vector3D newDirection( newDirectionHep.x(), newDirectionHep.y(), newDirectionHep.z() );
- // check if a bending correction has to be applied
- if (m_bendingCorrection) {
- ATH_MSG_WARNING("BendingCorrection not implemented!! (McMaterialEffectsEngine.cxx)");
- // const Trk::TrackingVolume* enclosingVolume = m_bendingCorrection ? m_layer->enclosingTrackingVolume() : 0;
- // if (enclosingVolume && enclosingVolume->bendingCorrector() ) {
- // double bendingCorrection = enclosingVolume->bendingCorrector()->bendingCorrection(*parm);
- // pathCorrection *= bendingCorrection;
- // }
- }
+ // recoil / save input momentum for validation
+ Amg::Vector3D inEl(pElectron*particleDir);
+ particleDir = (particleDir*pElectron- gammaE*newDirection).unit();
+
+ // -------> create the brem photon <--------------------
+ ISF::ISFParticle *bremPhoton = new ISF::ISFParticle( vertex,
+ gammaE*newDirection,
+ 0, //!< mass
+ 0, //!< charge
+ 22, //!< pdg code
+ time, //!< time
+ *parent );
- //--------------------------------------------------------------------------------------------------
- if (msgLvl(MSG::VERBOSE) && int(dir)){
- const Trk::TrackingVolume* enclosingVolume = m_layer->enclosingTrackingVolume();
- std::string volumeName = enclosingVolume ? enclosingVolume->volumeName() : "Unknown";
- double layerR = m_layer->surfaceRepresentation().bounds().r();
- double layerZ = m_layer->surfaceRepresentation().center().z();
- ATH_MSG_VERBOSE( " [M] Material effects update() on layer with index "<< m_layer->layerIndex() );
- ATH_MSG_VERBOSE( " -> path/X0 on layer [r/z] = [ " << layerR << " / " << layerZ << " ] :"
- << pathCorrection*m_matProp->thicknessInX0() );
- ATH_MSG_VERBOSE( " -> path correctin factor = " << pathCorrection );
+ // register TruthIncident
+ ISF::ISFParticleVector children(1, bremPhoton);
+ ISF::ISFTruthIncident truth( const_cast<ISF::ISFParticle&>(*parent),
+ children,
+ m_processCode,
+ parent->nextGeoID(),
+ ISF::fPrimarySurvives );
+ m_truthRecordSvc->registerTruthIncident(truth);
+
+ if (m_validationMode) {
+ ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
+ validInfo->setProcess(m_processCode);
+ if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
+ else validInfo->setGeneration(1); // assume parent is a primary track
+ bremPhoton->setUserInformation(validInfo);
+
+ // save interaction info
+ if ( m_validationTool ) {
+ Amg::Vector3D* nMom = new Amg::Vector3D(inEl);
+ m_validationTool->saveISFVertexInfo(m_processCode, vertex,*parent,inEl+gammaE*newDirection,nMom,children);
+ delete nMom;
+ }
+
+
}
- //--------------------------------------------------------------------------------------------------
+
+ return bremPhoton;
+
+}
+
+Trk::ExtrapolationCode iFatras::McMaterialEffectsEngine::processMaterialOnLayer( const ISF::ISFParticle* isp,
+ Trk::ExCellCharged& eCell,
+ Trk::PropDirection dir,
+ float& mFraction) const
+{
+
+ const Trk::TrackParameters* parm = eCell.leadParameters;
+
+ // path correction
+ double pathCorrection = fabs(m_layer->surfaceRepresentation().pathCorrection(parm->position(),dir*(parm->momentum())));
// figure out if particle stopped in the layer and recalculate path limit
int iStatus = 0;
- double dX0 = (1.-matFraction)*pathCorrection*m_matProp->thicknessInX0();
+ float dX0 = (1.-mFraction)*pathCorrection*m_thicknessInX0;
+ float dL0 = (1.-mFraction)*pathCorrection*m_thicknessInL0;
+
+ if ( eCell.materialLimitX0>0 && eCell.materialProcess<100 && eCell.materialX0+dX0>= eCell.materialLimitX0) { // elmg. interaction
+ float x0rem = eCell.materialLimitX0 - eCell.materialX0;
+ dL0 *= x0rem>0 ? x0rem/dX0 : 1.;
+ if ( x0rem>0 ) dX0 = x0rem;
+ iStatus = 1;
+ } else if ( eCell.materialLimitL0>0 && eCell.materialProcess>100 && eCell.materialL0+dL0 >= eCell.materialLimitL0 ) { // hadronic interaction
+ float l0rem = eCell.materialLimitX0 - eCell.materialL0 ;
+ dX0 *= l0rem>0 ? l0rem/dL0 : 1.;
+ if ( l0rem>0 ) dL0 = l0rem;
+ iStatus = 2;
+ }
+ // material update
+ eCell.materialX0 += dX0;
+ eCell.materialL0 += dL0;
+ eCell.zOaTrX += m_matProp->zOverAtimesRho()*dX0;
+ eCell.zX += m_matProp->averageZ()*dX0;
// get the kinematics
- double p = (tParameters->momentum()).mag();
+ double p = parm->momentum().mag();
+ double m = m_particleMasses.mass[eCell.pHypothesis];
+ double E = sqrt(p*p+m*m);
- matFraction += dX0/pathCorrection/m_matProp->thicknessInX0();
+ // radiation and ionization preceed the presampled interaction (if any)
- if ( isp->charge()!=0 ) {
- AmgVector(5) updatedParameters(tParameters->parameters());
- if ( m_eLoss ) ionize(*tParameters, updatedParameters, dX0, dir, eCell.pHypothesis );
+ if ( m_eLoss ) {
+ // the updatedParameters - first a deep copy
+ AmgVector(5) uParameters = parm->parameters();
- if ( m_ms ) {
+ // smeared/presampled energy loss
+ Trk::EnergyLoss* eloss = m_eLossSampler->energyLoss(*m_matProp, p, dX0/m_matProp->thicknessInX0(), dir, eCell.pHypothesis);
- double sigmaMSproj = (m_use_msUpdator && m_msUpdator ) ?
- sqrt(m_msUpdator->sigmaSquare(*m_matProp, p, dX0/m_matProp->thicknessInX0(), eCell.pHypothesis)) : msSigma(dX0,p,eCell.pHypothesis);
-
- multipleScatteringUpdate(*tParameters,updatedParameters,sigmaMSproj, pathCorrection);
+ if (eCell.pHypothesis==Trk::electron && m_createBremPhoton) {
+
+ // ionization update
+
+ double newP = E+eloss->meanIoni()>m ? sqrt((E+eloss->meanIoni())*(E+eloss->meanIoni())-m*m) : 0.5*m_minimumMomentum;
+ uParameters[Trk::qOverP] = parm->charge()/newP;
+
+ // radiation
+ if (newP>m_minimumMomentum)
+ radiate(isp,uParameters,eCell,dX0,mFraction,dir,pathCorrection*m_thicknessInX0); // mFraction used to estimate material thickness for brem photons
+
+ // save the actual radiation loss
+ float nqOp = uParameters[Trk::qOverP];
+ float radLoss = fabs(1./nqOp) - newP;
+ eloss->update(0.,0.,radLoss-eloss->meanRad(),eloss->meanRad()-radLoss);
+
+ } else {
+
+ // calculate the new momentum
+ double newP = E+eloss->deltaE()>m ? sqrt((E+eloss->deltaE())*(E+eloss->deltaE())-m*m) : 0.5*m_minimumMomentum;
+ uParameters[Trk::qOverP] = parm->charge()/newP;
+
+ }
+
+ // TODO straggling
+
+ if (m_validationMode && eloss) {
+ if (eCell.eLoss) {
+ eCell.eLoss->update(*eloss);
+ delete eloss;
+ } else
+ eCell.eLoss=eloss;
}
+
+
+ const Trk::TrackParameters* upd = parm->associatedSurface().createTrackParameters(uParameters[0],
+ uParameters[1],
+ uParameters[2],
+ uParameters[3],
+ uParameters[4]);
+ //eCell.updateLeadParameters(upd);
+ eCell.leadParameters = upd;
+
+ if (upd->momentum().mag() < m_minimumMomentum ) {
+
+ // save info for locally created particle
+ if (m_validationMode && isp!=m_isp) {
+ ATH_MSG_VERBOSE( " saving interaction info for locally produced particle " << isp->pdgCode() );
+ m_validationTool->saveISFParticleInfo(*isp,eCell,Trk::ExtrapolationCode::SuccessMaterialLimit);
+ }
+
+ if (isp!=m_isp) { delete isp; }
+ return Trk::ExtrapolationCode::FailureUpdateKill;
+ }
+ }
+
+ mFraction += dX0/pathCorrection/m_thicknessInX0;
+
+ if ( m_ms ) {
+ AmgVector(5) updatedParameters=eCell.leadParameters->parameters();
+
+ double simTheta = m_msSampler->simTheta(*m_matProp, p, dX0/m_thicknessInX0, eCell.pHypothesis);
+ //do the update -> You need 2 evaluation of simTheta. The second one is used to calculate deltaphi in multipleScatteringUpdate
+ multipleScatteringUpdate(*(eCell.leadParameters), updatedParameters, simTheta,
+ m_msSampler->simTheta(*m_matProp, p, dX0/m_thicknessInX0, eCell.pHypothesis));
+
+
// use the manipulator to update the track parameters -------> get rid of 0!
- tParameters = tParameters->associatedSurface().createTrackParameters(updatedParameters[0],
- updatedParameters[1],
- updatedParameters[2],
- updatedParameters[3],
- updatedParameters[4],
- 0);
+ const Trk::TrackParameters* upd = eCell.leadParameters->associatedSurface().createTrackParameters(updatedParameters[0],
+ updatedParameters[1],
+ updatedParameters[2],
+ updatedParameters[3],
+ updatedParameters[4]);
+
+ //eCell.updateLeadParameters(upd);
+ eCell.leadParameters=upd;
+ if (upd->momentum().mag() < m_minimumMomentum ) {
+ if (isp!=m_isp) { delete isp;}
+ return Trk::ExtrapolationCode::FailureUpdateKill;
+ }
+ }
+
+ if ( iStatus==1 || iStatus==2 ) { // interaction with particle stopping
+
+ ISF::ISFParticleVector childs = m_samplingTool->interact(isp,eCell); // TODO pass material when needed
+
+ // save info for locally created particles
+ if (m_validationMode && isp!=m_isp) {
+ ATH_MSG_VERBOSE( " saving interaction info for locally produced particle " << isp->pdgCode() );
+ m_validationTool->saveISFParticleInfo(*isp,eCell,Trk::ExtrapolationCode::SuccessMaterialLimit);
+ }
+
+ // loop over childrens and decide about their fate
- if ((tParameters->momentum()).mag() < m_minimumMomentum ) {
- if (isp!=m_isp) { delete isp; delete tParameters; }
- return 0;
+ for (unsigned ic=0; ic<childs.size(); ic++) {
+ double mom = childs[ic]->momentum().mag();
+
+ if (mom<m_minimumMomentum) continue;
+ Trk::ParticleHypothesis pHypothesis = m_pdgToParticleHypothesis.convert(childs[ic]->pdgCode(),childs[ic]->charge());
+ Trk::PathLimit pLim = m_samplingTool->sampleProcess(mom,childs[ic]->charge(),pHypothesis);
+ double pathLimit = ( !m_particleDecayHelper.empty()) ? m_particleDecayHelper->freePath(*childs[ic]) : -1.;
+
+ // TODO sample decays and save the material collection & path limits at the exit from the layer (ISFFatrasParticle ?)
+
+ // material fraction : flip if direction of propagation changed
+ float ci = m_layer->surfaceRepresentation().normal().dot( parm->momentum().unit() );
+ float co = m_layer->surfaceRepresentation().normal().dot( childs[ic]->momentum().unit() );
+ float remMat = ci*co <0 ? (1-mFraction) : mFraction;
+
+ // save presampled information
+ ISF::ParticleUserInformation* childInfo = childs[ic]->getUserInformation();
+ if (!childInfo) {
+ childInfo = new ISF::ParticleUserInformation();
+ childs[ic]->setUserInformation(childInfo);
+ }
+ childInfo->setMaterialLimit(pLim.process,pLim.x0Max,0.);
+ //std::cout <<"child info saved on isp:"<<childs[ic]<<","<<childs[ic]->getUserInformation()->materialLimit()<< std::endl;
+
+ // configure child cell and loop back
+ if (childs[ic]->charge()!=0) {
+ Trk::CurvilinearParameters cParm(childs[ic]->position(),childs[ic]->momentum(),childs[ic]->charge());
+ Trk::ExtrapolationCell< Trk::TrackParameters > ecc(cParm);
+
+ if (pathLimit>0) {
+ ecc.addConfigurationMode(Trk::ExtrapolationMode::StopWithPathLimit);
+ ecc.pathLimit = pathLimit;
+ }
+
+ if (pLim.x0Max>0) {
+ if (pLim.process==121) {
+ ecc.addConfigurationMode(Trk::ExtrapolationMode::StopWithMaterialLimitL0);
+ ecc.materialLimitL0 = pLim.x0Max;
+ } else {
+ ecc.addConfigurationMode(Trk::ExtrapolationMode::StopWithMaterialLimitX0);
+ ecc.materialLimitX0 = pLim.x0Max;
+ }
+ ecc.materialProcess = pLim.process;
+ }
+ ecc.leadLayer = m_layer;
+
+ // loop back
+ Trk::ExtrapolationCode exitCode=processMaterialOnLayer(childs[ic],ecc,dir,remMat);
+
+ if (exitCode == Trk::ExtrapolationCode::SuccessMaterialLimit )
+ ATH_MSG_VERBOSE(" child particle:"<< ic << "stopped in layer" );
+ else
+ ATH_MSG_VERBOSE(" child particle:"<< ic << "leaving layer" );
+
+ } else {
+ Trk::NeutralCurvilinearParameters nParm(childs[ic]->position(),childs[ic]->momentum(),childs[ic]->charge());
+ Trk::ExtrapolationCell< Trk::NeutralParameters > enc(nParm);
+
+ if (pathLimit>0) {
+ enc.addConfigurationMode(Trk::ExtrapolationMode::StopWithPathLimit);
+ enc.pathLimit = pathLimit;
+ }
+
+ if (pLim.x0Max>0) {
+ if (pLim.process==121) {
+ enc.addConfigurationMode(Trk::ExtrapolationMode::StopWithMaterialLimitL0);
+ enc.materialLimitL0 = pLim.x0Max;
+ } else {
+ enc.addConfigurationMode(Trk::ExtrapolationMode::StopWithMaterialLimitX0);
+ enc.materialLimitX0 = pLim.x0Max;
+ }
+ enc.materialProcess = pLim.process;
+ }
+ enc.leadLayer = m_layer;
+
+ // loop back
+ Trk::ExtrapolationCode exitCode=processMaterialOnLayer(childs[ic],enc,dir,remMat);
+
+ if (exitCode == Trk::ExtrapolationCode::SuccessMaterialLimit )
+ ATH_MSG_VERBOSE(" child particle:"<< ic << "stopped in layer" );
+ else
+ ATH_MSG_VERBOSE(" child particle:"<< ic << "leaving layer" );
+ }
}
+ if (childs.size()>0 || true ) { // TODO: pass the interaction status
+ if (isp!=m_isp) { delete isp; }
+ return Trk::ExtrapolationCode::SuccessMaterialLimit;
+ }
}
- // register particle if not in the stack already
- if (tParameters && isp!=m_isp ) {
- ISF::ISFParticle* regisp=new ISF::ISFParticle(isp->position(),tParameters->momentum(),isp->mass(),isp->charge(),
+ // register particle if not in the stack already
+ if (isp!=m_isp ) {
+ ISF::ISFParticle* regisp=new ISF::ISFParticle(isp->position(),parm->momentum(),isp->mass(),isp->charge(),
isp->pdgCode(),isp->timeStamp(),*m_isp,isp->barcode());
+ // add presampled process info
+ if (isp->getUserInformation() && isp->getUserInformation()->materialLimit()) {
+ ISF::MaterialPathInfo* matLim = isp->getUserInformation()->materialLimit();
+ ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
+ validInfo->setMaterialLimit(matLim->process,matLim->dMax,matLim->process==121 ? eCell.materialL0 : eCell.materialX0);
+ regisp->setUserInformation(validInfo);
+ //std::cout <<"collected material at layer exit:"<<regisp->getUserInformation()->materialLimit()->dCollected<<std::endl;
+ }
+ // validation mode
+ if (m_validationMode) {
+ ISF::ParticleUserInformation* validInfo = regisp->getUserInformation();
+ if (!validInfo){
+ validInfo = new ISF::ParticleUserInformation();
+ regisp->setUserInformation(validInfo);
+ }
+ if (isp->getUserInformation()) validInfo->setProcess(isp->getUserInformation()->process());
+ else validInfo->setProcess(-2); // signal problem in the validation chain
+ if (isp->getUserInformation()) validInfo->setGeneration(isp->getUserInformation()->generation());
+ else validInfo->setGeneration(-1); // signal problem in the validation chain
+ }
m_particleBroker->push(regisp, m_isp);
}
-
if (isp!=m_isp) delete isp;
+
+ return Trk::ExtrapolationCode::InProgress;
+}
+
+Trk::ExtrapolationCode iFatras::McMaterialEffectsEngine::processMaterialOnLayer( const ISF::ISFParticle* isp,
+ Trk::ExCellNeutral& eCell,
+ Trk::PropDirection dir,
+ float& mFraction) const
+{
+
+ const Trk::NeutralParameters* parm = eCell.leadParameters;
+
+ // path correction
+ double pathCorrection = fabs(m_layer->surfaceRepresentation().pathCorrection(parm->position(),dir*(parm->momentum())));
- if (tParameters) std::cout <<"layer "<< tParameters->position()<< " returns :"<< isp <<" with momentum "<< tParameters->momentum().mag() << std::endl;
- else std::cout <<"layer returns :"<< 0 << std::endl;
+ // figure out if particle stopped in the layer and recalculate path limit
+ int iStatus = 0;
+ float dX0 = (1.-mFraction)*pathCorrection*m_thicknessInX0;
+ float dL0 = (1.-mFraction)*pathCorrection*m_thicknessInL0;
- return (&parm);
+ if ( eCell.materialLimitX0>0 && eCell.materialProcess<100 && eCell.materialX0+dX0>= eCell.materialLimitX0) { // elmg. interaction
+ float x0rem = eCell.materialLimitX0 - eCell.materialX0;
+ dL0 *= x0rem>0 ? x0rem/dX0 : 1.;
+ if ( x0rem>0 ) dX0 = x0rem;
+ iStatus = 1;
+ } else if ( eCell.materialLimitL0>0 && eCell.materialProcess>100 && eCell.materialL0+dL0 >= eCell.materialLimitL0 ) { // hadronic interaction
+ float l0rem = eCell.materialLimitX0 - eCell.materialL0 ;
+ dX0 *= l0rem>0 ? l0rem/dL0 : 1.;
+ if ( l0rem>0 ) dL0 = l0rem;
+ iStatus = 2;
+ }
+ // material update
+ eCell.materialX0 += dX0;
+ eCell.materialL0 += dL0;
+ eCell.zOaTrX += m_matProp->zOverAtimesRho()*dX0;
+ eCell.zX += m_matProp->averageZ()*dX0;
-}
-
-void iFatras::McMaterialEffectsEngine::ionize(const Trk::TrackParameters& parm,
- AmgVector(5)& updatedPar,
- double dInX0,
- Trk::PropDirection /*dir*/,
- Trk::ParticleHypothesis particle ) const
-{
- double p = parm.momentum().mag();
- double m = m_particleMasses.mass[particle];
- double E = sqrt(p*p+m*m);
+ // get the kinematics
+ //double p = parm->momentum().mag();
+ //double m = m_particleMasses.mass[eCell.pHypothesis];
+ //double E = sqrt(p*p+m*m);
+
+ if ( iStatus==1 || iStatus==2 ) { // interaction with particle stopping
+
+ ISF::ISFParticleVector childs = m_samplingTool->interact(isp,eCell); // TODO pass material when needed
+
+ // save info for locally created particles
+ if (m_validationMode && isp!=m_isp) {
+ ATH_MSG_VERBOSE( " saving interaction info for locally produced particle " << isp->pdgCode() );
+ m_validationTool->saveISFParticleInfo(*isp,eCell,Trk::ExtrapolationCode::SuccessMaterialLimit);
+ }
+
+ // loop over childrens and decide about their fate
+
+ for (unsigned ic=0; ic<childs.size(); ic++) {
+ double mom = childs[ic]->momentum().mag();
+
+ if (mom<m_minimumMomentum) continue;
+ Trk::ParticleHypothesis pHypothesis = m_pdgToParticleHypothesis.convert(childs[ic]->pdgCode(),childs[ic]->charge());
+ Trk::PathLimit pLim = m_samplingTool->sampleProcess(mom,childs[ic]->charge(),pHypothesis);
+ double pathLimit = ( !m_particleDecayHelper.empty()) ? m_particleDecayHelper->freePath(*childs[ic]) : -1.;
+
+ // TODO sample decays and save the material collection & path limits at the exit from the layer (ISFFatrasParticle ?)
+
+ // material fraction : flip if direction of propagation changed
+ float ci = m_layer->surfaceRepresentation().normal().dot( parm->momentum().unit() );
+ float co = m_layer->surfaceRepresentation().normal().dot( childs[ic]->momentum().unit() );
+ float remMat = ci*co <0 ? (1-mFraction) : mFraction;
+
+ // save presampled information
+ ISF::ParticleUserInformation* childInfo = childs[ic]->getUserInformation();
+ if (!childInfo) {
+ childInfo = new ISF::ParticleUserInformation();
+ childs[ic]->setUserInformation(childInfo);
+ }
+ childInfo->setMaterialLimit(pLim.process,pLim.x0Max,0.);
+
+ // configure child cell and loop back
+ if (childs[ic]->charge()!=0) {
+ Trk::CurvilinearParameters cParm(childs[ic]->position(),childs[ic]->momentum(),childs[ic]->charge());
+ Trk::ExtrapolationCell< Trk::TrackParameters > ecc(cParm);
+
+ if (pathLimit>0) {
+ ecc.addConfigurationMode(Trk::ExtrapolationMode::StopWithPathLimit);
+ ecc.pathLimit = pathLimit;
+ }
+
+ if (pLim.x0Max>0) {
+ if (pLim.process==121) {
+ ecc.addConfigurationMode(Trk::ExtrapolationMode::StopWithMaterialLimitL0);
+ ecc.materialLimitL0 = pLim.x0Max;
+ } else {
+ ecc.addConfigurationMode(Trk::ExtrapolationMode::StopWithMaterialLimitX0);
+ ecc.materialLimitX0 = pLim.x0Max;
+ }
+ ecc.materialProcess = pLim.process;
+ }
+
+ // loop back
+ Trk::ExtrapolationCode exitCode = processMaterialOnLayer(childs[ic],ecc,dir,remMat);
- // the following formulas are imported from STEP
- // preparation of kinetic constants
-
- double me = m_particleMasses.mass[Trk::electron];
- double mfrac = me/m;
- double beta = p/E;
- double gamma = E/m;
- double Ionization = 0.;
-
- //Ionization - Bethe-Bloch
- double I = 16.e-6 * std::pow(m_matProp->averageZ(),0.9); //16 eV * Z**0.9 - bring to MeV
-
- //K/A*Z = 0.5 * 30.7075MeV/(g/mm2) * Z/A * rho[g/mm3] / scale to mm by this
- double kaz = 0.5*30.7075*m_matProp->zOverAtimesRho();
-
- if (particle == Trk::electron){
- // for electrons use slightly different BetheBloch adaption
- // see Stampfer, et al, "Track Fitting With Energy Loss", Comp. Pyhs. Comm. 79 (1994), 157-164
- Ionization = -kaz*(2.*log(2.*me/I)+3.*log(gamma) - 1.95);
+ if (exitCode == Trk::ExtrapolationCode::SuccessMaterialLimit )
+ ATH_MSG_VERBOSE(" child particle:"<< ic << "stopped in layer" );
+ else
+ ATH_MSG_VERBOSE(" child particle:"<< ic << "leaving layer" );
+
+ } else {
+ Trk::NeutralCurvilinearParameters nParm(childs[ic]->position(),childs[ic]->momentum(),childs[ic]->charge());
+ Trk::ExtrapolationCell< Trk::NeutralParameters > enc(nParm);
+
+ if (pathLimit>0) {
+ enc.addConfigurationMode(Trk::ExtrapolationMode::StopWithPathLimit);
+ enc.pathLimit = pathLimit;
+ }
+
+ if (pLim.x0Max>0) {
+ if (pLim.process==121) {
+ enc.addConfigurationMode(Trk::ExtrapolationMode::StopWithMaterialLimitL0);
+ enc.materialLimitL0 = pLim.x0Max;
+ } else {
+ enc.addConfigurationMode(Trk::ExtrapolationMode::StopWithMaterialLimitX0);
+ enc.materialLimitX0 = pLim.x0Max;
+ }
+ enc.materialProcess = pLim.process;
+ }
+
+ // loop back
+ Trk::ExtrapolationCode exitCode = processMaterialOnLayer(childs[ic],enc,dir,remMat);
+
+ if (exitCode == Trk::ExtrapolationCode::SuccessMaterialLimit )
+ ATH_MSG_VERBOSE(" child particle:"<< ic << "stopped in layer" );
+ else
+ ATH_MSG_VERBOSE(" child particle:"<< ic << "leaving layer" );
+ }
+ }
+
+ if (childs.size()>0 || true) { // TODO : pass the interaction status
+ if (isp!=m_isp) { delete isp; }
+ return Trk::ExtrapolationCode::SuccessMaterialLimit;
+ }
}
- else {
- double eta2 = beta*gamma; eta2 *= eta2;
- // density effect, only valid for high energies (gamma > 10 -> p > 1GeV for muons)
- double delta = 0.;
- if (gamma > 10.) {
- double eplasma = 28.816e-6 * sqrt(1000.*m_matProp->zOverAtimesRho());
- delta = 2.*log(eplasma/I) + log(eta2) - 1.;
+
+ // register particle if not in the stack already
+ if (isp!=m_isp ) {
+ ISF::ISFParticle* regisp=new ISF::ISFParticle(isp->position(),parm->momentum(),isp->mass(),isp->charge(),
+ isp->pdgCode(),isp->timeStamp(),*m_isp,isp->barcode());
+ // add presampled process info
+ if (isp->getUserInformation() && isp->getUserInformation()->materialLimit()) {
+ ISF::MaterialPathInfo* matLim = isp->getUserInformation()->materialLimit();
+ ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
+ validInfo->setMaterialLimit(matLim->process,matLim->dMax,matLim->process==121 ? eCell.materialL0 : eCell.materialX0);
+ regisp->setUserInformation(validInfo);
+ }
+ // validation mode
+ if (m_validationMode) {
+ ISF::ParticleUserInformation* validInfo = regisp->getUserInformation();
+ if (!validInfo){
+ validInfo = new ISF::ParticleUserInformation();
+ regisp->setUserInformation(validInfo);
+ }
+ if (isp->getUserInformation()) validInfo->setProcess(isp->getUserInformation()->process());
+ else validInfo->setProcess(-2); // signal problem in the validation chain
+ if (isp->getUserInformation()) validInfo->setGeneration(isp->getUserInformation()->generation());
+ else validInfo->setGeneration(-1); // signal problem in the validation chain
}
- // tmax - cut off energy
- double tMax = 2.*eta2*me/(1.+2.*gamma*mfrac+mfrac*mfrac);
- // divide by beta^2 for non-electrons
- kaz /= beta*beta;
- Ionization = -kaz*(log(2.*me*eta2*tMax/(I*I)) - 2.*(beta*beta) - delta);
+ m_particleBroker->push(regisp, m_isp);
}
-
- double energyLoss = Ionization*dInX0*m_matProp->x0();
- double newP = sqrt(fmax(0.,(E+energyLoss)*(E+energyLoss)-m*m));
-
- newP = newP <= 0.5*m_minimumMomentum ? 0.5*m_minimumMomentum : newP; // arbitrary regularization;
+ if (isp!=m_isp) delete isp;
- (updatedPar)[Trk::qOverP] = parm.charge()/newP;
-
+ return Trk::ExtrapolationCode::InProgress;
}
-double iFatras::McMaterialEffectsEngine::msSigma(double dInX0,double p,Trk::ParticleHypothesis particle) const {
-
- double m = m_particleMasses.mass[particle];
- double E = sqrt(p*p+m*m);
- double beta = p/E;
- // PDG Particle Review, Phys.Rev.D86,010001(2012), chapter 30.3, page 328
+// radiative effects
+void iFatras::McMaterialEffectsEngine::radiate(const ISF::ISFParticle* parent, AmgVector(5)& parm ,
+ Trk::ExCellCharged& eCell, float pathLim, float mFr,
+ Trk::PropDirection dir, float refX) const {
+
+ // sample energy loss and free path independently
+ double path = 0.;
+ double p = 1./ fabs(parm[Trk::qOverP]);
- double sigmaMSproj = 13.6/beta/p*sqrt(dInX0)*(1.+0.038*log(dInX0));
+ Amg::Vector3D eDir = eCell.leadParameters->momentum().unit();
+ Amg::Vector3D ePos = eCell.leadParameters->position();
- return sigmaMSproj;
-}
-
-void iFatras::McMaterialEffectsEngine::multipleScatteringUpdate(const Trk::TrackParameters& pars,
- AmgVector(5)& parameters,
- double sigmaMSproj, double pathCorrection) const
-{
- //block of initial values
- //------------------------------------//
- double d = 0.1*m_matProp->thickness(); //layer thickness, cm
- d*=pathCorrection;
- double X0 = 0.1*m_matProp->x0(); //layer rad len, cm
- double Z = m_matProp->averageZ(); //charge layer material
- //double AN = 2*Z; //mass number of the material (an approximation, bad for heavy nuclei)
- double p = m_isp->momentum().mag()*0.001; //particle momentum, GeV
- double m = m_isp->mass()*0.001; //particle mass, GeV
- double dOverX0 = d/X0;
- double E=sqrt(p*p+m*m);
- double beta=p/E;
- double scale = 0.860176; ///magic
-
- std::cout<<"McMaterialEffectsEngine::multipleScatteringUpdate: Material pars (d,beta,p,dOverX0,Z): "<<std::endl;
- std::cout<<d<<' '<<beta<<' '<<p<<' '<<dOverX0<<' '<<Z<<std::endl;
- //-----------------------------------//
-
- double theta = simulate_theta_space(beta, p, dOverX0, Z, scale);
- //theta*=m_projectionFactor*sigmaMSproj;
-
- double psi = 2*M_PI*CLHEP::RandFlat::shoot(m_randomEngine);
-
- CLHEP::Hep3Vector parsMomHep( pars.momentum().x(), pars.momentum().y(), pars.momentum().z() );
- CLHEP::Hep3Vector newDirection( parsMomHep.unit().x(), parsMomHep.unit().y(), parsMomHep.unit().z());
- double x = -newDirection.y();
- double y = newDirection.x();
- double z = 0.;
- // if it runs along the z axis - no good ==> take the x axis
- if (newDirection.z()*newDirection.z() > 0.999999) {
- x=1.; y=0.;
- }
- // deflector direction
- CLHEP::Hep3Vector deflector(x,y,z);
- // rotate the new direction for scattering
- newDirection.rotate(theta, deflector);
- // and arbitrarily in psi
- newDirection.rotate(psi, parsMomHep);
- // assign the new values
- parameters[Trk::phi] = newDirection.phi();
- parameters[Trk::theta] = newDirection.theta();
-
-}
+ while ( path < pathLim && p>m_minimumMomentum ) {
-double iFatras::McMaterialEffectsEngine::simulate_theta_space(double beta, double p,double dOverX0,double Z, double scale) const{
- double * scattering_params;
- // Decide which mixture is best
- if (dOverX0/(beta*beta)>0.6/pow(Z,0.6)){ //Gaussian
- // Gaussian mixture or pure Gaussian
- if (dOverX0/(beta*beta)>10){
- scattering_params=get_gaussian(beta,p,dOverX0,scale); // Get parameters
- std::cout<<"McMaterialEffectsEngine::multipleScatteringUpdate: using pure_gaussian"<<std::endl;
- }
- else{
- scattering_params=get_gaussmix(beta,p,dOverX0,Z,scale); // Get parameters
- std::cout<<"McMaterialEffectsEngine::multipleScatteringUpdate: using gaussian_mixture"<<std::endl;
- }
- return sim_gaussmix(scattering_params); // Simulate
- }
- //Semigaussian mixture
- scattering_params = get_semigauss(beta,p,dOverX0,Z,scale); // Get parameters
- std::cout<<"McMaterialEffectsEngine::multipleScatteringUpdate: using semi_gaussian mixture"<<std::endl;
- return sim_semigauss(scattering_params); // Simulate
-}
+ double rndx = CLHEP::RandFlat::shoot(m_randomEngine);
-double * iFatras::McMaterialEffectsEngine::get_gaussian(double beta, double p,double dOverX0, double scale) const{
- double * scattering_params = new double[4];
- scattering_params[0]=0.015/beta/p*sqrt(dOverX0)*scale; //Total standard deviation of mixture
- // scattering_params[0]=1.0;
- scattering_params[1]=1.0; //Variance of core
- scattering_params[2]=1.0; //Variance of tails
- scattering_params[3]=0.5; //Mixture weight of tail component
- return scattering_params;
-}
+ // sample visible fraction of the mother momentum taken according to 1/f
+
+ double eps = fmin(10.,m_minimumMomentum)/p;
-double * iFatras::McMaterialEffectsEngine::get_gaussmix(double beta, double p,double dOverX0,double Z, double scale) const{
- double * scattering_params = new double[4];
- scattering_params[0]=0.015/beta/p*sqrt(dOverX0)*scale; //Total standard deviation of mixture
-// scattering_params[0]=1.0;
- double d1=log(dOverX0/(beta*beta));
- double d2=log(pow(Z,2.0/3.0)*dOverX0/(beta*beta));
- double epsi;
- double var1=(-1.843e-3*d1+3.347e-2)*d1+8.471e-1; //Variance of core
- if (d2<0.5)
- epsi=(6.096e-4*d2+6.348e-3)*d2+4.841e-2;
- else
- epsi=(-5.729e-3*d2+1.106e-1)*d2-1.908e-2;
- scattering_params[1]=var1; //Variance of core
- scattering_params[2]=(1-(1-epsi)*var1)/epsi; //Variance of tails
- scattering_params[3]=epsi; //Mixture weight of tail component
- return scattering_params;
-}
+ double z = pow(eps,pow(rndx,exp(1.1))); // adjustment here ?
-double * iFatras::McMaterialEffectsEngine::get_semigauss(double beta,double p,double dOverX0,double Z, double scale) const{
- double * scattering_params = new double[6];
- double N=dOverX0*1.587E7*pow(Z,1.0/3.0)/(beta*beta)/(Z+1)/log(287/sqrt(Z));
- scattering_params[4]=0.015/beta/p*sqrt(dOverX0)*scale; //Total standard deviation of mixture
-// double eps = 0.125*log10(10*dOverX0);
-// scattering_params[4]=0.0175/beta/p*sqrt(dOverX0)*(1+eps);
- //scattering_params[4]= 0.0192/beta/p*sqrt(dOverX0)*(1.+0.038*log(dOverX0));
- //scattering_params[4]=1.0;
- double rho=41000/pow(Z,2.0/3.0);
- double b=rho/sqrt(N*(log(rho)-0.5));
- double n=pow(Z,0.1)*log(N);
- double var1=(5.783E-4*n+3.803E-2)*n+1.827E-1;
- double a=(((-4.590E-5*n+1.330E-3)*n-1.355E-2)*n+9.828E-2)*n+2.822E-1;
- double epsi = (1-var1)/(a*a*(log(b/a)-0.5)-var1);
- scattering_params[3]=(epsi>0) ? epsi : 0.0;//Mixture weight of tail component
- scattering_params[0]=a; //Parameter 1 of tails
- scattering_params[1]=b; //Parameter 2 of tails
- scattering_params[2]=var1; //Variance of core
- scattering_params[5]=N; //Average number of scattering processes
- return scattering_params;
-}
+ // convert into scaling factor for mother momentum
+ z = (1.- z);
-double iFatras::McMaterialEffectsEngine::sim_gaussmix(double * scattering_params) const{
- double sigma_tot = scattering_params[0];
- double var1 = scattering_params[1];
- double var2 = scattering_params[2];
- double epsi = scattering_params[3];
- bool ind = (double)rand()/RAND_MAX>epsi;
- double u=(double)rand()/RAND_MAX;
- if(ind)
- return sqrt(var1)*sqrt(-2*log(u))*sigma_tot;
- else
- return sqrt(var2)*sqrt(-2*log(u))*sigma_tot;
-}
+ // turn into path
+ //double dx = -log(z)*0.65; // adjust for mean of exp(-x)
+ double dx = -0.7*log(z); // adjust for mean of exp(-x)
+
+ // sample free path (in X0)
+ //double rnd = CLHEP::RandFlat::shoot(m_randomEngine);
+ //double dx = -log(rnd)/13.;
+
+ // resolve the case when there is not enough material left in the layer
+ if ( path+dx > pathLim ) {
+ double rndp = CLHEP::RandFlat::shoot(m_randomEngine);
+ if (rndp > (pathLim-path)/dx){
+ (parm)[Trk::qOverP] = (parm)[Trk::qOverP]>0 ? 1/p : -1/p;
+ mFr += (pathLim-path)/refX;
+ path = pathLim;
+ break; // radiation loop ended
+ }
+ path += dx*rndp;
+ mFr += dx*rndp/refX;
+
+ } else {
+ path+=dx;
+ mFr += dx/refX;
+ }
+ if ( p*(1-z) > m_minimumBremPhotonMomentum ) {
+
+ double deltaP = (1-z)*p;
+ recordBremPhotonLay(parent,eCell.time,p,deltaP,ePos,eDir,eCell.pHypothesis,dir,mFr);
+ p *=z ;
+
+ // std::cout <<"brem photon emitted, momentum update:"<< p<<","<<path<<","<<pathLim<<",mat.fraction:"<<mFr<<std::endl;
+ }
+ }
-double iFatras::McMaterialEffectsEngine::sim_semigauss(double * scattering_params) const{
- double a = scattering_params[0];
- double b = scattering_params[1];
- double var1 = scattering_params[2];
- double epsi = scattering_params[3];
- double sigma_tot = scattering_params[4];
- bool ind=(double)rand()/RAND_MAX>epsi;
- double u=(double)rand()/RAND_MAX;
- if(ind)
- return sqrt(var1)*sqrt(-2*log(u))*sigma_tot;
- else
- return a*b*sqrt((1-u)/(u*b*b+a*a))*sigma_tot;
+ (parm)[Trk::qOverP] = (parm)[Trk::qOverP] > 0 ? 1/p : -1./p;
+ (parm)[Trk::theta] = eDir.theta();
+ (parm)[Trk::phi] = eDir.phi();
+ return;
}
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsUpdator.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsUpdator.cxx
index 949c0320a11fcd53ffc70398770a987309f85455..eaca07196a0a8383fe0df9196d07aa5fe1d7af75 100755
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsUpdator.cxx
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/McMaterialEffectsUpdator.cxx
@@ -21,7 +21,7 @@
#include "ISF_Event/ISFTruthIncident.h"
#include "ISF_Event/ParticleClipboard.h"
#include "ISF_Event/ParticleUserInformation.h"
-#include "ISF_FatrasInterfaces/IParticleDecayer.h"
+#include "ISF_FatrasInterfaces/IParticleDecayHelper.h"
// iFatras
#include "ISF_FatrasInterfaces/IHadronicInteractionProcessor.h"
#include "ISF_FatrasInterfaces/IProcessSamplingTool.h"
@@ -140,7 +140,7 @@ iFatras::McMaterialEffectsUpdator::McMaterialEffectsUpdator(const std::string& t
declareProperty("ProcessSamplingTool" , m_samplingTool);
declareProperty("PhotonConversionTool" , m_conversionTool);
// tool handle for the particle decayer
- declareProperty( "ParticleDecayer" , m_particleDecayer );
+ declareProperty( "ParticleDecayHelper" , m_particleDecayer );
// MC Truth Properties
declareProperty("BremProcessCode" , m_processCode, "MCTruth Physics Process Code");
// the steering --------------------------------------------------------------
@@ -485,9 +485,20 @@ const Trk::TrackParameters* iFatras::McMaterialEffectsUpdator::updateInLay(const
if (parm && isp!=m_isp ) {
ISF::ISFParticle* regisp=new ISF::ISFParticle(isp->position(),parm->momentum(),isp->mass(),isp->charge(),
isp->pdgCode(),isp->timeStamp(),*m_isp,isp->barcode());
- // in the validation mode, add process info
- if (m_validationMode) {
+ // add presampled process info
+ if (isp->getUserInformation() && isp->getUserInformation()->materialLimit()) {
+ ISF::MaterialPathInfo* matLim = isp->getUserInformation()->materialLimit();
ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
+ validInfo->setMaterialLimit(matLim->process,matLim->dMax,matLim->process==121 ? pathLim.l0Collected : pathLim.x0Collected);
+ regisp->setUserInformation(validInfo);
+ }
+ // in the validation mode, add process info
+ if (m_validationMode) {
+ ISF::ParticleUserInformation* validInfo = regisp->getUserInformation();
+ if (!validInfo) {
+ validInfo = new ISF::ParticleUserInformation();
+ regisp->setUserInformation(validInfo);
+ }
if (isp->getUserInformation()) validInfo->setProcess(isp->getUserInformation()->process());
else validInfo->setProcess(-2); // signal problem in the validation chain
if (isp->getUserInformation()) validInfo->setGeneration(isp->getUserInformation()->generation());
@@ -575,7 +586,8 @@ const Trk::TrackParameters* iFatras::McMaterialEffectsUpdator::updateInLay(const
if ( iStatus==1 || iStatus==2 ) { // interaction
ISF::ISFParticleVector childs = iStatus==1 ? interactLay(isp,timeLim.time,*parm,particle,pathLim.process) :
- m_hadIntProcessor->doHadIntOnLayer(isp,timeLim.time,*parm,m_extMatProp, particle);
+ m_hadIntProcessor->doHadIntOnLayer(isp, timeLim.time, parm->position(), parm->momentum(),
+ ( m_extMatProp ? &m_extMatProp->material() : 0), particle);
// save info for locally created particles
if (m_validationMode && childs.size()>0 && isp!=m_isp) {
@@ -615,15 +627,25 @@ const Trk::TrackParameters* iFatras::McMaterialEffectsUpdator::updateInLay(const
if (parm && isp!=m_isp ) {
ISF::ISFParticle* regisp=new ISF::ISFParticle(isp->position(),parm->momentum(),isp->mass(),isp->charge(),
isp->pdgCode(),isp->timeStamp(),*m_isp,isp->barcode());
+ // add presampled process info
+ if (isp->getUserInformation() && isp->getUserInformation()->materialLimit()) {
+ ISF::MaterialPathInfo* matLim = isp->getUserInformation()->materialLimit();
+ ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
+ validInfo->setMaterialLimit(matLim->process,matLim->dMax,matLim->process==121 ? pathLim.l0Collected : pathLim.x0Collected);
+ regisp->setUserInformation(validInfo);
+ }
// in the validation mode, add process info
if (m_validationMode) {
- ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
+ ISF::ParticleUserInformation* validInfo = regisp->getUserInformation();
+ if (!validInfo){
+ validInfo = new ISF::ParticleUserInformation();
+ regisp->setUserInformation(validInfo);
+ }
if (isp->getUserInformation()) validInfo->setProcess(isp->getUserInformation()->process());
else validInfo->setProcess(-2); // signal problem in the validation chain
if (isp->getUserInformation()) validInfo->setGeneration(isp->getUserInformation()->generation());
else validInfo->setGeneration(-1); // signal problem in the validation chain
- regisp->setUserInformation(validInfo);
- }
+ }
m_particleBroker->push(regisp, m_isp);
}
@@ -860,7 +882,7 @@ const Trk::TrackParameters* iFatras::McMaterialEffectsUpdator::update(double tim
// Hadronic Interaction section ----------------------------------------------------------------------------------------------
// ST add hadronic interaction for secondaries
- if ( particle == Trk::pion && m_hadInt && m_hadIntProcessor->hadronicInteraction(parm, p, E,matprop, pathCorrection,particle)){
+ if ( particle == Trk::pion && m_hadInt && m_hadIntProcessor->hadronicInteraction(parm.position(), parm.momentum(), p, E, parm.charge(), matprop, pathCorrection, particle)){
ATH_MSG_VERBOSE( " [+] Hadronic interaction killed initial hadron ... stop simulation, tackle childs." );
return 0;
}
@@ -1500,6 +1522,7 @@ void iFatras::McMaterialEffectsUpdator::recordBremPhotonLay(const ISF::ISFPartic
timeLim.time, //!< time
*parent );
+
// in the validation mode, add process info
if (m_validationMode) {
ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
@@ -1612,7 +1635,7 @@ const Trk::TrackParameters* iFatras::McMaterialEffectsUpdator::interact(double
if ( process == 201 ) { // decay
// update parent before decay
- ISF::ISFParticleVector childVector = m_particleDecayer->decayParticle(*parent,position,momentum,parent->pdgCode(),time);
+ ISF::ISFParticleVector childVector = m_particleDecayer->decayParticle(*parent,position,momentum,time);
for (unsigned int i=0; i<childVector.size(); i++) {
// in the validation mode, add process info
@@ -1714,13 +1737,14 @@ const Trk::TrackParameters* iFatras::McMaterialEffectsUpdator::interact(double
if (process==14) { // photon conversion
- const Trk::TrackParameters* parm = new Trk::CurvilinearParameters(position,momentum,parent->charge());
+ const Trk::NeutralParameters* parm = new Trk::NeutralCurvilinearParameters(position,momentum,parent->charge());
bool cStat = m_conversionTool->doConversion(time, *parm);
- if ( cStat ) delete parm;
-
- return (cStat? 0 : parm);
+ if (!cStat) ATH_MSG_WARNING( "Conversion failed, killing photon anyway ");
+
+ // kill the mother particle
+ delete parm; return 0;
}
if (process==121) { // hadronic interaction
@@ -1729,7 +1753,7 @@ const Trk::TrackParameters* iFatras::McMaterialEffectsUpdator::interact(double
const Trk::TrackParameters* parm = new Trk::CurvilinearParameters(position,momentum,parent->charge());
- bool recHad = m_hadIntProcessor->doHadronicInteraction(time,*parm,m_extMatProp, particle, true);
+ bool recHad = m_hadIntProcessor->doHadronicInteraction(time, position, momentum, m_extMatProp, particle, true);
// eventually : bool recHad = m_hadIntProcessor->recordHadState( time, p, position, pDir, particle);
// kill the track if interaction recorded --------------------------
@@ -1823,7 +1847,29 @@ ISF::ISFParticleVector iFatras::McMaterialEffectsUpdator::interactLay(const ISF
if (process==14) { // photon conversion
- return m_conversionTool->doConversionOnLayer(parent, time, parm);
+ Trk::NeutralCurvilinearParameters neu(position,momentum,parent->charge());
+ childVector=m_conversionTool->doConversionOnLayer(parent, time, neu);
+
+ // validation mode
+ if (m_validationMode && m_validationTool) {
+
+ // add process info for children
+ for (unsigned int i=0; i<childVector.size(); i++) {
+ ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
+ validInfo->setProcess(process);
+ if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
+ else validInfo->setGeneration(1); // assume parent is a primary track
+ childVector[i]->setUserInformation(validInfo);
+ }
+ // save interaction info
+ if ( m_validationTool ) {
+ Amg::Vector3D* nMom = 0;
+ m_validationTool->saveISFVertexInfo(process, position,*parent,momentum,nMom,childVector);
+ delete nMom;
+ }
+ }
+
+ return childVector;
}
if (process==121) { // hadronic interaction
@@ -1832,12 +1878,13 @@ ISF::ISFParticleVector iFatras::McMaterialEffectsUpdator::interactLay(const ISF
const Trk::CurvilinearParameters parm(position,momentum,parent->charge());
- return ( m_hadIntProcessor->doHadIntOnLayer(parent, time,parm,m_extMatProp, particle) );
+ return ( m_hadIntProcessor->doHadIntOnLayer(parent, time, position, momentum,
+ m_extMatProp? &m_extMatProp->material() : 0, particle) );
}
if ( process == 201 ) { // decay
- childVector = m_particleDecayer->decayParticle(*parent,parm.position(),parm.momentum(),parent->pdgCode(),time);
+ childVector = m_particleDecayer->decayParticle(*parent,parm.position(),parm.momentum(),time);
// in the validation mode, add process info
if (m_validationMode) {
for (unsigned int i=0; i<childVector.size(); i++) {
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/MultipleScatteringSamplerGaussianMixture.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/MultipleScatteringSamplerGaussianMixture.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c7e84d57ec2262447e8d34f73aa33f0f9aa18c9e
--- /dev/null
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/MultipleScatteringSamplerGaussianMixture.cxx
@@ -0,0 +1,166 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// MultipleScatteringSamplerGaussianMixture.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+// Trk include
+#include "TrkGeometry/MaterialProperties.h"
+
+#include "ISF_FatrasTools/MultipleScatteringSamplerGaussianMixture.h"
+
+#include "AtlasCLHEP_RandomGenerators/RandGaussZiggurat.h"
+
+// static particle masses
+Trk::ParticleMasses iFatras::MultipleScatteringSamplerGaussianMixture::s_particleMasses;
+// static doubles
+double iFatras::MultipleScatteringSamplerGaussianMixture::s_main_RutherfordScott = 13.6*Gaudi::Units::MeV;
+double iFatras::MultipleScatteringSamplerGaussianMixture::s_log_RutherfordScott = 0.038;
+
+double iFatras::MultipleScatteringSamplerGaussianMixture::s_main_RossiGreisen = 17.5*Gaudi::Units::MeV;
+double iFatras::MultipleScatteringSamplerGaussianMixture::s_log_RossiGreisen = 0.125;
+
+// ============================= Gaussian mixture model =============
+double iFatras::MultipleScatteringSamplerGaussianMixture::s_gausMixSigma1_a0 = 8.471e-1;
+double iFatras::MultipleScatteringSamplerGaussianMixture::s_gausMixSigma1_a1 = 3.347e-2;
+double iFatras::MultipleScatteringSamplerGaussianMixture::s_gausMixSigma1_a2 = -1.843e-3;
+
+double iFatras::MultipleScatteringSamplerGaussianMixture::s_gausMixEpsilon_a0 = 4.841e-2;
+double iFatras::MultipleScatteringSamplerGaussianMixture::s_gausMixEpsilon_a1 = 6.348e-3;
+double iFatras::MultipleScatteringSamplerGaussianMixture::s_gausMixEpsilon_a2 = 6.096e-4;
+
+double iFatras::MultipleScatteringSamplerGaussianMixture::s_gausMixEpsilon_b0 = -1.908e-2;
+double iFatras::MultipleScatteringSamplerGaussianMixture::s_gausMixEpsilon_b1 = 1.106e-1;
+double iFatras::MultipleScatteringSamplerGaussianMixture::s_gausMixEpsilon_b2 = -5.729e-3;
+
+double iFatras::MultipleScatteringSamplerGaussianMixture::s_projectionFactor = sqrt(2.);
+
+// constructor
+iFatras::MultipleScatteringSamplerGaussianMixture::MultipleScatteringSamplerGaussianMixture(const std::string& t, const std::string& n, const IInterface* p) :
+ AthAlgTool(t,n,p),
+ m_log_include(true),
+ m_optGaussianMixtureG4(true),
+ m_rndGenSvc("AtRndmGenSvc", n),
+ m_randomEngine(0),
+ m_randomEngineName("TrkExRnd")
+{
+ declareInterface<IMultipleScatteringSampler>(this);
+ // multiple scattering parameters
+ declareProperty("MultipleScatteringLogarithmicTermOn", m_log_include);
+ declareProperty("G4OptimisedGaussianMixtureModel", m_optGaussianMixtureG4);
+ // random service for Gaussian mixture model
+ declareProperty("RandomNumberService", m_rndGenSvc, "Name of the random number service");
+ declareProperty("RandomStreamName", m_randomEngineName, "Name of the random number stream");
+}
+
+// destructor
+iFatras::MultipleScatteringSamplerGaussianMixture::~MultipleScatteringSamplerGaussianMixture()
+{}
+
+// Athena standard methods
+// initialize
+StatusCode iFatras::MultipleScatteringSamplerGaussianMixture::initialize()
+{
+ // get the random generator service
+ if (m_rndGenSvc.retrieve().isFailure()) {
+ ATH_MSG_WARNING( "Could not retrieve " << m_rndGenSvc );
+ } else {
+ ATH_MSG_VERBOSE( "Successfully retrieved " << m_rndGenSvc );
+ }
+
+ //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;
+ } else
+ ATH_MSG_VERBOSE ("Successfully got random engine '" << m_randomEngineName << "'");
+
+ ATH_MSG_VERBOSE( "Gaussian mixture model optimised = " << m_optGaussianMixtureG4 );
+
+ ATH_MSG_INFO( "initialize() successful" );
+ return StatusCode::SUCCESS;
+}
+
+// finalize
+StatusCode iFatras::MultipleScatteringSamplerGaussianMixture::finalize()
+{
+ ATH_MSG_INFO( "finalize() successful" );
+ return StatusCode::SUCCESS;
+}
+
+
+double iFatras::MultipleScatteringSamplerGaussianMixture::simTheta(const Trk::MaterialProperties& mat,
+ double p,
+ double pathcorrection,
+ Trk::ParticleHypothesis particle,
+ double ) const
+{
+ if (mat.thicknessInX0()<=0. || particle==Trk::geantino) return 0.;
+
+ // make sure the path correction is positive to avoid a floating point exception
+ pathcorrection *= pathcorrection < 0. ? (-1.) : (1) ;
+
+ // scale the path length to the radiation length
+ double t = pathcorrection * mat.thicknessInX0();
+
+ // kinematics (relativistic)
+ double m = s_particleMasses.mass[particle];
+ double E = sqrt(p*p + m*m);
+ double beta = p/E;
+
+ double sigma2(0.);
+
+ if (particle != Trk::electron) {
+
+ // the highland formula
+ sigma2 = s_main_RutherfordScott/(beta*p);
+
+ if (m_log_include)
+ sigma2 *= (1.+s_log_RutherfordScott*log(t));
+
+ sigma2 *= (sigma2*t);
+ }
+
+ else {
+
+ // Electron multiple scattering effects - see Highland NIM 129 (1975) p497-499
+ // (Highland extension to the Rossi-Greisen formulation)
+ // NOTE: The formula can be extended by replacing the momentum^2 term with pi * pf
+ sigma2 = s_main_RossiGreisen / ( beta * p );
+ sigma2 *= (sigma2*t);
+
+ if ( m_log_include ) {
+ double factor = 1. + s_log_RossiGreisen * log10( 10. * t );
+ factor *= factor;
+ sigma2 *= factor;
+ }
+ }
+
+ // d_0'
+ double dprime = t/(beta*beta);
+ double log_dprime = log(dprime);
+ // d_0''
+ double log_dprimeprime = log(std::pow(mat.averageZ(),2.0/3.0) * dprime);
+ // get epsilon
+ double epsilon = log_dprimeprime < 0.5 ?
+ s_gausMixEpsilon_a0 + s_gausMixEpsilon_a1*log_dprimeprime + s_gausMixEpsilon_a2*log_dprimeprime*log_dprimeprime :
+ s_gausMixEpsilon_b0 + s_gausMixEpsilon_b1*log_dprimeprime + s_gausMixEpsilon_b2*log_dprimeprime*log_dprimeprime;
+ // the standard sigma
+ double sigma1square = s_gausMixSigma1_a0 + s_gausMixSigma1_a1*log_dprime + s_gausMixSigma1_a2*log_dprime*log_dprime;
+ // G4 optimised / native double Gaussian model
+ if (!m_optGaussianMixtureG4) sigma2 = 225.*dprime/(p*p);
+ // throw the random number core/tail
+ if ( CLHEP::RandFlat::shoot(m_randomEngine) < epsilon) {
+ sigma2 *= (1.-(1.-epsilon)*sigma1square)/epsilon;
+ }
+
+ return s_projectionFactor*sqrt(sigma2)*CLHEP::RandGaussZiggurat::shoot(m_randomEngine);
+
+}
+
+
+
+
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/MultipleScatteringSamplerGeneralMixture.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/MultipleScatteringSamplerGeneralMixture.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b53a4ed8053ce488ffc498eabd3f7b776478fa8b
--- /dev/null
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/MultipleScatteringSamplerGeneralMixture.cxx
@@ -0,0 +1,224 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// MultipleScatteringSamplerGeneralMixture.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+// Trk include
+#include "TrkGeometry/MaterialProperties.h"
+
+#include "ISF_FatrasTools/MultipleScatteringSamplerGeneralMixture.h"
+#include "AtlasCLHEP_RandomGenerators/RandGaussZiggurat.h"
+
+// static particle masses
+Trk::ParticleMasses iFatras::MultipleScatteringSamplerGeneralMixture::s_particleMasses;
+// static doubles
+
+double iFatras::MultipleScatteringSamplerGeneralMixture::s_main_RossiGreisen = 17.5*Gaudi::Units::MeV;
+double iFatras::MultipleScatteringSamplerGeneralMixture::s_log_RossiGreisen = 0.125;
+
+double iFatras::MultipleScatteringSamplerGeneralMixture::s_projectionFactor = sqrt(2.);
+
+// ============================= General mixture model =============
+double iFatras::MultipleScatteringSamplerGeneralMixture::s_genMixScale = 0.608236; //numerically derived scaling factor
+
+// constructor
+iFatras::MultipleScatteringSamplerGeneralMixture::MultipleScatteringSamplerGeneralMixture(const std::string& t, const std::string& n, const IInterface* p) :
+ AthAlgTool(t,n,p),
+ m_log_include(true),
+ m_rndGenSvc("AtRndmGenSvc", n),
+ m_randomEngine(0),
+ m_randomEngineName("TrkExRnd")
+{
+ declareInterface<IMultipleScatteringSampler>(this);
+ // multiple scattering parameters
+ declareProperty("MultipleScatteringLogarithmicTermOn", m_log_include);
+ // random service for Gaussian mixture model
+ declareProperty("RandomNumberService", m_rndGenSvc, "Name of the random number service");
+ declareProperty("RandomStreamName", m_randomEngineName, "Name of the random number stream");
+}
+
+// destructor
+iFatras::MultipleScatteringSamplerGeneralMixture::~MultipleScatteringSamplerGeneralMixture()
+{}
+
+// Athena standard methods
+// initialize
+StatusCode iFatras::MultipleScatteringSamplerGeneralMixture::initialize()
+{
+ // get the random generator service
+ if (m_rndGenSvc.retrieve().isFailure()) {
+ ATH_MSG_WARNING( "Could not retrieve " << m_rndGenSvc );
+ } else {
+ ATH_MSG_VERBOSE( "Successfully retrieved " << m_rndGenSvc );
+ }
+
+ //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;
+ } else
+ ATH_MSG_VERBOSE ("Successfully got random engine '" << m_randomEngineName << "'");
+
+
+ ATH_MSG_INFO( "initialize() successful" );
+ return StatusCode::SUCCESS;
+}
+
+// finalize
+StatusCode iFatras::MultipleScatteringSamplerGeneralMixture::finalize()
+{
+ ATH_MSG_INFO( "finalize() successful" );
+ return StatusCode::SUCCESS;
+}
+
+
+double iFatras::MultipleScatteringSamplerGeneralMixture::simTheta(const Trk::MaterialProperties& mat,
+ double p,
+ double pathcorrection,
+ Trk::ParticleHypothesis particle,
+ double ) const
+{
+ if (mat.thicknessInX0()<=0. || particle==Trk::geantino) return 0.;
+
+ // make sure the path correction is positive to avoid a floating point exception
+ pathcorrection *= pathcorrection < 0. ? (-1.) : (1) ;
+
+ // scale the path length to the radiation length
+ double dOverX0 = pathcorrection * mat.thicknessInX0();
+
+ // kinematics (relativistic)
+ double m = s_particleMasses.mass[particle];
+ double E = sqrt(p*p + m*m);
+ double beta = p/E;
+
+ //material properties
+ double Z = mat.averageZ(); //charge layer material
+
+ double sigma2(0.);
+ double theta(0.);
+
+
+ if (particle != Trk::electron) {
+ //----------------------------------------------------------------------------------------------//
+ //see Mixture models of multiple scattering: computation and simulation. - R.Frühwirth, M. Liendl. -
+ //Computer Physics Communications 141 (2001) 230–246
+ //----------------------------------------------------------------------------------------------//
+ double * scattering_params;
+ // Decide which mixture is best
+ if (dOverX0/(beta*beta)>0.6/pow(Z,0.6)){ //Gaussian
+ // Gaussian mixture or pure Gaussian
+ if (dOverX0/(beta*beta)>10){
+ scattering_params=getGaussian(beta,p,dOverX0,s_genMixScale); // Get parameters
+ //std::cout<<"MultipleScatteringSamplerGeneralMixture::multipleScatteringUpdate: using pure_gaussian"<<std::endl;
+ }
+ else{
+ scattering_params=getGaussmix(beta,p,dOverX0,Z,s_genMixScale); // Get parameters
+ //std::cout<<"MultipleScatteringSamplerGeneralMixture::multipleScatteringUpdate: using gaussian_mixture"<<std::endl;
+ }
+ theta = simGaussmix(scattering_params); // Simulate
+ }
+ else{
+ //Semigaussian mixture
+ scattering_params = getSemigauss(beta,p,dOverX0,Z,s_genMixScale); // Get parameters
+ //std::cout<<"MultipleScatteringSamplerGeneralMixture::multipleScatteringUpdate: using semi_gaussian mixture"<<std::endl;
+ theta = simSemigauss(scattering_params); // Simulate
+ }
+ }
+
+ else {
+ // Electron multiple scattering effects - see Highland NIM 129 (1975) p497-499
+ // (Highland extension to the Rossi-Greisen formulation)
+ // NOTE: The formula can be extended by replacing the momentum^2 term with pi * pf
+ sigma2 = s_main_RossiGreisen / ( beta * p );
+ sigma2 *= (sigma2*dOverX0);
+
+ if ( m_log_include ) {
+ double factor = 1. + s_log_RossiGreisen * log10( 10. * dOverX0 );
+ factor *= factor;
+ sigma2 *= factor;
+ }
+
+ theta=sqrt(sigma2)*CLHEP::RandGaussZiggurat::shoot(m_randomEngine);
+ }
+
+ return theta*s_projectionFactor;
+}
+
+double * iFatras::MultipleScatteringSamplerGeneralMixture::getGaussian(double beta, double p,double dOverX0, double scale) const{
+ double * scattering_params = new double[4];
+ scattering_params[0]=15./beta/p*sqrt(dOverX0)*scale; //Total standard deviation of mixture
+ scattering_params[1]=1.0; //Variance of core
+ scattering_params[2]=1.0; //Variance of tails
+ scattering_params[3]=0.5; //Mixture weight of tail component
+ return scattering_params;
+}
+
+double * iFatras::MultipleScatteringSamplerGeneralMixture::getGaussmix(double beta, double p,double dOverX0,double Z, double scale) const{
+ double * scattering_params = new double[4];
+ scattering_params[0]=15./beta/p*sqrt(dOverX0)*scale; //Total standard deviation of mixture
+ double d1=log(dOverX0/(beta*beta));
+ double d2=log(pow(Z,2.0/3.0)*dOverX0/(beta*beta));
+ double epsi;
+ double var1=(-1.843e-3*d1+3.347e-2)*d1+8.471e-1; //Variance of core
+ if (d2<0.5)
+ epsi=(6.096e-4*d2+6.348e-3)*d2+4.841e-2;
+ else
+ epsi=(-5.729e-3*d2+1.106e-1)*d2-1.908e-2;
+ scattering_params[1]=var1; //Variance of core
+ scattering_params[2]=(1-(1-epsi)*var1)/epsi; //Variance of tails
+ scattering_params[3]=epsi; //Mixture weight of tail component
+ return scattering_params;
+}
+
+double * iFatras::MultipleScatteringSamplerGeneralMixture::getSemigauss(double beta,double p,double dOverX0,double Z, double scale) const{
+ double * scattering_params = new double[6];
+ double N=dOverX0*1.587E7*pow(Z,1.0/3.0)/(beta*beta)/(Z+1)/log(287/sqrt(Z));
+ scattering_params[4]=15./beta/p*sqrt(dOverX0)*scale; //Total standard deviation of mixture
+ double rho=41000/pow(Z,2.0/3.0);
+ double b=rho/sqrt(N*(log(rho)-0.5));
+ double n=pow(Z,0.1)*log(N);
+ double var1=(5.783E-4*n+3.803E-2)*n+1.827E-1;
+ double a=(((-4.590E-5*n+1.330E-3)*n-1.355E-2)*n+9.828E-2)*n+2.822E-1;
+ double epsi = (1-var1)/(a*a*(log(b/a)-0.5)-var1);
+ scattering_params[3]=(epsi>0) ? epsi : 0.0;//Mixture weight of tail component
+ scattering_params[0]=a; //Parameter 1 of tails
+ scattering_params[1]=b; //Parameter 2 of tails
+ scattering_params[2]=var1; //Variance of core
+ scattering_params[5]=N; //Average number of scattering processes
+ return scattering_params;
+}
+
+double iFatras::MultipleScatteringSamplerGeneralMixture::simGaussmix(double * scattering_params) const{
+ double sigma_tot = scattering_params[0];
+ double var1 = scattering_params[1];
+ double var2 = scattering_params[2];
+ double epsi = scattering_params[3];
+ bool ind = CLHEP::RandFlat::shoot(m_randomEngine)>epsi;
+ double u=CLHEP::RandFlat::shoot(m_randomEngine);
+ if(ind)
+ return sqrt(var1)*sqrt(-2*log(u))*sigma_tot;
+ else
+ return sqrt(var2)*sqrt(-2*log(u))*sigma_tot;
+}
+
+double iFatras::MultipleScatteringSamplerGeneralMixture::simSemigauss(double * scattering_params) const{
+ double a = scattering_params[0];
+ double b = scattering_params[1];
+ double var1 = scattering_params[2];
+ double epsi = scattering_params[3];
+ double sigma_tot = scattering_params[4];
+ bool ind=CLHEP::RandFlat::shoot(m_randomEngine)>epsi;
+ double u=CLHEP::RandFlat::shoot(m_randomEngine);
+ if(ind)
+ return sqrt(var1)*sqrt(-2*log(u))*sigma_tot;
+ else
+ return a*b*sqrt((1-u)/(u*b*b+a*a))*sigma_tot;
+}
+
+
+
+
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/MultipleScatteringSamplerHighland.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/MultipleScatteringSamplerHighland.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..163fa899838871503b8624170f0ea105303c3276
--- /dev/null
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/MultipleScatteringSamplerHighland.cxx
@@ -0,0 +1,139 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// MultipleScatteringSamplerHighland.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+// Trk include
+#include "TrkGeometry/MaterialProperties.h"
+#include "TrkExUtils/MaterialInteraction.h"
+
+#include "ISF_FatrasTools/MultipleScatteringSamplerHighland.h"
+
+#include "AtlasCLHEP_RandomGenerators/RandGaussZiggurat.h"
+
+// static particle masses
+Trk::ParticleMasses iFatras::MultipleScatteringSamplerHighland::s_particleMasses;
+// static doubles
+double iFatras::MultipleScatteringSamplerHighland::s_main_RutherfordScott = 13.6*Gaudi::Units::MeV;
+double iFatras::MultipleScatteringSamplerHighland::s_log_RutherfordScott = 0.038;
+
+double iFatras::MultipleScatteringSamplerHighland::s_main_RossiGreisen = 17.5*Gaudi::Units::MeV;
+double iFatras::MultipleScatteringSamplerHighland::s_log_RossiGreisen = 0.125;
+
+double iFatras::MultipleScatteringSamplerHighland::s_projectionFactor = sqrt(2.);
+
+Trk::MaterialInteraction m_matInt;
+
+// constructor
+iFatras::MultipleScatteringSamplerHighland::MultipleScatteringSamplerHighland(const std::string& t, const std::string& n, const IInterface* p) :
+ AthAlgTool(t,n,p),
+ m_log_include(true),
+ m_rndGenSvc("AtRndmGenSvc", n),
+ m_randomEngine(0),
+ m_randomEngineName("TrkExRnd")
+{
+ declareInterface<IMultipleScatteringSampler>(this);
+ // multiple scattering parameters
+ declareProperty("MultipleScatteringLogarithmicTermOn", m_log_include);
+ // random service
+ declareProperty("RandomNumberService", m_rndGenSvc, "Name of the random number service");
+ declareProperty("RandomStreamName", m_randomEngineName, "Name of the random number stream");
+}
+
+// destructor
+iFatras::MultipleScatteringSamplerHighland::~MultipleScatteringSamplerHighland()
+{}
+
+// Athena standard methods
+// initialize
+StatusCode iFatras::MultipleScatteringSamplerHighland::initialize()
+{
+
+ // get the random generator service
+ if (m_rndGenSvc.retrieve().isFailure()) {
+ ATH_MSG_WARNING( "Could not retrieve " << m_rndGenSvc );
+ } else {
+ ATH_MSG_VERBOSE( "Successfully retrieved " << m_rndGenSvc );
+ }
+
+ //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;
+ } else
+ ATH_MSG_VERBOSE ("Successfully got random engine '" << m_randomEngineName << "'");
+
+ ATH_MSG_INFO( "initialize() successful" );
+ return StatusCode::SUCCESS;
+}
+
+// finalize
+StatusCode iFatras::MultipleScatteringSamplerHighland::finalize()
+{
+ ATH_MSG_INFO( "finalize() successful" );
+ return StatusCode::SUCCESS;
+}
+
+
+double iFatras::MultipleScatteringSamplerHighland::simTheta(const Trk::MaterialProperties& mat,
+ double p,
+ double pathcorrection,
+ Trk::ParticleHypothesis particle,
+ double ) const
+{
+ if (mat.thicknessInX0()<=0. || particle==Trk::geantino) return 0.;
+
+ // make sure the path correction is positive to avoid a floating point exception
+ pathcorrection *= pathcorrection < 0. ? (-1.) : (1) ;
+
+ // scale the path length to the radiation length
+ double t = pathcorrection * mat.thicknessInX0();
+
+ // kinematics (relativistic)
+ double m = s_particleMasses.mass[particle];
+ double E = sqrt(p*p + m*m);
+ double beta = p/E;
+
+ double sigma2(0.);
+
+ double sigma = m_matInt.sigmaMS(t, p, beta);
+ sigma2 = sigma*sigma;
+
+ // Code below will not be used if the parameterization of TrkUtils is used
+ if (particle != Trk::electron) {
+
+ // the highland formula
+ sigma2 = s_main_RutherfordScott/(beta*p);
+
+ if (m_log_include)
+ sigma2 *= (1.+s_log_RutherfordScott*log(t));
+
+ sigma2 *= (sigma2*t);
+ }
+
+ else {
+
+ // Electron multiple scattering effects - see Highland NIM 129 (1975) p497-499
+ // (Highland extension to the Rossi-Greisen formulation)
+ // NOTE: The formula can be extended by replacing the momentum^2 term with pi * pf
+ sigma2 = s_main_RossiGreisen / ( beta * p );
+ sigma2 *= (sigma2*t);
+
+ if ( m_log_include ) {
+ double factor = 1. + s_log_RossiGreisen * log10( 10. * t );
+ factor *= factor;
+ sigma2 *= factor;
+ }
+ }
+
+ return s_projectionFactor*sqrt(sigma2)*CLHEP::RandGaussZiggurat::shoot(m_randomEngine);
+
+}
+
+
+
+
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/PhotonConversionTool.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/PhotonConversionTool.cxx
index cc5ad052f737d6bf05b306bfac04f74fa5aabd7b..ccd4af9bf42a64f3a741caf5381893055c22daab 100755
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/PhotonConversionTool.cxx
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/PhotonConversionTool.cxx
@@ -29,11 +29,11 @@
#include "TrkExInterfaces/IEnergyLossUpdator.h"
#include "TrkExInterfaces/ITimedExtrapolator.h"
#include "TrkExInterfaces/IMultipleScatteringUpdator.h"
-#include "TrkParameters/TrackParameters.h"
#include "TrkSurfaces/Surface.h"
#include "TrkGeometry/Layer.h"
#include "TrkGeometry/MaterialProperties.h"
#include "TrkVolumes/CylinderVolumeBounds.h"
+#include "TrkNeutralParameters/NeutralParameters.h"
// CLHEP
#include "CLHEP/Units/SystemOfUnits.h"
#include "CLHEP/Matrix/Vector.h"
@@ -252,10 +252,10 @@ void iFatras::PhotonConversionTool::recordChilds(double time,
*parent );
// in the validation mode, add process info
if (m_validationMode) {
- ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
+ ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
validInfo->setProcess(14);
if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
- else validInfo->setGeneration(1); // assume parent is a primary track
+ else validInfo->setGeneration(1); // assume parent is a primary track
ch1->setUserInformation(validInfo);
}
children[ichild] = ch1;
@@ -274,10 +274,10 @@ void iFatras::PhotonConversionTool::recordChilds(double time,
// in the validation mode, add process info
if (m_validationMode) {
- ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
+ ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
validInfo->setProcess(14);
- if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
- else validInfo->setGeneration(1); // assume parent is a primary track
+ if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
+ else validInfo->setGeneration(1); // assume parent is a primary track
ch2->setUserInformation(validInfo);
}
children[ichild] = ch2;
@@ -348,13 +348,13 @@ ISF::ISFParticleVector iFatras::PhotonConversionTool::getChilds(const ISF::ISFPa
pdg1,
time,
*parent );
- if (m_validationMode) {
- ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
- validInfo->setProcess(14);
- if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
- else validInfo->setGeneration(1); // assume parent is a primary track
- ch1->setUserInformation(validInfo);
- }
+ //if (m_validationMode) {
+ // ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
+ // validInfo->setProcess(14);
+ // if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
+ // else validInfo->setGeneration(1); // assume parent is a primary track
+ // ch1->setUserInformation(validInfo);
+ //}
children[0] = ch1;
ISF::ISFParticle* ch2 = new ISF::ISFParticle( vertex,
@@ -365,28 +365,28 @@ ISF::ISFParticleVector iFatras::PhotonConversionTool::getChilds(const ISF::ISFPa
time,
*parent );
- if (m_validationMode) {
- ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
- validInfo->setProcess(14);
- if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
- else validInfo->setGeneration(1); // assume parent is a primary track
- ch2->setUserInformation(validInfo);
- }
+ //if (m_validationMode) {
+ // ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
+ // validInfo->setProcess(14);
+ // if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
+ // else validInfo->setGeneration(1); // assume parent is a primary track
+ // ch2->setUserInformation(validInfo);
+ //}
children[1] = ch2;
// register TruthIncident
- ISF::ISFTruthIncident truth( const_cast<ISF::ISFParticle&>(*parent),
- children,
- m_processCode,
- parent->nextGeoID(),
- ISF::fKillsPrimary );
- m_truthRecordSvc->registerTruthIncident( truth);
+ //ISF::ISFTruthIncident truth( const_cast<ISF::ISFParticle&>(*parent),
+ // children,
+ // m_processCode,
+ // parent->nextGeoID(),
+ // ISF::fKillsPrimary );
+ //m_truthRecordSvc->registerTruthIncident( truth);
// save info for validation
- if (m_validationMode && m_validationTool) {
- Amg::Vector3D* nPrim=0;
- m_validationTool->saveISFVertexInfo(14,vertex,*parent,parent->momentum(),nPrim,children);
- }
+ //if (m_validationMode && m_validationTool) {
+ // Amg::Vector3D* nPrim=0;
+ // m_validationTool->saveISFVertexInfo(14,vertex,*parent,parent->momentum(),nPrim,children);
+ //}
return children;
}
@@ -599,7 +599,7 @@ Amg::Vector3D iFatras::PhotonConversionTool::childDirection(const Amg::Vector3D&
}
-bool iFatras::PhotonConversionTool::doConversion(double time, const Trk::TrackParameters& parm,
+bool iFatras::PhotonConversionTool::doConversion(double time, const Trk::NeutralParameters& parm,
const Trk::ExtendedMaterialProperties* /*extMatProp*/) const {
double p = parm.momentum().mag();
@@ -627,7 +627,7 @@ bool iFatras::PhotonConversionTool::doConversion(double time, const Trk::TrackPa
/** interface for processing of the presampled nuclear interactions on layer*/
ISF::ISFParticleVector iFatras::PhotonConversionTool::doConversionOnLayer(const ISF::ISFParticle* parent,
- double time, const Trk::TrackParameters& parm,
+ double time, const Trk::NeutralParameters& parm,
const Trk::ExtendedMaterialProperties* /*ematprop*/) const {
double p = parm.momentum().mag();
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/PhysicsValidationTool.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/PhysicsValidationTool.cxx
index d253484d56785ff407d46b13c81ac3f97bb1701d..14db928181f1e399cb43680cffce6818365793d9 100755
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/PhysicsValidationTool.cxx
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/PhysicsValidationTool.cxx
@@ -76,6 +76,10 @@ iFatras::PhysicsValidationTool::initialize()
m_particles->Branch("pph" , &m_pph , "pph/F" );
m_particles->Branch("p" , &m_p , "p/F" );
m_particles->Branch("eloss" , &m_eloss , "eloss/F" );
+ m_particles->Branch("ionloss" , &m_ionloss , "ionloss/F" );
+ m_particles->Branch("radloss" , &m_radloss , "radloss/F" );
+ m_particles->Branch("zOaTr" , &m_zOaTr , "zOaTr/F" );
+ m_particles->Branch("wZ" , &m_wZ , "wZ/F" );
m_particles->Branch("X0" , &m_X0 , "X0/F" );
m_particles->Branch("dt" , &m_dt , "dt/F" );
m_particles->Branch("thIn" , &m_thIn , "thIn/F" );
@@ -130,6 +134,83 @@ StatusCode iFatras::PhysicsValidationTool::finalize()
* ==> see headerfile
*=======================================================================*/
+/** new transport tool */
+void iFatras::PhysicsValidationTool::saveISFParticleInfo(const ISF::ISFParticle& isp,
+ const Trk::ExtrapolationCell<Trk::TrackParameters>& ec,
+ Trk::ExtrapolationCode ecode) {
+
+ m_p = isp.momentum().mag();
+ m_ionloss = ec.eLoss ? ec.eLoss->meanIoni() : 0.;
+ m_radloss = ec.eLoss ? ec.eLoss->meanRad() : 0.; // average estimate
+ m_zOaTr = ec.zOaTrX;
+ m_wZ = ec.zX;
+
+ if (ecode == Trk::ExtrapolationCode::SuccessBoundaryReached) { // surviving particle
+ m_scEnd = 0;
+ m_eloss = ec.endParameters->momentum().mag()-m_p;
+ //m_dt = time - isp.timeStamp();
+ m_dt = ec.pathLength/CLHEP::c_light; // TODO get timing info from the cell
+
+ m_thEnd = ec.endParameters->position().theta();
+ m_phEnd = ec.endParameters->position().phi();
+ m_dEnd = ec.endParameters->position().mag();
+
+ } else {
+
+ m_scEnd = -1; // undefined
+ if (ecode == Trk::ExtrapolationCode::SuccessPathLimit) m_scEnd = 201; // decay
+ else if (ecode == Trk::ExtrapolationCode::SuccessMaterialLimit ) m_scEnd = ec.materialProcess; // material interaction
+
+ m_eloss = -m_p;
+ m_thEnd = -10.; // TODO : retrieve position of stopped particle
+ m_phEnd = -10.; // TODO : retrieve position of stopped particle
+ m_dEnd = -1.; // TODO : retrieve position of stopped particle
+ }
+ m_X0 = ec.materialX0;
+
+ saveInfo(isp);
+
+}
+
+/** new transport tool */
+void iFatras::PhysicsValidationTool::saveISFParticleInfo(const ISF::ISFParticle& isp,
+ const Trk::ExtrapolationCell<Trk::NeutralParameters>& ec,
+ Trk::ExtrapolationCode ecode) {
+
+ m_p = isp.momentum().mag();
+ m_ionloss = 0.;
+ m_radloss = 0.;
+ m_zOaTr = ec.zOaTrX;
+ m_wZ = ec.zX;
+
+ if (ecode == Trk::ExtrapolationCode::SuccessBoundaryReached) { // surviving particle
+ m_scEnd = 0;
+ m_eloss = ec.endParameters->momentum().mag()-m_p;
+ //m_dt = time - isp.timeStamp();
+ m_dt = ec.pathLength/CLHEP::c_light; // TODO get timing info from the cell
+
+ m_thEnd = ec.endParameters->position().theta();
+ m_phEnd = ec.endParameters->position().phi();
+ m_dEnd = ec.endParameters->position().mag();
+
+ } else {
+
+ m_scEnd = -1; // undefined
+ if (ecode == Trk::ExtrapolationCode::SuccessPathLimit) m_scEnd = 201; // decay
+ else if (ecode == Trk::ExtrapolationCode::SuccessMaterialLimit ) m_scEnd = ec.materialProcess; // material interaction
+
+ m_eloss = -m_p;
+ m_thEnd = -10.; // TODO : retrieve position of stopped particle
+ m_phEnd = -10.; // TODO : retrieve position of stopped particle
+ m_dEnd = -1.; // TODO : retrieve position of stopped particle
+ }
+ m_X0 = ec.materialX0;
+
+ saveInfo(isp);
+
+}
+
+/** old transport tool/prompt decay */
void iFatras::PhysicsValidationTool::saveISFParticleInfo(const ISF::ISFParticle& isp,
int endProcess,
const Trk::TrackParameters* ePar,
@@ -137,17 +218,22 @@ void iFatras::PhysicsValidationTool::saveISFParticleInfo(const ISF::ISFParticle&
m_pdg = isp.pdgCode();
m_scIn = isp.getUserInformation()? isp.getUserInformation()->process() : 0; // assume primary track
- m_scEnd = endProcess;
m_gen = isp.getUserInformation()? isp.getUserInformation()->generation() : 0; // assume primary track
m_pth=isp.momentum().theta();
m_pph=isp.momentum().phi();
- m_p =isp.momentum().mag();
m_geoID = isp.nextGeoID();
- m_eloss = ePar? ePar->momentum().mag()-m_p : -m_p;
- m_dt = time - isp.timeStamp();
m_dIn = isp.position().mag();
m_thIn = isp.position().theta();
m_phIn = isp.position().phi();
+ m_p = isp.momentum().mag();
+ m_scEnd = endProcess;
+ m_eloss = ePar? ePar->momentum().mag()-m_p : -m_p;
+ m_dt = time - isp.timeStamp();
+
+ m_ionloss = 0.; // n.a.
+ m_radloss = 0.; // n.a.
+ m_zOaTr = 0.; // n.a.
+
if (ePar) {
//error propagation
@@ -163,6 +249,25 @@ void iFatras::PhysicsValidationTool::saveISFParticleInfo(const ISF::ISFParticle&
m_dEnd = -1.; // particle stopped and position unknown
}
m_X0 = dX0;
+
+ m_particles->Fill();
+
+}
+
+void iFatras::PhysicsValidationTool::saveInfo(const ISF::ISFParticle& isp) {
+
+ // ISF particle info ( common )
+
+ m_pdg = isp.pdgCode();
+ m_scIn = isp.getUserInformation()? isp.getUserInformation()->process() : 0; // assume primary track
+ m_gen = isp.getUserInformation()? isp.getUserInformation()->generation() : 0; // assume primary track
+ m_pth=isp.momentum().theta();
+ m_pph=isp.momentum().phi();
+ m_geoID = isp.nextGeoID();
+ m_dIn = isp.position().mag();
+ m_thIn = isp.position().theta();
+ m_phIn = isp.position().phi();
+
m_particles->Fill();
}
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/ProcessSamplingTool.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/ProcessSamplingTool.cxx
index b3fceeab2e09b743b90821c29ebd22713e034741..8830505feb3e79687607aa33029de4c4070e16a4 100755
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/ProcessSamplingTool.cxx
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/ProcessSamplingTool.cxx
@@ -13,6 +13,14 @@
#include "TrkParameters/TrackParameters.h"
// ISF includes
#include "ISF_Event/ISFParticle.h"
+#include "ISF_Event/ITruthIncident.h"
+#include "ISF_Event/ISFTruthIncident.h"
+#include "ISF_Interfaces/ITruthSvc.h"
+// iFatras
+#include "ISF_FatrasInterfaces/IHadronicInteractionProcessor.h"
+#include "ISF_FatrasInterfaces/IPhysicsValidationTool.h"
+#include "ISF_FatrasInterfaces/IPhotonConversionTool.h"
+#include "ISF_FatrasInterfaces/IParticleDecayHelper.h"
// CLHEP
#include "CLHEP/Units/SystemOfUnits.h"
#include "CLHEP/Random/RandFlat.h"
@@ -28,13 +36,29 @@ iFatras::ProcessSamplingTool::ProcessSamplingTool( const std::string& t,
const std::string& n,
const IInterface* p )
: AthAlgTool(t,n,p),
- m_rndGenSvc("AtDSFMTGenSvc", n),
- m_randomEngine(0),
- m_randomEngineName("FatrasRnd")
+ m_rndGenSvc("AtDSFMTGenSvc", n),
+ m_randomEngine(0),
+ m_randomEngineName("FatrasRnd"),
+ m_truthRecordSvc("ISF_ValidationTruthService", n),
+ m_hadInt(true),
+ m_hadIntProcessor("iFatras::G4HadIntProcessor/FatrasG4HadIntProcessor"),
+ m_particleDecayer(),
+ m_conversionTool("iFatras::McConversionCreator/FatrasConversionCreator"),
+ m_validationMode(false),
+ m_validationTool("")
{
declareInterface<IProcessSamplingTool>(this);
// service handles
- declareProperty( "RandomNumberService", m_rndGenSvc );
+ declareProperty( "RandomNumberService" , m_rndGenSvc );
+ // hadronic interactions
+ declareProperty( "HadronicInteraction" , m_hadInt);
+ declareProperty( "HadronicInteractionProcessor" , m_hadIntProcessor);
+ declareProperty( "ParticleDecayHelper" , m_particleDecayer );
+ declareProperty( "PhotonConversionTool" , m_conversionTool);
+ // validation mode ------------------------------------------------------------
+ declareProperty("ValidationMode" , m_validationMode);
+ declareProperty("PhysicsValidationTool" , m_validationTool);
+ declareProperty("TruthRecordSvc" , m_truthRecordSvc , "ISF Particle Truth Svc");
}
/*=========================================================================
@@ -65,6 +89,43 @@ iFatras::ProcessSamplingTool::initialize()
return StatusCode::FAILURE;
}
+ if (m_truthRecordSvc.retrieve().isFailure()){
+ ATH_MSG_FATAL( "Could not retrieve " << m_truthRecordSvc );
+ return StatusCode::FAILURE;
+ }
+
+ // retrieve the photon conversion tool
+ if (m_conversionTool.retrieve().isFailure()){
+ ATH_MSG_FATAL( "Could not retrieve " << m_conversionTool );
+ return StatusCode::FAILURE;
+ } else
+ ATH_MSG_VERBOSE( "Successfully retrieved " << m_conversionTool );
+
+ // retrieve the hadronic interaction tool
+ if (m_hadInt) {
+ if (m_hadIntProcessor.retrieve().isFailure()){
+ ATH_MSG_FATAL( "Could not retrieve " << m_hadIntProcessor );
+ return StatusCode::FAILURE;
+ } else
+ ATH_MSG_VERBOSE( "Successfully retrieved " << m_hadIntProcessor );
+ }
+
+ // Particle decayer
+ if (m_particleDecayer.retrieve().isFailure()){
+ ATH_MSG_FATAL( "Could not retrieve " << m_particleDecayer );
+ return StatusCode::FAILURE;
+ } else
+ ATH_MSG_VERBOSE( "Successfully retrieved " << m_particleDecayer );
+
+ // retrieve the physics validation tool
+ if (m_validationMode) {
+ if (m_validationTool.retrieve().isFailure()){
+ ATH_MSG_FATAL( "Could not retrieve " << m_validationTool );
+ return StatusCode::FAILURE;
+ } else
+ ATH_MSG_VERBOSE( "Successfully retrieved " << m_validationTool );
+ }
+
return StatusCode::SUCCESS;
}
@@ -78,7 +139,7 @@ StatusCode iFatras::ProcessSamplingTool::finalize()
return StatusCode::SUCCESS;
}
-Trk::PathLimit iFatras::ProcessSamplingTool::sampleProcess(double momentum,double charge,Trk::ParticleHypothesis particle)
+Trk::PathLimit iFatras::ProcessSamplingTool::sampleProcess(double momentum,double charge,Trk::ParticleHypothesis particle) const
{
int process=0;
double x0Max = -1.;
@@ -144,28 +205,189 @@ Trk::PathLimit iFatras::ProcessSamplingTool::sampleProcess(double momentum,doubl
}
// presumably here we are left with hadrons only
+ if (m_hadInt) {
+
+ // sample free path in terms of nuclear interaction length
+ double al = 1.; // scaling here
+
+ /*
+
+ if ( particle == Trk::pion || particle == Trk::kaon || particle == Trk::pi0 || particle == Trk::k0) {
+ al *= 1./(1.+ exp(-0.5*(momentum-270.)*(momentum-270.)/60./60.));
+ }
+ if ( particle == Trk::proton || particle == Trk::neutron ) al *=0.7;
+ if ( particle == Trk::pion || particle == Trk::pi0) al *=0.9;
+ */
+
+ x0Max = -log(rndx)*al ;
+
+ process = 121;
+
+ //std::cout <<"hadronic path limit:"<<momentum<<","<<al<<","<< x0Max << std::endl;
+ return Trk::PathLimit(x0Max,process);
+ }
+
+
+ return Trk::PathLimit(x0Max,process);
+
+ // return Trk::PathLimit(x0Max,process);
+}
+
+
+ISF::ISFParticleVector iFatras::ProcessSamplingTool::interact(const ISF::ISFParticle* parent,
+ Trk::ExCellCharged& eCell,
+ const Trk::Material* mat) const
+{
+ ISF::ISFParticleVector childVector(0);
+
+ int process = eCell.materialProcess;
+ if ( process==0 ) return childVector;
+
+ Trk::ParticleHypothesis particle = eCell.pHypothesis;
+ const Trk::TrackParameters* parm=eCell.leadParameters;
+
+ double mass = m_particleMasses.mass[particle];
+ double p = parm->momentum().mag();
+ //double E = sqrt( p*p + mass*mass);
+ Amg::Vector3D position=parm->position();
+ Amg::Vector3D momentum=parm->momentum();
+
+ if ( process == 5 ) { // positron annihilation
+
+ double fmin = mass/p;
+
+ double fr = 1.-pow(fmin,CLHEP::RandFlat::shoot(m_randomEngine));
+
+ // double cTh = 1.-fmin/(1.-fr)/fr;
- // sample free path in terms of nuclear interaction length
- double al = 1.; // scaling here
+ // first implementation: ctH=1
- /*
+ ISF::ISFParticle* child1 = new ISF::ISFParticle( position,
+ fr*momentum,
+ 0.,
+ 0.,
+ 22,
+ eCell.time,
+ *parent );
- if ( particle == Trk::pion || particle == Trk::kaon || particle == Trk::pi0 || particle == Trk::k0) {
- al *= 1./(1.+ exp(-0.5*(momentum-270.)*(momentum-270.)/60./60.));
+ ISF::ISFParticle* child2 = new ISF::ISFParticle( position,
+ (1-fr)*momentum,
+ 0.,
+ 0.,
+ 22,
+ eCell.time,
+ *parent );
+
+ childVector.push_back(child1);
+ childVector.push_back(child2);
+
+ } else if (process==121) { // hadronic interaction
+
+ //const Amg::Vector3D pDir = momentum.unit();
+
+ // const Trk::CurvilinearParameters parm(position,momentum,parent->charge());
+
+ childVector = m_hadIntProcessor->doHadIntOnLayer(parent, eCell.time, position, momentum, mat, particle);
+
+ } else if ( process == 201 ) { // decay
+ childVector = m_particleDecayer->decayParticle(*parent, position, momentum, eCell.time);
+ }
+
+ // validation mode
+ if (m_validationMode) {
+
+ // add process info for children
+ for (unsigned int i=0; i<childVector.size(); i++) {
+ ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
+ validInfo->setProcess(process);
+ if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
+ else validInfo->setGeneration(1); // assume parent is a primary track
+ childVector[i]->setUserInformation(validInfo);
+ }
+ // save interaction info
+ if ( m_validationTool ) {
+ Amg::Vector3D* nMom = 0;
+ m_validationTool->saveISFVertexInfo(process, position,*parent,momentum,nMom,childVector);
+ delete nMom;
+ }
}
- if ( particle == Trk::proton || particle == Trk::neutron ) al *=0.7;
- if ( particle == Trk::pion || particle == Trk::pi0) al *=0.9;
- */
- x0Max = -log(rndx)*al ;
+ // register TruthIncident
+ ISF::ISFTruthIncident truth( const_cast<ISF::ISFParticle&>(*parent),
+ childVector,
+ process,
+ parent->nextGeoID(), // inherits from the parent
+ ISF::fKillsPrimary );
+ m_truthRecordSvc->registerTruthIncident( truth);
+
+ //
+ return childVector;
+}
+
+ISF::ISFParticleVector iFatras::ProcessSamplingTool::interact(const ISF::ISFParticle* parent,
+ Trk::ExCellNeutral& eCell,
+ const Trk::Material* mat) const
+{
+ ISF::ISFParticleVector childVector(0);
+
+ int process = eCell.materialProcess;
+ Trk::ParticleHypothesis particle = eCell.pHypothesis;
+ const Trk::NeutralParameters* parm=eCell.leadParameters;
+
+ if ( process==0 ) return childVector;
+ //double mass = m_particleMasses.mass[particle];
+ //double p = parm.momentum().mag();
+ /*double E = sqrt( p*p + mass*mass);*/
+ Amg::Vector3D position=parm->position();
+ Amg::Vector3D momentum=parm->momentum();
- process = 121;
+ if (process==14) { // photon conversion
- //std::cout <<"hadronic path limit:"<<momentum<<","<<al<<","<< x0Max << std::endl;
+ childVector = m_conversionTool->doConversionOnLayer(parent, eCell.time, *parm);
+
+ } else if (process==121) { // hadronic interaction
+
+ //const Amg::Vector3D pDir = momentum.unit();
+
+ //const Trk::CurvilinearParameters parm(position,momentum,parent->charge());
+
+ childVector = m_hadIntProcessor->doHadIntOnLayer(parent, eCell.time, position, momentum, mat, particle);
+
+ } else if ( process == 201 ) { // decay
+
+ childVector = m_particleDecayer->decayParticle(*parent, position, momentum, eCell.time);
+ }
+
+ // validation mode
+ if (m_validationMode) {
+
+ // add process info for children
+ for (unsigned int i=0; i<childVector.size(); i++) {
+ ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
+ validInfo->setProcess(process);
+ if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
+ else validInfo->setGeneration(1); // assume parent is a primary track
+ childVector[i]->setUserInformation(validInfo);
+ }
+ // save interaction info
+ if ( m_validationTool ) {
+ Amg::Vector3D* nMom = 0;
+ m_validationTool->saveISFVertexInfo(process, position,*parent,momentum,nMom,childVector);
+ delete nMom;
+ }
+ }
+
+ // register TruthIncident
+ ISF::ISFTruthIncident truth( const_cast<ISF::ISFParticle&>(*parent),
+ childVector,
+ process,
+ parent->nextGeoID(), // inherits from the parent
+ ISF::fKillsPrimary );
+ m_truthRecordSvc->registerTruthIncident( truth);
+
+ //
+ return childVector;
- return Trk::PathLimit(x0Max,process);
-
- // return Trk::PathLimit(x0Max,process);
}
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/TransportEngine.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/TransportEngine.cxx
index 1eeb264c0b025e2ea49f6c71f342a465f9fdcbbb..14ff7365b71b313f2237d4f0601f93bb2f5043be 100644
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/TransportEngine.cxx
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/TransportEngine.cxx
@@ -14,12 +14,16 @@
#include "TrkParameters/TrackParameters.h"
#include "TrkNeutralParameters/NeutralParameters.h"
#include "TrkExInterfaces/IExtrapolationEngine.h"
+#include "TrkExInterfaces/ITimedExtrapolator.h"
+
#include "TrkGeometry/Layer.h"
+#include "TrkGeometry/TrackingVolume.h"
// iFatras includes
#include "ISF_FatrasInterfaces/ISimHitCreator.h"
#include "ISF_FatrasInterfaces/IParticleDecayHelper.h"
#include "ISF_FatrasInterfaces/IProcessSamplingTool.h"
+#include "ISF_FatrasInterfaces/IPhysicsValidationTool.h"
// ISF includes
#include "ISF_Event/ISFParticle.h"
@@ -30,6 +34,9 @@
#include "CLHEP/Units/SystemOfUnits.h"
#include "CLHEP/Random/RandFlat.h"
+// AtlasDetDescr
+#include "AtlasDetDescr/AtlasRegion.h"
+
#include <iostream>
//#include "AtlasDetDescr/AtlasDetectorID.h"
@@ -50,10 +57,15 @@ iFatras::TransportEngine::TransportEngine( const std::string& t,
m_trackFilter(""),
m_neutralHadronFilter(""),
m_photonFilter(""),
- m_samplingTool("")
+ m_samplingTool(""),
+ m_validationMode(false),
+ m_validationTool("")
{
- declareInterface<ITransportTool>(this);
+ declareInterface<ISF::IParticleProcessor>(this);
+ // validation output section
+ declareProperty( "ValidationMode", m_validationMode );
+ declareProperty( "PhysicsValidationTool", m_validationTool );
// tool handle for the particle decayer
declareProperty( "ParticleDecayHelper", m_particleDecayHelper );
// tool handle for the track creator
@@ -64,7 +76,7 @@ iFatras::TransportEngine::TransportEngine( const std::string& t,
declareProperty( "TrackFilter", m_trackFilter );
declareProperty( "NeutralFilter", m_neutralHadronFilter );
declareProperty( "PhotonFilter", m_photonFilter );
- //
+ // tools handles ISF Framework
declareProperty( "ProcessSamplingTool", m_samplingTool );
// service handles
declareProperty( "RandomNumberService", m_rndGenSvc );
@@ -102,6 +114,8 @@ iFatras::TransportEngine::initialize()
return StatusCode::FAILURE;
if (retrieveTool<iFatras::IProcessSamplingTool>(m_samplingTool).isFailure())
return StatusCode::FAILURE;
+ if (m_validationMode && retrieveTool<iFatras::IPhysicsValidationTool>(m_validationTool).isFailure())
+ return StatusCode::FAILURE;
if ( m_rndGenSvc.retrieve().isFailure() ){
ATH_MSG_FATAL( "Could not retrieve " << m_rndGenSvc );
@@ -136,8 +150,8 @@ ISF::ISFParticle* iFatras::TransportEngine::process( const ISF::ISFParticle& isp
// copy the current particle onto the particle clipboard
ISF::ParticleClipboard::getInstance().setParticle(isp);
- ATH_MSG_INFO ("[ fatras transport ] processing particle " << isp.pdgCode() );
-
+ ATH_MSG_DEBUG ("[ fatras transport ] processing particle " << isp.pdgCode() );
+
// [0] ============================================================================================
// - INPUT PARTICLE PREPARATION
//
@@ -174,13 +188,14 @@ ISF::ISFParticle* iFatras::TransportEngine::process( const ISF::ISFParticle& isp
ATH_MSG_VERBOSE( "[ fatras transport ] Determined processor w/o filter.");
// now transport the particle - if it passes the filter
- if ( filter && !filter->passFilter(isp) ) {
- ATH_MSG_DEBUG( "[ fatras transport ] Filter not passed, ignore particle.");
- ATH_MSG_DEBUG( " -> for more information, turn on debug/verbose output of " << filter->name() );
- return 0;
- }
+ //if ( filter && !filter->passFilter(isp) ) {
+ // ATH_MSG_DEBUG( "[ fatras transport ] Filter not passed, ignore particle.");
+ // ATH_MSG_DEBUG( " -> for more information, turn on debug/verbose output of " << filter->name() );
+ // return 0;
+ //}
- ATH_MSG_INFO( "[ fatras transport ] The StackParticle passed filter - starting transport.");
+ //ATH_MSG_DEBUG( "[ fatras transport ] The StackParticle passed filter - starting transport.");
+ ATH_MSG_DEBUG( "[ fatras transport ] starting transport.");
// [a] DECAY section ------------------------------------------------------------------------------
//
@@ -192,7 +207,13 @@ ISF::ISFParticle* iFatras::TransportEngine::process( const ISF::ISFParticle& isp
if ( pathLimit > 0. && pathLimit < 0.01 ) {
if (!m_particleDecayHelper.empty()) {
ATH_MSG_VERBOSE( "[ fatras transport ] Decay is triggered for input particle.");
- m_particleDecayHelper->decay(isp);
+ m_particleDecayHelper->decay(isp,isp.position(),isp.momentum(),isp.timeStamp());
+ }
+
+ // validation mode - for all particle registered into stack
+ if ( m_validationMode && m_validationTool ) {
+ Trk::CurvilinearParameters inputPar(isp.position(),isp.momentum(),isp.charge());
+ m_validationTool->saveISFParticleInfo(isp,201,&inputPar,isp.timeStamp(),0.);
}
return 0;
}
@@ -202,20 +223,24 @@ ISF::ISFParticle* iFatras::TransportEngine::process( const ISF::ISFParticle& isp
//
double materialLimitX0 = -1.;
double materialLimitL0 = -1.;
- // now sample the path limit with the path sampling tool
- // -> the outcome is a material path limit
- Trk::PathLimit matLim(-1.,0);
- if (absPdg!=999 && pHypothesis<99) {
- // sample the material process
- matLim = m_samplingTool->sampleProcess(isp.momentum().mag(),isp.charge(),pHypothesis);
- // now assign the limit to L0 or X0 (121 is the code used inside the sampling tool)
- if (matLim.process == 121 ){
- materialLimitL0 = matLim.x0Max;
- } else
- materialLimitX0 = matLim.x0Max;
- ATH_MSG_VERBOSE( "[ fatras transport ] Particle free path : " << matLim.x0Max << " and process : " << matLim.process);
+ int materialProcess = -1;
+ double materialX0=0.;
+ double materialL0=0.;
+ // now sample the path limit if not done already
+ ISF::MaterialPathInfo* matLimit = isp.getUserInformation() ? isp.getUserInformation()->materialLimit() : 0;
+ if (matLimit) {
+ materialProcess = matLimit->process;
+ if (materialProcess==121) materialLimitL0 = matLimit->dMax;
+ else materialLimitX0 = matLimit->dMax;
+ if (materialProcess==121) materialL0=matLimit->dCollected;
+ else materialX0 = matLimit->dCollected;
+ } else if (absPdg!=999 && pHypothesis<99) { // need to resample
+ Trk::PathLimit matLim = m_samplingTool->sampleProcess(isp.momentum().mag(),isp.charge(),pHypothesis);
+ materialProcess = matLim.process;
+ if (materialProcess==121) materialLimitL0 = matLim.x0Max;
+ else materialLimitX0 = matLim.x0Max;
}
-
+
// [1] ============================================================================================
// - TRANSPORT THROUGH THE DETECTOR
//
@@ -237,9 +262,20 @@ ISF::ISFParticle* iFatras::TransportEngine::process( const ISF::ISFParticle& isp
pathLimit,
materialLimitX0,
materialLimitL0);
- ecc.materialProcess = matLim.process;
+ ecc.materialProcess = materialProcess;
+ if (ecc.materialLimitL0 != -1)
+ ecc.materialL0 = materialL0;
+ ecc.addConfigurationMode(Trk::ExtrapolationMode::StopWithMaterialLimitL0);
+ if (ecc.materialLimitX0 != -1)
+ ecc.materialX0 = materialX0;
+ ecc.addConfigurationMode(Trk::ExtrapolationMode::StopWithMaterialLimitX0);
+
// do the transport using the new engine
- Trk::ExtrapolationCode eCode = m_extrapolationEngine->extrapolate(ecc);
+ Trk::ExtrapolationCode eCode = m_extrapolationEngine->extrapolate(ecc);
+
+ // validation info
+ if (m_validationMode ) m_validationTool->saveISFParticleInfo(isp,ecc,eCode);
+
// resolve the different success states
if (eCode.isSuccess()) {
// this is a success - let's find out which one
@@ -249,11 +285,20 @@ ISF::ISFParticle* iFatras::TransportEngine::process( const ISF::ISFParticle& isp
eCode,
ecc.endParameters->position(),
ecc.endParameters->momentum(),
- ecc.endParameters->charge(),
- ecc.endParameters->position().mag()/CLHEP::c_light); //!< TODO update with pathLength
+ ecc.pathLength/CLHEP::c_light,
+ ecc.nextGeometrySignature);
// memory cleanup
- delete ecc.endParameters;
+ // delete ecc.endParameters;
// no hit creation -> return can be done now
+ if (rParticle && m_validationMode) {
+ // 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
+ rParticle->setUserInformation(validInfo);
+ }
+
return rParticle;
} else {
ATH_MSG_WARNING( "Error code " << eCode.toString() << " running the Extrapolator for neutral particle!" );
@@ -277,21 +322,32 @@ ISF::ISFParticle* iFatras::TransportEngine::process( const ISF::ISFParticle& isp
materialLimitX0,
materialLimitL0);
ecc.addConfigurationMode(Trk::ExtrapolationMode::CollectSensitive);
- ecc.materialProcess = matLim.process;
+ ecc.materialProcess = materialProcess;
+ if (ecc.materialLimitL0 != -1)
+ ecc.materialL0 = materialL0;
+ ecc.addConfigurationMode(Trk::ExtrapolationMode::StopWithMaterialLimitL0);
+ if (ecc.materialLimitX0 != -1)
+ ecc.materialX0 = materialX0;
+ ecc.addConfigurationMode(Trk::ExtrapolationMode::StopWithMaterialLimitX0);
// do the transport using the new engine
// [B - 1] extrapoalte
Trk::ExtrapolationCode eCode = m_extrapolationEngine->extrapolate(ecc);
+
+ // validation info
+ if (m_validationMode ) m_validationTool->saveISFParticleInfo(isp,ecc,eCode);
+
// resolve the different success states
if (eCode.isSuccess()) {
// this is a success - let's find out which one
ATH_MSG_VERBOSE( "[ fatras transport ] Successfully run extrapolation for charged particle.");
- // either give a particle back (at boundary) or the sub tools took care of it
+
+ // either give a particle back (at boundary) or the sub tools took care of it
rParticle = handleExtrapolationResult(isp,
eCode,
ecc.endParameters->position(),
ecc.endParameters->momentum(),
- ecc.endParameters->charge(),
- ecc.endParameters->position().mag()/CLHEP::c_light); //!< TODO update with pathLength
+ ecc.pathLength/CLHEP::c_light,
+ ecc.nextGeometrySignature);
// memory cleanup
delete ecc.endParameters;
} else {
@@ -323,6 +379,16 @@ ISF::ISFParticle* iFatras::TransportEngine::process( const ISF::ISFParticle& isp
}
}
}
+
+ if (rParticle && m_validationMode) {
+ // 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
+ rParticle->setUserInformation(validInfo);
+ }
+
// return what you have to the ISF
return rParticle;
}
@@ -333,30 +399,34 @@ ISF::ISFParticle* iFatras::TransportEngine::handleExtrapolationResult(const ISF:
Trk::ExtrapolationCode eCode,
const Amg::Vector3D& position,
const Amg::Vector3D& momentum,
- double charge,
- double stime)
+ double stime,
+ Trk::GeometrySignature nextGeoID)
{
+
// the return particle
ISF::ISFParticle* rParticle = 0;
// switch the extrapolation code returns
if (eCode == Trk::ExtrapolationCode::SuccessBoundaryReached){
ATH_MSG_VERBOSE( "[ fatras transport ] Successfully reached detector boundary with the particle -> return particle at boundary.");
// create the updated particle at the end of processing step
+ AtlasDetDescr::AtlasRegion geoID=AtlasDetDescr::AtlasRegion(5);
+ if (nextGeoID<99) geoID = AtlasDetDescr::AtlasRegion(nextGeoID);
rParticle = new ISF::ISFParticle(position,
momentum,
isp.mass(),
isp.charge(),
- stime+isp.timeStamp(),
isp.pdgCode(),
+ stime+isp.timeStamp(),
isp,
isp.barcode());
+ rParticle->setNextGeoID(geoID);
} else if (eCode == Trk::ExtrapolationCode::SuccessPathLimit){
ATH_MSG_VERBOSE( "[ fatras transport ] Successfully reached path limit for the particle -> decay & return 0.");
// call the decayer tool
- m_particleDecayHelper->decay(isp);
+ m_particleDecayHelper->decay(isp,position,momentum,stime+isp.timeStamp());
} else if (eCode == Trk::ExtrapolationCode::SuccessMaterialLimit )
ATH_MSG_VERBOSE( "[ fatras transport ] Successfully reached material limit for the particle -> return 0.");
+
return rParticle;
}
-
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/TransportTool.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/TransportTool.cxx
index 46be49de12d7a554275be585d72e020f73b04b56..980323fadf6ff8e6c4c8ab789d31b20615f7310c 100755
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/TransportTool.cxx
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/TransportTool.cxx
@@ -62,7 +62,7 @@ iFatras::TransportTool::TransportTool( const std::string& t,
m_errorPropagation(false),
m_hitsOff(false)
{
- declareInterface<ITransportTool>(this);
+ declareInterface<ISF::IParticleProcessor>(this);
// validation output section
declareProperty( "ValidationOutput", m_validationOutput );
@@ -108,7 +108,7 @@ iFatras::TransportTool::initialize()
if (retrieveTool<iFatras::ISimHitCreator>(m_simHitCreatorID).isFailure())
return StatusCode::FAILURE;
if (retrieveTool<iFatras::ISimHitCreator>(m_simHitCreatorMS).isFailure())
- return StatusCode::FAILURE;
+ return StatusCode::FAILURE;
if (retrieveTool<Trk::ITimedExtrapolator>(m_extrapolator).isFailure())
return StatusCode::FAILURE;
if (retrieveTool<iFatras::IParticleDecayHelper>(m_particleDecayHelper).isFailure())
@@ -245,7 +245,7 @@ ISF::ISFParticle* iFatras::TransportTool::process( const ISF::ISFParticle& isp)
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);
+ m_particleDecayHelper->decay(isp,isp.position(),isp.momentum(),isp.timeStamp());
}
// validation mode - for all particle registered into stack
@@ -257,9 +257,15 @@ ISF::ISFParticle* iFatras::TransportTool::process( const ISF::ISFParticle& isp)
return 0;
}
- // presample interactions
+ // presample interactions if not done already
Trk::PathLimit pathLim(-1.,0);
- if (absPdg!=999 && pHypothesis<99) pathLim = m_samplingTool->sampleProcess(mom,isp.charge(),pHypothesis);
+ ISF::MaterialPathInfo* matLimit = isp.getUserInformation() ? isp.getUserInformation()->materialLimit() : 0;
+ 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(isp.momentum().mag(),isp.charge(),pHypothesis);
+ }
// use extrapolation with path limit - automatic exit at subdetector boundary
// NB: don't delete eParameters, it's memory is managed inside the extrapolator
@@ -369,7 +375,7 @@ ISF::ISFParticle* iFatras::TransportTool::process( const ISF::ISFParticle& isp)
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);
+ m_particleDecayHelper->decay(*uisp,uisp->position(),uisp->momentum(),uisp->timeStamp());
}
delete uisp;
return 0;
diff --git a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/components/ISF_FatrasTools_entries.cxx b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/components/ISF_FatrasTools_entries.cxx
index 18734e63a1f14ba488c57259b62ca9d4b0dd912d..3f455cd6b4385c7b07355e5e0bf03aa050d21334 100755
--- a/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/components/ISF_FatrasTools_entries.cxx
+++ b/Simulation/ISF/ISF_Fatras/ISF_FatrasTools/src/components/ISF_FatrasTools_entries.cxx
@@ -5,19 +5,24 @@
#include "ISF_FatrasTools/McMaterialEffectsEngine.h"
#include "ISF_FatrasTools/HadIntProcessorParametric.h"
#include "ISF_FatrasTools/McEnergyLossUpdator.h"
-#include "ISF_FatrasTools/McBetheHeitlerEnergyLossUpdator.h"
#include "ISF_FatrasTools/PhotonConversionTool.h"
#include "ISF_FatrasTools/ProcessSamplingTool.h"
#include "ISF_FatrasTools/PhysicsValidationTool.h"
+#include "ISF_FatrasTools/MultipleScatteringSamplerHighland.h"
+#include "ISF_FatrasTools/MultipleScatteringSamplerGaussianMixture.h"
+#include "ISF_FatrasTools/MultipleScatteringSamplerGeneralMixture.h"
+
DECLARE_NAMESPACE_TOOL_FACTORY( iFatras , TransportTool )
DECLARE_NAMESPACE_TOOL_FACTORY( iFatras , TransportEngine )
DECLARE_NAMESPACE_TOOL_FACTORY( iFatras , McMaterialEffectsUpdator )
DECLARE_NAMESPACE_TOOL_FACTORY( iFatras , McMaterialEffectsEngine )
DECLARE_NAMESPACE_TOOL_FACTORY( iFatras , HadIntProcessorParametric )
DECLARE_NAMESPACE_TOOL_FACTORY( iFatras , McEnergyLossUpdator )
-DECLARE_NAMESPACE_TOOL_FACTORY( iFatras , McBetheHeitlerEnergyLossUpdator )
DECLARE_NAMESPACE_TOOL_FACTORY( iFatras , PhotonConversionTool )
+DECLARE_NAMESPACE_TOOL_FACTORY( iFatras , MultipleScatteringSamplerGaussianMixture )
+DECLARE_NAMESPACE_TOOL_FACTORY( iFatras , MultipleScatteringSamplerGeneralMixture )
+DECLARE_NAMESPACE_TOOL_FACTORY( iFatras , MultipleScatteringSamplerHighland )
DECLARE_NAMESPACE_TOOL_FACTORY( iFatras , ProcessSamplingTool )
DECLARE_NAMESPACE_TOOL_FACTORY( iFatras , PhysicsValidationTool )