ATLAS Offline Software
Public Member Functions | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
PunchThroughG4Tool Class Reference

Punch through calculations. More...

#include <PunchThroughG4Tool.h>

Inheritance diagram for PunchThroughG4Tool:
Collaboration diagram for PunchThroughG4Tool:

Public Member Functions

 PunchThroughG4Tool (const std::string &, const std::string &, const IInterface *)
 
virtual ~PunchThroughG4Tool ()=default
 
virtual StatusCode initialize () override
 AlgTool initialize method. More...
 
virtual StatusCode finalize () override
 AlgTool finalize method. More...
 
StatusCode initializePhysics () override
 
virtual std::vector< std::map< std::string, double > > computePunchThroughParticles (const G4FastTrack &fastTrack, CLHEP::HepRandomEngine *rndmEngine, double punchThroughProbability, double punchThroughClassifierRand) override
 interface function: fill a vector with the punch-through particles More...
 
virtual void createAllSecondaryTracks (G4ParticleTable &ptable, G4FastStep &fastStep, const G4Track &g4PrimaryTrack, std::vector< std::map< std::string, double >> &secKinematicsMapVect, G4TrackVector &secTrackCont, const std::vector< double > &caloMSVars) override
 create all secondary tracks from kinematics map More...
 
virtual std::vector< double > getCaloMSVars () override
 

Private Member Functions

void checkParticleTable (G4ParticleTable &ptable, int secondarySignedPDG)
 
StatusCode registerPunchThroughParticle (G4ParticleTable &ptable, int pdg, bool doAntiparticle=false, double minEnergy=0., int maxNumParticles=-1, double numParticlesFactor=1., double energyFactor=1., double posAngleFactor=1., double momAngleFactor=1.)
 registers a type of punch-through particles which will be simulated More...
 
StatusCode initializeRegisterPunchThroughParticles ()
 initialize register all the punch-through particles which will be simulated More...
 
StatusCode registerCorrelation (int pdgID1, int pdgID2, double minCorrEnergy=0., double fullCorrEnergy=0.)
 register a correlation for the two given types of punch-through particles with a given energy threshold above which we will have full correlation More...
 
StatusCode initializeRegisterCorrelations ()
 initialize register all correlations between particles More...
 
std::unique_ptr< PunchThroughPDFCreatorreadLookuptablePDF (int pdgID, TFile *fileLookupTable, const std::string &folderName)
 reads out the lookuptable for the given type of particle More...
 
StatusCode checkCaloMSBoundaries (const std::vector< std::pair< double, double >> *rzMS, const std::vector< std::pair< double, double >> *rzCalo)
 Check calo-MS boundaries. More...
 
int getAllParticles (const G4Track &g4PrimaryTrack, std::vector< std::map< std::string, double >> &secKinematicsMapVect, CLHEP::HepRandomEngine *rndmEngine, int pdg, double interpEnergy, double interpEta, int numParticles=-1)
 create the right number of punch-through particles for the given pdg and return the number of particles which was created. More...
 
G4ThreeVector punchTroughPosPropagator (double theta, double phi, double R1, double R2, double z1, double z2) const
 get particle through the calorimeter More...
 
std::vector< std::map< std::string, double > > checkEnergySumFromSecondaries (double mainEnergyInit, std::vector< std::map< std::string, double >> &secKinematicsMapVect)
 check the energies satisfying energy condition More...
 
G4Track * createSecondaryTrack (G4ParticleTable &ptable, G4FastStep &fastStep, double currentTime, int secondarySignedPDG, double energy, double theta, double phi, double momTheta, double momPhi, const std::vector< double > &caloMSVars)
 create secondary track for each given the kinematics More...
 
int getCorrelatedParticles (const G4Track &g4PrimaryTrack, std::vector< std::map< std::string, double >> &secKinematicsMapVect, int pdg, int corrParticles, CLHEP::HepRandomEngine *rndmEngine, double interpEnergy, double interpEta)
 get the right number of particles for the given pdg while considering the correlation to an other particle type, which has already created 'corrParticles' number of particles More...
 
std::map< std::string, double > getOneParticleKinematics (CLHEP::HepRandomEngine *rndmEngine, int secondaryPDG, float initParticleTheta, float initParticlePhi, double interpEnergy, double interpEta) const
 create exactly one punch-through particle with the given pdg and the given max energy More...
 
double getFloatAfterPatternInStr (const char *str, const char *pattern)
 get the floating point number in a string, after the given pattern More...
 
std::vector< double > inversePCA (int pcaCdfIterator, std::vector< double > &variables) const
 
double interpolateEnergy (const double &energy, CLHEP::HepRandomEngine *rndmEngine) const
 
double interpolateEta (const double &eta, CLHEP::HepRandomEngine *rndmEngine) const
 
std::vector< std::map< std::string, std::string > > getInfoMap (const std::string &mainNode, const std::string &xmlFilePath)
 
int passedParamIterator (int pid, double eta, const std::vector< std::map< std::string, std::string >> &mapvect) const
 
StatusCode initializeInverseCDF (const std::string &quantileTransformerConfigFile)
 
StatusCode initializeInversePCA (const std::string &inversePCAConfigFile)
 

Static Private Member Functions

static double inverseCdfTransform (double variable, const std::map< double, double > &inverse_cdf_map)
 
static std::vector< double > dotProduct (const std::vector< std::vector< double >> &m, const std::vector< double > &v)
 
static double normal_cdf (double x)
 
static std::map< double, double > getVariableCDFmappings (xmlNodePtr &nodeParent)
 

Private Attributes

std::vector< double > m_energyPoints
 energy and eta points in param More...
 
std::vector< double > m_etaPoints
 
double m_R1 {0.}
 calo-MS borders More...
 
double m_R2 {0.}
 
double m_z1 {0.}
 
double m_z2 {0.}
 
TFile * m_fileLookupTable {nullptr}
 ROOT objects. More...
 
std::map< int, PunchThroughParticle * > m_particles
 needed to initially create punch-through particles with the right distributions More...
 
StringProperty m_filenameLookupTable {this, "FilenameLookupTable", "CaloPunchThroughParametrisation.root", "holds the filename of the lookup table"}
 
StringProperty m_filenameInverseCDF {this, "FilenameInverseCdf", "", "holds the filename of inverse quantile transformer config"}
 
StringProperty m_filenameInversePCA {this, "FilenameInversePca", "", "holds the filename of inverse PCA config"}
 
IntegerArrayProperty m_pdgInitiators {this, "PunchThroughInitiators", {}, "vector of punch-through initiator pgds"}
 
IntegerArrayProperty m_initiatorsMinEnergy {this, "InitiatorsMinEnergy", {}, "vector of punch-through initiator min energies to create punch through"}
 
DoubleArrayProperty m_initiatorsEtaRange {this, "InitiatorsEtaRange", {}, "vector of min and max abs eta range to allow punch through initiators"}
 
IntegerArrayProperty m_punchThroughParticles {this, "PunchThroughParticles", {}, "vector of pdgs of the particles produced in punch-throughs"}
 
BooleanArrayProperty m_doAntiParticles {this, "DoAntiParticles", {}, "vector of bools to determine if anti-particles are created for each punch-through particle type"}
 
IntegerArrayProperty m_correlatedParticle {this, "CorrelatedParticle", {}, "holds the pdg of the correlated particle for each given pdg"}
 
DoubleArrayProperty m_minCorrEnergy {this, "MinCorrelationEnergy", {}, "holds the energy threshold below which no particle correlation is computed"}
 
DoubleArrayProperty m_fullCorrEnergy {this, "FullCorrelationEnergy", {}, "holds the energy threshold above which a particle correlation is fully developed"}
 
DoubleArrayProperty m_posAngleFactor {this, "ScalePosDeflectionAngles", {}, "tuning parameter to scale the position deflection angles"}
 
DoubleArrayProperty m_momAngleFactor {this, "ScaleMomDeflectionAngles", {}, "tuning parameter to scale the momentum deflection angles"}
 
DoubleArrayProperty m_minEnergy {this, "MinEnergy", {}, "punch-through particles minimum energies"}
 
IntegerArrayProperty m_maxNumParticles {this, "MaxNumParticles", {}, "maximum number of punch-through particles for each particle type"}
 
DoubleArrayProperty m_numParticlesFactor {this, "NumParticlesFactor", {}, "scale the number of punch-through particles"}
 
DoubleArrayProperty m_energyFactor {this, "EnergyFactor", {}, "scale the energy of the punch-through particles"}
 
ServiceHandle< ISF::IGeoIDSvcm_geoIDSvc {this, "GeoIDSvc", "ISF::GeoIDSvc"}
 
ServiceHandle< IEnvelopeDefSvcm_envDefSvc {this, "EnvelopeDefSvc", "AtlasGeometry_EnvelopeDefSvc"}
 
DoubleProperty m_beamPipe {this, "BeamPipeRadius", 500.}
 beam pipe radius More...
 
std::vector< std::vector< std::vector< double > > > m_inverse_PCA_matrix
 pca vectors More...
 
std::vector< std::vector< double > > m_PCA_means
 
std::vector< std::map< std::string, std::string > > m_xml_info_pca
 infoMaps More...
 
std::vector< std::map< std::string, std::string > > m_xml_info_cdf
 
std::vector< std::map< double, double > > m_variable0_inverse_cdf
 (vector of map) for CDF mappings More...
 
std::vector< std::map< double, double > > m_variable1_inverse_cdf
 
std::vector< std::map< double, double > > m_variable2_inverse_cdf
 
std::vector< std::map< double, double > > m_variable3_inverse_cdf
 
std::vector< std::map< double, double > > m_variable4_inverse_cdf
 

Detailed Description

Punch through calculations.

Author
Elmar Ritsch Elmar.nosp@m..Rit.nosp@m.sch@c.nosp@m.ern..nosp@m.ch @maintainer/updater Thomas Carter thoma.nosp@m.s.mi.nosp@m.chael.nosp@m..car.nosp@m.ter@c.nosp@m.ern..nosp@m.ch @maintainer/updater Firdaus Soberi firda.nosp@m.us.s.nosp@m.oberi.nosp@m.@cer.nosp@m.n.ch

Definition at line 50 of file PunchThroughG4Tool.h.

Constructor & Destructor Documentation

◆ PunchThroughG4Tool()

PunchThroughG4Tool::PunchThroughG4Tool ( const std::string &  type,
const std::string &  name,
const IInterface *  parent 
)

Definition at line 39 of file PunchThroughG4Tool.cxx.

40  : base_class(type,name,parent)
41 {
42 }

◆ ~PunchThroughG4Tool()

virtual PunchThroughG4Tool::~PunchThroughG4Tool ( )
virtualdefault

Member Function Documentation

◆ checkCaloMSBoundaries()

StatusCode PunchThroughG4Tool::checkCaloMSBoundaries ( const std::vector< std::pair< double, double >> *  rzMS,
const std::vector< std::pair< double, double >> *  rzCalo 
)
private

Check calo-MS boundaries.

Definition at line 205 of file PunchThroughG4Tool.cxx.

207 {
208  bool found1, found2;
209  found1=false; found2=false;
210 
211  // procedure below is to find intersection of calo and MS boundary (right at the edge of calo and into MS)
212  // at least 2 distinct calo-MS boundary points must be found to proceed
213  for ( const auto & [r_tempCalo , z_tempCalo] : *rzCalo )
214  {
215  // calo point must be more than beampipe radius
216  if (r_tempCalo> m_beamPipe)
217  {
218  for ( const auto & [r_tempMS , z_tempMS] : *rzMS )
219  {
220  // first point, if R and Z of calo matches R and Z of MS
221  if (r_tempCalo==r_tempMS && z_tempCalo==z_tempMS && found1==false )
222  {
223  m_R1=r_tempMS;
224  m_z1=std::fabs(z_tempMS);
225  found1=true;
226  continue;
227  }
228  // second point, if R and Z of calo matches R and Z of MS
229  else if (r_tempCalo==r_tempMS && z_tempCalo==z_tempMS && r_tempCalo!=m_R1 && std::fabs(z_tempCalo)!=m_z1)
230  {
231  m_R2=r_tempMS;
232  m_z2=std::fabs(z_tempMS);
233  found2=true;
234  }
235  }
236 
237  if (found1==true && found2==true) break;
238  }
239  }
240 
241  //in case geometry description changes
242  if (found1 == false){
243  ATH_MSG_ERROR ("[PunchThroughG4Tool] first coordinate of calo-MS border not found");
244  return StatusCode::FAILURE;
245  }
246  if (found2 == false){
247  ATH_MSG_ERROR ("[PunchThroughG4Tool] second coordinate of calo-MS border not found; first one is: R1 ="<<m_R1<<" z1 ="<<m_z1);
248  return StatusCode::FAILURE;
249  }
250 
251  //now order the found values
252  if (m_R1>m_R2) { std::swap(m_R1, m_R2); } //m_R1 - smaller one
253  if (m_z1<m_z2) { std::swap(m_z1, m_z2); } //m_z1 - bigger one
254 
255  if (m_R1==m_R2 || m_z1==m_z2){
256  ATH_MSG_ERROR ("[PunchThroughG4Tool] Bug in propagation calculation! R1="<<m_R1<<" R2 = "<<m_R2<<" z1="<<m_z1<<" z2= "<<m_z2 );
257  return StatusCode::FAILURE;
258  }
259 
260  ATH_MSG_DEBUG ("[PunchThroughG4Tool] calo-MS boundary coordinates: R1="<<m_R1<<" R2 = "<<m_R2<<" z1="<<m_z1<<" z2= "<<m_z2);
261 
262  return StatusCode::SUCCESS;
263 }

◆ checkEnergySumFromSecondaries()

std::vector< std::map< std::string, double > > PunchThroughG4Tool::checkEnergySumFromSecondaries ( double  mainEnergyInit,
std::vector< std::map< std::string, double >> &  secKinematicsMapVect 
)
private

check the energies satisfying energy condition

Definition at line 712 of file PunchThroughG4Tool.cxx.

