ATLAS Offline Software
McMaterialEffectsUpdator.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // class header
6 #include "McMaterialEffectsUpdator.h"
7 
8 // Gaudi Kernel
9 #include "GaudiKernel/IRndmGenSvc.h"
10 #include "GaudiKernel/RndmGenerators.h"
11 #include "GaudiKernel/DataSvc.h"
12 #include "GaudiKernel/SmartDataPtr.h"
13 // ISF includes
14 #include "ISF_Interfaces/IParticleBroker.h"
15 #include "ISF_Interfaces/ITruthSvc.h"
16 #include "ISF_Event/ISFParticle.h"
17 #include "ISF_Event/ISFTruthIncident.h"
18 #include "ISF_Event/ParticleClipboard.h"
19 #include "ISF_Event/ParticleUserInformation.h"
20 #include "ISF_FatrasInterfaces/IParticleDecayHelper.h"
21 #include "TruthUtils/HepMCHelpers.h"
22 // iFatras
23 #include "ISF_FatrasInterfaces/IHadronicInteractionProcessor.h"
24 #include "ISF_FatrasInterfaces/IProcessSamplingTool.h"
25 #include "ISF_FatrasInterfaces/IPhysicsValidationTool.h"
26 #include "ISF_FatrasInterfaces/IPhotonConversionTool.h"
27 // Trk inlcude
28 #include "TrkExInterfaces/IEnergyLossUpdator.h"
29 #include "TrkExInterfaces/IMultipleScatteringUpdator.h"
30 #include "TrkParameters/TrackParameters.h"
31 #include "TrkEventPrimitives/ParamDefs.h"
32 #include "TrkSurfaces/Surface.h"
33 #include "TrkGeometry/Layer.h"
34 #include "TrkGeometry/TrackingGeometry.h"
35 #include "TrkGeometry/TrackingVolume.h"
36 #include "TrkGeometry/MaterialProperties.h"
37 #include "TrkVolumes/CylinderVolumeBounds.h"
38 #include "TrkMaterialOnTrack/EnergyLoss.h"
39 #include "TrkMaterialOnTrack/MaterialEffectsOnTrack.h"
40 // CLHEP
41 #include "CLHEP/Units/SystemOfUnits.h"
42 #include "CLHEP/Matrix/Vector.h"
43 #include "CLHEP/Random/RandFlat.h"
44 #include "CLHEP/Random/RandLandau.h"
45 #include "CLHEP/Random/RandGaussZiggurat.h"
46 // Validation mode - TTree includes
47 #include "TTree.h"
48 #include "GaudiKernel/ITHistSvc.h"
49 // STD
50 #include <cmath>
51 
52 // constructor
53 iFatras::McMaterialEffectsUpdator::McMaterialEffectsUpdator(const std::string& t, const std::string& n, const IInterface* p) :
54  base_class(t,n,p),
55  m_eLoss(true),
56  m_eLossUpdator("iFatras::McEnergyLossUpdator/FatrasEnergyLossUpdator"),
57  m_ms(true),
58  m_msUpdator("Trk::MultipleScatteringUpdator/AtlasMultipleScatteringUpdator"),
59  m_conversionTool("iFatras::McConversionCreator/FatrasConversionCreator"),
60  m_processCode(3),
61  m_hadInt(true),
62  //m_hadIntProcessor("iFatras::HadIntProcessorParametric/FatrasHadIntProcessorParametric"),
63  m_hadIntProcessor("iFatras::G4HadIntProcessor/FatrasG4HadIntProcessor"),
64  m_particleDecayer(),
65  m_minimumMomentum(50.*CLHEP::MeV),
66  m_minimumBremPhotonMomentum(50.*CLHEP::MeV),
67  m_parametricScattering(false),
68  m_use_msUpdator(false),
69  m_uniformHertzDipoleAngle(false),
70  m_referenceMaterial(true),
71  m_bendingCorrection(false),
72  m_rndGenSvc("AtDSFMTGenSvc", n),
73  m_randomEngine(nullptr),
74  m_randomEngineName("FatrasRnd"),
75  m_recordedBremPhotons(0),
76  m_currentSample(-1),
77  m_recordEnergyDeposition(false),
78  m_layerIndexCaloSampleMapName("LayerIndexCaloSampleMap"),
79  m_layerIndexCaloSampleMap(nullptr),
80  m_projectionFactor(sqrt(2.)),
81  m_validationMode(false),
82  m_validationTool(""),
83  m_validationTreeName("FatrasMaterialEffects"),
84  m_validationTreeDescription("Validation output from the McMaterialEffectsUpdator"),
85  m_validationTreeFolder("/val/FatrasSimulationMaterial"),
86  m_validationTree(nullptr),
87  m_layerIndex(0),
88  m_tInX0(0.),
89  m_thetaMSproj(0.),
90  m_thetaMSphi(0.),
91  m_thetaMStheta(0.),
92  m_deltaP(0.),
93  m_deltaPsigma(0.),
94  m_bremValidation(false),
95  m_bremValidationTreeName("FatrasBremPhotons"),
96  m_bremValidationTreeDescription("Validation output from the McMaterialEffectsUpdator"),
97  m_bremValidationTreeFolder("/val/FatrasBremPhotons"),
98  m_bremValidationTree(nullptr),
99  m_bremPointX(0.),
100  m_bremPointY(0.),
101  m_bremPointR(0.),
102  m_bremPointZ(0.),
103  m_bremMotherEnergy(0.),
104  m_bremPhotonEnergy(0.),
105  m_bremPhotonAngle(0.),
106  m_edValidation(false),
107  m_edValidationTreeName("FatrasEnergyInCaloDeposit"),
108  m_edValidationTreeDescription("Validation output from the McMaterialEffectUpdator"),
109  m_edValidationTreeFolder("/val/FatrasEnergyInCaloDeposit"),
110  m_edValidationTree(nullptr),
111  m_edLayerIntersectX(0.),
112  m_edLayerIntersectY(0.),
113  m_edLayerIntersectZ(0.),
114  m_edLayerIntersectR(0.),
115  m_edLayerEnergyDeposit(0.),
116  m_edLayerSample(0.),
117  m_oneOverThree(1./3.),
118  m_particleBroker("ISF_ParticleParticleBroker", n),
119  m_truthRecordSvc("ISF_TruthRecordSvc", n),
120  m_isp(nullptr),
121  m_layer(nullptr),
122  m_matProp(nullptr)
123 {
124  // the tool parameters -----------------------------------------------------
125  declareProperty("EnergyLoss" , m_eLoss);
126  declareProperty("EnergyLossUpdator" , m_eLossUpdator);
127  declareProperty("MultipleScattering" , m_ms);
128  declareProperty("MultipleScatteringUpdator" , m_msUpdator);
129  declareProperty("ProcessSamplingTool" , m_samplingTool);
130  declareProperty("PhotonConversionTool" , m_conversionTool);
131  // tool handle for the particle decayer
132  declareProperty( "ParticleDecayHelper" , m_particleDecayer );
133  // MC Truth Properties
134  declareProperty("BremProcessCode" , m_processCode, "MCTruth Physics Process Code");
135  // the steering --------------------------------------------------------------
136  declareProperty("HadronicInteraction" , m_hadInt);
137  declareProperty("HadronicInteractionProcessor" , m_hadIntProcessor);
138  // the steering --------------------------------------------------------------
139  declareProperty("ParametericScattering" , m_parametricScattering);
140  declareProperty("UseMultipleScatteringUpdator" , m_use_msUpdator);
141  declareProperty("MomentumCut" , m_minimumMomentum);
142  declareProperty("MinimumBremPhotonMomentum" , m_minimumBremPhotonMomentum);
143  // the steering --------------------------------------------------------------
144  declareProperty("ReferenceMaterial" , m_referenceMaterial);
145  declareProperty("BendingCorrection" , m_bendingCorrection);
146  // validation mode ------------------------------------------------------------
147  declareProperty("RecordEnergyDeposition" , m_recordEnergyDeposition);
148  declareProperty("LayerIndexCaloSampleMap" , m_layerIndexCaloSampleMapName);
149  // validation mode ------------------------------------------------------------
150  declareProperty("ValidationMode" , m_validationMode);
151  declareProperty("PhysicsValidationTool" , m_validationTool);
152  declareProperty("BremPhotonValidation" , m_bremValidation);
153  declareProperty("EnergyDepositValidation" , m_edValidation);
154  declareProperty("RandomNumberService" , m_rndGenSvc , "Random number generator");
155  declareProperty("RandomStreamName" , m_randomEngineName , "Name of the random number stream");
156 
157  // ISF Services and Tools
158  declareProperty("ParticleBroker" , m_particleBroker , "ISF Particle Broker Svc");
159  declareProperty("TruthRecordSvc" , m_truthRecordSvc , "ISF Particle Truth Svc");
160 }
161 
162 // destructor
163 iFatras::McMaterialEffectsUpdator::~McMaterialEffectsUpdator()
164 = default;
165 
166 // Athena standard methods
167 // initialize
168 StatusCode iFatras::McMaterialEffectsUpdator::initialize()
169 {
170 
171  ATH_MSG_DEBUG( "initialize()" );
172 
173  // retrieve the process sampling tool
174  if (m_samplingTool.retrieve().isFailure()){
175  ATH_MSG_FATAL( "Could not retrieve " << m_samplingTool );
176  return StatusCode::FAILURE;
177  } else
178  ATH_MSG_VERBOSE( "Successfully retrieved " << m_samplingTool );
179 
180  // retrieve the energy loss updator
181  if (m_eLoss){
182  if (m_eLossUpdator.retrieve().isFailure()){
183  ATH_MSG_FATAL( "Could not retrieve " << m_eLossUpdator );
184  return StatusCode::FAILURE;
185  } else
186  ATH_MSG_VERBOSE( "Successfully retrieved " << m_eLossUpdator );
187  } else {
188  m_eLossUpdator.disable();
189  }
190 
191  // retrieve the multiple scattering updator tool
192  if (m_ms && m_use_msUpdator){
193  if (m_msUpdator.retrieve().isFailure()){
194  ATH_MSG_FATAL( "Could not retrieve " << m_msUpdator );
195  return StatusCode::FAILURE;
196  } else
197  ATH_MSG_VERBOSE( "Successfully retrieved " << m_msUpdator );
198  } else {
199  m_msUpdator.disable();
200  }
201 
202  // retrieve the photon conversion tool
203  if (m_conversionTool.retrieve().isFailure()){
204  ATH_MSG_FATAL( "Could not retrieve " << m_conversionTool );
205  return StatusCode::FAILURE;
206  } else
207  ATH_MSG_VERBOSE( "Successfully retrieved " << m_conversionTool );
208 
209  // retrieve the hadronic interaction tool
210  if (m_hadInt){
211  if (m_hadIntProcessor.retrieve().isFailure()){
212  ATH_MSG_FATAL( "Could not retrieve " << m_hadIntProcessor );
213  return StatusCode::FAILURE;
214  } else
215  ATH_MSG_VERBOSE( "Successfully retrieved " << m_hadIntProcessor );
216  } else {
217  m_hadIntProcessor.disable();
218  }
219 
220  // get the random generator serice
221  if (m_rndGenSvc.retrieve().isFailure()){
222  ATH_MSG_FATAL( "Could not retrieve " << m_rndGenSvc );
223  return StatusCode::FAILURE;
224  } else
225  ATH_MSG_VERBOSE( "Successfully retrieved " << m_rndGenSvc );
226 
227  //Get own engine with own seeds:
228  m_randomEngine = m_rndGenSvc->GetEngine(m_randomEngineName);
229  if (!m_randomEngine) {
230  ATH_MSG_FATAL( "Could not get random engine '" << m_randomEngineName << "'" );
231  return StatusCode::FAILURE;
232  }
233 
234  //We can use conditions when the key is not empty
235  m_useConditions=!m_trackingGeometryReadKey.key().empty();
236 
237  // get the TrackingGeometrySvc
238  if (!m_useConditions) {
239  if (m_trackingGeometrySvc.retrieve().isSuccess()) {
240  ATH_MSG_DEBUG("Successfully retrieved " << m_trackingGeometrySvc);
241  m_trackingGeometryName = m_trackingGeometrySvc->trackingGeometryName();
242  } else {
243  ATH_MSG_WARNING("Couldn't retrieve " << m_trackingGeometrySvc << ". ");
244  ATH_MSG_WARNING(" -> Trying to retrieve default '"
245  << m_trackingGeometryName << "' from DetectorStore.");
246  }
247  }
248 
249  // get the tracking geometry for layer lookup
250  // init the TrackingGeometryReadKey
251  ATH_CHECK(m_trackingGeometryReadKey.initialize(m_useConditions));
252 
253  // Particle decayer
254  if (m_particleDecayer.retrieve().isFailure()){
255  ATH_MSG_FATAL( "Could not retrieve " << m_particleDecayer );
256  return StatusCode::FAILURE;
257  } else
258  ATH_MSG_VERBOSE( "Successfully retrieved " << m_particleDecayer );
259 
260  ITHistSvc* tHistSvc = nullptr;
261  if (m_validationMode || m_bremValidation || m_edValidation){
262  // now register the Tree
263  if (service("THistSvc",tHistSvc).isFailure())
264  ATH_MSG_ERROR( "initialize() Could not find Hist Service -> Switching ValidationMode Off !" );
265  }
266 
267  // ISF Services
268  if (m_particleBroker.retrieve().isFailure()){
269  ATH_MSG_FATAL( "Could not retrieve " << m_particleBroker );
270  return StatusCode::FAILURE;
271  }
272  if (m_truthRecordSvc.retrieve().isFailure()){
273  ATH_MSG_FATAL( "Could not retrieve " << m_truthRecordSvc );
274  return StatusCode::FAILURE;
275  }
276 
277  // the validation setup -------------------------------- PART 1: General ----------------------------------
278 
279  // retrieve the physics validation tool
280  ATH_CHECK( m_validationTool.retrieve( DisableTool{ m_validationTool.empty() || !m_validationMode } ) );
281 
282  if (m_validationMode){
283  ATH_MSG_VERBOSE( "Booking material validation TTree ... " );
284 
285  // create the new Tree
286  m_validationTree = new TTree(m_validationTreeName.c_str(), m_validationTreeDescription.c_str());
287 
288  // counter for material effects
289  m_validationTree->Branch("LayerIndex" , &m_layerIndex, "layerIdx/I");
290  m_validationTree->Branch("PathInX0" , &m_tInX0, "tInX0/F");
291  m_validationTree->Branch("ThetaMS" , &m_thetaMSproj, "thetaMS/F");
292  m_validationTree->Branch("ThetaMSphi" , &m_thetaMSphi, "thetaMSphi/F");
293  m_validationTree->Branch("ThetaMStheta" , &m_thetaMStheta, "thetaMStheta/F");
294  m_validationTree->Branch("DeltaP" , &m_deltaP, "deltaP/F");
295  m_validationTree->Branch("DeltaPsigma" , &m_deltaPsigma, "deltaPsigma/F");
296 
297 
298  if ((tHistSvc->regTree(m_validationTreeFolder, m_validationTree)).isFailure()) {
299  ATH_MSG_ERROR( "initialize() Could not register the validation Tree -> Switching ValidationMode Off !" );
300  delete m_validationTree; m_validationTree = nullptr;
301  } else
302  ATH_MSG_INFO( "TTree for MaterialEffects validation booked." );
303 
304  }
305  // the validation setup -------------------------------- PART 2: Brem Photons -----------------------------
306  if (m_bremValidation){
307 
308  ATH_MSG_VERBOSE( "Booking bremstrahlung validation TTree ... " );
309 
310  // create the new Tree
311  m_bremValidationTree = new TTree(m_bremValidationTreeName.c_str(), m_bremValidationTreeDescription.c_str());
312 
313  // counter for bremstrahlung
314  m_bremValidationTree->Branch("BremPositionX" , &m_bremPointX, "bremX/F");
315  m_bremValidationTree->Branch("BremPositionY" , &m_bremPointY, "bremY/F");
316  m_bremValidationTree->Branch("BremPositionR" , &m_bremPointR, "bremR/F");
317  m_bremValidationTree->Branch("BremPositionZ" , &m_bremPointZ, "bremZ/F");
318  m_bremValidationTree->Branch("BremMotherEnergy" , &m_bremMotherEnergy, "bremMotherE/F");
319  m_bremValidationTree->Branch("BremPhotonEnergy" , &m_bremPhotonEnergy, "bremPhotonE/F");
320  m_bremValidationTree->Branch("BremPhotonAngle" , &m_bremPhotonAngle, "bremPhotondA/F");
321 
322  if ((tHistSvc->regTree(m_bremValidationTreeFolder, m_bremValidationTree)).isFailure()) {
323  ATH_MSG_ERROR("initialize() Could not register the validation Tree -> Switching ValidationMode Off !" );
324  delete m_bremValidationTree; m_bremValidationTree = nullptr;
325  } else
326  ATH_MSG_INFO( "TTree for Bremsstrahlung validation booked." );
327 
328  } // ------------- end of validation mode -----------------------------------------------------------------
329 
330  // the validation setup -------------------------------- PART 3: Brem Photons -----------------------------
331  if (m_edValidation){
332 
333  ATH_MSG_VERBOSE( "Booking Energy deposition validation TTree ... " );
334 
335  // create the new Tree
336  m_edValidationTree = new TTree(m_edValidationTreeName.c_str(), m_edValidationTreeDescription.c_str());
337 
338  // counter for boundary surfaces
339  m_edValidationTree->Branch("EdepositPositionX" , &m_edLayerIntersectX, "edX/F");
340  m_edValidationTree->Branch("EdepositPositionY" , &m_edLayerIntersectY, "edY/F");
341  m_edValidationTree->Branch("EdepositPositionR" , &m_edLayerIntersectR, "edR/F");
342  m_edValidationTree->Branch("EdepositPositionZ" , &m_edLayerIntersectZ, "edZ/F");
343  m_edValidationTree->Branch("Edeposit" , &m_edLayerEnergyDeposit, "ed/F");
344  m_edValidationTree->Branch("EdepositLayerSample" , &m_edLayerSample, "edLayerSample/F");
345 
346  if ((tHistSvc->regTree(m_edValidationTreeFolder, m_edValidationTree)).isFailure()) {
347  ATH_MSG_ERROR("initialize() Could not register the validation Tree -> Switching ValidationMode Off !" );
348  delete m_edValidationTree; m_edValidationTree = nullptr;
349  } else
350  ATH_MSG_INFO( "TTree for Energy deposition validation booked." );
351 
352 
353 
354  } // ------------- end of validation mode -----------------------------------------------------------------
355  ATH_MSG_DEBUG( "finalize() successful" );
356  return StatusCode::SUCCESS;
357 }
358 
359 // finalize
360 StatusCode iFatras::McMaterialEffectsUpdator::finalize()
361 {
362 
363  ATH_MSG_INFO( " ---------- Statistics output -------------------------- " );
364  ATH_MSG_INFO( " Minimum energy cut for brem photons : " << m_minimumBremPhotonMomentum );
365  ATH_MSG_INFO( " Brem photons (above cut, recorded) : " << m_recordedBremPhotons );
366 
367  ATH_MSG_DEBUG( "finalize() successful" );
368  return StatusCode::SUCCESS;
369 }
370 
371 std::unique_ptr<Trk::TrackParameters>
372 iFatras::McMaterialEffectsUpdator::update(
373  const Trk::TrackParameters* parm,
374  const Trk::Layer& lay,
375  Trk::TimeLimit& timeLim,
376  Trk::PathLimit& pathLim,
377  Trk::GeometrySignature /*geoID*/,
378  Trk::PropDirection dir,
379  Trk::ParticleHypothesis particle) const
380 {
381 
382  // the call should work for particles created on the layer, too
383  // + additional parameter for material scaling
384 
385  double traversedMatFraction = 0.;
386 
387  m_layer = &lay;
388 
389  // get parent particle
390  const ISF::ISFParticle *parent = ISF::ParticleClipboard::getInstance().getParticle();
391  // something is seriously wrong if there is no parent particle
392  assert(parent);
393 
394  m_isp = parent;
395 
396  // get the material properties
397  m_matProp = nullptr;
398 
399  if (m_referenceMaterial){
400  // very first go : get the Surface
401  const Trk::Surface &parmSurface = parm->associatedSurface();
402  // only go on if the parmSurface has an associatedDetectorElement
403  if (parmSurface.associatedDetectorElement()){
404  // first get the LayerMaterialProperties
405  const Trk::LayerMaterialProperties* layerMaterial = m_layer->layerMaterialProperties();
406  // assign the material properties
407  m_matProp = layerMaterial ? layerMaterial->fullMaterial( parm->position() ) : nullptr;
408  }
409  }
410 
411 
412  return updateInLay(parent,parm,traversedMatFraction,timeLim, pathLim,dir, particle);
413 }
414 
415 std::unique_ptr<Trk::TrackParameters>
416 iFatras::McMaterialEffectsUpdator::updateInLay(
417  const ISF::ISFParticle* isp,
418  const Trk::TrackParameters* inParm,
419  double& matFraction,
420  Trk::TimeLimit& timeLim,
421  Trk::PathLimit& pathLim,
422  Trk::PropDirection dir,
423  Trk::ParticleHypothesis particle) const
424 {
425 
426  if (!inParm)
427  {
428  ATH_MSG_WARNING("Trk::TrackParameters parm is null -- will not proceed. Returning nullptr.");
429  return nullptr;
430  }
431 
432  std::unique_ptr<Trk::TrackParameters> currPar = inParm->uniqueClone();
433  // recalculate if missing
434  m_matProp = m_matProp ? m_matProp : m_layer->fullUpdateMaterialProperties(*currPar);
435  if (!m_matProp)
436  {
437  ATH_MSG_WARNING("Something went wrong -- m_matProp is missing but Trk::TrackParameters is not null! Returning nullptr.");
438  return nullptr;
439  }
440  double pathCorrection =
441  m_layer->surfaceRepresentation().normal().dot(currPar->momentum()) != 0
442  ? std::abs(m_layer->surfaceRepresentation().pathCorrection(currPar->position(), dir * (currPar->momentum())))
443  : 1.;
444 
445  // check if a bending correction has to be applied
446  if (m_bendingCorrection) {
447  ATH_MSG_WARNING("BendingCorrection not implemented!! (McMaterialEffectsUpdator.cxx)");
448  }
449 
450  //--------------------------------------------------------------------------------------------------
451  if (msgLvl(MSG::VERBOSE) && dir != Trk::PropDirection::anyDirection){
452  const Trk::TrackingVolume* enclosingVolume = m_layer->enclosingTrackingVolume();
453  std::string volumeName = enclosingVolume ? enclosingVolume->volumeName() : "Unknown";
454  double layerR = m_layer->surfaceRepresentation().bounds().r();
455  double layerZ = m_layer->surfaceRepresentation().center().z();
456  ATH_MSG_VERBOSE( " [M] Material effects update() on layer with index "<< m_layer->layerIndex() );
457  ATH_MSG_VERBOSE( " -> path/X0 on layer [r/z] = [ " << layerR << " / " << layerZ << " ] :"
458  << pathCorrection*m_matProp->thicknessInX0() );
459  ATH_MSG_VERBOSE( " -> path correctin factor = " << pathCorrection );
460  }
461  //--------------------------------------------------------------------------------------------------
462 
463  // validation mode ---------------------------------------------------------------
464  if (m_validationMode) {
465  m_tInX0 = (1-matFraction)*pathCorrection * m_matProp->thicknessInX0();
466  m_layerIndex = m_layer->layerIndex().value();
467  m_validationTree->Fill();
468  }
469  // -------------------------------------------------------------------------------
470 
471  // prepare material collection
472  const Trk::MaterialProperties* extMatProp = nullptr;
473  double dInL0 = 0.;
474  extMatProp = dynamic_cast<const Trk::MaterialProperties*>(m_matProp);
475  dInL0 = extMatProp ? (1-matFraction)*pathCorrection*extMatProp->thicknessInL0() :
476  (1-matFraction)*pathCorrection*m_matProp->thicknessInX0()/0.37/m_matProp->averageZ();
477 
478  // figure out if particle stopped in the layer and recalculate path limit
479  int iStatus = 0;
480  double dX0 = (1.-matFraction)*pathCorrection*m_matProp->thicknessInX0();
481 
482  if (particle == Trk::geantino || particle == Trk::nonInteracting || dX0 == 0.) {
483  // non-interacting - pass the input back after checks
484  pathLim.updateMat(dX0, m_matProp->averageZ(), dInL0);
485  // register particle if not in the stack already
486  if (isp != m_isp) {
487  ISF::TruthBinding *regTruthBinding{};
488  if (isp->getTruthBinding()) {
489  regTruthBinding = new ISF::TruthBinding(*(isp->getTruthBinding()));
490  }
491  else {
492  ATH_MSG_WARNING("Incomming ISParticle had no TruthBinding " << *isp);
493  regTruthBinding = new ISF::TruthBinding(nullptr, nullptr, nullptr);
494  }
495  HepMcParticleLink *regHMPL{};
496  if (isp->getParticleLink()) {
497  regHMPL = new HepMcParticleLink(*(isp->getParticleLink()));
498  }
499  else {
500  ATH_MSG_WARNING("Incomming ISParticle had no ParticleLink: " << *isp);
501  regHMPL = new HepMcParticleLink(isp->id(), 0, HepMcParticleLink::IS_POSITION, HepMcParticleLink::IS_ID);
502  }
503  ISF::ISFParticle* regisp = new ISF::ISFParticle(isp->position(),
504  currPar->momentum(),
505  isp->mass(),
506  isp->charge(),
507  isp->pdgCode(),
508  HepMC::status(isp),
509  isp->timeStamp(),
510  *m_isp,
511  HepMC::uniqueID(isp),
512  HepMC::barcode(isp), // FIXME barcode-based
513  regTruthBinding,
514  regHMPL
515  );
516  // add presampled process info
517  if (isp->getUserInformation() && isp->getUserInformation()->materialLimit()) {
518  const ISF::MaterialPathInfo* matLim = isp->getUserInformation()->materialLimit();
519  ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
520  validInfo->setMaterialLimit(
521  matLim->process, matLim->dMax, matLim->process == 121 ? pathLim.l0Collected : pathLim.x0Collected);
522  regisp->setUserInformation(validInfo);
523  }
524  // in the validation mode, add process info
525  if (m_validationMode) {
526  ISF::ParticleUserInformation* validInfo = regisp->getUserInformation();
527  if (!validInfo) {
528  validInfo = new ISF::ParticleUserInformation();
529  regisp->setUserInformation(validInfo);
530  }
531  if (isp->getUserInformation()) {
532  validInfo->setProcess(isp->getUserInformation()->process());
533  } else {
534  validInfo->setProcess(-2); // signal problem in the validation chain}
535  }
536  if (isp->getUserInformation()) {
537  validInfo->setGeneration(isp->getUserInformation()->generation());
538  } else {
539  validInfo->setGeneration(-1); // signal problem in the validation chain}
540  }
541  regisp->setUserInformation(validInfo);
542  }
543  // Check that the returned ISFParticle has a valid TruthBinding
544  // before pushing into the particle broker
545  if (!regisp->getTruthBinding()) {
546  ATH_MSG_ERROR("Could not retrieve TruthBinding from non-interacting ISFParticle "<< *regisp);
547  }
548  m_particleBroker->push(regisp, m_isp);
549  }
550  if (isp != m_isp) {
551  delete isp;
552  return nullptr;
553  }
554  // non-interacting - pass the input back
555  return currPar;
556  } else {
557  if (pathLim.x0Max > 0 && pathLim.process < 100 && pathLim.x0Collected + dX0 >= pathLim.x0Max) { // elmg. interaction
558  float x0rem = pathLim.x0Max - pathLim.x0Collected;
559  dInL0 *= x0rem > 0 ? x0rem / dX0 : 1.;
560  if (x0rem > 0)
561  dX0 = x0rem;
562  iStatus = 1;
563  } else if (pathLim.x0Max > 0 && extMatProp && pathLim.process > 100 &&
564  pathLim.l0Collected + dInL0 >= pathLim.x0Max) { // hadronic interaction
565  float l0rem = pathLim.x0Max - pathLim.l0Collected;
566  dX0 *= l0rem > 0 ? l0rem / dInL0 : 1.;
567  if (l0rem > 0)
568  dInL0 = l0rem;
569  iStatus = 2;
570  }
571  pathLim.updateMat(dX0, m_matProp->averageZ(), dInL0);
572  }
573 
574  // get the kinematics
575  double p = currPar->momentum().mag();
576 
577  // radiation and ionization preceed the presampled interaction (if any)
578  // the updatedParameters - first a deep copy
579  AmgVector(5) updatedParameters(currPar->parameters());
580  std::unique_ptr<Trk::TrackParameters> updated = nullptr;
581 
582  if (m_eLoss && particle == Trk::electron && dX0 > 0.) {
583 
584  double matTot = (1 - matFraction) * pathCorrection * m_matProp->thicknessInX0();
585  Amg::Vector3D edir = currPar->momentum().unit();
586  radiate(isp, updatedParameters, currPar->position(), edir, timeLim, dX0, matFraction, matTot, dir, particle);
587  // update parent parameters:
588  updated = currPar->associatedSurface().createUniqueTrackParameters(updatedParameters[0],
589  updatedParameters[1],
590  updatedParameters[2],
591  updatedParameters[3],
592  updatedParameters[4],
593  std::nullopt);
594 
595  currPar = std::move(updated);
596 
597  if (currPar->momentum().mag() < m_minimumMomentum) {
598  if (isp != m_isp) {
599  delete isp;
600  }
601  return nullptr;
602  }
603  } else {
604  matFraction += dX0 / pathCorrection / m_matProp->thicknessInX0();
605  }
606 
607  if (isp->charge() != 0 && dX0 > 0.) {
608  if (m_eLoss) {
609  ionize(*currPar, updatedParameters, dX0, dir, particle);
610  }
611 
612  if (m_ms && m_matProp->thicknessInX0() > 0.) {
613  double sigmaMSproj =
614  (m_use_msUpdator && m_msUpdator)
615  ? sqrt(m_msUpdator->sigmaSquare(*m_matProp, p, dX0 / m_matProp->thicknessInX0(), particle))
616  : msSigma(dX0, p, particle);
617 
618  multipleScatteringUpdate(*currPar, updatedParameters, sigmaMSproj);
619  }
620 
621  // use the manipulator to update the track parameters -------> get rid of 0!
622  updated = currPar->associatedSurface().createUniqueTrackParameters(updatedParameters[0],
623  updatedParameters[1],
624  updatedParameters[2],
625  updatedParameters[3],
626  updatedParameters[4],
627  std::nullopt);
628  currPar = std::move(updated);
629 
630  if (currPar->momentum().mag() < m_minimumMomentum) {
631  if (isp != m_isp) {
632  delete isp;
633  }
634  return nullptr;
635  }
636  }
637 
638  if (iStatus == 1 || iStatus == 2) { // interaction
639  ISF::ISFParticleVector childs;
640 
641  if (iStatus == 1) {
642  childs = interactLay(isp, timeLim.time, *currPar, particle, pathLim.process); // Registers TruthIncident internally
643  } else {
644  if (extMatProp) {
645  childs = m_hadIntProcessor->doHadIntOnLayer(
646  isp, timeLim.time, currPar->position(), currPar->momentum(), &extMatProp->material(), particle);
647  } else {
648  childs = m_hadIntProcessor->doHadIntOnLayer(
649  isp, timeLim.time, currPar->position(), currPar->momentum(), nullptr, particle);
650  }
651  }
652  // save info for locally created particles
653  if (m_validationMode && m_validationTool.isEnabled() && !childs.empty() && isp != m_isp) {
654  ATH_MSG_VERBOSE(" saving interaction info for locally produced particle " << isp->pdgCode());
655  m_validationTool->saveISFParticleInfo(*isp, pathLim.process, currPar.get(), timeLim.time, pathLim.x0Max);
656  }
657 
658  // loop over childrens and decide about their fate
659 
660  for (unsigned ic = 0; ic < childs.size(); ic++) {
661  double mom = childs[ic]->momentum().mag();
662 
663  if (mom < m_minimumMomentum) {
664  continue;
665  }
666  Trk::ParticleHypothesis pHypothesis =
667  m_pdgToParticleHypothesis.convert(childs[ic]->pdgCode(), childs[ic]->charge());
668  auto cparm = std::make_unique<Trk::CurvilinearParameters>(childs[ic]->position(), childs[ic]->momentum(), childs[ic]->charge());
669  Trk::PathLimit pLim = m_samplingTool->sampleProcess(m_randomEngine, mom, childs[ic]->charge(), pHypothesis);
670 
671  // TODO sample decays and save the material collection & path limits at the exit from the layer
672  // (ISFFatrasParticle ?)
673 
674  // material fraction : flip if direction of propagation changed
675  double ci = m_layer->surfaceRepresentation().normal().dot(currPar->momentum().unit());
676  double co = m_layer->surfaceRepresentation().normal().dot(childs[ic]->momentum().unit());
677  double remMat = ci * co < 0 ? (1 - matFraction) : matFraction;
678 
679  // loop back
680  std::unique_ptr<Trk::TrackParameters> uPar =
681  updateInLay(childs[ic], cparm.get(), remMat, timeLim, pLim, dir, pHypothesis);
682  if (uPar)
683  ATH_MSG_VERBOSE("Check this: parameters should be dummy here " << isp->pdgCode() << ","
684  << uPar->position());
685  }
686  if (!childs.empty()) { // assume that interaction processing failed if no children
687  if (isp!=m_isp) {
688  delete isp;
689  }
690  return nullptr;
691  }
692  }
693 
694  // register particle if not in the stack already
695  if (isp != m_isp) {
696  ISF::TruthBinding *regTruthBinding{};
697  if (isp->getTruthBinding()) {
698  regTruthBinding = new ISF::TruthBinding(*(isp->getTruthBinding()));
699  }
700  else {
701  ATH_MSG_WARNING("Incomming ISParticle had no TruthBinding " << *isp);
702  regTruthBinding = new ISF::TruthBinding(nullptr, nullptr, nullptr);
703  }
704  HepMcParticleLink *regHMPL{};
705  if (isp->getParticleLink()) {
706  regHMPL = new HepMcParticleLink(*(isp->getParticleLink()));
707  }
708  else {
709  ATH_MSG_WARNING("Incomming ISParticle had no ParticleLink: " << *isp);
710  regHMPL = new HepMcParticleLink(isp->id(), 0, HepMcParticleLink::IS_POSITION, HepMcParticleLink::IS_ID);
711  }
712  ISF::ISFParticle* regisp = new ISF::ISFParticle(isp->position(),
713  currPar->momentum(),
714  isp->mass(),
715  isp->charge(),
716  isp->pdgCode(),
717  HepMC::status(isp),
718  isp->timeStamp(),
719  *m_isp,
720  HepMC::uniqueID(isp),
721  HepMC::barcode(isp), // FIXME barcode-based
722  regTruthBinding,
723  regHMPL
724  );
725  // add presampled process info
726  if (isp->getUserInformation() && isp->getUserInformation()->materialLimit()) {
727  const ISF::MaterialPathInfo* matLim = isp->getUserInformation()->materialLimit();
728  ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
729  validInfo->setMaterialLimit(
730  matLim->process, matLim->dMax, matLim->process == 121 ? pathLim.l0Collected : pathLim.x0Collected);
731  regisp->setUserInformation(validInfo);
732  }
733  // in the validation mode, add process info
734  if (m_validationMode) {
735  ISF::ParticleUserInformation* validInfo = regisp->getUserInformation();
736  if (!validInfo) {
737  validInfo = new ISF::ParticleUserInformation();
738  regisp->setUserInformation(validInfo);
739  }
740  if (isp->getUserInformation())
741  validInfo->setProcess(isp->getUserInformation()->process());
742  else
743  validInfo->setProcess(-2); // signal problem in the validation chain
744  if (isp->getUserInformation())
745  validInfo->setGeneration(isp->getUserInformation()->generation());
746  else
747  validInfo->setGeneration(-1); // signal problem in the validation chain
748  }
749  // Check that the returned ISFParticle has a valid TruthBinding
750  // before pushing into the particle broker
751  if (!regisp->getTruthBinding()) {
752  ATH_MSG_ERROR("Could not retrieve TruthBinding from non-interacting ISFParticle "<< *regisp);
753  }
754  m_particleBroker->push(regisp, m_isp);
755  }
756 
757  if (isp!=m_isp)
758  {
759  delete isp;
760  return nullptr;
761  }
762 
763  return currPar;
764 }
765 
766 void
767 iFatras::McMaterialEffectsUpdator::ionize(const Trk::TrackParameters& parm,
768  AmgVector(5) & updatedPar,
769  double dInX0,
770  Trk::PropDirection /*dir*/,
771  Trk::ParticleHypothesis particle) const
772 {
773 
774  double p = parm.momentum().mag();
775  double m = Trk::ParticleMasses::mass[particle];
776  double E = sqrt(p*p+m*m);
777 
778  // the following formulas are imported from STEP
779  // preparation of kinetic constants
780 
781  double me = Trk::ParticleMasses::mass[Trk::electron];
782  double mfrac = me/m;
783  double beta = p/E;
784  double gamma = E/m;
785  double Ionization = 0.;
786 
787  //Ionization - Bethe-Bloch
788  double I = 16.e-6 * std::pow(m_matProp->averageZ(),0.9); //16 eV * Z**0.9 - bring to MeV
789 
790  //K/A*Z = 0.5 * 30.7075MeV/(g/mm2) * Z/A * rho[g/mm3] / scale to mm by this
791  double kaz = 0.5*30.7075*m_matProp->zOverAtimesRho();
792 
793  if (particle == Trk::electron){
794  // for electrons use slightly different BetheBloch adaption
795  // see Stampfer, et al, "Track Fitting With Energy Loss", Comp. Pyhs. Comm. 79 (1994), 157-164
796  Ionization = -kaz*(2.*log(2.*me/I)+3.*log(gamma) - 1.95);
797  }
798  else {
799  double eta2 = beta*gamma; eta2 *= eta2;
800  // density effect, only valid for high energies (gamma > 10 -> p > 1GeV for muons)
801  double delta = 0.;
802  if (gamma > 10.) {
803  double eplasma = 28.816e-6 * sqrt(1000.*m_matProp->zOverAtimesRho());
804  delta = 2.*log(eplasma/I) + log(eta2) - 1.;
805  }
806  // tmax - cut off energy
807  double tMax = 2.*eta2*me/(1.+2.*gamma*mfrac+mfrac*mfrac);
808  // divide by beta^2 for non-electrons
809  kaz /= beta*beta;
810  Ionization = -kaz*(log(2.*me*eta2*tMax/(I*I)) - 2.*(beta*beta) - delta);
811  }
812 
813  double energyLoss = Ionization*dInX0*m_matProp->x0();
814  double newP = sqrt(fmax(0.,(E+energyLoss)*(E+energyLoss)-m*m));
815 
816  newP = newP <= 0.5*m_minimumMomentum ? 0.5*m_minimumMomentum : newP; // arbitrary regularization;
817 
818  (updatedPar)[Trk::qOverP] = parm.charge()/newP;
819 
820 }
821 
822 // radiative effects
823 void iFatras::McMaterialEffectsUpdator::radiate(const ISF::ISFParticle* parent, AmgVector(5)& parm ,
824  const Amg::Vector3D& pos, Amg::Vector3D& dir,
825  Trk::TimeLimit timeLim, double pathLim, double& matFraction, double matTot,
826  Trk::PropDirection pdir, Trk::ParticleHypothesis particle) const {
827 
828  // sample energy loss and free path independently
829  double path = 0.;
830  double p = 1./ fabs(parm[Trk::qOverP]);
831 
832 
833  while ( path < pathLim && p>m_minimumMomentum ) {
834 
835  double rndx = CLHEP::RandFlat::shoot(m_randomEngine);
836 
837  // sample visible fraction of the mother momentum taken according to 1/f
838 
839  double eps = fmin(10.,m_minimumMomentum)/p;
840 
841  double z = pow(eps,pow(rndx,exp(1.1))); // adjustment here ?
842 
843  // convert into scaling factor for mother momentum
844  z = (1.- z);
845 
846  // turn into path
847  double dx = -log(z)*0.65; // adjust for mean of exp(-x)
848 
849  // resolve the case when there is not enough material left in the layer
850  if ( path+dx > pathLim ) {
851  double rndp = CLHEP::RandFlat::shoot(m_randomEngine);
852  if (rndp > (pathLim-path)/dx){
853  (parm)[Trk::qOverP] = (parm)[Trk::qOverP]>0 ? 1/p : -1/p;
854  path = pathLim;
855  matFraction = 1.;
856  return; // radiation loop ended
857  }
858  path += dx*rndp;
859  matFraction += dx*rndp/matTot;
860 
861  } else {
862  path+=dx;
863  matFraction += dx/matTot;
864  }
865 
866  /*
867  // sample free path
868  double dx = -log(rndx)/2.; // can be adjusted here
869  //double dx = log(1./rndx)/2.; // can be adjusted here
870  // there may not be enough material to proceed - sample
871  //std::cout <<"brem in layer:"<< path<<"," << dx <<","<< pathLim<< std::endl;
872  if ( path+dx > pathLim ) {
873  double rndp = CLHEP::RandFlat::shoot(m_randomEngine);
874  if (rndp > (pathLim-path)/dx){
875  (parm)[Trk::qOverP] = (parm)[Trk::qOverP]>0 ? 1/p : -1/p;
876  path = pathLim;
877  matFraction = 1.;
878  return; // radiation loop ended
879  }
880  path += dx*rndp;
881  matFraction += dx*rndp/matTot;
882  } else {
883  path+=dx;
884  matFraction += dx/matTot;
885  }
886 
887  // radiate
888  double z = exp(-dx);
889  // resample according to 1/f
890  //double rr = CLHEP::RandFlat::shoot(m_randomEngine);
891  //if (rr*fxx > 1/z ) z=1./rr/fxx ;
892 
893  */
894 
895  // additional adjustement
896  //double rndp = CLHEP::RandFlat::shoot(m_randomEngine);
897  //if ( p*(1-z) > m_minimumBremPhotonMomentum && rndp>0.5 ) {
898  if ( p*(1-z) > m_minimumBremPhotonMomentum ) {
899 
900  double deltaP = (1-z)*p;
901  recordBremPhotonLay(parent,timeLim,p,deltaP, pos,dir,matFraction,pdir,particle);
902  // recordEnergyDeposit(-deltaP,p,particle);
903  p *=z ;
904  }
905 
906  }
907 
908  (parm)[Trk::qOverP] = (parm)[Trk::qOverP] > 0 ? 1/p : -1./p;
909  (parm)[Trk::theta] = dir.theta();
910  (parm)[Trk::phi] = dir.phi();
911 }
912 
913 
914 // actual update mehtod
915 std::unique_ptr<Trk::TrackParameters>
916 iFatras::McMaterialEffectsUpdator::update(
917  double time,
918  const Trk::TrackParameters& parm,
919  const Trk::MaterialProperties& matprop,
920  double pathCorrection,
921  Trk::PropDirection dir,
922  Trk::ParticleHypothesis particle,
923  Trk::MaterialUpdateMode) const
924 {
925  // no material properties - pass them back
926  if (particle==Trk::geantino || (!m_ms && !m_eLoss) ) return nullptr;
927 
928  // get the kinematics
929  double p = parm.momentum().mag();
930  double m = Trk::ParticleMasses::mass[particle];
931  double E = sqrt(p*p+m*m);
932 
933  // Hadronic Interaction section ----------------------------------------------------------------------------------------------
934  // ST add hadronic interaction for secondaries
935  if (particle == Trk::pion && m_hadInt &&
936  m_hadIntProcessor->hadronicInteraction(parm.position(),
937  parm.momentum(),
938  p,
939  E,
940  parm.charge(),
941  matprop,
942  pathCorrection,
943  particle)) {
944  ATH_MSG_VERBOSE(" [+] Hadronic interaction killed initial hadron ... stop "
945  "simulation, tackle childs.");
946  return nullptr;
947  }
948 
949  if (particle == Trk::neutron || particle == Trk::pi0 || particle == Trk::k0) {
950  return parm.uniqueClone();
951  }
952 
953  // the updatedParameters - first a copy
954  AmgVector(5) updatedParameters(parm.parameters());
955 
956  double newP = p;
957  if (m_eLoss){
958  // get the momentum change plus the according sigma
959  Trk::EnergyLoss sampledEnergyLoss( m_eLossUpdator->energyLoss(matprop, p, pathCorrection, dir, particle));
960 
961  double energyLoss = sampledEnergyLoss.deltaE();
962  // protection against NaN
963  if (E + energyLoss < m) {
964  ATH_MSG_VERBOSE(
965  " [+] particle momentum fell under momentum cut - stop simulation");
966  return nullptr;
967  }
968  // smear the mometnum change with the sigma
969  newP = sqrt((E + energyLoss) * (E + energyLoss) - m * m);
970 
971  if (m_recordEnergyDeposition && m_currentSample > 0){
972  if (m_edValidationTree){
973  // fill the variables
974  m_edLayerIntersectX = parm.position().x();
975  m_edLayerIntersectY = parm.position().y();
976  m_edLayerIntersectZ = parm.position().z();
977  m_edLayerIntersectR = parm.position().perp();
978  // the layer and the deposited energy
979  m_edLayerSample = m_currentSample;
980  m_edLayerEnergyDeposit = -1*sampledEnergyLoss.deltaE();
981  // and fill it
982  m_edValidationTree->Fill();
983  }
984  }
985 
986  // the deltaP
987  double deltaP = newP - p;
988 
989  // validation
990  if (m_bremValidation && particle==Trk::electron) m_bremMotherEnergy = p;
991  // record the brem
992  if (particle == Trk::electron && deltaP < -1.*m_minimumBremPhotonMomentum) {
993  Amg::Vector3D recoilDir = parm.momentum().unit(); // to be used for recoil
994  recordBremPhoton(time,p,-deltaP,
995  parm.position(),
996  recoilDir,
997  Trk::electron);
998  }
999  // bail out if the update drops below the momentum cut
1000  if ( newP < m_minimumMomentum) {
1001  ATH_MSG_VERBOSE( " [+] particle momentum fell under momentum cut - stop simulation" );
1002  return nullptr;
1003  }
1004  ATH_MSG_VERBOSE( " [+] delta(P) = " << deltaP );
1005  // continue
1006  // TOM - The random number generator for energy loss is now controlled by MC versions of the energy loss updator
1007  //double sigmaDeltaP = (p+deltaP)*deltaPsigmaQoP.second;
1008  (updatedParameters)[Trk::qOverP] = parm.charge()/(p+deltaP);
1009  // for the validation
1010  m_deltaP = deltaP;
1011  }
1012 
1013  if (m_ms && matprop.x0()>0 ){
1014  // get the projected scattering angle
1015  double sigmaMSproj = (m_use_msUpdator && m_msUpdator) ? sqrt(m_msUpdator->sigmaSquare(matprop, p, pathCorrection, particle)) :
1016  msSigma(pathCorrection*matprop.thicknessInX0(),p,particle);
1017  // do the update
1018  multipleScatteringUpdate(parm,updatedParameters,sigmaMSproj);
1019  }
1020 
1021  // validation mode ---------------------------------------------------------------
1022  if (m_validationMode) {
1023  m_tInX0 = pathCorrection * matprop.thicknessInX0();
1024  m_validationTree->Fill();
1025  }
1026  // -------------------------------------------------------------------------------
1027 
1028  std::optional<AmgSymMatrix(5)> errorMatrix =
1029  (parm.covariance()) ? std::optional<AmgSymMatrix(5)>(*parm.covariance())
1030  : std::nullopt;
1031 
1032  return parm.associatedSurface()
1033  .createUniqueTrackParameters(updatedParameters[0],
1034  updatedParameters[1],
1035  updatedParameters[2],
1036  updatedParameters[3],
1037  updatedParameters[4],
1038  errorMatrix);
1039 }
1040 
1041 std::unique_ptr<Trk::TrackParameters>
1042 iFatras::McMaterialEffectsUpdator::update(
1043  double /*time*/,
1044  const Trk::TrackParameters* parm,
1045  const Trk::MaterialEffectsOnTrack& meff,
1046  Trk::ParticleHypothesis particle,
1047  Trk::MaterialUpdateMode) const
1048 
1049 {
1050 
1051  if (particle != Trk::muon) {
1052  ATH_MSG_WARNING( " [ ---- ] Only implemented for muons yet. No parameterisation for other particles." );
1053  parm->uniqueClone();
1054  }
1055 
1056  // the updatedParameters - first a copy
1057  AmgVector(5) updatedParameters(parm->parameters());
1058 
1059  // get the path correciton --- is 1. for dynamic layers
1060  double pathCorrection = 1.;
1061 
1062  // get the kinematics
1063  double p = parm->momentum().mag();
1064  double m = Trk::ParticleMasses::mass[particle];
1065  double E = sqrt(p*p+m*m);
1066  // the updatedParameters - first a d
1067 
1068 
1069  // MaterialEffectsOnTrack object should contain
1070  // a) multiple scattering part
1071  // b) energy loss part for Landau simulation
1072 
1073  // (a)
1074  if (m_ms){
1075  //If the meff has scattering angles use those, otherwise use MEffUpdator
1076  if(!meff.scatteringAngles()){
1077  //Trick to keep using existing MultipleScatteringUpdator interface
1078  //and create a dummy materialProperties with the properties we are interested in
1079  Trk::MaterialProperties mprop(meff.thicknessInX0(),1.,10.e-10,10.e-10,10.e-10,10.e-10);
1080 
1081  // get the projected scattering angle
1082  double sigmaMSproj = (m_use_msUpdator && m_msUpdator) ? sqrt(m_msUpdator->sigmaSquare(mprop, p, pathCorrection, particle)) :
1083  msSigma(pathCorrection*mprop.thicknessInX0(),p,particle);
1084  // do the update
1085  multipleScatteringUpdate(*parm,updatedParameters,sigmaMSproj);
1086  }
1087 
1088 
1089  }
1090 
1091  // (b)
1092  if (m_eLoss && meff.energyLoss()){
1093  // energy loss update
1094  double energyLossMPV = meff.energyLoss()->deltaE();
1095  double energyLossSigma = meff.energyLoss()->sigmaDeltaE();
1096 
1097  double simulatedEnergyLoss = -1.*(energyLossMPV+energyLossSigma*CLHEP::RandLandau::shoot(m_randomEngine));
1098 
1099  if (m_recordEnergyDeposition){
1100 
1101  ATH_MSG_VERBOSE( " [+] try to record deposited energy (if mapped with a calo sample) " );
1102  const EventContext& ctx = Gaudi::Hive::currentContext(); // FIXME This is slow, but can't change the interface.
1103  const Trk::TrackingGeometry* pTrackingGeometry = trackingGeometry(ctx);
1104  if (pTrackingGeometry) {
1105  const Trk::LayerIndexSampleMap* lism = layerIndexSampleMap();
1106 
1107  if (lism){
1108  // get the Volume and then get the layer
1109  const Trk::TrackingVolume* currentVolume = pTrackingGeometry->lowestTrackingVolume(parm->position());
1110  const Trk::Layer* associatedLayer = currentVolume ? currentVolume->associatedLayer(parm->position()) : nullptr;
1111 
1112  // only go on if you have found an associated Layer && the guy has an index
1113  if (associatedLayer && associatedLayer->layerIndex().value()){
1114  // now try to find it
1115  Trk::LayerIndexSampleMap::const_iterator layerIndexCaloSample = lism->find(associatedLayer->layerIndex());
1116  // layerindex could be found
1117  if (layerIndexCaloSample != lism->end()){
1118  // the validation section
1119  if (m_edValidationTree){
1120  // fill the variables
1121  const Amg::Vector3D& edeposition = parm->position();
1122  m_edLayerIntersectX = edeposition.x();
1123  m_edLayerIntersectY = edeposition.y();
1124  m_edLayerIntersectZ = edeposition.z();
1125  m_edLayerIntersectR = edeposition.perp();
1126  // the layer and the deposited energy
1127  m_edLayerSample = layerIndexCaloSample->second;
1128  m_edLayerEnergyDeposit = -1*simulatedEnergyLoss;
1129  // and fill it
1130  m_edValidationTree->Fill();
1131  } // end of validation section
1132  } // end of calo sample found
1133  } // end of associatedLayer
1134  } // end of lism
1135  } // else -> trackingGeometry exists
1136  } // end of recordEnergyDeposition
1137 
1138  double newP = sqrt((E+simulatedEnergyLoss)*(E+simulatedEnergyLoss)-m*m);
1139 
1140  // set the new momentum
1141  updatedParameters[Trk::qOverP] = parm->charge()/newP;
1142 
1143  }
1144 
1145  // use the manipulator to update the track parameters -------> get ridd of 0!
1146  return parm->associatedSurface()
1147  .createUniqueTrackParameters(updatedParameters[0],
1148  updatedParameters[1],
1149  updatedParameters[2],
1150  updatedParameters[3],
1151  updatedParameters[4],
1152  std::nullopt);
1153 }
1154 
1155 double iFatras::McMaterialEffectsUpdator::msSigma(double dInX0,double p,Trk::ParticleHypothesis particle) const {
1156 
1157  double m = Trk::ParticleMasses::mass[particle];
1158  double E = sqrt(p*p+m*m);
1159  double beta = p/E;
1160 
1161  // PDG Particle Review, Phys.Rev.D86,010001(2012), chapter 30.3, page 328
1162 
1163  double sigmaMSproj = 13.6/beta/p*sqrt(dInX0)*(1.+0.038*log(dInX0));
1164 
1165  return sigmaMSproj;
1166 }
1167 
1168 void iFatras::McMaterialEffectsUpdator::multipleScatteringUpdate(const Trk::TrackParameters& pars,
1169  AmgVector(5)& parameters,
1170  double sigmaMSproj) const
1171 {
1172 
1173  // parametric scattering - independent in x/y
1174  if (m_parametricScattering){
1175  // the initial values
1176  double theta = parameters[Trk::theta];
1177  double phi = parameters[Trk::phi];
1178  double sinTheta = (theta*theta > 10e-10) ? sin(theta) : 1.;
1179 
1180  // sample them in an independent way
1181  double deltaTheta = m_projectionFactor*sigmaMSproj*CLHEP::RandGaussZiggurat::shoot(m_randomEngine);
1182  double deltaPhi = m_projectionFactor*sigmaMSproj*CLHEP::RandGaussZiggurat::shoot(m_randomEngine)/sinTheta;
1183 
1184  phi += deltaPhi;
1185  if (phi >= M_PI) phi -= M_PI;
1186  else if (phi < -M_PI) phi += M_PI;
1187  if (theta > M_PI) theta -= M_PI;
1188 
1189  ATH_MSG_VERBOSE( " [+] deltaPhi / deltaTheta = " << deltaPhi << " / " << deltaTheta );
1190 
1191  // assign the new values
1192  parameters[Trk::phi] = phi;
1193  parameters[Trk::theta] = fabs(theta + deltaTheta);
1194 
1195  // for the validation
1196  m_thetaMStheta = deltaTheta;
1197  m_thetaMSphi = deltaPhi;
1198 
1199  } else {
1200  // scale the projected out of the plane
1201  double thetaMs = m_projectionFactor*sigmaMSproj*CLHEP::RandGaussZiggurat::shoot(m_randomEngine);
1202  double psi = 2.*M_PI*CLHEP::RandFlat::shoot(m_randomEngine);
1203  // more complex but "more true" -
1204  CLHEP::Hep3Vector parsMomHep( pars.momentum().x(), pars.momentum().y(), pars.momentum().z() );
1205  CLHEP::Hep3Vector newDirectionHep( parsMomHep.unit().x(), parsMomHep.unit().y(), parsMomHep.unit().z());
1206  double x = -newDirectionHep.y();
1207  double y = newDirectionHep.x();
1208  double z = 0.;
1209  // if it runs along the z axis - no good ==> take the x axis
1210  if (newDirectionHep.z()*newDirectionHep.z() > 0.999999) {
1211  x = 1.;
1212  y = 0.;
1213  }
1214  // deflector direction
1215  CLHEP::Hep3Vector deflector(x,y,z);
1216  // rotate the new direction for scattering
1217  newDirectionHep.rotate(thetaMs, deflector);
1218  // and arbitrarily in psi
1219  newDirectionHep.rotate(psi, parsMomHep);
1220  // assign the new values
1221  parameters[Trk::phi] = newDirectionHep.phi();
1222  parameters[Trk::theta] = newDirectionHep.theta();
1223 
1224  if (m_validationMode){
1225  m_thetaMStheta = pars.parameters()[Trk::theta] - parameters[Trk::theta];
1226  m_thetaMSphi = pars.parameters()[Trk::phi] - parameters[Trk::phi];
1227  }
1228  }
1229 
1230 }
1231 
1232 
1233 void iFatras::McMaterialEffectsUpdator::recordBremPhoton(double time,
1234  double pElectron,
1235  double gammaE,
1236  const Amg::Vector3D& vertex,
1237  Amg::Vector3D& particleDir,
1238  Trk::ParticleHypothesis particle ) const
1239 {
1240  // get parent particle
1241  // @TODO: replace by Fatras internal bookkeeping
1242  const ISF::ISFParticle *parent = ISF::ParticleClipboard::getInstance().getParticle();
1243  // something is seriously wrong if there is no parent particle
1244  assert(parent);
1245 
1246  // statistics
1247  ++m_recordedBremPhotons;
1248  // ------------------------------------------------------
1249  // simple approach
1250  // (a) simulate theta uniform within the opening angle of the relativistic Hertz dipole
1251  // theta_max = 1/gamma
1252  // (b)Following the Geant4 approximation from L. Urban -> encapsulate that later
1253  // the azimutal angle
1254 
1255  double psi = 2.*M_PI*CLHEP::RandFlat::shoot(m_randomEngine);
1256 
1257  // the start of the equation
1258  double theta = 0.;
1259  double m = Trk::ParticleMasses::mass[particle];
1260  double E = sqrt(pElectron*pElectron + m*m);
1261 
1262  if (m_uniformHertzDipoleAngle) {
1263  // the simplest simulation
1264  theta = m/E * CLHEP::RandFlat::shoot(m_randomEngine);
1265  } else {
1266  // ----->
1267  theta = m/E;
1268  // follow
1269  double a = 0.625; // 5/8
1270 
1271  double r1 = CLHEP::RandFlat::shoot(m_randomEngine);
1272  double r2 = CLHEP::RandFlat::shoot(m_randomEngine);
1273  double r3 = CLHEP::RandFlat::shoot(m_randomEngine);
1274 
1275  double u = -log(r2*r3)/a;
1276 
1277  theta *= (r1 < 0.25 ) ? u : u*m_oneOverThree; // 9./(9.+27) = 0.25
1278  }
1279 
1280  ATH_MSG_VERBOSE( " [ brem ] Simulated angle to electron = " << theta << "." );
1281 
1282  /*
1283  // more complex but "more true"
1284  Amg::Vector3D newDirection(particleDir);
1285  double x = -newDirection.y();
1286  double y = newDirection.x();
1287  double z = 0.;
1288  // if it runs along the z axis - no good ==> take the x axis
1289  if (newDirection.z()*newDirection.z() > 0.999999)
1290  x = 1.;
1291  // deflector direction
1292  Amg::Vector3D deflector(x,y,z);
1293  // rotate the new direction for scattering
1294  newDirection.rotate(theta, deflector);
1295  // and arbitrarily in psi
1296  newDirection.rotate(psi,particleDir);
1297  */
1298  // resample
1299  //double rand = CLHEP::RandFlat::shoot(m_randomEngine);
1300  //double eps = 0.001;
1301  //theta = eps/(1.-rand*(1-eps)) - eps;
1302 
1303  double th = particleDir.theta()-theta;
1304  double ph = particleDir.phi();
1305  if ( th<0.) { th *=-1; ph += M_PI; }
1306  CLHEP::Hep3Vector newDirectionHep( sin(th)*cos(ph),sin(th)*sin(ph),cos(th) );
1307  CLHEP::Hep3Vector particleDirHep( particleDir.x(), particleDir.y(), particleDir.z() );
1308  newDirectionHep.rotate(psi,particleDirHep);
1309  Amg::Vector3D newDirection( newDirectionHep.x(), newDirectionHep.y(), newDirectionHep.z() );
1310 
1311  // recoil / save input momentum for validation
1312  Amg::Vector3D inEl(pElectron*particleDir);
1313  particleDir = (particleDir*pElectron- gammaE*newDirection).unit();
1314 
1315  //std::cout <<"brem opening angle:in:out:"<< cos(theta) <<","<<newDirection*particleDir<< std::endl;
1316 
1317  // -------> create the brem photon <--------------------
1318  const int status = 1 + HepMC::SIM_STATUS_THRESHOLD;
1319  const int id = HepMC::UNDEFINED_ID; // This will be set if the child particle is saved to the GenEvent
1320  ISF::ISFParticle *bremPhoton = new ISF::ISFParticle( vertex,
1321  gammaE*newDirection,
1322  0,
1323  0,
1324  MC::PHOTON,
1325  status,
1326  time,
1327  *parent,
1328  id
1329  );
1330 
1331  // in the validation mode, add process info
1332  if (m_validationMode) {
1333  ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
1334  validInfo->setProcess(m_processCode);
1335  if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
1336  else validInfo->setGeneration(1); // assume parent is a primary
1337  bremPhoton->setUserInformation(validInfo);
1338  }
1339 
1340  // register TruthIncident
1341  ISF::ISFParticleVector children(1, bremPhoton);
1342  ISF::ISFTruthIncident truth( const_cast<ISF::ISFParticle&>(*parent),
1343  children,
1344  m_processCode ,
1345  parent->nextGeoID(),
1346  ISF::fPrimarySurvives );
1347  m_truthRecordSvc->registerTruthIncident( truth);
1348  // At this point we need to update the properties of the
1349  // ISFParticles produced in the interaction
1350  truth.updateParentAfterIncidentProperties();
1351  truth.updateChildParticleProperties();
1352 
1353  // Check that the new/updated ISFParticles have a valid TruthBinding before pushing into the particle broker
1354  if (!parent->getTruthBinding()) {
1355  ATH_MSG_ERROR("Could not retrieve TruthBinding from parent ISFParticle "<< *parent);
1356  }
1357  for (auto *childParticle : children) {
1358  if (!childParticle->getTruthBinding()) {
1359  ATH_MSG_ERROR("Could not retrieve TruthBinding from child ISFParticle "<< *childParticle);
1360  }
1361  m_particleBroker->push(childParticle, parent);
1362  }
1363 
1364  // save info for validation
1365  if (m_validationMode && m_validationTool.isEnabled()) {
1366  Amg::Vector3D* nMom = new Amg::Vector3D((pElectron-gammaE)*particleDir);
1367  m_validationTool->saveISFVertexInfo(3,vertex,*parent,inEl,nMom,children);
1368  delete nMom;
1369  }
1370 
1371 
1372  // if validation is turned on - go for it
1373  if (m_bremValidationTree){
1374  //
1375  ATH_MSG_VERBOSE( " [ brem ] Filling bremsstrahlung validation TTree ..." );
1376  // spatical information
1377  m_bremPointX = float(vertex.x());
1378  m_bremPointY = float(vertex.y());
1379  m_bremPointR = float(vertex.perp());
1380  m_bremPointZ = float(vertex.z());
1381  // photon energyEnergy
1382  m_bremPhotonEnergy = float(gammaE);
1383  m_bremPhotonAngle = float(theta);
1384  // fill the tree
1385  m_bremValidationTree->Fill();
1386  } else {
1387  ATH_MSG_VERBOSE( " [ brem ] No TTree booked ! " );
1388  }
1389 
1390 }
1391 
1392 void iFatras::McMaterialEffectsUpdator::recordBremPhotonLay(const ISF::ISFParticle* parent,
1393  Trk::TimeLimit timeLim,
1394  double pElectron,
1395  double gammaE,
1396  const Amg::Vector3D& vertex,
1397  Amg::Vector3D& particleDir,
1398  double matFraction,
1399  Trk::PropDirection dir,
1400  Trk::ParticleHypothesis particle ) const
1401 {
1402  // statistics
1403  ++m_recordedBremPhotons;
1404  // ------------------------------------------------------
1405  // simple approach
1406  // (a) simulate theta uniform within the opening angle of the relativistic Hertz dipole
1407  // theta_max = 1/gamma
1408  // (b)Following the Geant4 approximation from L. Urban -> encapsulate that later
1409  // the azimutal angle
1410 
1411  double psi = 2.*M_PI*CLHEP::RandFlat::shoot(m_randomEngine);
1412 
1413  // the start of the equation
1414  double theta = 0.;
1415  double m = Trk::ParticleMasses::mass[particle];
1416  double E = sqrt(pElectron*pElectron + m*m);
1417 
1418  if (m_uniformHertzDipoleAngle) {
1419  // the simplest simulation
1420  theta = m/E * CLHEP::RandFlat::shoot(m_randomEngine);
1421  } else {
1422  // ----->
1423  theta = m/E;
1424  // follow
1425  double a = 0.625; // 5/8
1426 
1427  double r1 = CLHEP::RandFlat::shoot(m_randomEngine);
1428  double r2 = CLHEP::RandFlat::shoot(m_randomEngine);
1429  double r3 = CLHEP::RandFlat::shoot(m_randomEngine);
1430 
1431  double u = -log(r2*r3)/a;
1432 
1433  theta *= (r1 < 0.25 ) ? u : u*m_oneOverThree; // 9./(9.+27) = 0.25
1434  }
1435 
1436  ATH_MSG_VERBOSE( " [ brem ] Simulated angle to electron = " << theta << "." );
1437 
1438  // more complex but "more true"
1439  CLHEP::Hep3Vector particleDirHep( particleDir.x(), particleDir.y(), particleDir.z() );
1440  CLHEP::Hep3Vector newDirectionHep( particleDirHep );
1441  double x = -newDirectionHep.y();
1442  double y = newDirectionHep.x();
1443  double z = 0.;
1444  // if it runs along the z axis - no good ==> take the x axis
1445  if (newDirectionHep.z()*newDirectionHep.z() > 0.999999)
1446  x = 1.;
1447  // deflector direction
1448  CLHEP::Hep3Vector deflectorHep(x,y,z);
1449  // rotate the new direction for scattering
1450  newDirectionHep.rotate(theta, deflectorHep);
1451  // and arbitrarily in psi
1452  newDirectionHep.rotate(psi,particleDirHep);
1453 
1454  // recoil
1455  Amg::Vector3D newDirection( newDirectionHep.x(), newDirectionHep.y(), newDirectionHep.z() );
1456  particleDir = (particleDir*pElectron- gammaE*newDirection).unit();
1457 
1458  // -------> create the brem photon <--------------------
1459  const int status = 1 + HepMC::SIM_STATUS_THRESHOLD;
1460  const int id = HepMC::UNDEFINED_ID; // This will be set if the child particle is saved to the GenEvent
1461  ISF::ISFParticle *bremPhoton = new ISF::ISFParticle( vertex,
1462  gammaE*newDirection,
1463  0,
1464  0,
1465  MC::PHOTON,
1466  status,
1467  timeLim.time,
1468  *parent,
1469  id
1470  );
1471 
1472 
1473  // in the validation mode, add process info
1474  if (m_validationMode) {
1475  ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
1476  validInfo->setProcess(m_processCode);
1477  if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
1478  else validInfo->setGeneration(1); // assume parent is a primary track
1479  bremPhoton->setUserInformation(validInfo);
1480  }
1481 
1482  // register TruthIncident
1483  ISF::ISFParticleVector children(1, bremPhoton);
1484  ISF::ISFTruthIncident truth( const_cast<ISF::ISFParticle&>(*parent),
1485  children,
1486  m_processCode ,
1487  parent->nextGeoID(),
1488  ISF::fPrimarySurvives );
1489  m_truthRecordSvc->registerTruthIncident( truth);
1490  // At this point we need to update the properties of the
1491  // ISFParticles produced in the interaction
1492  truth.updateParentAfterIncidentProperties();
1493  truth.updateChildParticleProperties();
1494 
1495  // Check that the new/updated ISFParticles have a valid TruthBinding
1496  if (!parent->getTruthBinding()) {
1497  ATH_MSG_ERROR("Could not retrieve TruthBinding from parent ISFParticle "<< *parent);
1498  }
1499  for (auto *childParticle : children) {
1500  if (!childParticle->getTruthBinding()) {
1501  ATH_MSG_ERROR("Could not retrieve TruthBinding from child ISFParticle "<< *childParticle);
1502  }
1503  }
1504 
1505  // save info for validation
1506  if (m_validationMode && m_validationTool.isEnabled()) {
1507  Amg::Vector3D* nMom = new Amg::Vector3D((pElectron-gammaE)*particleDir);
1508  m_validationTool->saveISFVertexInfo(3,vertex,*parent,pElectron*particleDir,nMom,children);
1509  delete nMom;
1510  }
1511 
1512 
1513  // layer update : don't push into particle stack untill destiny resolved
1514  Trk::PathLimit pLim = m_samplingTool->sampleProcess(m_randomEngine, bremPhoton->momentum().mag(),0.,Trk::photon);
1515 
1516  // material fraction : flip if direction of propagation changed
1517  double ci = m_layer->surfaceRepresentation().normal().dot(particleDir);
1518  double co = m_layer->surfaceRepresentation().normal().dot(newDirection);
1519 
1520  double remMat = ci*co <0 ? (1-matFraction) : matFraction;
1521 
1522  // decide if photon leaves the layer
1523  // amount of material to traverse
1524  double pathCorrection = m_layer->surfaceRepresentation().normal().dot(bremPhoton->momentum()) !=0 ?
1525  fabs(m_layer->surfaceRepresentation().pathCorrection(bremPhoton->position(),bremPhoton->momentum())) : 1.;
1526 
1527  double mTot = m_matProp->thicknessInX0()*pathCorrection*(1.-remMat);
1528 
1529  if (mTot < pLim.x0Max) { // release
1530 
1531  m_particleBroker->push( bremPhoton, parent);
1532 
1533  //std::cout <<"brem photon:"<<bremPhoton<<" registered for extrapolation from "<< vertex << "; presampled pathLim, layer dx:"<< pLim.x0Max<<","<< mTot <<std::endl;
1534 
1535  } else { // interaction within the layer
1536  auto cparm = std::make_unique<Trk::CurvilinearParameters>(bremPhoton->position(),bremPhoton->momentum(),bremPhoton->charge());
1537  std::unique_ptr<Trk::TrackParameters> uPar = updateInLay(bremPhoton,cparm.get(),remMat,timeLim, pLim, dir, Trk::photon);
1538  if (uPar) { ATH_MSG_VERBOSE( "Check this: parameters should be dummy here (brem.photon) " <<","<<uPar->position() ); }
1539  }
1540 
1541  // if validation is turned on - go for it
1542  if (m_bremValidationTree){
1543  //
1544  ATH_MSG_VERBOSE( " [ brem ] Filling bremsstrahlung validation TTree ..." );
1545  // spatical information
1546  m_bremPointX = float(vertex.x());
1547  m_bremPointY = float(vertex.y());
1548  m_bremPointR = float(vertex.perp());
1549  m_bremPointZ = float(vertex.z());
1550  // photon energyEnergy
1551  m_bremPhotonEnergy = float(gammaE);
1552  m_bremPhotonAngle = float(theta);
1553  // fill the tree
1554  m_bremValidationTree->Fill();
1555  } else {
1556  ATH_MSG_VERBOSE( " [ brem ] No TTree booked ! " );
1557  }
1558 
1559 }
1560 
1561 
1562 const Trk::LayerIndexSampleMap* iFatras::McMaterialEffectsUpdator::layerIndexSampleMap() const
1563 {
1564  // no layer index map --- retrieve it
1565  if (!m_layerIndexCaloSampleMap){
1566  // try to retrieve it from detector store
1567  if ((detStore()->retrieve(m_layerIndexCaloSampleMap, m_layerIndexCaloSampleMapName).isFailure())){
1568  ATH_MSG_WARNING( "Could not retrieve '" << m_layerIndexCaloSampleMapName << "' from DetectorStore." );
1569  ATH_MSG_WARNING( " -> switching recording of energy deposit off." );
1570  }
1571 
1572  if (m_layerIndexCaloSampleMap)
1573  ATH_MSG_VERBOSE( "Successfully retrieved '" << m_layerIndexCaloSampleMapName << "' with entries:"
1574  << m_layerIndexCaloSampleMap->size() );
1575  }
1576  return m_layerIndexCaloSampleMap;
1577 }
1578 
1579 std::unique_ptr<Trk::TrackParameters>
1580 iFatras::McMaterialEffectsUpdator::interact(double time,
1581  const Amg::Vector3D& position,
1582  const Amg::Vector3D& momentum,
1583  Trk::ParticleHypothesis particle,
1584  int process,
1585  const Trk::Material* extMatProp) const
1586 {
1587  // Responsible for registering TruthIncidents
1588  if (process == 0)
1589  return nullptr;
1590 
1591  // get parent particle
1592  // @TODO: replace by Fatras internal bookkeeping
1593  const ISF::ISFParticle *parent = ISF::ParticleClipboard::getInstance().getParticle();
1594  // something is seriously wrong if there is no parent particle
1595  assert(parent);
1596 
1597  double mass = Trk::ParticleMasses::mass[particle];
1598  double p = momentum.mag();
1599  /*double E = sqrt( p*p + mass*mass);*/
1600 
1601  if ( process == 201 ) { // decay
1602  // update parent before decay
1603  ISF::ISFParticleVector childVector = m_particleDecayer->decayParticle(*parent,position,momentum,time);
1604 
1605  // register TruthIncident
1606  ISF::ISFTruthIncident truth( const_cast<ISF::ISFParticle&>(*parent),
1607  childVector,
1608  process,
1609  parent->nextGeoID(), // inherits from the parent
1610  ISF::fKillsPrimary );
1611  m_truthRecordSvc->registerTruthIncident( truth);
1612  // At this point we need to update the properties of the
1613  // ISFParticles produced in the interaction
1614  truth.updateChildParticleProperties();
1615 
1616  for (unsigned int i=0; i<childVector.size(); i++) {
1617  // in the validation mode, add process info
1618  if (m_validationMode) {
1619  ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
1620  validInfo->setProcess(process);
1621  if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
1622  else validInfo->setGeneration(1); // assume parent is a primary track
1623  childVector[i]->setUserInformation(validInfo);
1624  }
1625  // register next geo (is current), next flavor can be defined by filter
1626  childVector[i]->setNextGeoID( parent->nextGeoID() );
1627  // feed it the particle broker with parent information
1628  m_particleBroker->push(childVector[i], parent);
1629  }
1630 
1631  // save info for validation
1632  if (m_validationMode && m_validationTool.isEnabled()) {
1633  Amg::Vector3D* nMom = nullptr;
1634  m_validationTool->saveISFVertexInfo(process,position,*parent,momentum,nMom,childVector);
1635  delete nMom;
1636  }
1637 
1638  return nullptr;
1639  }
1640 
1641  if ( process == 5 ) { // positron annihilation
1642 
1643  ISF::ISFParticleVector children(2);
1644 
1645  double fmin = mass/p;
1646 
1647  double fr = 1.-pow(fmin,CLHEP::RandFlat::shoot(m_randomEngine));
1648 
1649  // double cTh = 1.-fmin/(1.-fr)/fr;
1650 
1651  // first implementation: ctH=1
1652  const int status = 1 + HepMC::SIM_STATUS_THRESHOLD;
1653  const int id = HepMC::UNDEFINED_ID; // This will be set if the child particle is saved to the GenEvent
1654 
1655  children[0] = new ISF::ISFParticle( position,
1656  fr*momentum,
1657  0.,
1658  0.,
1659  MC::PHOTON,
1660  status,
1661  time,
1662  *parent,
1663  id
1664  );
1665 
1666  children[1] = new ISF::ISFParticle( position,
1667  (1-fr)*momentum,
1668  0.,
1669  0.,
1670  MC::PHOTON,
1671  status,
1672  time,
1673  *parent,
1674  id
1675  );
1676 
1677  // register TruthIncident
1678  ISF::ISFTruthIncident truth( const_cast<ISF::ISFParticle&>(*parent),
1679  children,
1680  process,
1681  parent->nextGeoID(), // inherits from the parent
1682  ISF::fPrimarySurvives );
1683  m_truthRecordSvc->registerTruthIncident( truth);
1684  // At this point we need to update the properties of the
1685  // ISFParticles produced in the interaction
1686  truth.updateParentAfterIncidentProperties();
1687  truth.updateChildParticleProperties();
1688 
1689  // Check that the new/updated ISFParticles have a valid TruthBinding
1690  if (!parent->getTruthBinding()) {
1691  ATH_MSG_ERROR("Could not retrieve TruthBinding from parent ISFParticle "<< *parent);
1692  }
1693  for (auto *childParticle : children) {
1694  if (!childParticle->getTruthBinding()) {
1695  ATH_MSG_ERROR("Could not retrieve TruthBinding from child ISFParticle "<< *childParticle);
1696  }
1697  }
1698 
1699  // in the validation mode, add process info
1700  if (m_validationMode) {
1701  ISF::ParticleUserInformation* validInfo1 = new ISF::ParticleUserInformation();
1702  validInfo1->setProcess(process);
1703  if (parent->getUserInformation()) validInfo1->setGeneration(parent->getUserInformation()->generation()+1);
1704  else validInfo1->setGeneration(1); // assume parent is a primary track
1705  children[0]->setUserInformation(validInfo1);
1706  ISF::ParticleUserInformation* validInfo2 = new ISF::ParticleUserInformation();
1707  validInfo2->setProcess(process);
1708  if (parent->getUserInformation()) validInfo2->setGeneration(parent->getUserInformation()->generation()+1);
1709  else validInfo2->setGeneration(1); // assume parent is a primary track
1710  children[1]->setUserInformation(validInfo2);
1711  }
1712 
1713  // push child particles onto stack
1714  m_particleBroker->push( children[0], parent);
1715  m_particleBroker->push( children[1], parent);
1716 
1717  // save info for validation
1718  if (m_validationMode && m_validationTool.isEnabled()) {
1719  Amg::Vector3D* nMom = nullptr;
1720  m_validationTool->saveISFVertexInfo(process,position,*parent,momentum,nMom,children);
1721  delete nMom;
1722  }
1723 
1724  // kill the track -----------------------------
1725  return nullptr;
1726  }
1727 
1728  if (process==14) { // photon conversion
1729 
1730  auto parm = std::make_unique<Trk::NeutralCurvilinearParameters>(position,momentum,parent->charge());
1731 
1732  bool cStat = m_conversionTool->doConversion(time, *parm); // Registers TruthIncident internally
1733 
1734  if (!cStat) ATH_MSG_WARNING( "Conversion failed, killing photon anyway ");
1735 
1736  // kill the mother particle
1737  return nullptr;
1738  }
1739 
1740  if (process==121) { // hadronic interaction
1741 
1742  //const Amg::Vector3D pDir = momentum.unit();
1743 
1744  auto parm = std::make_unique<Trk::CurvilinearParameters>(position,momentum,parent->charge());
1745 
1746  bool recHad = m_hadIntProcessor->doHadronicInteraction(time, position, momentum, extMatProp, particle, true); // Registers TruthIncident internally
1747  // eventually : bool recHad = m_hadIntProcessor->recordHadState( time, p, position, pDir, particle);
1748 
1749  // kill the track if interaction recorded --------------------------
1750  if (recHad) {
1751  return nullptr;
1752  } else {
1753  return parm;
1754  }
1755  }
1756 
1757  return nullptr;
1758 }
1759 
1760 ISF::ISFParticleVector iFatras::McMaterialEffectsUpdator::interactLay(const ISF::ISFParticle* parent,
1761  double time,
1762  const Trk::TrackParameters& parm,
1763  Trk::ParticleHypothesis particle,
1764  int process,
1765  const Trk::MaterialProperties* extMatProp) const {
1766  ISF::ISFParticleVector childVector(0);
1767 
1768  if ( process==0 ) return childVector;
1769 
1770  double mass = Trk::ParticleMasses::mass[particle];
1771  double p = parm.momentum().mag();
1772  /*double E = sqrt( p*p + mass*mass);*/
1773  const Amg::Vector3D& position=parm.position();
1774  const Amg::Vector3D& momentum=parm.momentum();
1775 
1776  if ( process == 5 ) { // positron annihilation
1777 
1778  ISF::ISFParticleVector children(2);
1779 
1780  double fmin = mass/p;
1781 
1782  double fr = 1.-pow(fmin,CLHEP::RandFlat::shoot(m_randomEngine));
1783 
1784  // double cTh = 1.-fmin/(1.-fr)/fr;
1785 
1786  // first implementation: ctH=1
1787  const int status = 1 + HepMC::SIM_STATUS_THRESHOLD;
1788  const int id = HepMC::UNDEFINED_ID; // This will be set if the child particle is saved to the GenEvent
1789 
1790  children[0] = new ISF::ISFParticle( position,
1791  fr*momentum,
1792  0.,
1793  0.,
1794  MC::PHOTON,
1795  status,
1796  time,
1797  *parent,
1798  id
1799  );
1800 
1801  children[1] = new ISF::ISFParticle( position,
1802  (1-fr)*momentum,
1803  0.,
1804  0.,
1805  MC::PHOTON,
1806  status,
1807  time,
1808  *parent,
1809  id
1810  );
1811 
1812  // in the validation mode, add process info
1813  if (m_validationMode) {
1814  ISF::ParticleUserInformation* validInfo1 = new ISF::ParticleUserInformation();
1815  validInfo1->setProcess(process);
1816  if (parent->getUserInformation()) validInfo1->setGeneration(parent->getUserInformation()->generation()+1);
1817  else validInfo1->setGeneration(1); // assume parent is a primary track
1818  children[0]->setUserInformation(validInfo1);
1819  ISF::ParticleUserInformation* validInfo2 = new ISF::ParticleUserInformation();
1820  validInfo2->setProcess(process);
1821  if (parent->getUserInformation()) validInfo2->setGeneration(parent->getUserInformation()->generation()+1);
1822  else validInfo2->setGeneration(-1); // assume parent is a primary track
1823  children[1]->setUserInformation(validInfo2);
1824  }
1825 
1826  // register TruthIncident
1827  ISF::ISFTruthIncident truth( const_cast<ISF::ISFParticle&>(*parent),
1828  children,
1829  process,
1830  parent->nextGeoID(), // inherits from the parent
1831  ISF::fPrimarySurvives );
1832  m_truthRecordSvc->registerTruthIncident( truth);
1833  // At this point we need to update the properties of the
1834  // ISFParticles produced in the interaction
1835  truth.updateParentAfterIncidentProperties();
1836  truth.updateChildParticleProperties();
1837 
1838  // save info for validation
1839  if (m_validationMode && m_validationTool.isEnabled()) {
1840  Amg::Vector3D* nMom = nullptr;
1841  m_validationTool->saveISFVertexInfo(process,position,*parent,momentum,nMom,children);
1842  delete nMom;
1843  }
1844 
1845  // Check that the new/updated ISFParticles have a valid TruthBinding
1846  if (!parent->getTruthBinding()) {
1847  ATH_MSG_ERROR("Could not retrieve TruthBinding from parent ISFParticle "<< *parent);
1848  }
1849  for (auto *childParticle : children) {
1850  if (!childParticle->getTruthBinding()) {
1851  ATH_MSG_ERROR("Could not retrieve TruthBinding from child ISFParticle "<< *childParticle);
1852  }
1853  }
1854 
1855  return children;
1856  }
1857 
1858  if (process==14) { // photon conversion
1859 
1860  Trk::NeutralCurvilinearParameters neu(position,momentum,parent->charge());
1861  childVector=m_conversionTool->doConversionOnLayer(parent, time, neu); // Registers TruthIncident internally
1862 
1863  // validation mode
1864  if (m_validationMode && m_validationTool.isEnabled()) {
1865 
1866  // add process info for children
1867  for (unsigned int i=0; i<childVector.size(); i++) {
1868  ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
1869  validInfo->setProcess(process);
1870  if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
1871  else validInfo->setGeneration(1); // assume parent is a primary track
1872  childVector[i]->setUserInformation(validInfo);
1873  }
1874  // save interaction info
1875  Amg::Vector3D* nMom = nullptr;
1876  m_validationTool->saveISFVertexInfo(process, position,*parent,momentum,nMom,childVector);
1877  delete nMom;
1878  }
1879 
1880  // Check that the new ISFParticles have a valid TruthBinding
1881  for (auto *childParticle : childVector) {
1882  if (!childParticle->getTruthBinding()) {
1883  ATH_MSG_ERROR("Could not retrieve TruthBinding from child ISFParticle "<< *childParticle);
1884  }
1885  }
1886 
1887  return childVector;
1888  }
1889 
1890  if (process==121) { // hadronic interaction
1891 
1892  return ( m_hadIntProcessor->doHadIntOnLayer(parent, time, position, momentum,
1893  extMatProp? &extMatProp->material() : nullptr, particle) );
1894 
1895  }
1896 
1897  if ( process == 201 ) { // decay
1898  childVector = m_particleDecayer->decayParticle(*parent,parm.position(),parm.momentum(),time);
1899  // in the validation mode, add process info
1900  if (m_validationMode) {
1901  for (unsigned int i=0; i<childVector.size(); i++) {
1902  ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
1903  validInfo->setProcess(process);
1904  if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
1905  else validInfo->setGeneration(1); // assume parent is a primary track
1906  childVector[i]->setUserInformation(validInfo);
1907  }
1908  }
1909 
1910  // register TruthIncident
1911  ISF::ISFTruthIncident truth( const_cast<ISF::ISFParticle&>(*parent),
1912  childVector,
1913  process,
1914  parent->nextGeoID(), // inherits from the parent
1915  ISF::fKillsPrimary );
1916  m_truthRecordSvc->registerTruthIncident( truth);
1917  // At this point we need to update the properties of the
1918  // ISFParticles produced in the interaction
1919  truth.updateChildParticleProperties();
1920 
1921  // save info for validation
1922  if (m_validationMode && m_validationTool.isEnabled()) {
1923  Amg::Vector3D* nMom = nullptr;
1924  m_validationTool->saveISFVertexInfo(process,parm.position(),*parent,parm.momentum(),nMom,childVector);
1925  delete nMom;
1926  }
1927  }
1928 
1929  // Check that the new ISFParticles have a valid TruthBinding
1930  for (auto *childParticle : childVector) {
1931  if (!childParticle->getTruthBinding()) {
1932  ATH_MSG_ERROR("Could not retrieve TruthBinding from child ISFParticle "<< *childParticle);
1933  }
1934  }
1935 
1936  return childVector;
1937 }
1938 
1939 const Trk::TrackingGeometry*
1940 iFatras::McMaterialEffectsUpdator::trackingGeometry(const EventContext& ctx) const
1941 {
1942  if (m_useConditions) {
1943  SG::ReadCondHandle<Trk::TrackingGeometry> handle(m_trackingGeometryReadKey, ctx);
1944  if (!handle.isValid()) {
1945  ATH_MSG_FATAL(
1946  "Could not retrieve TrackingGeometry from Conditions Store.");
1947  throw iFatras::McMaterialEffectsUpdatorException();
1948  }
1949  return handle.cptr();
1950  } else {
1951  const Trk::TrackingGeometry* trackingGeometry = nullptr;
1952  if (detStore()
1953  ->retrieve(trackingGeometry, m_trackingGeometryName)
1954  .isFailure()) {
1955  ATH_MSG_FATAL("Could not retrieve TrackingGeometry from DetectorStore.");
1956  throw iFatras::McMaterialEffectsUpdatorException();
1957  }
1958  return trackingGeometry;
1959  }
1960 }
1961 
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
iFatras::McMaterialEffectsUpdator::m_processCode
int m_processCode
MCTruth process code for TruthIncidents created by this tool.
Definition: McMaterialEffectsUpdator.h:232
Trk::anyDirection
@ anyDirection
Definition: PropDirection.h:22
make_hlt_rep.pars
pars
Definition: make_hlt_rep.py:90
ISF::ISFTruthIncident::updateChildParticleProperties
void updateChildParticleProperties()
Update the id and particleLink properties of the child particles (to be called after registerTruthInc...
Definition: ISFTruthIncident.cxx:239
iFatras::McMaterialEffectsUpdator::recordBremPhoton
void recordBremPhoton(double time, double pElectron, double gammaE, const Amg::Vector3D &vertex, Amg::Vector3D &particleDir, Trk::ParticleHypothesis particle) const
the helper function for a brem photon record
Definition: McMaterialEffectsUpdator.cxx:1233
EnergyLoss.h
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
Trk::ParticleSwitcher::particle
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:76
iFatras::McMaterialEffectsUpdator::msSigma
double msSigma(double dInX0, double p, Trk::ParticleHypothesis particle) const
handle the Energy loss
Definition: McMaterialEffectsUpdator.cxx:1155
Trk::k0
@ k0
Definition: ParticleHypothesis.h:35
Trk::Surface::associatedDetectorElement
const TrkDetElementBase * associatedDetectorElement() const
return associated Detector Element
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
athena.path
path
python interpreter configuration --------------------------------------—
Definition: athena.py:126
IParticleDecayHelper.h
TrackParameters.h
Trk::PathLimit::process
int process
Definition: HelperStructs.h:39
ParticleGun_SamplingFraction.eta2
eta2
Definition: ParticleGun_SamplingFraction.py:96
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
SG::ReadCondHandle
Definition: ReadCondHandle.h:44
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
iFatras::McMaterialEffectsUpdator::m_samplingTool
ToolHandle< iFatras::IProcessSamplingTool > m_samplingTool
MCTruth process sampling.
Definition: McMaterialEffectsUpdator.h:235
Trk::ParametersBase::charge
double charge() const
Returns the charge.
Surface.h
Trk::LayerMaterialProperties::fullMaterial
virtual const MaterialProperties * fullMaterial(const Amg::Vector3D &gp) const =0
Return method for full material description of the Layer.
iFatras::McMaterialEffectsUpdator::multipleScatteringUpdate
void multipleScatteringUpdate(const Trk::TrackParameters &parm, AmgVector(5) &parameters, double sigmaMSproj) const
the private multiple Scattering update method, thetaMs is the projected random number
Definition: McMaterialEffectsUpdator.cxx:1168
MaterialProperties.h
iFatras::McMaterialEffectsUpdator::~McMaterialEffectsUpdator
virtual ~McMaterialEffectsUpdator()
Destructor.
Trk::ParametersBase::position
const Amg::Vector3D & position() const
Access method for the position.
ISF::ParticleUserInformation::setProcess
void setProcess(int proc)
Definition: ParticleUserInformation.h:90
Trk::PathLimit::updateMat
void updateMat(float dX0, float Z, float dL0)
collected material update
Definition: HelperStructs.h:51
ISF::TruthBinding
Definition: TruthBinding.h:18
Trk::ParametersBase::associatedSurface
virtual const Surface & associatedSurface() const override=0
Access to the Surface associated to the Parameters.
Trk::ParametersBase::uniqueClone
std::unique_ptr< ParametersBase< DIM, T > > uniqueClone() const
clone method for polymorphic deep copy returning unique_ptr; it is not overriden, but uses the existi...
Definition: ParametersBase.h:97
IEnergyLossUpdator.h
Trk::pi0
@ pi0
Definition: ParticleHypothesis.h:34
xAOD::deltaPhi
setSAddress setEtaMS setDirPhiMS setDirZMS setBarrelRadius setEndcapAlpha setEndcapRadius setInterceptInner setEtaMap setEtaBin setIsTgcFailure setDeltaPt deltaPhi
Definition: L2StandAloneMuon_v1.cxx:160
ISF::ParticleClipboard::getInstance
static ParticleClipboard & getInstance()
get the singleton instance
Definition: ParticleClipboard.h:61
iFatras::McMaterialEffectsUpdator::recordBremPhotonLay
void recordBremPhotonLay(const ISF::ISFParticle *parent, Trk::TimeLimit time, double pElectron, double gammaE, const Amg::Vector3D &vertex, Amg::Vector3D &particleDir, double matFraction, Trk::PropDirection dir, Trk::ParticleHypothesis particle) const
the helper function for a brem photon record
Definition: McMaterialEffectsUpdator.cxx:1392
IParticleBroker.h
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:71
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
iFatras::McMaterialEffectsUpdator::m_layerIndexCaloSampleMapName
std::string m_layerIndexCaloSampleMapName
name to record it
Definition: McMaterialEffectsUpdator.h:277
python.SystemOfUnits.MeV
int MeV
Definition: SystemOfUnits.py:154
Trk::MaterialProperties::thicknessInX0
float thicknessInX0() const
Return the radiationlength fraction.
ParticleClipboard.h
ISF::ISFParticle
Definition: ISFParticle.h:42
iFatras::McMaterialEffectsUpdator::trackingGeometry
virtual const Trk::TrackingGeometry * trackingGeometry(const EventContext &ctx) const
Definition: McMaterialEffectsUpdator.cxx:1940
Trk::EnergyLoss::sigmaDeltaE
double sigmaDeltaE() const
returns the symmatric error
M_PI
#define M_PI
Definition: ActiveFraction.h:11
SG::ReadCondHandle::isValid
bool isValid()
Definition: ReadCondHandle.h:206
ISF::MaterialPathInfo::dMax
float dMax
Definition: ParticleUserInformation.h:29
ISFTruthIncident.h
ISF::ISFParticle::pdgCode
int pdgCode() const
PDG value.
Trk::MaterialProperties::thicknessInL0
float thicknessInL0() const
Return the nuclear interaction length fraction.
Layer.h
Trk::MaterialEffectsBase::thicknessInX0
double thicknessInX0() const
returns the actually traversed material .
Trk::MaterialUpdateMode
MaterialUpdateMode
This is a steering enum to force the material update it can be: (1) addNoise (-1) removeNoise Second ...
Definition: MaterialUpdateMode.h:18
Trk::MaterialProperties::x0
float x0() const
Return the radiation length.
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
iFatras::McMaterialEffectsUpdator::m_parametricScattering
bool m_parametricScattering
describe deflection parametric/do real deflection
Definition: McMaterialEffectsUpdator.h:251
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
IPhysicsValidationTool.h
Trk::MaterialProperties::material
const Material & material() const
Return the stored Material.
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
x
#define x
iFatras::McMaterialEffectsUpdator::m_randomEngineName
std::string m_randomEngineName
Name of the random number stream.
Definition: McMaterialEffectsUpdator.h:269
iFatras::McMaterialEffectsUpdator::m_recordEnergyDeposition
bool m_recordEnergyDeposition
for deposition methods
Definition: McMaterialEffectsUpdator.h:276
ParamDefs.h
ISF::ISFParticle::id
int id() const
unique ID
SUSY_SimplifiedModel_PostInclude.process
string process
Definition: SUSY_SimplifiedModel_PostInclude.py:42
Trk::u
@ u
Enums for curvilinear frames.
Definition: ParamDefs.h:83
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition: EventPrimitives.h:52
dqt_zlumi_pandas.mass
mass
Definition: dqt_zlumi_pandas.py:170
Trk::LayerIndexSampleMap
std::map< Trk::LayerIndex, int > LayerIndexSampleMap
Definition: LayerIndexSampleMap.h:29
ISF::ISFParticle::position
const Amg::Vector3D & position() const
The current position of the ISFParticle.
iFatras::McMaterialEffectsUpdator::update
std::unique_ptr< Trk::TrackParameters > update(const Trk::TrackParameters *parm, const Trk::Layer &sf, Trk::TimeLimit &time, Trk::PathLimit &path, Trk::GeometrySignature geoID, Trk::PropDirection dir=Trk::alongMomentum, Trk::ParticleHypothesis particle=Trk::pion) const
Updator interface (full update for a layer): A unique_ptr to Trk::TrackParameters is returned it coul...
Definition: McMaterialEffectsUpdator.cxx:372
IPhotonConversionTool.h
iFatras::McMaterialEffectsUpdator::m_eLoss
bool m_eLoss
IEnergyLossUpdator.
Definition: McMaterialEffectsUpdator.h:221
MaterialEffectsOnTrack.h
Trk::ParticleHypothesis
ParticleHypothesis
Definition: ParticleHypothesis.h:25
Trk::MaterialEffectsOnTrack
represents the full description of deflection and e-loss of a track in material.
Definition: MaterialEffectsOnTrack.h:40
CylinderVolumeBounds.h
iFatras::McMaterialEffectsUpdator::updateInLay
std::unique_ptr< Trk::TrackParameters > updateInLay(const ISF::ISFParticle *isp, const Trk::TrackParameters *parm, double &matFraction, Trk::TimeLimit &time, Trk::PathLimit &path, Trk::PropDirection dir=Trk::alongMomentum, Trk::ParticleHypothesis particle=Trk::pion) const
Definition: McMaterialEffectsUpdator.cxx:416
iFatras::McMaterialEffectsUpdator::m_rndGenSvc
ServiceHandle< IAtRndmGenSvc > m_rndGenSvc
Random Generator service
Definition: McMaterialEffectsUpdator.h:266
Trk::PropDirection
PropDirection
Definition: PropDirection.h:19
ISFParticle.h
Trk::GeometrySignature
GeometrySignature
Definition: GeometrySignature.h:24
ParticleGun_EoverP_Config.mom
mom
Definition: ParticleGun_EoverP_Config.py:63
python.TriggerHandler.th
th
Definition: TriggerHandler.py:296
Trk::PathLimit::x0Collected
float x0Collected
Definition: HelperStructs.h:36
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
ParticleGun_EoverP_Config.momentum
momentum
Definition: ParticleGun_EoverP_Config.py:63
TrigVtx::gamma
@ gamma
Definition: TrigParticleTable.h:26
ISF::ISFTruthIncident
Definition: ISFTruthIncident.h:35
lumiFormat.i
int i
Definition: lumiFormat.py:92
iFatras::McMaterialEffectsUpdator::m_msUpdator
ToolHandle< Trk::IMultipleScatteringUpdator > m_msUpdator
Definition: McMaterialEffectsUpdator.h:226
Trk::Surface::createUniqueTrackParameters
virtual ChargedTrackParametersUniquePtr createUniqueTrackParameters(double l1, double l2, double phi, double theat, double qop, std::optional< AmgSymMatrix(5)> cov=std::nullopt) const =0
Use the Surface as a ParametersBase constructor, from local parameters - charged.
z
#define z
beamspotman.n
n
Definition: beamspotman.py:731
Trk::theta
@ theta
Definition: ParamDefs.h:72
iFatras::McMaterialEffectsUpdator::m_hadInt
bool m_hadInt
hadronic interaction setting
Definition: McMaterialEffectsUpdator.h:238
Trk::TrackingGeometry
Definition: TrackingGeometry.h:67
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
Trk::electron
@ electron
Definition: ParticleHypothesis.h:27
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
HepMC::barcode
int barcode(const T *p)
Definition: Barcode.h:16
AmgVector
AmgVector(4) T2BSTrackFilterTool
Definition: T2BSTrackFilterTool.cxx:114
Trk::LayerMaterialProperties
Definition: LayerMaterialProperties.h:62
CLHEP
STD'S.
Definition: IAtRndmGenSvc.h:19
ISF::ParticleUserInformation
Definition: ParticleUserInformation.h:52
Trk::pion
@ pion
Definition: ParticleHypothesis.h:29
ISF::ISFParticleVector
std::vector< ISF::ISFParticle * > ISFParticleVector
ISFParticle vector.
Definition: ISFParticleContainer.h:26
iFatras::McMaterialEffectsUpdator::McMaterialEffectsUpdator
McMaterialEffectsUpdator(const std::string &, const std::string &, const IInterface *)
AlgTool constructor for McMaterialEffectsUpdator.
Definition: McMaterialEffectsUpdator.cxx:53
Trk::geantino
@ geantino
Definition: ParticleHypothesis.h:26
Trk::EnergyLoss::deltaE
double deltaE() const
returns the
HepMC::uniqueID
int uniqueID(const T &p)
Definition: MagicNumbers.h:113
iFatras::McMaterialEffectsUpdatorException
Exception to be thrown when TrackingGeometry not found.
Definition: McMaterialEffectsUpdator.h:56
test_pyathena.parent
parent
Definition: test_pyathena.py:15
iFatras::McMaterialEffectsUpdator::interact
std::unique_ptr< Trk::TrackParameters > interact(double time, const Amg::Vector3D &position, const Amg::Vector3D &momentum, Trk::ParticleHypothesis particle, int process, const Trk::Material *extMatProp=0) const
Definition: McMaterialEffectsUpdator.cxx:1580
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
iFatras::McMaterialEffectsUpdator::layerIndexSampleMap
const Trk::LayerIndexSampleMap * layerIndexSampleMap() const
return the TrackingGeometry used
Definition: McMaterialEffectsUpdator.cxx:1562
IHadronicInteractionProcessor.h
Trk::ParametersBase
Definition: ParametersBase.h:55
Trk::CurvilinearParametersT
Definition: CurvilinearParametersT.h:48
Trk::TimeLimit
Definition: HelperStructs.h:58
Trk::muon
@ muon
Definition: ParticleHypothesis.h:28
Trk::neutron
@ neutron
Definition: ParticleHypothesis.h:33
Trk::LayerIndex::value
int value() const
layerIndex expressed in an integer
Definition: LayerIndex.h:71
beamspotman.dir
string dir
Definition: beamspotman.py:623
HepMC::UNDEFINED_ID
constexpr int UNDEFINED_ID
Definition: MagicNumbers.h:55
ISF::ParticleUserInformation::setGeneration
void setGeneration(int gen)
Definition: ParticleUserInformation.h:92
iFatras::McMaterialEffectsUpdator::m_ms
bool m_ms
IMultipleScatteringUpdator.
Definition: McMaterialEffectsUpdator.h:225
grepfile.ic
int ic
Definition: grepfile.py:33
iFatras::McMaterialEffectsUpdator::m_minimumMomentum
double m_minimumMomentum
Minimum momentum cut.
Definition: McMaterialEffectsUpdator.h:245
Trk::ParticleMasses::mass
constexpr double mass[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:53
iFatras::McMaterialEffectsUpdator::finalize
StatusCode finalize()
AlgTool finalize method.
Definition: McMaterialEffectsUpdator.cxx:360
iFatras::McMaterialEffectsUpdator::initialize
StatusCode initialize()
AlgTool initailize method.
Definition: McMaterialEffectsUpdator.cxx:168
python.PyKernel.detStore
detStore
Definition: PyKernel.py:41
ParticleUserInformation.h
Trk::TrackingVolume::volumeName
const std::string & volumeName() const
Returns the VolumeName - for debug reason, might be depreciated later.
IProcessSamplingTool.h
IMultipleScatteringUpdator.h
ISF::ISFParticle::timeStamp
double timeStamp() const
Timestamp of the ISFParticle.
HepMC::SIM_STATUS_THRESHOLD
constexpr int SIM_STATUS_THRESHOLD
Constant definiting the status threshold for simulated particles, eg. can be used to separate generat...
Definition: MagicNumbers.h:38
iFatras::McMaterialEffectsUpdator::m_bendingCorrection
bool m_bendingCorrection
Switch to use bending correction.
Definition: McMaterialEffectsUpdator.h:263
iFatras::McMaterialEffectsUpdator::m_minimumBremPhotonMomentum
double m_minimumBremPhotonMomentum
Minimum momentum for brem photons.
Definition: McMaterialEffectsUpdator.h:248
ISF::ParticleUserInformation::setMaterialLimit
void setMaterialLimit(int process, float x0lim, float x0coll)
Definition: ParticleUserInformation.h:94
VP1PartSpect::E
@ E
Definition: VP1PartSpectFlags.h:21
iFatras::McMaterialEffectsUpdator::m_edValidation
bool m_edValidation
Definition: McMaterialEffectsUpdator.h:349
ISF::ISFParticle::momentum
const Amg::Vector3D & momentum() const
The current momentum vector of the ISFParticle.
iFatras::McMaterialEffectsUpdator::m_particleBroker
ServiceHandle< ISF::IParticleBroker > m_particleBroker
Definition: McMaterialEffectsUpdator.h:375
Trk::nonInteracting
@ nonInteracting
Definition: ParticleHypothesis.h:25
Trk::EnergyLoss
This class describes energy loss material effects in the ATLAS tracking EDM.
Definition: EnergyLoss.h:34
charge
double charge(const T &p)
Definition: AtlasPID.h:494
iFatras::McMaterialEffectsUpdator::m_conversionTool
ToolHandle< iFatras::IPhotonConversionTool > m_conversionTool
IPhotonConversionTool.
Definition: McMaterialEffectsUpdator.h:229
ISF::ISFParticle::setUserInformation
void setUserInformation(ParticleUserInformation *userInfo)
Trk::PathLimit
Definition: HelperStructs.h:34
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
ISF::ISFParticle::charge
double charge() const
charge of the particle
iFatras::McMaterialEffectsUpdator::m_truthRecordSvc
ServiceHandle< ISF::ITruthSvc > m_truthRecordSvc
Definition: McMaterialEffectsUpdator.h:376
Trk::PathLimit::l0Collected
float l0Collected
Definition: HelperStructs.h:37
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
ISF::ParticleClipboard::getParticle
const ISF::ISFParticle * getParticle() const
get the particle from the clipboard
Definition: ParticleClipboard.h:72
ISF::ISFParticle::getParticleLink
const HepMcParticleLink * getParticleLink() const
HepMcParticleLink accessors.
Definition: ISFParticle.h:172
ISF::MaterialPathInfo
Definition: ParticleUserInformation.h:28
iFatras::McMaterialEffectsUpdator::m_bremValidation
bool m_bremValidation
Definition: McMaterialEffectsUpdator.h:328
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
iFatras::McMaterialEffectsUpdator::m_validationMode
bool m_validationMode
Definition: McMaterialEffectsUpdator.h:307
Trk::vertex
@ vertex
Definition: MeasurementType.h:21
Trk::MaterialProperties
Definition: MaterialProperties.h:40
TrackingVolume.h
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
ISF::ParticleUserInformation::materialLimit
const MaterialPathInfo * materialLimit() const
Definition: ParticleUserInformation.h:88
a
TList * a
Definition: liststreamerinfos.cxx:10
Trk::MaterialEffectsOnTrack::energyLoss
const EnergyLoss * energyLoss() const
returns the energy loss object.
y
#define y
CaloSwCorrections.time
def time(flags, cells_name, *args, **kw)
Definition: CaloSwCorrections.py:242
iFatras::McMaterialEffectsUpdator::m_use_msUpdator
bool m_use_msUpdator
switch between MSUpdator and local parametrization
Definition: McMaterialEffectsUpdator.h:254
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
ISF::ISFTruthIncident::updateParentAfterIncidentProperties
void updateParentAfterIncidentProperties()
Update the id and particleLink properties of the parentAfterIncident (to be called after registerTrut...
Definition: ISFTruthIncident.cxx:211
iFatras::McMaterialEffectsUpdator::radiate
void radiate(const ISF::ISFParticle *parent, AmgVector(5) &updatedParameters, const Amg::Vector3D &pos, Amg::Vector3D &dir, Trk::TimeLimit time, double dX0, double &matFraction, double matTot, Trk::PropDirection pdir, Trk::ParticleHypothesis particle) const
Definition: McMaterialEffectsUpdator.cxx:823
unit
const PlainObject unit() const
This is a plugin that makes Eigen look like CLHEP & defines some convenience methods.
Definition: AmgMatrixBasePlugin.h:20
makeTRTBarrelCans.dx
tuple dx
Definition: makeTRTBarrelCans.py:20
Trk::MaterialEffectsOnTrack::scatteringAngles
const ScatteringAngles * scatteringAngles() const
returns the MCS-angles object.
iFatras::McMaterialEffectsUpdator::m_eLossUpdator
ToolHandle< Trk::IEnergyLossUpdator > m_eLossUpdator
Definition: McMaterialEffectsUpdator.h:222
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:73
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
McMaterialEffectsUpdator.h
Trk::TrackingGeometry::lowestTrackingVolume
const TrackingVolume * lowestTrackingVolume(const Amg::Vector3D &gp) const
return the lowest tracking Volume
Definition: TrackingGeometry.cxx:53
Trk::TimeLimit::time
float time
Definition: HelperStructs.h:60
python.DecayParser.children
children
Definition: DecayParser.py:32
Trk::photon
@ photon
Definition: ParticleHypothesis.h:32
iFatras::McMaterialEffectsUpdator::ionize
void ionize(const Trk::TrackParameters &parm, AmgVector(5) &updatedPar, double dInX0, Trk::PropDirection pdir, Trk::ParticleHypothesis particle) const
Definition: McMaterialEffectsUpdator.cxx:767
physics_parameters.parameters
parameters
Definition: physics_parameters.py:144
declareProperty
#define declareProperty(n, p, h)
Definition: BaseFakeBkgTool.cxx:15
ISF::ParticleUserInformation::process
int process() const
Definition: ParticleUserInformation.h:84
TrackingGeometry.h
Trk::phi
@ phi
Definition: ParamDefs.h:81
merge.status
status
Definition: merge.py:17
Trk::PathLimit::x0Max
float x0Max
Definition: HelperStructs.h:35
ITruthSvc.h
Trk::Material
Definition: Material.h:116
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:14
I
#define I(x, y, z)
Definition: MD5.cxx:116
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
MuonParameters::beta
@ beta
Definition: MuonParamDefs.h:144
HepMC::status
int status(const T &p)
Definition: MagicNumbers.h:130
ISF::fPrimarySurvives
@ fPrimarySurvives
Definition: ISFTruthIncident.h:24
Trk::Surface
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:75
iFatras::McMaterialEffectsUpdator::interactLay
ISF::ISFParticleVector interactLay(const ISF::ISFParticle *isp, double time, const Trk::TrackParameters &parm, Trk::ParticleHypothesis particle, int process, const Trk::MaterialProperties *extMatProp=0) const
Definition: McMaterialEffectsUpdator.cxx:1760
grepfile.fr
fr
Definition: grepfile.py:32
iFatras::McMaterialEffectsUpdator::m_validationTool
ToolHandle< IPhysicsValidationTool > m_validationTool
Definition: McMaterialEffectsUpdator.h:308
Trk::TrackingVolume
Definition: TrackingVolume.h:121
ISF::ParticleUserInformation::generation
int generation() const
Definition: ParticleUserInformation.h:86
Trk::TrackingVolume::associatedLayer
const Layer * associatedLayer(const Amg::Vector3D &gp) const
Return the associated Layer.
Definition: TrackingVolume.cxx:569
iFatras::McMaterialEffectsUpdator::m_hadIntProcessor
ToolHandle< IHadronicInteractionProcessor > m_hadIntProcessor
Definition: McMaterialEffectsUpdator.h:239
iFatras::McMaterialEffectsUpdator::m_particleDecayer
ToolHandle< IParticleDecayHelper > m_particleDecayer
Particle Decay.
Definition: McMaterialEffectsUpdator.h:242
readCCLHist.float
float
Definition: readCCLHist.py:83
ISF::MaterialPathInfo::process
int process
Definition: ParticleUserInformation.h:31
HepMCHelpers.h
ISF::ISFParticle::getUserInformation
const ParticleUserInformation * getUserInformation() const
get/set ParticleUserInformation
ISF::ISFParticle::getTruthBinding
const TruthBinding * getTruthBinding() const
pointer to the simulation truth - optional, can be 0
iFatras::McMaterialEffectsUpdator::m_referenceMaterial
bool m_referenceMaterial
Switch to use reference material.
Definition: McMaterialEffectsUpdator.h:260
ISF::fKillsPrimary
@ fKillsPrimary
Definition: ISFTruthIncident.h:25
Trk::Layer
Definition: Layer.h:73
SG::ReadCondHandle::cptr
const_pointer_type cptr()
Definition: ReadCondHandle.h:67
Trk::Layer::layerIndex
const LayerIndex & layerIndex() const
get the layerIndex
ISF::ISFParticle::mass
double mass() const
mass of the particle
Generated on Fri Jul 5 2024 21:21:23 for ATLAS Offline Software by   doxygen 1.8.18