ATLAS Offline Software
PhotonConversionTool.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 "PhotonConversionTool.h"
7 
8 // Gaudi Kernel
9 #include "GaudiKernel/RndmGenerators.h"
10 #include "GaudiKernel/DataSvc.h"
11 #include "GaudiKernel/SmartDataPtr.h"
12 #include "StoreGate/StoreGateSvc.h"
13 // ISF includes
14 #include "ISF_Event/ISFParticle.h"
15 #include "ISF_Event/ISFParticleContainer.h"
16 #include "ISF_Event/ParticleClipboard.h"
17 #include "ISF_Event/ParticleUserInformation.h"
18 #include "ISF_Interfaces/IParticleBroker.h"
19 #include "ISF_Interfaces/ITruthSvc.h"
20 #include "ISF_Event/ISFTruthIncident.h"
21 #include "TruthUtils/MagicNumbers.h"
22 // iFatras
23 #include "ISF_FatrasInterfaces/IPhysicsValidationTool.h"
24 // Trk inlcude
25 #include "TrkEventPrimitives/PdgToParticleHypothesis.h"
26 #include "TrkExInterfaces/IEnergyLossUpdator.h"
27 #include "TrkExInterfaces/ITimedExtrapolator.h"
28 #include "TrkExInterfaces/IMultipleScatteringUpdator.h"
29 #include "TrkSurfaces/Surface.h"
30 #include "TrkGeometry/Layer.h"
31 #include "TrkGeometry/MaterialProperties.h"
32 #include "TrkVolumes/CylinderVolumeBounds.h"
33 #include "TrkNeutralParameters/NeutralParameters.h"
34 // CLHEP
35 #include "CLHEP/Units/SystemOfUnits.h"
36 #include "CLHEP/Matrix/Vector.h"
37 #include "CLHEP/Random/RandFlat.h"
38 #include "CLHEP/Vector/LorentzVector.h"
39 // Validation mode - TTree includes
40 #include "TTree.h"
41 #include "GaudiKernel/ITHistSvc.h"
42 // STD
43 #include <cmath>
44 
45 #include "TrkGeometry/TrackingGeometry.h"
46 #include "TrkGeometry/TrackingVolume.h"
47 
48 #include "TrkEventPrimitives/ParticleHypothesis.h"
49 
50 
52 namespace {
53  constexpr double s_alpha = 1./137.;
54  constexpr double s_oneOverThree = 1./3.;
55 }
56 
57 const Trk::PdgToParticleHypothesis iFatras::PhotonConversionTool::s_pdgToHypo;
58 
59 // constructor
60 iFatras::PhotonConversionTool::PhotonConversionTool(const std::string& t, const std::string& n, const IInterface* p) :
61  base_class(t,n,p),
62  m_particleBroker("ISF_ParticleParticleBroker", n),
63  m_truthRecordSvc("ISF_TruthRecordSvc", n),
64  m_processCode(14),
65  m_referenceMaterial(true),
66  m_minChildEnergy(50.*CLHEP::MeV),
67  m_childEnergyScaleFactor(1.),
68  m_conversionProbScaleFactor(0.98),
69  m_rndGenSvc("AtDSFMTGenSvc", n),
70  m_randomEngine(nullptr),
71  m_randomEngineName("FatrasRnd"),
72  m_recordedConversions(0),
73  m_droppedConversions(0),
74  m_validationMode(false),
75  m_validationTreeName("FatrasPhotonConversions"),
76  m_validationTreeDescription("Validation output from the PhotonConversionTool"),
77  m_validationTreeFolder("/val/FatrasPhotonConversions"),
78  m_validationTree(nullptr),
79  m_validationTool(""),
80  m_conversionPointX(0.),
81  m_conversionPointY(0.),
82  m_conversionPointR(0.),
83  m_conversionPointZ(0.),
84  m_conversionPhotonEnergy(0.),
85  m_conversionChildEnergy(0.),
86  m_conversionChildAngle(0.)
87 {
88  // ISF Services and Tools
89  declareProperty("ParticleBroker" , m_particleBroker , "ISF ParticleBroker Svc" );
90  declareProperty("TruthRecordSvc" , m_truthRecordSvc , "ISF Particle Truth Svc" );
91  // MC Truth Properties
92  declareProperty("PhysicsProcessCode" , m_processCode , "MCTruth Physics Process Code");
93  // the steering - general things
94  declareProperty("ReferenceMaterial" , m_referenceMaterial);
95  declareProperty("MinimumChildEnergy" , m_minChildEnergy);
96  // NB: ChildEnergyScaling=1 corresponds to energy sharing formulas of GEANT4
97  declareProperty("ChildEnergyScaling" , m_childEnergyScaleFactor);
98  declareProperty("ConversionProbScaling" , m_conversionProbScaleFactor);
99  declareProperty("RandomNumberService" , m_rndGenSvc , "Random number generator");
100  declareProperty("RandomStreamName" , m_randomEngineName , "Name of the random number stream");
101  // validation --------------------------------------------------------------
102  declareProperty("ValidationMode" , m_validationMode);
103  declareProperty("PhysicsValidationTool" , m_validationTool);
104 }
105 
106 // destructor
107 iFatras::PhotonConversionTool::~PhotonConversionTool()
108 = default;
109 
110 // Athena standard methods
111 // initialize
112 StatusCode iFatras::PhotonConversionTool::initialize()
113 {
114  ATH_MSG_DEBUG( "initialize()" );
115 
116  // ISF Services
117  if (m_particleBroker.retrieve().isFailure()){
118  ATH_MSG_FATAL( "Could not retrieve ISF Particle Broker: " << m_particleBroker );
119  return StatusCode::FAILURE;
120  }
121  if (m_truthRecordSvc.retrieve().isFailure()){
122  ATH_MSG_FATAL( "Could not retrieve ISF Truth Record Service: " << m_truthRecordSvc );
123  return StatusCode::FAILURE;
124  }
125 
126  // the random service
127  if ( m_rndGenSvc.retrieve().isFailure() ){
128  ATH_MSG_FATAL( "Could not retrieve " << m_rndGenSvc );
129  return StatusCode::FAILURE;
130  }
131  //Get own engine with own seeds:
132  m_randomEngine = m_rndGenSvc->GetEngine(m_randomEngineName);
133  if (!m_randomEngine) {
134  ATH_MSG_FATAL( "Could not get random engine '" << m_randomEngineName << "'" );
135  return StatusCode::FAILURE;
136  }
137 
138  // the validation setup ----------------------------------------------------------------------------------
139  ATH_CHECK( m_validationTool.retrieve( DisableTool{ m_validationTool.empty() || !m_validationMode } ) );
140  if (m_validationMode){
141  ATH_MSG_VERBOSE( "Booking conversion validation TTree ... " );
142 
143  // create the new Tree
144  m_validationTree = new TTree(m_validationTreeName.c_str(), m_validationTreeDescription.c_str());
145 
146  // conversion positions and information
147  m_validationTree->Branch("ConversionPositionX" , &m_conversionPointX, "convX/F");
148  m_validationTree->Branch("ConversionPositionY" , &m_conversionPointY, "convY/F");
149  m_validationTree->Branch("ConversionPositionR" , &m_conversionPointR, "convR/F");
150  m_validationTree->Branch("ConversionPositionZ" , &m_conversionPointZ, "convZ/F");
151  m_validationTree->Branch("ConversionPhotonEnergy" , &m_conversionPhotonEnergy, "convPhotonE/F");
152  m_validationTree->Branch("ConversionChildEnergy" , &m_conversionChildEnergy, "convChildE/F");
153  m_validationTree->Branch("ConversionChildAngle " , &m_conversionChildAngle, "convChildA/F");
154 
155  // now register the Tree
156  SmartIF<ITHistSvc> tHistSvc{Gaudi::svcLocator()->service("THistSvc")};
157  if (!tHistSvc) {
158  ATH_MSG_ERROR( "initialize() Could not find Hist Service -> Switching ValidationMode Off !" );
159  delete m_validationTree; m_validationTree = nullptr;
160  }
161  if ((tHistSvc->regTree(m_validationTreeFolder, m_validationTree)).isFailure()) {
162  ATH_MSG_ERROR( "initialize() Could not register the validation Tree -> Switching ValidationMode Off !" );
163  delete m_validationTree; m_validationTree = nullptr;
164  }
165 
166  } // ------------- end of validation mode -----------------------------------------------------------------
167  ATH_MSG_DEBUG( "finalize() successful" );
168 
169  return StatusCode::SUCCESS;
170 }
171 
172 // finalize
173 StatusCode iFatras::PhotonConversionTool::finalize()
174 {
175 
176  ATH_MSG_INFO( " ---------- Statistics output -------------------------- " );
177  ATH_MSG_INFO( " Minimum energy cut of conversions into e+e- : " << m_minChildEnergy << " [MeV] " );
178  ATH_MSG_INFO( " Conversions into e+e- (above cut, recorded) : " << m_recordedConversions );
179  ATH_MSG_INFO( " Conversions into e+e- (below cut, dropped) : " << m_droppedConversions );
180  ATH_MSG_DEBUG( "finalize() successful" );
181 
182  return StatusCode::SUCCESS;
183 }
184 
185 // ----------------- private helper methods ----------------------------------------------------
186 void iFatras::PhotonConversionTool::recordChilds(double time,
187  const Amg::Vector3D& vertex,
188  const Amg::Vector3D& photonDirection,
189  double childEnergy, double photonEnergy,
190  const Amg::Vector3D& childDirection,
191  Trk::ParticleHypothesis childType) const
192 {
193  // get parent particle
194  // @TODO: replace by Fatras internal bookkeeping
195  const ISF::ISFParticle *parent = ISF::ParticleClipboard::getInstance().getParticle();
196  // something is seriously wrong if there is no parent particle
197  assert(parent);
198 
199  // calculate the child momentum
200  double p1 = sqrt(childEnergy*childEnergy-Trk::ParticleMasses::mass[childType]*Trk::ParticleMasses::mass[childType]);
201 
202  // now properly : energy-momentum conservation
203  CLHEP::HepLorentzVector vtmp;
204  CLHEP::Hep3Vector photonDirectionHep( photonDirection.x(), photonDirection.y(), photonDirection.z() );
205  CLHEP::Hep3Vector childDirectionHep( childDirection.x(), childDirection.y(), childDirection.z() );
206  vtmp = CLHEP::HepLorentzVector(photonEnergy*photonDirectionHep, photonEnergy)
207  - CLHEP::HepLorentzVector(p1*childDirectionHep, childEnergy);
208 
209  double p2 = vtmp.vect().mag();
210  //Amg::Vector3D secondChildDir(vtmp.vect().unit());
211 
212  // charge sampling
213  double charge1, charge2;
214  charge1 = charge2 = 0.;
215  if (CLHEP::RandFlat::shoot(m_randomEngine)>0.5) {
216  charge1 = -1.;
217  charge2 = 1.;
218  }
219  else {
220  charge1 = 1.;
221  charge2 = -1.;
222  }
223 
224  // add the new secondary states to the ISF particle stack
225  double mass = Trk::ParticleMasses::mass[childType];
226  int pdg1 = s_pdgToHypo.convert(childType, charge1, false);
227  int pdg2 = s_pdgToHypo.convert(childType, charge2, false);
228 
229  // remove soft children
230  int nchild = 0;
231  if ( p1 > m_minChildEnergy ) nchild++;
232  if ( p2 > m_minChildEnergy ) nchild++;
233 
234  ISF::ISFParticleVector children(nchild);
235 
236  int ichild = 0;
237  const int status = 1 + HepMC::SIM_STATUS_THRESHOLD;
238  const int id = HepMC::UNDEFINED_ID; // This will be set if the child particle is saved to the GenEvent
239  if ( p1 > m_minChildEnergy ) {
240  ISF::ISFParticle* ch1 = new ISF::ISFParticle( vertex,
241  p1*childDirection,
242  mass,
243  charge1,
244  pdg1,
245  status,
246  time,
247  *parent,
248  id
249  );
250  // in the validation mode, add process info
251  if (m_validationMode) {
252  ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
253  validInfo->setProcess(14);
254  if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
255  else validInfo->setGeneration(1); // assume parent is a primary track
256  ch1->setUserInformation(validInfo);
257  }
258  children[ichild] = ch1;
259  ichild++;
260  }
261 
262  if ( p2 > m_minChildEnergy ) {
263  ISF::ISFParticle* ch2 = new ISF::ISFParticle( vertex,
264  p2*childDirection,
265  mass,
266  charge2,
267  pdg2,
268  status,
269  time,
270  *parent,
271  id
272  );
273 
274  // in the validation mode, add process info
275  if (m_validationMode) {
276  ISF::ParticleUserInformation* validInfo = new ISF::ParticleUserInformation();
277  validInfo->setProcess(14);
278  if (parent->getUserInformation()) validInfo->setGeneration(parent->getUserInformation()->generation()+1);
279  else validInfo->setGeneration(1); // assume parent is a primary track
280  ch2->setUserInformation(validInfo);
281  }
282  children[ichild] = ch2;
283  }
284 
285  // register TruthIncident
286  ISF::ISFTruthIncident truth( const_cast<ISF::ISFParticle&>(*parent),
287  children,
288  m_processCode,
289  parent->nextGeoID(),
290  ISF::fKillsPrimary );
291  m_truthRecordSvc->registerTruthIncident( truth);
292  // At this point we need to update the properties of the
293  // ISFParticles produced in the interaction
294  truth.updateChildParticleProperties();
295 
296  // push onto ParticleStack
297  if (!children.empty() ) {
298  for (auto *childParticle : children) {
299  //Check that the new ISFParticles have a valid TruthBinding
300  if (!childParticle->getTruthBinding()) {
301  ATH_MSG_ERROR("Could not retrieve TruthBinding from child ISFParticle "<< *childParticle);
302  }
303  m_particleBroker->push(childParticle, parent);
304  }
305  }
306 
307  // save info for validation
308  if (m_validationMode && m_validationTool.isEnabled()) {
309  Amg::Vector3D* nPrim=nullptr;
310  m_validationTool->saveISFVertexInfo(14,vertex,*parent,parent->momentum(),nPrim,children);
311  }
312 
313 }
314 
315 ISF::ISFParticleVector iFatras::PhotonConversionTool::getChilds(const ISF::ISFParticle* parent,
316  double time,
317  const Amg::Vector3D& vertex,
318  const Amg::Vector3D& photonDirection,
319  double childEnergy, double photonEnergy,
320  const Amg::Vector3D& childDirection,
321  Trk::ParticleHypothesis childType) const
322 {
323  // Called by PhotonConversionTool::doConversionOnLayer
324  // calculate the child momentum
325  double p1 = sqrt(childEnergy*childEnergy-Trk::ParticleMasses::mass[childType]*Trk::ParticleMasses::mass[childType]);
326 
327  // now properly : energy-momentum conservation
328  CLHEP::HepLorentzVector vtmp;
329  CLHEP::Hep3Vector photonDirectionHep( photonDirection.x(), photonDirection.y(), photonDirection.z() );
330  CLHEP::Hep3Vector childDirectionHep( childDirection.x(), childDirection.y(), childDirection.z() );
331  vtmp = CLHEP::HepLorentzVector(photonEnergy*photonDirectionHep, photonEnergy)
332  - CLHEP::HepLorentzVector(p1*childDirectionHep, childEnergy);
333 
334  double p2 = vtmp.vect().mag();
335  //Amg::Vector3D secondChildDir(vtmp.vect().unit());
336 
337  // charge sampling
338  double charge1, charge2;
339  charge1 = charge2 = 0.;
340  if (CLHEP::RandFlat::shoot(m_randomEngine)>0.5) {
341  charge1 = -1.;
342  charge2 = 1.;
343  }
344  else {
345  charge1 = 1.;
346  charge2 = -1.;
347  }
348 
349  double mass = Trk::ParticleMasses::mass[childType];
350  int pdg1 = s_pdgToHypo.convert(childType, charge1, false);
351  int pdg2 = s_pdgToHypo.convert(childType, charge2, false);
352 
353  // removal of soft children to be done in layer mat updator
354  const int status = 1 + HepMC::SIM_STATUS_THRESHOLD;
355  const int id = HepMC::UNDEFINED_ID; // This will be set if the child particle is saved to the GenEvent
356  std::unique_ptr<ISF::ISFParticle> ch1(new ISF::ISFParticle(vertex,
357  p1*childDirection,
358  mass,
359  charge1,
360  pdg1,
361  status,
362  time,
363  *parent,
364  id
365  ));
366 
367  std::unique_ptr<ISF::ISFParticle> ch2(new ISF::ISFParticle(vertex,
368  p2*childDirection,
369  mass,
370  charge2,
371  pdg2,
372  status,
373  time,
374  *parent,
375  id
376  ));
377 
378  ISF::ISFParticleVector children{ch1.release(),
379  ch2.release()};
380 
381  // register TruthIncident
382  ISF::ISFTruthIncident truth( const_cast<ISF::ISFParticle&>(*parent),
383  children,
384  m_processCode,
385  parent->nextGeoID(),
386  ISF::fKillsPrimary );
387  m_truthRecordSvc->registerTruthIncident( truth);
388  // At this point we need to update the properties of the
389  // ISFParticles produced in the interaction
390  truth.updateChildParticleProperties();
391 
392  // Check that the new ISFParticles have a valid TruthBinding
393  for (auto *childParticle : children) {
394  if (!childParticle->getTruthBinding()) {
395  ATH_MSG_ERROR("Could not retrieve TruthBinding from child ISFParticle "<< *childParticle);
396  }
397  }
398 
399  return children;
400 }
401 
402 bool iFatras::PhotonConversionTool::pairProduction(const Trk::MaterialProperties& mprop,
403  double pathCorrection,
404  double p) const
405 {
406 
407  // use for the moment only Al data - Yung Tsai - Rev.Mod.Particle Physics Vol. 74, No.4, October 1974
408  // optainef from a fit given in the momentum range 100 10 6 2 1 0.6 0.4 0.2 0.1 GeV
409 
411  // Double_t fitFunction(Double_t *x, Double_t *par) {
412  // return par[0] + par[1]*pow(x[0],par[2]);
413  // }
414 
415  // EXT PARAMETER STEP FIRST
416  // NO. NAME VALUE ERROR SIZE DERIVATIVE
417  // 1 p0 -7.01612e-03 8.43478e-01 1.62766e-04 1.11914e-05
418  // 2 p1 7.69040e-02 1.00059e+00 8.90718e-05 -8.41167e-07
419  // 3 p2 -6.07682e-01 5.13256e+00 6.07228e-04 -9.44448e-07
420 
421  double p0 = -7.01612e-03;
422  double p1 = 7.69040e-02;
423  double p2 = -6.07682e-01;
424  // calculate xi
425  double xi = p0 + p1*pow(p/1000.,p2);
426  // now calculate what's left
427  double attenuation = exp( -7.777e-01*pathCorrection*mprop.thicknessInX0()*(1.-xi) );
428 
429  ATH_MSG_VERBOSE( "[ conv ] Propability of conversion = " << (1.-attenuation) * 100. << " %." );
430  /* ST do when the actual conversion happens */
431  /*
432  // record for the validation mode
433  if (m_validationMode){
434  // spatical information
435  m_conversionPointX = float(vertex.x());
436  m_conversionPointY = float(vertex.y());
437  m_conversionPointR = float(vertex.perp());
438  m_conversionPointZ = float(vertex.z());
439  // photon energyEnergy
440  m_conversionPhotonEnergy = float(p);
441  }
442  */
443 
444  return m_conversionProbScaleFactor*CLHEP::RandFlat::shoot(m_randomEngine) > attenuation;
445 
446 }
447 
448 
449 double iFatras::PhotonConversionTool::childEnergyFraction(const Trk::MaterialProperties& mprop,
450  double gammaMom) const {
451 
452  // the fraction
453  double epsilon0 = Trk::ParticleMasses::mass[Trk::electron]/gammaMom;
454  // some needed manipolations
455  double Z = mprop.averageZ();
456  double oneOverZpow = 1./pow(Z,s_oneOverThree);
457  double alphaZsquare = (s_alpha*s_alpha*Z*Z);
458  // now f(Z) - from Z and s_alpha
459  double fZ = alphaZsquare*(1./(1.+alphaZsquare)+0.20206-0.0369*alphaZsquare+0.0083*alphaZsquare*alphaZsquare);
460  double FZ = (8./3)*log(Z)+8.*fZ;
461  // delta_max
462  double deltaMax = exp((42.24-FZ)*.1195)-0.952;
463  // delta_min
464  double deltaMin = 4.*epsilon0*136.*oneOverZpow;
465  // the minimum fraction
466  double epsilon1 = 0.5-0.5*sqrt(1.-deltaMin/deltaMax);
467  double epsilonMin = epsilon1 > epsilon0 ? epsilon1 : epsilon0;
468  // calculate phi1 / phi2 - calculate from deltaMin
469  double Phi1 = phi1(deltaMin);
470  double Phi2 = phi2(deltaMin);
471  // then calculate F10/F20
472  double F10 = 3.*Phi1 - Phi2 - FZ;
473  double F20 = 1.5*Phi1 + 0.5*Phi2 - FZ;
474  // and finally calucate N1, N2
475  double N1 = (0.25-epsilonMin+epsilonMin*epsilonMin)*F10;
476  double N2 = 1.5*F20;
477  // ------------ decide which one to take
478  if ( CLHEP::RandFlat::shoot(m_randomEngine) < N1/(N1+N2) ) {
479  // sample from f1,g1 distribution
480  for ( ; ; ){
481  double epsilon = 0.5 - (0.5 - epsilonMin)*pow(CLHEP::RandFlat::shoot(m_randomEngine),s_oneOverThree);
482  // prepare for the rejection check
483  double delta = 136.*epsilon0*oneOverZpow/(epsilon-epsilon*epsilon);
484  double F1 = 3.*phi1(delta)-phi2(delta)-FZ;
485  // reject ? - or redo the exercise
486  if (F1/F10 > CLHEP::RandFlat::shoot(m_randomEngine)) return m_childEnergyScaleFactor*epsilon;
487  }
488  } else {
489  // sample from f2,g2 distribution
490  for ( ; ; ){
491  double epsilon = epsilonMin + (0.5-epsilonMin)*CLHEP::RandFlat::shoot(m_randomEngine);
492  // prepare for the rejection check
493  double delta = 136.*epsilon0*oneOverZpow/(epsilon-epsilon*epsilon);
494  double F2 = 1.5*phi1(delta)-0.5*phi2(delta)-FZ;
495  // reject ? - or redo the exercise
496  if (F2/F20 > CLHEP::RandFlat::shoot(m_randomEngine)) return m_childEnergyScaleFactor*epsilon;
497  }
498  }
499 
500 }
501 
502 
503 
504 double iFatras::PhotonConversionTool::childEnergyFraction(double gammaMom) const {
505 
506  // the fraction
507  double epsilon0 = Trk::ParticleMasses::mass[Trk::electron]/gammaMom;
508  // some needed manipolations
509  //double Z = mprop.averageZ();
510  double Z = 13.;
511  double oneOverZpow = 1./pow(Z,s_oneOverThree);
512  double alphaZsquare = (s_alpha*s_alpha*Z*Z);
513  // now f(Z) - from Z and s_alpha
514  double fZ = alphaZsquare*(1./(1.+alphaZsquare)+0.20206-0.0369*alphaZsquare+0.0083*alphaZsquare*alphaZsquare);
515  double FZ = (8./3)*log(Z)+8.*fZ;
516  // delta_max
517  double deltaMax = exp((42.24-FZ)*.1195)-0.952;
518  // delta_min
519  double deltaMin = 4.*epsilon0*136.*oneOverZpow;
520  // the minimum fraction
521  double epsilon1 = 0.5-0.5*sqrt(1.-deltaMin/deltaMax);
522  double epsilonMin = epsilon1 > epsilon0 ? epsilon1 : epsilon0;
523  // calculate phi1 / phi2 - calculate from deltaMin
524  double Phi1 = phi1(deltaMin);
525  double Phi2 = phi2(deltaMin);
526  // then calculate F10/F20
527  double F10 = 3.*Phi1 - Phi2 - FZ;
528  double F20 = 1.5*Phi1 + 0.5*Phi2 - FZ;
529  // and finally calucate N1, N2
530  double N1 = (0.25-epsilonMin+epsilonMin*epsilonMin)*F10;
531  double N2 = 1.5*F20;
532  // ------------ decide which one to take
533  if ( CLHEP::RandFlat::shoot(m_randomEngine) < N1/(N1+N2)) {
534  // sample from f1,g1 distribution
535  for ( ; ; ){
536  double epsilon = 0.5 - (0.5 - epsilonMin)*pow(CLHEP::RandFlat::shoot(m_randomEngine),s_oneOverThree);
537  // prepare for the rejection check
538  double delta = 136.*epsilon0*oneOverZpow/(epsilon-epsilon*epsilon);
539  double F1 = 3.*phi1(delta)-phi2(delta)-FZ;
540  // reject ? - or redo the exercise
541  if (F1/F10 > CLHEP::RandFlat::shoot(m_randomEngine)) return m_childEnergyScaleFactor*epsilon;
542  }
543  } else {
544  // sample from f2,g2 distribution
545  for ( ; ; ){
546  double epsilon = epsilonMin + (0.5-epsilonMin)*CLHEP::RandFlat::shoot(m_randomEngine);
547  // prepare for the rejection check
548  double delta = 136.*epsilon0*oneOverZpow/(epsilon-epsilon*epsilon);
549  double F2 = 1.5*phi1(delta)-0.5*phi2(delta)-FZ;
550  // reject ? - or redo the exercise
551  if (F2/F20 > CLHEP::RandFlat::shoot(m_randomEngine)) return m_childEnergyScaleFactor*epsilon;
552  }
553  }
554 
555 }
556 
557 
558 Amg::Vector3D iFatras::PhotonConversionTool::childDirection(const Amg::Vector3D& gammaMom,
559  double childE) const
560 {
561  // --------------------------------------------------
562  // Following the Geant4 approximation from L. Urban
563  // the azimutal angle
564  double psi = 2.*M_PI*CLHEP::RandFlat::shoot(m_randomEngine);
565 
566  // the start of the equation
567  double theta = Trk::ParticleMasses::mass[Trk::electron]/childE;
568  // follow
569  double a = 0.625; // 5/8
570  //double d = 27.;
571 
572  double r1 = CLHEP::RandFlat::shoot(m_randomEngine);
573  double r2 = CLHEP::RandFlat::shoot(m_randomEngine);
574  double r3 = CLHEP::RandFlat::shoot(m_randomEngine);
575 
576  double u = -log(r2*r3)/a;
577 
578  theta *= (r1 < 0.25 ) ? u : u*s_oneOverThree; // 9./(9.+27) = 0.25
579 
580  ATH_MSG_VERBOSE( "[ conv ] Simulated angle to photon = " << theta << "." );
581 
582  // more complex but "more true"
583  CLHEP::Hep3Vector gammaMomHep( gammaMom.x(), gammaMom.y(), gammaMom.z() );
584  CLHEP::Hep3Vector newDirectionHep(gammaMomHep.unit());
585  double x = -newDirectionHep.y();
586  double y = newDirectionHep.x();
587  double z = 0.;
588  // if it runs along the z axis - no good ==> take the x axis
589  if (newDirectionHep.z()*newDirectionHep.z() > 0.999999)
590  x = 1.;
591  // deflector direction
592  CLHEP::Hep3Vector deflectorHep(x,y,z);
593  // rotate the new direction for scattering
594  newDirectionHep.rotate(theta, deflectorHep);
595  // and arbitrarily in psi
596  newDirectionHep.rotate(psi,gammaMomHep);
597 
598  // record for the validation mode
599  if (m_validationMode){
600  // child energy and direction
601  m_conversionChildEnergy = float(childE);
602  m_conversionChildAngle = float(theta);
603  }
604 
605  // assign the new values
606  Amg::Vector3D newDirection( newDirectionHep.x(), newDirectionHep.y(), newDirectionHep.z() );
607  return newDirection;
608 
609 }
610 
611 bool iFatras::PhotonConversionTool::doConversion(double time, const Trk::NeutralParameters& parm,
612  const Trk::ExtendedMaterialProperties* /*extMatProp*/) const {
613 
614  // Called by McMaterialEffectsUpdator::interact
615  double p = parm.momentum().mag();
616 
617  // get the energy
618  double efrac = childEnergyFraction(p);
619  double childEnergy = p*efrac;
620 
621  // count the conversion
622  ++m_recordedConversions;
623 
624  // now get the deflection
625  Amg::Vector3D childDir(childDirection(parm.momentum(), childEnergy));
626  // verbose output
627  ATH_MSG_VERBOSE( "[ conv ] Child energy simulated as : " << childEnergy << " MeV" );
628  // calculate the second child direction
629  recordChilds(time,
630  parm.position(),
631  parm.momentum().unit(),
632  childEnergy, p,
633  childDir,
634  Trk::electron); // Registers TruthIncident internally
635  // fill the TTree ----------------------------
636  if (m_validationTree)
637  m_validationTree->Fill();
638 
639  return true;
640 }
641 
643 ISF::ISFParticleVector iFatras::PhotonConversionTool::doConversionOnLayer(const ISF::ISFParticle* parent,
644  double time, const Trk::NeutralParameters& parm,
645  const Trk::ExtendedMaterialProperties* /*ematprop*/) const {
646  // Called by McMaterialEffectsUpdator::interactLay
647  double p = parm.momentum().mag();
648 
649  // get the energy
650  double efrac = childEnergyFraction(p);
651  double childEnergy = p*efrac;
652 
653  // count the conversion
654  ++m_recordedConversions;
655 
656  // now get the deflection
657  Amg::Vector3D childDir(childDirection(parm.momentum(), childEnergy));
658  // verbose output
659  ATH_MSG_VERBOSE( "[ conv ] Child energy simulated as : " << childEnergy << " MeV" );
660  // fill the TTree ----------------------------
661  if (m_validationTree)
662  m_validationTree->Fill();
663  // calculate the second child direction and return
664  return getChilds(parent,
665  time,
666  parm.position(),
667  parm.momentum().unit(),
668  childEnergy, p,
669  childDir,
670  Trk::electron); // Registers TruthIncident internally
671 
672 }
673 
iFatras::PhotonConversionTool::childEnergyFraction
double childEnergyFraction(const Trk::MaterialProperties &mprop, double gammaMom) const
simulate the child energy
Definition: PhotonConversionTool.cxx:449
iFatras::PhotonConversionTool::doConversionOnLayer
ISF::ISFParticleVector doConversionOnLayer(const ISF::ISFParticle *parent, double time, const Trk::NeutralParameters &parm, const Trk::ExtendedMaterialProperties *ematprop=0) const
interface for processing of the presampled nuclear interactions on layer
Definition: PhotonConversionTool.cxx:643
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
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::PhotonConversionTool::doConversion
bool doConversion(double time, const Trk::NeutralParameters &parm, const Trk::ExtendedMaterialProperties *extMatProp=0) const
interface for processing of the presampled pair production
Definition: PhotonConversionTool.cxx:611
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
iFatras::PhotonConversionTool::getChilds
ISF::ISFParticleVector getChilds(const ISF::ISFParticle *parent, double time, const Amg::Vector3D &vertex, const Amg::Vector3D &photonDirection, double childEnergy, double photonEnergy, const Amg::Vector3D &childDirection, Trk::ParticleHypothesis childType) const
get childs - for layer update
Definition: PhotonConversionTool.cxx:315
dqt_zlumi_pandas.N2
N2
Definition: dqt_zlumi_pandas.py:325
iFatras::PhotonConversionTool::finalize
StatusCode finalize()
AlgTool finalize method.
Definition: PhotonConversionTool.cxx:173
ITimedExtrapolator.h
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
Surface.h
Base_Fragment.mass
mass
Definition: Sherpa_i/share/common/Base_Fragment.py:59
MaterialProperties.h
Trk::ParametersBase::position
const Amg::Vector3D & position() const
Access method for the position.
ISF::ParticleUserInformation::setProcess
void setProcess(int proc)
Definition: ParticleUserInformation.h:90
IEnergyLossUpdator.h
ISF::ParticleClipboard::getInstance
static ParticleClipboard & getInstance()
get the singleton instance
Definition: ParticleClipboard.h:61
iFatras::PhotonConversionTool::m_processCode
int m_processCode
MCTruth process code for TruthIncidents created by this tool.
Definition: PhotonConversionTool.h:114
IParticleBroker.h
python.SystemOfUnits.MeV
int MeV
Definition: SystemOfUnits.py:154
Trk::MaterialProperties::thicknessInX0
float thicknessInX0() const
Return the radiationlength fraction.
TRTCalib_cfilter.p1
p1
Definition: TRTCalib_cfilter.py:130
ParticleClipboard.h
ISF::ISFParticle
Definition: ISFParticle.h:42
M_PI
#define M_PI
Definition: ActiveFraction.h:11
ISFTruthIncident.h
Layer.h
iFatras::PhotonConversionTool::recordChilds
void recordChilds(double time, const Amg::Vector3D &vertex, const Amg::Vector3D &photonDirection, double childEnergy, double photonEnergy, const Amg::Vector3D &childDirection, Trk::ParticleHypothesis childType=Trk::electron) const
record childs - create interface already for mu-/mu+ pair production
Definition: PhotonConversionTool.cxx:186
NeutralParameters.h
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
iFatras::PhotonConversionTool::m_childEnergyScaleFactor
double m_childEnergyScaleFactor
Definition: PhotonConversionTool.h:121
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
IPhysicsValidationTool.h
iFatras::PhotonConversionTool::m_referenceMaterial
bool m_referenceMaterial
Switch to use reference material.
Definition: PhotonConversionTool.h:117
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
x
#define x
PdgToParticleHypothesis.h
Trk::u
@ u
Enums for curvilinear frames.
Definition: ParamDefs.h:77
ISFParticleContainer.h
Trk::ParticleHypothesis
ParticleHypothesis
Definition: ParticleHypothesis.h:25
CylinderVolumeBounds.h
TRTCalib_cfilter.p2
p2
Definition: TRTCalib_cfilter.py:131
ISFParticle.h
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
python.PhysicalConstants.epsilon0
int epsilon0
Definition: PhysicalConstants.py:97
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
ISF::ISFTruthIncident
Definition: ISFTruthIncident.h:35
z
#define z
beamspotman.n
n
Definition: beamspotman.py:731
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
CLHEP
STD'S.
Definition: IAtRndmGenSvc.h:19
ISF::ParticleUserInformation
Definition: ParticleUserInformation.h:52
iFatras::PhotonConversionTool::m_rndGenSvc
ServiceHandle< IAtRndmGenSvc > m_rndGenSvc
Random Generator service
Definition: PhotonConversionTool.h:126
ISF::ISFParticleVector
std::vector< ISF::ISFParticle * > ISFParticleVector
ISFParticle vector.
Definition: ISFParticleContainer.h:26
iFatras::PhotonConversionTool::initialize
StatusCode initialize()
AlgTool initailize method.
Definition: PhotonConversionTool.cxx:112
test_pyathena.parent
parent
Definition: test_pyathena.py:15
dqt_zlumi_pandas.N1
int N1
Definition: dqt_zlumi_pandas.py:322
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
iFatras::PhotonConversionTool::m_randomEngineName
std::string m_randomEngineName
Name of the random number stream.
Definition: PhotonConversionTool.h:128
Trk::ParametersBase
Definition: ParametersBase.h:55
ParticleHypothesis.h
iFatras::PhotonConversionTool::childDirection
Amg::Vector3D childDirection(const Amg::Vector3D &gammaMom, double childE) const
simulate the one child direction - the second one is given clear then
Definition: PhotonConversionTool.cxx:558
HepMC::UNDEFINED_ID
constexpr int UNDEFINED_ID
Definition: MagicNumbers.h:56
ISF::ParticleUserInformation::setGeneration
void setGeneration(int gen)
Definition: ParticleUserInformation.h:92
Trk::ParticleMasses::mass
constexpr double mass[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:53
ParticleUserInformation.h
IMultipleScatteringUpdator.h
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:46
MagicNumbers.h
Trk::PdgToParticleHypothesis
Definition: PdgToParticleHypothesis.h:29
iFatras::PhotonConversionTool::~PhotonConversionTool
virtual ~PhotonConversionTool()
Destructor.
ISF::ISFParticle::setUserInformation
void setUserInformation(ParticleUserInformation *userInfo)
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
PhotonConversionTool.h
iFatras::PhotonConversionTool::m_conversionProbScaleFactor
double m_conversionProbScaleFactor
Definition: PhotonConversionTool.h:122
ISF::ParticleClipboard::getParticle
const ISF::ISFParticle * getParticle() const
get the particle from the clipboard
Definition: ParticleClipboard.h:72
iFatras::PhotonConversionTool::pairProduction
bool pairProduction(const Trk::MaterialProperties &mprop, double pathCorrection, double p) const
interface for processing of the pair production
Definition: PhotonConversionTool.cxx:402
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
Trk::vertex
@ vertex
Definition: MeasurementType.h:21
iFatras::PhotonConversionTool::s_pdgToHypo
static const Trk::PdgToParticleHypothesis s_pdgToHypo
Definition: PhotonConversionTool.h:153
Trk::MaterialProperties
Definition: MaterialProperties.h:40
TrackingVolume.h
iFatras::PhotonConversionTool::m_minChildEnergy
double m_minChildEnergy
The cut from which on the child products are followed.
Definition: PhotonConversionTool.h:120
a
TList * a
Definition: liststreamerinfos.cxx:10
y
#define y
Trk::MaterialProperties::averageZ
float averageZ() const
Returns the average Z of the material.
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
python.DecayParser.children
children
Definition: DecayParser.py:32
iFatras::PhotonConversionTool::m_validationMode
bool m_validationMode
Definition: PhotonConversionTool.h:135
TrackingGeometry.h
merge.status
status
Definition: merge.py:17
ITruthSvc.h
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
iFatras::PhotonConversionTool::m_truthRecordSvc
ServiceHandle< ISF::ITruthSvc > m_truthRecordSvc
Definition: PhotonConversionTool.h:111
iFatras::PhotonConversionTool::m_particleBroker
ServiceHandle< ISF::IParticleBroker > m_particleBroker
ISF services & Tools.
Definition: PhotonConversionTool.h:110
StoreGateSvc.h
readCCLHist.float
float
Definition: readCCLHist.py:83
iFatras::PhotonConversionTool::m_validationTool
ToolHandle< IPhysicsValidationTool > m_validationTool
Definition: PhotonConversionTool.h:141
TRTCalib_cfilter.p0
p0
Definition: TRTCalib_cfilter.py:129
ISF::fKillsPrimary
@ fKillsPrimary
Definition: ISFTruthIncident.h:25
iFatras::PhotonConversionTool::PhotonConversionTool
PhotonConversionTool(const std::string &, const std::string &, const IInterface *)
AlgTool constructor for PhotonConversionTool.
Definition: PhotonConversionTool.cxx:60
Generated on Sun Dec 22 2024 21:16:10 for ATLAS Offline Software by   doxygen 1.8.18