diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/ISF_PunchThroughTools/PDFcreator.h b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/ISF_PunchThroughTools/PDFcreator.h
new file mode 100644
index 0000000000000000000000000000000000000000..9cd5083a324b11a7d0fa1a0d251ce037db00fd84
--- /dev/null
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/ISF_PunchThroughTools/PDFcreator.h
@@ -0,0 +1,56 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// PDFcreator.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#ifndef ISF_PUNCHTHROUGHTOOLS_SRC_PDFCREATOR_H
+#define ISF_PUNCHTHROUGHTOOLS_SRC_PDFCREATOR_H
+
+// Athena Base
+#include "AthenaKernel/IAtRndmGenSvc.h"
+
+// ROOT includes
+#include "TFile.h"
+#include "TF1.h"
+#include "TH2F.h"
+
+namespace ISF
+{
+ /** @class PDFcreator
+
+ Creates random numbers with a distribution given as ROOT TF1.
+ The TF1 function parameters will be retrieved from a histogram given by addPar.
+
+ @author Elmar Ritsch <Elmar.Ritsch@cern.ch>
+ */
+
+ class PDFcreator
+ {
+
+ public:
+ /** construct the class with a given TF1 and a random engine */
+ PDFcreator(TF1 *pdf, CLHEP::HepRandomEngine *engine):m_randomEngine(engine), m_pdf(pdf), m_randmin(0), m_randmax(0) {} ;
+
+ ~PDFcreator() { };
+
+ /** all following is used to set up the class */
+ void addPar(TH1 *hist) { m_par.push_back(hist); }; //!< get the function's parameter from the given histogram
+ void setRange( TH1 *min, TH1 *max) { m_randmin = min; m_randmax = max; }; //!< get the function's range from the given histograms
+
+ /** get the random value with this methode, by providing the input parameters */
+ double getRand( std::vector<double> inputPar, bool discrete = false, double randMin = 0., double randMax = 0.);
+
+ private:
+ CLHEP::HepRandomEngine *m_randomEngine; //!< Random Engine
+ TF1 *m_pdf; //!< the probability density function
+ TH1 *m_randmin; //!< the fit-functions minimum values
+ TH1 *m_randmax; //!< the fit-functions maximum values
+ std::vector<TH1*> m_par; /*!< the 'histograms' which hold the parameters
+ for the above PDF */
+ };
+}
+
+#endif
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/ISF_PunchThroughTools/PunchThroughParticle.h b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/ISF_PunchThroughTools/PunchThroughParticle.h
new file mode 100644
index 0000000000000000000000000000000000000000..f213e4e925715272301ebffd1350643779a786e3
--- /dev/null
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/ISF_PunchThroughTools/PunchThroughParticle.h
@@ -0,0 +1,112 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// PunchThroughParticle.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#ifndef ISF_PUNCHTHROUGHTOOLS_PUNCHTHROUGHPARTICLE_H
+#define ISF_PUNCHTHROUGHTOOLS_PUNCHTHROUGHPARTICLE_H
+
+// ISF
+#include "PDFcreator.h"
+
+// forward declarations
+class TH2F;
+
+/*-------------------------------------------------------------------------
+ * class PunchThroughParticle
+ *-------------------------------------------------------------------------*/
+namespace ISF
+{
+ /** @class PunchThroughParticle
+
+ This class holds information for different properties of a punch-through
+ particle (energy, theta, phi) distributions
+
+ @author Elmar Ritsch <Elmar.Ritsch@cern.ch>
+ */
+
+ class PunchThroughParticle
+ {
+ public:
+ PunchThroughParticle(int pdg, bool doAnti = false); //!< set this particle's pdg code and if anti-particle should be done or not
+ ~PunchThroughParticle();
+
+ /** set methods */
+ void setMinEnergy(double minEnergy); //!< the minimum energy with which these particles should be created
+ void setMaxNumParticles(int maxNum); //!< the maximum number of particles which will be created
+ void setNumParticlesFactor(double numFactor); //!< to scale the number of punch-through particles
+ void setEnergyFactor(double energyFactor); //!< to scale the energy of created particles
+ void setPosAngleFactor(double momAngleFactor); //!< to scale the position deflection angles
+ void setMomAngleFactor(double momAngleFactor); //!< to scale the momentum deviation of created particles
+ void setNumParticlesPDF(PDFcreator *pdf); //!< set the PDFcreator for the number of exit particles distribution
+ void setCorrelation(int corrPdg, TH2F *histLowE, TH2F *histHighE,
+ double minCorrE = 0., double fullCorrE = 0., double lowE = 0.,
+ double midE = 0., double upperE = 0.); /*!< set the correlated particle type + correlation histograms
+ + full energy correlation threshold
+ + histogram domains histLowE:[lowE,midE] histHighE:[midE,upperE]*/
+ void setExitEnergyPDF(PDFcreator *pdf); //!< set the PDFcreator for the energy distribution
+ void setExitDeltaThetaPDF(PDFcreator *pdf); //!< set the PDFcreator for the deltaTheta distribution
+ void setExitDeltaPhiPDF(PDFcreator *pdf); //!< set the PDFcreator for the deltaPhi distribution
+ void setMomDeltaThetaPDF(PDFcreator *pdf); //!< set the PDFcreator for the momentumDeltaTheta distribution
+ void setMomDeltaPhiPDF(PDFcreator *pdf); //!< set the PDFcreator for the momentumDeltaPhi distribution
+
+
+ /** get-access methods */
+ inline int getId() { return m_pdgId; }; //!< get this particle's pdg code
+ inline bool getdoAnti() { return m_doAnti; }; //!< get if anti-particles should be done or not
+ inline double getMinEnergy() { return m_minEnergy; }; //!< the minimum energy with which these particles should be created
+ inline double getNumParticlesFactor() { return m_numParticlesFactor; }; //!< to scale the number of punch-through particles
+ inline double getEnergyFactor() { return m_energyFactor; }; //!< to scale the energy of created particles
+ inline double getPosAngleFactor() { return m_posAngleFactor; }; //!< to scale the position deviation angles
+ inline double getMomAngleFactor() { return m_momAngleFactor; }; //!< to scale the momentum deviation of created particles
+
+ inline int getMaxNumParticles() { return m_maxNum; }; //!< the maximum number of particles which will be created
+ inline int getCorrelatedPdg() { return m_corrPdg; }; //!< the correlated particle id
+ inline int getMinCorrelationEnergy() { return m_corrMinEnergy; }; //!< minimum energy for correlation
+ inline int getFullCorrelationEnergy() { return m_corrFullEnergy; }; //!< the full correlation energy
+ inline TH2F *getCorrelationLowEHist() { return m_histCorrLowE; }; //!< the correlation histogram (low energies)
+ inline TH2F *getCorrelationHighEHist() { return m_histCorrHighE; }; //!< the correlation histogram (high energies)
+ inline double *getCorrelationHistDomains(){ return m_corrHistDomains; }; //!< correlation histogram domains [lowE,midE,upperE]
+ inline PDFcreator *getNumParticlesPDF() { return m_pdfNumParticles; }; //!< distribution parameters for the number of particles
+ inline PDFcreator *getExitEnergyPDF() { return m_pdfExitEnergy; }; //!< distribution parameters for the exit energy
+ inline PDFcreator *getExitDeltaThetaPDF() { return m_pdfExitDeltaTheta; }; //!< distribution parameters for the exit delta theta
+ inline PDFcreator *getExitDeltaPhiPDF() { return m_pdfExitDeltaPhi; }; //!< distribution parameters for the exit delta phi
+ inline PDFcreator *getMomDeltaThetaPDF() { return m_pdfMomDeltaTheta; }; //!< distribution parameters for the momentum delta theta
+ inline PDFcreator *getMomDeltaPhiPDF() { return m_pdfMomDeltaPhi; }; //!< distribution parameters for the momentum delta phi
+
+ private:
+ int m_pdgId; //!< the pdg-id of this particle
+ bool m_doAnti; //!< do also anti-particles?
+
+ /** some cut-parameters which will be set via python */
+ double m_minEnergy; //!< the minimum energy with which these particles should be created
+ int m_maxNum; //!< the maximum number of particles which will be created
+
+ /** some tuning-parameters which will be set via python */
+ double m_numParticlesFactor; //!< scale the number of particles created
+ double m_energyFactor; //!< scale the energy of this particle type
+ double m_posAngleFactor; //!< scale the position deflection angles
+ double m_momAngleFactor; //!< scale the momentum deviation
+
+ /** all following stores the right distributions for all properties of the punch-through particles */
+ int m_corrPdg; //!< correlation to any other punch-through particle type
+ double m_corrMinEnergy; //!< below this energy threshold, no particle correlation is computed
+ double m_corrFullEnergy; /*!< holds the energy threshold above which
+ a particle correlation is fully developed */
+
+ TH2F* m_histCorrLowE; //!< low energy correlation histogram (x:this particle, y:the correlated particle)
+ TH2F* m_histCorrHighE; //!< high energy correlation histogram (x:this particle, y:the correlated particle)
+ double *m_corrHistDomains; //!< correlation histogram domains
+ PDFcreator* m_pdfNumParticles; //!< number of punch-through particles
+ PDFcreator* m_pdfExitEnergy; //!< energy of punch-through particles
+ PDFcreator* m_pdfExitDeltaTheta; //!< theta deviation of punch-through particles
+ PDFcreator* m_pdfExitDeltaPhi; //!< phi deviation of punch-through particles
+ PDFcreator* m_pdfMomDeltaTheta; //!< delta theta angle of punch-through particle momentum
+ PDFcreator* m_pdfMomDeltaPhi; //!< delta phi angle of punch-through particle momentum
+ };
+}
+
+#endif
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/ISF_PunchThroughTools/PunchThroughTool.h b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/ISF_PunchThroughTools/PunchThroughTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..62150add0eaacca5b650b96f33798945091fc532
--- /dev/null
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/ISF_PunchThroughTools/PunchThroughTool.h
@@ -0,0 +1,183 @@
+/*
+ Copyright (C) 2002-2017 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>
+
+// Athena Base
+#include "AthenaBaseComps/AthAlgTool.h"
+
+#include "BarcodeInterfaces/Barcode.h"
+#include "BarcodeInterfaces/PhysicsProcessCode.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+
+/*-------------------------------------------------------------------------
+ * Forward declarations
+ *-------------------------------------------------------------------------*/
+
+class IPartPropSvc;
+class IAtRndmGenSvc;
+class TFile;
+class IEnvelopeDefSvc;
+
+namespace CLHEP {
+ class HepRandomEngine;
+}
+
+namespace HepPDT {
+ class ParticleDataTable;
+}
+
+namespace HepMC {
+ class GenEvent;
+}
+
+namespace Barcode {
+ class IBarcodeSvc;
+}
+
+namespace ISF {
+
+ class PunchThroughParticle;
+ class PDFcreator;
+ class IGeoIDSvc;
+}
+
+namespace ISF {
+
+ class PunchThroughTool : virtual public IPunchThroughTool,
+ public AthAlgTool
+ {
+ public:
+ /** Constructor */
+ PunchThroughTool(const std::string&,const std::string&,const IInterface*);
+
+ /** Destructor */
+ virtual ~PunchThroughTool ();
+
+ /** AlgTool initialize method */
+ virtual StatusCode initialize();
+ /** AlgTool finalize method */
+ virtual StatusCode finalize ();
+ /** interface function: fill a vector with the punch-through particles */
+ const ISF::ISFParticleContainer* computePunchThroughParticles(const ISF::ISFParticle &isfp) const;
+
+ private:
+ /*---------------------------------------------------------------------
+ * Private member functions
+ *---------------------------------------------------------------------*/
+ /** registers a type of punch-through particles which will be simulated */
+ 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.);
+ /** register a correlation for the two given types of punch-through particles
+ with a given energy threshold above which we will have full correlation */
+ StatusCode registerCorrelation(int pdgID1, int pdgID2,double minCorrEnergy = 0., double fullCorrEnergy = 0.);
+
+ /** reads out the lookuptable for the given type of particle */
+ PDFcreator *readLookuptablePDF(int pdgID, std::string folderName);
+
+ /** 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 */
+ int getAllParticles(int pdg, int numParticles = -1) 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 */
+ int getCorrelatedParticles(int doPdg, int corrParticles ) const;
+
+ /** create exactly one punch-through particle with the given pdg and the given max energy */
+ ISF::ISFParticle *getOneParticle(int pdg, double maxEnergy) const;
+
+ /** create a ISF Particle state at the MS entrace containing a particle with the given properties */
+ ISF::ISFParticle *createExitPs(int PDGcode, double energy, double theta, double phi, double momTheta, double momPhi) const;
+
+ /** get the floating point number in a string, after the given pattern */
+ double getFloatAfterPatternInStr(const char *str, const char *pattern);
+ /** get particle through the calorimeter */
+ Amg::Vector3D propagator(double theta, double phi) const;
+
+ /*---------------------------------------------------------------------
+ * Private members
+ *---------------------------------------------------------------------*/
+
+ /** initial particle properties */
+ mutable const ISF::ISFParticle* m_initPs; //!< the incoming particle
+ mutable double m_initEnergy; //!< the incoming particle's energy
+ mutable double m_initEta; //!< the incoming particle's eta
+ mutable double m_initTheta; //!< the incoming particle's theta
+ mutable double m_initPhi; //!< the incoming particle's phi
+
+ /** calo-MS borders */
+ double R1, R2, z1, z2;
+
+ /** the returned vector of ISFParticles */
+ mutable ISF::ISFParticleContainer *m_isfpCont;
+
+ /** parent event */
+ mutable HepMC::GenEvent* m_parentGenEvt; //!< all newly created particles/vertices will have this common parent
+
+ /** ParticleDataTable needed to get connection pdg_code <-> charge */
+ const HepPDT::ParticleDataTable* m_particleDataTable;
+
+ /** random generator service */
+ ServiceHandle<IAtRndmGenSvc> m_randomSvc; //!< Random Svc
+ std::string m_randomEngineName; //!< Name of the random number stream
+ CLHEP::HepRandomEngine* m_randomEngine; //!< Random Engine
+
+ /** ROOT objects */
+ TFile* m_fileLookupTable; //!< the punch-through lookup table file
+ /** general information of the punch-through particles which will be created */
+ mutable std::map<int, bool> m_doAntiParticleMap; /*!< stores information, if anti-particles are
+ created for any given PDG id */
+ /** needed to create punch-through particles with the right distributions */
+ mutable std::map<int, PunchThroughParticle*> m_particles; //!< store all punch-through information for each particle id
+
+ /** barcode steering */
+
+ mutable Barcode::PhysicsProcessCode m_processCode;
+ mutable Barcode::ParticleBarcode m_primBC;
+ mutable Barcode::ParticleBarcode m_secBC;
+
+ /*---------------------------------------------------------------------
+ * Properties
+ *---------------------------------------------------------------------*/
+
+ /** Properties */
+ std::string m_filenameLookupTable; //!< holds the filename of the lookup table (property)
+ std::vector<int> m_pdgInitiators; //!< vector of punch-through initiator pgds
+ std::vector<int> m_punchThroughParticles; //!< vector of pdgs of the particles produced in punch-throughs
+ std::vector<bool> m_doAntiParticles; //!< vector of bools to determine if anti-particles are created for each punch-through particle type
+ std::vector<int> m_correlatedParticle; //!< holds the pdg of the correlated particle for each given pdg
+ std::vector<double> m_minCorrEnergy; //!< holds the energy threshold below which no particle correlation is computed
+ std::vector<double> m_fullCorrEnergy; //!< holds the energy threshold above which a particle correlation is fully developed
+ std::vector<double> m_posAngleFactor; //!< tuning parameter to scale the position deflection angles
+ std::vector<double> m_momAngleFactor; //!< tuning parameter to scale the momentum deflection angles
+ std::vector<double> m_minEnergy; //!< punch-through particles minimum energies
+ std::vector<int> m_maxNumParticles; //!< maximum number of punch-through particles for each particle type
+ std::vector<double> m_numParticlesFactor; //!< scale the number of punch-through particles
+ std::vector<double> m_energyFactor; //!< scale the energy of the punch-through particles
+
+ /*---------------------------------------------------------------------
+ * ServiceHandles
+ *---------------------------------------------------------------------*/
+ ServiceHandle<IPartPropSvc> m_particlePropSvc; //!< particle properties svc
+ ServiceHandle<IGeoIDSvc> m_geoIDSvc;
+ ServiceHandle<Barcode::IBarcodeSvc> m_barcodeSvc;
+ ServiceHandle<IEnvelopeDefSvc> m_envDefSvc;
+
+ /** beam pipe radius */
+ double m_beamPipe;
+ };
+}
+
+#endif
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/cmt/requirements b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/cmt/requirements
new file mode 100644
index 0000000000000000000000000000000000000000..2ced03874e05630f45318bd246a8442988e3a7ed
--- /dev/null
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/cmt/requirements
@@ -0,0 +1,49 @@
+package ISF_PunchThroughTools
+
+manager Andreas Salzburger <Andreas.Salzburger@cern.ch>
+manager Elmar Ritsch <Elmar.Ritsch@cern.ch>
+manager Wolfgang Lukas <Wolfgang.Lukas@cern.ch>
+
+public
+use AtlasPolicy AtlasPolicy-*
+use BarcodeInterfaces BarcodeInterfaces-* Simulation/Barcode
+
+########## Control #############################################################
+
+use AthenaKernel AthenaKernel-* Control
+use AthenaBaseComps AthenaBaseComps-* Control
+use AtlasROOT AtlasROOT-* External
+
+use ISF_FastCaloSimInterfaces ISF_FastCaloSimInterfaces-* Simulation/ISF/ISF_FastCaloSim
+
+########## Tracking ############################################################
+
+use GeoPrimitives GeoPrimitives-* DetectorDescription
+
+private
+
+use DataModel DataModel-* Control
+use GaudiInterface GaudiInterface-* External
+use SubDetectorEnvelopes SubDetectorEnvelopes-* AtlasGeometryCommon
+########## External ############################################################
+use AtlasCLHEP AtlasCLHEP-* External
+use AtlasHepMC AtlasHepMC-* External
+
+########## ISF #################################################################
+use ISF_Event ISF_Event-* Simulation/ISF/ISF_Core
+use ISF_Interfaces ISF_Interfaces-* Simulation/ISF/ISF_Core
+
+##############################################
+use HepPDT v* LCG_Interfaces
+use PathResolver PathResolver-* Tools
+
+declare_runtime_extras extras = ../Data/*.root
+
+public
+library ISF_PunchThroughTools *.cxx components/*.cxx
+apply_pattern component_library
+
+# use the following to compile with debug information
+#private
+#macro cppdebugflags '$(cppdebugflags_s)'
+#macro_remove componentshr_linkopts "-Wl,-s"
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PDFcreator.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PDFcreator.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..d16d6a0f335f2ed9253878a467c84a77d98908df
--- /dev/null
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PDFcreator.cxx
@@ -0,0 +1,164 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// PDFcreator.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+// class header
+#include "ISF_PunchThroughTools/PDFcreator.h"
+
+// Random number generators
+#include "CLHEP/Random/RandFlat.h"
+
+// ROOT
+#include "TFile.h"
+#include "TH2F.h"
+#include "TAxis.h"
+#include "TRandom.h"
+
+/*=========================================================================
+* DESCRIPTION OF FUNCTION:
+* ==> see headerfile
+*=======================================================================*/
+double ISF::PDFcreator::getRand(std::vector<double> inputParameters, bool discrete, double randMin, double randMax)
+{
+// Since the histograms which hold the function's fit parameters are binned
+// we have to choose which bin we use for the current energy & eta regime.
+// It's very unlikely to exactly hit a bin's center, therefore we randomly
+// choose either the lower next or the upper next bin.
+
+ int numInputPar = inputParameters.size();
+ TAxis **axis = new TAxis*[numInputPar];
+ Int_t *chosenBin = new Int_t[numInputPar];
+
+
+// get the axis from the first (could be any other too, but the first one
+// always exists) histogram (which holds the function's first fit parameter)
+ axis[0] = m_par[0]->GetXaxis();
+ axis[1] = m_par[0]->GetYaxis();
+
+// loop over all input inputParameters (e.g. eta & energy)
+ for (int inputPar=0; inputPar<numInputPar; inputPar++) {
+ // store the axis for the current inputParameter
+ TAxis *curaxis = axis[inputPar];
+ double curvalue = inputParameters[inputPar];
+
+ // get the id of the bin corresponding to the current value
+ // (use FindFixBin to avoid the axis from being rebinned)
+ Int_t bin = curaxis->FindFixBin(curvalue);
+ // get the center of the closest bin to the input inputParameter
+ double closestCenter = curaxis->GetBinCenter(bin);
+ // get the bins edge closest to the current value
+ // and find out if the next closest bin has id bin+1 or bin-1
+ double closestEdge = (curvalue <= closestCenter) ? curaxis->GetBinLowEdge(bin) : curaxis->GetBinUpEdge(bin);
+ int binDelta = (curvalue <= closestCenter) ? -1 : +1;
+
+ // Calculate the 'normalised distance' between the input-value,
+ // the closest bin's center and this bin's closest edge
+ // distance == 0 means that the input value lies exactly on its closest bin's center.
+ // distance == 0.5 means that the input value lies exactly on the edge to the next bin
+ double distance = (curvalue - closestCenter) / (closestEdge - closestCenter) / 2.;
+
+ // now randomly choose if we take the closest bin or the one next to the closest one.
+ // with this, we kind of linearly interpolate between two bins if an input value
+ // betweeen these two bins does occure many times.
+ double rand = CLHEP::RandFlat::shoot(m_randomEngine);
+ chosenBin[inputPar] = ( rand >= distance ) ? bin : bin+binDelta;
+
+ // check if we have chosen an over- or underflow-bin
+ // if this has happened, choose the last bin still in range
+
+ if ( chosenBin[inputPar] <= 0 ) chosenBin[inputPar] = 1;
+ else if ( chosenBin[inputPar] > curaxis->GetNbins() ) chosenBin[inputPar] = curaxis->GetNbins();
+ }
+
+ // now get the bin number (since all histograms are binned
+ // in exactly the same way it should be the same for all
+ // of them)
+ Int_t bin = 0;
+ if (numInputPar == 1) bin = m_par[0]->GetBin( chosenBin[0] );
+ else if (numInputPar ==2 ) bin = m_par[0]->GetBin( chosenBin[0], chosenBin[1] );
+ else if (numInputPar == 3) bin = m_par[0]->GetBin( chosenBin[0], chosenBin[1], chosenBin[2] );
+ // TODO: implement case of >3 input parameters
+
+ // free some memory
+ delete [] axis;
+
+ // fill in the parameters into the distribution-function
+ // (e.g. for the current energy & eta regime)
+ for (int funcPar = 0; funcPar < m_pdf->GetNpar(); funcPar++) {
+ double value = m_par[funcPar]->GetBinContent(bin);
+ m_pdf->SetParameter(funcPar, value);
+ }
+
+ // now finally we can determine the random number from the given distribution
+ //
+ // set the minimum and maximum random values
+ double xMin = ( (randMin != 0.) && (randMin < m_randmax->GetBinContent(bin)) )
+ ? randMin : m_randmin->GetBinContent(bin);
+ double xMax = ( (randMax != 0.) && (randMax <= m_randmax->GetBinContent(bin)) )
+ ? randMax : m_randmax->GetBinContent(bin);
+
+ // free some more memory
+ delete [] chosenBin;
+
+ // (1.) for discrete distributions
+ // - this is implemented separately from the continuous case,
+ // since it significantly speeds up the creation of random
+ // numbers from a distribution with a high negative gradient
+ // - for this the sum( f(x_i) ) has to be normalized !
+ if (discrete) {
+ // get a flat random value between 0. and 1.
+ double flat = CLHEP::RandFlat::shoot(m_randomEngine);
+
+ double sum = 0.;
+
+ // set the minimum x value
+ xMin = lround(xMin);
+
+ // now determine the corresponding value from the distribution
+ for (int rand=xMin; rand<lround(xMax); rand++) {
+ // sum up the probabilities for each current bin
+ sum += m_pdf->Eval(rand);
+ // if the sum gets larger than the random value
+ // -> the current 'rand' is the value from the distribution
+ // -> return this value
+ if (sum >= flat) return ((double)rand);
+ }
+
+ // This point is reached if the the x value corresponding to the random
+ // value is the maximum possible value in the probability distribution
+ // or if the distribution is not normalized.
+ // -> return the maximum value
+ return (double)m_pdf->GetXmax();
+ }
+
+
+ // (2.) for continuous distributions
+ // * algorithmus: rejection sampling
+ else {
+ // find max and min values in the PDF
+ double yMin = m_pdf->GetMinimum(xMin, xMax);
+ //if (yMin < 0.) yMin = 0.; // in case the PDF contains negative values
+ double yMax = m_pdf->GetMaximum(xMin, xMax);
+
+ // now get a random value rx with the right distribution
+ // - use the count variable just for backup to prevent
+ // an endless loop
+ for (int count=0; count < 1e7; count++) {
+ // first choose a random value within [yMin, yMax]
+ double ry = CLHEP::RandFlat::shoot(m_randomEngine)*(yMax-yMin) + yMin;
+ // then choose a random value within [xMin, xMax]
+ double rx = CLHEP::RandFlat::shoot(m_randomEngine)*(xMax-xMin) + xMin;
+ // calculate y=pdf(rx)
+ double y = m_pdf->Eval(rx);
+ // if ry<=pdf(rx) is fulfilled, we have our random value
+ if ( ry <= y ) return rx;
+ }
+ }
+
+ // this point will only be reached, if something went wrong
+ return 0.;
+}
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PunchThroughParticle.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PunchThroughParticle.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..e521f618a701c527e9df49ec6560ab72a8696a83
--- /dev/null
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PunchThroughParticle.cxx
@@ -0,0 +1,182 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// PunchThroughParticle.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+// class header
+#include "ISF_PunchThroughTools/PunchThroughParticle.h"
+
+// ROOT
+#include "TH2F.h"
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+ISF::PunchThroughParticle::PunchThroughParticle(int pdg, bool doAnti):
+ m_pdgId(pdg),
+ m_doAnti(doAnti),
+ m_minEnergy(0.),
+ m_maxNum(-1),
+ m_numParticlesFactor(1.),
+ m_energyFactor(1.),
+ m_posAngleFactor(1.),
+ m_momAngleFactor(1.),
+ m_corrPdg(0),
+ m_corrMinEnergy(0),
+ m_corrFullEnergy(0),
+ m_histCorrLowE(0),
+ m_histCorrHighE(0),
+ m_corrHistDomains(0),
+ m_pdfNumParticles(0), //does this number ever change?
+ m_pdfExitEnergy(0),
+ m_pdfExitDeltaTheta(0),
+ m_pdfExitDeltaPhi(0),
+ m_pdfMomDeltaTheta(0),
+ m_pdfMomDeltaPhi(0)
+{ }
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+void ISF::PunchThroughParticle::setMinEnergy(double minEnergy)
+{
+ m_minEnergy = minEnergy;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+void ISF::PunchThroughParticle::setMaxNumParticles(int maxNum)
+{
+ m_maxNum = maxNum;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+void ISF::PunchThroughParticle::setNumParticlesFactor(double numFactor)
+{
+ m_numParticlesFactor = numFactor;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+void ISF::PunchThroughParticle::setEnergyFactor(double energyFactor)
+{
+ m_energyFactor = energyFactor;
+}
+
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+void ISF::PunchThroughParticle::setPosAngleFactor(double posAngleFactor)
+{
+ m_posAngleFactor = posAngleFactor;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+void ISF::PunchThroughParticle::setMomAngleFactor(double momAngleFactor)
+{
+ m_momAngleFactor = momAngleFactor;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+void ISF::PunchThroughParticle::setNumParticlesPDF(PDFcreator *pdf)
+{
+ m_pdfNumParticles = pdf;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+void ISF::PunchThroughParticle::setCorrelation(int corrPdg,
+ TH2F *histLowE, TH2F *histHighE,
+ double minCorrelationEnergy,
+ double fullCorrelationEnergy,
+ double lowE,
+ double midE,
+ double upperE)
+{
+ m_corrPdg = corrPdg;
+ m_histCorrLowE = histLowE;
+ m_histCorrHighE = histHighE;
+ m_corrMinEnergy = minCorrelationEnergy;
+ m_corrFullEnergy = fullCorrelationEnergy;
+
+ m_corrHistDomains = new double [3];
+ m_corrHistDomains[0] = lowE;
+ m_corrHistDomains[1] = midE;
+ m_corrHistDomains[2] = upperE;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+void ISF::PunchThroughParticle::setExitEnergyPDF(PDFcreator *pdf)
+{
+ m_pdfExitEnergy = pdf;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+void ISF::PunchThroughParticle::setExitDeltaThetaPDF(PDFcreator *pdf)
+{
+ m_pdfExitDeltaTheta = pdf;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+void ISF::PunchThroughParticle::setExitDeltaPhiPDF(PDFcreator *pdf)
+{
+ m_pdfExitDeltaPhi = pdf;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+void ISF::PunchThroughParticle::setMomDeltaThetaPDF(PDFcreator *pdf)
+{
+ m_pdfMomDeltaTheta = pdf;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+void ISF::PunchThroughParticle::setMomDeltaPhiPDF(PDFcreator *pdf)
+{
+ m_pdfMomDeltaPhi = pdf;
+}
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PunchThroughTool.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PunchThroughTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ca6eadb52854c0aa0585b51bc80db46c85700021
--- /dev/null
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/PunchThroughTool.cxx
@@ -0,0 +1,1060 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// PunchThroughTool.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+// class header
+#include "ISF_PunchThroughTools/PunchThroughTool.h"
+
+// standard C++ libraries
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <algorithm>
+#include <vector>
+
+// standard C libraries
+#include <math.h>
+
+// Control
+#include "DataModel/DataVector.h"
+
+// Gaudi & StoreGate
+#include "GaudiKernel/IPartPropSvc.h"
+
+// HepMC
+#include "HepMC/SimpleVector.h"
+#include "HepMC/GenVertex.h"
+#include "HepMC/GenParticle.h"
+#include "HepMC/GenEvent.h"
+#include "HepPDT/ParticleDataTable.hh"
+
+// CLHEP
+#include "CLHEP/Random/RandFlat.h"
+
+// ROOT
+#include "TFile.h"
+#include "TH2F.h"
+#include "TAxis.h"
+#include "TH1.h"
+#include "TRandom.h"
+#include "TMath.h"
+
+// PathResolver
+#include "PathResolver/PathResolver.h"
+
+//ISF
+#include "ISF_Event/ISFParticle.h"
+#include "ISF_PunchThroughTools/PDFcreator.h"
+#include "ISF_PunchThroughTools/PunchThroughParticle.h"
+#include "ISF_Interfaces/IGeoIDSvc.h"
+
+//Barcode
+#include "BarcodeInterfaces/IBarcodeSvc.h"
+
+//Geometry
+#include "SubDetectorEnvelopes/IEnvelopeDefSvc.h"
+
+//Amg
+#include "GeoPrimitives/GeoPrimitivesHelpers.h"
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+ISF::PunchThroughTool::PunchThroughTool( const std::string& type,
+ const std::string& name,
+ const IInterface* parent )
+ : AthAlgTool(type, name, parent),
+ m_parentGenEvt(0),
+ m_particleDataTable(0),
+ m_randomSvc("AtDSFMTGenSvc", name),
+ m_randomEngineName("ISFRnd"),
+ m_randomEngine(0),
+ m_filenameLookupTable("CaloPunchThroughParametrisation.root"),
+ m_pdgInitiators(),
+ m_punchThroughParticles(),
+ m_doAntiParticles(),
+ m_correlatedParticle(),
+ m_minCorrEnergy(),
+ m_fullCorrEnergy(),
+ m_posAngleFactor(),
+ m_momAngleFactor(),
+ m_minEnergy(),
+ m_maxNumParticles(),
+ m_numParticlesFactor(),
+ m_energyFactor(),
+ m_particlePropSvc("PartPropSvc",name),
+ m_geoIDSvc("iGeant4::G4PolyconeGeoIDSvc/ISF_G4PolyconeGeoIDSvc",name),
+ m_barcodeSvc("BarcodeSvc",name),
+ m_envDefSvc("AtlasGeometry_EnvelopeDefSvc",name),
+ m_beamPipe(500.)
+{
+ // declare this as an interface
+ declareInterface<IPunchThroughTool>(this);
+
+ // property definition and initialization - allows to set variables from python
+
+ declareProperty( "FilenameLookupTable", m_filenameLookupTable );
+ declareProperty( "PunchThroughInitiators", m_pdgInitiators );
+ declareProperty( "PunchThroughParticles", m_punchThroughParticles );
+ declareProperty( "DoAntiParticles", m_doAntiParticles );
+ declareProperty( "CorrelatedParticle", m_correlatedParticle );
+ declareProperty( "MinCorrelationEnergy", m_minCorrEnergy );
+ declareProperty( "FullCorrelationEnergy", m_fullCorrEnergy );
+ declareProperty( "ScalePosDeflectionAngles", m_posAngleFactor );
+ declareProperty( "ScaleMomDeflectionAngles", m_momAngleFactor );
+ declareProperty( "MinEnergy", m_minEnergy );
+ declareProperty( "MaxNumParticles", m_maxNumParticles );
+ declareProperty( "NumParticlesFactor", m_numParticlesFactor );
+ declareProperty( "EnergyFactor", m_energyFactor );
+ declareProperty( "RandomNumberService", m_randomSvc, "Random number generator" );
+ declareProperty( "RandomStreamName", m_randomEngineName, "Name of the random number stream");
+ declareProperty( "BarcodeSvc", m_barcodeSvc );
+ declareProperty( "EnvelopeDefSvc", m_envDefSvc );
+ declareProperty( "BeamPipeRadius", m_beamPipe );
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+ISF::PunchThroughTool::~PunchThroughTool()
+{}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+StatusCode ISF::PunchThroughTool::initialize()
+{
+ ATH_MSG_INFO( "initialize()" );
+
+ // resolving lookuptable file
+ std::string resolvedFileName = PathResolver::find_file (m_filenameLookupTable, "DATAPATH");
+ if (resolvedFileName != "")
+ {
+ ATH_MSG_INFO( "[ punchthrough ] Parameterisation file found: " << resolvedFileName );
+ }
+ else
+ {
+ ATH_MSG_ERROR( "[ punchthrough ] Parameterisation file not found" );
+ return StatusCode::FAILURE;
+ }
+
+ // open the LookupTable file
+ m_fileLookupTable = new TFile( resolvedFileName.c_str(), "READ");
+ if (!m_fileLookupTable)
+ {
+ ATH_MSG_WARNING("[ punchthrough ] unable to open the lookup-tabel for the punch-through simulation (file does not exist)");
+ return StatusCode::FAILURE;
+ }
+
+ if (!m_fileLookupTable->IsOpen())
+ {
+ ATH_MSG_WARNING("[ punchthrough ] unable to open the lookup-tabel for the punch-through simulation (wrong or empty file?)");
+ return StatusCode::FAILURE;
+ }
+
+ // retrieve the ParticleProperties handle
+ if ( m_particlePropSvc.retrieve().isFailure() )
+ {
+ ATH_MSG_FATAL( "[ punchthrough ] Can not retrieve " << m_particlePropSvc << " ! Abort. " );
+ return StatusCode::FAILURE;
+ }
+
+ // and the particle data table
+ m_particleDataTable = m_particlePropSvc->PDT();
+ if (m_particleDataTable==0)
+ {
+ ATH_MSG_FATAL( " [ punchthrough ] Could not get ParticleDataTable! Cannot associate pdg code with charge! Abort. " );
+ return StatusCode::FAILURE;
+ }
+
+ // Random number service
+ if ( m_randomSvc.retrieve().isFailure() )
+ {
+ ATH_MSG_ERROR( "[ punchthrough ] Could not retrieve " << m_randomSvc );
+ return StatusCode::FAILURE;
+ }
+
+ // Get own engine with own seeds:
+ m_randomEngine = m_randomSvc->GetEngine(m_randomEngineName);
+ if (!m_randomEngine)
+ {
+ ATH_MSG_ERROR( "[ punchthrough ] Could not get random engine '" << m_randomEngineName << "'" );
+ return StatusCode::FAILURE;
+ }
+
+ // Geometry identifier service
+ if ( !m_geoIDSvc.empty() && m_geoIDSvc.retrieve().isFailure())
+ {
+ ATH_MSG_FATAL ( "[ punchthrough ] Could not retrieve GeometryIdentifier Service. Abort");
+ return StatusCode::FAILURE;
+ }
+
+ //barcode service
+
+ if (m_barcodeSvc.retrieve().isFailure() )
+ {
+ ATH_MSG_ERROR( "[ punchthrough ] Could not retrieve " << m_barcodeSvc );
+ return StatusCode::FAILURE;
+ }
+
+ //envelope definition service
+ if (m_envDefSvc.retrieve().isFailure() )
+ {
+ ATH_MSG_ERROR( "[ punchthrough ] Could not retrieve " << m_envDefSvc );
+ return StatusCode::FAILURE;
+ }
+
+//--------------------------------------------------------------------------------
+ // register all the punch-through particles which will be simulated
+ for ( unsigned int num = 0; num < m_punchThroughParticles.size(); num++ )
+ {
+ int pdg = m_punchThroughParticles[num];
+ // if no information is given on the creation of anti-particles -> do not simulate anti-particles
+ bool doAnti = ( num < m_doAntiParticles.size() ) ? m_doAntiParticles[num] : false;
+ // if no information is given on the minimum energy -> take 50. MeV as default
+ double minEnergy = ( num < m_minEnergy.size() ) ? m_minEnergy[num] : 50.;
+ // if no information is given on the maximum number of punch-through particles -> take -1 as default
+ int maxNum = ( num < m_minEnergy.size() ) ? m_maxNumParticles[num] : -1;
+ // if no information is given on the scale factor for the number of particles -> take 1. as defaulft
+ double numFactor = ( num < m_numParticlesFactor.size() ) ? m_numParticlesFactor[num] : 1.;
+ // if no information is given on the position angle factor -> take 1.
+ double posAngleFactor = ( num < m_posAngleFactor.size() ) ? m_posAngleFactor[num] : 1.;
+ // if no information is given on the momentum angle factor -> take 1.
+ double momAngleFactor = ( num < m_momAngleFactor.size() ) ? m_momAngleFactor[num] : 1.;
+ // if no information is given on the scale factor for the energy -> take 1. as default
+ double energyFactor = ( num < m_energyFactor.size() ) ? m_energyFactor[num] : 1.;
+
+ // register the particle
+ ATH_MSG_VERBOSE("[ punchthrough ] registering punch-through particle type with pdg = " << pdg );
+ if ( registerParticle( pdg, doAnti, minEnergy, maxNum, numFactor, energyFactor, posAngleFactor, momAngleFactor )
+ != StatusCode::SUCCESS)
+ {
+ ATH_MSG_ERROR("[ punchthrough ] unable to register punch-through particle type with pdg = " << pdg);
+ }
+ }
+
+ // TODO: implement punch-through parameters for different m_pdgInitiators
+ // currently m_pdgInitiators is only used to filter out particles
+
+ // check if more correlations were given than particle types were registered
+ unsigned int numCorrelations = m_correlatedParticle.size();
+ if ( numCorrelations > m_punchThroughParticles.size() )
+ {
+ ATH_MSG_WARNING("[ punchthrough ] more punch-through particle correlations are given, than punch-through particle types are registered (skipping the last ones)");
+ numCorrelations = m_punchThroughParticles.size();
+ }
+
+ // now register correlation between particles
+ for ( unsigned int num = 0; num < numCorrelations; num++ )
+ {
+ int pdg1 = m_punchThroughParticles[num];
+ int pdg2 = m_correlatedParticle[num];
+ double fullCorrEnergy = ( num < m_fullCorrEnergy.size() ) ? m_fullCorrEnergy[num] : 0.;
+ double minCorrEnergy = ( num < m_minCorrEnergy.size() ) ? m_minCorrEnergy[num] : 0.;
+
+ // if correlatedParticle==0 is given -> no correlation
+ if ( ! pdg2) continue;
+ // register it
+ if ( registerCorrelation(pdg1, pdg2, minCorrEnergy, fullCorrEnergy) != StatusCode::SUCCESS )
+ {
+ ATH_MSG_ERROR("[ punchthrough ] unable to register punch-through particle correlation for pdg1=" << pdg1 << " pdg2=" << pdg2 );
+ }
+ }
+
+ // get the calo-MS border coordinates. Look at calo and MS geometry definitions, if same R and Z -> boundary surface
+
+ RZPairVector* rzMS = &(m_envDefSvc->getMuonRZValues());
+ RZPairVector* rzCalo = &(m_envDefSvc->getCaloRZValues());
+
+ bool found1, found2;
+ found1=false; found2=false;
+
+ for ( size_t i=0; i<rzCalo->size();++i)
+ {
+ double r_tempCalo = rzCalo->at(i).first;
+ double z_tempCalo = rzCalo->at(i).second;
+
+ if (r_tempCalo> m_beamPipe)
+ {
+ for ( size_t j=0; j<rzMS->size();++j)
+ {
+ double r_tempMS =rzMS->at(j).first;
+ double z_tempMS =rzMS->at(j).second;
+
+ if (r_tempCalo==r_tempMS && z_tempCalo==z_tempMS && found1==false )
+ {
+ R1=r_tempMS;
+ z1=fabs(z_tempMS);
+ found1=true;
+ continue;
+ }
+ else if (r_tempCalo==r_tempMS && z_tempCalo==z_tempMS && r_tempCalo!=R1 && fabs(z_tempCalo)!=z1)
+ {
+ R2=r_tempMS;
+ z2=fabs(z_tempMS);
+ found2=true;
+ }
+ }
+
+ if (found1==true && found2==true) break;
+ }
+ }
+
+ //in case geometry description changes
+ if (found1 == false) ATH_MSG_ERROR ("first coordinate of calo-MS border not found");
+ if (found2 == false) ATH_MSG_ERROR ("second coordinate of calo-MS border not found; first one is: R1 ="<<R1<<" z1 ="<<z1);
+
+ //now order the found values
+ double r_temp, z_temp;
+ if (R1>R2) { r_temp=R1; R1=R2; R2=r_temp; } //R1 - smaller one
+ if (z1<z2) { z_temp=z1; z1=z2; z2=z_temp; } //z1 - bigger one
+
+ if (R1==R2 || z1==z2) ATH_MSG_ERROR ("[punch-though] Bug in propagation calculation! R1="<<R1<<" R2 = "<<R2<<" z1="<<z1<<" z2= "<<z2 );
+ else ATH_MSG_DEBUG ("calo-MS boundary coordinates: R1="<<R1<<" R2 = "<<R2<<" z1="<<z1<<" z2= "<<z2);
+
+ ATH_MSG_INFO( "punchthrough initialization is successful" );
+ return StatusCode::SUCCESS;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+StatusCode ISF::PunchThroughTool::finalize()
+{
+ ATH_MSG_INFO( "[punchthrough] finalize() successful" );
+
+ // close the file with the lookuptable
+ m_fileLookupTable->Close();
+
+ return StatusCode::SUCCESS;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+const ISF::ISFParticleContainer* ISF::PunchThroughTool::computePunchThroughParticles(const ISF::ISFParticle &isfp) const
+{
+ ATH_MSG_DEBUG( "[ punchthrough ] starting punch-through simulation");
+
+ // reset the output particle collection
+ m_isfpCont = new ISF::ISFParticleContainer();
+
+ // reset the parent GenEvent
+ m_parentGenEvt = 0;
+
+ // store the initial particle state locally
+ m_initPs = &isfp;
+
+ ATH_MSG_VERBOSE("[ GeoIDSvc ] position of the input particle: r"<<m_initPs->position().perp()<<" z= "<<m_initPs->position().z() );
+ //if not on ID surface - don't simulate
+
+ if ( m_geoIDSvc->inside(m_initPs->position(),AtlasDetDescr::fAtlasID) != 1)
+ {
+ ATH_MSG_DEBUG("[ GeoIDSvc ] input particle position is not on reference surface -> no punch-through simulation");
+ return 0;
+ }
+
+ //check if it points to the calorimeter - if not, don't simulate
+
+ if ( m_geoIDSvc->identifyNextGeoID(*m_initPs) != AtlasDetDescr::fAtlasCalo)
+ {
+ ATH_MSG_VERBOSE ("[ GeoIDSvc ] input particle doesn't point to calorimeter"<< "Next GeoID: "<<m_geoIDSvc->identifyNextGeoID(*m_initPs) );
+ return 0;
+ }
+
+ // check if the particle's pdg is registered as a punch-through-causing type
+ {
+ std::vector<int>::const_iterator pdgIt = m_pdgInitiators.begin();
+ std::vector<int>::const_iterator pdgItEnd = m_pdgInitiators.end();
+ // loop over all known punch-through initiators
+ for ( ; pdgIt != pdgItEnd; ++pdgIt)
+ { if (abs(m_initPs->pdgCode()) == *pdgIt) break; }
+
+ // particle will not cause punch-through -> bail out
+ if (pdgIt == pdgItEnd)
+ {
+ ATH_MSG_DEBUG("[ punchthrough ] particle is not registered as punch-through initiator. Dropping it in the calo.");
+ return 0;
+ }
+ }
+
+ //if initial particle is on ID surface, points to the calorimeter and is a punch-through initiator
+
+ // store initial particle mass
+ double mass = m_initPs->mass();
+ //store its barcode
+ m_primBC = m_initPs->barcode();
+ // now store some parameters from the given particle locally
+ // in this class
+ // -> the energy
+ m_initEnergy = sqrt( m_initPs->momentum().mag2() + mass*mass );
+ // -> get geometrical properties
+ // -> the pseudorapidity eta
+ m_initEta = m_initPs->position().eta();
+ // -> the angle theta
+ m_initTheta = m_initPs->position().theta();
+ // -> the azimutal angle phi
+ m_initPhi = m_initPs->position().phi();
+
+ // this is the place where the magic is done:
+ // test for each registered punch-through pdg if a punch-through
+ // occures and create these particles
+ // -> therefore loop over all registered pdg ids
+ std::map<int,PunchThroughParticle*>::const_iterator it;
+ // to keep track of the correlated particles which were already simulated:
+ // first int is pdg, second int is number of particles created
+ std::map<int, int> corrPdgNumDone;
+
+ // loop over all particle pdgs
+ for ( it = m_particles.begin(); it != m_particles.end(); ++it )
+ {
+ // the pdg that is currently treated
+ int doPdg = it->first;
+ // get the current particle's correlated pdg
+ int corrPdg = it->second->getCorrelatedPdg();
+
+ // if there is a correlated particle type to this one
+ if (corrPdg)
+ {
+ // find out if the current pdg was already simulated
+ std::map<int,int>::iterator pos = corrPdgNumDone.find(doPdg);
+ // if the current pdg was not simulated yet, find out if
+ // it's correlated one was simulated
+ if ( pos == corrPdgNumDone.end() ) pos = corrPdgNumDone.find(corrPdg);
+
+ // neither this nor the correlated particle type was simulated
+ // so far:
+ if ( pos == corrPdgNumDone.end() )
+ {
+ // -> roll a dice if we create this particle or its correlated one
+ if ( CLHEP::RandFlat::shoot(m_randomEngine) > 0.5 ) doPdg = corrPdg;
+ // now create the particles with the given pdg and note how many
+ // particles of this pdg are created
+ corrPdgNumDone[doPdg] = getAllParticles(doPdg);
+ }
+
+ // one of the two correlated particle types was already simulated
+ // 'pos' points to the already simulated one
+ else
+ {
+ // get the pdg of the already simulated particle and the number
+ // of these particles that were created
+ int donePdg = pos->first;
+ int doneNumPart = pos->second;
+ // set the pdg of the particle type that will be done
+ if (donePdg == doPdg) doPdg = corrPdg;
+
+ // now create the correlated particles
+ getCorrelatedParticles(doPdg, doneNumPart);
+ // note: no need to take note, that this particle type is now simulated,
+ // since this is the second of two correlated particles, which is
+ // simulated and we do not have correlations of more than two particles.
+ }
+
+ // if no correlation for this particle
+ // -> directly create all particles with the current pdg
+ }
+ else getAllParticles(doPdg);
+
+ } // for-loop over all particle pdgs
+
+ if (m_isfpCont->size() > 0) ATH_MSG_DEBUG( "[ punchthrough ] returning ISFparticle vector , size: "<<m_isfpCont->size() );
+
+ return m_isfpCont;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+int ISF::PunchThroughTool::getAllParticles(int pdg, int numParticles) const
+{
+ // first check if the ISF particle vector already exists
+ if (!m_isfpCont) m_isfpCont = new ISFParticleContainer();
+
+ // get the current particle
+ PunchThroughParticle *p = m_particles[pdg];
+
+ // if no number of particles (=-1) was handed over as an argument
+ // -> get the number of particles from the pdf
+ if ( numParticles < 0 )
+ {
+ // prepare the function arguments for the PDFcreator class
+ std::vector<double> parameters;
+ parameters.push_back( m_initEnergy );
+ parameters.push_back( fabs(m_initEta) );
+ // the maximum number of particles which should be produced
+ // if no maximum number is given, this is -1
+ int maxParticles = p->getMaxNumParticles();
+
+ // get the right number of punch-through particles
+ // and ensure that we do not create too many particles
+ do
+ {
+ numParticles = lround( p->getNumParticlesPDF()->getRand(parameters, true) );
+ // scale the number of particles if requested
+ numParticles = lround( numParticles *= p->getNumParticlesFactor() );
+ }
+ while ( (maxParticles >= 0.) && (numParticles > maxParticles) );
+ }
+
+ ATH_MSG_VERBOSE("[ punchthrough ] adding " << numParticles << " punch-through particles with pdg " << pdg);
+
+ // now create the exact number of particles which was just computed before
+ double energyRest = m_initEnergy;
+ double minEnergy = p->getMinEnergy();
+ int numCreated = 0;
+
+ for ( numCreated = 0; (numCreated < numParticles) && (energyRest > minEnergy); numCreated++ )
+ {
+ // create one particle which fullfills the right energy distribution
+ ISF::ISFParticle *par = getOneParticle(pdg, energyRest);
+
+ // if something went wrong
+ if (!par)
+ {
+ ATH_MSG_ERROR("[ punchthrough ] something went wrong while creating punch-through particles");
+ return StatusCode::FAILURE;
+ }
+
+ // get the energy of the particle which was just created
+ double restMass = m_particleDataTable->particle(abs(pdg))->mass();
+ double curEnergy = sqrt(par->momentum().mag2() + restMass*restMass);
+
+ // calculate the maximum energy to be available for all
+ // following punch-through particles created
+ energyRest -= curEnergy;
+
+ // add this ISFparticle to the vector
+ m_isfpCont->push_back( par );
+ }
+
+ // the number of particles which was created is numCreated
+ return (numCreated);
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+int ISF::PunchThroughTool::getCorrelatedParticles(int pdg, int corrParticles) const
+{
+ // get the PunchThroughParticle class
+ PunchThroughParticle *p = m_particles[pdg];
+
+ // (1.) decide if we do correlation or not
+ double rand = CLHEP::RandFlat::shoot(m_randomEngine)
+ *(p->getFullCorrelationEnergy()-p->getMinCorrelationEnergy())
+ + p->getMinCorrelationEnergy();
+ if ( m_initEnergy < rand )
+ {
+ // here we do not do correlation
+ return getAllParticles(pdg);
+ }
+
+ // (2.) if this point is reached, we do correlation
+ // decide which 2d correlation histogram to use
+ double *histDomains = p->getCorrelationHistDomains();
+ TH2F *hist2d = 0;
+ // compute the center values of the lowE and highE
+ // correlation histogram domains
+ if ( m_initEnergy < histDomains[1])
+ {
+ // initial energy lower than border between lowEnergy and highEnergy histogram domain
+ // --> choose lowEnergy correlation histogram
+ hist2d = p->getCorrelationLowEHist();
+ }
+ else
+ {
+ double rand = CLHEP::RandFlat::shoot(m_randomEngine)*(histDomains[2]-histDomains[1])
+ + histDomains[1];
+ hist2d = ( m_initEnergy < rand) ? p->getCorrelationLowEHist()
+ : p->getCorrelationHighEHist();
+ }
+
+ // get the correlation 2d histogram
+
+ // now find out where on the x-axis the the bin for number of
+ // correlated particles is
+ Int_t xbin = hist2d->GetXaxis()->FindFixBin(corrParticles);
+ int numParticles = 0;
+ int maxParticles = p->getMaxNumParticles();
+ // now the distribution along the y-axis is a PDF for the number
+ // of 'pdg' particles
+ do
+ {
+ double rand = CLHEP::RandFlat::shoot(m_randomEngine);
+ double sum = 0.;
+ for ( int ybin = 1; ybin <= hist2d->GetNbinsY(); ybin++ )
+ {
+ sum += hist2d->GetBinContent(xbin, ybin);
+ // check if we choose the current bin or not
+ if ( sum >= rand )
+ {
+ numParticles = ybin - 1;
+ break;
+ }
+ }
+ // scale the number of particles is requested
+ numParticles = lround( numParticles * p->getNumParticlesFactor() );
+ }
+ while ( (maxParticles >= 0.) && (numParticles > maxParticles) );
+
+ // finally create this exact number of particles
+ return getAllParticles(pdg, numParticles);
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+ISF::ISFParticle *ISF::PunchThroughTool::getOneParticle(int pdg, double maxEnergy) const
+{
+ // get a local copy of the needed punch-through particle class
+ PunchThroughParticle *p = m_particles[pdg];
+
+ // (1.) decide if we create a particle or an anti-particle
+ int anti = 1;
+ if ( p->getdoAnti() )
+ {
+ // get a random-value
+ double rand = CLHEP::RandFlat::shoot(m_randomEngine);
+ // 50/50 chance to be a particle or its anti-particle
+ if (rand > 0.5) anti = -1;
+ }
+
+ // (2.) get the right punch-through distributions
+ // prepare the function arguments for the PDFcreator class
+ std::vector<double> parInitEnergyEta;
+ parInitEnergyEta.push_back( m_initEnergy );
+ parInitEnergyEta.push_back( fabs(m_initEta) );
+
+ // (2.1) get the energy
+ double energy = p->getExitEnergyPDF()->getRand(
+ parInitEnergyEta, false, p->getMinEnergy(), maxEnergy/p->getEnergyFactor() );
+ energy *= p->getEnergyFactor(); // scale the energy if requested
+
+ // set up the new parameters for the next PDFcreator::getRand calls
+ std::vector<double> parExitEnergyInitEta;
+ parExitEnergyInitEta.push_back( energy );
+ parExitEnergyInitEta.push_back( fabs(m_initEta) );
+
+ // (2.2) get the particles delta theta relative to the incoming particle
+ double theta = 0;
+ // loop to keep theta within range [0,PI]
+ do
+ {
+ // get random value
+ double deltaTheta = p->getExitDeltaThetaPDF()->getRand(
+ parExitEnergyInitEta, false );
+ // decide if delta positive/negative
+ deltaTheta *= ( CLHEP::RandFlat::shoot(m_randomEngine) > 0.5 ) ? 1. : -1.;
+ // calculate the exact theta value of the later created
+ // punch-through particle
+ theta = m_initTheta + deltaTheta*p->getPosAngleFactor();
+ }
+ while ( (theta > M_PI) || (theta < 0.) );
+
+ // (2.3) get the particle's delta phi relative to the incoming particle
+
+ double deltaPhi = p->getExitDeltaPhiPDF()->getRand(
+ parExitEnergyInitEta, false );
+ deltaPhi *= ( CLHEP::RandFlat::shoot(m_randomEngine) > 0.5 ) ? 1. : -1.;
+
+ // keep phi within range [-PI,PI]
+ double phi = m_initPhi + deltaPhi*p->getPosAngleFactor();
+ while ( fabs(phi) > 2*M_PI) phi /= 2.;
+ while (phi > M_PI) phi -= 2*M_PI;
+ while (phi < -M_PI) phi += 2*M_PI;
+
+ // set up the new parameters for the next PDFcreator::getRand calls
+ std::vector<double> parExitEnergyExitEta;
+ parExitEnergyExitEta.push_back( energy );
+ parExitEnergyExitEta.push_back( -log(tan(theta/2.)) );
+
+ // (2.4) get the particle momentum delta theta, relative to its position
+ //
+ // loop to keep momTheta within range [0,PI]
+
+ double momTheta = 0.;
+ do
+ {
+ // get random value
+ double momDeltaTheta = p->getMomDeltaThetaPDF()->getRand(
+ parExitEnergyExitEta, false );
+ // decide if delta positive/negative
+ momDeltaTheta *= ( CLHEP::RandFlat::shoot(m_randomEngine) > 0.5 ) ? 1. : -1.;
+ // calculate the exact momentum theta value of the later created
+ // punch-through particle
+ momTheta = theta + momDeltaTheta*p->getMomAngleFactor();
+ }
+ while ( (momTheta > M_PI) || (momTheta < 0.) );
+
+ // (2.5) get the particle momentum delta phi, relative to its position
+
+ double momDeltaPhi = p->getMomDeltaPhiPDF()->getRand(
+ parExitEnergyExitEta, false );
+ momDeltaPhi *= ( CLHEP::RandFlat::shoot(m_randomEngine) > 0.5 ) ? 1. : -1.;
+
+ double momPhi = phi + momDeltaPhi*p->getMomAngleFactor();
+ // keep momPhi within range [-PI,PI]
+ while ( fabs(momPhi) > 2*M_PI) momPhi /= 2.;
+ while (momPhi > M_PI) momPhi -= 2*M_PI;
+ while (momPhi < -M_PI) momPhi += 2*M_PI;
+
+ //assign barcodes to the produced particles
+ m_processCode = 0;
+ m_secBC = m_barcodeSvc->newSecondary( m_primBC, m_processCode);
+
+ // (**) finally create the punch-through particle as a ISFParticle
+
+ ISF::ISFParticle *par = createExitPs( pdg*anti, energy, theta, phi, momTheta, momPhi);
+
+ return par;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+StatusCode
+ISF::PunchThroughTool::registerParticle(int pdg, bool doAntiparticle,
+ double minEnergy, int maxNumParticles, double numParticlesFactor,
+ double energyFactor, double posAngleFactor, double momAngleFactor)
+{
+ // read in the data needed to construct the distributions for the number of punch-through particles
+
+ // (1.) get the distribution function for the number of punch-through particles
+ PDFcreator *pdf_num = readLookuptablePDF(pdg, "NumExitPDG");
+ if (!pdf_num ) return StatusCode::FAILURE; // return error if something went wrong
+
+ // (2.) get the PDF for the punch-through energy
+ PDFcreator *pdf_energy = readLookuptablePDF(pdg, "ExitEnergyPDG");
+ if (!pdf_energy)
+ {
+ delete pdf_num;
+ return StatusCode::FAILURE; // return error if something went wrong
+ }
+
+ // (3.) get the PDF for the punch-through particles difference in
+ // theta compared to the incoming particle
+ PDFcreator *pdf_theta = readLookuptablePDF(pdg, "ExitDeltaThetaPDG");
+ if (!pdf_theta)
+ {
+ delete pdf_num; delete pdf_energy;
+ return StatusCode::FAILURE;
+ }
+
+ // (4.) get the PDF for the punch-through particles difference in
+ // phi compared to the incoming particle
+ PDFcreator *pdf_phi = readLookuptablePDF(pdg, "ExitDeltaPhiPDG");
+ if (!pdf_phi)
+ {
+ delete pdf_num; delete pdf_energy; delete pdf_theta;
+ return StatusCode::FAILURE;
+ }
+
+ // (5.) get the PDF for the punch-through particle momentum delta theta angle
+ PDFcreator *pdf_momTheta = readLookuptablePDF(pdg, "MomDeltaThetaPDG");
+ if (!pdf_momTheta)
+ {
+ delete pdf_num; delete pdf_energy; delete pdf_theta; delete pdf_phi;
+ return StatusCode::FAILURE;
+ }
+
+ // (6.) get the PDF for the punch-through particle momentum delta phi angle
+ PDFcreator *pdf_momPhi = readLookuptablePDF(pdg, "MomDeltaPhiPDG");
+ if (!pdf_momPhi)
+ {
+ delete pdf_num; delete pdf_energy; delete pdf_theta; delete pdf_phi; delete pdf_momTheta;
+ return StatusCode::FAILURE;
+ }
+
+ // (7.) now finally store all this in the right std::map
+ PunchThroughParticle *particle = new PunchThroughParticle(pdg, doAntiparticle);
+ particle->setNumParticlesPDF(pdf_num);
+ particle->setExitEnergyPDF(pdf_energy);
+ particle->setExitDeltaThetaPDF(pdf_theta);
+ particle->setExitDeltaPhiPDF(pdf_phi);
+ particle->setMomDeltaThetaPDF(pdf_momTheta);
+ particle->setMomDeltaPhiPDF(pdf_momPhi);
+
+ // (8.) set some additional particle and simulation properties
+ double restMass = m_particleDataTable->particle(abs(pdg))->mass();
+ minEnergy = ( minEnergy > restMass ) ? minEnergy : restMass;
+ particle->setMinEnergy(minEnergy);
+ particle->setMaxNumParticles(maxNumParticles);
+ particle->setNumParticlesFactor(numParticlesFactor);
+ particle->setEnergyFactor(energyFactor);
+ particle->setPosAngleFactor(posAngleFactor);
+ particle->setMomAngleFactor(momAngleFactor);
+
+ // (9.) insert this PunchThroughParticle instance into the std::map class member
+ m_particles[pdg] = particle;
+
+ return StatusCode::SUCCESS;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+StatusCode ISF::PunchThroughTool::registerCorrelation(int pdgID1, int pdgID2,
+ double minCorrEnergy, double fullCorrEnergy)
+{
+ // find the given pdgs in the registered particle ids
+ std::map<int, PunchThroughParticle*>::iterator location1 = m_particles.find(pdgID1);
+ std::map<int, PunchThroughParticle*>::iterator location2 = m_particles.find(pdgID2);
+
+ // if at least one of the given pdgs was not registered yet -> return an error
+ if ( (location1 == m_particles.end()) || (location2 == m_particles.end()) )
+ return StatusCode::FAILURE;
+
+ // now look for the correlation histograms
+ std::stringstream name;
+ name << "NumExitCorrelations/x_PDG" << abs(pdgID1) << "__y_PDG" << abs(pdgID2) << "__lowE";
+ TH2F *hist1_2_lowE = (TH2F*)m_fileLookupTable->Get(name.str().c_str());
+ name.str("");
+ name << "NumExitCorrelations/x_PDG" << abs(pdgID1) << "__y_PDG" << abs(pdgID2) << "__highE";
+ TH2F *hist1_2_highE = (TH2F*)m_fileLookupTable->Get(name.str().c_str());
+ name.str("");
+ name << "NumExitCorrelations/x_PDG" << abs(pdgID2) << "__y_PDG" << abs(pdgID1) << "__lowE";
+ TH2F *hist2_1_lowE = (TH2F*)m_fileLookupTable->Get(name.str().c_str());
+ name.str("");
+ name << "NumExitCorrelations/x_PDG" << abs(pdgID2) << "__y_PDG" << abs(pdgID1) << "__highE";
+ TH2F *hist2_1_highE = (TH2F*)m_fileLookupTable->Get(name.str().c_str());
+ // check if the histograms exist
+ if ( (!hist1_2_lowE) || (!hist2_1_lowE) || (!hist1_2_highE) || (!hist2_1_highE) )
+ {
+ ATH_MSG_ERROR("[ punchthrough ] unable to retrieve the correlation data for PDG IDs " << pdgID1 << " and " << pdgID2);
+ return StatusCode::FAILURE;
+ }
+
+ // TODO: if only one of the two histograms exists, create the other one
+ // by mirroring the data
+
+ double lowE = getFloatAfterPatternInStr( hist1_2_lowE->GetTitle(), "elow_");
+ double midE = getFloatAfterPatternInStr( hist1_2_lowE->GetTitle(), "ehigh_");
+ //TODO: check if the same:
+ // double midE = getFloatAfterPatternInStr( hist1_2_lowE->GetTitle(), "elow_");
+ double upperE = getFloatAfterPatternInStr( hist1_2_highE->GetTitle(), "ehigh_");
+ // now store the correlation either way id1->id2 and id2->id1
+ m_particles[pdgID1]->setCorrelation(pdgID2, hist2_1_lowE, hist2_1_highE,
+ minCorrEnergy, fullCorrEnergy,
+ lowE, midE, upperE);
+
+ m_particles[pdgID2]->setCorrelation(pdgID1, hist1_2_lowE, hist1_2_highE,
+ minCorrEnergy, fullCorrEnergy,
+ lowE, midE, upperE);
+ return StatusCode::SUCCESS;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *======================================================================*/
+
+ISF::PDFcreator *ISF::PunchThroughTool::readLookuptablePDF(int pdg, std::string folderName)
+{
+ // will hold the PDFcreator class which will be returned at the end
+ // this will store the probability density functions for the number of punch-through particles
+ // (as functions of energy & eta of the incoming particle)
+ PDFcreator *pdf = 0;
+
+ // (1.) retrieve the distribution function
+ std::stringstream name;
+ name << folderName << pdg << "/function";
+ TF1 *func = (TF1*)m_fileLookupTable->Get(name.str().c_str());
+ if (! func )
+ {
+ ATH_MSG_ERROR("[ punchthrough ] unable to retrieve the PDF (" << folderName << ") from the lookuptable" );
+ return 0;
+ }
+ // store this function in the PDFcreator class
+ pdf = new PDFcreator(func, m_randomEngine);
+
+ // (2.) get this function maximum and minimum
+ std::stringstream namemin, namemax;
+ namemin << folderName << pdg << "/randmin";
+ TH1 *randmin = (TH1*)m_fileLookupTable->Get(namemin.str().c_str());
+ namemax << folderName << pdg << "/randmax";
+ TH1 *randmax = (TH1*)m_fileLookupTable->Get(namemax.str().c_str());
+
+ if ( !(randmin && randmax) )
+ {
+ ATH_MSG_ERROR("[ punchthrough ] unable to retrieve the PDF boundaries (" << folderName << ") from the lookuptable" );
+ delete pdf;
+ return 0;
+ }
+
+ // store this histogram in the PDFcreator class
+ pdf->setRange( randmin, randmax);
+
+ // (3.) get the parameters for the just received distribution function (stored in TH2F histograms)
+ for (int par = 0; par < func->GetNpar(); par++)
+ {
+ // prepare the name of the histogram in the file
+ name.str("");
+ name << folderName << pdg<< "/parameter" << par;
+
+ // get the histogram from the file
+ TH1 *hist = (TH1*)m_fileLookupTable->Get(name.str().c_str());
+
+ // if histogram does not exist -> error!
+ if(! hist)
+ {
+ ATH_MSG_ERROR( "[ punchthrough ] unable to retrieve all parameters for the distribution ("<< folderName << pdg << ")" );
+ delete pdf;
+ return 0;
+ }
+
+ // add this histogram to the PDFcreator class
+ pdf->addPar( hist );
+ }
+
+ return pdf;
+}
+
+/* =========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=========================================================================*/
+
+ISF::ISFParticle* ISF::PunchThroughTool::createExitPs( int pdg,
+ double energy, double theta, double phi,double momTheta, double momPhi) const
+{
+ // the intersection point with Calo-MS surface
+
+ Amg::Vector3D pos = propagator(theta,phi);
+
+ // set up the real punch-through particle at this position
+ // set up the momentum vector of this particle as a GlobalMomentum
+ // by using the given energy and mass of the particle and also using
+ // the given theta and phi
+
+ Amg::Vector3D mom;
+ double mass = m_particleDataTable->particle(abs(pdg))->mass();
+ Amg::setRThetaPhi( mom, sqrt(energy*energy - mass*mass), momTheta, momPhi);
+
+ double charge = m_particleDataTable->particle(abs(pdg))->charge();
+ // since the PDT table only has abs(PID) values for the charge
+ charge *= (pdg > 0.) ? 1. : -1.;
+
+ double pTime = 0; /** @TODO: fix */
+
+ ISF::ISFParticle* finalPar = new ISF::ISFParticle (pos, mom, mass, charge, pdg, pTime, *m_initPs, m_secBC);
+ finalPar->setNextGeoID( AtlasDetDescr::fAtlasMS);
+
+ // return the punch-through particle
+ return finalPar;
+}
+
+/*=========================================================================
+ * DESCRIPTION OF FUNCTION:
+ * ==> see headerfile
+ *=======================================================================*/
+
+double ISF::PunchThroughTool::getFloatAfterPatternInStr(const char *cstr, const char *cpattern)
+{
+ double num = 0.;
+
+ std::string str( cstr);
+ std::string pattern( cpattern);
+ size_t pos = str.find(cpattern);
+
+ if ( pos == std::string::npos)
+ {
+ ATH_MSG_WARNING("[ punchthrough ] unable to retrieve floating point number from string");
+ return -999999999999.;
+ }
+
+ std::istringstream iss( cstr+pos+pattern.length());
+ iss >> std::dec >> num;
+
+ return num;
+}
+
+Amg::Vector3D ISF::PunchThroughTool::propagator(double theta,double phi) const
+{
+ // phi, theta - direction of the punch-through particle coming into calo
+ //particle propagates inside the calorimeter along the straight line
+ //coordinates of this particles when exiting the calo (on calo-MS boundary)
+
+ double x, y, z, r;
+
+ // cylinders border angles
+ double theta1 = atan (R1/z1);
+ double theta2 = atan (R1/z2);
+ double theta3 = atan (R2/z2);
+ //where is the particle
+
+ if (theta >= 0 && theta < theta1)
+ {
+ z = z1;
+ r = fabs (z1*tan(theta));
+ }
+ else if (theta >= theta1 && theta < theta2)
+ {
+ z = R1/tan(theta);
+ r = R1;
+ }
+ else if (theta >= theta2 && theta < theta3)
+ {
+ z = z2;
+ r = fabs(z2*tan(theta));;
+ }
+ else if (theta >= theta3 && theta < (TMath::Pi()-theta3) )
+ {
+ z = R2/tan(theta);
+ r = R2;
+ }
+ else if (theta >= (TMath::Pi()-theta3) && theta < (TMath::Pi()-theta2) )
+ {
+ z = -z2;
+ r = fabs(z2*tan(theta));
+ }
+ else if (theta >= (TMath::Pi()-theta2) && theta < (TMath::Pi()-theta1) )
+ {
+ z = R1/tan(theta);
+ r = R1;
+ }
+ else if (theta >= (TMath::Pi()-theta1) && theta <= TMath::Pi() )
+ {
+ z = -z1;
+ r = fabs(z1*tan(theta));
+ }
+
+ //parallel universe
+ else
+ {
+ ATH_MSG_WARNING ( "Given theta angle is incorrect, setting particle position to (0, 0, 0)");
+ x = 0.0; y = 0.0; z = 0.0; r = 0.0;
+ }
+
+ x = r*cos(phi);
+ y = r*sin(phi);
+ Amg::Vector3D pos(x, y, z);
+
+ ATH_MSG_DEBUG("position of produced punch-through particle: x = "<< x <<" y = "<< y <<" z = "<< z<<" r = "<< pos.perp() <<"sqrt(x^2+y^2) = "<< sqrt(x*x+y*y) );
+ ATH_MSG_DEBUG("GeoID thinks: Calo: "<< m_geoIDSvc->inside(pos, AtlasDetDescr::fAtlasCalo) <<" MS: "<< m_geoIDSvc->inside(pos,AtlasDetDescr::fAtlasMS));
+
+ return pos;
+}
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/components/ISF_PunchThroughTools_entries.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/components/ISF_PunchThroughTools_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..d360f7ce0f163c34b948ac306e24a8a5252b7d38
--- /dev/null
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/components/ISF_PunchThroughTools_entries.cxx
@@ -0,0 +1,8 @@
+#include "GaudiKernel/DeclareFactoryEntries.h"
+#include "ISF_PunchThroughTools/PunchThroughTool.h"
+
+DECLARE_NAMESPACE_TOOL_FACTORY( ISF,PunchThroughTool)
+DECLARE_FACTORY_ENTRIES( ISF_PunchThroughTools ) {
+ DECLARE_NAMESPACE_TOOL( ISF,PunchThroughTool)
+}
+
diff --git a/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/components/ISF_PunchThroughTools_load.cxx b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/components/ISF_PunchThroughTools_load.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..428a49e65990035e8a7c73da3843dd7b59500657
--- /dev/null
+++ b/Simulation/ISF/ISF_FastCaloSim/ISF_PunchThroughTools/src/components/ISF_PunchThroughTools_load.cxx
@@ -0,0 +1,5 @@
+
+#include "GaudiKernel/LoadFactoryEntries.h"
+
+LOAD_FACTORY_ENTRIES( ISF_PunchThroughTools )
+