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