712  {
713  //initialize variables
714  double energy;
715  double totEnergySecondaries = 0;
716 
717  // Check energy conservation
718  // Step 1: sum all of the energies from the vector of secondaries (map) together, to initialize procedure later
719  for (std::size_t i = 0; i < secKinematicsMapVect.size(); i++) {
720  energy = secKinematicsMapVect[i].at("energy");
721  totEnergySecondaries += energy;
722  }
723 
724  // Step 2: Sort the vector based on "energy" in descending order using lambda expression
725  std::sort(
726  secKinematicsMapVect.begin(),
727  secKinematicsMapVect.end(),
728  [](const std::map<std::string, double>& lhs, const std::map<std::string, double>& rhs) {
729  return lhs.at("energy") > rhs.at("energy");
730  }
731  );
732 
733  // Just for printing here
734  if(totEnergySecondaries > mainEnergyInit){
735  ATH_MSG_DEBUG("[PunchThroughG4Tool::checkEnergySumFromSecondaries] Case where energy of created secondaries more than parent track identified! ");
736  ATH_MSG_DEBUG("[PunchThroughG4Tool::checkEnergySumFromSecondaries] ==> TotalSecondariesEnergy = " << totEnergySecondaries << ", ParentEnergy = "<< mainEnergyInit);
737  }
738 
739  // Step 3: Adjust the energy and remove secondary with lowest energy to ensure total energy of secondaries <= mainEnergyInit
740  while (totEnergySecondaries > mainEnergyInit && !secKinematicsMapVect.empty()) {
741  // Get the last map in the vector (which has the lowest energy due to sorting)
742  std::map<std::string, double> lastMap = secKinematicsMapVect.back();
743 
744  // Get the energy value of the last map
745  double lastEnergy = lastMap.at("energy");
746 
747  // Remove the last map from the vector
748  secKinematicsMapVect.pop_back();
749 
750  // Subtract the removed energy from totEnergySecondaries
751  totEnergySecondaries -= lastEnergy;
752  }
753 
754  //finally return the secKinematicsMapVect after the checks
755  return secKinematicsMapVect;
756 }

◆ checkParticleTable()

void PunchThroughG4Tool::checkParticleTable ( G4ParticleTable &  ptable,
int  secondarySignedPDG 
)
private

Definition at line 277 of file PunchThroughG4Tool.cxx.

277  {
278  // Get Geant4 particle definition
279  G4ParticleDefinition* secG4Particle = ptable.FindParticle(secondarySignedPDG);
280  double mass = secG4Particle->GetPDGMass();
281  ATH_MSG_DEBUG("[PunchThroughG4Tool::checkParticleTable] secondarySignedPDG = " << secondarySignedPDG << ", mass = "<< mass);
282 }

◆ computePunchThroughParticles()

std::vector< std::map< std::string, double > > PunchThroughG4Tool::computePunchThroughParticles ( const G4FastTrack &  fastTrack,
CLHEP::HepRandomEngine *  rndmEngine,
double  punchThroughProbability,
double  punchThroughClassifierRand 
)
overridevirtual

interface function: fill a vector with the punch-through particles

Definition at line 284 of file PunchThroughG4Tool.cxx.

285 {
286  ATH_MSG_DEBUG( "[PunchThroughG4Tool] starting punch-through simulation");
287 
288  // holder to all secondary particle kinematics
289  std::vector<std::map<std::string, double>> secKinematicsMapVect;
290 
291  // Get primary Geant4 track (main particle producing punchthrough)
292  const G4Track * g4PrimaryTrack = fastTrack.GetPrimaryTrack();
293 
294  // Get Geant4 particle definition
295  const G4ParticleDefinition * mainG4Particle = g4PrimaryTrack->GetDefinition();
296 
297  // Get primary Geant4 particle properties
298  int pdgID = mainG4Particle->GetPDGEncoding();
299  float mainPartMass = mainG4Particle->GetPDGMass();
300  float mainMomMag2 = g4PrimaryTrack->GetMomentum().mag2();
301  float mainPartEta = g4PrimaryTrack->GetPosition().eta();
302  // get primary track's current position and momentum direction
303  const G4ThreeVector mainMomentumDir = g4PrimaryTrack->GetMomentumDirection();
304  const G4ThreeVector mainMomentum = g4PrimaryTrack->GetMomentum();
305  const G4ThreeVector mainPosition = g4PrimaryTrack->GetPosition();
306 
307  // calculate incoming energy and eta
308  const double initEnergy = std::sqrt( mainMomMag2 + mainPartMass*mainPartMass );
309  const double initEta = mainPartEta;
310 
311  ATH_MSG_DEBUG("[PunchThroughG4Tool] position of the input particle: r"<<mainPosition.perp()<<" z= "<<mainPosition.z() );
312 
313  // check if it points to the calorimeter - if not, don't simulate
314  // first convert to Amg::Vector3D
315  const Amg::Vector3D amg3DPosVec(mainPosition.x(), mainPosition.y(), mainPosition.z());
316  const Amg::Vector3D amg3DMomDirVec(mainMomentumDir.x(), mainMomentumDir.y(), mainMomentumDir.z());
317  const Amg::Vector3D amg3DMomVec(mainMomentum.x(), mainMomentum.y(), mainMomentum.z());
318  // find nextGeoID to get AtlasDetDescr::AtlasRegion
319  AtlasDetDescr::AtlasRegion nextGeoID = m_geoIDSvc->identifyNextGeoID(amg3DPosVec,amg3DMomVec);
320  if ( nextGeoID != AtlasDetDescr::fAtlasCalo)
321  {
322  ATH_MSG_DEBUG("[PunchThroughG4Tool](GeoIDSvc) input particle doesn't point to calorimeter"<< "Next GeoID: "<< nextGeoID );
323  //return nullptr;
324  return {};
325  }
326 
327  // check if the particle's pdg is registered as a punch-through-causing type
328  {
329  std::vector<int>::const_iterator pdgIt = m_pdgInitiators.begin();
330  std::vector<int>::const_iterator pdgItEnd = m_pdgInitiators.end();
331 
332  std::vector<int>::const_iterator minEnergyIt = m_initiatorsMinEnergy.begin();
333  // loop over all known punch-through initiators
334  for ( ; pdgIt != pdgItEnd; ++pdgIt, ++minEnergyIt)
335  {
336  if (std::abs(pdgID) == *pdgIt){
337  if(initEnergy < *minEnergyIt){
338  ATH_MSG_DEBUG("[PunchThroughG4Tool] particle does not meet initiator min energy requirement. Dropping it in the calo.");
339  //return nullptr;
340  return {};
341  }
342  break;
343  }
344  }
345 
346  // particle will not cause punch-through -> bail out
347  if (pdgIt == pdgItEnd)
348  {
349  ATH_MSG_DEBUG("[PunchThroughG4Tool] particle is not registered as punch-through initiator. Dropping it in the calo.");
350  //return nullptr;
351  return {};
352  }
353  }
354 
355  if(initEta < m_initiatorsEtaRange.value().at(0) || initEta > m_initiatorsEtaRange.value().at(1) ){
356  ATH_MSG_DEBUG("[PunchThroughG4Tool] particle does not meet initiator eta range requirement. Dropping it in the calo.");
357  //return nullptr;
358  return {};
359  }
360 
361  ATH_MSG_DEBUG("[PunchThroughG4Tool] punchThroughProbability output: " << punchThroughProbability << " RandFlat: " << punchThroughClassifierRand );
362 
363  //Safety check: If probability < random number then don't simulate punch through
364  if( punchThroughClassifierRand > punchThroughProbability){
365  ATH_MSG_DEBUG("[PunchThroughG4Tool] particle not classified to create punch through. Dropping it in the calo.");
366  //return nullptr;
367  return {};
368  }
369 
370  //if initial particle is on ID surface, points to the calorimeter, is a punch-through initiator, meets initiator min enery and eta range
371 
372  // this is the place where the magic is done:
373  // test for each registered punch-through pdg if a punch-through
374  // occures and create these particles
375  // -> therefore loop over all registered pdg ids
376  // to keep track of the correlated particles which were already simulated:
377  // first int is pdg, second int is number of particles created
378 
379  // interpolate energy and eta
380  const double interpEnergy = interpolateEnergy(initEnergy, rndmEngine);
381  const double interpEta = interpolateEta(initEta, rndmEngine);
382 
383  std::map<int, int> corrPdgNumDone;
384 
385  int maxTries = 10;
386  int nTries = 0;
387 
388  // loop over all particle pdgs
389  while(secKinematicsMapVect.empty() && nTries < maxTries) { //ensure we always create at least one punch through particle, maxTries to catch very rare cases
390  // loop over all particle pdgs
391  for (const auto& currentParticle : m_particles)
392  {
393  // the pdg that is currently treated
394  int doPdg = currentParticle.first;
395  // get the current particle's correlated pdg
396  int corrPdg = currentParticle.second->getCorrelatedPdg();
397 
398  // if there is a correlated particle type to this one
399  if (corrPdg)
400  {
401  // find out if the current pdg was already simulated
402  std::map<int,int>::iterator itPdgPos = corrPdgNumDone.find(doPdg);
403  // if the current pdg was not simulated yet, find out if
404  // it's correlated one was simulated
405  if ( itPdgPos == corrPdgNumDone.end() ) itPdgPos = corrPdgNumDone.find(corrPdg);
406 
407  // neither this nor the correlated particle type was simulated
408  // so far:
409  if ( itPdgPos == corrPdgNumDone.end() )
410  {
411  // -> roll a dice if we create this particle or its correlated one
412  if ( CLHEP::RandFlat::shoot(rndmEngine) > 0.5 ) doPdg = corrPdg;
413  // now create the particles with the given pdg and note how many
414  // particles of this pdg are created
415  corrPdgNumDone[doPdg] = getAllParticles(*g4PrimaryTrack, secKinematicsMapVect, rndmEngine, doPdg, interpEnergy, interpEta);
416  }
417 
418  // one of the two correlated particle types was already simulated
419  // 'itPdgPos' points to the already simulated one
420  else
421  {
422  // get the pdg of the already simulated particle and the number
423  // of these particles that were created
424  const int donePdg = itPdgPos->first;
425  const int doneNumPart = itPdgPos->second;
426  // set the pdg of the particle type that will be done
427  if (donePdg == doPdg) doPdg = corrPdg;
428 
429  // now create the correlated particles
430  getCorrelatedParticles(*g4PrimaryTrack, secKinematicsMapVect, doPdg, doneNumPart, rndmEngine, interpEnergy, interpEta);
431  // note: no need to take note, that this particle type is now simulated,
432  // since this is the second of two correlated particles, which is
433  // simulated and we do not have correlations of more than two particles.
434  }
435 
436  // if no correlation for this particle
437  // -> directly create all particles with the current pdg
438  }
439  else{
440  getAllParticles(*g4PrimaryTrack, secKinematicsMapVect, rndmEngine, doPdg, interpEnergy, interpEta);
441  }
442  } // for-loop over all particle pdgs
443  }
444 
445  return secKinematicsMapVect;
446 }

◆ createAllSecondaryTracks()

void PunchThroughG4Tool::createAllSecondaryTracks ( G4ParticleTable &  ptable,
G4FastStep &  fastStep,
const G4Track &  g4PrimaryTrack,
std::vector< std::map< std::string, double >> &  secKinematicsMapVect,
G4TrackVector &  secTrackCont,
const std::vector< double > &  caloMSVars 
)
overridevirtual

create all secondary tracks from kinematics map

Definition at line 758 of file PunchThroughG4Tool.cxx.

758  {
759  // Get Geant4 particle definition
760  const G4ParticleDefinition* mainG4Particle = g4PrimaryTrack.GetDefinition();
761  float mainMomMag2 = g4PrimaryTrack.GetMomentum().mag2();
762  float mainPartMass = mainG4Particle->GetPDGMass();
763  double mainEnergyInit = std::sqrt(mainMomMag2 + mainPartMass * mainPartMass);
764 
765  // Energy conservation: call checkEnergySumFromSecondaries to check/update secKinematicsMapVect
766  secKinematicsMapVect = checkEnergySumFromSecondaries(mainEnergyInit, secKinematicsMapVect);
767 
768  // set number of secondary tracks to be produced
769  int numSecondaries = secKinematicsMapVect.size();
770 
771  if(numSecondaries>0){ //safety
772  // set number of secondary tracks from numSecondaries above
773  fastStep.SetNumberOfSecondaryTracks(numSecondaries);
774 
775  // Get current (global) time from original main G4Track
776  double currentTime = g4PrimaryTrack.GetGlobalTime();
777 
778  // create the secondaries given the number of secondaries to be created
779  for(int i=0;i < numSecondaries; i++){
780  int anti = (int)(secKinematicsMapVect[i].at("anti"));
781  int secondaryPDG = (int)(secKinematicsMapVect[i].at("secondaryPDG"));
782  int signedPDG = secondaryPDG*anti;
783  double energy = secKinematicsMapVect[i].at("energy");
784  double theta = secKinematicsMapVect[i].at("theta");
785  double phi = secKinematicsMapVect[i].at("phi");
786  double momTheta = secKinematicsMapVect[i].at("momTheta");
787  double momPhi = secKinematicsMapVect[i].at("momPhi");
788 
789  // (**) finally create the punch-through particle as a G4Track (secondary particle track)
790  ATH_MSG_DEBUG("[PunchThroughG4Tool::createAllSecondaryTracks] createSecondaryTrack input parameters: currentTime = " << currentTime << " anti? = "<< anti <<" signedPDG = "<< signedPDG <<" energy = "<< energy <<" theta = "<< theta <<" phi = "<< phi <<" momTheta = "<< momTheta << " momPhi " << momPhi);
791  G4Track* newSecTrack = createSecondaryTrack( ptable, fastStep, currentTime, signedPDG, energy, theta, phi, momTheta, momPhi, caloMSVars);
792 
793  // if something went wrong and the secondary track is simply not created
794  if (!newSecTrack)
795  {
796  G4Exception("[PunchThroughG4Tool::createAllSecondaryTracks]", "ExceptionError", FatalException, "something went wrong while creating punch-through particle tracks");
797  return;
798  }
799 
800  // add this secondary track to the g4trackvector
801  secTrackCont.push_back( newSecTrack );
802  }
803  }
804 
805  //printout: fastStep.DumpInfo()
806 }

◆ createSecondaryTrack()

G4Track * PunchThroughG4Tool::createSecondaryTrack ( G4ParticleTable &  ptable,
G4FastStep &  fastStep,
double  currentTime,
int  secondarySignedPDG,
double  energy,
double  theta,
double  phi,
double  momTheta,
double  momPhi,
const std::vector< double > &  caloMSVars 
)
private

create secondary track for each given the kinematics

Definition at line 808 of file PunchThroughG4Tool.cxx.

