PunchThroughTool.h

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
 
#ifndef ISF_PUNCHTHROUGHTOOLS_SRC_PUNCHTHROUGHTOOL_H
#define ISF_PUNCHTHROUGHTOOLS_SRC_PUNCHTHROUGHTOOL_H 1
 
#include "ISF_FastCaloSimInterfaces/IPunchThroughTool.h"
#include <string>
 
#include "ISF_FastCaloSimInterfaces/IPunchThroughClassifier.h"
 
// Athena Base
#include "AthenaBaseComps/AthAlgTool.h"
 
//Geometry
#include "SubDetectorEnvelopes/IEnvelopeDefSvc.h"
 
#include "ISF_Interfaces/IGeoIDSvc.h"
 
// Gaudi & StoreGate
#include "GaudiKernel/IPartPropSvc.h"
 
#include "GeoPrimitives/GeoPrimitives.h"
 
#include "ISF_Event/ISFParticleContainer.h"
 
#include "AtlasHepMC/GenEvent_fwd.h"
 
//libXML
#include <libxml/xmlmemory.h>
#include <libxml/parser.h>
#include <libxml/tree.h>
#include <libxml/xmlreader.h>
#include <libxml/xpath.h>
#include <libxml/xpathInternals.h>
 
/*-------------------------------------------------------------------------
 *  Forward declarations
 *-------------------------------------------------------------------------*/
 
class TFile;
 
namespace HepPDT {
  class ParticleDataTable;
}
 
namespace ISF {
  class PunchThroughParticle;
  class PDFcreator;
}
 
namespace ISF {
 
  class PunchThroughTool : public extends<AthAlgTool, IPunchThroughTool>
  {
  public:
    PunchThroughTool(const std::string&,const std::string&,const IInterface*);

    virtual ~PunchThroughTool () = default;

    virtual StatusCode initialize();
    virtual StatusCode finalize  ();
    const ISF::ISFParticleVector* computePunchThroughParticles(const ISF::ISFParticle &isfp, const TFCSSimulationState& simulstate, CLHEP::HepRandomEngine* rndmEngine) const;
 
  private:
    /*---------------------------------------------------------------------
     *  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.);
    StatusCode registerCorrelation(int pdgID1, int pdgID2,double minCorrEnergy = 0., double fullCorrEnergy = 0.);

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

    int getAllParticles(const ISF::ISFParticle &isfp, ISFParticleVector& isfpCont, CLHEP::HepRandomEngine* rndmEngine, int pdg, double interpEnergy, double interpEta, int numParticles = -1) const;

    int getCorrelatedParticles(const ISF::ISFParticle &isfp, ISFParticleVector& isfpCont, int doPdg, int corrParticles, CLHEP::HepRandomEngine* rndmEngine, double interpEnergy, double interpEta) const;

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

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

    double getFloatAfterPatternInStr(const char *str, const char *pattern);
    Amg::Vector3D propagator(double theta, double phi) const;
 
    //apply the inverse PCA transform
    std::vector<double> inversePCA(int pcaCdfIterator, std::vector<double> &variables) const;
 
    //apply the inverse CDF trainsform
    static double inverseCdfTransform(double variable, const std::map<double, double>& inverse_cdf_map) ;
 
    //dot product between matrix and vector, used to inverse PCA
    static std::vector<double> dotProduct(const std::vector<std::vector<double>> &m, const std::vector<double> &v) ;
 
    //returns normal cdf based on normal gaussian value
    static double normal_cdf(double x) ;
 
    //apply energy interpolation
    double interpolateEnergy(const double &energy, CLHEP::HepRandomEngine* rndmEngine) const;
 
    //apply eta interpolation
    double interpolateEta(const double &eta, CLHEP::HepRandomEngine* rndmEngine) const;
 
    //get the infoMap from xml file based on the xmlpathname and also name of mainNode
    std::vector<std::map<std::string,std::string>> getInfoMap(const std::string& mainNode, const std::string &xmlFilePath);
 
    //decide the pca / cdf part to read based on pdgId and eta
    int passedParamIterator(int pid, double eta, const std::vector<std::map<std::string,std::string>> &mapvect) const;
 
    //load inverse quantile transformer from XML
    StatusCode initializeInverseCDF(const std::string & quantileTransformerConfigFile);
 
    //get CDF mapping for individual XML node
    static std::map<double, double> getVariableCDFmappings(xmlNodePtr& nodeParent);
 
    //load inverse PCA from XML
    StatusCode initializeInversePCA(const std::string & inversePCAConfigFile);
 
 
    /*---------------------------------------------------------------------
     *  Private members
     *---------------------------------------------------------------------*/

     std::vector<double>                 m_energyPoints;
     std::vector<double>                 m_etaPoints;

    double                               m_R1{0.};
    double                               m_R2{0.};
    double                               m_z1{0.};
    double                               m_z2{0.};

    const HepPDT::ParticleDataTable*    m_particleDataTable{nullptr};

    TFile*                              m_fileLookupTable{nullptr};   

    std::map<int, PunchThroughParticle*> m_particles;       
 
    /*---------------------------------------------------------------------
     *  Properties
     *---------------------------------------------------------------------*/

    StringProperty m_filenameLookupTable{this, "FilenameLookupTable", "CaloPunchThroughParametrisation.root", "holds the filename of the lookup table"};
    StringProperty m_filenameInverseCDF{this, "FilenameInverseCdf", "", "holds the filename of inverse quantile transformer config"};
    StringProperty m_filenameInversePCA{this, "FilenameInversePca", "",  "holds the filename of inverse PCA config"};
 
    PublicToolHandle<IPunchThroughClassifier> m_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"};
 
    
 
    /*---------------------------------------------------------------------
     *  ServiceHandles
     *---------------------------------------------------------------------*/
    ServiceHandle<IPartPropSvc>          m_particlePropSvc{this, "PartPropSvc", "PartPropSvc", "particle properties svc"};
    ServiceHandle<IGeoIDSvc>             m_geoIDSvc{this, "GeoIDSvc", "ISF::GeoIDSvc"};
    ServiceHandle<IEnvelopeDefSvc>       m_envDefSvc{this, "EnvelopeDefSvc", "AtlasGeometry_EnvelopeDefSvc"};

    DoubleProperty                                  m_beamPipe{this, "BeamPipeRadius", 500.};

    std::vector<std::vector<std::vector<double>>> m_inverse_PCA_matrix;
    std::vector<std::vector<double>> m_PCA_means;

    std::vector<std::map<std::string, std::string>> m_xml_info_pca;
    std::vector<std::map<std::string, std::string>> m_xml_info_cdf;

    std::vector<std::map<double, double>>  m_variable0_inverse_cdf;
    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;
  };
}
 
#endif