Par01/src/Par01EMShowerModel.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file Par01/src/Par01EMShowerModel.cc
0027 /// \brief Implementation of the Par01EMShowerModel class
0028 //
0029 //
0030 //
0031 #include "Par01EMShowerModel.hh"
0032
0033 #include "Par01EnergySpot.hh"
0034
0035 #include "G4Electron.hh"
0036 #include "G4Gamma.hh"
0037 #include "G4NistManager.hh"
0038 #include "G4PhysicalConstants.hh"
0039 #include "G4Positron.hh"
0040 #include "G4SystemOfUnits.hh"
0041 #include "G4TouchableHandle.hh"
0042 #include "G4TransportationManager.hh"
0043 #include "G4VSensitiveDetector.hh"
0044 #include "Randomize.hh"
0045
0046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0047
0048 Par01EMShowerModel::Par01EMShowerModel(G4String modelName, G4Region* envelope)
0049   : G4VFastSimulationModel(modelName, envelope)
0050 {
0051   fFakeStep = new G4Step();
0052   fFakePreStepPoint = fFakeStep->GetPreStepPoint();
0053   fFakePostStepPoint = fFakeStep->GetPostStepPoint();
0054   fTouchableHandle = new G4TouchableHistory();
0055   fpNavigator = new G4Navigator();
0056   fNaviSetup = false;
0057   fCsI = nullptr;
0058 }
0059
0060 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0061
0062 Par01EMShowerModel::Par01EMShowerModel(G4String modelName) : G4VFastSimulationModel(modelName)
0063 {
0064   fFakeStep = new G4Step();
0065   fFakePreStepPoint = fFakeStep->GetPreStepPoint();
0066   fFakePostStepPoint = fFakeStep->GetPostStepPoint();
0067   fTouchableHandle = new G4TouchableHistory();
0068   fpNavigator = new G4Navigator();
0069   fNaviSetup = false;
0070   fCsI = nullptr;
0071 }
0072
0073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0074
0075 Par01EMShowerModel::~Par01EMShowerModel()
0076 {
0077   delete fFakeStep;
0078   delete fpNavigator;
0079 }
0080
0081 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0082
0083 G4bool Par01EMShowerModel::IsApplicable(const G4ParticleDefinition& particleType)
0084 {
0085   return &particleType == G4Electron::ElectronDefinition()
0086          || &particleType == G4Positron::PositronDefinition()
0087          || &particleType == G4Gamma::GammaDefinition();
0088 }
0089
0090 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0091
0092 G4bool Par01EMShowerModel::ModelTrigger(const G4FastTrack& fastTrack)
0093 {
0094   // Applies the parameterisation above 100 MeV:
0095   return fastTrack.GetPrimaryTrack()->GetKineticEnergy() > 100 * MeV;
0096 }
0097
0098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0099
0100 void Par01EMShowerModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
0101 {
0102   // Kill the parameterised particle:
0103   fastStep.KillPrimaryTrack();
0104   fastStep.ProposePrimaryTrackPathLength(0.0);
0105   fastStep.ProposeTotalEnergyDeposited(fastTrack.GetPrimaryTrack()->GetKineticEnergy());
0106
0107   // split into "energy spots" energy according to the shower shape:
0108   Explode(fastTrack);
0109
0110   // and put those energy spots into the crystals:
0111   BuildDetectorResponse();
0112 }
0113
0114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0115
0116 void Par01EMShowerModel::Explode(const G4FastTrack& fastTrack)
0117 {
0118   //-----------------------------------------------------
0119   //
0120   //-----------------------------------------------------
0121
0122   // Reduced quantities:
0123   // -- critical energy in CsI:
0124   G4double Ec = 800 * MeV / (54. + 1.2);  // 54 = mean Z of CsI
0125   G4double Energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy();
0126   G4double y = Energy / Ec;
0127
0128   // compute value of parameter "a" of longitudinal profile, b assumed = 0.5
0129   G4double a, tmax, b(0.5), C;
0130   if (fastTrack.GetPrimaryTrack()->GetDefinition() == G4Gamma::GammaDefinition())
0131     C = 0.5;
0132   else
0133     C = -0.5;
0134   tmax = 1.0 * (std::log(y) + C);
0135   a = 1.0 + b * tmax;
0136
0137   // t : reduced quantity = z/X0:
0138   G4double t, bt;
0139   if (fCsI == nullptr) fCsI = G4NistManager::Instance()->FindOrBuildMaterial("G4_CESIUM_IODIDE");
0140   G4double X0 = fCsI->GetRadlen();
0141   // Moliere radius:
0142   G4double Es = 21 * MeV;
0143   G4double Rm = X0 * Es / Ec;
0144
0145   // axis of the shower, in global reference frame:
0146   G4ThreeVector xShower, yShower, zShower;
0147   zShower = fastTrack.GetPrimaryTrack()->GetMomentumDirection();
0148   xShower = zShower.orthogonal();
0149   yShower = zShower.cross(xShower);
0150   // starting point of the shower:
0151   G4ThreeVector sShower = fastTrack.GetPrimaryTrack()->GetPosition();
0152
0153   // We shoot 100 spots of energy:
0154   G4int nSpots = 100;
0155   G4double deposit = Energy / double(nSpots);
0156   Par01EnergySpot eSpot;
0157   eSpot.SetEnergy(deposit);
0158   G4ThreeVector ePoint;
0159   G4double z, r, phi;
0160
0161   feSpotList.clear();
0162   for (int i = 0; i < nSpots; i++) {
0163     // Longitudinal profile:
0164     // -- shoot z according to Gamma distribution:
0165     bt = G4RandGamma::shoot(a, 1.0);
0166     t = bt / b;
0167     z = t * X0;
0168
0169     // transverse profile:
0170     // we set 90% of energy in one Rm,
0171     // the rest between 1 and 3.5 Rm:
0172     G4double xr = G4UniformRand();
0173     if (xr < 0.9)
0174       r = xr / 0.9 * Rm;
0175     else
0176       r = ((xr - 0.9) / 0.1 * 2.5 + 1.0) * Rm;
0177     phi = G4UniformRand() * twopi;
0178
0179     // build the position:
0180     ePoint = sShower + z * zShower + r * std::cos(phi) * xShower + r * std::sin(phi) * yShower;
0181
0182     // and the energy spot:
0183     eSpot.SetPosition(ePoint);
0184
0185     // Records the eSpot:
0186     feSpotList.push_back(eSpot);
0187   }
0188 }
0189
0190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0191
0192 void Par01EMShowerModel::BuildDetectorResponse()
0193 {
0194   // Does the assignation of the energy spots to the sensitive volumes:
0195   for (size_t i = 0; i < feSpotList.size(); i++) {
0196     // Draw the energy spot:
0197     //      feSpotList[i].Draw();
0198     //      feSpotList[i].Print();
0199
0200     // "converts" the energy spot into the fake
0201     // G4Step to pass to sensitive detector:
0202     AssignSpotAndCallHit(feSpotList[i]);
0203   }
0204 }
0205
0206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0207
0208 void Par01EMShowerModel::AssignSpotAndCallHit(const Par01EnergySpot& eSpot)
0209 {
0210   //
0211   // "converts" the energy spot into the fake
0212   // G4Step to pass to sensitive detector:
0213   //
0214   FillFakeStep(eSpot);
0215
0216   //
0217   // call sensitive part: taken/adapted from the stepping:
0218   // Send G4Step information to Hit/Dig if the volume is sensitive
0219   //
0220   G4VPhysicalVolume* pCurrentVolume = fFakeStep->GetPreStepPoint()->GetPhysicalVolume();
0221   G4VSensitiveDetector* pSensitive;
0222
0223   if (pCurrentVolume != nullptr) {
0224     pSensitive = pCurrentVolume->GetLogicalVolume()->GetSensitiveDetector();
0225     if (pSensitive != nullptr) {
0226       pSensitive->Hit(fFakeStep);
0227     }
0228   }
0229 }
0230
0231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0232
0233 void Par01EMShowerModel::FillFakeStep(const Par01EnergySpot& eSpot)
0234 {
0235   //-----------------------------------------------------------
0236   // find in which volume the spot is.
0237   //-----------------------------------------------------------
0238   if (!fNaviSetup) {
0239     fpNavigator->SetWorldVolume(G4TransportationManager::GetTransportationManager()
0240                                   ->GetNavigatorForTracking()
0241                                   ->GetWorldVolume());
0242     fpNavigator->LocateGlobalPointAndUpdateTouchableHandle(
0243       eSpot.GetPosition(), G4ThreeVector(0., 0., 0.), fTouchableHandle, false);
0244     fNaviSetup = true;
0245   }
0246   else {
0247     fpNavigator->LocateGlobalPointAndUpdateTouchableHandle(
0248       eSpot.GetPosition(), G4ThreeVector(0., 0., 0.), fTouchableHandle);
0249   }
0250   //--------------------------------------
0251   // Fills attribute of the G4Step needed
0252   // by our sensitive detector:
0253   //-------------------------------------
0254   // set touchable volume at PreStepPoint:
0255   fFakePreStepPoint->SetTouchableHandle(fTouchableHandle);
0256   // set total energy deposit:
0257   fFakeStep->SetTotalEnergyDeposit(eSpot.GetEnergy());
0258 }