810 {
811  //initialize to nullptr
812  G4Track *newSecTrack = nullptr;
813 
814  G4ParticleDefinition* secG4Particle = ptable.FindParticle(secondarySignedPDG);
815  double mass = secG4Particle->GetPDGMass();
816 
817  // the intersection point with Calo-MS surface
818  G4ThreeVector posVec = punchTroughPosPropagator(theta, phi, caloMSVars[0], caloMSVars[1], caloMSVars[2], caloMSVars[3]);
819 
820  // set up the real punch-through particle at this position
821  // set up the momentum vector of this particle as a GlobalMomentum
822  // by using the given energy and mass of the particle and also using
823  // the given theta and phi
824  G4ThreeVector momVec;
825  double momMag = std::sqrt(energy*energy - mass*mass);
826  momVec.setRThetaPhi(momMag,momTheta,momPhi);
827  ATH_MSG_DEBUG("[PunchThroughG4Tool]::createSecondaryTrack] setRThetaPhi pre input parameters: energy = " << energy <<" mass = "<< mass);
828  ATH_MSG_DEBUG("[PunchThroughG4Tool]::createSecondaryTrack] setRThetaPhi input parameters: std::sqrt(energy*energy - mass*mass) = " << std::sqrt(energy*energy - mass*mass) << " momTheta = "<< momTheta <<" momPhi = "<< momPhi);
829 
830  // create dynamic particle
831  G4DynamicParticle dynParticle(secG4Particle,momVec);
832  newSecTrack = fastStep.CreateSecondaryTrack(dynParticle, posVec, currentTime, false); // use position in global coordinates instead
833 
834  return newSecTrack;
835 }

◆ dotProduct()

std::vector< double > PunchThroughG4Tool::dotProduct ( const std::vector< std::vector< double >> &  m,
const std::vector< double > &  v 
)
staticprivate

Definition at line 841 of file PunchThroughG4Tool.cxx.

842 {
843  std::vector<double> result;
844  result.reserve(m.size());
845  for (const auto& r : m){
846  result.push_back(std::inner_product(v.begin(), v.end(), r.begin(), 0.0));
847  }
848 
849  return result;
850 }

◆ finalize()

StatusCode PunchThroughG4Tool::finalize ( )
overridevirtual

AlgTool finalize method.

Definition at line 265 of file PunchThroughG4Tool.cxx.

266 {
267  ATH_MSG_DEBUG( "[PunchThroughG4Tool] finalize() starting" );
268  for(auto & each : m_particles) {
269  delete each.second;
270  }
271 
272  ATH_MSG_DEBUG( "[PunchThroughG4Tool] finalize() successful" );
273 
274  return StatusCode::SUCCESS;
275 }

◆ getAllParticles()

int PunchThroughG4Tool::getAllParticles ( const G4Track &  g4PrimaryTrack,
std::vector< std::map< std::string, double >> &  secKinematicsMapVect,
CLHEP::HepRandomEngine *  rndmEngine,
int  pdg,
double  interpEnergy,
double  interpEta,
int  numParticles = -1 
)
private

create the right number of punch-through particles for the given pdg and return the number of particles which was created.

