ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Cerenkov.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Cerenkov.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 //
29 // Cerenkov Radiation Class Implementation
31 //
32 // File: G4Cerenkov.cc
33 // Description: Discrete Process -- Generation of Cerenkov Photons
34 // Version: 2.1
35 // Created: 1996-02-21
36 // Author: Juliet Armstrong
37 // Updated: 2007-09-30 by Peter Gumplinger
38 // > change inheritance to G4VDiscreteProcess
39 // GetContinuousStepLimit -> GetMeanFreePath (StronglyForced)
40 // AlongStepDoIt -> PostStepDoIt
41 // 2005-08-17 by Peter Gumplinger
42 // > change variable name MeanNumPhotons -> MeanNumberOfPhotons
43 // 2005-07-28 by Peter Gumplinger
44 // > add G4ProcessType to constructor
45 // 2001-09-17, migration of Materials to pure STL (mma)
46 // 2000-11-12 by Peter Gumplinger
47 // > add check on CerenkovAngleIntegrals->IsFilledVectorExist()
48 // in method GetAverageNumberOfPhotons
49 // > and a test for MeanNumberOfPhotons <= 0.0 in DoIt
50 // 2000-09-18 by Peter Gumplinger
51 // > change: aSecondaryPosition=x0+rand*aStep.GetDeltaPosition();
52 // aSecondaryTrack->SetTouchable(0);
53 // 1999-10-29 by Peter Gumplinger
54 // > change: == into <= in GetContinuousStepLimit
55 // 1997-08-08 by Peter Gumplinger
56 // > add protection against /0
57 // > G4MaterialPropertiesTable; new physics/tracking scheme
58 //
59 // mail: gum@triumf.ca
60 //
62 
63 #include "G4ios.hh"
64 #include "G4PhysicalConstants.hh"
65 #include "G4SystemOfUnits.hh"
66 #include "G4Poisson.hh"
67 #include "G4EmProcessSubType.hh"
68 
69 #include "G4LossTableManager.hh"
70 #include "G4MaterialCutsCouple.hh"
71 #include "G4ParticleDefinition.hh"
72 
73 #include "G4Cerenkov.hh"
74 
76  : G4VProcess(processName, type),
77  fTrackSecondariesFirst(false),
78  fMaxBetaChange(0.0),
79  fMaxPhotons(0),
80  fStackingFlag(true),
81  fNumPhotons(0)
82 {
84 
85  thePhysicsTable = nullptr;
86 
87  if (verboseLevel>0) {
88  G4cout << GetProcessName() << " is created " << G4endl;
89  }
90 }
91 
93 {
94  if (thePhysicsTable != nullptr) {
96  delete thePhysicsTable;
97  }
98 }
99 
101 {
102  return (aParticleType.GetPDGCharge() != 0.0 &&
103  aParticleType.GetPDGMass() != 0.0 &&
104  aParticleType.GetParticleName() != "chargedgeantino" &&
105  !aParticleType.IsShortLived() ) ? true : false;
106 }
107 
109 {
110  fTrackSecondariesFirst = state;
111 }
112 
114 {
116 }
117 
119 {
120  fMaxPhotons = NumPhotons;
121 }
122 
124 {
126 }
127 
128 // PostStepDoIt
129 // -------------
130 //
132 G4Cerenkov::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
133 
134 // This routine is called for each tracking Step of a charged particle
135 // in a radiator. A Poisson-distributed number of photons is generated
136 // according to the Cerenkov formula, distributed evenly along the track
137 // segment and uniformly azimuth w.r.t. the particle direction. The
138 // parameters are then transformed into the Master Reference System, and
139 // they are added to the particle change.
140 
141 {
143  // Should we ensure that the material is dispersive?
145 
146  aParticleChange.Initialize(aTrack);
147 
148  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
149  const G4Material* aMaterial = aTrack.GetMaterial();
150 
151  G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
152  G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
153 
154  G4ThreeVector x0 = pPreStepPoint->GetPosition();
155  G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
156  G4double t0 = pPreStepPoint->GetGlobalTime();
157 
158  G4MaterialPropertiesTable* aMaterialPropertiesTable =
159  aMaterial->GetMaterialPropertiesTable();
160  if (!aMaterialPropertiesTable) return pParticleChange;
161 
162  G4MaterialPropertyVector* Rindex =
163  aMaterialPropertiesTable->GetProperty(kRINDEX);
164  if (!Rindex) return pParticleChange;
165 
166  // particle charge
167  G4double charge = aParticle->GetDefinition()->GetPDGCharge();
168 
169  // particle beta
170  G4double beta = (pPreStepPoint->GetBeta() + pPostStepPoint->GetBeta())*0.5;
171 
172  //fNumPhotons = 0; // in PostStepGetPhysicalInteractionLength()
173 
174  G4double MeanNumberOfPhotons =
175  GetAverageNumberOfPhotons(charge,beta,aMaterial,Rindex);
176 
177  if (MeanNumberOfPhotons <= 0.0) {
178 
179  // return unchanged particle and no secondaries
180 
182 
183  return pParticleChange;
184 
185  }
186 
187  G4double step_length = aStep.GetStepLength();
188 
189  MeanNumberOfPhotons = MeanNumberOfPhotons * step_length;
190 
191  fNumPhotons = (G4int) G4Poisson(MeanNumberOfPhotons);
192 
193  if ( fNumPhotons <= 0 || !fStackingFlag ) {
194 
195  // return unchanged particle and no secondaries
196 
198 
199  return pParticleChange;
200 
201  }
202 
204 
206 
208  if (aTrack.GetTrackStatus() == fAlive )
210  }
211 
213 
214  G4double Pmin = Rindex->GetMinLowEdgeEnergy();
215  G4double Pmax = Rindex->GetMaxLowEdgeEnergy();
216  G4double dp = Pmax - Pmin;
217 
218  G4double nMax = Rindex->GetMaxValue();
219 
220  G4double BetaInverse = 1./beta;
221 
222  G4double maxCos = BetaInverse / nMax;
223  G4double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos);
224 
225  G4double beta1 = pPreStepPoint ->GetBeta();
226  G4double beta2 = pPostStepPoint->GetBeta();
227 
228  G4double MeanNumberOfPhotons1 =
229  GetAverageNumberOfPhotons(charge,beta1,aMaterial,Rindex);
230  G4double MeanNumberOfPhotons2 =
231  GetAverageNumberOfPhotons(charge,beta2,aMaterial,Rindex);
232 
233  for (G4int i = 0; i < fNumPhotons; i++) {
234 
235  // Determine photon energy
236 
237  G4double rand;
238  G4double sampledEnergy, sampledRI;
239  G4double cosTheta, sin2Theta;
240 
241  // sample an energy
242 
243  do {
244  rand = G4UniformRand();
245  sampledEnergy = Pmin + rand * dp;
246  sampledRI = Rindex->Value(sampledEnergy);
247  cosTheta = BetaInverse / sampledRI;
248 
249  sin2Theta = (1.0 - cosTheta)*(1.0 + cosTheta);
250  rand = G4UniformRand();
251 
252  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
253  } while (rand*maxSin2 > sin2Theta);
254 
255  // Generate random position of photon on cone surface
256  // defined by Theta
257 
258  rand = G4UniformRand();
259 
260  G4double phi = twopi*rand;
261  G4double sinPhi = std::sin(phi);
262  G4double cosPhi = std::cos(phi);
263 
264  // calculate x,y, and z components of photon energy
265  // (in coord system with primary particle direction
266  // aligned with the z axis)
267 
268  G4double sinTheta = std::sqrt(sin2Theta);
269  G4double px = sinTheta*cosPhi;
270  G4double py = sinTheta*sinPhi;
271  G4double pz = cosTheta;
272 
273  // Create photon momentum direction vector
274  // The momentum direction is still with respect
275  // to the coordinate system where the primary
276  // particle direction is aligned with the z axis
277 
278  G4ParticleMomentum photonMomentum(px, py, pz);
279 
280  // Rotate momentum direction back to global reference
281  // system
282 
283  photonMomentum.rotateUz(p0);
284 
285  // Determine polarization of new photon
286 
287  G4double sx = cosTheta*cosPhi;
288  G4double sy = cosTheta*sinPhi;
289  G4double sz = -sinTheta;
290 
291  G4ThreeVector photonPolarization(sx, sy, sz);
292 
293  // Rotate back to original coord system
294 
295  photonPolarization.rotateUz(p0);
296 
297  // Generate a new photon:
298 
299  G4DynamicParticle* aCerenkovPhoton =
301 
302  aCerenkovPhoton->SetPolarization(photonPolarization.x(),
303  photonPolarization.y(),
304  photonPolarization.z());
305 
306  aCerenkovPhoton->SetKineticEnergy(sampledEnergy);
307 
308  // Generate new G4Track object:
309 
310  G4double NumberOfPhotons, N;
311 
312  do {
313  rand = G4UniformRand();
314  NumberOfPhotons = MeanNumberOfPhotons1 - rand *
315  (MeanNumberOfPhotons1-MeanNumberOfPhotons2);
316  N = G4UniformRand() *
317  std::max(MeanNumberOfPhotons1,MeanNumberOfPhotons2);
318  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
319  } while (N > NumberOfPhotons);
320 
321  G4double delta = rand * aStep.GetStepLength();
322 
323  G4double deltaTime = delta / (pPreStepPoint->GetVelocity()+
324  rand*(pPostStepPoint->GetVelocity()-
325  pPreStepPoint->GetVelocity())*0.5);
326 
327  G4double aSecondaryTime = t0 + deltaTime;
328 
329  G4ThreeVector aSecondaryPosition = x0 + rand * aStep.GetDeltaPosition();
330 
331  G4Track* aSecondaryTrack =
332  new G4Track(aCerenkovPhoton,aSecondaryTime,aSecondaryPosition);
333 
334  aSecondaryTrack->SetTouchableHandle(
336 
337  aSecondaryTrack->SetParentID(aTrack.GetTrackID());
338 
339  aParticleChange.AddSecondary(aSecondaryTrack);
340  }
341 
342  if (verboseLevel>0) {
343  G4cout <<"\n Exiting from G4Cerenkov::DoIt -- NumberOfSecondaries = "
345  }
346 
347  return pParticleChange;
348 }
349 
350 // BuildThePhysicsTable for the Cerenkov process
351 // ---------------------------------------------
352 //
353 
355 {
356  if (thePhysicsTable) return;
357 
358  const G4MaterialTable* theMaterialTable=
360  G4int numOfMaterials = G4Material::GetNumberOfMaterials();
361 
362  // create new physics table
363 
364  thePhysicsTable = new G4PhysicsTable(numOfMaterials);
365 
366  // loop for materials
367 
368  for (G4int i=0 ; i < numOfMaterials; i++) {
369 
370  G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector = 0;
371 
372  // Retrieve vector of refraction indices for the material
373  // from the material's optical properties table
374 
375  G4Material* aMaterial = (*theMaterialTable)[i];
376 
377  G4MaterialPropertiesTable* aMaterialPropertiesTable =
378  aMaterial->GetMaterialPropertiesTable();
379 
380  if (aMaterialPropertiesTable) {
381  aPhysicsOrderedFreeVector = new G4PhysicsOrderedFreeVector();
382  G4MaterialPropertyVector* theRefractionIndexVector =
383  aMaterialPropertiesTable->GetProperty(kRINDEX);
384 
385  if (theRefractionIndexVector) {
386 
387  // Retrieve the first refraction index in vector
388  // of (photon energy, refraction index) pairs
389 
390  G4double currentRI = (*theRefractionIndexVector)[0];
391 
392  if (currentRI > 1.0) {
393 
394  // Create first (photon energy, Cerenkov Integral)
395  // pair
396 
397  G4double currentPM = theRefractionIndexVector->Energy(0);
398  G4double currentCAI = 0.0;
399 
400  aPhysicsOrderedFreeVector->InsertValues(currentPM , currentCAI);
401 
402  // Set previous values to current ones prior to loop
403 
404  G4double prevPM = currentPM;
405  G4double prevCAI = currentCAI;
406  G4double prevRI = currentRI;
407 
408  // loop over all (photon energy, refraction index)
409  // pairs stored for this material
410 
411  for (size_t ii = 1;
412  ii < theRefractionIndexVector->GetVectorLength();
413  ++ii) {
414  currentRI = (*theRefractionIndexVector)[ii];
415  currentPM = theRefractionIndexVector->Energy(ii);
416 
417  currentCAI = 0.5*(1.0/(prevRI*prevRI) +
418  1.0/(currentRI*currentRI));
419 
420  currentCAI = prevCAI + (currentPM - prevPM) * currentCAI;
421 
422  aPhysicsOrderedFreeVector->
423  InsertValues(currentPM, currentCAI);
424 
425  prevPM = currentPM;
426  prevCAI = currentCAI;
427  prevRI = currentRI;
428  }
429 
430  }
431  }
432  }
433 
434  // The Cerenkov integral for a given material
435  // will be inserted in thePhysicsTable
436  // according to the position of the material in
437  // the material table.
438 
439  thePhysicsTable->insertAt(i,aPhysicsOrderedFreeVector);
440 
441  }
442 }
443 
444 // GetMeanFreePath
445 // ---------------
446 //
447 
449  G4double,
451 {
452  return 1.;
453 }
454 
456  const G4Track& aTrack,
457  G4double,
459 {
460  *condition = NotForced;
461  G4double StepLimit = DBL_MAX;
462  fNumPhotons = 0;
463 
464  const G4Material* aMaterial = aTrack.GetMaterial();
465  G4int materialIndex = aMaterial->GetIndex();
466 
467  // If Physics Vector is not defined no Cerenkov photons
468  // this check avoid string comparison below
469 
470  if(!(*thePhysicsTable)[materialIndex]) { return StepLimit; }
471 
472  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
473  const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
474 
475  G4double kineticEnergy = aParticle->GetKineticEnergy();
476  const G4ParticleDefinition* particleType = aParticle->GetDefinition();
477  G4double mass = particleType->GetPDGMass();
478 
479  // particle beta
480  G4double beta = aParticle->GetTotalMomentum() /
481  aParticle->GetTotalEnergy();
482  // particle gamma
483  G4double gamma = aParticle->GetTotalEnergy()/mass;
484 
485  G4MaterialPropertiesTable* aMaterialPropertiesTable =
486  aMaterial->GetMaterialPropertiesTable();
487 
488  G4MaterialPropertyVector* Rindex = NULL;
489 
490  if (aMaterialPropertiesTable)
491  Rindex = aMaterialPropertiesTable->GetProperty(kRINDEX);
492 
493  G4double nMax;
494  if (Rindex) {
495  nMax = Rindex->GetMaxValue();
496  } else {
497  return StepLimit;
498  }
499 
500  G4double BetaMin = 1./nMax;
501  if ( BetaMin >= 1. ) return StepLimit;
502 
503  G4double GammaMin = 1./std::sqrt(1.-BetaMin*BetaMin);
504 
505  if (gamma < GammaMin ) return StepLimit;
506 
507  G4double kinEmin = mass*(GammaMin-1.);
508 
509  G4double RangeMin = G4LossTableManager::Instance()->GetRange(particleType,
510  kinEmin,
511  couple);
513  kineticEnergy,
514  couple);
515 
516  G4double Step = Range - RangeMin;
517 
518  // If the step is smaller than 1e-16 mm, it may happen that the particle
519  // does not move. See bug 1992.
520  // 2019-03-11: change to 1e-15
521  if (Step < 1.e-15*mm) return StepLimit;
522 
523  if (Step < StepLimit) StepLimit = Step;
524  // If user has defined an average maximum number of photons to
525  // be generated in a Step, then calculate the Step length for
526  // that number of photons.
527 
528  if (fMaxPhotons > 0) {
529 
530  // particle charge
531  const G4double charge = aParticle->GetDefinition()->GetPDGCharge();
532 
533  G4double MeanNumberOfPhotons =
534  GetAverageNumberOfPhotons(charge,beta,aMaterial,Rindex);
535 
536  Step = 0.;
537  if (MeanNumberOfPhotons > 0.0) Step = fMaxPhotons / MeanNumberOfPhotons;
538 
539  if (Step > 0. && Step < StepLimit) StepLimit = Step;
540  }
541 
542  // If user has defined an maximum allowed change in beta per step
543  if (fMaxBetaChange > 0.) {
544 
545  G4double dedx = G4LossTableManager::Instance()->GetDEDX(particleType,
546  kineticEnergy,
547  couple);
548 
549  G4double deltaGamma = gamma - 1./std::sqrt(1.-beta*beta*
550  (1.-fMaxBetaChange)*
551  (1.-fMaxBetaChange));
552 
553  Step = mass * deltaGamma / dedx;
554 
555  if (Step > 0. && Step < StepLimit) StepLimit = Step;
556 
557  }
558 
559  *condition = StronglyForced;
560  return StepLimit;
561 }
562 
563 // GetAverageNumberOfPhotons
564 // -------------------------
565 // This routine computes the number of Cerenkov photons produced per
566 // GEANT-unit (millimeter) in the current medium.
567 // ^^^^^^^^^^
568 
569 G4double
571  const G4double beta,
572  const G4Material* aMaterial,
573  G4MaterialPropertyVector* Rindex) const
574 {
575  const G4double Rfact = 369.81/(eV * cm);
576 
577  if(beta <= 0.0)return 0.0;
578 
579  G4double BetaInverse = 1./beta;
580 
581  // Vectors used in computation of Cerenkov Angle Integral:
582  // - Refraction Indices for the current material
583  // - new G4PhysicsOrderedFreeVector allocated to hold CAI's
584 
585  G4int materialIndex = aMaterial->GetIndex();
586 
587  // Retrieve the Cerenkov Angle Integrals for this material
588 
589  G4PhysicsOrderedFreeVector* CerenkovAngleIntegrals =
590  (G4PhysicsOrderedFreeVector*)((*thePhysicsTable)(materialIndex));
591 
592  if(!(CerenkovAngleIntegrals->IsFilledVectorExist()))return 0.0;
593 
594  // Min and Max photon energies
595  G4double Pmin = Rindex->GetMinLowEdgeEnergy();
596  G4double Pmax = Rindex->GetMaxLowEdgeEnergy();
597 
598  // Min and Max Refraction Indices
599  G4double nMin = Rindex->GetMinValue();
600  G4double nMax = Rindex->GetMaxValue();
601 
602  // Max Cerenkov Angle Integral
603  G4double CAImax = CerenkovAngleIntegrals->GetMaxValue();
604 
605  G4double dp, ge;
606 
607  // If n(Pmax) < 1/Beta -- no photons generated
608 
609  if (nMax < BetaInverse) {
610  dp = 0.0;
611  ge = 0.0;
612  }
613 
614  // otherwise if n(Pmin) >= 1/Beta -- photons generated
615 
616  else if (nMin > BetaInverse) {
617  dp = Pmax - Pmin;
618  ge = CAImax;
619  }
620 
621  // If n(Pmin) < 1/Beta, and n(Pmax) >= 1/Beta, then
622  // we need to find a P such that the value of n(P) == 1/Beta.
623  // Interpolation is performed by the GetEnergy() and
624  // Value() methods of the G4MaterialPropertiesTable and
625  // the GetValue() method of G4PhysicsVector.
626 
627  else {
628  Pmin = Rindex->GetEnergy(BetaInverse);
629  dp = Pmax - Pmin;
630 
631  // need boolean for current implementation of G4PhysicsVector
632  // ==> being phased out
633  G4bool isOutRange;
634  G4double CAImin = CerenkovAngleIntegrals->GetValue(Pmin, isOutRange);
635  ge = CAImax - CAImin;
636 
637  if (verboseLevel>0) {
638  G4cout << "CAImin = " << CAImin << G4endl;
639  G4cout << "ge = " << ge << G4endl;
640  }
641  }
642 
643  // Calculate number of photons
644  G4double NumPhotons = Rfact * charge/eplus * charge/eplus *
645  (dp - ge * BetaInverse*BetaInverse);
646 
647  return NumPhotons;
648 }
649 
651 {
652  G4int PhysicsTableSize = thePhysicsTable->entries();
654 
655  for (G4int i = 0 ; i < PhysicsTableSize ; i++ ) {
657  v->DumpValues();
658  }
659 }
660