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