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::shared_ptr< GenDatam_gendata {nullptr}
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.}
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< 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 46 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 282 of file PunchThroughTool.cxx.

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

◆ 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 1250 of file PunchThroughTool.cxx.

1252{
1253 // the intersection point with Calo-MS surface
1254
1256
1257 // set up the real punch-through particle at this position
1258 // set up the momentum vector of this particle as a GlobalMomentum
1259 // by using the given energy and mass of the particle and also using
1260 // the given theta and phi
1261
1263 double mass = m_gendata->particleMass(std::abs(pdg)).value();
1264 Amg::setRThetaPhi( mom, std::sqrt(energy*energy - mass*mass), momTheta, momPhi);
1265 ATH_MSG_DEBUG("setRThetaPhi pre input parameters: energy = "<< energy <<" mass = "<< mass);
1266 ATH_MSG_DEBUG("setRThetaPhi input parameters: std::sqrt(energy*energy - mass*mass) = "<< std::sqrt(energy*energy - mass*mass) <<" momTheta = "<< momTheta <<" momPhi = "<< momPhi);
1267
1268
1269 double charge = MC::charge(std::abs(pdg));
1270 // since the PDT table only has abs(PID) values for the charge
1271 charge *= (pdg > 0.) ? 1. : -1.;
1272
1273 const double pTime = 0;
1274 const int status = 1 + HepMC::SIM_STATUS_THRESHOLD;
1275 const int id = HepMC::UNDEFINED_ID;
1276 // NB we are not considering the possibility that the punch-through
1277 // particle is the incoming particle having survived an interaction.
1278 ISF::ISFParticle* finalPar = new ISF::ISFParticle ( pos, mom, mass, charge, pdg, status, pTime, isfp, id);
1280
1281 // return the punch-through particle
1282 return finalPar;
1283}
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
std::shared_ptr< GenData > m_gendata
Amg::Vector3D propagator(double theta, double phi) const
get particle through the calorimeter
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...
double charge(const T &p)
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 723 of file PunchThroughTool.cxx.

724{
725 std::vector<double> result;
726 result.reserve(m.size());
727 for (const auto& r : m){
728 result.push_back(std::inner_product(v.begin(), v.end(), r.begin(), 0.0));
729 }
730
731 return result;
732}
int r
Definition globals.cxx:22

◆ finalize()

StatusCode ISF::PunchThroughTool::finalize ( )
virtual

AlgTool finalize method.

Definition at line 265 of file PunchThroughTool.cxx.

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

◆ 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 442 of file PunchThroughTool.cxx.

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

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

1291{
1292 double num = 0.;
1293
1294 const std::string_view str( cstr);
1295 const std::string_view pattern( cpattern);
1296 const size_t pos = str.find(pattern);
1297
1298 if ( pos == std::string::npos)
1299 {
1300 ATH_MSG_WARNING("[ punchthrough ] unable to retrieve floating point number from string");
1301 return -999999999999.;
1302 }
1303 const std::string_view substring = str.substr(pos+pattern.length());
1304 std::from_chars(substring.data(), substring.data() + substring.size(), num);
1305 return num;
1306}
#define ATH_MSG_WARNING(x)

◆ getInfoMap()

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

Definition at line 791 of file PunchThroughTool.cxx.

792{
793 std::vector<InfoMap> xml_info;
794
795 for (const XMLCoreNode* node : doc.get_children (mainNode + "/info/item")) {
796 xml_info.emplace_back(*node);
797 }
798
799 return xml_info;
800}

◆ 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 586 of file PunchThroughTool.cxx.

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

