ATLAS Offline Software
Loading...
Searching...
No Matches
ISF::PunchThroughTool Class Reference

#include <PunchThroughTool.h>

Inheritance diagram for ISF::PunchThroughTool:

Classes

struct  InfoMap

Public Member Functions

 PunchThroughTool (const std::string &, const std::string &, const IInterface *)
 Constructor.
virtual ~PunchThroughTool ()=default
 Destructor.
virtual StatusCode initialize ()
 AlgTool initialize method.
virtual StatusCode finalize ()
 AlgTool finalize method.
const ISF::ISFParticleVectorcomputePunchThroughParticles (const ISF::ISFParticle &isfp, const TFCSSimulationState &simulstate, CLHEP::HepRandomEngine *rndmEngine) const
 interface function: fill a vector with the punch-through particles

Private Member Functions

StatusCode registerParticle (int pdgID, 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
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
std::unique_ptr< ISF::PDFcreatorreadLookuptablePDF (int pdgID, const std::string &folderName)
 reads out the lookuptable for the given type of particle
int getAllParticles (const ISF::ISFParticle &isfp, ISFParticleVector &isfpCont, CLHEP::HepRandomEngine *rndmEngine, int pdg, double interpEnergy, double interpEta, int numParticles=-1) const
 create the right number of punch-through particles for the given pdg and return the number of particles which was created.
int getCorrelatedParticles (const ISF::ISFParticle &isfp, ISFParticleVector &isfpCont, int doPdg, int corrParticles, CLHEP::HepRandomEngine *rndmEngine, double interpEnergy, double interpEta) const
 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
ISF::ISFParticlegetOneParticle (const ISF::ISFParticle &isfp, int pdg, CLHEP::HepRandomEngine *rndmEngine, double interpEnergy, double interpEta) const
 create exactly one punch-through particle with the given pdg and the given max energy
ISF::ISFParticlecreateExitPs (const ISF::ISFParticle &isfp, int PDGcode, double energy, double theta, double phi, double momTheta, double momPhi) const
 create a ISF Particle state at the MS entrace containing a particle with the given properties
double getFloatAfterPatternInStr (const char *str, const char *pattern)
 get the floating point number in a string, after the given pattern
Amg::Vector3D propagator (double theta, double phi) const
 get particle through the calorimeter
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< InfoMapgetInfoMap (const std::string &mainNode, const XMLCoreNode &doc)
int passedParamIterator (int pid, double eta, const std::vector< InfoMap > &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 (const XMLCoreNode *node)

Private Attributes

std::vector< double > m_energyPoints
 energy and eta points in param
std::vector< double > m_etaPoints
double m_R1 {0.}
 calo-MS borders
double m_R2 {0.}
double m_z1 {0.}
double m_z2 {0.}
const HepPDT::ParticleDataTable * m_particleDataTable {nullptr}
 ParticleDataTable needed to get connection pdg_code <-> charge.
TFile * m_fileLookupTable {nullptr}
 ROOT objects.
std::map< int, PunchThroughParticle * > m_particles
 needed to create punch-through particles with the right distributions
StringProperty m_filenameLookupTable {this, "FilenameLookupTable", "CaloPunchThroughParametrisation.root", "holds the filename of the lookup table"}
 Properties.
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"}
PublicToolHandle< IPunchThroughClassifierm_punchThroughClassifier {this, "PunchThroughClassifier", "ISF_PunchThroughClassifier", ""}
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< IPartPropSvc > m_particlePropSvc {this, "PartPropSvc", "PartPropSvc", "particle properties svc"}
ServiceHandle< IGeoIDSvcm_geoIDSvc {this, "GeoIDSvc", "ISF::GeoIDSvc"}
ServiceHandle< IEnvelopeDefSvcm_envDefSvc {this, "EnvelopeDefSvc", "AtlasGeometry_EnvelopeDefSvc"}
DoubleProperty m_beamPipe {this, "BeamPipeRadius", 500.}
 beam pipe radius
std::vector< std::vector< std::vector< double > > > m_inverse_PCA_matrix
 pca vectors
std::vector< std::vector< double > > m_PCA_means
std::vector< InfoMapm_xml_info_pca
 infoMaps
std::vector< InfoMapm_xml_info_cdf
std::vector< std::map< double, double > > m_variable0_inverse_cdf
 (vector of map) for CDF mappings
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

Definition at line 49 of file PunchThroughTool.h.

Constructor & Destructor Documentation

◆ PunchThroughTool()

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

Constructor.

Definition at line 65 of file PunchThroughTool.cxx.

68: base_class(type, name, parent)
69{
70}

◆ ~PunchThroughTool()

virtual ISF::PunchThroughTool::~PunchThroughTool ( )
virtualdefault

Destructor.

Member Function Documentation

◆ computePunchThroughParticles()

const ISF::ISFParticleVector * ISF::PunchThroughTool::computePunchThroughParticles ( const ISF::ISFParticle & isfp,
const TFCSSimulationState & simulstate,
CLHEP::HepRandomEngine * rndmEngine ) const

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

Definition at line 291 of file PunchThroughTool.cxx.

292{
293 ATH_MSG_DEBUG( "[ punchthrough ] starting punch-through simulation");
294
295 // create output particle collection
296 auto isfpCont = std::make_unique<ISF::ISFParticleVector>();
297
298 ATH_MSG_VERBOSE("[ punchthrough ] position of the input particle: r"<<isfp.position().perp()<<" z= "<<isfp.position().z() );
299
300 //check if it points to the calorimeter - if not, don't simulate
301
302 if ( m_geoIDSvc->identifyNextGeoID(isfp) != AtlasDetDescr::fAtlasCalo)
303 {
304 ATH_MSG_VERBOSE ("[ GeoIDSvc ] input particle doesn't point to calorimeter"<< "Next GeoID: "<<m_geoIDSvc->identifyNextGeoID(isfp) );
305 return nullptr;
306 }
307
308
309 // check if the particle's pdg is registered as a punch-through-causing type
310 {
311 std::vector<int>::const_iterator pdgIt = m_pdgInitiators.begin();
312 std::vector<int>::const_iterator pdgItEnd = m_pdgInitiators.end();
313
314 std::vector<int>::const_iterator minEnergyIt = m_initiatorsMinEnergy.begin();
315 // loop over all known punch-through initiators
316 for ( ; pdgIt != pdgItEnd; ++pdgIt, ++minEnergyIt)
317 {
318 if (std::abs(isfp.pdgCode()) == *pdgIt){
319 if(std::sqrt( isfp.momentum().mag2() + isfp.mass()*isfp.mass() ) < *minEnergyIt){
320 ATH_MSG_DEBUG("[ punchthrough ] particle does not meet initiator min energy requirement. Dropping it in the calo.");
321 return nullptr;
322 }
323 break;
324 }
325 }
326
327 // particle will not cause punch-through -> bail out
328 if (pdgIt == pdgItEnd)
329 {
330 ATH_MSG_DEBUG("[ punchthrough ] particle is not registered as punch-through initiator. Dropping it in the calo.");
331 return nullptr;
332 }
333 }
334
335 if(isfp.position().eta() < m_initiatorsEtaRange.value().at(0) || isfp.position().eta() > m_initiatorsEtaRange.value().at(1) ){
336 ATH_MSG_DEBUG("[ punchthrough ] particle does not meet initiator eta range requirement. Dropping it in the calo.");
337 return nullptr;
338 }
339
340 //Calculate probability of punch through using punchThroughClassifier
341 double punchThroughProbability = m_punchThroughClassifier->computePunchThroughProbability(isfp, simulstate);
342
343 //Draw random number to compare to probability
344 double punchThroughClassifierRand = CLHEP::RandFlat::shoot(rndmEngine);
345
346 ATH_MSG_DEBUG("[ punchthrough ] punchThroughProbability output: " << punchThroughProbability << " RandFlat: " << punchThroughClassifierRand );
347
348 //If probability < random number then don't simulate punch through
349 if( punchThroughClassifierRand > punchThroughProbability){
350 ATH_MSG_DEBUG("[ punchthrough ] particle not classified to create punch through. Dropping it in the calo.");
351 return nullptr;
352 }
353
354 //if initial particle is on ID surface, points to the calorimeter, is a punch-through initiator, meets initiator min enery and eta range
355
356 // this is the place where the magic is done:
357 // test for each registered punch-through pdg if a punch-through
358 // occures and create these particles
359 // -> therefore loop over all registered pdg ids
360 // to keep track of the correlated particles which were already simulated:
361 // first int is pdg, second int is number of particles created
362
363 // calculate incoming energy and eta
364 const double initEnergy = std::sqrt( isfp.momentum().mag2() + isfp.mass()*isfp.mass() );
365 const double initEta = isfp.position().eta();
366
367 // interpolate energy and eta
368 const double interpEnergy = interpolateEnergy(initEnergy, rndmEngine);
369 const double interpEta = interpolateEta(initEta, rndmEngine);
370
371 std::map<int, int> corrPdgNumDone;
372
373 int maxTries = 10;
374 int nTries = 0;
375
376 // loop over all particle pdgs
377 while(isfpCont->empty() && nTries < maxTries) { //ensure we always create at least one punch through particle, maxTries to catch very rare cases
378
379 for (const auto& currentParticle : m_particles)
380 {
381 // the pdg that is currently treated
382 int doPdg = currentParticle.first;
383 // get the current particle's correlated pdg
384 int corrPdg = currentParticle.second->getCorrelatedPdg();
385
386 // if there is a correlated particle type to this one
387 if (corrPdg)
388 {
389 // find out if the current pdg was already simulated
390 std::map<int,int>::iterator pos = corrPdgNumDone.find(doPdg);
391 // if the current pdg was not simulated yet, find out if
392 // it's correlated one was simulated
393 if ( pos == corrPdgNumDone.end() ) pos = corrPdgNumDone.find(corrPdg);
394
395 // neither this nor the correlated particle type was simulated
396 // so far:
397 if ( pos == corrPdgNumDone.end() )
398 {
399 // -> roll a dice if we create this particle or its correlated one
400 if ( CLHEP::RandFlat::shoot(rndmEngine) > 0.5 ) doPdg = corrPdg;
401 // now create the particles with the given pdg and note how many
402 // particles of this pdg are created
403 corrPdgNumDone[doPdg] = getAllParticles(isfp, *isfpCont, rndmEngine, doPdg, interpEnergy, interpEta);
404 }
405
406 // one of the two correlated particle types was already simulated
407 // 'pos' points to the already simulated one
408 else
409 {
410 // get the pdg of the already simulated particle and the number
411 // of these particles that were created
412 const int donePdg = pos->first;
413 const int doneNumPart = pos->second;
414 // set the pdg of the particle type that will be done
415 if (donePdg == doPdg) doPdg = corrPdg;
416
417 // now create the correlated particles
418 getCorrelatedParticles(isfp, *isfpCont, doPdg, doneNumPart, rndmEngine, interpEnergy, interpEta);
419 // note: no need to take note, that this particle type is now simulated,
420 // since this is the second of two correlated particles, which is
421 // simulated and we do not have correlations of more than two particles.
422 }
423
424 // if no correlation for this particle
425 // -> directly create all particles with the current pdg
426 }
427 else getAllParticles(isfp, *isfpCont, rndmEngine, doPdg, interpEnergy, interpEta);
428
429 } // for-loop over all particle pdgs
430
431 nTries++;
432 }
433
434 if (!isfpCont->empty()) ATH_MSG_DEBUG( "[ punchthrough ] returning ISFparticle vector , size: "<<isfpCont->size() );
435
436 for (ISF::ISFParticle *particle : *isfpCont) {
437 ATH_MSG_DEBUG("codes of produced punch through particle: pdg = "<< particle->pdgCode());
438 Amg::Vector3D position = particle->position();
439 ATH_MSG_DEBUG("position of produced punch-through particle: x = "<< position.x() <<" y = "<< position.y() <<" z = "<< position.z());
440 Amg::Vector3D momentum = particle->momentum();
441 ATH_MSG_DEBUG("momentum of produced punch-through particle: px = "<< momentum.x() <<" py = "<< momentum.x() <<" pz = "<< momentum.x() <<" e = "<< particle->ekin() << " mass = " << particle->mass());
442 }
443
444 return isfpCont.release();
445}
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_DEBUG(x)
const Amg::Vector3D & momentum() const
The current momentum vector of the ISFParticle.
const Amg::Vector3D & position() const
The current position of the ISFParticle.
int pdgCode() const
PDG value.
double mass() const
mass of the particle
PublicToolHandle< IPunchThroughClassifier > m_punchThroughClassifier
int getAllParticles(const ISF::ISFParticle &isfp, ISFParticleVector &isfpCont, CLHEP::HepRandomEngine *rndmEngine, int pdg, double interpEnergy, double interpEta, int numParticles=-1) const
create the right number of punch-through particles for the given pdg and return the number of particl...
DoubleArrayProperty m_initiatorsEtaRange
std::map< int, PunchThroughParticle * > m_particles
needed to create punch-through particles with the right distributions
IntegerArrayProperty m_initiatorsMinEnergy
ServiceHandle< IGeoIDSvc > m_geoIDSvc
int getCorrelatedParticles(const ISF::ISFParticle &isfp, ISFParticleVector &isfpCont, int doPdg, int corrParticles, CLHEP::HepRandomEngine *rndmEngine, double interpEnergy, double interpEta) const
get the right number of particles for the given pdg while considering the correlation to an other par...
IntegerArrayProperty m_pdgInitiators
double interpolateEta(const double &eta, CLHEP::HepRandomEngine *rndmEngine) const
double interpolateEnergy(const double &energy, CLHEP::HepRandomEngine *rndmEngine) const
Eigen::Matrix< double, 3, 1 > Vector3D
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses

◆ createExitPs()

ISF::ISFParticle * ISF::PunchThroughTool::createExitPs ( const ISF::ISFParticle & isfp,
int PDGcode,
double energy,
double theta,
double phi,
double momTheta,
double momPhi ) const
private

create a ISF Particle state at the MS entrace containing a particle with the given properties

@TODO: fix

Definition at line 1259 of file PunchThroughTool.cxx.

1261{
1262 // the intersection point with Calo-MS surface
1263
1265
1266 // set up the real punch-through particle at this position
1267 // set up the momentum vector of this particle as a GlobalMomentum
1268 // by using the given energy and mass of the particle and also using
1269 // the given theta and phi
1270
1272 double mass = m_particleDataTable->particle(std::abs(pdg))->mass();
1273 Amg::setRThetaPhi( mom, std::sqrt(energy*energy - mass*mass), momTheta, momPhi);
1274 ATH_MSG_DEBUG("setRThetaPhi pre input parameters: energy = "<< energy <<" mass = "<< mass);
1275 ATH_MSG_DEBUG("setRThetaPhi input parameters: std::sqrt(energy*energy - mass*mass) = "<< std::sqrt(energy*energy - mass*mass) <<" momTheta = "<< momTheta <<" momPhi = "<< momPhi);
1276
1277
1278 double charge = m_particleDataTable->particle(std::abs(pdg))->charge();
1279 // since the PDT table only has abs(PID) values for the charge
1280 charge *= (pdg > 0.) ? 1. : -1.;
1281
1282 const double pTime = 0;
1283 const int status = 1 + HepMC::SIM_STATUS_THRESHOLD;
1284 const int id = HepMC::UNDEFINED_ID;
1285 // NB we are not considering the possibility that the punch-through
1286 // particle is the incoming particle having survived an interaction.
1287 ISF::ISFParticle* finalPar = new ISF::ISFParticle ( pos, mom, mass, charge, pdg, status, pTime, isfp, id);
1289
1290 // return the punch-through particle
1291 return finalPar;
1292}
Scalar phi() const
phi method
Scalar theta() const
theta method
double charge(const T &p)
Definition AtlasPID.h:997
void setNextGeoID(AtlasDetDescr::AtlasRegion geoID)
register the next AtlasDetDescr::AtlasRegion
Amg::Vector3D propagator(double theta, double phi) const
get particle through the calorimeter
const HepPDT::ParticleDataTable * m_particleDataTable
ParticleDataTable needed to get connection pdg_code <-> charge.
void setRThetaPhi(Amg::Vector3D &v, double r, double theta, double phi)
sets radius, the theta and phi angle of a vector.
constexpr int UNDEFINED_ID
constexpr int SIM_STATUS_THRESHOLD
Constant definiting the status threshold for simulated particles, eg. can be used to separate generat...
status
Definition merge.py:16

◆ dotProduct()

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

Definition at line 732 of file PunchThroughTool.cxx.

733{
734 std::vector<double> result;
735 result.reserve(m.size());
736 for (const auto& r : m){
737 result.push_back(std::inner_product(v.begin(), v.end(), r.begin(), 0.0));
738 }
739
740 return result;
741}
int r
Definition globals.cxx:22

◆ finalize()

StatusCode ISF::PunchThroughTool::finalize ( )
virtual

AlgTool finalize method.

Definition at line 274 of file PunchThroughTool.cxx.

275{
276 ATH_MSG_DEBUG( "[punchthrough] finalize() starting" );
277 for(auto & each : m_particles) {
278 delete each.second;
279 }
280
281 ATH_MSG_DEBUG( "[punchthrough] finalize() successful" );
282
283 return StatusCode::SUCCESS;
284}

◆ getAllParticles()

int ISF::PunchThroughTool::getAllParticles ( const ISF::ISFParticle & isfp,
ISFParticleVector & isfpCont,
CLHEP::HepRandomEngine * rndmEngine,
int pdg,
double interpEnergy,
double interpEta,
int numParticles = -1 ) const
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 451 of file PunchThroughTool.cxx.

452{
453
454 // get the current particle
455 PunchThroughParticle *p = m_particles.at(pdg);
456
457 // if no number of particles (=-1) was handed over as an argument
458 // -> get the number of particles from the pdf
459 if ( numParticles < 0 )
460 {
461 // prepare the function arguments for the PDFcreator class
462 std::vector<int> parameters;
463 parameters.push_back( std::round(interpEnergy) );
464 parameters.push_back( std::round(interpEta*100) );
465 // the maximum number of particles which should be produced
466 // if no maximum number is given, this is -1
467 int maxParticles = p->getMaxNumParticles();
468
469 // get the right number of punch-through particles
470 // and ensure that we do not create too many particles
471 do
472 {
473 numParticles = int( p->getNumParticlesPDF()->getRand(rndmEngine, parameters) );
474
475 // scale the number of particles if requested
476 numParticles = lround( numParticles *= p->getNumParticlesFactor() );
477 }
478 while ( (maxParticles >= 0.) && (numParticles > maxParticles) );
479 }
480
481 ATH_MSG_VERBOSE("[ punchthrough ] adding " << numParticles << " punch-through particles with pdg " << pdg);
482
483 // now create the exact number of particles which was just computed before
484 double energyRest = std::sqrt( isfp.momentum().mag2() + isfp.mass()*isfp.mass() );
485 double minEnergy = p->getMinEnergy();
486 int numCreated = 0;
487
488 for ( numCreated = 0; (numCreated < numParticles) && (energyRest > minEnergy); numCreated++ )
489 {
490 // create one particle which fullfills the right energy distribution
491 ISF::ISFParticle *par = getOneParticle(isfp, pdg, rndmEngine, interpEnergy, interpEta);
492
493 // if something went wrong
494 if (!par)
495 {
496 ATH_MSG_ERROR("[ punchthrough ] something went wrong while creating punch-through particles");
497 return 0;
498 }
499
500 // get the energy of the particle which was just created
501 const double restMass = m_particleDataTable->particle(std::abs(pdg))->mass();
502 double curEnergy = std::sqrt(par->momentum().mag2() + restMass*restMass);
503
504 // calculate the maximum energy to be available for all
505 // following punch-through particles created
506 energyRest -= curEnergy;
507
508 // add this ISFparticle to the vector
509 isfpCont.push_back( par );
510 }
511
512 // the number of particles which was created is numCreated
513 return (numCreated);
514}
#define ATH_MSG_ERROR(x)
ISF::ISFParticle * getOneParticle(const ISF::ISFParticle &isfp, int pdg, CLHEP::HepRandomEngine *rndmEngine, double interpEnergy, double interpEta) const
create exactly one punch-through particle with the given pdg and the given max energy
constexpr uint8_t maxParticles()

◆ getCorrelatedParticles()

int ISF::PunchThroughTool::getCorrelatedParticles ( const ISF::ISFParticle & isfp,
ISFParticleVector & isfpCont,
int doPdg,
int corrParticles,
CLHEP::HepRandomEngine * rndmEngine,
double interpEnergy,
double interpEta ) const
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 521 of file PunchThroughTool.cxx.

522{
523 // get the PunchThroughParticle class
524 PunchThroughParticle *p = m_particles.at(pdg);
525
526 const double initEnergy = std::sqrt( isfp.momentum().mag2() + isfp.mass()*isfp.mass() );
527
528 // (1.) decide if we do correlation or not
529 double rand = CLHEP::RandFlat::shoot(rndmEngine)
530 *(p->getFullCorrelationEnergy()-p->getMinCorrelationEnergy())
531 + p->getMinCorrelationEnergy();
532 if ( initEnergy < rand )
533 {
534 // here we do not do correlation
535 return getAllParticles(isfp, isfpCont, rndmEngine, pdg, interpEnergy, interpEta);
536 }
537
538 // (2.) if this point is reached, we do correlation
539 // decide which 2d correlation histogram to use
540 double *histDomains = p->getCorrelationHistDomains();
541 TH2F *hist2d = nullptr;
542 // compute the center values of the lowE and highE
543 // correlation histogram domains
544 if ( initEnergy < histDomains[1])
545 {
546 // initial energy lower than border between lowEnergy and highEnergy histogram domain
547 // --> choose lowEnergy correlation histogram
548 hist2d = p->getCorrelationLowEHist();
549 }
550 else
551 {
552 double rand = CLHEP::RandFlat::shoot(rndmEngine)*(histDomains[2]-histDomains[1])
553 + histDomains[1];
554 hist2d = ( initEnergy < rand) ? p->getCorrelationLowEHist()
555 : p->getCorrelationHighEHist();
556 }
557
558 // get the correlation 2d histogram
559
560 // now find out where on the x-axis the the bin for number of
561 // correlated particles is
562 Int_t xbin = hist2d->GetXaxis()->FindFixBin(corrParticles);
563 int numParticles = 0;
564 int maxParticles = p->getMaxNumParticles();
565 // now the distribution along the y-axis is a PDF for the number
566 // of 'pdg' particles
567 do
568 {
569 double rand = CLHEP::RandFlat::shoot(rndmEngine);
570 double sum = 0.;
571 for ( int ybin = 1; ybin <= hist2d->GetNbinsY(); ybin++ )
572 {
573 sum += hist2d->GetBinContent(xbin, ybin);
574 // check if we choose the current bin or not
575 if ( sum >= rand )
576 {
577 numParticles = ybin - 1;
578 break;
579 }
580 }
581 // scale the number of particles is requested
582 numParticles = lround( numParticles * p->getNumParticlesFactor() );
583 }
584 while ( (maxParticles >= 0.) && (numParticles > maxParticles) );
585
586 // finally create this exact number of particles
587 return getAllParticles(isfp, isfpCont, rndmEngine, pdg, interpEnergy, interpEta, numParticles);
588}
TH2F(name, title, nxbins, bins_par2, bins_par3, bins_par4, bins_par5=None, bins_par6=None, path='', **kwargs)

◆ getFloatAfterPatternInStr()

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

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

Definition at line 1299 of file PunchThroughTool.cxx.

1300{
1301 double num = 0.;
1302
1303 const std::string_view str( cstr);
1304 const std::string_view pattern( cpattern);
1305 const size_t pos = str.find(pattern);
1306
1307 if ( pos == std::string::npos)
1308 {
1309 ATH_MSG_WARNING("[ punchthrough ] unable to retrieve floating point number from string");
1310 return -999999999999.;
1311 }
1312 const std::string_view substring = str.substr(pos+pattern.length());
1313 std::from_chars(substring.data(), substring.data() + substring.size(), num);
1314 return num;
1315}
#define ATH_MSG_WARNING(x)

◆ getInfoMap()

auto ISF::PunchThroughTool::getInfoMap ( const std::string & mainNode,
const XMLCoreNode & doc )
private

Definition at line 800 of file PunchThroughTool.cxx.

801{
802 std::vector<InfoMap> xml_info;
803
804 for (const XMLCoreNode* node : doc.get_children (mainNode + "/info/item")) {
805 xml_info.emplace_back(*node);
806 }
807
808 return xml_info;
809}

◆ getOneParticle()

ISF::ISFParticle * ISF::PunchThroughTool::getOneParticle ( const ISF::ISFParticle & isfp,
int pdg,
CLHEP::HepRandomEngine * rndmEngine,
double interpEnergy,
double interpEta ) const
private

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

Definition at line 595 of file PunchThroughTool.cxx.

596{
597 // get a local copy of the needed punch-through particle class
598 PunchThroughParticle *p = m_particles.at(pdg);
599
600 // (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
601 int pcaCdfIterator = passedParamIterator(pdg, interpEta*100, m_xml_info_pca); //pca and cdf info should be of same size
602
603 ATH_MSG_DEBUG("[ punchthrough ] passedPCAIterator ==> pcaCdfIterator = "<< pcaCdfIterator <<" , pdg = "<< pdg <<" , interpEnergy = "<< interpEnergy <<" MeV, interpEta(*100) = "<< interpEta*100);
604
605 // (1.) decide if we create a particle or an anti-particle
606 int anti = 1;
607 if ( p->getdoAnti() )
608 {
609 // get a random-value
610 double rand = CLHEP::RandFlat::shoot(rndmEngine);
611 // 50/50 chance to be a particle or its anti-particle
612 if (rand > 0.5) anti = -1;
613 }
614
615 // (2.) get the right punch-through distributions
616 // prepare the function arguments for the PDFcreator class
617 std::vector<int> parInitEnergyEta;
618 parInitEnergyEta.push_back( std::round(interpEnergy) );
619 parInitEnergyEta.push_back( std::round(interpEta*100) );
620
621 //initialise variables to store punch through particle kinematics
622 double energy = 0.;
623 double deltaTheta = 0.;
624 double deltaPhi = 0.;
625 double momDeltaTheta = 0.;
626 double momDeltaPhi = 0.;
627
628 double principal_component_0 = 0.;
629 double principal_component_1 = 0.;
630 double principal_component_2 = 0.;
631 double principal_component_3 = 0.;
632 double principal_component_4 = 0.;
633 std::vector<double> transformed_variables;
634
635
636 principal_component_0 = p->getPCA0PDF()->getRand(rndmEngine, parInitEnergyEta);
637 principal_component_1 = p->getPCA1PDF()->getRand(rndmEngine, parInitEnergyEta);
638 principal_component_2 = p->getPCA2PDF()->getRand(rndmEngine, parInitEnergyEta);
639 principal_component_3 = p->getPCA3PDF()->getRand(rndmEngine, parInitEnergyEta);
640 principal_component_4 = p->getPCA4PDF()->getRand(rndmEngine, parInitEnergyEta);
641
642 ATH_MSG_DEBUG("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 );
643
644 std::vector<double> principal_components {
645 principal_component_0,
646 principal_component_1,
647 principal_component_2,
648 principal_component_3,
649 principal_component_4
650 };
651
652 transformed_variables = inversePCA(pcaCdfIterator,principal_components);
653
654 energy = inverseCdfTransform(transformed_variables.at(0), m_variable0_inverse_cdf[pcaCdfIterator]);
655 deltaTheta = inverseCdfTransform(transformed_variables.at(1), m_variable1_inverse_cdf[pcaCdfIterator]);
656 deltaPhi = inverseCdfTransform(transformed_variables.at(2), m_variable2_inverse_cdf[pcaCdfIterator]);
657 momDeltaTheta = inverseCdfTransform(transformed_variables.at(3), m_variable3_inverse_cdf[pcaCdfIterator]);
658 momDeltaPhi = inverseCdfTransform(transformed_variables.at(4), m_variable4_inverse_cdf[pcaCdfIterator]);
659
660 ATH_MSG_DEBUG("Transformed punch through kinematics: energy = "<< energy <<" MeV deltaTheta = "<< deltaTheta <<" deltaPhi = "<< deltaPhi <<" momDeltaTheta = "<< momDeltaTheta <<" momDeltaPhi = "<< momDeltaPhi );
661
662 energy *= p->getEnergyFactor(); // scale the energy if requested
663 if (energy < p->getMinEnergy()) {
664 energy = p->getMinEnergy() + 10;
665 }
666
667
668 // (2.2) get the particles delta theta relative to the incoming particle
669 double theta = 0;
670 // loop to keep theta within range [0,PI]
671 do
672 {
673 // decide if delta positive/negative
674 deltaTheta *= ( CLHEP::RandFlat::shoot(rndmEngine) > 0.5 ) ? 1. : -1.;
675 // calculate the exact theta value of the later created
676 // punch-through particle
677 theta = isfp.position().theta() + deltaTheta*p->getPosAngleFactor();
678
679 }
680 while ( (theta > M_PI) || (theta < 0.) );
681 // (2.3) get the particle's delta phi relative to the incoming particle
682
683 deltaPhi *= ( CLHEP::RandFlat::shoot(rndmEngine) > 0.5 ) ? 1. : -1.;
684
685 // keep phi within range [-PI,PI]
686 double phi = isfp.position().phi() + deltaPhi*p->getPosAngleFactor();
687 while ( std::fabs(phi) > 2*M_PI) phi /= 2.;
688 while (phi > M_PI) phi -= 2*M_PI;
689 while (phi < -M_PI) phi += 2*M_PI;
690
691 // (2.4) get the particle momentum delta theta, relative to its position
692 //
693 // loop to keep momTheta within range [0,PI]
694
695 double momTheta = 0.;
696 do
697 {
698 // decide if delta positive/negative
699 momDeltaTheta *= ( CLHEP::RandFlat::shoot(rndmEngine) > 0.5 ) ? 1. : -1.;
700 // calculate the exact momentum theta value of the later created
701 // punch-through particle
702 momTheta = theta + momDeltaTheta*p->getMomAngleFactor();
703
704 }
705 while ( (momTheta > M_PI) || (momTheta < 0.) );
706
707 // (2.5) get the particle momentum delta phi, relative to its position
708
709 momDeltaPhi *= ( CLHEP::RandFlat::shoot(rndmEngine) > 0.5 ) ? 1. : -1.;
710
711 double momPhi = phi + momDeltaPhi*p->getMomAngleFactor();
712 // keep momPhi within range [-PI,PI]
713 while ( std::fabs(momPhi) > 2*M_PI) momPhi /= 2.;
714 while (momPhi > M_PI) momPhi -= 2*M_PI;
715 while (momPhi < -M_PI) momPhi += 2*M_PI;
716
717 // (**) finally create the punch-through particle as a ISFParticle
718
719 ATH_MSG_DEBUG("createExitPs input parameters: doAnti? = "<< pdg*anti <<" energy = "<< energy <<" theta = "<< theta <<" phi = "<< phi <<" momTheta = "<< momTheta << " momPhi " << momPhi );
720
721
722 ISF::ISFParticle *par = createExitPs( isfp, pdg*anti, energy, theta, phi, momTheta, momPhi);
723
724 return par;
725}
#define M_PI
Scalar deltaPhi(const MatrixBase< Derived > &vec) const
std::vector< double > inversePCA(int pcaCdfIterator, std::vector< double > &variables) const
std::vector< std::map< double, double > > m_variable1_inverse_cdf
static double inverseCdfTransform(double variable, const std::map< double, double > &inverse_cdf_map)
std::vector< std::map< double, double > > m_variable3_inverse_cdf
std::vector< std::map< double, double > > m_variable2_inverse_cdf
std::vector< std::map< double, double > > m_variable0_inverse_cdf
(vector of map) for CDF mappings
int passedParamIterator(int pid, double eta, const std::vector< InfoMap > &mapvect) const
std::vector< InfoMap > m_xml_info_pca
infoMaps
std::vector< std::map< double, double > > m_variable4_inverse_cdf
ISF::ISFParticle * createExitPs(const ISF::ISFParticle &isfp, int PDGcode, double energy, double theta, double phi, double momTheta, double momPhi) const
create a ISF Particle state at the MS entrace containing a particle with the given properties

◆ getVariableCDFmappings()

std::map< double, double > ISF::PunchThroughTool::getVariableCDFmappings ( const XMLCoreNode * node)
staticprivate

Definition at line 888 of file PunchThroughTool.cxx.

888 {
889
890 std::map<double, double> mappings;
891
892 if (node) {
893 for (const XMLCoreNode* child : node->get_children ("CDFmap")) {
894 double ref = child->get_double_attrib ("ref");
895 double quant = child->get_double_attrib ("quant");
896 mappings[ref] = quant;
897 }
898 }
899
900 return mappings;
901}
const boost::regex ref(r_ef)
static int quant(double min, double max, unsigned nSteps, double val)
std::vector< const XMLCoreNode * > get_children(const std::string &path="*") const
Return all children matching a pattern.

◆ initialize()

StatusCode ISF::PunchThroughTool::initialize ( )
virtual

AlgTool initialize method.

Definition at line 77 of file PunchThroughTool.cxx.

78{
79 ATH_MSG_DEBUG( "initialize()" );
80
81 // initialise punch through classifier
82 if (m_punchThroughClassifier.retrieve().isFailure() )
83 {
84 ATH_MSG_ERROR (m_punchThroughClassifier.propertyName() << ": Failed to retrieve tool " << m_punchThroughClassifier.type());
85 return StatusCode::FAILURE;
86 }
87
88 // resolving lookuptable file
89 std::string resolvedFileName = PathResolverFindCalibFile (m_filenameLookupTable);
90 if (resolvedFileName.empty()) {
91 ATH_MSG_ERROR( "[ punchthrough ] Parametrisation file '" << m_filenameLookupTable << "' not found" );
92 return StatusCode::FAILURE;
93 }
94 ATH_MSG_INFO( "[ punchthrough ] Parametrisation file found: " << resolvedFileName );
95
96 // open the LookupTable file
97 m_fileLookupTable = new TFile( resolvedFileName.c_str(), "READ");
98 if (!m_fileLookupTable) {
99 ATH_MSG_WARNING("[ punchthrough ] unable to open the lookup-table for the punch-through simulation (file does not exist)");
100 return StatusCode::FAILURE;
101 }
102
103 if (!m_fileLookupTable->IsOpen()) {
104 ATH_MSG_WARNING("[ punchthrough ] unable to open the lookup-table for the punch-through simulation (wrong or empty file?)");
105 return StatusCode::FAILURE;
106 }
107
108 //retrieve inverse CDF config file
110 {
111 ATH_MSG_WARNING("[ punchthrough ] unable to open or read the inverse CDF config");
112 return StatusCode::FAILURE;
113 }
114
115 //retrieve inverse PCA config file
117 {
118 ATH_MSG_WARNING("[ punchthrough ] unable to open or read the inverse PCA config");
119 return StatusCode::FAILURE;
120 }
121
122 //check first the size of infoMap for both PCA and CDF, they should be equal
123 if (!(m_xml_info_pca.size() == m_xml_info_cdf.size()))
124 {
125 ATH_MSG_WARNING("[ punchthrough ] size of infoMap for PCA and CDF differs! Something is wrong with input xml files.");
126 return StatusCode::FAILURE;
127 }
128
129 // retrieve the ParticleProperties handle
130 ATH_CHECK( m_particlePropSvc.retrieve() );
131
132 // and the particle data table
135 {
136 ATH_MSG_FATAL( " [ punchthrough ] Could not get ParticleDataTable! Cannot associate pdg code with charge! Abort. " );
137 return StatusCode::FAILURE;
138 }
139
140 // Geometry identifier service
141 if ( !m_geoIDSvc.empty() && m_geoIDSvc.retrieve().isFailure())
142 {
143 ATH_MSG_FATAL ( "[ punchthrough ] Could not retrieve GeometryIdentifier Service. Abort");
144 return StatusCode::FAILURE;
145 }
146
147 //envelope definition service
148 if (m_envDefSvc.retrieve().isFailure() )
149 {
150 ATH_MSG_ERROR( "[ punchthrough ] Could not retrieve " << m_envDefSvc );
151 return StatusCode::FAILURE;
152 }
153
154 //--------------------------------------------------------------------------------
155 // register all the punch-through particles which will be simulated
156 for ( unsigned int num = 0; num < m_punchThroughParticles.size(); num++ )
157 {
158 const int pdg = m_punchThroughParticles[num];
159 // if no information is given on the creation of anti-particles -> do not simulate anti-particles
160 const bool doAnti = ( num < m_doAntiParticles.size() ) ? m_doAntiParticles[num] : false;
161 // if no information is given on the minimum energy -> take 50. MeV as default
162 const double minEnergy = ( num < m_minEnergy.size() ) ? m_minEnergy[num] : 50.;
163 // if no information is given on the maximum number of punch-through particles -> take -1 as default
164 const int maxNum = ( num < m_minEnergy.size() ) ? m_maxNumParticles[num] : -1;
165 // if no information is given on the scale factor for the number of particles -> take 1. as defaulft
166 const double numFactor = ( num < m_numParticlesFactor.size() ) ? m_numParticlesFactor[num] : 1.;
167 // if no information is given on the position angle factor -> take 1.
168 const double posAngleFactor = ( num < m_posAngleFactor.size() ) ? m_posAngleFactor[num] : 1.;
169 // if no information is given on the momentum angle factor -> take 1.
170 const double momAngleFactor = ( num < m_momAngleFactor.size() ) ? m_momAngleFactor[num] : 1.;
171 // if no information is given on the scale factor for the energy -> take 1. as default
172 const double energyFactor = ( num < m_energyFactor.size() ) ? m_energyFactor[num] : 1.;
173
174 // register the particle
175 ATH_MSG_VERBOSE("[ punchthrough ] registering punch-through particle type with pdg = " << pdg );
176 if ( registerParticle( pdg, doAnti, minEnergy, maxNum, numFactor, energyFactor, posAngleFactor, momAngleFactor )
177 != StatusCode::SUCCESS)
178 {
179 ATH_MSG_ERROR("[ punchthrough ] unable to register punch-through particle type with pdg = " << pdg);
180 }
181 }
182
183 // TODO: implement punch-through parameters for different m_pdgInitiators
184 // currently m_pdgInitiators is only used to filter out particles
185
186 // check if more correlations were given than particle types were registered
187 unsigned int numCorrelations = m_correlatedParticle.size();
188 if ( numCorrelations > m_punchThroughParticles.size() )
189 {
190 ATH_MSG_WARNING("[ punchthrough ] more punch-through particle correlations are given, than punch-through particle types are registered (skipping the last ones)");
191 numCorrelations = m_punchThroughParticles.size();
192 }
193
194 // now register correlation between particles
195 for ( unsigned int num = 0; num < numCorrelations; num++ )
196 {
197 const int pdg1 = m_punchThroughParticles[num];
198 const int pdg2 = m_correlatedParticle[num];
199 const double fullCorrEnergy = ( num < m_fullCorrEnergy.size() ) ? m_fullCorrEnergy[num] : 0.;
200 const double minCorrEnergy = ( num < m_minCorrEnergy.size() ) ? m_minCorrEnergy[num] : 0.;
201
202 // if correlatedParticle==0 is given -> no correlation
203 if ( ! pdg2) continue;
204 // register it
205 if ( registerCorrelation(pdg1, pdg2, minCorrEnergy, fullCorrEnergy) != StatusCode::SUCCESS )
206 {
207 ATH_MSG_ERROR("[ punchthrough ] unable to register punch-through particle correlation for pdg1=" << pdg1 << " pdg2=" << pdg2 );
208 }
209 }
210
211 // get the calo-MS border coordinates. Look at calo and MS geometry definitions, if same R and Z -> boundary surface
212
213 const RZPairVector* rzMS = &(m_envDefSvc->getMuonRZBoundary());
214 const RZPairVector* rzCalo = &(m_envDefSvc->getCaloRZBoundary());
215
216 bool found1, found2;
217 found1=false; found2=false;
218
219 for ( size_t i=0; i<rzCalo->size();++i)
220 {
221 const double r_tempCalo = rzCalo->at(i).first;
222 const double z_tempCalo = rzCalo->at(i).second;
223
224 if (r_tempCalo> m_beamPipe)
225 {
226 for ( size_t j=0; j<rzMS->size();++j)
227 {
228 const double r_tempMS =rzMS->at(j).first;
229 const double z_tempMS =rzMS->at(j).second;
230
231 if (r_tempCalo==r_tempMS && z_tempCalo==z_tempMS && found1==false )
232 {
233 m_R1=r_tempMS;
234 m_z1=std::fabs(z_tempMS);
235 found1=true;
236 continue;
237 }
238 else if (r_tempCalo==r_tempMS && z_tempCalo==z_tempMS && r_tempCalo!=m_R1 && std::fabs(z_tempCalo)!=m_z1)
239 {
240 m_R2=r_tempMS;
241 m_z2=std::fabs(z_tempMS);
242 found2=true;
243 }
244 }
245
246 if (found1==true && found2==true) break;
247 }
248 }
249
250 //in case geometry description changes
251 if (found1 == false) ATH_MSG_ERROR ("first coordinate of calo-MS border not found");
252 if (found2 == false) ATH_MSG_ERROR ("second coordinate of calo-MS border not found; first one is: R1 ="<<m_R1<<" z1 ="<<m_z1);
253
254 //now order the found values
255 double r_temp, z_temp;
256 if (m_R1>m_R2) { r_temp=m_R1; m_R1=m_R2; m_R2=r_temp; } //m_R1 - smaller one
257 if (m_z1<m_z2) { z_temp=m_z1; m_z1=m_z2; m_z2=z_temp; } //m_z1 - bigger one
258
259 if (m_R1==m_R2 || m_z1==m_z2) ATH_MSG_ERROR ("[punch-though] Bug in propagation calculation! R1="<<m_R1<<" R2 = "<<m_R2<<" z1="<<m_z1<<" z2= "<<m_z2 );
260 else ATH_MSG_DEBUG ("calo-MS boundary coordinates: R1="<<m_R1<<" R2 = "<<m_R2<<" z1="<<m_z1<<" z2= "<<m_z2);
261
262 // close the file with the lookuptable
263 m_fileLookupTable->Close();
264
265 ATH_MSG_INFO( "punchthrough initialization is successful" );
266 return StatusCode::SUCCESS;
267}
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_FATAL(x)
#define ATH_MSG_INFO(x)
std::string PathResolverFindCalibFile(const std::string &logical_file_name)
std::vector< RZPair > RZPairVector
Definition RZPair.h:18
StringProperty m_filenameInverseCDF
IntegerArrayProperty m_correlatedParticle
DoubleArrayProperty m_minCorrEnergy
ServiceHandle< IPartPropSvc > m_particlePropSvc
StringProperty m_filenameLookupTable
Properties.
StringProperty m_filenameInversePCA
BooleanArrayProperty m_doAntiParticles
DoubleArrayProperty m_energyFactor
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...
DoubleProperty m_beamPipe
beam pipe radius
DoubleArrayProperty m_momAngleFactor
DoubleArrayProperty m_fullCorrEnergy
std::vector< InfoMap > m_xml_info_cdf
StatusCode registerParticle(int pdgID, 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
DoubleArrayProperty m_numParticlesFactor
StatusCode initializeInverseCDF(const std::string &quantileTransformerConfigFile)
DoubleArrayProperty m_minEnergy
IntegerArrayProperty m_maxNumParticles
double m_R1
calo-MS borders
StatusCode initializeInversePCA(const std::string &inversePCAConfigFile)
TFile * m_fileLookupTable
ROOT objects.
ServiceHandle< IEnvelopeDefSvc > m_envDefSvc
DoubleArrayProperty m_posAngleFactor
IntegerArrayProperty m_punchThroughParticles
float j(const xAOD::IParticle &, const xAOD::TrackMeasurementValidation &hit, const Eigen::Matrix3d &jab_inv)

◆ initializeInverseCDF()

StatusCode ISF::PunchThroughTool::initializeInverseCDF ( const std::string & quantileTransformerConfigFile)
private

Definition at line 862 of file PunchThroughTool.cxx.

862 {
863 //parse xml that contains config for inverse CDF for each of punch through particle kinematics
864 XMLCoreParser p;
865 std::unique_ptr<XMLCoreNode> doc = p.parse (inverseCdfConfigFile);
866
867 ATH_MSG_INFO( "[ punchthrough ] Loading inverse CDF: " << inverseCdfConfigFile);
868
869 //check info first
870 m_xml_info_cdf = getInfoMap("CDFMappings", *doc);
871
872 for (unsigned int i = 0; i < m_xml_info_cdf.size(); i++) {
873 ATH_MSG_DEBUG( "[PunchThroughG4Tool] m_xml_info_cdf[" << i << "].name = " << m_xml_info_cdf[i].name);
874
875 if (const XMLCoreNode* node = doc->get_child ("CDFMappings/" + m_xml_info_cdf[i].name)) {
876 m_variable0_inverse_cdf.push_back(getVariableCDFmappings (node->get_child("variable0")));
877 m_variable1_inverse_cdf.push_back(getVariableCDFmappings (node->get_child("variable1")));
878 m_variable2_inverse_cdf.push_back(getVariableCDFmappings (node->get_child("variable2")));
879 m_variable3_inverse_cdf.push_back(getVariableCDFmappings (node->get_child("variable3")));
880 m_variable4_inverse_cdf.push_back(getVariableCDFmappings (node->get_child("variable4")));
881 }
882 }
883
884 return StatusCode::SUCCESS;
885}
static std::map< double, double > getVariableCDFmappings(const XMLCoreNode *node)
std::vector< InfoMap > getInfoMap(const std::string &mainNode, const XMLCoreNode &doc)

◆ initializeInversePCA()

StatusCode ISF::PunchThroughTool::initializeInversePCA ( const std::string & inversePCAConfigFile)
private

Definition at line 820 of file PunchThroughTool.cxx.

820 {
821
822 XMLCoreParser p;
823 std::unique_ptr<XMLCoreNode> doc = p.parse (inversePCAConfigFile);
824
825 ATH_MSG_INFO( "[ punchthrough ] Loading inversePCA: " << inversePCAConfigFile);
826
827 //check info first
828 m_xml_info_pca = getInfoMap("PCAinverse", *doc);
829
830 auto getRow = [] (const XMLCoreNode* node,
831 const std::string& attrib,
832 int imax)
833 -> std::vector<double>
834 {
835 std::vector<double> row;
836 for (int i = 0; i <= imax; ++i) {
837 // Dynamically create property name
838 std::string propName = attrib + std::to_string(i);
839 row.push_back (node->get_double_attrib (propName));
840 }
841 return row;
842 };
843
844 for (unsigned int i = 0; i < m_xml_info_pca.size(); i++) {
845 std::vector<std::vector<double>> PCA_matrix;
846
847 if (const XMLCoreNode* node = doc->get_child ("PCAinverse/" + m_xml_info_pca[i].name)) {
848 for (const XMLCoreNode* matrix : node->get_children ("PCAmatrix")) {
849 PCA_matrix.push_back(getRow(matrix, "comp_", 4));
850 }
851 for (const XMLCoreNode* means : node->get_children ("PCAmeans")) {
852 m_PCA_means.push_back(getRow(means, "mean_", 4));
853 }
854 }
855
856 m_inverse_PCA_matrix.push_back(std::move(PCA_matrix));
857 }
858
859 return StatusCode::SUCCESS;
860}
int imax(int i, int j)
std::vector< std::vector< std::vector< double > > > m_inverse_PCA_matrix
pca vectors
std::vector< std::vector< double > > m_PCA_means
double get_double_attrib(const std::string &name) const
Retrieve the value of an attribute as a double.
row
Appending html table to final .html summary file.

◆ interpolateEnergy()

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

Definition at line 918 of file PunchThroughTool.cxx.

918 {
919
920 ATH_MSG_DEBUG("[ punchthrough ] interpolating incoming energy: " << energy);
921
922 std::string energyPointsString;
923 for (auto element:m_energyPoints){
924 energyPointsString += std::to_string(element) + " ";
925 }
926
927 ATH_MSG_DEBUG("[ punchthrough ] available energy points: " << energyPointsString);
928
929 auto const upperEnergy = std::upper_bound(m_energyPoints.begin(), m_energyPoints.end(), energy);
930
931 if(upperEnergy == m_etaPoints.end()){ //if no energy greater than input energy, choose greatest energy
932 ATH_MSG_DEBUG("[ punchthrough ] incoming energy > largest energy point, returning greatest energy point: " << m_energyPoints.back());
933 return m_energyPoints.back();
934 }
935 else if(upperEnergy == m_etaPoints.begin()){ //if smallest energy greater than input energy, choose smallest energy
936 ATH_MSG_DEBUG("[ punchthrough ] incoming energy < smallest energy point, returning smallest energy point: " << *upperEnergy);
937 return *upperEnergy;
938 }
939
940 ATH_MSG_DEBUG("[ punchthrough ] energy points upper_bound: "<< *upperEnergy);
941
942 double randomShoot = CLHEP::RandFlat::shoot(rndmEngine);
943
944 ATH_MSG_DEBUG("[ punchthrough ] Shooting random number: "<< randomShoot);
945
946 double midPoint = *std::prev(upperEnergy)*M_SQRT2;
947
948 if(energy < midPoint){ //if energy smaller than mid point in log(energy)
949
950 double distance = std::abs(energy - *std::prev(upperEnergy))/((midPoint) - *std::prev(upperEnergy));
951
952 ATH_MSG_DEBUG( "[ punchthrough ] incoming energy is closest to prev(upper_bound) in log(energy), distance: " << distance );
953
954 if(randomShoot < distance){
955 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot < distance, returning upper_bound " << *upperEnergy );
956 return *upperEnergy;
957 }
958 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot > distance, returning prev(upper_bound) " << *std::prev(upperEnergy) );
959
960 return *std::prev(upperEnergy);
961 }
962 else if(energy > midPoint){ //if energy greater than mid point in log(energy)
963
964 double distance = std::abs(energy - *upperEnergy)/((*upperEnergy - midPoint));
965
966 ATH_MSG_DEBUG( "[ punchthrough ] incoming energy is closest to upper_bound in log(energy), distance: " << distance );
967
968 if(randomShoot < distance){
969 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot < distance, returning prev(upper_bound) " << *std::prev(upperEnergy) );
970 return *std::prev(upperEnergy);
971 }
972 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot > distance, returning upper_bound " << *upperEnergy );
973 return *upperEnergy;
974 }
975
976 return *upperEnergy;
977}
std::vector< double > m_energyPoints
energy and eta points in param
std::vector< double > m_etaPoints
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space

◆ interpolateEta()

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

Definition at line 979 of file PunchThroughTool.cxx.

979 {
980
981 double absEta = std::abs(eta);
982
983 ATH_MSG_DEBUG("[ punchthrough ] interpolating incoming abs(eta): " << absEta);
984
985 std::string etaPointsString;
986 for (auto element:m_etaPoints){
987 etaPointsString += std::to_string(element) + " ";
988 }
989
990 ATH_MSG_DEBUG("[ punchthrough ] available eta points: " << etaPointsString);
991
992 auto const upperEta = std::upper_bound(m_etaPoints.begin(), m_etaPoints.end(), absEta);
993
994 if(upperEta == m_etaPoints.end()){
995 ATH_MSG_DEBUG("[ punchthrough ] incoming abs(eta) > largest eta point, returning greatest eta point: " << m_etaPoints.back());
996 return m_etaPoints.back();
997 }
998
999
1000 ATH_MSG_DEBUG("[ punchthrough ] eta points upper_bound: "<< *upperEta);
1001
1002 double randomShoot = CLHEP::RandFlat::shoot(rndmEngine);
1003
1004 ATH_MSG_DEBUG("[ punchthrough ] Shooting random number: "<< randomShoot);
1005
1006 if(std::abs(absEta - *upperEta) < std::abs(absEta - *std::prev(upperEta))){
1007
1008 double distance = std::abs(absEta - *upperEta)/((*upperEta - *std::prev(upperEta))/2);
1009
1010 ATH_MSG_DEBUG( "[ punchthrough ] abs(eta) is closer to eta points upper_bound, distance: " << distance );
1011
1012 if(randomShoot > distance){
1013 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot > distance, returning upper_bound " << *upperEta );
1014 return *upperEta;
1015 }
1016
1017 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot < distance, returning prev(upper_bound) " << *std::prev(upperEta) );
1018
1019 return *std::prev(upperEta);
1020 }
1021 else if(std::abs(absEta - *std::prev(upperEta)) < std::abs(absEta - *upperEta)){
1022
1023 if(std::prev(upperEta) == m_etaPoints.begin()){
1024 ATH_MSG_DEBUG( "[ punchthrough ] prev of upper bound is begin, returning that: " << *std::prev(upperEta) );
1025 return *std::prev(upperEta);
1026 }
1027
1028 double distance = std::abs(absEta - *std::prev(upperEta))/((*std::prev(upperEta) - *std::prev(std::prev(upperEta)))/2);
1029
1030 ATH_MSG_DEBUG( "[ punchthrough ] abs(eta) is closer to eta points prev(upper_bound), distance: " << distance );
1031
1032 if(randomShoot > distance){
1033 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot > distance, returning prev(prev(upper_bound)) " << *std::prev(std::prev(upperEta)) );
1034
1035 return *std::prev(std::prev(upperEta));
1036 }
1037 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot < distance, returning prev(upper_bound) " << *std::prev(upperEta) );
1038
1039 return *std::prev(upperEta);
1040 }
1041
1042 return *std::prev(upperEta);
1043}
Scalar eta() const
pseudorapidity method
bool absEta(const xAOD::TauJet &tau, float &out)

◆ inverseCdfTransform()

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

Definition at line 903 of file PunchThroughTool.cxx.

903 {
904
905 double norm_cdf = normal_cdf(variable);
906
907 auto upper = inverse_cdf_map.upper_bound(norm_cdf);
908 auto lower = upper--;
909
910 double m = (upper->second - lower->second)/(upper->first - lower->first);
911 double c = lower->second - m * lower->first;
912 double transformed = m * norm_cdf + c;
913
914 return transformed;
915
916}
int upper(int c)
static double normal_cdf(double x)

◆ inversePCA()

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

Definition at line 811 of file PunchThroughTool.cxx.

812{
813 std::vector<double> transformed_variables = dotProduct(m_inverse_PCA_matrix[pcaCdfIterator], variables);
814
815 std::transform (transformed_variables.begin(), transformed_variables.end(), m_PCA_means[pcaCdfIterator].begin(), transformed_variables.begin(), std::plus<double>()); // + means
816
817 return transformed_variables;
818}
static std::vector< double > dotProduct(const std::vector< std::vector< double > > &m, const std::vector< double > &v)

◆ normal_cdf()

double ISF::PunchThroughTool::normal_cdf ( double x)
staticprivate

Definition at line 727 of file PunchThroughTool.cxx.

727 {
728
729 return 0.5 * TMath::Erfc(-x * M_SQRT1_2);
730}
#define x

◆ passedParamIterator()

int ISF::PunchThroughTool::passedParamIterator ( int pid,
double eta,
const std::vector< InfoMap > & mapvect ) const
private

Definition at line 767 of file PunchThroughTool.cxx.

768{
769 //convert the pid to absolute value and string for query
770 int pidStrSingle = std::abs(pid);
771 //STEP 1
772 //filter items matching pid first
773
774 for (int i = 0; const InfoMap& m : mapvect) {
775 ++i;
776 const std::vector<int>& v = m.pidStr;
777 if(std::find(v.begin(), v.end(),pidStrSingle)==v.end()) continue;
778 assert(m.etaMaxs.size() == m.etaMins.size());
779 for (unsigned int j = 0; j < m.etaMins.size(); j++){ // assume size etaMinsVect == etaMaxsVect
780 double etaMinToCompare = m.etaMins[j];
781 double etaMaxToCompare = m.etaMaxs[j];
782 if((eta >= etaMinToCompare) && (eta < etaMaxToCompare)){
783 //PASS CONDITION
784 //then choose the passing one and note it's iterator
785 return i-1; //in case more than 1 match (ambiguous case)
786 }
787 }
788 }
789 return 0;
790}

◆ propagator()

Amg::Vector3D ISF::PunchThroughTool::propagator ( double theta,
double phi ) const
private

get particle through the calorimeter

Definition at line 1317 of file PunchThroughTool.cxx.

1318{
1319 // phi, theta - direction of the punch-through particle coming into calo
1320 //particle propagates inside the calorimeter along the straight line
1321 //coordinates of this particles when exiting the calo (on calo-MS boundary)
1322
1323 double x, y, z, r;
1324
1325 // cylinders border angles
1326 const double theta1 = atan (m_R1/m_z1);
1327 const double theta2 = atan (m_R1/m_z2);
1328 const double theta3 = atan (m_R2/m_z2);
1329 //where is the particle
1330
1331 if (theta >= 0 && theta < theta1)
1332 {
1333 z = m_z1;
1334 r = std::fabs (m_z1*tan(theta));
1335 }
1336 else if (theta >= theta1 && theta < theta2)
1337 {
1338 z = m_R1/tan(theta);
1339 r = m_R1;
1340 }
1341 else if (theta >= theta2 && theta < theta3)
1342 {
1343 z = m_z2;
1344 r = std::fabs(m_z2*tan(theta));;
1345 }
1346 else if (theta >= theta3 && theta < (TMath::Pi()-theta3) )
1347 {
1348 z = m_R2/tan(theta);
1349 r = m_R2;
1350 }
1351 else if (theta >= (TMath::Pi()-theta3) && theta < (TMath::Pi()-theta2) )
1352 {
1353 z = -m_z2;
1354 r = std::fabs(m_z2*tan(theta));
1355 }
1356 else if (theta >= (TMath::Pi()-theta2) && theta < (TMath::Pi()-theta1) )
1357 {
1358 z = m_R1/tan(theta);
1359 r = m_R1;
1360 }
1361 else if (theta >= (TMath::Pi()-theta1) && theta <= TMath::Pi() )
1362 {
1363 z = -m_z1;
1364 r = std::fabs(m_z1*tan(theta));
1365 }
1366
1367 //parallel universe
1368 else
1369 {
1370 ATH_MSG_WARNING ( "Given theta angle is incorrect, setting particle position to (0, 0, 0)");
1371 x = 0.0; y = 0.0; z = 0.0; r = 0.0;
1372 }
1373
1374 x = r*cos(phi);
1375 y = r*sin(phi);
1376 Amg::Vector3D pos(x, y, z);
1377
1378 ATH_MSG_DEBUG("position of produced punch-through particle: x = "<< x <<" y = "<< y <<" z = "<< z<<" r = "<< pos.perp() <<"std::sqrt(x^2+y^2) = "<< std::sqrt(x*x+y*y) );
1379 ATH_MSG_DEBUG("GeoID thinks: Calo: "<< m_geoIDSvc->inside(pos, AtlasDetDescr::fAtlasCalo) <<" MS: "<< m_geoIDSvc->inside(pos,AtlasDetDescr::fAtlasMS));
1380
1381 return pos;
1382}
#define y
#define z

◆ readLookuptablePDF()

std::unique_ptr< ISF::PDFcreator > ISF::PunchThroughTool::readLookuptablePDF ( int pdgID,
const std::string & folderName )
private

reads out the lookuptable for the given type of particle

Definition at line 1183 of file PunchThroughTool.cxx.

1184{
1185
1186 // will hold the PDFcreator class which will be returned at the end
1187 // this will store the distributions for the punch through particles
1188 // (as map of energy & eta of the incoming particle)
1189 //PDFcreator *pdf = new PDFcreator();
1190 std::unique_ptr<ISF::PDFcreator> pdf = std::make_unique<ISF::PDFcreator>();
1191
1192 //Get directory object
1193 std::stringstream dirName;
1194 dirName << folderName << pdg;
1195 pdf->setName(dirName.str().c_str());
1196
1197 TDirectory * dir = (TDirectory*)m_fileLookupTable->Get(dirName.str().c_str());
1198 if(! dir)
1199 {
1200 ATH_MSG_ERROR( "[ punchthrough ] unable to retrieve directory object ("<< folderName << pdg << ")" );
1201 return nullptr;
1202 }
1203
1204
1205
1206 //Get list of all objects in directory
1207 TIter keyList(dir->GetListOfKeys());
1208 TKey *key;
1209
1210 while ((key = (TKey*)keyList())) {
1211
1212 //Get histogram object from key and its name
1213 TH1* hist = nullptr;
1214
1215 std::string histName;
1216 if(strcmp(key->GetClassName(), "TH1F") == 0){
1217 hist = (TH1*)key->ReadObj();
1218 histName = hist->GetName();
1219 }
1220
1221 //extract energy and eta from hist name 6 and 1 to position delimeters correctly
1222 std::string strEnergy = histName.substr( histName.find_first_of('E') + 1, histName.find_first_of('_')-histName.find_first_of('E') - 1 );
1223 histName.erase(0, histName.find_first_of('_') + 1);
1224 std::string strEtaMin = histName.substr( histName.find("etaMin") + 6, histName.find_first_of('_') - histName.find("etaMin") - 6 );
1225 histName.erase(0, histName.find('_') + 1);
1226 std::string strEtaMax = histName.substr( histName.find("etaMax") + 6, histName.length());
1227
1228 //create integers to store in map
1229 const int energy = std::stoi(strEnergy);
1230 const int etaMin = std::stoi(strEtaMin);
1231
1232 //Add entry to pdf map
1233 pdf->addToEnergyEtaHist1DMap(energy, etaMin, hist);
1234
1235 //create doubles to store energy and eta points for interpolation
1236 const double energyDbl = static_cast<double>(energy);
1237 const double etaDbl = static_cast<double>(etaMin)/100.;
1238
1239 //create vectors to store the eta and energy points, this allows us to interpolate
1240 if (std::find(m_energyPoints.begin(), m_energyPoints.end(), energyDbl) == m_energyPoints.end()) {
1241 m_energyPoints.push_back(energyDbl);
1242 }
1243 if (std::find(m_etaPoints.begin(), m_etaPoints.end(), etaDbl) == m_etaPoints.end()) {
1244 m_etaPoints.push_back(etaDbl);
1245 }
1246
1247 }
1248
1249
1250
1251 return pdf;
1252}

◆ registerCorrelation()

StatusCode ISF::PunchThroughTool::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 1128 of file PunchThroughTool.cxx.

1130{
1131 // find the given pdgs in the registered particle ids
1132 std::map<int, PunchThroughParticle*>::iterator location1 = m_particles.find(pdgID1);
1133 std::map<int, PunchThroughParticle*>::iterator location2 = m_particles.find(pdgID2);
1134
1135 // if at least one of the given pdgs was not registered yet -> return an error
1136 if ( (location1 == m_particles.end()) || (location2 == m_particles.end()) )
1137 return StatusCode::FAILURE;
1138
1139 // now look for the correlation histograms
1140 std::stringstream name;
1141 name << "NumExitCorrelations/x_PDG" << std::abs(pdgID1) << "__y_PDG" << std::abs(pdgID2) << "__lowE";
1142 TH2F *hist1_2_lowE = (TH2F*)m_fileLookupTable->Get(name.str().c_str());
1143 name.str("");
1144 name << "NumExitCorrelations/x_PDG" << std::abs(pdgID1) << "__y_PDG" << std::abs(pdgID2) << "__highE";
1145 TH2F *hist1_2_highE = (TH2F*)m_fileLookupTable->Get(name.str().c_str());
1146 name.str("");
1147 name << "NumExitCorrelations/x_PDG" << std::abs(pdgID2) << "__y_PDG" << std::abs(pdgID1) << "__lowE";
1148 TH2F *hist2_1_lowE = (TH2F*)m_fileLookupTable->Get(name.str().c_str());
1149 name.str("");
1150 name << "NumExitCorrelations/x_PDG" << std::abs(pdgID2) << "__y_PDG" << std::abs(pdgID1) << "__highE";
1151 TH2F *hist2_1_highE = (TH2F*)m_fileLookupTable->Get(name.str().c_str());
1152 // check if the histograms exist
1153 if ( (!hist1_2_lowE) || (!hist2_1_lowE) || (!hist1_2_highE) || (!hist2_1_highE) )
1154 {
1155 ATH_MSG_ERROR("[ punchthrough ] unable to retrieve the correlation data for PDG IDs " << pdgID1 << " and " << pdgID2);
1156 return StatusCode::FAILURE;
1157 }
1158
1159 // TODO: if only one of the two histograms exists, create the other one
1160 // by mirroring the data
1161
1162 const double lowE = getFloatAfterPatternInStr( hist1_2_lowE->GetTitle(), "elow_");
1163 const double midE = getFloatAfterPatternInStr( hist1_2_lowE->GetTitle(), "ehigh_");
1164 //TODO: check if the same:
1165 // double midE = getFloatAfterPatternInStr( hist1_2_lowE->GetTitle(), "elow_");
1166 const double upperE = getFloatAfterPatternInStr( hist1_2_highE->GetTitle(), "ehigh_");
1167 // now store the correlation either way id1->id2 and id2->id1
1168 m_particles[pdgID1]->setCorrelation(pdgID2, hist2_1_lowE, hist2_1_highE,
1169 minCorrEnergy, fullCorrEnergy,
1170 lowE, midE, upperE);
1171
1172 m_particles[pdgID2]->setCorrelation(pdgID1, hist1_2_lowE, hist1_2_highE,
1173 minCorrEnergy, fullCorrEnergy,
1174 lowE, midE, upperE);
1175 return StatusCode::SUCCESS;
1176}
double getFloatAfterPatternInStr(const char *str, const char *pattern)
get the floating point number in a string, after the given pattern

◆ registerParticle()

StatusCode ISF::PunchThroughTool::registerParticle ( int pdgID,
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 1051 of file PunchThroughTool.cxx.

1054{
1055 // read in the data needed to construct the distributions for the number of punch-through particles
1056
1057 // (1.) get the distribution function for the number of punch-through particles
1058 std::unique_ptr<ISF::PDFcreator> pdf_num(readLookuptablePDF(pdg, "FREQ_PDG"));
1059 if (!pdf_num ) return StatusCode::FAILURE; // return error if something went wrong
1060
1061 // (2.) get the PDF for the punch-through energy
1062 std::unique_ptr<PDFcreator> pdf_pca0 (readLookuptablePDF(pdg, "PCA0_PDG"));
1063 if (!pdf_pca0)
1064 {
1065 return StatusCode::FAILURE; // return error if something went wrong
1066 }
1067
1068 // (3.) get the PDF for the punch-through particles difference in
1069 // theta compared to the incoming particle
1070 std::unique_ptr<PDFcreator> pdf_pca1 (readLookuptablePDF(pdg, "PCA1_PDG"));
1071 if (!pdf_pca1)
1072 {
1073 return StatusCode::FAILURE;
1074 }
1075
1076 // (4.) get the PDF for the punch-through particles difference in
1077 // phi compared to the incoming particle
1078 std::unique_ptr<PDFcreator> pdf_pca2 (readLookuptablePDF(pdg, "PCA2_PDG"));
1079 if (!pdf_pca2)
1080 {
1081 return StatusCode::FAILURE;
1082 }
1083
1084 // (5.) get the PDF for the punch-through particle momentum delta theta angle
1085 std::unique_ptr<PDFcreator> pdf_pca3 (readLookuptablePDF(pdg, "PCA3_PDG"));
1086 if (!pdf_pca3)
1087 {
1088 return StatusCode::FAILURE;
1089 }
1090
1091 // (6.) get the PDF for the punch-through particle momentum delta phi angle
1092 std::unique_ptr<PDFcreator> pdf_pca4 (readLookuptablePDF(pdg, "PCA4_PDG"));
1093 if (!pdf_pca4)
1094 {
1095 return StatusCode::FAILURE;
1096 }
1097
1098 // (7.) now finally store all this in the right std::map
1099 PunchThroughParticle *particle = new PunchThroughParticle(pdg, doAntiparticle);
1100 particle->setNumParticlesPDF(std::move(pdf_num));
1101 particle->setPCA0PDF(std::move(pdf_pca0));
1102 particle->setPCA1PDF(std::move(pdf_pca1));
1103 particle->setPCA2PDF(std::move(pdf_pca2));
1104 particle->setPCA3PDF(std::move(pdf_pca3));
1105 particle->setPCA4PDF(std::move(pdf_pca4));
1106
1107 // (8.) set some additional particle and simulation properties
1108 const double restMass = m_particleDataTable->particle(std::abs(pdg))->mass();
1109 minEnergy = ( minEnergy > restMass ) ? minEnergy : restMass;
1110 particle->setMinEnergy(minEnergy);
1111 particle->setMaxNumParticles(maxNumParticles);
1112 particle->setNumParticlesFactor(numParticlesFactor);
1113 particle->setEnergyFactor(energyFactor);
1114 particle->setPosAngleFactor(posAngleFactor);
1115 particle->setMomAngleFactor(momAngleFactor);
1116
1117 // (9.) insert this PunchThroughParticle instance into the std::map class member
1118 m_particles[pdg] = particle;
1119
1120 return StatusCode::SUCCESS;
1121}
std::unique_ptr< ISF::PDFcreator > readLookuptablePDF(int pdgID, const std::string &folderName)
reads out the lookuptable for the given type of particle

Member Data Documentation

◆ m_beamPipe

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

beam pipe radius

Definition at line 207 of file PunchThroughTool.h.

207{this, "BeamPipeRadius", 500.};

◆ m_correlatedParticle

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

Definition at line 187 of file PunchThroughTool.h.

187{this, "CorrelatedParticle", {}, "holds the pdg of the correlated particle for each given pdg"};

◆ m_doAntiParticles

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

Definition at line 186 of file PunchThroughTool.h.

186{this, "DoAntiParticles", {}, "vector of bools to determine if anti-particles are created for each punch-through particle type"};

◆ m_energyFactor

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

Definition at line 195 of file PunchThroughTool.h.

195{this, "EnergyFactor", {}, "scale the energy of the punch-through particles"};

◆ m_energyPoints

std::vector<double> ISF::PunchThroughTool::m_energyPoints
private

energy and eta points in param

Definition at line 154 of file PunchThroughTool.h.

◆ m_envDefSvc

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

Definition at line 204 of file PunchThroughTool.h.

204{this, "EnvelopeDefSvc", "AtlasGeometry_EnvelopeDefSvc"};

◆ m_etaPoints

std::vector<double> ISF::PunchThroughTool::m_etaPoints
private

Definition at line 155 of file PunchThroughTool.h.

◆ m_fileLookupTable

TFile* ISF::PunchThroughTool::m_fileLookupTable {nullptr}
private

ROOT objects.

the punch-through lookup table file

Definition at line 167 of file PunchThroughTool.h.

167{nullptr};

◆ m_filenameInverseCDF

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

Definition at line 178 of file PunchThroughTool.h.

178{this, "FilenameInverseCdf", "", "holds the filename of inverse quantile transformer config"};

◆ m_filenameInversePCA

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

Definition at line 179 of file PunchThroughTool.h.

179{this, "FilenameInversePca", "", "holds the filename of inverse PCA config"};

◆ m_filenameLookupTable

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

Properties.

Definition at line 177 of file PunchThroughTool.h.

177{this, "FilenameLookupTable", "CaloPunchThroughParametrisation.root", "holds the filename of the lookup table"};

◆ m_fullCorrEnergy

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

Definition at line 189 of file PunchThroughTool.h.

189{this, "FullCorrelationEnergy", {}, "holds the energy threshold above which a particle correlation is fully developed"};

◆ m_geoIDSvc

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

Definition at line 203 of file PunchThroughTool.h.

203{this, "GeoIDSvc", "ISF::GeoIDSvc"};

◆ m_initiatorsEtaRange

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

Definition at line 184 of file PunchThroughTool.h.

184{this, "InitiatorsEtaRange", {}, "vector of min and max abs eta range to allow punch through initiators"};

◆ m_initiatorsMinEnergy

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

Definition at line 183 of file PunchThroughTool.h.

183{this, "InitiatorsMinEnergy", {}, "vector of punch-through initiator min energies to create punch through"};

◆ m_inverse_PCA_matrix

std::vector<std::vector<std::vector<double> > > ISF::PunchThroughTool::m_inverse_PCA_matrix
private

pca vectors

Definition at line 210 of file PunchThroughTool.h.

◆ m_maxNumParticles

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

Definition at line 193 of file PunchThroughTool.h.

193{this, "MaxNumParticles", {}, "maximum number of punch-through particles for each particle type"};

◆ m_minCorrEnergy

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

Definition at line 188 of file PunchThroughTool.h.

188{this, "MinCorrelationEnergy", {}, "holds the energy threshold below which no particle correlation is computed"};

◆ m_minEnergy

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

Definition at line 192 of file PunchThroughTool.h.

192{this, "MinEnergy", {}, "punch-through particles minimum energies"};

◆ m_momAngleFactor

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

Definition at line 191 of file PunchThroughTool.h.

191{this, "ScaleMomDeflectionAngles", {}, "tuning parameter to scale the momentum deflection angles"};

◆ m_numParticlesFactor

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

Definition at line 194 of file PunchThroughTool.h.

194{this, "NumParticlesFactor", {}, "scale the number of punch-through particles"};

◆ m_particleDataTable

const HepPDT::ParticleDataTable* ISF::PunchThroughTool::m_particleDataTable {nullptr}
private

ParticleDataTable needed to get connection pdg_code <-> charge.

Definition at line 164 of file PunchThroughTool.h.

164{nullptr};

◆ m_particlePropSvc

ServiceHandle<IPartPropSvc> ISF::PunchThroughTool::m_particlePropSvc {this, "PartPropSvc", "PartPropSvc", "particle properties svc"}
private

Definition at line 202 of file PunchThroughTool.h.

202{this, "PartPropSvc", "PartPropSvc", "particle properties svc"};

◆ m_particles

std::map<int, PunchThroughParticle*> ISF::PunchThroughTool::m_particles
private

needed to create punch-through particles with the right distributions

store all punch-through information for each particle id

Definition at line 170 of file PunchThroughTool.h.

◆ m_PCA_means

std::vector<std::vector<double> > ISF::PunchThroughTool::m_PCA_means
private

Definition at line 211 of file PunchThroughTool.h.

◆ m_pdgInitiators

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

Definition at line 182 of file PunchThroughTool.h.

182{this, "PunchThroughInitiators", {}, "vector of punch-through initiator pgds"};

◆ m_posAngleFactor

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

Definition at line 190 of file PunchThroughTool.h.

190{this, "ScalePosDeflectionAngles", {}, "tuning parameter to scale the position deflection angles"};

◆ m_punchThroughClassifier

PublicToolHandle<IPunchThroughClassifier> ISF::PunchThroughTool::m_punchThroughClassifier {this, "PunchThroughClassifier", "ISF_PunchThroughClassifier", ""}
private

Definition at line 181 of file PunchThroughTool.h.

181{this, "PunchThroughClassifier", "ISF_PunchThroughClassifier", ""};

◆ m_punchThroughParticles

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

Definition at line 185 of file PunchThroughTool.h.

185{this, "PunchThroughParticles", {}, "vector of pdgs of the particles produced in punch-throughs"};

◆ m_R1

double ISF::PunchThroughTool::m_R1 {0.}
private

calo-MS borders

Definition at line 158 of file PunchThroughTool.h.

158{0.};

◆ m_R2

double ISF::PunchThroughTool::m_R2 {0.}
private

Definition at line 159 of file PunchThroughTool.h.

159{0.};

◆ m_variable0_inverse_cdf

std::vector<std::map<double, double> > ISF::PunchThroughTool::m_variable0_inverse_cdf
private

(vector of map) for CDF mappings

Definition at line 218 of file PunchThroughTool.h.

◆ m_variable1_inverse_cdf

std::vector<std::map<double, double> > ISF::PunchThroughTool::m_variable1_inverse_cdf
private

Definition at line 219 of file PunchThroughTool.h.

◆ m_variable2_inverse_cdf

std::vector<std::map<double, double> > ISF::PunchThroughTool::m_variable2_inverse_cdf
private

Definition at line 220 of file PunchThroughTool.h.

◆ m_variable3_inverse_cdf

std::vector<std::map<double, double> > ISF::PunchThroughTool::m_variable3_inverse_cdf
private

Definition at line 221 of file PunchThroughTool.h.

◆ m_variable4_inverse_cdf

std::vector<std::map<double, double> > ISF::PunchThroughTool::m_variable4_inverse_cdf
private

Definition at line 222 of file PunchThroughTool.h.

◆ m_xml_info_cdf

std::vector<InfoMap> ISF::PunchThroughTool::m_xml_info_cdf
private

Definition at line 215 of file PunchThroughTool.h.

◆ m_xml_info_pca

std::vector<InfoMap> ISF::PunchThroughTool::m_xml_info_pca
private

infoMaps

Definition at line 214 of file PunchThroughTool.h.

◆ m_z1

double ISF::PunchThroughTool::m_z1 {0.}
private

Definition at line 160 of file PunchThroughTool.h.

160{0.};

◆ m_z2

double ISF::PunchThroughTool::m_z2 {0.}
private

Definition at line 161 of file PunchThroughTool.h.

161{0.};

The documentation for this class was generated from the following files: