| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 /// \file Par01/src/Par01EMShowerModel.cc
27 /// \brief Implementation of the Par01EMShowerModel class
28 //
29 //
30 //
31 #include "Par01EMShowerModel.hh"
32
33 #include "Par01EnergySpot.hh"
34
35 #include "G4Electron.hh"
36 #include "G4Gamma.hh"
37 #include "G4NistManager.hh"
38 #include "G4PhysicalConstants.hh"
39 #include "G4Positron.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4TouchableHandle.hh"
42 #include "G4TransportationManager.hh"
43 #include "G4VSensitiveDetector.hh"
44 #include "Randomize.hh"
45
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47
48 Par01EMShowerModel::Par01EMShowerModel(G4String modelName, G4Region* envelope)
49 : G4VFastSimulationModel(modelName, envelope)
50 {
51 fFakeStep = new G4Step();
52 fFakePreStepPoint = fFakeStep->GetPreStepPoint();
53 fFakePostStepPoint = fFakeStep->GetPostStepPoint();
54 fTouchableHandle = new G4TouchableHistory();
55 fpNavigator = new G4Navigator();
56 fNaviSetup = false;
57 fCsI = nullptr;
58 }
59
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61
62 Par01EMShowerModel::Par01EMShowerModel(G4String modelName) : G4VFastSimulationModel(modelName)
63 {
64 fFakeStep = new G4Step();
65 fFakePreStepPoint = fFakeStep->GetPreStepPoint();
66 fFakePostStepPoint = fFakeStep->GetPostStepPoint();
67 fTouchableHandle = new G4TouchableHistory();
68 fpNavigator = new G4Navigator();
69 fNaviSetup = false;
70 fCsI = nullptr;
71 }
72
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74
75 Par01EMShowerModel::~Par01EMShowerModel()
76 {
77 delete fFakeStep;
78 delete fpNavigator;
79 }
80
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82
83 G4bool Par01EMShowerModel::IsApplicable(const G4ParticleDefinition& particleType)
84 {
85 return &particleType == G4Electron::ElectronDefinition()
86 || &particleType == G4Positron::PositronDefinition()
87 || &particleType == G4Gamma::GammaDefinition();
88 }
89
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91
92 G4bool Par01EMShowerModel::ModelTrigger(const G4FastTrack& fastTrack)
93 {
94 // Applies the parameterisation above 100 MeV:
95 return fastTrack.GetPrimaryTrack()->GetKineticEnergy() > 100 * MeV;
96 }
97
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
99
100 void Par01EMShowerModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
101 {
102 // Kill the parameterised particle:
103 fastStep.KillPrimaryTrack();
104 fastStep.ProposePrimaryTrackPathLength(0.0);
105 fastStep.ProposeTotalEnergyDeposited(fastTrack.GetPrimaryTrack()->GetKineticEnergy());
106
107 // split into "energy spots" energy according to the shower shape:
108 Explode(fastTrack);
109
110 // and put those energy spots into the crystals:
111 BuildDetectorResponse();
112 }
113
114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
115
116 void Par01EMShowerModel::Explode(const G4FastTrack& fastTrack)
117 {
118 //-----------------------------------------------------
119 //
120 //-----------------------------------------------------
121
122 // Reduced quantities:
123 // -- critical energy in CsI:
124 G4double Ec = 800 * MeV / (54. + 1.2); // 54 = mean Z of CsI
125 G4double Energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy();
126 G4double y = Energy / Ec;
127
128 // compute value of parameter "a" of longitudinal profile, b assumed = 0.5
129 G4double a, tmax, b(0.5), C;
130 if (fastTrack.GetPrimaryTrack()->GetDefinition() == G4Gamma::GammaDefinition())
131 C = 0.5;
132 else
133 C = -0.5;
134 tmax = 1.0 * (std::log(y) + C);
135 a = 1.0 + b * tmax;
136
137 // t : reduced quantity = z/X0:
138 G4double t, bt;
139 if (fCsI == nullptr) fCsI = G4NistManager::Instance()->FindOrBuildMaterial("G4_CESIUM_IODIDE");
140 G4double X0 = fCsI->GetRadlen();
141 // Moliere radius:
142 G4double Es = 21 * MeV;
143 G4double Rm = X0 * Es / Ec;
144
145 // axis of the shower, in global reference frame:
146 G4ThreeVector xShower, yShower, zShower;
147 zShower = fastTrack.GetPrimaryTrack()->GetMomentumDirection();
148 xShower = zShower.orthogonal();
149 yShower = zShower.cross(xShower);
150 // starting point of the shower:
151 G4ThreeVector sShower = fastTrack.GetPrimaryTrack()->GetPosition();
152
153 // We shoot 100 spots of energy:
154 G4int nSpots = 100;
155 G4double deposit = Energy / double(nSpots);
156 Par01EnergySpot eSpot;
157 eSpot.SetEnergy(deposit);
158 G4ThreeVector ePoint;
159 G4double z, r, phi;
160
161 feSpotList.clear();
162 for (int i = 0; i < nSpots; i++) {
163 // Longitudinal profile:
164 // -- shoot z according to Gamma distribution:
165 bt = G4RandGamma::shoot(a, 1.0);
166 t = bt / b;
167 z = t * X0;
168
169 // transverse profile:
170 // we set 90% of energy in one Rm,
171 // the rest between 1 and 3.5 Rm:
172 G4double xr = G4UniformRand();
173 if (xr < 0.9)
174 r = xr / 0.9 * Rm;
175 else
176 r = ((xr - 0.9) / 0.1 * 2.5 + 1.0) * Rm;
177 phi = G4UniformRand() * twopi;
178
179 // build the position:
180 ePoint = sShower + z * zShower + r * std::cos(phi) * xShower + r * std::sin(phi) * yShower;
181
182 // and the energy spot:
183 eSpot.SetPosition(ePoint);
184
185 // Records the eSpot:
186 feSpotList.push_back(eSpot);
187 }
188 }
189
190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
191
192 void Par01EMShowerModel::BuildDetectorResponse()
193 {
194 // Does the assignation of the energy spots to the sensitive volumes:
195 for (size_t i = 0; i < feSpotList.size(); i++) {
196 // Draw the energy spot:
197 // feSpotList[i].Draw();
198 // feSpotList[i].Print();
199
200 // "converts" the energy spot into the fake
201 // G4Step to pass to sensitive detector:
202 AssignSpotAndCallHit(feSpotList[i]);
203 }
204 }
205
206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
207
208 void Par01EMShowerModel::AssignSpotAndCallHit(const Par01EnergySpot& eSpot)
209 {
210 //
211 // "converts" the energy spot into the fake
212 // G4Step to pass to sensitive detector:
213 //
214 FillFakeStep(eSpot);
215
216 //
217 // call sensitive part: taken/adapted from the stepping:
218 // Send G4Step information to Hit/Dig if the volume is sensitive
219 //
220 G4VPhysicalVolume* pCurrentVolume = fFakeStep->GetPreStepPoint()->GetPhysicalVolume();
221 G4VSensitiveDetector* pSensitive;
222
223 if (pCurrentVolume != nullptr) {
224 pSensitive = pCurrentVolume->GetLogicalVolume()->GetSensitiveDetector();
225 if (pSensitive != nullptr) {
226 pSensitive->Hit(fFakeStep);
227 }
228 }
229 }
230
231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
232
233 void Par01EMShowerModel::FillFakeStep(const Par01EnergySpot& eSpot)
234 {
235 //-----------------------------------------------------------
236 // find in which volume the spot is.
237 //-----------------------------------------------------------
238 if (!fNaviSetup) {
239 fpNavigator->SetWorldVolume(G4TransportationManager::GetTransportationManager()
240 ->GetNavigatorForTracking()
241 ->GetWorldVolume());
242 fpNavigator->LocateGlobalPointAndUpdateTouchableHandle(
243 eSpot.GetPosition(), G4ThreeVector(0., 0., 0.), fTouchableHandle, false);
244 fNaviSetup = true;
245 }
246 else {
247 fpNavigator->LocateGlobalPointAndUpdateTouchableHandle(
248 eSpot.GetPosition(), G4ThreeVector(0., 0., 0.), fTouchableHandle);
249 }
250 //--------------------------------------
251 // Fills attribute of the G4Step needed
252 // by our sensitive detector:
253 //-------------------------------------
254 // set touchable volume at PreStepPoint:
255 fFakePreStepPoint->SetTouchableHandle(fTouchableHandle);
256 // set total energy deposit:
257 fFakeStep->SetTotalEnergyDeposit(eSpot.GetEnergy());
258 }
259