879 {
880
881 std::map<double, double> mappings;
882
883 if (node) {
884 for (const XMLCoreNode* child : node->get_children ("CDFmap")) {
885 double ref = child->get_double_attrib ("ref");
886 double quant = child->get_double_attrib ("quant");
887 mappings[ref] = quant;
888 }
889 }
890
891 return mappings;
892}
const std::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 m_gendata = std::make_shared<GenData>();
130
131 // Geometry identifier service
132 if ( !m_geoIDSvc.empty() && m_geoIDSvc.retrieve().isFailure())
133 {
134 ATH_MSG_FATAL ( "[ punchthrough ] Could not retrieve GeometryIdentifier Service. Abort");
135 return StatusCode::FAILURE;
136 }
137
138 //envelope definition service
139 if (m_envDefSvc.retrieve().isFailure() )
140 {
141 ATH_MSG_ERROR( "[ punchthrough ] Could not retrieve " << m_envDefSvc );
142 return StatusCode::FAILURE;
143 }
144
145 //--------------------------------------------------------------------------------
146 // register all the punch-through particles which will be simulated
147 for ( unsigned int num = 0; num < m_punchThroughParticles.size(); num++ )
148 {
149 const int pdg = m_punchThroughParticles[num];
150 // if no information is given on the creation of anti-particles -> do not simulate anti-particles
151 const bool doAnti = ( num < m_doAntiParticles.size() ) ? m_doAntiParticles[num] : false;
152 // if no information is given on the minimum energy -> take 50. MeV as default
153 const double minEnergy = ( num < m_minEnergy.size() ) ? m_minEnergy[num] : 50.;
154 // if no information is given on the maximum number of punch-through particles -> take -1 as default
155 const int maxNum = ( num < m_minEnergy.size() ) ? m_maxNumParticles[num] : -1;
156 // if no information is given on the scale factor for the number of particles -> take 1. as defaulft
157 const double numFactor = ( num < m_numParticlesFactor.size() ) ? m_numParticlesFactor[num] : 1.;
158 // if no information is given on the position angle factor -> take 1.
159 const double posAngleFactor = ( num < m_posAngleFactor.size() ) ? m_posAngleFactor[num] : 1.;
160 // if no information is given on the momentum angle factor -> take 1.
161 const double momAngleFactor = ( num < m_momAngleFactor.size() ) ? m_momAngleFactor[num] : 1.;
162 // if no information is given on the scale factor for the energy -> take 1. as default
163 const double energyFactor = ( num < m_energyFactor.size() ) ? m_energyFactor[num] : 1.;
164
165 // register the particle
166 ATH_MSG_VERBOSE("[ punchthrough ] registering punch-through particle type with pdg = " << pdg );
167 if ( registerParticle( pdg, doAnti, minEnergy, maxNum, numFactor, energyFactor, posAngleFactor, momAngleFactor )
168 != StatusCode::SUCCESS)
169 {
170 ATH_MSG_ERROR("[ punchthrough ] unable to register punch-through particle type with pdg = " << pdg);
171 }
172 }
173
174 // TODO: implement punch-through parameters for different m_pdgInitiators
175 // currently m_pdgInitiators is only used to filter out particles
176
177 // check if more correlations were given than particle types were registered
178 unsigned int numCorrelations = m_correlatedParticle.size();
179 if ( numCorrelations > m_punchThroughParticles.size() )
180 {
181 ATH_MSG_WARNING("[ punchthrough ] more punch-through particle correlations are given, than punch-through particle types are registered (skipping the last ones)");
182 numCorrelations = m_punchThroughParticles.size();
183 }
184
185 // now register correlation between particles
186 for ( unsigned int num = 0; num < numCorrelations; num++ )
187 {
188 const int pdg1 = m_punchThroughParticles[num];
189 const int pdg2 = m_correlatedParticle[num];
190 const double fullCorrEnergy = ( num < m_fullCorrEnergy.size() ) ? m_fullCorrEnergy[num] : 0.;
191 const double minCorrEnergy = ( num < m_minCorrEnergy.size() ) ? m_minCorrEnergy[num] : 0.;
192
193 // if correlatedParticle==0 is given -> no correlation
194 if ( ! pdg2) continue;
195 // register it
196 if ( registerCorrelation(pdg1, pdg2, minCorrEnergy, fullCorrEnergy) != StatusCode::SUCCESS )
197 {
198 ATH_MSG_ERROR("[ punchthrough ] unable to register punch-through particle correlation for pdg1=" << pdg1 << " pdg2=" << pdg2 );
199 }
200 }
201
202 // get the calo-MS border coordinates. Look at calo and MS geometry definitions, if same R and Z -> boundary surface
203
204 const RZPairVector* rzMS = &(m_envDefSvc->getMuonRZBoundary());
205 const RZPairVector* rzCalo = &(m_envDefSvc->getCaloRZBoundary());
206
207 bool found1, found2;
208 found1=false; found2=false;
209
210 for ( size_t i=0; i<rzCalo->size();++i)
211 {
212 const double r_tempCalo = rzCalo->at(i).first;
213 const double z_tempCalo = rzCalo->at(i).second;
214
215 if (r_tempCalo> m_beamPipe)
216 {
217 for ( size_t j=0; j<rzMS->size();++j)
218 {
219 const double r_tempMS =rzMS->at(j).first;
220 const double z_tempMS =rzMS->at(j).second;
221
222 if (r_tempCalo==r_tempMS && z_tempCalo==z_tempMS && found1==false )
223 {
224 m_R1=r_tempMS;
225 m_z1=std::fabs(z_tempMS);
226 found1=true;
227 continue;
228 }
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) ATH_MSG_ERROR ("first coordinate of calo-MS border not found");
243 if (found2 == false) ATH_MSG_ERROR ("second coordinate of calo-MS border not found; first one is: R1 ="<<m_R1<<" z1 ="<<m_z1);
244
245 //now order the found values
246 double r_temp, z_temp;
247 if (m_R1>m_R2) { r_temp=m_R1; m_R1=m_R2; m_R2=r_temp; } //m_R1 - smaller one
248 if (m_z1<m_z2) { z_temp=m_z1; m_z1=m_z2; m_z2=z_temp; } //m_z1 - bigger one
249
250 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 );
251 else ATH_MSG_DEBUG ("calo-MS boundary coordinates: R1="<<m_R1<<" R2 = "<<m_R2<<" z1="<<m_z1<<" z2= "<<m_z2);
252
253 // close the file with the lookuptable
254 m_fileLookupTable->Close();
255
256 ATH_MSG_INFO( "punchthrough initialization is successful" );
257 return StatusCode::SUCCESS;
258}
#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
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 853 of file PunchThroughTool.cxx.

