| 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 //
>> 27 //
26 ////////////////////////////////////////////// 28 ////////////////////////////////////////////////////////////////////////
27 // Cerenkov Radiation Class Implementation 29 // Cerenkov Radiation Class Implementation
28 ////////////////////////////////////////////// 30 ////////////////////////////////////////////////////////////////////////
29 // 31 //
30 // File: G4Cerenkov.cc 32 // File: G4Cerenkov.cc
31 // Description: Discrete Process -- Generation 33 // Description: Discrete Process -- Generation of Cerenkov Photons
32 // Version: 2.1 34 // Version: 2.1
33 // Created: 1996-02-21 35 // Created: 1996-02-21
34 // Author: Juliet Armstrong 36 // Author: Juliet Armstrong
35 // Updated: 2007-09-30 by Peter Gumplinger 37 // Updated: 2007-09-30 by Peter Gumplinger
36 // > change inheritance to G4VDis 38 // > change inheritance to G4VDiscreteProcess
37 // GetContinuousStepLimit -> GetM 39 // GetContinuousStepLimit -> GetMeanFreePath (StronglyForced)
38 // AlongStepDoIt -> PostStepDoIt 40 // AlongStepDoIt -> PostStepDoIt
39 // 2005-08-17 by Peter Gumplinger 41 // 2005-08-17 by Peter Gumplinger
40 // > change variable name MeanNum 42 // > change variable name MeanNumPhotons -> MeanNumberOfPhotons
41 // 2005-07-28 by Peter Gumplinger 43 // 2005-07-28 by Peter Gumplinger
42 // > add G4ProcessType to constru 44 // > add G4ProcessType to constructor
43 // 2001-09-17, migration of Mater 45 // 2001-09-17, migration of Materials to pure STL (mma)
44 // 2000-11-12 by Peter Gumplinger 46 // 2000-11-12 by Peter Gumplinger
45 // > add check on CerenkovAngleIn 47 // > add check on CerenkovAngleIntegrals->IsFilledVectorExist()
46 // in method GetAverageNumberOfPh 48 // in method GetAverageNumberOfPhotons
47 // > and a test for MeanNumberOfP 49 // > and a test for MeanNumberOfPhotons <= 0.0 in DoIt
48 // 2000-09-18 by Peter Gumplinger 50 // 2000-09-18 by Peter Gumplinger
49 // > change: aSecondaryPosition=x 51 // > change: aSecondaryPosition=x0+rand*aStep.GetDeltaPosition();
50 // aSecondaryTrack->Set 52 // aSecondaryTrack->SetTouchable(0);
51 // 1999-10-29 by Peter Gumplinger 53 // 1999-10-29 by Peter Gumplinger
52 // > change: == into <= in GetCon 54 // > change: == into <= in GetContinuousStepLimit
53 // 1997-08-08 by Peter Gumplinger 55 // 1997-08-08 by Peter Gumplinger
54 // > add protection against /0 56 // > add protection against /0
55 // > G4MaterialPropertiesTable; n 57 // > G4MaterialPropertiesTable; new physics/tracking scheme
56 // 58 //
>> 59 // mail: gum@triumf.ca
>> 60 //
57 ////////////////////////////////////////////// 61 ////////////////////////////////////////////////////////////////////////
58 62
59 #include "G4Cerenkov.hh" <<
60 <<
61 #include "G4ios.hh" 63 #include "G4ios.hh"
>> 64 #include "G4PhysicalConstants.hh"
>> 65 #include "G4SystemOfUnits.hh"
>> 66 #include "G4Poisson.hh"
>> 67 #include "G4EmProcessSubType.hh"
>> 68
62 #include "G4LossTableManager.hh" 69 #include "G4LossTableManager.hh"
63 #include "G4Material.hh" <<
64 #include "G4MaterialCutsCouple.hh" 70 #include "G4MaterialCutsCouple.hh"
65 #include "G4MaterialPropertiesTable.hh" <<
66 #include "G4OpticalParameters.hh" <<
67 #include "G4OpticalPhoton.hh" <<
68 #include "G4ParticleDefinition.hh" 71 #include "G4ParticleDefinition.hh"
69 #include "G4ParticleMomentum.hh" <<
70 #include "G4PhysicalConstants.hh" <<
71 #include "G4PhysicsFreeVector.hh" <<
72 #include "G4Poisson.hh" <<
73 #include "G4SystemOfUnits.hh" <<
74 #include "G4ThreeVector.hh" <<
75 #include "Randomize.hh" <<
76 #include "G4PhysicsModelCatalog.hh" <<
77 72
78 //....oooOO0OOooo........oooOO0OOooo........oo << 73 #include "G4Cerenkov.hh"
>> 74
79 G4Cerenkov::G4Cerenkov(const G4String& process 75 G4Cerenkov::G4Cerenkov(const G4String& processName, G4ProcessType type)
80 : G4VProcess(processName, type) << 76 : G4VProcess(processName, type),
81 , fNumPhotons(0) << 77 fTrackSecondariesFirst(false),
>> 78 fMaxBetaChange(0.0),
>> 79 fMaxPhotons(0),
>> 80 fStackingFlag(true),
>> 81 fNumPhotons(0)
82 { 82 {
83 secID = G4PhysicsModelCatalog::GetModelID("m <<
84 SetProcessSubType(fCerenkov); 83 SetProcessSubType(fCerenkov);
85 84
86 thePhysicsTable = nullptr; 85 thePhysicsTable = nullptr;
87 86
88 if(verboseLevel > 0) << 87 if (verboseLevel>0) {
89 { << 88 G4cout << GetProcessName() << " is created " << G4endl;
90 G4cout << GetProcessName() << " is created <<
91 } 89 }
92 Initialise(); <<
93 } 90 }
94 91
95 //....oooOO0OOooo........oooOO0OOooo........oo <<
96 G4Cerenkov::~G4Cerenkov() 92 G4Cerenkov::~G4Cerenkov()
97 { 93 {
98 if(thePhysicsTable != nullptr) << 94 if (thePhysicsTable != nullptr) {
99 { << 95 thePhysicsTable->clearAndDestroy();
100 thePhysicsTable->clearAndDestroy(); << 96 delete thePhysicsTable;
101 delete thePhysicsTable; <<
102 } 97 }
103 } 98 }
104 99
105 void G4Cerenkov::ProcessDescription(std::ostre << 100 G4bool G4Cerenkov::IsApplicable(const G4ParticleDefinition& aParticleType)
106 { 101 {
107 out << "The Cerenkov effect simulates optica << 102 return (aParticleType.GetPDGCharge() != 0.0 &&
108 out << "passage of charged particles through << 103 aParticleType.GetPDGMass() != 0.0 &&
109 out << "to have the property RINDEX (refract << 104 aParticleType.GetParticleName() != "chargedgeantino" &&
110 G4VProcess::DumpInfo(); << 105 !aParticleType.IsShortLived() ) ? true : false;
>> 106 }
111 107
112 G4OpticalParameters* params = G4OpticalParam << 108 void G4Cerenkov::SetTrackSecondariesFirst(const G4bool state)
113 out << "Maximum beta change per step: " << p << 109 {
114 out << "Maximum photons per step: " << param << 110 fTrackSecondariesFirst = state;
115 out << "Track secondaries first: " <<
116 << params->GetCerenkovTrackSecondariesFi <<
117 out << "Stack photons: " << params->GetCeren <<
118 out << "Verbose level: " << params->GetCeren <<
119 } 111 }
120 112
121 //....oooOO0OOooo........oooOO0OOooo........oo << 113 void G4Cerenkov::SetMaxBetaChangePerStep(const G4double value)
122 G4bool G4Cerenkov::IsApplicable(const G4Partic <<
123 { 114 {
124 return (aParticleType.GetPDGCharge() != 0.0 << 115 fMaxBetaChange = value*CLHEP::perCent;
125 aParticleType.GetPDGMass() != 0.0 && <<
126 aParticleType.GetParticleName() != " <<
127 !aParticleType.IsShortLived()) <<
128 ? true <<
129 : false; <<
130 } 116 }
131 117
132 //....oooOO0OOooo........oooOO0OOooo........oo << 118 void G4Cerenkov::SetMaxNumPhotonsPerStep(const G4int NumPhotons)
133 void G4Cerenkov::Initialise() <<
134 { 119 {
135 G4OpticalParameters* params = G4OpticalParam << 120 fMaxPhotons = NumPhotons;
136 SetMaxBetaChangePerStep(params->GetCerenkovM <<
137 SetMaxNumPhotonsPerStep(params->GetCerenkovM <<
138 SetTrackSecondariesFirst(params->GetCerenkov <<
139 SetStackPhotons(params->GetCerenkovStackPhot <<
140 SetVerboseLevel(params->GetCerenkovVerboseLe <<
141 } 121 }
142 122
143 //....oooOO0OOooo........oooOO0OOooo........oo <<
144 void G4Cerenkov::BuildPhysicsTable(const G4Par 123 void G4Cerenkov::BuildPhysicsTable(const G4ParticleDefinition&)
145 { 124 {
146 if(thePhysicsTable) << 125 if (!thePhysicsTable) BuildThePhysicsTable();
147 return; <<
148 <<
149 const G4MaterialTable* theMaterialTable = G4 <<
150 std::size_t numOfMaterials = G4 <<
151 <<
152 thePhysicsTable = new G4PhysicsTable(numOfMa <<
153 <<
154 // loop over materials <<
155 for(std::size_t i = 0; i < numOfMaterials; + <<
156 { <<
157 G4PhysicsFreeVector* cerenkovIntegral = nu <<
158 <<
159 // Retrieve vector of refraction indices f <<
160 // from the material's optical properties <<
161 G4Material* aMaterial = (*theMate <<
162 G4MaterialPropertiesTable* MPT = aMaterial <<
163 <<
164 if(MPT) <<
165 { <<
166 cerenkovIntegral <<
167 G4MaterialPropertyVector* refractiveInde <<
168 <<
169 if(refractiveIndex) <<
170 { <<
171 // Retrieve the first refraction index <<
172 // of (photon energy, refraction index <<
173 G4double currentRI = (*refractiveIndex <<
174 if(currentRI > 1.0) <<
175 { <<
176 // Create first (photon energy, Cere <<
177 G4double currentPM = refractiveInde <<
178 G4double currentCAI = 0.0; <<
179 <<
180 cerenkovIntegral->InsertValues(curre <<
181 <<
182 // Set previous values to current on <<
183 G4double prevPM = currentPM; <<
184 G4double prevCAI = currentCAI; <<
185 G4double prevRI = currentRI; <<
186 <<
187 // loop over all (photon energy, ref <<
188 // pairs stored for this material <<
189 for(std::size_t ii = 1; ii < refract <<
190 { <<
191 currentRI = (*refractiveIndex)[ii <<
192 currentPM = refractiveIndex->Ener <<
193 currentCAI = prevCAI + (currentPM <<
194 (1.0 / (p <<
195 1.0 / (c <<
196 <<
197 cerenkovIntegral->InsertValues(cur <<
198 <<
199 prevPM = currentPM; <<
200 prevCAI = currentCAI; <<
201 prevRI = currentRI; <<
202 } <<
203 } <<
204 } <<
205 } <<
206 <<
207 // The Cerenkov integral for a given mater <<
208 // thePhysicsTable according to the positi <<
209 // the material table. <<
210 thePhysicsTable->insertAt(i, cerenkovInteg <<
211 } <<
212 } 126 }
213 127
214 //....oooOO0OOooo........oooOO0OOooo........oo << 128 // PostStepDoIt
215 G4VParticleChange* G4Cerenkov::PostStepDoIt(co << 129 // -------------
216 co << 130 //
>> 131 G4VParticleChange*
>> 132 G4Cerenkov::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
>> 133
217 // This routine is called for each tracking St 134 // This routine is called for each tracking Step of a charged particle
218 // in a radiator. A Poisson-distributed number 135 // in a radiator. A Poisson-distributed number of photons is generated
219 // according to the Cerenkov formula, distribu 136 // according to the Cerenkov formula, distributed evenly along the track
220 // segment and uniformly azimuth w.r.t. the pa 137 // segment and uniformly azimuth w.r.t. the particle direction. The
221 // parameters are then transformed into the Ma 138 // parameters are then transformed into the Master Reference System, and
222 // they are added to the particle change. 139 // they are added to the particle change.
223 140
224 { 141 {
>> 142 ////////////////////////////////////////////////////
>> 143 // Should we ensure that the material is dispersive?
>> 144 ////////////////////////////////////////////////////
>> 145
225 aParticleChange.Initialize(aTrack); 146 aParticleChange.Initialize(aTrack);
226 147
227 const G4DynamicParticle* aParticle = aTrack. 148 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
228 const G4Material* aMaterial = aTrack. << 149 const G4Material* aMaterial = aTrack.GetMaterial();
229 150
230 G4StepPoint* pPreStepPoint = aStep.GetPreSt 151 G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
231 G4StepPoint* pPostStepPoint = aStep.GetPostS 152 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
232 153
233 G4ThreeVector x0 = pPreStepPoint->GetPositio 154 G4ThreeVector x0 = pPreStepPoint->GetPosition();
234 G4ThreeVector p0 = aStep.GetDeltaPosition(). 155 G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
235 G4double t0 = pPreStepPoint->GetGlobalT << 156 G4double t0 = pPreStepPoint->GetGlobalTime();
236 157
237 G4MaterialPropertiesTable* MPT = aMaterial-> << 158 G4MaterialPropertiesTable* aMaterialPropertiesTable =
238 if(!MPT) << 159 aMaterial->GetMaterialPropertiesTable();
239 return pParticleChange; << 160 if (!aMaterialPropertiesTable) return pParticleChange;
240 << 161
241 G4MaterialPropertyVector* Rindex = MPT->GetP << 162 G4MaterialPropertyVector* Rindex =
242 if(!Rindex) << 163 aMaterialPropertiesTable->GetProperty(kRINDEX);
243 return pParticleChange; << 164 if (!Rindex) return pParticleChange;
244 165
>> 166 // particle charge
245 G4double charge = aParticle->GetDefinition() 167 G4double charge = aParticle->GetDefinition()->GetPDGCharge();
246 168
247 G4double beta1 = pPreStepPoint->GetBeta(); << 169 // particle beta
248 G4double beta2 = pPostStepPoint->GetBeta(); << 170 G4double beta = (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta())*0.5;
249 G4double beta = (beta1 + beta2) * 0.5; <<
250 171
251 G4double MeanNumberOfPhotons = << 172 //fNumPhotons = 0; // in PostStepGetPhysicalInteractionLength()
252 GetAverageNumberOfPhotons(charge, beta, aM << 173
253 G4double MeanNumberOfPhotons1 = << 174 G4double MeanNumberOfPhotons =
254 GetAverageNumberOfPhotons(charge, beta1, a << 175 GetAverageNumberOfPhotons(charge,beta,aMaterial,Rindex);
255 G4double MeanNumberOfPhotons2 = << 176
256 GetAverageNumberOfPhotons(charge, beta2, a << 177 if (MeanNumberOfPhotons <= 0.0) {
>> 178
>> 179 // return unchanged particle and no secondaries
>> 180
>> 181 aParticleChange.SetNumberOfSecondaries(0);
>> 182
>> 183 return pParticleChange;
>> 184
>> 185 }
>> 186
>> 187 G4double step_length = aStep.GetStepLength();
>> 188
>> 189 MeanNumberOfPhotons = MeanNumberOfPhotons * step_length;
>> 190
>> 191 fNumPhotons = (G4int) G4Poisson(MeanNumberOfPhotons);
>> 192
>> 193 if ( fNumPhotons <= 0 || !fStackingFlag ) {
>> 194
>> 195 // return unchanged particle and no secondaries
>> 196
>> 197 aParticleChange.SetNumberOfSecondaries(0);
>> 198
>> 199 return pParticleChange;
257 200
258 if(MeanNumberOfPhotons <= 0.0) <<
259 { <<
260 // return unchanged particle and no second <<
261 aParticleChange.SetNumberOfSecondaries(0); <<
262 return pParticleChange; <<
263 } <<
264 <<
265 MeanNumberOfPhotons *= aStep.GetStepLength() <<
266 fNumPhotons = (G4int) G4Poisson(Mean <<
267 <<
268 // third condition added to prevent infinite <<
269 // see bugzilla 2555 <<
270 if(fNumPhotons <= 0 || !fStackingFlag || <<
271 std::max(MeanNumberOfPhotons1, MeanNumber <<
272 { <<
273 // return unchanged particle and no second <<
274 aParticleChange.SetNumberOfSecondaries(0); <<
275 return pParticleChange; <<
276 } 201 }
277 202
278 //////////////////////////////////////////// 203 ////////////////////////////////////////////////////////////////
>> 204
279 aParticleChange.SetNumberOfSecondaries(fNumP 205 aParticleChange.SetNumberOfSecondaries(fNumPhotons);
280 206
281 if(fTrackSecondariesFirst) << 207 if (fTrackSecondariesFirst) {
282 { << 208 if (aTrack.GetTrackStatus() == fAlive )
283 if(aTrack.GetTrackStatus() == fAlive) << 209 aParticleChange.ProposeTrackStatus(fSuspend);
284 aParticleChange.ProposeTrackStatus(fSusp <<
285 } 210 }
286 211
287 //////////////////////////////////////////// 212 ////////////////////////////////////////////////////////////////
288 G4double Pmin = Rindex->Energy(0); <<
289 G4double Pmax = Rindex->GetMaxEnergy(); <<
290 G4double dp = Pmax - Pmin; <<
291 213
292 G4double nMax = Rindex->GetMaxValue() << 214 G4double Pmin = Rindex->GetMinLowEdgeEnergy();
293 G4double BetaInverse = 1. / beta; << 215 G4double Pmax = Rindex->GetMaxLowEdgeEnergy();
>> 216 G4double dp = Pmax - Pmin;
294 217
295 G4double maxCos = BetaInverse / nMax; << 218 G4double nMax = Rindex->GetMaxValue();
>> 219
>> 220 G4double BetaInverse = 1./beta;
>> 221
>> 222 G4double maxCos = BetaInverse / nMax;
296 G4double maxSin2 = (1.0 - maxCos) * (1.0 + m 223 G4double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos);
297 224
298 for(G4int i = 0; i < fNumPhotons; ++i) << 225 G4double beta1 = pPreStepPoint ->GetBeta();
299 { << 226 G4double beta2 = pPostStepPoint->GetBeta();
300 // Determine photon energy << 227
301 G4double rand; << 228 G4double MeanNumberOfPhotons1 =
302 G4double sampledEnergy, sampledRI; << 229 GetAverageNumberOfPhotons(charge,beta1,aMaterial,Rindex);
303 G4double cosTheta, sin2Theta; << 230 G4double MeanNumberOfPhotons2 =
304 << 231 GetAverageNumberOfPhotons(charge,beta2,aMaterial,Rindex);
305 // sample an energy << 232
306 do << 233 for (G4int i = 0; i < fNumPhotons; i++) {
307 { << 234
308 rand = G4UniformRand(); << 235 // Determine photon energy
309 sampledEnergy = Pmin + rand * dp; << 236
310 sampledRI = Rindex->Value(sampledEne << 237 G4double rand;
311 cosTheta = BetaInverse / sampledRI; << 238 G4double sampledEnergy, sampledRI;
312 << 239 G4double cosTheta, sin2Theta;
313 sin2Theta = (1.0 - cosTheta) * (1.0 + co << 240
314 rand = G4UniformRand(); << 241 // sample an energy
315 << 242
316 // Loop checking, 07-Aug-2015, Vladimir << 243 do {
317 } while(rand * maxSin2 > sin2Theta); << 244 rand = G4UniformRand();
318 << 245 sampledEnergy = Pmin + rand * dp;
319 // Create photon momentum direction vector << 246 sampledRI = Rindex->Value(sampledEnergy);
320 // with respect to the coordinate system w << 247 cosTheta = BetaInverse / sampledRI;
321 // direction is aligned with the z axis << 248
322 rand = G4UniformRand(); << 249 sin2Theta = (1.0 - cosTheta)*(1.0 + cosTheta);
323 G4double phi = twopi * rand; << 250 rand = G4UniformRand();
324 G4double sinPhi = std::sin(phi); << 251
325 G4double cosPhi = std::cos(phi); << 252 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
326 G4double sinTheta = std::sqrt(sin2Theta); << 253 } while (rand*maxSin2 > sin2Theta);
327 G4ParticleMomentum photonMomentum(sinTheta << 254
328 cosTheta << 255 // Generate random position of photon on cone surface
329 << 256 // defined by Theta
330 // Rotate momentum direction back to globa << 257
331 photonMomentum.rotateUz(p0); << 258 rand = G4UniformRand();
332 << 259
333 // Determine polarization of new photon << 260 G4double phi = twopi*rand;
334 G4ThreeVector photonPolarization(cosTheta << 261 G4double sinPhi = std::sin(phi);
335 -sinTheta << 262 G4double cosPhi = std::cos(phi);
336 << 263
337 // Rotate back to original coord system << 264 // calculate x,y, and z components of photon energy
338 photonPolarization.rotateUz(p0); << 265 // (in coord system with primary particle direction
339 << 266 // aligned with the z axis)
340 // Generate a new photon: << 267
341 auto aCerenkovPhoton = << 268 G4double sinTheta = std::sqrt(sin2Theta);
342 new G4DynamicParticle(G4OpticalPhoton::O << 269 G4double px = sinTheta*cosPhi;
343 << 270 G4double py = sinTheta*sinPhi;
344 aCerenkovPhoton->SetPolarization(photonPol << 271 G4double pz = cosTheta;
345 aCerenkovPhoton->SetKineticEnergy(sampledE << 272
346 << 273 // Create photon momentum direction vector
347 G4double NumberOfPhotons, N; << 274 // The momentum direction is still with respect
348 << 275 // to the coordinate system where the primary
349 do << 276 // particle direction is aligned with the z axis
350 { << 277
351 rand = G4UniformRand(); << 278 G4ParticleMomentum photonMomentum(px, py, pz);
352 NumberOfPhotons = MeanNumberOfPhotons1 - << 279
353 rand * (MeanNumberOfPh << 280 // Rotate momentum direction back to global reference
354 N = << 281 // system
355 G4UniformRand() * std::max(MeanNumberO << 282
356 // Loop checking, 07-Aug-2015, Vladimir << 283 photonMomentum.rotateUz(p0);
357 } while(N > NumberOfPhotons); << 284
358 << 285 // Determine polarization of new photon
359 G4double delta = rand * aStep.GetStepLengt << 286
360 G4double deltaTime = << 287 G4double sx = cosTheta*cosPhi;
361 delta / << 288 G4double sy = cosTheta*sinPhi;
362 (pPreStepPoint->GetVelocity() + << 289 G4double sz = -sinTheta;
363 rand * (pPostStepPoint->GetVelocity() - << 290
364 0.5); << 291 G4ThreeVector photonPolarization(sx, sy, sz);
365 << 292
366 G4double aSecondaryTime = t0 + de << 293 // Rotate back to original coord system
367 G4ThreeVector aSecondaryPosition = x0 + ra << 294
368 << 295 photonPolarization.rotateUz(p0);
369 // Generate new G4Track object: << 296
370 G4Track* aSecondaryTrack = << 297 // Generate a new photon:
371 new G4Track(aCerenkovPhoton, aSecondaryT << 298
372 << 299 G4DynamicParticle* aCerenkovPhoton =
373 aSecondaryTrack->SetTouchableHandle( << 300 new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(),photonMomentum);
374 aStep.GetPreStepPoint()->GetTouchableHan << 301
375 aSecondaryTrack->SetParentID(aTrack.GetTra << 302 aCerenkovPhoton->SetPolarization(photonPolarization.x(),
376 aSecondaryTrack->SetCreatorModelID(secID); << 303 photonPolarization.y(),
377 aParticleChange.AddSecondary(aSecondaryTra << 304 photonPolarization.z());
378 } << 305
379 << 306 aCerenkovPhoton->SetKineticEnergy(sampledEnergy);
380 if(verboseLevel > 1) << 307
381 { << 308 // Generate new G4Track object:
382 G4cout << "\n Exiting from G4Cerenkov::DoI << 309
383 << aParticleChange.GetNumberOfSecon << 310 G4double NumberOfPhotons, N;
>> 311
>> 312 do {
>> 313 rand = G4UniformRand();
>> 314 NumberOfPhotons = MeanNumberOfPhotons1 - rand *
>> 315 (MeanNumberOfPhotons1-MeanNumberOfPhotons2);
>> 316 N = G4UniformRand() *
>> 317 std::max(MeanNumberOfPhotons1,MeanNumberOfPhotons2);
>> 318 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
>> 319 } while (N > NumberOfPhotons);
>> 320
>> 321 G4double delta = rand * aStep.GetStepLength();
>> 322
>> 323 G4double deltaTime = delta / (pPreStepPoint->GetVelocity()+
>> 324 rand*(pPostStepPoint->GetVelocity()-
>> 325 pPreStepPoint->GetVelocity())*0.5);
>> 326
>> 327 G4double aSecondaryTime = t0 + deltaTime;
>> 328
>> 329 G4ThreeVector aSecondaryPosition = x0 + rand * aStep.GetDeltaPosition();
>> 330
>> 331 G4Track* aSecondaryTrack =
>> 332 new G4Track(aCerenkovPhoton,aSecondaryTime,aSecondaryPosition);
>> 333
>> 334 aSecondaryTrack->SetTouchableHandle(
>> 335 aStep.GetPreStepPoint()->GetTouchableHandle());
>> 336
>> 337 aSecondaryTrack->SetParentID(aTrack.GetTrackID());
>> 338
>> 339 aParticleChange.AddSecondary(aSecondaryTrack);
>> 340 }
>> 341
>> 342 if (verboseLevel>0) {
>> 343 G4cout <<"\n Exiting from G4Cerenkov::DoIt -- NumberOfSecondaries = "
>> 344 << aParticleChange.GetNumberOfSecondaries() << G4endl;
384 } 345 }
385 346
386 return pParticleChange; 347 return pParticleChange;
387 } 348 }
388 349
389 //....oooOO0OOooo........oooOO0OOooo........oo << 350 // BuildThePhysicsTable for the Cerenkov process
390 void G4Cerenkov::PreparePhysicsTable(const G4P << 351 // ---------------------------------------------
>> 352 //
>> 353
>> 354 void G4Cerenkov::BuildThePhysicsTable()
391 { 355 {
392 Initialise(); << 356 if (thePhysicsTable) return;
>> 357
>> 358 const G4MaterialTable* theMaterialTable=
>> 359 G4Material::GetMaterialTable();
>> 360 G4int numOfMaterials = G4Material::GetNumberOfMaterials();
>> 361
>> 362 // create new physics table
>> 363
>> 364 thePhysicsTable = new G4PhysicsTable(numOfMaterials);
>> 365
>> 366 // loop for materials
>> 367
>> 368 for (G4int i=0 ; i < numOfMaterials; i++) {
>> 369
>> 370 G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector = 0;
>> 371
>> 372 // Retrieve vector of refraction indices for the material
>> 373 // from the material's optical properties table
>> 374
>> 375 G4Material* aMaterial = (*theMaterialTable)[i];
>> 376
>> 377 G4MaterialPropertiesTable* aMaterialPropertiesTable =
>> 378 aMaterial->GetMaterialPropertiesTable();
>> 379
>> 380 if (aMaterialPropertiesTable) {
>> 381 aPhysicsOrderedFreeVector = new G4PhysicsOrderedFreeVector();
>> 382 G4MaterialPropertyVector* theRefractionIndexVector =
>> 383 aMaterialPropertiesTable->GetProperty(kRINDEX);
>> 384
>> 385 if (theRefractionIndexVector) {
>> 386
>> 387 // Retrieve the first refraction index in vector
>> 388 // of (photon energy, refraction index) pairs
>> 389
>> 390 G4double currentRI = (*theRefractionIndexVector)[0];
>> 391
>> 392 if (currentRI > 1.0) {
>> 393
>> 394 // Create first (photon energy, Cerenkov Integral)
>> 395 // pair
>> 396
>> 397 G4double currentPM = theRefractionIndexVector->Energy(0);
>> 398 G4double currentCAI = 0.0;
>> 399
>> 400 aPhysicsOrderedFreeVector->InsertValues(currentPM , currentCAI);
>> 401
>> 402 // Set previous values to current ones prior to loop
>> 403
>> 404 G4double prevPM = currentPM;
>> 405 G4double prevCAI = currentCAI;
>> 406 G4double prevRI = currentRI;
>> 407
>> 408 // loop over all (photon energy, refraction index)
>> 409 // pairs stored for this material
>> 410
>> 411 for (size_t ii = 1;
>> 412 ii < theRefractionIndexVector->GetVectorLength();
>> 413 ++ii) {
>> 414 currentRI = (*theRefractionIndexVector)[ii];
>> 415 currentPM = theRefractionIndexVector->Energy(ii);
>> 416
>> 417 currentCAI = 0.5*(1.0/(prevRI*prevRI) +
>> 418 1.0/(currentRI*currentRI));
>> 419
>> 420 currentCAI = prevCAI + (currentPM - prevPM) * currentCAI;
>> 421
>> 422 aPhysicsOrderedFreeVector->
>> 423 InsertValues(currentPM, currentCAI);
>> 424
>> 425 prevPM = currentPM;
>> 426 prevCAI = currentCAI;
>> 427 prevRI = currentRI;
>> 428 }
>> 429
>> 430 }
>> 431 }
>> 432 }
>> 433
>> 434 // The Cerenkov integral for a given material
>> 435 // will be inserted in thePhysicsTable
>> 436 // according to the position of the material in
>> 437 // the material table.
>> 438
>> 439 thePhysicsTable->insertAt(i,aPhysicsOrderedFreeVector);
>> 440
>> 441 }
393 } 442 }
394 443
395 //....oooOO0OOooo........oooOO0OOooo........oo << 444 // GetMeanFreePath
396 G4double G4Cerenkov::GetMeanFreePath(const G4T << 445 // ---------------
397 G4ForceCo << 446 //
>> 447
>> 448 G4double G4Cerenkov::GetMeanFreePath(const G4Track&,
>> 449 G4double,
>> 450 G4ForceCondition*)
398 { 451 {
399 return 1.; 452 return 1.;
400 } 453 }
401 454
402 //....oooOO0OOooo........oooOO0OOooo........oo <<
403 G4double G4Cerenkov::PostStepGetPhysicalIntera 455 G4double G4Cerenkov::PostStepGetPhysicalInteractionLength(
404 const G4Track& aTrack, G4double, G4ForceCond << 456 const G4Track& aTrack,
>> 457 G4double,
>> 458 G4ForceCondition* condition)
405 { 459 {
406 *condition = NotForced; << 460 *condition = NotForced;
407 G4double StepLimit = DBL_MAX; 461 G4double StepLimit = DBL_MAX;
408 fNumPhotons = 0; << 462 fNumPhotons = 0;
409 463
410 const G4Material* aMaterial = aTrack.GetMate 464 const G4Material* aMaterial = aTrack.GetMaterial();
411 std::size_t materialIndex = aMaterial->Get << 465 G4int materialIndex = aMaterial->GetIndex();
412 466
413 // If Physics Vector is not defined no Ceren 467 // If Physics Vector is not defined no Cerenkov photons
414 if(!(*thePhysicsTable)[materialIndex]) << 468 // this check avoid string comparison below
415 { << 469
416 return StepLimit; << 470 if(!(*thePhysicsTable)[materialIndex]) { return StepLimit; }
417 } <<
418 471
419 const G4DynamicParticle* aParticle = aTrack. 472 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
420 const G4MaterialCutsCouple* couple = aTrack. 473 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
421 474
422 G4double kineticEnergy = a << 475 G4double kineticEnergy = aParticle->GetKineticEnergy();
423 const G4ParticleDefinition* particleType = a 476 const G4ParticleDefinition* particleType = aParticle->GetDefinition();
424 G4double mass = p << 477 G4double mass = particleType->GetPDGMass();
425 478
426 G4double beta = aParticle->GetTotalMomentum << 479 // particle beta
427 G4double gamma = aParticle->GetTotalEnergy() << 480 G4double beta = aParticle->GetTotalMomentum() /
>> 481 aParticle->GetTotalEnergy();
>> 482 // particle gamma
>> 483 G4double gamma = aParticle->GetTotalEnergy()/mass;
428 484
429 G4MaterialPropertiesTable* aMaterialProperti 485 G4MaterialPropertiesTable* aMaterialPropertiesTable =
430 aMaterial->GetMaterialPropertiesTable(); << 486 aMaterial->GetMaterialPropertiesTable();
431 487
432 G4MaterialPropertyVector* Rindex = nullptr; << 488 G4MaterialPropertyVector* Rindex = NULL;
433 489
434 if(aMaterialPropertiesTable) << 490 if (aMaterialPropertiesTable)
435 Rindex = aMaterialPropertiesTable->GetProp << 491 Rindex = aMaterialPropertiesTable->GetProperty(kRINDEX);
436 492
437 G4double nMax; 493 G4double nMax;
438 if(Rindex) << 494 if (Rindex) {
439 { << 495 nMax = Rindex->GetMaxValue();
440 nMax = Rindex->GetMaxValue(); << 496 } else {
441 } << 497 return StepLimit;
442 else << 498 }
443 { << 499
444 return StepLimit; << 500 G4double BetaMin = 1./nMax;
445 } << 501 if ( BetaMin >= 1. ) return StepLimit;
446 << 502
447 G4double BetaMin = 1. / nMax; << 503 G4double GammaMin = 1./std::sqrt(1.-BetaMin*BetaMin);
448 if(BetaMin >= 1.) << 504
449 return StepLimit; << 505 if (gamma < GammaMin ) return StepLimit;
450 << 506
451 G4double GammaMin = 1. / std::sqrt(1. - Beta << 507 G4double kinEmin = mass*(GammaMin-1.);
452 if(gamma < GammaMin) << 508
453 return StepLimit; << 509 G4double RangeMin = G4LossTableManager::Instance()->GetRange(particleType,
454 << 510 kinEmin,
455 G4double kinEmin = mass * (GammaMin - 1.); << 511 couple);
456 G4double RangeMin = << 512 G4double Range = G4LossTableManager::Instance()->GetRange(particleType,
457 G4LossTableManager::Instance()->GetRange(p << 513 kineticEnergy,
458 G4double Range = G4LossTableManager::Instanc << 514 couple);
459 particleType, kineticEnergy, couple); << 515
460 G4double Step = Range - RangeMin; 516 G4double Step = Range - RangeMin;
461 517
462 // If the step is smaller than G4ThreeVector << 518 // If the step is smaller than 1e-16 mm, it may happen that the particle
463 // that the particle does not move. See bug << 519 // does not move. See bug 1992.
464 static const G4double minAllowedStep = G4Thr << 520 // 2019-03-11: change to 1e-15
465 if(Step < minAllowedStep) << 521 if (Step < 1.e-15*mm) return StepLimit;
466 return StepLimit; << 522
467 << 523 if (Step < StepLimit) StepLimit = Step;
468 if(Step < StepLimit) << 524 // If user has defined an average maximum number of photons to
469 StepLimit = Step; << 525 // be generated in a Step, then calculate the Step length for
470 << 526 // that number of photons.
471 // If user has defined an average maximum nu << 527
472 // a Step, then calculate the Step length fo << 528 if (fMaxPhotons > 0) {
473 if(fMaxPhotons > 0) << 529
474 { << 530 // particle charge
475 const G4double charge = aParticle->GetDefi << 531 const G4double charge = aParticle->GetDefinition()->GetPDGCharge();
476 G4double MeanNumberOfPhotons = << 532
477 GetAverageNumberOfPhotons(charge, beta, << 533 G4double MeanNumberOfPhotons =
478 Step = 0.; << 534 GetAverageNumberOfPhotons(charge,beta,aMaterial,Rindex);
479 if(MeanNumberOfPhotons > 0.0) << 535
480 Step = fMaxPhotons / MeanNumberOfPhotons << 536 Step = 0.;
481 if(Step > 0. && Step < StepLimit) << 537 if (MeanNumberOfPhotons > 0.0) Step = fMaxPhotons / MeanNumberOfPhotons;
482 StepLimit = Step; << 538
>> 539 if (Step > 0. && Step < StepLimit) StepLimit = Step;
483 } 540 }
484 541
485 // If user has defined an maximum allowed ch 542 // If user has defined an maximum allowed change in beta per step
486 if(fMaxBetaChange > 0.) << 543 if (fMaxBetaChange > 0.) {
487 { << 544
488 G4double dedx = G4LossTableManager::Instan << 545 G4double dedx = G4LossTableManager::Instance()->GetDEDX(particleType,
489 particleType, kineticEnergy, couple); << 546 kineticEnergy,
490 G4double deltaGamma = << 547 couple);
491 gamma - 1. / std::sqrt(1. - beta * beta << 548
492 (1. - fMax << 549 G4double deltaGamma = gamma - 1./std::sqrt(1.-beta*beta*
493 << 550 (1.-fMaxBetaChange)*
494 Step = mass * deltaGamma / dedx; << 551 (1.-fMaxBetaChange));
495 if(Step > 0. && Step < StepLimit) << 552
496 StepLimit = Step; << 553 Step = mass * deltaGamma / dedx;
>> 554
>> 555 if (Step > 0. && Step < StepLimit) StepLimit = Step;
>> 556
497 } 557 }
498 558
499 *condition = StronglyForced; 559 *condition = StronglyForced;
500 return StepLimit; 560 return StepLimit;
501 } 561 }
502 562
503 //....oooOO0OOooo........oooOO0OOooo........oo << 563 // GetAverageNumberOfPhotons
504 G4double G4Cerenkov::GetAverageNumberOfPhotons << 564 // -------------------------
505 const G4double charge, const G4double beta, <<
506 G4MaterialPropertyVector* Rindex) const <<
507 // This routine computes the number of Cerenko 565 // This routine computes the number of Cerenkov photons produced per
508 // Geant4-unit (millimeter) in the current med << 566 // GEANT-unit (millimeter) in the current medium.
>> 567 // ^^^^^^^^^^
>> 568
>> 569 G4double
>> 570 G4Cerenkov::GetAverageNumberOfPhotons(const G4double charge,
>> 571 const G4double beta,
>> 572 const G4Material* aMaterial,
>> 573 G4MaterialPropertyVector* Rindex) const
509 { 574 {
510 constexpr G4double Rfact = 369.81 / (eV * cm << 575 const G4double Rfact = 369.81/(eV * cm);
511 if(beta <= 0.0) << 576
512 return 0.0; << 577 if(beta <= 0.0)return 0.0;
513 G4double BetaInverse = 1. / beta; << 578
>> 579 G4double BetaInverse = 1./beta;
514 580
515 // Vectors used in computation of Cerenkov A 581 // Vectors used in computation of Cerenkov Angle Integral:
516 // - Refraction Indices for the current mat 582 // - Refraction Indices for the current material
517 // - new G4PhysicsFreeVector allocated to h << 583 // - new G4PhysicsOrderedFreeVector allocated to hold CAI's
518 std::size_t materialIndex = aMaterial->GetIn << 584
>> 585 G4int materialIndex = aMaterial->GetIndex();
>> 586
>> 587 // Retrieve the Cerenkov Angle Integrals for this material
>> 588
>> 589 G4PhysicsOrderedFreeVector* CerenkovAngleIntegrals =
>> 590 (G4PhysicsOrderedFreeVector*)((*thePhysicsTable)(materialIndex));
519 591
520 // Retrieve the Cerenkov Angle Integrals for << 592 if(!(CerenkovAngleIntegrals->IsFilledVectorExist()))return 0.0;
521 G4PhysicsVector* CerenkovAngleIntegrals = (( <<
522 593
523 std::size_t length = CerenkovAngleIntegrals- << 594 // Min and Max photon energies
524 if(0 == length) << 595 G4double Pmin = Rindex->GetMinLowEdgeEnergy();
525 return 0.0; << 596 G4double Pmax = Rindex->GetMaxLowEdgeEnergy();
526 <<
527 // Min and Max photon energies <<
528 G4double Pmin = Rindex->Energy(0); <<
529 G4double Pmax = Rindex->GetMaxEnergy(); <<
530 597
531 // Min and Max Refraction Indices << 598 // Min and Max Refraction Indices
532 G4double nMin = Rindex->GetMinValue(); << 599 G4double nMin = Rindex->GetMinValue();
533 G4double nMax = Rindex->GetMaxValue(); 600 G4double nMax = Rindex->GetMaxValue();
534 601
535 // Max Cerenkov Angle Integral << 602 // Max Cerenkov Angle Integral
536 G4double CAImax = (*CerenkovAngleIntegrals)[ << 603 G4double CAImax = CerenkovAngleIntegrals->GetMaxValue();
537 604
538 G4double dp, ge; 605 G4double dp, ge;
539 // If n(Pmax) < 1/Beta -- no photons generat <<
540 if(nMax < BetaInverse) <<
541 { <<
542 dp = 0.0; <<
543 ge = 0.0; <<
544 } <<
545 // otherwise if n(Pmin) >= 1/Beta -- photons <<
546 else if(nMin > BetaInverse) <<
547 { <<
548 dp = Pmax - Pmin; <<
549 ge = CAImax; <<
550 } <<
551 // If n(Pmin) < 1/Beta, and n(Pmax) >= 1/Bet <<
552 // that the value of n(P) == 1/Beta. Interpo <<
553 // GetEnergy() and Value() methods of the G4 <<
554 // the Value() method of G4PhysicsVector. <<
555 else <<
556 { <<
557 Pmin = Rindex->GetEnergy(BetaInverse); <<
558 dp = Pmax - Pmin; <<
559 <<
560 G4double CAImin = CerenkovAngleIntegrals-> <<
561 ge = CAImax - CAImin; <<
562 <<
563 if(verboseLevel > 1) <<
564 { <<
565 G4cout << "CAImin = " << CAImin << G4end <<
566 } <<
567 } <<
568 <<
569 // Calculate number of photons <<
570 G4double NumPhotons = Rfact * charge / eplus <<
571 (dp - ge * BetaInverse <<
572 606
573 return NumPhotons; << 607 // If n(Pmax) < 1/Beta -- no photons generated
574 } <<
575 608
576 //....oooOO0OOooo........oooOO0OOooo........oo << 609 if (nMax < BetaInverse) {
577 void G4Cerenkov::SetTrackSecondariesFirst(cons << 610 dp = 0.0;
578 { << 611 ge = 0.0;
579 fTrackSecondariesFirst = state; << 612 }
580 G4OpticalParameters::Instance()->SetCerenkov << 613
581 fTrackSecondariesFirst); << 614 // otherwise if n(Pmin) >= 1/Beta -- photons generated
582 } << 615
>> 616 else if (nMin > BetaInverse) {
>> 617 dp = Pmax - Pmin;
>> 618 ge = CAImax;
>> 619 }
>> 620
>> 621 // If n(Pmin) < 1/Beta, and n(Pmax) >= 1/Beta, then
>> 622 // we need to find a P such that the value of n(P) == 1/Beta.
>> 623 // Interpolation is performed by the GetEnergy() and
>> 624 // Value() methods of the G4MaterialPropertiesTable and
>> 625 // the GetValue() method of G4PhysicsVector.
>> 626
>> 627 else {
>> 628 Pmin = Rindex->GetEnergy(BetaInverse);
>> 629 dp = Pmax - Pmin;
>> 630
>> 631 // need boolean for current implementation of G4PhysicsVector
>> 632 // ==> being phased out
>> 633 G4bool isOutRange;
>> 634 G4double CAImin = CerenkovAngleIntegrals->GetValue(Pmin, isOutRange);
>> 635 ge = CAImax - CAImin;
>> 636
>> 637 if (verboseLevel>0) {
>> 638 G4cout << "CAImin = " << CAImin << G4endl;
>> 639 G4cout << "ge = " << ge << G4endl;
>> 640 }
>> 641 }
>> 642
>> 643 // Calculate number of photons
>> 644 G4double NumPhotons = Rfact * charge/eplus * charge/eplus *
>> 645 (dp - ge * BetaInverse*BetaInverse);
583 646
584 //....oooOO0OOooo........oooOO0OOooo........oo << 647 return NumPhotons;
585 void G4Cerenkov::SetMaxBetaChangePerStep(const <<
586 { <<
587 fMaxBetaChange = value * CLHEP::perCent; <<
588 G4OpticalParameters::Instance()->SetCerenkov <<
589 } 648 }
590 649
591 //....oooOO0OOooo........oooOO0OOooo........oo <<
592 void G4Cerenkov::SetMaxNumPhotonsPerStep(const <<
593 { <<
594 fMaxPhotons = NumPhotons; <<
595 G4OpticalParameters::Instance()->SetCerenkov <<
596 } <<
597 <<
598 void G4Cerenkov::SetStackPhotons(const G4bool <<
599 { <<
600 fStackingFlag = stackingFlag; <<
601 G4OpticalParameters::Instance()->SetCerenkov <<
602 } <<
603 <<
604 //....oooOO0OOooo........oooOO0OOooo........oo <<
605 void G4Cerenkov::DumpPhysicsTable() const 650 void G4Cerenkov::DumpPhysicsTable() const
606 { 651 {
607 G4cout << "Dump Physics Table!" << G4endl; << 652 G4int PhysicsTableSize = thePhysicsTable->entries();
608 for(std::size_t i = 0; i < thePhysicsTable-> << 653 G4PhysicsOrderedFreeVector *v;
609 { << 654
610 (*thePhysicsTable)[i]->DumpValues(); << 655 for (G4int i = 0 ; i < PhysicsTableSize ; i++ ) {
>> 656 v = (G4PhysicsOrderedFreeVector*)(*thePhysicsTable)[i];
>> 657 v->DumpValues();
611 } 658 }
612 } 659 }
613 660
614 //....oooOO0OOooo........oooOO0OOooo........oo <<
615 void G4Cerenkov::SetVerboseLevel(G4int verbose <<
616 { <<
617 verboseLevel = verbose; <<
618 G4OpticalParameters::Instance()->SetCerenkov <<
619 } <<
620 661