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 92046 2015-08-14 07:22:07Z 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 G4double charge = aParticle->GetDefinition()->GetPDGCharge(); >> 206 >> 207 // particle beta >> 208 G4double beta = (pPreStepPoint ->GetBeta() + >> 209 pPostStepPoint->GetBeta())*0.5; >> 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 G4double beta1 = pPreStepPoint ->GetBeta(); >> 263 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 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko >> 290 } while (rand*maxSin2 > sin2Theta); >> 291 >> 292 // Generate random position of photon on cone surface >> 293 // defined by Theta >> 294 >> 295 rand = G4UniformRand(); >> 296 >> 297 G4double phi = twopi*rand; >> 298 G4double sinPhi = std::sin(phi); >> 299 G4double cosPhi = std::cos(phi); >> 300 >> 301 // calculate x,y, and z components of photon energy >> 302 // (in coord system with primary particle direction >> 303 // aligned with the z axis) >> 304 >> 305 G4double sinTheta = std::sqrt(sin2Theta); >> 306 G4double px = sinTheta*cosPhi; >> 307 G4double py = sinTheta*sinPhi; >> 308 G4double pz = cosTheta; >> 309 >> 310 // Create photon momentum direction vector >> 311 // The momentum direction is still with respect >> 312 // to the coordinate system where the primary >> 313 // particle direction is aligned with the z axis >> 314 >> 315 G4ParticleMomentum photonMomentum(px, py, pz); >> 316 >> 317 // Rotate momentum direction back to global reference >> 318 // system >> 319 >> 320 photonMomentum.rotateUz(p0); >> 321 >> 322 // Determine polarization of new photon >> 323 >> 324 G4double sx = cosTheta*cosPhi; >> 325 G4double sy = cosTheta*sinPhi; >> 326 G4double sz = -sinTheta; >> 327 >> 328 G4ThreeVector photonPolarization(sx, sy, sz); >> 329 >> 330 // Rotate back to original coord system >> 331 >> 332 photonPolarization.rotateUz(p0); >> 333 >> 334 // Generate a new photon: >> 335 >> 336 G4DynamicParticle* aCerenkovPhoton = >> 337 new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(), >> 338 photonMomentum); >> 339 aCerenkovPhoton->SetPolarization >> 340 (photonPolarization.x(), >> 341 photonPolarization.y(), >> 342 photonPolarization.z()); >> 343 >> 344 aCerenkovPhoton->SetKineticEnergy(sampledEnergy); >> 345 >> 346 // Generate new G4Track object: >> 347 >> 348 G4double NumberOfPhotons, N; >> 349 >> 350 do { >> 351 rand = G4UniformRand(); >> 352 NumberOfPhotons = MeanNumberOfPhotons1 - rand * >> 353 (MeanNumberOfPhotons1-MeanNumberOfPhotons2); >> 354 N = G4UniformRand() * >> 355 std::max(MeanNumberOfPhotons1,MeanNumberOfPhotons2); >> 356 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko >> 357 } while (N > NumberOfPhotons); >> 358 >> 359 G4double delta = rand * aStep.GetStepLength(); >> 360 >> 361 G4double deltaTime = delta / (pPreStepPoint->GetVelocity()+ >> 362 rand*(pPostStepPoint->GetVelocity()- >> 363 pPreStepPoint->GetVelocity())*0.5); >> 364 >> 365 G4double aSecondaryTime = t0 + deltaTime; >> 366 >> 367 G4ThreeVector aSecondaryPosition = >> 368 x0 + rand * aStep.GetDeltaPosition(); >> 369 >> 370 G4Track* aSecondaryTrack = >> 371 new G4Track(aCerenkovPhoton,aSecondaryTime,aSecondaryPosition); >> 372 >> 373 aSecondaryTrack->SetTouchableHandle( >> 374 aStep.GetPreStepPoint()->GetTouchableHandle()); >> 375 >> 376 aSecondaryTrack->SetParentID(aTrack.GetTrackID()); >> 377 >> 378 aParticleChange.AddSecondary(aSecondaryTrack); >> 379 } >> 380 >> 381 if (verboseLevel>0) { >> 382 G4cout <<"\n Exiting from G4Cerenkov::DoIt -- NumberOfSecondaries = " >> 383 << aParticleChange.GetNumberOfSecondaries() << G4endl; >> 384 } >> 385 >> 386 return pParticleChange; 589 } 387 } 590 388 591 //....oooOO0OOooo........oooOO0OOooo........oo << 389 // BuildThePhysicsTable for the Cerenkov process 592 void G4Cerenkov::SetMaxNumPhotonsPerStep(const << 390 // --------------------------------------------- >> 391 // >> 392 >> 393 void G4Cerenkov::BuildThePhysicsTable() 593 { 394 { 594 fMaxPhotons = NumPhotons; << 395 if (thePhysicsTable) return; 595 G4OpticalParameters::Instance()->SetCerenkov << 396 >> 397 const G4MaterialTable* theMaterialTable= >> 398 G4Material::GetMaterialTable(); >> 399 G4int numOfMaterials = G4Material::GetNumberOfMaterials(); >> 400 >> 401 // create new physics table >> 402 >> 403 thePhysicsTable = new G4PhysicsTable(numOfMaterials); >> 404 >> 405 // loop for materials >> 406 >> 407 for (G4int i=0 ; i < numOfMaterials; i++) >> 408 { >> 409 G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector = 0; >> 410 >> 411 // Retrieve vector of refraction indices for the material >> 412 // from the material's optical properties table >> 413 >> 414 G4Material* aMaterial = (*theMaterialTable)[i]; >> 415 >> 416 G4MaterialPropertiesTable* aMaterialPropertiesTable = >> 417 aMaterial->GetMaterialPropertiesTable(); >> 418 >> 419 if (aMaterialPropertiesTable) { >> 420 >> 421 aPhysicsOrderedFreeVector = new G4PhysicsOrderedFreeVector(); >> 422 G4MaterialPropertyVector* theRefractionIndexVector = >> 423 aMaterialPropertiesTable->GetProperty("RINDEX"); >> 424 >> 425 if (theRefractionIndexVector) { >> 426 >> 427 // Retrieve the first refraction index in vector >> 428 // of (photon energy, refraction index) pairs >> 429 >> 430 G4double currentRI = (*theRefractionIndexVector)[0]; >> 431 >> 432 if (currentRI > 1.0) { >> 433 >> 434 // Create first (photon energy, Cerenkov Integral) >> 435 // pair >> 436 >> 437 G4double currentPM = theRefractionIndexVector-> >> 438 Energy(0); >> 439 G4double currentCAI = 0.0; >> 440 >> 441 aPhysicsOrderedFreeVector-> >> 442 InsertValues(currentPM , currentCAI); >> 443 >> 444 // Set previous values to current ones prior to loop >> 445 >> 446 G4double prevPM = currentPM; >> 447 G4double prevCAI = currentCAI; >> 448 G4double prevRI = currentRI; >> 449 >> 450 // loop over all (photon energy, refraction index) >> 451 // pairs stored for this material >> 452 >> 453 for (size_t ii = 1; >> 454 ii < theRefractionIndexVector->GetVectorLength(); >> 455 ++ii) >> 456 { >> 457 currentRI = (*theRefractionIndexVector)[ii]; >> 458 currentPM = theRefractionIndexVector->Energy(ii); >> 459 >> 460 currentCAI = 0.5*(1.0/(prevRI*prevRI) + >> 461 1.0/(currentRI*currentRI)); >> 462 >> 463 currentCAI = prevCAI + >> 464 (currentPM - prevPM) * currentCAI; >> 465 >> 466 aPhysicsOrderedFreeVector-> >> 467 InsertValues(currentPM, currentCAI); >> 468 >> 469 prevPM = currentPM; >> 470 prevCAI = currentCAI; >> 471 prevRI = currentRI; >> 472 } >> 473 >> 474 } >> 475 } >> 476 } >> 477 >> 478 // The Cerenkov integral for a given material >> 479 // will be inserted in thePhysicsTable >> 480 // according to the position of the material in >> 481 // the material table. >> 482 >> 483 thePhysicsTable->insertAt(i,aPhysicsOrderedFreeVector); >> 484 >> 485 } 596 } 486 } 597 487 598 void G4Cerenkov::SetStackPhotons(const G4bool << 488 // GetMeanFreePath >> 489 // --------------- >> 490 // >> 491 >> 492 G4double G4Cerenkov::GetMeanFreePath(const G4Track&, >> 493 G4double, >> 494 G4ForceCondition*) 599 { 495 { 600 fStackingFlag = stackingFlag; << 496 return 1.; 601 G4OpticalParameters::Instance()->SetCerenkov << 602 } 497 } 603 498 604 //....oooOO0OOooo........oooOO0OOooo........oo << 499 G4double G4Cerenkov::PostStepGetPhysicalInteractionLength( 605 void G4Cerenkov::DumpPhysicsTable() const << 500 const G4Track& aTrack, 606 { << 501 G4double, 607 G4cout << "Dump Physics Table!" << G4endl; << 502 G4ForceCondition* condition) 608 for(std::size_t i = 0; i < thePhysicsTable-> << 503 { 609 { << 504 *condition = NotForced; 610 (*thePhysicsTable)[i]->DumpValues(); << 505 G4double StepLimit = DBL_MAX; 611 } << 506 >> 507 const G4Material* aMaterial = aTrack.GetMaterial(); >> 508 G4int materialIndex = aMaterial->GetIndex(); >> 509 >> 510 // If Physics Vector is not defined no Cerenkov photons >> 511 // this check avoid string comparison below >> 512 if(!(*thePhysicsTable)[materialIndex]) { return StepLimit; } >> 513 >> 514 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); >> 515 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple(); >> 516 >> 517 G4double kineticEnergy = aParticle->GetKineticEnergy(); >> 518 const G4ParticleDefinition* particleType = aParticle->GetDefinition(); >> 519 G4double mass = particleType->GetPDGMass(); >> 520 >> 521 // particle beta >> 522 G4double beta = aParticle->GetTotalMomentum() / >> 523 aParticle->GetTotalEnergy(); >> 524 // particle gamma >> 525 G4double gamma = aParticle->GetTotalEnergy()/mass; >> 526 >> 527 G4MaterialPropertiesTable* aMaterialPropertiesTable = >> 528 aMaterial->GetMaterialPropertiesTable(); >> 529 >> 530 G4MaterialPropertyVector* Rindex = NULL; >> 531 >> 532 if (aMaterialPropertiesTable) >> 533 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX"); >> 534 >> 535 G4double nMax; >> 536 if (Rindex) { >> 537 nMax = Rindex->GetMaxValue(); >> 538 } else { >> 539 return StepLimit; >> 540 } >> 541 >> 542 G4double BetaMin = 1./nMax; >> 543 if ( BetaMin >= 1. ) return StepLimit; >> 544 >> 545 G4double GammaMin = 1./std::sqrt(1.-BetaMin*BetaMin); >> 546 >> 547 if (gamma < GammaMin ) return StepLimit; >> 548 >> 549 G4double kinEmin = mass*(GammaMin-1.); >> 550 >> 551 G4double RangeMin = G4LossTableManager::Instance()-> >> 552 GetRange(particleType, >> 553 kinEmin, >> 554 couple); >> 555 G4double Range = G4LossTableManager::Instance()-> >> 556 GetRange(particleType, >> 557 kineticEnergy, >> 558 couple); >> 559 >> 560 G4double Step = Range - RangeMin; >> 561 if (Step < 1.*um ) return StepLimit; >> 562 >> 563 if (Step > 0. && Step < StepLimit) StepLimit = Step; >> 564 >> 565 // If user has defined an average maximum number of photons to >> 566 // be generated in a Step, then calculate the Step length for >> 567 // that number of photons. >> 568 >> 569 if (fMaxPhotons > 0) { >> 570 >> 571 // particle charge >> 572 const G4double charge = aParticle-> >> 573 GetDefinition()->GetPDGCharge(); >> 574 >> 575 G4double MeanNumberOfPhotons = >> 576 GetAverageNumberOfPhotons(charge,beta,aMaterial,Rindex); >> 577 >> 578 Step = 0.; >> 579 if (MeanNumberOfPhotons > 0.0) Step = fMaxPhotons / >> 580 MeanNumberOfPhotons; >> 581 >> 582 if (Step > 0. && Step < StepLimit) StepLimit = Step; >> 583 } >> 584 >> 585 // If user has defined an maximum allowed change in beta per step >> 586 if (fMaxBetaChange > 0.) { >> 587 >> 588 G4double dedx = G4LossTableManager::Instance()-> >> 589 GetDEDX(particleType, >> 590 kineticEnergy, >> 591 couple); >> 592 >> 593 G4double deltaGamma = gamma - >> 594 1./std::sqrt(1.-beta*beta* >> 595 (1.-fMaxBetaChange)* >> 596 (1.-fMaxBetaChange)); >> 597 >> 598 Step = mass * deltaGamma / dedx; >> 599 >> 600 if (Step > 0. && Step < StepLimit) StepLimit = Step; >> 601 >> 602 } >> 603 >> 604 *condition = StronglyForced; >> 605 return StepLimit; 612 } 606 } 613 607 614 //....oooOO0OOooo........oooOO0OOooo........oo << 608 // GetAverageNumberOfPhotons 615 void G4Cerenkov::SetVerboseLevel(G4int verbose << 609 // ------------------------- >> 610 // This routine computes the number of Cerenkov photons produced per >> 611 // GEANT-unit (millimeter) in the current medium. >> 612 // ^^^^^^^^^^ >> 613 >> 614 G4double >> 615 G4Cerenkov::GetAverageNumberOfPhotons(const G4double charge, >> 616 const G4double beta, >> 617 const G4Material* aMaterial, >> 618 G4MaterialPropertyVector* Rindex) const 616 { 619 { 617 verboseLevel = verbose; << 620 const G4double Rfact = 369.81/(eV * cm); 618 G4OpticalParameters::Instance()->SetCerenkov << 621 >> 622 if(beta <= 0.0)return 0.0; >> 623 >> 624 G4double BetaInverse = 1./beta; >> 625 >> 626 // Vectors used in computation of Cerenkov Angle Integral: >> 627 // - Refraction Indices for the current material >> 628 // - new G4PhysicsOrderedFreeVector allocated to hold CAI's >> 629 >> 630 G4int materialIndex = aMaterial->GetIndex(); >> 631 >> 632 // Retrieve the Cerenkov Angle Integrals for this material >> 633 >> 634 G4PhysicsOrderedFreeVector* CerenkovAngleIntegrals = >> 635 (G4PhysicsOrderedFreeVector*)((*thePhysicsTable)(materialIndex)); >> 636 >> 637 if(!(CerenkovAngleIntegrals->IsFilledVectorExist()))return 0.0; >> 638 >> 639 // Min and Max photon energies >> 640 G4double Pmin = Rindex->GetMinLowEdgeEnergy(); >> 641 G4double Pmax = Rindex->GetMaxLowEdgeEnergy(); >> 642 >> 643 // Min and Max Refraction Indices >> 644 G4double nMin = Rindex->GetMinValue(); >> 645 G4double nMax = Rindex->GetMaxValue(); >> 646 >> 647 // Max Cerenkov Angle Integral >> 648 G4double CAImax = CerenkovAngleIntegrals->GetMaxValue(); >> 649 >> 650 G4double dp, ge; >> 651 >> 652 // If n(Pmax) < 1/Beta -- no photons generated >> 653 >> 654 if (nMax < BetaInverse) { >> 655 dp = 0; >> 656 ge = 0; >> 657 } >> 658 >> 659 // otherwise if n(Pmin) >= 1/Beta -- photons generated >> 660 >> 661 else if (nMin > BetaInverse) { >> 662 dp = Pmax - Pmin; >> 663 ge = CAImax; >> 664 } >> 665 >> 666 // If n(Pmin) < 1/Beta, and n(Pmax) >= 1/Beta, then >> 667 // we need to find a P such that the value of n(P) == 1/Beta. >> 668 // Interpolation is performed by the GetEnergy() and >> 669 // Value() methods of the G4MaterialPropertiesTable and >> 670 // the GetValue() method of G4PhysicsVector. >> 671 >> 672 else { >> 673 Pmin = Rindex->GetEnergy(BetaInverse); >> 674 dp = Pmax - Pmin; >> 675 >> 676 // need boolean for current implementation of G4PhysicsVector >> 677 // ==> being phased out >> 678 G4bool isOutRange; >> 679 G4double CAImin = CerenkovAngleIntegrals-> >> 680 GetValue(Pmin, isOutRange); >> 681 ge = CAImax - CAImin; >> 682 >> 683 if (verboseLevel>0) { >> 684 G4cout << "CAImin = " << CAImin << G4endl; >> 685 G4cout << "ge = " << ge << G4endl; >> 686 } >> 687 } >> 688 >> 689 // Calculate number of photons >> 690 G4double NumPhotons = Rfact * charge/eplus * charge/eplus * >> 691 (dp - ge * BetaInverse*BetaInverse); >> 692 >> 693 return NumPhotons; 619 } 694 } 620 695