also create these particles with the right distributions (energy, theta, phi). if a second argument is given, create exactly this number of particles (also with the right energy,theta,phi distributions

Definition at line 448 of file PunchThroughG4Tool.cxx.

449 {
450  // Get initial Geant4 position theta and phi
451  float initParticleTheta = g4PrimaryTrack.GetPosition().theta();
452  float initParticlePhi = g4PrimaryTrack.GetPosition().phi();
453 
454  // get the current particle
456  double minAllowedEnergy = p->getMinEnergy();
457 
458  // if no number of particles (=-1) was handed over as an argument
459  // -> get the number of particles from the pdf
460  if ( numParticles < 0 )
461  {
462  // prepare the function arguments for the PunchThroughPDFCreator class
463  std::vector<int> parameters;
464  parameters.push_back( std::round(interpEnergy) );
465  parameters.push_back( std::round(interpEta*100) );
466  // the maximum number of particles which should be produced
467  // if no maximum number is given, this is -1
468  int maxParticles = p->getMaxNumParticles();
469 
470  // get the right number of punch-through particles
471  // and ensure that we do not create too many particles
472  do
473  {
474  numParticles = int( p->getNumParticlesPDF()->getRand(rndmEngine, parameters) );
475 
476  // scale the number of particles if requested
477  numParticles = lround( numParticles *= p->getNumParticlesFactor() );
478  }
479  while ( (maxParticles >= 0.) && (numParticles > maxParticles) );
480  }
481 
482  ATH_MSG_DEBUG("[PunchThroughG4Tool] adding " << numParticles << " punch-through particles with pdg code: " << pdg);
483 
484  // get how many secondary particles to be created
485  for ( int numCreated = 0; numCreated < numParticles; numCreated++ )
486  {
487  // create one particle which fullfills the right energy distribution
488  std::map<std::string, double> secondaryKinematicsMap = getOneParticleKinematics(rndmEngine, pdg, initParticleTheta, initParticlePhi, interpEnergy, interpEta);
489 
490  //safety to cater minimum energy as rest mass
491  double energy = secondaryKinematicsMap.at("energy");
492  if(energy > minAllowedEnergy){
493  // push to the vector of maps
494  secKinematicsMapVect.push_back(secondaryKinematicsMap);
495  }
496  }
497 
498  // Get the size of the vector
499  std::size_t numSecondaries = secKinematicsMapVect.size();
500 
501  // the actual (pre) number of particles which was created is numSecondaries, which is incremented in the loop
502  return (numSecondaries);
503 }

◆ getCaloMSVars()

std::vector< double > PunchThroughG4Tool::getCaloMSVars ( )
overridevirtual

Definition at line 197 of file PunchThroughG4Tool.cxx.

197  {
198  //m_R1, m_R2, m_z1, m_z2
199  std::vector<double> caloMSVarsVect = {m_R1, m_R2, m_z1, m_z2};
200  ATH_MSG_DEBUG("[PunchThroughG4Tool] getCaloMSVars()==> m_R1 = "<< caloMSVarsVect[0] << ", m_R2 = " << caloMSVarsVect[1] << ", m_z1 = " << caloMSVarsVect[2] << ", m_z2 = " << caloMSVarsVect[3]);
201 
202  return caloMSVarsVect;
203 }

◆ getCorrelatedParticles()

int PunchThroughG4Tool::getCorrelatedParticles ( const G4Track &  g4PrimaryTrack,
std::vector< std::map< std::string, double >> &  secKinematicsMapVect,
int  pdg,
int  corrParticles,
CLHEP::HepRandomEngine *  rndmEngine,
double  interpEnergy,
double  interpEta 
)
private

get the right number of particles for the given pdg while considering the correlation to an other particle type, which has already created 'corrParticles' number of particles

Definition at line 505 of file PunchThroughG4Tool.cxx.

506 {
507  // Get Geant4 particle definition
508  const G4ParticleDefinition * mainG4Particle = g4PrimaryTrack.GetDefinition();
509  // Get Geant4 particle pdgID
510  float mainMomMag2 = g4PrimaryTrack.GetMomentum().mag2();
511  float mainPartMass = mainG4Particle->GetPDGMass();
512 
513  // get the PunchThroughParticle class
515 
516  const double initEnergy = std::sqrt( mainMomMag2 + mainPartMass*mainPartMass );
517 
518  // (1.) decide if we do correlation or not
519  double rand = CLHEP::RandFlat::shoot(rndmEngine)
520  *(p->getFullCorrelationEnergy() - p->getMinCorrelationEnergy())
521  + p->getMinCorrelationEnergy();
522 
523  // if initEnergy less than random corr energy (meaning threshold for min energy for doing correlation not passed)
524  if ( initEnergy < rand )
525  {
526  // here we do not do correlation
527  return getAllParticles(g4PrimaryTrack, secKinematicsMapVect, rndmEngine, pdg, interpEnergy, interpEta);
528  }
529 
530  // (2.) if this point is reached, we do correlation
531  // decide which 2d correlation histogram to use
532  double *histDomains = p->getCorrelationHistDomains();
533  TH2F *hist2d = nullptr;
534  // compute the center values of the lowE and highE
535  // correlation histogram domains
536  if ( initEnergy < histDomains[1])
537  {
538  // initial energy lower than border between lowEnergy and highEnergy histogram domain
539  // --> choose lowEnergy correlation histogram
540  hist2d = p->getCorrelationLowEHist();
541  }
542  else
543  {
544  double rand = CLHEP::RandFlat::shoot(rndmEngine)*(histDomains[2]-histDomains[1]) + histDomains[1];
545  hist2d = ( initEnergy < rand) ? p->getCorrelationLowEHist() : p->getCorrelationHighEHist();
546  }
547 
548  // get the correlation 2d histogram
549  // now find out where on the x-axis the the bin for number of
550  // correlated particles is
551  Int_t xbin = hist2d->GetXaxis()->FindFixBin(corrParticles);
552  int numParticles = 0;
553  int maxParticles = p->getMaxNumParticles();
554  // now the distribution along the y-axis is a PDF for the number
555  // of 'pdg' particles
556  do
557  {
558  double rand = CLHEP::RandFlat::shoot(rndmEngine);
559  double sum = 0.;
560  for ( int ybin = 1; ybin <= hist2d->GetNbinsY(); ybin++ )
561  {
562  sum += hist2d->GetBinContent(xbin, ybin);
563  // check if we choose the current bin or not
564  if ( sum >= rand )
565  {
566  numParticles = ybin - 1;
567  break;
568  }
569  }
570  // scale the number of particles is requested
571  numParticles = lround( numParticles * p->getNumParticlesFactor() );
572  }
573  while ( (maxParticles >= 0.) && (numParticles > maxParticles) );
574 
575  // finally create this exact number of particles
576  return getAllParticles(g4PrimaryTrack, secKinematicsMapVect, rndmEngine, pdg, interpEnergy, interpEta, numParticles);
577 }

◆ getFloatAfterPatternInStr()

double PunchThroughG4Tool::getFloatAfterPatternInStr ( const char *  str,
const char *  pattern 
)
private

get the floating point number in a string, after the given pattern

Definition at line 1440 of file PunchThroughG4Tool.cxx.

1441 {
1442  double num = 0.;
1443 
1444  const std::string_view str( cstr);
1445  const std::string_view pattern( cpattern);
1446  const size_t pos = str.find(pattern);
1447 
1448  if ( pos == std::string::npos)
1449  {
1450  ATH_MSG_WARNING("[PunchThroughG4Tool] unable to retrieve floating point number from string");
1451  return -999999999999.;
1452  }
1453  const std::string_view substring = str.substr(pos+pattern.length());
1454  std::from_chars(substring.data(), substring.data() + substring.size(), num);
1455  return num;
1456 }

◆ getInfoMap()

std::vector< std::map< std::string, std::string > > PunchThroughG4Tool::getInfoMap ( const std::string &  mainNode,
const std::string &  xmlFilePath 
)
private

Definition at line 905 of file PunchThroughG4Tool.cxx.

905  {
906  // Initialize pointers
907  xmlDocPtr doc;
908  xmlChar* xmlBuff = nullptr;
909 
910  std::vector<std::map<std::string, std::string>> xml_info;
911  doc = xmlParseFile( xmlFilePath.c_str() );
912 
913  //check info first
914  for( xmlNodePtr nodeRoot = doc->children; nodeRoot != nullptr; nodeRoot = nodeRoot->next) {
915  if (xmlStrEqual( nodeRoot->name, BAD_CAST mainNode.c_str() )) {
916  for( xmlNodePtr nodeRootChild = nodeRoot->children; nodeRootChild != nullptr; nodeRootChild = nodeRootChild->next ) {
917  if (xmlStrEqual( nodeRootChild->name, BAD_CAST "info" )) {
918  if (nodeRootChild->children != NULL) {
919  for( xmlNodePtr infoNode = nodeRootChild->children; infoNode != nullptr; infoNode = infoNode->next) {
920  if(xmlStrEqual( infoNode->name, BAD_CAST "item" )){
921  std::map<std::string, std::string> xml_info_item;
922 
923  if ((xmlBuff = xmlGetProp(infoNode, BAD_CAST "name")) != nullptr) {
924  xml_info_item.insert({"name", reinterpret_cast<const char*>(xmlBuff)});
925  }
926  if ((xmlBuff = xmlGetProp(infoNode, BAD_CAST "etaMins")) != nullptr) {
927  xml_info_item.insert({"etaMins", reinterpret_cast<const char*>(xmlBuff)});
928  }
929  if ((xmlBuff = xmlGetProp(infoNode, BAD_CAST "etaMaxs")) != nullptr) {
930  xml_info_item.insert({"etaMaxs", reinterpret_cast<const char*>(xmlBuff)});
931  }
932  if ((xmlBuff = xmlGetProp(infoNode, BAD_CAST "pidStr")) != nullptr) {
933  xml_info_item.insert({"pidStr", reinterpret_cast<const char*>(xmlBuff)});
934  }
935 
936  xml_info.push_back(xml_info_item);
937  }
938  }
939  }
940  }
941  }
942  }
943  }
944 
945  // free memory when done
946  xmlFreeDoc(doc);
947 
948  return xml_info;
949 }

◆ getOneParticleKinematics()

std::map< std::string, double > PunchThroughG4Tool::getOneParticleKinematics ( CLHEP::HepRandomEngine *  rndmEngine,
int  secondaryPDG,
float  initParticleTheta,
float  initParticlePhi,
double  interpEnergy,
double  interpEta 
) const
private

create exactly one punch-through particle with the given pdg and the given max energy

Definition at line 579 of file PunchThroughG4Tool.cxx.

580 {
581  // get a local copy of the needed punch-through particle class
582  PunchThroughParticle *p = m_particles.at(secondaryPDG);
583 
584  // (0.) get the pca / cdf group based on pdgId and eta, eta times 100, e.g eta -4 to 4 is from eta -400 to 400
585  int pcaCdfIterator = passedParamIterator(secondaryPDG, interpEta*100, m_xml_info_pca); //pca and cdf info should be of same size
586 
587  ATH_MSG_DEBUG("[PunchThroughG4Tool] passedPCAIterator ==> pcaCdfIterator = "<< pcaCdfIterator <<" , pdg = "<< secondaryPDG <<" , interpEnergy = "<< interpEnergy <<" MeV, interpEta(*100) = "<< interpEta*100);
588 
589  // (1.) decide if we create a particle or an anti-particle,
590  // this is from m_doAntiParticles properties
591  int anti = 1;
592  if ( p->getdoAnti() )
593  {
594  // get a random-value
595  double rand = CLHEP::RandFlat::shoot(rndmEngine);
596  // 50/50 chance to be a particle or its anti-particle
597  if (rand > 0.5) anti = -1;
598  }
599 
600  // (2.) get the right punch-through distributions
601  // prepare the function arguments for the PunchThroughPDFCreator class
602  std::vector<int> parInitEnergyEta;
603  parInitEnergyEta.push_back( std::round(interpEnergy) );
604  parInitEnergyEta.push_back( std::round(interpEta*100) );
605 
606  //initialise variables to store punch through particle kinematics
607  double energy = 0.;
608  double deltaTheta = 0.;
609  double deltaPhi = 0.;
610  double momDeltaTheta = 0.;
611  double momDeltaPhi = 0.;
612 
613  double principal_component_0 = 0.;
614  double principal_component_1 = 0.;
615  double principal_component_2 = 0.;
616  double principal_component_3 = 0.;
617  double principal_component_4 = 0.;
618  std::vector<double> transformed_variables;
619 
620  principal_component_0 = p->getPCA0PDF()->getRand(rndmEngine, parInitEnergyEta);
621  principal_component_1 = p->getPCA1PDF()->getRand(rndmEngine, parInitEnergyEta);
622  principal_component_2 = p->getPCA2PDF()->getRand(rndmEngine, parInitEnergyEta);
623  principal_component_3 = p->getPCA3PDF()->getRand(rndmEngine, parInitEnergyEta);
624  principal_component_4 = p->getPCA4PDF()->getRand(rndmEngine, parInitEnergyEta);
625 
626  ATH_MSG_DEBUG("[PunchThroughG4Tool] Drawn punch through kinematics PCA components: PCA0 = "<< principal_component_0 <<" PCA1 = "<< principal_component_1 <<" PCA2 = "<< principal_component_2 <<" PCA3 = "<< principal_component_3 <<" PCA4 = "<< principal_component_4 );
627 
628  std::vector<double> principal_components {
629  principal_component_0,
630  principal_component_1,
631  principal_component_2,
632  principal_component_3,
633  principal_component_4
634  };
635 
636  transformed_variables = inversePCA(pcaCdfIterator,principal_components);
637 
638  energy = inverseCdfTransform(transformed_variables.at(0), m_variable0_inverse_cdf[pcaCdfIterator]);
639  deltaTheta = inverseCdfTransform(transformed_variables.at(1), m_variable1_inverse_cdf[pcaCdfIterator]);
640  deltaPhi = inverseCdfTransform(transformed_variables.at(2), m_variable2_inverse_cdf[pcaCdfIterator]);
641  momDeltaTheta = inverseCdfTransform(transformed_variables.at(3), m_variable3_inverse_cdf[pcaCdfIterator]);
642  momDeltaPhi = inverseCdfTransform(transformed_variables.at(4), m_variable4_inverse_cdf[pcaCdfIterator]);
643 
644  ATH_MSG_DEBUG("[PunchThroughG4Tool] Transformed punch through kinematics: energy = "<< energy <<" MeV deltaTheta = "<< deltaTheta <<" deltaPhi = "<< deltaPhi <<" momDeltaTheta = "<< momDeltaTheta <<" momDeltaPhi = "<< momDeltaPhi );
645 
646  energy *= p->getEnergyFactor(); // scale the energy if requested
647 
648  // safety: if energy less than minimum energy, add a small energy (10MeV to it to continue)
649  if (energy < p->getMinEnergy()) {
650  energy = p->getMinEnergy() + 10;
651  }
652 
653  // (2.2) get the particles DeltaTheta relative to the incoming particle
654  double theta = 0;
655  // loop to keep theta within range [0,PI]
656  do
657  {
658  // decide if delta positive/negative
659  deltaTheta *= ( CLHEP::RandFlat::shoot(rndmEngine) > 0.5 ) ? 1. : -1.;
660 
661  // calculate the exact theta value of the later created
662  // punch-through particle
663  theta = initParticleTheta + deltaTheta*p->getPosAngleFactor();
664  }
665  while ( (theta > M_PI) || (theta < 0.) );
666 
667  // (2.3) get the particle's delta phi relative to the incoming particle
668  deltaPhi *= ( CLHEP::RandFlat::shoot(rndmEngine) > 0.5 ) ? 1. : -1.;
669 
670  // keep phi within range [-PI,PI]
671  double phi = initParticlePhi + deltaPhi*p->getPosAngleFactor();
672  while ( std::fabs(phi) > 2*M_PI) phi /= 2.;
673  while (phi > M_PI) phi -= 2*M_PI;
674  while (phi < -M_PI) phi += 2*M_PI;
675 
676  // (2.4) get the particle momentum delta theta, relative to its position
677  //
678  // loop to keep momTheta within range [0,PI]
679  double momTheta = 0.;
680  do
681  {
682  // decide if delta positive/negative
683  momDeltaTheta *= ( CLHEP::RandFlat::shoot(rndmEngine) > 0.5 ) ? 1. : -1.;
684  // calculate the exact momentum theta value of the later created
685  // punch-through particle
686  momTheta = theta + momDeltaTheta*p->getMomAngleFactor();
687  }
688  while ( (momTheta > M_PI) || (momTheta < 0.) );
689 
690  // (2.5) get the particle momentum delta phi, relative to its position
691  momDeltaPhi *= ( CLHEP::RandFlat::shoot(rndmEngine) > 0.5 ) ? 1. : -1.;
692 
693  double momPhi = phi + momDeltaPhi*p->getMomAngleFactor();
694  // keep momPhi within range [-PI,PI]
695  while ( std::fabs(momPhi) > 2*M_PI) momPhi /= 2.;
696  while (momPhi > M_PI) momPhi -= 2*M_PI;
697  while (momPhi < -M_PI) momPhi += 2*M_PI;
698 
699  // store the kinematics in a map and return it
700  std::map<std::string, double> secondaryKinematicsMap;
701  secondaryKinematicsMap.insert({ "anti" , anti });
702  secondaryKinematicsMap.insert({ "secondaryPDG" , secondaryPDG });
703  secondaryKinematicsMap.insert({ "energy" , energy });
704  secondaryKinematicsMap.insert({ "theta" , theta });
705  secondaryKinematicsMap.insert({ "phi" , phi });
706  secondaryKinematicsMap.insert({ "momTheta" , momTheta });
707  secondaryKinematicsMap.insert({ "momPhi" , momPhi });
708 
709  return secondaryKinematicsMap;
710 }

◆ getVariableCDFmappings()

std::map< double, double > PunchThroughG4Tool::getVariableCDFmappings ( xmlNodePtr &  nodeParent)
staticprivate

Definition at line 1088 of file PunchThroughG4Tool.cxx.

1088  {
1089  std::map<double, double> mappings;
1090  xmlChar* xmlBuff = nullptr;
1091  double ref = -1;
1092  double quant = -1;
1093  for( xmlNodePtr node = nodeParent->children; node != nullptr; node = node->next ) {
1094  //Get min and max values that we normalise values to
1095  if (xmlStrEqual( node->name, BAD_CAST "CDFmap" )) {
1096  if ((xmlBuff = xmlGetProp(node, BAD_CAST "ref")) != nullptr) {
1097  ref = atof( reinterpret_cast<const char*> (xmlBuff) );
1098  }
1099  if ((xmlBuff = xmlGetProp(node, BAD_CAST "quant")) != nullptr) {
1100  quant = atof( reinterpret_cast<const char*> (xmlBuff) );
1101  }
1102  mappings.insert(std::pair<double, double>(ref, quant) );
1103  }
1104  }
1105 
1106  return mappings;
1107 }

◆ initialize()

StatusCode PunchThroughG4Tool::initialize ( )
overridevirtual

AlgTool initialize method.

Definition at line 45 of file PunchThroughG4Tool.cxx.

45  {
46  ATH_MSG_DEBUG("[PunchThroughG4Tool] ==>" << name() << "::initialize()");
47 
48  //retrieve inverse CDF config file
50  {
51  ATH_MSG_WARNING("[PunchThroughG4Tool] unable to open or read the inverse CDF config");
52  return StatusCode::FAILURE;
53  }
54 
55  //retrieve inverse PCA config file
57  {
58  ATH_MSG_WARNING("[PunchThroughG4Tool] unable to open or read the inverse PCA config");
59  return StatusCode::FAILURE;
60  }
61 
62  //check first the size of infoMap for both PCA and CDF, they should be equal
63  if (!(m_xml_info_pca.size() == m_xml_info_cdf.size()))
64  {
65  ATH_MSG_WARNING("[PunchThroughG4Tool] size of infoMap for PCA and CDF differs! Something is wrong with input xml files.");
66  return StatusCode::FAILURE;
67  }
68 
69  ATH_MSG_INFO( "[PunchThroughG4Tool] initialization is successful" );
70  return StatusCode::SUCCESS;
71 }

◆ initializeInverseCDF()

StatusCode PunchThroughG4Tool::initializeInverseCDF ( const std::string &  quantileTransformerConfigFile)
private

Definition at line 1023 of file PunchThroughG4Tool.cxx.

1023  {
1024  std::map<double, double> variable0_inverse_cdf_row;
1025  std::map<double, double> variable1_inverse_cdf_row;
1026  std::map<double, double> variable2_inverse_cdf_row;
1027  std::map<double, double> variable3_inverse_cdf_row;
1028  std::map<double, double> variable4_inverse_cdf_row;
1029 
1030  // Initialize pointers
1031  xmlDocPtr doc;
1032 
1033  //parse xml that contains config for inverse CDF for each of punch through particle kinematics
1034 
1035  doc = xmlParseFile( inverseCdfConfigFile.c_str() );
1036 
1037  ATH_MSG_DEBUG( "[PunchThroughG4Tool] Loading inverse CDF: " << inverseCdfConfigFile);
1038 
1039  //check info first
1040  m_xml_info_cdf = getInfoMap("CDFMappings",inverseCdfConfigFile);
1041 
1042  //do the saving
1043  for (unsigned int i = 0; i < m_xml_info_cdf.size(); i++) {
1044  ATH_MSG_DEBUG( "[PunchThroughG4Tool] m_xml_info_cdf[" << i << "].at('name') = " << m_xml_info_cdf[i].at("name"));
1045 
1046  for( xmlNodePtr nodeRoot = doc->children; nodeRoot != nullptr; nodeRoot = nodeRoot->next) {
1047  if (xmlStrEqual( nodeRoot->name, BAD_CAST "CDFMappings" )) {
1048  for( xmlNodePtr typeMappings = nodeRoot->children; typeMappings != nullptr; typeMappings = typeMappings->next ) {
1049  if (xmlStrEqual( typeMappings->name, BAD_CAST m_xml_info_cdf[i].at("name").c_str() )) {
1050  if (typeMappings->children != NULL) {
1051  for( xmlNodePtr nodeMappings = typeMappings->children; nodeMappings != nullptr; nodeMappings = nodeMappings->next) {
1052 
1053  if (xmlStrEqual( nodeMappings->name, BAD_CAST "variable0" )) {
1054  variable0_inverse_cdf_row = getVariableCDFmappings(nodeMappings);
1055  }
1056  else if (xmlStrEqual( nodeMappings->name, BAD_CAST "variable1" )) {
1057  variable1_inverse_cdf_row = getVariableCDFmappings(nodeMappings);
1058  }
1059  else if (xmlStrEqual( nodeMappings->name, BAD_CAST "variable2" )) {
1060  variable2_inverse_cdf_row = getVariableCDFmappings(nodeMappings);
1061  }
1062  else if (xmlStrEqual( nodeMappings->name, BAD_CAST "variable3" )) {
1063  variable3_inverse_cdf_row = getVariableCDFmappings(nodeMappings);
1064  }
1065  else if (xmlStrEqual( nodeMappings->name, BAD_CAST "variable4" )) {
1066  variable4_inverse_cdf_row = getVariableCDFmappings(nodeMappings);
1067  }
1068  }
1069 
1070  }
1071  }
1072  }
1073  }
1074  }
1075  m_variable0_inverse_cdf.push_back(variable0_inverse_cdf_row);
1076  m_variable1_inverse_cdf.push_back(variable1_inverse_cdf_row);
1077  m_variable2_inverse_cdf.push_back(variable2_inverse_cdf_row);
1078  m_variable3_inverse_cdf.push_back(variable3_inverse_cdf_row);
1079  m_variable4_inverse_cdf.push_back(variable4_inverse_cdf_row);
1080  }
1081 
1082  // free memory when done
1083  xmlFreeDoc(doc);
1084 
1085  return StatusCode::SUCCESS;
1086 }

◆ initializeInversePCA()

StatusCode PunchThroughG4Tool::initializeInversePCA ( const std::string &  inversePCAConfigFile)
private

Definition at line 960 of file PunchThroughG4Tool.cxx.

960  {
961  // Initialize pointers
962  xmlDocPtr doc;
963  xmlChar* xmlBuff = nullptr;
964 
965  doc = xmlParseFile( inversePCAConfigFile.c_str() );
966 
967  ATH_MSG_DEBUG( "[PunchThroughG4Tool] Loading inversePCA: " << inversePCAConfigFile);
968 
969  //check info first
970  m_xml_info_pca = getInfoMap("PCAinverse",inversePCAConfigFile);
971 
972  //do the saving
973  for (unsigned int i = 0; i < m_xml_info_pca.size(); i++) {
974  std::vector<std::vector<double>> PCA_matrix;
975  ATH_MSG_DEBUG( "[PunchThroughG4Tool] m_xml_info_pca[" << i << "].at('name') = " << m_xml_info_pca[i].at("name"));
976 
977  for( xmlNodePtr nodeRoot = doc->children; nodeRoot != nullptr; nodeRoot = nodeRoot->next) {
978  if (xmlStrEqual( nodeRoot->name, BAD_CAST "PCAinverse" )) {
979  for( xmlNodePtr nodePCAinverse = nodeRoot->children; nodePCAinverse != nullptr; nodePCAinverse = nodePCAinverse->next ) {
980 
981  if (xmlStrEqual( nodePCAinverse->name, BAD_CAST m_xml_info_pca[i].at("name").c_str() )) {
982  if (nodePCAinverse->children != NULL) {
983  for( xmlNodePtr pcaNode = nodePCAinverse->children; pcaNode != nullptr; pcaNode = pcaNode->next) {
984 
985  if (xmlStrEqual( pcaNode->name, BAD_CAST "PCAmatrix" )) {
986  std::vector<double> PCA_matrix_row;
987  for (int i = 0; i <= 4; ++i) {
988  std::string propName = "comp_" + std::to_string(i); // Dynamically create property name
989  if ((xmlBuff = xmlGetProp(pcaNode, BAD_CAST propName.c_str())) != nullptr) {
990  PCA_matrix_row.push_back(atof(reinterpret_cast<const char*>(xmlBuff))); // Convert and push to the row
991  }
992  }
993  PCA_matrix.push_back(PCA_matrix_row);
994  }
995  else if (xmlStrEqual( pcaNode->name, BAD_CAST "PCAmeans" )) {
996  std::vector<double> PCA_means_row;
997  for (int i = 0; i <= 4; ++i) {
998  std::string propName = "mean_" + std::to_string(i); // Dynamically create property name
999  if ((xmlBuff = xmlGetProp(pcaNode, BAD_CAST propName.c_str())) != nullptr) {
1000  PCA_means_row.push_back(atof(reinterpret_cast<const char*>(xmlBuff))); // Convert and push to the row
1001  }
1002  }
1003  m_PCA_means.push_back(PCA_means_row);
1004  }
1005 
1006  }
1007 
1008  }
1009  }
1010 
1011  }
1012  }
1013  }
1014  m_inverse_PCA_matrix.push_back(PCA_matrix);
1015  }
1016 
1017  // free memory when done
1018  xmlFreeDoc(doc);
1019 
1020  return StatusCode::SUCCESS;
1021 }

◆ initializePhysics()

StatusCode PunchThroughG4Tool::initializePhysics ( )
override

Definition at line 105 of file PunchThroughG4Tool.cxx.

105  {
106  // resolving lookuptable file
107 
108  ATH_MSG_INFO("[PunchThroughG4Tool] PunchThroughG4Tool::initializePhysics() called");
109  std::string resolvedFileName = PathResolverFindCalibFile (m_filenameLookupTable);
110  if (resolvedFileName.empty()) {
111  ATH_MSG_ERROR( "[PunchThroughG4Tool] Parametrisation file '" << m_filenameLookupTable << "' not found" );
112  return StatusCode::FAILURE;
113  }
114  ATH_MSG_DEBUG( "[PunchThroughG4Tool] Parametrisation file found: " << resolvedFileName );
115  // open the LookupTable file
116  m_fileLookupTable = new TFile( resolvedFileName.c_str(), "READ");
117  if (!m_fileLookupTable) {
118  ATH_MSG_WARNING("[PunchThroughG4Tool] unable to open the lookup-table for the punch-through simulation (file does not exist)");
119  return StatusCode::FAILURE;
120  }
121 
122  if (!m_fileLookupTable->IsOpen()) {
123  ATH_MSG_WARNING("[PunchThroughG4Tool] unable to open the lookup-table for the punch-through simulation (wrong or empty file?)");
124  return StatusCode::FAILURE;
125  }
126 
127  //--------------------------------------------------------------------------------
128  // register all the punch-through particles which will be simulated
130 
131  //--------------------------------------------------------------------------------
132  // register all the punch-through correlations between particles to be simulated
134 
135  //--------------------------------------------------------------------------------
136  // close the file with the lookuptable
137  m_fileLookupTable->Close();
138 
139  // Geometry identifier service
140  if ( !m_geoIDSvc.empty() && m_geoIDSvc.retrieve().isFailure())
141  {
142  ATH_MSG_FATAL ( "[PunchThroughG4Tool] Could not retrieve GeometryIdentifier Service. Abort");
143  return StatusCode::FAILURE;
144  }
145 
146  //envelope definition service
147  if (m_envDefSvc.retrieve().isFailure() )
148  {
149  ATH_MSG_ERROR( "[PunchThroughG4Tool] Could not retrieve " << m_envDefSvc );
150  return StatusCode::FAILURE;
151  }
152 
153  // direct-initialization of the Muon and Calo boundaries (giving R, Z coordinates)
154  // can also do like below (RZPairVector is vector<pair<double, double>>)
155  // const RZPairVector* rzMS = &(m_envDefSvc->getMuonRZBoundary());
156  const std::vector<std::pair<double, double>>* rzMS = &(m_envDefSvc->getMuonRZBoundary());
157  const std::vector<std::pair<double, double>>* rzCalo = &(m_envDefSvc->getCaloRZBoundary());
158  ATH_CHECK(checkCaloMSBoundaries(rzMS, rzCalo));
159  return StatusCode::SUCCESS;
160 }

◆ initializeRegisterCorrelations()

StatusCode PunchThroughG4Tool::initializeRegisterCorrelations ( )
private

initialize register all correlations between particles

Definition at line 162 of file PunchThroughG4Tool.cxx.

162  {
163  //--------------------------------------------------------------------------------
164  // TODO: implement punch-through parameters for different m_pdgInitiators
165  // currently m_pdgInitiators is only used to filter out particles,
166  // inside computePunchThroughParticles
167 
168  // check if more correlations were given than particle types were registered
169  unsigned int numCorrelations = m_correlatedParticle.size();
170  if ( numCorrelations > m_punchThroughParticles.size() )
171  {
172  ATH_MSG_WARNING("[PunchThroughG4Tool] more punch-through particle correlations are given, than punch-through particle types are registered (skipping the last ones)");
173  numCorrelations = m_punchThroughParticles.size();
174  }
175 
176  // now register correlation between particles
177  for ( unsigned int num = 0; num < numCorrelations; num++ )
178  {
179  const int pdg1 = m_punchThroughParticles[num];
180  const int pdg2 = m_correlatedParticle[num];
181  const double fullCorrEnergy = ( num < m_fullCorrEnergy.size() ) ? m_fullCorrEnergy[num] : 0.;
182  const double minCorrEnergy = ( num < m_minCorrEnergy.size() ) ? m_minCorrEnergy[num] : 0.;
183 
184  // if correlatedParticle==0 is given -> no correlation
185  if ( ! pdg2) continue;
186  // register it
187  if ( registerCorrelation(pdg1, pdg2, minCorrEnergy, fullCorrEnergy) != StatusCode::SUCCESS )
188  {
189  ATH_MSG_ERROR("[PunchThroughG4Tool] unable to register punch-through particle correlation for pdg1=" << pdg1 << " pdg2=" << pdg2 );
190  return StatusCode::FAILURE;
191  }
192  }
193  // if all goes well
194  return StatusCode::SUCCESS;
195 }

◆ initializeRegisterPunchThroughParticles()

StatusCode PunchThroughG4Tool::initializeRegisterPunchThroughParticles ( )
private

initialize register all the punch-through particles which will be simulated

Definition at line 73 of file PunchThroughG4Tool.cxx.

73  {
74  G4ParticleTable *ptable = G4ParticleTable::GetParticleTable();
75  for ( unsigned int num = 0; num < m_punchThroughParticles.size(); num++ )
76  {
77  const int pdg = m_punchThroughParticles[num];
78  // if no information is given on the creation of anti-particles -> do not simulate anti-particles
79  const bool doAnti = ( num < m_doAntiParticles.size() ) ? m_doAntiParticles[num] : false;
80  // if no information is given on the minimum energy -> take 50. MeV as default
81  const double minEnergy = ( num < m_minEnergy.size() ) ? m_minEnergy[num] : 50.;
82  // if no information is given on the maximum number of punch-through particles -> take -1 as default
83  const int maxNum = ( num < m_minEnergy.size() ) ? m_maxNumParticles[num] : -1;
84  // if no information is given on the scale factor for the number of particles -> take 1. as defaulft
85  const double numFactor = ( num < m_numParticlesFactor.size() ) ? m_numParticlesFactor[num] : 1.;
86  // if no information is given on the position angle factor -> take 1.
87  const double posAngleFactor = ( num < m_posAngleFactor.size() ) ? m_posAngleFactor[num] : 1.;
88  // if no information is given on the momentum angle factor -> take 1.
89  const double momAngleFactor = ( num < m_momAngleFactor.size() ) ? m_momAngleFactor[num] : 1.;
90  // if no information is given on the scale factor for the energy -> take 1. as default
91  const double energyFactor = ( num < m_energyFactor.size() ) ? m_energyFactor[num] : 1.;
92 
93  // register the particle
94  ATH_MSG_VERBOSE("VERBOSE: [PunchThroughG4Tool] registering punch-through particle type with pdg = " << pdg );
95  if (registerPunchThroughParticle( *ptable, pdg, doAnti, minEnergy, maxNum, numFactor, energyFactor, posAngleFactor, momAngleFactor ) != StatusCode::SUCCESS)
96  {
97  ATH_MSG_ERROR("[PunchThroughG4Tool] unable to register punch-through particle type with pdg = " << pdg);
98  return StatusCode::FAILURE;
99  }
100  }
101  // if all goes well
102  return StatusCode::SUCCESS;
103 }

◆ interpolateEnergy()

double PunchThroughG4Tool::interpolateEnergy ( const double &  energy,
CLHEP::HepRandomEngine *  rndmEngine 
) const
private

Definition at line 1124 of file PunchThroughG4Tool.cxx.

1124  {
1125 
1126  ATH_MSG_DEBUG("[PunchThroughG4Tool] interpolating incoming energy: " << energy);
1127 
1128  std::string energyPointsString;
1129  for (auto element:m_energyPoints){
1130  energyPointsString += std::to_string(element) + " ";
1131  }
1132 
1133  ATH_MSG_DEBUG("[PunchThroughG4Tool] available energy points: " << energyPointsString);
1134 
1135  auto const upperEnergy = std::upper_bound(m_energyPoints.begin(), m_energyPoints.end(), energy);
1136 
1137  if(upperEnergy == m_etaPoints.end()){ //if no energy greater than input energy, choose greatest energy
1138  ATH_MSG_DEBUG("[PunchThroughG4Tool] incoming energy > largest energy point, returning greatest energy point: " << m_energyPoints.back());
1139  return m_energyPoints.back();
1140  }
1141  else if(upperEnergy == m_etaPoints.begin()){ //if smallest energy greater than input energy, choose smallest energy
1142  ATH_MSG_DEBUG("[PunchThroughG4Tool] incoming energy < smallest energy point, returning smallest energy point: " << *upperEnergy);
1143  return *upperEnergy;
1144  }
1145 
1146  ATH_MSG_DEBUG("[PunchThroughG4Tool] energy points upper_bound: "<< *upperEnergy);
1147 
1148  double randomShoot = CLHEP::RandFlat::shoot(rndmEngine);
1149 
1150  ATH_MSG_DEBUG("[PunchThroughG4Tool] Shooting random number: "<< randomShoot);
1151 
1152  double midPoint = *std::prev(upperEnergy)*M_SQRT2;
1153 
1154  if(energy < midPoint){ //if energy smaller than mid point in log(energy)
1155 
1156  double distance = std::abs(energy - *std::prev(upperEnergy))/((midPoint) - *std::prev(upperEnergy));
1157 
1158  ATH_MSG_DEBUG( "[PunchThroughG4Tool] incoming energy is closest to prev(upper_bound) in log(energy), distance: " << distance );
1159 
1160  if(randomShoot < distance){
1161  ATH_MSG_DEBUG( "[PunchThroughG4Tool] randomShoot < distance, returning upper_bound " << *upperEnergy );
1162  return *upperEnergy;
1163  }
1164  ATH_MSG_DEBUG( "[PunchThroughG4Tool] randomShoot > distance, returning prev(upper_bound) " << *std::prev(upperEnergy) );
1165 
1166  return *std::prev(upperEnergy);
1167  }
1168  else if(energy > midPoint){ //if energy greater than mid point in log(energy)
1169 
1170  double distance = std::abs(energy - *upperEnergy)/((*upperEnergy - midPoint));
1171 
1172  ATH_MSG_DEBUG( "[PunchThroughG4Tool] incoming energy is closest to upper_bound in log(energy), distance: " << distance );
1173 
1174  if(randomShoot < distance){
1175  ATH_MSG_DEBUG( "[PunchThroughG4Tool] randomShoot < distance, returning prev(upper_bound) " << *std::prev(upperEnergy) );
1176  return *std::prev(upperEnergy);
1177  }
1178  ATH_MSG_DEBUG( "[PunchThroughG4Tool] randomShoot > distance, returning upper_bound " << *upperEnergy );
1179  return *upperEnergy;
1180  }
1181 
1182  return *upperEnergy;
1183 }

◆ interpolateEta()

double PunchThroughG4Tool::interpolateEta ( const double &  eta,
CLHEP::HepRandomEngine *  rndmEngine 
) const
private

Definition at line 1185 of file PunchThroughG4Tool.cxx.

1185  {
1186 
1187  double absEta = std::abs(eta);
1188 
1189  ATH_MSG_DEBUG("[PunchThroughG4Tool] interpolating incoming abs(eta): " << absEta);
1190 
1191  std::string etaPointsString;
1192  for (auto element:m_etaPoints){
1193  etaPointsString += std::to_string(element) + " ";
1194  }
1195 
1196  ATH_MSG_DEBUG("[PunchThroughG4Tool] available eta points: " << etaPointsString);
1197 
1198  auto const upperEta = std::upper_bound(m_etaPoints.begin(), m_etaPoints.end(), absEta);
1199 
1200  if(upperEta == m_etaPoints.end()){
1201  ATH_MSG_DEBUG("[PunchThroughG4Tool] incoming abs(eta) > largest eta point, returning greatest eta point: " << m_etaPoints.back());
1202  return m_etaPoints.back();
1203  }
1204 
1205 
1206  ATH_MSG_DEBUG("[PunchThroughG4Tool] eta points upper_bound: "<< *upperEta);
1207 
1208  double randomShoot = CLHEP::RandFlat::shoot(rndmEngine);
1209 
1210  ATH_MSG_DEBUG("[PunchThroughG4Tool] Shooting random number: "<< randomShoot);
1211 
1212  if(std::abs(absEta - *upperEta) < std::abs(absEta - *std::prev(upperEta))){
1213 
1214  double distance = std::abs(absEta - *upperEta)/((*upperEta - *std::prev(upperEta))/2);
1215 
1216  ATH_MSG_DEBUG( "[PunchThroughG4Tool] abs(eta) is closer to eta points upper_bound, distance: " << distance );
1217 
1218  if(randomShoot > distance){
1219  ATH_MSG_DEBUG( "[PunchThroughG4Tool] randomShoot > distance, returning upper_bound " << *upperEta );
1220  return *upperEta;
1221  }
1222 
1223  ATH_MSG_DEBUG( "[PunchThroughG4Tool] randomShoot < distance, returning prev(upper_bound) " << *std::prev(upperEta) );
1224 
1225  return *std::prev(upperEta);
1226  }
1227  else if(std::abs(absEta - *std::prev(upperEta)) < std::abs(absEta - *upperEta)){
1228 
1229  if(std::prev(upperEta) == m_etaPoints.begin()){
1230  ATH_MSG_DEBUG( "[PunchThroughG4Tool] prev of upper bound is begin, returning that: " << *std::prev(upperEta) );
1231  return *std::prev(upperEta);
1232  }
1233 
1234  double distance = std::abs(absEta - *std::prev(upperEta))/((*std::prev(upperEta) - *std::prev(std::prev(upperEta)))/2);
1235 
1236  ATH_MSG_DEBUG( "[PunchThroughG4Tool] abs(eta) is closer to eta points prev(upper_bound), distance: " << distance );
1237 
1238  if(randomShoot > distance){
1239  ATH_MSG_DEBUG( "[PunchThroughG4Tool] randomShoot > distance, returning prev(prev(upper_bound)) " << *std::prev(std::prev(upperEta)) );
1240 
1241  return *std::prev(std::prev(upperEta));
1242  }
1243  ATH_MSG_DEBUG( "[PunchThroughG4Tool] randomShoot < distance, returning prev(upper_bound) " << *std::prev(upperEta) );
1244 
1245  return *std::prev(upperEta);
1246  }
1247 
1248  return *std::prev(upperEta);
1249 }

◆ inverseCdfTransform()

double PunchThroughG4Tool::inverseCdfTransform ( double  variable,
const std::map< double, double > &  inverse_cdf_map 
)
staticprivate

Definition at line 1109 of file PunchThroughG4Tool.cxx.

1109  {
1110 
1111  double norm_cdf = normal_cdf(variable);
1112 
1113  auto upper = inverse_cdf_map.upper_bound(norm_cdf);
1114  auto lower = upper--;
1115 
1116  double m = (upper->second - lower->second)/(upper->first - lower->first);
1117  double c = lower->second - m * lower->first;
1118  double transformed = m * norm_cdf + c;
1119 
1120  return transformed;
1121 
1122 }

◆ inversePCA()

std::vector< double > PunchThroughG4Tool::inversePCA ( int  pcaCdfIterator,
std::vector< double > &  variables 
) const
private

Definition at line 951 of file PunchThroughG4Tool.cxx.

952 {
953  std::vector<double> transformed_variables = dotProduct(m_inverse_PCA_matrix[pcaCdfIterator], variables);
954 
955  std::transform (transformed_variables.begin(), transformed_variables.end(), m_PCA_means[pcaCdfIterator].begin(), transformed_variables.begin(), std::plus<double>()); // + means
956 
957  return transformed_variables;
958 }

◆ normal_cdf()

double PunchThroughG4Tool::normal_cdf ( double  x)
staticprivate

Definition at line 837 of file PunchThroughG4Tool.cxx.

837  {
838  return 0.5 * TMath::Erfc(-x * M_SQRT1_2);
839 }

◆ passedParamIterator()

int PunchThroughG4Tool::passedParamIterator ( int  pid,
double  eta,
const std::vector< std::map< std::string, std::string >> &  mapvect 
) const
private

Definition at line 876 of file PunchThroughG4Tool.cxx.

877 {
878  //convert the pid to absolute value and string for query
879  int pidStrSingle = std::abs(pid);
880  //STEP 1
881  //filter items matching pid first
882 
883  for (unsigned int i = 0; i < mapvect.size(); i++){
884  const std::string &pidStr = mapvect[i].at("pidStr");
885  auto v = str_to_list<int>(pidStr);
886  if(std::find(v.begin(), v.end(),pidStrSingle)==v.end()) continue;
887  const std::string &etaMinsStr = mapvect[i].at("etaMins");
888  const std::string &etaMaxsStr = mapvect[i].at("etaMaxs");
889  std::vector<double> etaMinsVect = str_to_list<double>(etaMinsStr);
890  std::vector<double> etaMaxsVect = str_to_list<double>(etaMaxsStr);
891  assert(etaMaxsVect.size() == etaMinsVect.size());
892  for (unsigned int j = 0; j < etaMinsVect.size(); j++){ // assume size etaMinsVect == etaMaxsVect
893  double etaMinToCompare = etaMinsVect[j];
894  double etaMaxToCompare = etaMaxsVect[j];
895  if((eta >= etaMinToCompare) && (eta < etaMaxToCompare)){
896  //PASS CONDITION
897  //then choose the passing one and note it's iterator
898  return (i); //in case more than 1 match (ambiguous case)
899  }
900  }
901  }
902  return 0;
903 }

◆ punchTroughPosPropagator()

G4ThreeVector PunchThroughG4Tool::punchTroughPosPropagator ( double  theta,
double  phi,
double  R1,
double  R2,
double  z1,
double  z2 
) const
private

get particle through the calorimeter

Definition at line 1459 of file PunchThroughG4Tool.cxx.

1460 {
1461  // phi, theta - direction of the punch-through particle coming into calo
1462  // particle propagates inside the calorimeter along the straight line
1463  // coordinates of this particles when exiting the calo (on calo-MS boundary)
1464 
1465  double x, y, z, r;
1466 
1467  // cylinders border angles
1468  const double theta1 = atan (R1 / z1);
1469  const double theta2 = atan (R1 / z2);
1470  const double theta3 = atan (R2 / z2);
1471  //where is the particle
1472 
1473  if (theta >= 0 && theta < theta1)
1474  {
1475  z = z1;
1476  r = std::fabs (z1 * tan(theta));
1477  }
1478  else if (theta >= theta1 && theta < theta2)
1479  {
1480  z = R1 / tan(theta);
1481  r = R1;
1482  }
1483  else if (theta >= theta2 && theta < theta3)
1484  {
1485  z = z2;
1486  r = std::fabs(z2 * tan(theta));;
1487  }
1488  else if (theta >= theta3 && theta < (TMath::Pi() - theta3) )
1489  {
1490  z = R2 / tan(theta);
1491  r = R2;
1492  }
1493  else if (theta >= (TMath::Pi() - theta3) && theta < (TMath::Pi() - theta2) )
1494  {
1495  z = -z2;
1496  r = std::fabs(z2 * tan(theta));
1497  }
1498  else if (theta >= (TMath::Pi() - theta2) && theta < (TMath::Pi() - theta1) )
1499  {
1500  z = R1 / tan(theta);
1501  r = R1;
1502  }
1503  else if (theta >= (TMath::Pi() - theta1) && theta <= TMath::Pi() )
1504  {
1505  z = -z1;
1506  r = std::fabs(z1 * tan(theta));
1507  }
1508 
1509  //parallel universe
1510  else
1511  {
1512  ATH_MSG_WARNING("[PunchThroughG4Tool::punchTroughPosPropagator] Given theta angle is incorrect, setting particle position to (0, 0, 0)");
1513  x = 0.0; y = 0.0; z = 0.0; r = 0.0;
1514  }
1515 
1516  x = r * cos(phi);
1517  y = r * sin(phi);
1518  G4ThreeVector posVec(x, y, z);
1519 
1520  ATH_MSG_DEBUG("[PunchThroughG4Tool::punchTroughPosPropagator] position of produced punch-through particle: x = " << x <<" y = "<< y <<" z = "<< z<<" r = "<< posVec.perp() <<"std::sqrt(x^2+y^2) = "<< std::sqrt(x * x + y * y));
1521 
1522  return posVec;
1523 }

◆ readLookuptablePDF()

std::unique_ptr< PunchThroughPDFCreator > PunchThroughG4Tool::readLookuptablePDF ( int  pdgID,
TFile *  fileLookupTable,
const std::string &  folderName 
)
private

reads out the lookuptable for the given type of particle

Definition at line 1373 of file PunchThroughG4Tool.cxx.

1374 {
1375 
1376  // will hold the PDFcreator class which will be returned at the end
1377  // this will store the distributions for the punch through particles
1378  // (as map of energy & eta of the incoming particle)
1379  //PunchThroughPDFCreator *pdf = new PunchThroughPDFCreator();
1380  std::unique_ptr<PunchThroughPDFCreator> pdf = std::make_unique<PunchThroughPDFCreator>();
1381 
1382  //Get directory object
1383  std::stringstream dirName;
1384  dirName << folderName << pdg;
1385  pdf->setName(dirName.str());
1386 
1387  TDirectory * dir = (TDirectory*)fileLookupTable->Get(dirName.str().c_str());
1388  if(! dir)
1389  {
1390  ATH_MSG_ERROR( "[PunchThroughG4Tool] unable to retrieve directory object ("<< folderName << pdg << ")" );
1391  return nullptr;
1392  }
1393 
1394  //Get list of all objects in directory
1395  TIter keyList(dir->GetListOfKeys());
1396  TKey *key;
1397 
1398  while ((key = (TKey*)keyList())) {
1399 
1400  //Get histogram object from key and its name
1401  TH1* hist = nullptr;
1402 
1403  std::string histName;
1404  if(strcmp(key->GetClassName(), "TH1F") == 0){
1405  hist = (TH1*)key->ReadObj();
1406  histName = hist->GetName();
1407  }
1408 
1409  //extract energy and eta from hist name 6 and 1 to position delimeters correctly
1410  std::string strEnergy = histName.substr( histName.find_first_of('E') + 1, histName.find_first_of('_')-histName.find_first_of('E') - 1 );
1411  histName.erase(0, histName.find_first_of('_') + 1);
1412  std::string strEtaMin = histName.substr( histName.find("etaMin") + 6, histName.find_first_of('_') - histName.find("etaMin") - 6 );
1413  histName.erase(0, histName.find('_') + 1);
1414  std::string strEtaMax = histName.substr( histName.find("etaMax") + 6, histName.length());
1415 
1416  //create integers to store in map
1417  const int energy = std::stoi(strEnergy);
1418  const int etaMin = std::stoi(strEtaMin);
1419 
1420  //Add entry to pdf map
1421  pdf->addToEnergyEtaHist1DMap(energy, etaMin, hist);
1422 
1423  //create doubles to store energy and eta points for interpolation
1424  const double energyDbl = static_cast<double>(energy);
1425  const double etaDbl = static_cast<double>(etaMin)/100.;
1426 
1427  //create vectors to store the eta and energy points, this allows us to interpolate
1428  if (std::find(m_energyPoints.begin(), m_energyPoints.end(), energyDbl) == m_energyPoints.end()) {
1429  m_energyPoints.push_back(energyDbl);
1430  }
1431  if (std::find(m_etaPoints.begin(), m_etaPoints.end(), etaDbl) == m_etaPoints.end()) {
1432  m_etaPoints.push_back(etaDbl);
1433  }
1434 
1435  }
1436 
1437  return pdf;
1438 }

◆ registerCorrelation()

StatusCode PunchThroughG4Tool::registerCorrelation ( int  pdgID1,
int  pdgID2,
double  minCorrEnergy = 0.,
double  fullCorrEnergy = 0. 
)
private

register a correlation for the two given types of punch-through particles with a given energy threshold above which we will have full correlation

Definition at line 1325 of file PunchThroughG4Tool.cxx.

1327 {
1328  // find the given pdgs in the registered particle ids
1331 
1332  // if at least one of the given pdgs was not registered yet -> return an error
1333  if ( (location1 == m_particles.end()) || (location2 == m_particles.end()) )
1334  return StatusCode::FAILURE;
1335 
1336  // now look for the correlation histograms
1337  std::stringstream name;
1338  name << "NumExitCorrelations/x_PDG" << std::abs(pdgID1) << "__y_PDG" << std::abs(pdgID2) << "__lowE";
1339  TH2F *hist1_2_lowE = (TH2F*)m_fileLookupTable->Get(name.str().c_str());
1340  name.str("");
1341  name << "NumExitCorrelations/x_PDG" << std::abs(pdgID1) << "__y_PDG" << std::abs(pdgID2) << "__highE";
1342  TH2F *hist1_2_highE = (TH2F*)m_fileLookupTable->Get(name.str().c_str());
1343  name.str("");
1344  name << "NumExitCorrelations/x_PDG" << std::abs(pdgID2) << "__y_PDG" << std::abs(pdgID1) << "__lowE";
1345  TH2F *hist2_1_lowE = (TH2F*)m_fileLookupTable->Get(name.str().c_str());
1346  name.str("");
1347  name << "NumExitCorrelations/x_PDG" << std::abs(pdgID2) << "__y_PDG" << std::abs(pdgID1) << "__highE";
1348  TH2F *hist2_1_highE = (TH2F*)m_fileLookupTable->Get(name.str().c_str());
1349  // check if the histograms exist
1350  if ( (!hist1_2_lowE) || (!hist2_1_lowE) || (!hist1_2_highE) || (!hist2_1_highE) )
1351  {
1352  ATH_MSG_ERROR("[PunchThroughG4Tool] unable to retrieve the correlation data for PDG IDs " << pdgID1 << " and " << pdgID2);
1353  return StatusCode::FAILURE;
1354  }
1355 
1356  // TODO: if only one of the two histograms exists, create the other one
1357  // by mirroring the data
1358  const double lowE = getFloatAfterPatternInStr( hist1_2_lowE->GetTitle(), "elow_");
1359  const double midE = getFloatAfterPatternInStr( hist1_2_lowE->GetTitle(), "ehigh_");
1360  //TODO: check if the same:
1361  // double midE = getFloatAfterPatternInStr( hist1_2_lowE->GetTitle(), "elow_");
1362  const double upperE = getFloatAfterPatternInStr( hist1_2_highE->GetTitle(), "ehigh_");
1363  // now store the correlation either way id1->id2 and id2->id1
1364  m_particles[pdgID1]->setCorrelation(pdgID2, hist2_1_lowE, hist2_1_highE,
1365  minCorrEnergy, fullCorrEnergy,
1366  lowE, midE, upperE);
1367  m_particles[pdgID2]->setCorrelation(pdgID1, hist1_2_lowE, hist1_2_highE,
1368  minCorrEnergy, fullCorrEnergy,
1369  lowE, midE, upperE);
1370  return StatusCode::SUCCESS;
1371 }

◆ registerPunchThroughParticle()

StatusCode PunchThroughG4Tool::registerPunchThroughParticle ( G4ParticleTable &  ptable,
int  pdg,
bool  doAntiparticle = false,
double  minEnergy = 0.,
int  maxNumParticles = -1,
double  numParticlesFactor = 1.,
double  energyFactor = 1.,
double  posAngleFactor = 1.,
double  momAngleFactor = 1. 
)
private

registers a type of punch-through particles which will be simulated

Definition at line 1252 of file PunchThroughG4Tool.cxx.

1255 {
1256  G4ParticleDefinition* secG4Particle = ptable.FindParticle(std::abs(pdg));
1257  double restMass = secG4Particle->GetPDGMass();
1258 
1259  // read in the data needed to construct the distributions for the number of punch-through particles
1260  // (1.) get the distribution function for the number of punch-through particles
1261  std::unique_ptr<PunchThroughPDFCreator> pdf_num(readLookuptablePDF(pdg, m_fileLookupTable, "FREQ_PDG"));
1262  if (!pdf_num ) return StatusCode::FAILURE; // return error if something went wrong
1263 
1264  // (2.) get the PDF for the punch-through energy
1265  std::unique_ptr<PunchThroughPDFCreator> pdf_pca0 (readLookuptablePDF(pdg, m_fileLookupTable, "PCA0_PDG"));
1266  if (!pdf_pca0)
1267  {
1268  return StatusCode::FAILURE; // return error if something went wrong
1269  }
1270 
1271  // (3.) get the PDF for the punch-through particles difference in
1272  // theta compared to the incoming particle
1273  std::unique_ptr<PunchThroughPDFCreator> pdf_pca1 (readLookuptablePDF(pdg, m_fileLookupTable, "PCA1_PDG"));
1274  if (!pdf_pca1)
1275  {
1276  return StatusCode::FAILURE;
1277  }
1278 
1279  // (4.) get the PDF for the punch-through particles difference in
1280  // phi compared to the incoming particle
1281  std::unique_ptr<PunchThroughPDFCreator> pdf_pca2 (readLookuptablePDF(pdg, m_fileLookupTable, "PCA2_PDG"));
1282  if (!pdf_pca2)
1283  {
1284  return StatusCode::FAILURE;
1285  }
1286 
1287  // (5.) get the PDF for the punch-through particle momentum delta theta angle
1288  std::unique_ptr<PunchThroughPDFCreator> pdf_pca3 (readLookuptablePDF(pdg, m_fileLookupTable, "PCA3_PDG"));
1289  if (!pdf_pca3)
1290  {
1291  return StatusCode::FAILURE;
1292  }
1293 
1294  // (6.) get the PDF for the punch-through particle momentum delta phi angle
1295  std::unique_ptr<PunchThroughPDFCreator> pdf_pca4 (readLookuptablePDF(pdg, m_fileLookupTable, "PCA4_PDG"));
1296  if (!pdf_pca4)
1297  {
1298  return StatusCode::FAILURE;
1299  }
1300 
1301  // (7.) now finally store all this in the right std::map
1302  PunchThroughParticle *particle = new PunchThroughParticle(pdg, doAntiparticle);
1303  particle->setNumParticlesPDF(std::move(pdf_num));
1304  particle->setPCA0PDF(std::move(pdf_pca0));
1305  particle->setPCA1PDF(std::move(pdf_pca1));
1306  particle->setPCA2PDF(std::move(pdf_pca2));
1307  particle->setPCA3PDF(std::move(pdf_pca3));
1308  particle->setPCA4PDF(std::move(pdf_pca4));
1309 
1310  // (8.) set some additional particle and simulation properties
1311  minEnergy = ( minEnergy > restMass ) ? minEnergy : restMass;
1312  particle->setMinEnergy(minEnergy);
1313  particle->setMaxNumParticles(maxNumParticles);
1314  particle->setNumParticlesFactor(numParticlesFactor);
1315  particle->setEnergyFactor(energyFactor);
1316  particle->setPosAngleFactor(posAngleFactor);
1317  particle->setMomAngleFactor(momAngleFactor);
1318 
1319  // (9.) insert this PunchThroughParticle instance into the std::map class member
1320  m_particles[pdg] = particle;
1321 
1322  return StatusCode::SUCCESS;
1323 }

Member Data Documentation

◆ m_beamPipe

DoubleProperty PunchThroughG4Tool::m_beamPipe {this, "BeamPipeRadius", 500.}
private

beam pipe radius

Definition at line 220 of file PunchThroughG4Tool.h.

◆ m_correlatedParticle

IntegerArrayProperty PunchThroughG4Tool::m_correlatedParticle {this, "CorrelatedParticle", {}, "holds the pdg of the correlated particle for each given pdg"}
private

Definition at line 203 of file PunchThroughG4Tool.h.

◆ m_doAntiParticles

BooleanArrayProperty PunchThroughG4Tool::m_doAntiParticles {this, "DoAntiParticles", {}, "vector of bools to determine if anti-particles are created for each punch-through particle type"}
private

Definition at line 202 of file PunchThroughG4Tool.h.

◆ m_energyFactor

DoubleArrayProperty PunchThroughG4Tool::m_energyFactor {this, "EnergyFactor", {}, "scale the energy of the punch-through particles"}
private

Definition at line 211 of file PunchThroughG4Tool.h.

◆ m_energyPoints

std::vector<double> PunchThroughG4Tool::m_energyPoints
private

energy and eta points in param

Definition at line 176 of file PunchThroughG4Tool.h.

◆ m_envDefSvc

ServiceHandle<IEnvelopeDefSvc> PunchThroughG4Tool::m_envDefSvc {this, "EnvelopeDefSvc", "AtlasGeometry_EnvelopeDefSvc"}
private

Definition at line 217 of file PunchThroughG4Tool.h.

◆ m_etaPoints

std::vector<double> PunchThroughG4Tool::m_etaPoints
private

Definition at line 177 of file PunchThroughG4Tool.h.

◆ m_fileLookupTable

TFile* PunchThroughG4Tool::m_fileLookupTable {nullptr}
private

ROOT objects.

the punch-through lookup table file

Definition at line 186 of file PunchThroughG4Tool.h.

◆ m_filenameInverseCDF

StringProperty PunchThroughG4Tool::m_filenameInverseCDF {this, "FilenameInverseCdf", "", "holds the filename of inverse quantile transformer config"}
private

Definition at line 195 of file PunchThroughG4Tool.h.

◆ m_filenameInversePCA

StringProperty PunchThroughG4Tool::m_filenameInversePCA {this, "FilenameInversePca", "", "holds the filename of inverse PCA config"}
private

Definition at line 196 of file PunchThroughG4Tool.h.

◆ m_filenameLookupTable

StringProperty PunchThroughG4Tool::m_filenameLookupTable {this, "FilenameLookupTable", "CaloPunchThroughParametrisation.root", "holds the filename of the lookup table"}
private

Definition at line 194 of file PunchThroughG4Tool.h.

◆ m_fullCorrEnergy

DoubleArrayProperty PunchThroughG4Tool::m_fullCorrEnergy {this, "FullCorrelationEnergy", {}, "holds the energy threshold above which a particle correlation is fully developed"}
private

Definition at line 205 of file PunchThroughG4Tool.h.

◆ m_geoIDSvc

ServiceHandle<ISF::IGeoIDSvc> PunchThroughG4Tool::m_geoIDSvc {this, "GeoIDSvc", "ISF::GeoIDSvc"}
private

Definition at line 216 of file PunchThroughG4Tool.h.

◆ m_initiatorsEtaRange

DoubleArrayProperty PunchThroughG4Tool::m_initiatorsEtaRange {this, "InitiatorsEtaRange", {}, "vector of min and max abs eta range to allow punch through initiators"}
private

Definition at line 200 of file PunchThroughG4Tool.h.

◆ m_initiatorsMinEnergy

IntegerArrayProperty PunchThroughG4Tool::m_initiatorsMinEnergy {this, "InitiatorsMinEnergy", {}, "vector of punch-through initiator min energies to create punch through"}
private

Definition at line 199 of file PunchThroughG4Tool.h.

◆ m_inverse_PCA_matrix

std::vector<std::vector<std::vector<double> > > PunchThroughG4Tool::m_inverse_PCA_matrix
private

pca vectors

Definition at line 223 of file PunchThroughG4Tool.h.

◆ m_maxNumParticles

IntegerArrayProperty PunchThroughG4Tool::m_maxNumParticles {this, "MaxNumParticles", {}, "maximum number of punch-through particles for each particle type"}
private

Definition at line 209 of file PunchThroughG4Tool.h.

◆ m_minCorrEnergy

DoubleArrayProperty PunchThroughG4Tool::m_minCorrEnergy {this, "MinCorrelationEnergy", {}, "holds the energy threshold below which no particle correlation is computed"}
private

Definition at line 204 of file PunchThroughG4Tool.h.

◆ m_minEnergy

DoubleArrayProperty PunchThroughG4Tool::m_minEnergy {this, "MinEnergy", {}, "punch-through particles minimum energies"}
private

Definition at line 208 of file PunchThroughG4Tool.h.

◆ m_momAngleFactor

DoubleArrayProperty PunchThroughG4Tool::m_momAngleFactor {this, "ScaleMomDeflectionAngles", {}, "tuning parameter to scale the momentum deflection angles"}
private

Definition at line 207 of file PunchThroughG4Tool.h.

◆ m_numParticlesFactor

DoubleArrayProperty PunchThroughG4Tool::m_numParticlesFactor {this, "NumParticlesFactor", {}, "scale the number of punch-through particles"}
private

Definition at line 210 of file PunchThroughG4Tool.h.

◆ m_particles

std::map<int, PunchThroughParticle*> PunchThroughG4Tool::m_particles
private

needed to initially create punch-through particles with the right distributions

store all punch-through information for each particle id

Definition at line 189 of file PunchThroughG4Tool.h.

◆ m_PCA_means

std::vector<std::vector<double> > PunchThroughG4Tool::m_PCA_means
private

Definition at line 224 of file PunchThroughG4Tool.h.

◆ m_pdgInitiators

IntegerArrayProperty PunchThroughG4Tool::m_pdgInitiators {this, "PunchThroughInitiators", {}, "vector of punch-through initiator pgds"}
private

Definition at line 198 of file PunchThroughG4Tool.h.

◆ m_posAngleFactor

DoubleArrayProperty PunchThroughG4Tool::m_posAngleFactor {this, "ScalePosDeflectionAngles", {}, "tuning parameter to scale the position deflection angles"}
private

Definition at line 206 of file PunchThroughG4Tool.h.

◆ m_punchThroughParticles

IntegerArrayProperty PunchThroughG4Tool::m_punchThroughParticles {this, "PunchThroughParticles", {}, "vector of pdgs of the particles produced in punch-throughs"}
private

Definition at line 201 of file PunchThroughG4Tool.h.

◆ m_R1

double PunchThroughG4Tool::m_R1 {0.}
private

calo-MS borders

Definition at line 180 of file PunchThroughG4Tool.h.

◆ m_R2

double PunchThroughG4Tool::m_R2 {0.}
private

Definition at line 181 of file PunchThroughG4Tool.h.

◆ m_variable0_inverse_cdf

std::vector<std::map<double, double> > PunchThroughG4Tool::m_variable0_inverse_cdf
private

(vector of map) for CDF mappings

Definition at line 231 of file PunchThroughG4Tool.h.

◆ m_variable1_inverse_cdf

std::vector<std::map<double, double> > PunchThroughG4Tool::m_variable1_inverse_cdf
private

Definition at line 232 of file PunchThroughG4Tool.h.

◆ m_variable2_inverse_cdf

std::vector<std::map<double, double> > PunchThroughG4Tool::m_variable2_inverse_cdf
private

Definition at line 233 of file PunchThroughG4Tool.h.

◆ m_variable3_inverse_cdf

std::vector<std::map<double, double> > PunchThroughG4Tool::m_variable3_inverse_cdf
private

Definition at line 234 of file PunchThroughG4Tool.h.

◆ m_variable4_inverse_cdf

std::vector<std::map<double, double> > PunchThroughG4Tool::m_variable4_inverse_cdf
private

Definition at line 235 of file PunchThroughG4Tool.h.

◆ m_xml_info_cdf

std::vector<std::map<std::string, std::string> > PunchThroughG4Tool::m_xml_info_cdf
private

Definition at line 228 of file PunchThroughG4Tool.h.

◆ m_xml_info_pca

std::vector<std::map<std::string, std::string> > PunchThroughG4Tool::m_xml_info_pca
private

infoMaps

Definition at line 227 of file PunchThroughG4Tool.h.

◆ m_z1

double PunchThroughG4Tool::m_z1 {0.}
private

Definition at line 182 of file PunchThroughG4Tool.h.

◆ m_z2

double PunchThroughG4Tool::m_z2 {0.}
private

Definition at line 183 of file PunchThroughG4Tool.h.


The documentation for this class was generated from the following files:
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
xmlChar
unsigned char xmlChar
Definition: TGoodRunsListWriter.h:28
PunchThroughParticle
Definition: G4Atlas/G4AtlasTools/src/PunchThroughParticle.h:31
mergePhysValFiles.pattern
pattern
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:25
beamspotman.r
def r
Definition: beamspotman.py:674
PunchThroughG4Tool::m_correlatedParticle
IntegerArrayProperty m_correlatedParticle
Definition: PunchThroughG4Tool.h:203
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:79
get_generator_info.result
result
Definition: get_generator_info.py:21
checkCoolLatestUpdate.variables
variables
Definition: checkCoolLatestUpdate.py:12
AddEmptyComponent.histName
string histName
Definition: AddEmptyComponent.py:64
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
PunchThroughG4Tool::m_filenameLookupTable
StringProperty m_filenameLookupTable
Definition: PunchThroughG4Tool.h:194
Base_Fragment.mass
mass
Definition: Sherpa_i/share/common/Base_Fragment.py:59
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:83
PunchThroughG4Tool::initializeInverseCDF
StatusCode initializeInverseCDF(const std::string &quantileTransformerConfigFile)
Definition: PunchThroughG4Tool.cxx:1023
xAOD::deltaPhi
setSAddress setEtaMS setDirPhiMS setDirZMS setBarrelRadius setEndcapAlpha setEndcapRadius setInterceptInner setEtaMap setEtaBin setIsTgcFailure setDeltaPt deltaPhi
Definition: L2StandAloneMuon_v1.cxx:161
PunchThroughG4Tool::m_fullCorrEnergy
DoubleArrayProperty m_fullCorrEnergy
Definition: PunchThroughG4Tool.h:205
plotmaker.hist
hist
Definition: plotmaker.py:148
PunchThroughG4Tool::inversePCA
std::vector< double > inversePCA(int pcaCdfIterator, std::vector< double > &variables) const
Definition: PunchThroughG4Tool.cxx:951
AtlasDetDescr::AtlasRegion
AtlasRegion
Definition: AtlasRegion.h:21
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:75
MuonGM::round
float round(const float toRound, const unsigned int decimals)
Definition: Mdt.cxx:27
PunchThroughG4Tool::initializeRegisterPunchThroughParticles
StatusCode initializeRegisterPunchThroughParticles()
initialize register all the punch-through particles which will be simulated
Definition: PunchThroughG4Tool.cxx:73
M_PI
#define M_PI
Definition: ActiveFraction.h:11
PunchThroughG4Tool::interpolateEnergy
double interpolateEnergy(const double &energy, CLHEP::HepRandomEngine *rndmEngine) const
Definition: PunchThroughG4Tool.cxx:1124
PunchThroughG4Tool::checkCaloMSBoundaries
StatusCode checkCaloMSBoundaries(const std::vector< std::pair< double, double >> *rzMS, const std::vector< std::pair< double, double >> *rzCalo)
Check calo-MS boundaries.
Definition: PunchThroughG4Tool.cxx:205
PunchThroughG4Tool::dotProduct
static std::vector< double > dotProduct(const std::vector< std::vector< double >> &m, const std::vector< double > &v)
Definition: PunchThroughG4Tool.cxx:841
upper
int upper(int c)
Definition: LArBadChannelParser.cxx:49
PunchThroughG4Tool::registerPunchThroughParticle
StatusCode registerPunchThroughParticle(G4ParticleTable &ptable, int pdg, bool doAntiparticle=false, double minEnergy=0., int maxNumParticles=-1, double numParticlesFactor=1., double energyFactor=1., double posAngleFactor=1., double momAngleFactor=1.)
registers a type of punch-through particles which will be simulated
Definition: PunchThroughG4Tool.cxx:1252
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
PrintTrkAnaSummary.dirName
dirName
Definition: PrintTrkAnaSummary.py:275
PunchThroughG4Tool::m_numParticlesFactor
DoubleArrayProperty m_numParticlesFactor
Definition: PunchThroughG4Tool.h:210
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
PunchThroughG4Tool::m_inverse_PCA_matrix
std::vector< std::vector< std::vector< double > > > m_inverse_PCA_matrix
pca vectors
Definition: PunchThroughG4Tool.h:223
python.TrigEgammaMonitorHelper.TH2F
def TH2F(name, title, nxbins, bins_par2, bins_par3, bins_par4, bins_par5=None, bins_par6=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:45
PunchThroughG4Tool::m_R1
double m_R1
calo-MS borders
Definition: PunchThroughG4Tool.h:180
x
#define x
PunchThroughG4Tool::getFloatAfterPatternInStr
double getFloatAfterPatternInStr(const char *str, const char *pattern)
get the floating point number in a string, after the given pattern
Definition: PunchThroughG4Tool.cxx:1440
PunchThroughG4Tool::getCorrelatedParticles
int getCorrelatedParticles(const G4Track &g4PrimaryTrack, std::vector< std::map< std::string, double >> &secKinematicsMapVect, int pdg, int corrParticles, CLHEP::HepRandomEngine *rndmEngine, double interpEnergy, double interpEta)
get the right number of particles for the given pdg while considering the correlation to an other par...
Definition: PunchThroughG4Tool.cxx:505
python.CaloAddPedShiftConfig.type
type
Definition: CaloAddPedShiftConfig.py:42
PunchThroughG4Tool::m_geoIDSvc
ServiceHandle< ISF::IGeoIDSvc > m_geoIDSvc
Definition: PunchThroughG4Tool.h:216
drawFromPickle.atan
atan
Definition: drawFromPickle.py:36
PunchThroughG4Tool::normal_cdf
static double normal_cdf(double x)
Definition: PunchThroughG4Tool.cxx:837
PunchThroughG4Tool::m_variable0_inverse_cdf
std::vector< std::map< double, double > > m_variable0_inverse_cdf
(vector of map) for CDF mappings
Definition: PunchThroughG4Tool.h:231
PunchThroughG4Tool::m_z1
double m_z1
Definition: PunchThroughG4Tool.h:182
PunchThroughG4Tool::registerCorrelation
StatusCode registerCorrelation(int pdgID1, int pdgID2, double minCorrEnergy=0., double fullCorrEnergy=0.)
register a correlation for the two given types of punch-through particles with a given energy thresho...
Definition: PunchThroughG4Tool.cxx:1325
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
PunchThroughG4Tool::m_doAntiParticles
BooleanArrayProperty m_doAntiParticles
Definition: PunchThroughG4Tool.h:202
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
convertTimingResiduals.sum
sum
Definition: convertTimingResiduals.py:55
LArG4FSStartPointFilter.rand
rand
Definition: LArG4FSStartPointFilter.py:80
ParticleGun_FastCalo_ChargeFlip_Config.energy
energy
Definition: ParticleGun_FastCalo_ChargeFlip_Config.py:78
lumiFormat.i
int i
Definition: lumiFormat.py:85
z
#define z
PunchThroughG4Tool::createSecondaryTrack
G4Track * createSecondaryTrack(G4ParticleTable &ptable, G4FastStep &fastStep, double currentTime, int secondarySignedPDG, double energy, double theta, double phi, double momTheta, double momPhi, const std::vector< double > &caloMSVars)
create secondary track for each given the kinematics
Definition: PunchThroughG4Tool.cxx:808
PunchThroughG4Tool::punchTroughPosPropagator
G4ThreeVector punchTroughPosPropagator(double theta, double phi, double R1, double R2, double z1, double z2) const
get particle through the calorimeter
Definition: PunchThroughG4Tool.cxx:1459
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
python.CaloAddPedShiftConfig.str
str
Definition: CaloAddPedShiftConfig.py:42
Amg::transform
Amg::Vector3D transform(Amg::Vector3D &v, Amg::Transform3D &tr)
Transform a point from a Trasformation3D.
Definition: GeoPrimitivesHelpers.h:156
ParticleGun_EoverP_Config.pid
pid
Definition: ParticleGun_EoverP_Config.py:62
test_pyathena.parent
parent
Definition: test_pyathena.py:15
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
PunchThroughG4Tool::m_initiatorsEtaRange
DoubleArrayProperty m_initiatorsEtaRange
Definition: PunchThroughG4Tool.h:200
python.selection.variable
variable
Definition: selection.py:33
PunchThroughG4Tool::m_filenameInversePCA
StringProperty m_filenameInversePCA
Definition: PunchThroughG4Tool.h:196
PunchThroughG4Tool::m_xml_info_cdf
std::vector< std::map< std::string, std::string > > m_xml_info_cdf
Definition: PunchThroughG4Tool.h:228
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
PunchThroughG4Tool::m_beamPipe
DoubleProperty m_beamPipe
beam pipe radius
Definition: PunchThroughG4Tool.h:220
WriteCalibToCool.swap
swap
Definition: WriteCalibToCool.py:94
PunchThroughG4Tool::initializeInversePCA
StatusCode initializeInversePCA(const std::string &inversePCAConfigFile)
Definition: PunchThroughG4Tool.cxx:960
beamspotman.dir
string dir
Definition: beamspotman.py:621
CxxUtils::atof
double atof(std::string_view str)
Converts a string into a double / float.
Definition: Control/CxxUtils/Root/StringUtils.cxx:91
trigbs_pickEvents.num
num
Definition: trigbs_pickEvents.py:76
PunchThroughG4Tool::m_initiatorsMinEnergy
IntegerArrayProperty m_initiatorsMinEnergy
Definition: PunchThroughG4Tool.h:199
PunchThroughG4Tool::m_PCA_means
std::vector< std::vector< double > > m_PCA_means
Definition: PunchThroughG4Tool.h:224
node::name
void name(const std::string &n)
Definition: node.h:37
CaloCellTimeCorrFiller.folderName
string folderName
Definition: CaloCellTimeCorrFiller.py:19
merge_scale_histograms.doc
string doc
Definition: merge_scale_histograms.py:9
PunchThroughG4Tool::m_variable3_inverse_cdf
std::vector< std::map< double, double > > m_variable3_inverse_cdf
Definition: PunchThroughG4Tool.h:234
PunchThroughG4Tool::getInfoMap
std::vector< std::map< std::string, std::string > > getInfoMap(const std::string &mainNode, const std::string &xmlFilePath)
Definition: PunchThroughG4Tool.cxx:905
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
PunchThroughG4Tool::m_minEnergy
DoubleArrayProperty m_minEnergy
Definition: PunchThroughG4Tool.h:208
ActsTrk::to_string
std::string to_string(const DetectorType &type)
Definition: GeometryDefs.h:34
L1CaloPhase1Monitoring.propName
propName
Definition: L1CaloPhase1Monitoring.py:558
PunchThroughG4Tool::interpolateEta
double interpolateEta(const double &eta, CLHEP::HepRandomEngine *rndmEngine) const
Definition: PunchThroughG4Tool.cxx:1185
FakeBkgTools::maxParticles
constexpr uint8_t maxParticles()
Definition: FakeBkgInternals.h:93
PunchThroughG4Tool::inverseCdfTransform
static double inverseCdfTransform(double variable, const std::map< double, double > &inverse_cdf_map)
Definition: PunchThroughG4Tool.cxx:1109
PunchThroughG4Tool::initializeRegisterCorrelations
StatusCode initializeRegisterCorrelations()
initialize register all correlations between particles
Definition: PunchThroughG4Tool.cxx:162
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
PunchThroughG4Tool::m_pdgInitiators
IntegerArrayProperty m_pdgInitiators
Definition: PunchThroughG4Tool.h:198
PunchThroughG4Tool::checkEnergySumFromSecondaries
std::vector< std::map< std::string, double > > checkEnergySumFromSecondaries(double mainEnergyInit, std::vector< std::map< std::string, double >> &secKinematicsMapVect)
check the energies satisfying energy condition
Definition: PunchThroughG4Tool.cxx:712
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:16
PathResolverFindCalibFile
std::string PathResolverFindCalibFile(const std::string &logical_file_name)
Definition: PathResolver.cxx:283
AtlasDetDescr::fAtlasCalo
@ fAtlasCalo
Definition: AtlasRegion.h:29
python.PyAthena.v
v
Definition: PyAthena.py:154
PunchThroughG4Tool::m_minCorrEnergy
DoubleArrayProperty m_minCorrEnergy
Definition: PunchThroughG4Tool.h:204
PunchThroughG4Tool::m_filenameInverseCDF
StringProperty m_filenameInverseCDF
Definition: PunchThroughG4Tool.h:195
PunchThroughG4Tool::getOneParticleKinematics
std::map< std::string, double > getOneParticleKinematics(CLHEP::HepRandomEngine *rndmEngine, int secondaryPDG, float initParticleTheta, float initParticlePhi, double interpEnergy, double interpEta) const
create exactly one punch-through particle with the given pdg and the given max energy
Definition: PunchThroughG4Tool.cxx:579
LArCellBinning.etaMin
etaMin
Definition: LArCellBinning.py:84
PunchThroughG4Tool::m_xml_info_pca
std::vector< std::map< std::string, std::string > > m_xml_info_pca
infoMaps
Definition: PunchThroughG4Tool.h:227
PunchThroughG4Tool::m_energyFactor
DoubleArrayProperty m_energyFactor
Definition: PunchThroughG4Tool.h:211
y
#define y
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
PunchThroughG4Tool::m_variable2_inverse_cdf
std::vector< std::map< double, double > > m_variable2_inverse_cdf
Definition: PunchThroughG4Tool.h:233
PunchThroughG4Tool::m_punchThroughParticles
IntegerArrayProperty m_punchThroughParticles
Definition: PunchThroughG4Tool.h:201
std::sort
void sort(typename std::reverse_iterator< DataModel_detail::iterator< DVL > > beg, typename std::reverse_iterator< DataModel_detail::iterator< DVL > > end, const Compare &comp)
Specialization of sort for DataVector/List.
Definition: DVL_algorithms.h:623
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
PunchThroughG4Tool::passedParamIterator
int passedParamIterator(int pid, double eta, const std::vector< std::map< std::string, std::string >> &mapvect) const
Definition: PunchThroughG4Tool.cxx:876
ref
const boost::regex ref(r_ef)
PunchThroughG4Tool::getVariableCDFmappings
static std::map< double, double > getVariableCDFmappings(xmlNodePtr &nodeParent)
Definition: PunchThroughG4Tool.cxx:1088
PunchThroughG4Tool::m_maxNumParticles
IntegerArrayProperty m_maxNumParticles
Definition: PunchThroughG4Tool.h:209
PunchThroughG4Tool::m_fileLookupTable
TFile * m_fileLookupTable
ROOT objects.
Definition: PunchThroughG4Tool.h:186
PunchThroughG4Tool::m_envDefSvc
ServiceHandle< IEnvelopeDefSvc > m_envDefSvc
Definition: PunchThroughG4Tool.h:217
PowhegPythia8EvtGen_jetjet.pdf
pdf
Definition: PowhegPythia8EvtGen_jetjet.py:4
TauGNNUtils::Variables::absEta
bool absEta(const xAOD::TauJet &tau, double &out)
Definition: TauGNNUtils.cxx:245
PunchThroughG4Tool::readLookuptablePDF
std::unique_ptr< PunchThroughPDFCreator > readLookuptablePDF(int pdgID, TFile *fileLookupTable, const std::string &folderName)
reads out the lookuptable for the given type of particle
Definition: PunchThroughG4Tool.cxx:1373
physics_parameters.parameters
parameters
Definition: physics_parameters.py:144
str
Definition: BTagTrackIpAccessor.cxx:11
PunchThroughG4Tool::m_momAngleFactor
DoubleArrayProperty m_momAngleFactor
Definition: PunchThroughG4Tool.h:207
PunchThroughG4Tool::m_R2
double m_R2
Definition: PunchThroughG4Tool.h:181
PunchThroughG4Tool::m_energyPoints
std::vector< double > m_energyPoints
energy and eta points in param
Definition: PunchThroughG4Tool.h:176
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
PunchThroughG4Tool::getAllParticles
int getAllParticles(const G4Track &g4PrimaryTrack, std::vector< std::map< std::string, double >> &secKinematicsMapVect, CLHEP::HepRandomEngine *rndmEngine, int pdg, double interpEnergy, double interpEta, int numParticles=-1)
create the right number of punch-through particles for the given pdg and return the number of particl...
Definition: PunchThroughG4Tool.cxx:448
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
python.compressB64.c
def c
Definition: compressB64.py:93
PunchThroughG4Tool::m_posAngleFactor
DoubleArrayProperty m_posAngleFactor
Definition: PunchThroughG4Tool.h:206
PunchThroughG4Tool::m_variable1_inverse_cdf
std::vector< std::map< double, double > > m_variable1_inverse_cdf
Definition: PunchThroughG4Tool.h:232
PunchThroughG4Tool::m_variable4_inverse_cdf
std::vector< std::map< double, double > > m_variable4_inverse_cdf
Definition: PunchThroughG4Tool.h:235
PunchThroughG4Tool::m_particles
std::map< int, PunchThroughParticle * > m_particles
needed to initially create punch-through particles with the right distributions
Definition: PunchThroughG4Tool.h:189
PunchThroughG4Tool::m_etaPoints
std::vector< double > m_etaPoints
Definition: PunchThroughG4Tool.h:177
node
Definition: node.h:21
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106
mapkey::key
key
Definition: TElectronEfficiencyCorrectionTool.cxx:37
PunchThroughG4Tool::m_z2
double m_z2
Definition: PunchThroughG4Tool.h:183