DDG4/plugins/Geant4P1ShowerModel.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013 //
0014 // Please note:
0015 //
0016 // These shower parametrizations come direct from the Geant4 example:
0017 // <geant4-source-dir>/examples/extended/parameterisations/Par01
0018 //
0019 // From the README:
0020 //        o Par01EMShowerModel which provides a crude
0021 //      parameterisation for e+/e-/gamma. This model
0022 //      is bound to the EM calorimeter.
0023 //
0024 //        o Par01PionShowerModel: an even more crude
0025 //      parameterisation for pi+/pi-. This model
0026 //      is bound to a ghost volume.
0027 //
0028 //==========================================================================
0029
0030 // Framework include files
0031 #include <DDG4/Geant4FastSimShowerModel.inl.h>
0032 #include <DDG4/Geant4FastSimSpot.h>
0033 #include <DDG4/Geant4Random.h>
0034
0035 // Geant4 include files
0036 #include <G4Gamma.hh>
0037 #include <G4SystemOfUnits.hh>
0038
0039 // C/C++ include files
0040
0041
0042 /// Namespace for the AIDA detector description toolkit
0043 namespace dd4hep  {
0044
0045   /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0046   namespace sim  {
0047
0048     ///===================================================================================================
0049     ///
0050     ///  Par01 electromagnetic shower model (e+, e-)
0051     ///
0052     ///===================================================================================================
0053
0054     /// Configuration structure for the fast simulation shower model Geant4FSShowerModel<par01_em_model>
0055     class par01_em_model  {
0056     public:
0057       G4FastSimHitMaker hitMaker          { };
0058       std::string       materialName      { };
0059       G4Material*       material          { nullptr };
0060       double            criticalEnergyRef { 800*MeV };
0061       double            criticalEnergy    { 800*MeV };
0062     };
0063
0064     /// Declare optional properties from embedded structure
0065     template <>
0066     void Geant4FSShowerModel<par01_em_model>::initialize()     {
0067       declareProperty("Material",       this->locals.materialName);
0068       declareProperty("CriticalEnergy", this->locals.criticalEnergyRef);
0069       this->m_applicablePartNames.emplace_back("e+");
0070       this->m_applicablePartNames.emplace_back("e-");
0071     }
0072
0073     /// Sensitive detector construction callback. Called at "ConstructSDandField()"
0074     template <>
0075     void Geant4FSShowerModel<par01_em_model>::constructSensitives(Geant4DetectorConstructionContext* ctxt)   {
0076       locals.material       = this->getMaterial(this->locals.materialName);
0077       locals.criticalEnergy = this->locals.criticalEnergyRef * (this->locals.material->GetZ() + 1.2);
0078       this->Geant4FastSimShowerModel::constructSensitives(ctxt);
0079     }
0080
0081     /// User callback to model the particle/energy shower
0082     template <>
0083     void Geant4FSShowerModel<par01_em_model>::modelShower(const G4FastTrack& track, G4FastStep& step)   {
0084       auto* primary = track.GetPrimaryTrack();
0085       // Kill the parameterised particle:
0086       this->killParticle(step, primary->GetKineticEnergy(), 0e0);
0087       //-----------------------------------------------------
0088       // split into "energy spots" energy according to the shower shape:
0089       // See Par01EMShowerModel::Explode(track);
0090       //-----------------------------------------------------
0091       G4FastHit hit;
0092       Geant4FastSimSpot spot(&hit, &track);
0093
0094       // Reduced quantities:   -- critical energy in material:
0095       G4double Ec      = this->locals.criticalEnergy;
0096       G4double Energy  = spot.kineticEnergy();
0097       G4double y       = Energy/Ec;
0098       // compute value of parameter "a" of longitudinal profile, b assumed = 0.5
0099       G4double b = 0.5, C = -0.5;
0100       if (spot.trackDefinition() == G4Gamma::GammaDefinition()) C = 0.5;
0101       G4double tmax    = 1.0 * (std::log(y) + C);
0102       G4double a       = 1.0 + b*tmax;
0103       // radiation length
0104       G4double X0      = this->locals.material->GetRadlen();
0105       // Moliere radius:
0106       G4double Es      = 21*MeV;
0107       G4double Rm      = X0*Es/Ec;
0108       // We shoot 100 spots of energy:
0109       G4int    nSpots  = 100;
0110       G4double deposit = Energy/double(nSpots);
0111
0112       // axis of the shower, in global reference frame:
0113       G4ThreeVector xShower, yShower, zShower;
0114       zShower = spot.particleLocalDirection();
0115       xShower = zShower.orthogonal();
0116       yShower = zShower.cross(xShower);
0117       // starting point of the shower:
0118       Geant4Random* rndm    = Geant4Random::instance();
0119       G4ThreeVector sShower = spot.particleLocalPosition();
0120       for (int i = 0; i < nSpots; i++)    {
0121     // Longitudinal profile: -- shoot z according to Gamma distribution:
0122     G4double bt  = rndm->gamma(a, 1e0);
0123     G4double t   = bt/b;                       // t : reduced quantity = z/X0:
0124     G4double z   = t*X0;
0125     // transverse profile: we set 90% of energy in one Rm, the rest between 1 and 3.5 Rm:
0126     G4double xr  = rndm->uniform(0e0,1e0);
0127     G4double phi = rndm->uniform(0e0,twopi);
0128     G4double r;
0129     if (xr < 0.9) r = xr/0.9*Rm;
0130     else r = ((xr - 0.9)/0.1*2.5 + 1.0)*Rm;
0131     // build the position:
0132     G4ThreeVector position = sShower + z*zShower + r*std::cos(phi)*xShower + r*std::sin(phi)*yShower;
0133     /// Process spot and call sensitive detector
0134     hit.SetPosition(position);
0135     hit.SetEnergy(deposit);
0136     this->locals.hitMaker.make(hit, track);
0137       }
0138     }
0139
0140     ///===================================================================================================
0141     ///
0142     ///  Par01 Pion shower model (pi+, pi-)
0143     ///
0144     ///===================================================================================================
0145
0146     /// Configuration structure for the fast simulation shower model Geant4FSShowerModel<par01_pion_model>
0147     class par01_pion_model  {
0148     public:
0149       G4FastSimHitMaker hitMaker { };
0150     };
0151
0152     /// Declare optional properties from embedded structure
0153     template <>
0154     void Geant4FSShowerModel<par01_pion_model>::initialize()     {
0155       this->m_applicablePartNames.emplace_back("pi+");
0156       this->m_applicablePartNames.emplace_back("pi-");
0157     }
0158
0159     /// User callback to determine if the shower creation should be triggered
0160     template <>
0161     bool Geant4FSShowerModel<par01_pion_model>::check_trigger(const G4FastTrack& /* track */)   {
0162       return true;
0163     }
0164
0165     /// User callback to model the particle/energy shower
0166     template <>
0167     void Geant4FSShowerModel<par01_pion_model>::modelShower(const G4FastTrack& track, G4FastStep& step)   {
0168       auto* primary = track.GetPrimaryTrack();
0169       // Kill the parameterised particle:
0170       this->killParticle(step, primary->GetKineticEnergy(), 0e0);
0171
0172       //-----------------------------------------------------
0173       // Non-physical shower generated: exp along z and transverse.
0174       //-----------------------------------------------------
0175       // center of the shower, we put at the middle of the ghost:
0176       G4ThreeVector pos = track.GetPrimaryTrackLocalPosition();
0177       G4ThreeVector dir = track.GetPrimaryTrackLocalDirection();
0178       G4double      distOut = track.GetEnvelopeSolid()->DistanceToOut(pos, dir);
0179       G4ThreeVector showerCenter = pos + (distOut/2.)*dir;
0180
0181       showerCenter = track.GetInverseAffineTransformation()->TransformPoint(showerCenter);
0182
0183       // axis of the shower, in global reference frame:
0184       G4ThreeVector zShower = primary->GetMomentumDirection();
0185       G4ThreeVector xShower = zShower.orthogonal();
0186       G4ThreeVector yShower = zShower.cross(xShower);
0187
0188       // shoot the energy spots:
0189       G4double      Energy  = primary->GetKineticEnergy();
0190       G4int         nSpot   = 50;
0191       G4double      deposit = Energy/double(nSpot);
0192       Geant4Random* rndm    = Geant4Random::instance();
0193       for (int i = 0; i < nSpot; i++)  {
0194     double z   = rndm->gauss(0, 20*cm);
0195     double r   = rndm->gauss(0, 10*cm);
0196     double phi = rndm->uniform(0e0, twopi);
0197     G4ThreeVector position = showerCenter + z*zShower + r*std::cos(phi)*xShower + r*std::sin(phi)*yShower;
0198     /// Process spot and call sensitive detector
0199     G4FastHit hit(position, deposit);
0200     this->locals.hitMaker.make(hit, track);
0201       }
0202     }
0203
0204     typedef Geant4FSShowerModel<par01_em_model>   Geant4Par01EMShowerModel;
0205     typedef Geant4FSShowerModel<par01_pion_model> Geant4Par01PionShowerModel;
0206   }
0207 }
0208
0209 #include <DDG4/Factories.h>
0210 DECLARE_GEANT4ACTION_NS(dd4hep::sim,Geant4Par01EMShowerModel)
0211 DECLARE_GEANT4ACTION_NS(dd4hep::sim,Geant4Par01PionShowerModel)