853 {
854 //parse xml that contains config for inverse CDF for each of punch through particle kinematics
855 XMLCoreParser p;
856 std::unique_ptr<XMLCoreNode> doc = p.parse (inverseCdfConfigFile);
857
858 ATH_MSG_INFO( "[ punchthrough ] Loading inverse CDF: " << inverseCdfConfigFile);
859
860 //check info first
861 m_xml_info_cdf = getInfoMap("CDFMappings", *doc);
862
863 for (unsigned int i = 0; i < m_xml_info_cdf.size(); i++) {
864 ATH_MSG_DEBUG( "[PunchThroughG4Tool] m_xml_info_cdf[" << i << "].name = " << m_xml_info_cdf[i].name);
865
866 if (const XMLCoreNode* node = doc->get_child ("CDFMappings/" + m_xml_info_cdf[i].name)) {
867 m_variable0_inverse_cdf.push_back(getVariableCDFmappings (node->get_child("variable0")));
868 m_variable1_inverse_cdf.push_back(getVariableCDFmappings (node->get_child("variable1")));
869 m_variable2_inverse_cdf.push_back(getVariableCDFmappings (node->get_child("variable2")));
870 m_variable3_inverse_cdf.push_back(getVariableCDFmappings (node->get_child("variable3")));
871 m_variable4_inverse_cdf.push_back(getVariableCDFmappings (node->get_child("variable4")));
872 }
873 }
874
875 return StatusCode::SUCCESS;
876}
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 811 of file PunchThroughTool.cxx.

