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