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