811 {
812
813 XMLCoreParser p;
814 std::unique_ptr<XMLCoreNode> doc = p.parse (inversePCAConfigFile);
815
816 ATH_MSG_INFO( "[ punchthrough ] Loading inversePCA: " << inversePCAConfigFile);
817
818 //check info first
819 m_xml_info_pca = getInfoMap("PCAinverse", *doc);
820
821 auto getRow = [] (const XMLCoreNode* node,
822 const std::string& attrib,
823 int imax)
824 -> std::vector<double>
825 {
826 std::vector<double> row;
827 for (int i = 0; i <= imax; ++i) {
828 // Dynamically create property name
829 std::string propName = attrib + std::to_string(i);
830 row.push_back (node->get_double_attrib (propName));
831 }
832 return row;
833 };
834
835 for (unsigned int i = 0; i < m_xml_info_pca.size(); i++) {
836 std::vector<std::vector<double>> PCA_matrix;
837
838 if (const XMLCoreNode* node = doc->get_child ("PCAinverse/" + m_xml_info_pca[i].name)) {
839 for (const XMLCoreNode* matrix : node->get_children ("PCAmatrix")) {
840 PCA_matrix.push_back(getRow(matrix, "comp_", 4));
841 }
842 for (const XMLCoreNode* means : node->get_children ("PCAmeans")) {
843 m_PCA_means.push_back(getRow(means, "mean_", 4));
844 }
845 }
846
847 m_inverse_PCA_matrix.push_back(std::move(PCA_matrix));
848 }
849
850 return StatusCode::SUCCESS;
851}
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(std::string_view 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 909 of file PunchThroughTool.cxx.

909 {
910
911 ATH_MSG_DEBUG("[ punchthrough ] interpolating incoming energy: " << energy);
912
913 std::string energyPointsString;
914 for (auto element:m_energyPoints){
915 energyPointsString += std::to_string(element) + " ";
916 }
917
918 ATH_MSG_DEBUG("[ punchthrough ] available energy points: " << energyPointsString);
919
920 auto const upperEnergy = std::upper_bound(m_energyPoints.begin(), m_energyPoints.end(), energy);
921
922 if(upperEnergy == m_etaPoints.end()){ //if no energy greater than input energy, choose greatest energy
923 ATH_MSG_DEBUG("[ punchthrough ] incoming energy > largest energy point, returning greatest energy point: " << m_energyPoints.back());
924 return m_energyPoints.back();
925 }
926 else if(upperEnergy == m_etaPoints.begin()){ //if smallest energy greater than input energy, choose smallest energy
927 ATH_MSG_DEBUG("[ punchthrough ] incoming energy < smallest energy point, returning smallest energy point: " << *upperEnergy);
928 return *upperEnergy;
929 }
930
931 ATH_MSG_DEBUG("[ punchthrough ] energy points upper_bound: "<< *upperEnergy);
932
933 double randomShoot = CLHEP::RandFlat::shoot(rndmEngine);
934
935 ATH_MSG_DEBUG("[ punchthrough ] Shooting random number: "<< randomShoot);
936
937 double midPoint = *std::prev(upperEnergy)*M_SQRT2;
938
939 if(energy < midPoint){ //if energy smaller than mid point in log(energy)
940
941 double distance = std::abs(energy - *std::prev(upperEnergy))/((midPoint) - *std::prev(upperEnergy));
942
943 ATH_MSG_DEBUG( "[ punchthrough ] incoming energy is closest to prev(upper_bound) in log(energy), distance: " << distance );
944
945 if(randomShoot < distance){
946 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot < distance, returning upper_bound " << *upperEnergy );
947 return *upperEnergy;
948 }
949 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot > distance, returning prev(upper_bound) " << *std::prev(upperEnergy) );
950
951 return *std::prev(upperEnergy);
952 }
953 else if(energy > midPoint){ //if energy greater than mid point in log(energy)
954
955 double distance = std::abs(energy - *upperEnergy)/((*upperEnergy - midPoint));
956
957 ATH_MSG_DEBUG( "[ punchthrough ] incoming energy is closest to upper_bound in log(energy), distance: " << distance );
958
959 if(randomShoot < distance){
960 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot < distance, returning prev(upper_bound) " << *std::prev(upperEnergy) );
961 return *std::prev(upperEnergy);
962 }
963 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot > distance, returning upper_bound " << *upperEnergy );
964 return *upperEnergy;
965 }
966
967 return *upperEnergy;
968}
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 970 of file PunchThroughTool.cxx.

970 {
971
972 double absEta = std::abs(eta);
973
974 ATH_MSG_DEBUG("[ punchthrough ] interpolating incoming abs(eta): " << absEta);
975
976 std::string etaPointsString;
977 for (auto element:m_etaPoints){
978 etaPointsString += std::to_string(element) + " ";
979 }
980
981 ATH_MSG_DEBUG("[ punchthrough ] available eta points: " << etaPointsString);
982
983 auto const upperEta = std::upper_bound(m_etaPoints.begin(), m_etaPoints.end(), absEta);
984
985 if(upperEta == m_etaPoints.end()){
986 ATH_MSG_DEBUG("[ punchthrough ] incoming abs(eta) > largest eta point, returning greatest eta point: " << m_etaPoints.back());
987 return m_etaPoints.back();
988 }
989
990
991 ATH_MSG_DEBUG("[ punchthrough ] eta points upper_bound: "<< *upperEta);
992
993 double randomShoot = CLHEP::RandFlat::shoot(rndmEngine);
994
995 ATH_MSG_DEBUG("[ punchthrough ] Shooting random number: "<< randomShoot);
996
997 if(std::abs(absEta - *upperEta) < std::abs(absEta - *std::prev(upperEta))){
998
999 double distance = std::abs(absEta - *upperEta)/((*upperEta - *std::prev(upperEta))/2);
1000
1001 ATH_MSG_DEBUG( "[ punchthrough ] abs(eta) is closer to eta points upper_bound, distance: " << distance );
1002
1003 if(randomShoot > distance){
1004 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot > distance, returning upper_bound " << *upperEta );
1005 return *upperEta;
1006 }
1007
1008 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot < distance, returning prev(upper_bound) " << *std::prev(upperEta) );
1009
1010 return *std::prev(upperEta);
1011 }
1012 else if(std::abs(absEta - *std::prev(upperEta)) < std::abs(absEta - *upperEta)){
1013
1014 if(std::prev(upperEta) == m_etaPoints.begin()){
1015 ATH_MSG_DEBUG( "[ punchthrough ] prev of upper bound is begin, returning that: " << *std::prev(upperEta) );
1016 return *std::prev(upperEta);
1017 }
1018
1019 double distance = std::abs(absEta - *std::prev(upperEta))/((*std::prev(upperEta) - *std::prev(std::prev(upperEta)))/2);
1020
1021 ATH_MSG_DEBUG( "[ punchthrough ] abs(eta) is closer to eta points prev(upper_bound), distance: " << distance );
1022
1023 if(randomShoot > distance){
1024 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot > distance, returning prev(prev(upper_bound)) " << *std::prev(std::prev(upperEta)) );
1025
1026 return *std::prev(std::prev(upperEta));
1027 }
1028 ATH_MSG_DEBUG( "[ punchthrough ] randomShoot < distance, returning prev(upper_bound) " << *std::prev(upperEta) );
1029
1030 return *std::prev(upperEta);
1031 }
1032
1033 return *std::prev(upperEta);
1034}
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 894 of file PunchThroughTool.cxx.

894 {
895
896 double norm_cdf = normal_cdf(variable);
897
898 auto upper = inverse_cdf_map.upper_bound(norm_cdf);
899 auto lower = upper--;
900
901 double m = (upper->second - lower->second)/(upper->first - lower->first);
902 double c = lower->second - m * lower->first;
903 double transformed = m * norm_cdf + c;
904
905 return transformed;
906
907}
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 802 of file PunchThroughTool.cxx.

803{
804 std::vector<double> transformed_variables = dotProduct(m_inverse_PCA_matrix[pcaCdfIterator], variables);
805
806 std::transform (transformed_variables.begin(), transformed_variables.end(), m_PCA_means[pcaCdfIterator].begin(), transformed_variables.begin(), std::plus<double>()); // + means
807
808 return transformed_variables;
809}
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 718 of file PunchThroughTool.cxx.

718 {
719
720 return 0.5 * TMath::Erfc(-x * M_SQRT1_2);
721}
#define x

◆ passedParamIterator()

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

Definition at line 758 of file PunchThroughTool.cxx.

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

◆ propagator()

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

get particle through the calorimeter

Definition at line 1308 of file PunchThroughTool.cxx.

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

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

◆ 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 1119 of file PunchThroughTool.cxx.

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

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

201{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 182 of file PunchThroughTool.h.

182{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 181 of file PunchThroughTool.h.

181{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 190 of file PunchThroughTool.h.

190{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 152 of file PunchThroughTool.h.

◆ m_envDefSvc

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

Definition at line 198 of file PunchThroughTool.h.

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

◆ m_etaPoints

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

Definition at line 153 of file PunchThroughTool.h.

◆ m_fileLookupTable

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

ROOT objects.

the punch-through lookup table file

Definition at line 162 of file PunchThroughTool.h.

162{nullptr};

◆ m_filenameInverseCDF

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

Definition at line 173 of file PunchThroughTool.h.

173{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 174 of file PunchThroughTool.h.

174{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 172 of file PunchThroughTool.h.

172{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 184 of file PunchThroughTool.h.

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

◆ m_gendata

std::shared_ptr<GenData> ISF::PunchThroughTool::m_gendata {nullptr}
private

Definition at line 63 of file PunchThroughTool.h.

63{nullptr};

◆ m_geoIDSvc

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

Definition at line 197 of file PunchThroughTool.h.

197{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 179 of file PunchThroughTool.h.

179{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 178 of file PunchThroughTool.h.

178{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 204 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 188 of file PunchThroughTool.h.

188{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 183 of file PunchThroughTool.h.

183{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 187 of file PunchThroughTool.h.

187{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 186 of file PunchThroughTool.h.

186{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 189 of file PunchThroughTool.h.

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

◆ 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 165 of file PunchThroughTool.h.

◆ m_PCA_means

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

Definition at line 205 of file PunchThroughTool.h.

◆ m_pdgInitiators

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

Definition at line 177 of file PunchThroughTool.h.

177{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 185 of file PunchThroughTool.h.

185{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 176 of file PunchThroughTool.h.

176{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 180 of file PunchThroughTool.h.

180{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 156 of file PunchThroughTool.h.

156{0.};

◆ m_R2

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

Definition at line 157 of file PunchThroughTool.h.

157{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 212 of file PunchThroughTool.h.

◆ m_variable1_inverse_cdf

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

Definition at line 213 of file PunchThroughTool.h.

◆ m_variable2_inverse_cdf

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

Definition at line 214 of file PunchThroughTool.h.

◆ m_variable3_inverse_cdf

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

Definition at line 215 of file PunchThroughTool.h.

◆ m_variable4_inverse_cdf

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

Definition at line 216 of file PunchThroughTool.h.

◆ m_xml_info_cdf

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

Definition at line 209 of file PunchThroughTool.h.

◆ m_xml_info_pca

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

infoMaps

Definition at line 208 of file PunchThroughTool.h.

◆ m_z1

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

Definition at line 158 of file PunchThroughTool.h.

158{0.};

◆ m_z2

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

Definition at line 159 of file PunchThroughTool.h.

159{0.};

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