Geant4  10.03
G4Scintillation.cc
Go to the documentation of this file.
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 // $Id: G4Scintillation.cc 98002 2016-06-30 13:03:36Z gcosmo $
28 //
30 // Scintillation Light Class Implementation
32 //
33 // File: G4Scintillation.cc
34 // Description: RestDiscrete Process - Generation of Scintillation Photons
35 // Version: 1.0
36 // Created: 1998-11-07
37 // Author: Peter Gumplinger
38 // Updated: 2010-10-20 Allow the scintillation yield to be a function
39 // of energy deposited by particle type
40 // Thanks to Zach Hartwig (Department of Nuclear
41 // Science and Engineeering - MIT)
42 // 2010-09-22 by Peter Gumplinger
43 // > scintillation rise time included, thanks to
44 // > Martin Goettlich/DESY
45 // 2005-08-17 by Peter Gumplinger
46 // > change variable name MeanNumPhotons -> MeanNumberOfPhotons
47 // 2005-07-28 by Peter Gumplinger
48 // > add G4ProcessType to constructor
49 // 2004-08-05 by Peter Gumplinger
50 // > changed StronglyForced back to Forced in GetMeanLifeTime
51 // 2002-11-21 by Peter Gumplinger
52 // > change to use G4Poisson for small MeanNumberOfPhotons
53 // 2002-11-07 by Peter Gumplinger
54 // > now allow for fast and slow scintillation component
55 // 2002-11-05 by Peter Gumplinger
56 // > now use scintillation constants from G4Material
57 // 2002-05-09 by Peter Gumplinger
58 // > use only the PostStepPoint location for the origin of
59 // scintillation photons when energy is lost to the medium
60 // by a neutral particle
61 // 2000-09-18 by Peter Gumplinger
62 // > change: aSecondaryPosition=x0+rand*aStep.GetDeltaPosition();
63 // aSecondaryTrack->SetTouchable(0);
64 // 2001-09-17, migration of Materials to pure STL (mma)
65 // 2003-06-03, V.Ivanchenko fix compilation warnings
66 //
67 // mail: gum@triumf.ca
68 //
70 
71 #include "G4ios.hh"
72 #include "globals.hh"
73 #include "G4PhysicalConstants.hh"
74 #include "G4SystemOfUnits.hh"
75 #include "G4ParticleTypes.hh"
76 #include "G4EmProcessSubType.hh"
78 
79 #include "G4Scintillation.hh"
80 
82 // Class Implementation
84 
86  // static data members
88 
89 
90 //G4bool G4Scintillation::fTrackSecondariesFirst = false;
91 //G4bool G4Scintillation::fFiniteRiseTime = false;
92 //G4double G4Scintillation::fYieldFactor = 1.0;
93 //G4double G4Scintillation::fExcitationRatio = 1.0;
94 //G4bool G4Scintillation::fScintillationByParticleType = false;
95 //G4bool G4Scintillation::fScintillationTrackInfo = false;
96 //G4bool G4Scintillation::fStackingFlag = true;
97 //G4EmSaturation* G4Scintillation::fEmSaturation = NULL;
98 
100  // Operators
102 
103 // G4Scintillation::operator=(const G4Scintillation &right)
104 // {
105 // }
106 
108  // Constructors
110 
112  G4ProcessType type)
113  : G4VRestDiscreteProcess(processName, type) ,
114  fTrackSecondariesFirst(false),
115  fFiniteRiseTime(false),
116  fYieldFactor(1.0),
117  fExcitationRatio(1.0),
118  fScintillationByParticleType(false),
119  fScintillationTrackInfo(false),
120  fStackingFlag(true),
121  fNumPhotons(0),
122  fEmSaturation(nullptr)
123 {
125 
126 #ifdef G4DEBUG_SCINTILLATION
127  ScintTrackEDep = 0.;
128  ScintTrackYield = 0.;
129 #endif
130 
131  fFastIntegralTable = nullptr;
132  fSlowIntegralTable = nullptr;
133 
134  if (verboseLevel>0) {
135  G4cout << GetProcessName() << " is created " << G4endl;
136  }
137 }
138 
140  // Destructors
142 
144 {
145  if (fFastIntegralTable != nullptr) {
147  delete fFastIntegralTable;
148  }
149  if (fSlowIntegralTable != nullptr) {
151  delete fSlowIntegralTable;
152  }
153 }
154 
156  // Methods
158 
160 {
161  if (fFastIntegralTable != nullptr) {
163  delete fFastIntegralTable;
164  fFastIntegralTable = nullptr;
165  }
166  if (fSlowIntegralTable != nullptr) {
168  delete fSlowIntegralTable;
169  fSlowIntegralTable = nullptr;
170  }
172 }
173 
174 
175 // AtRestDoIt
176 // ----------
177 //
179 G4Scintillation::AtRestDoIt(const G4Track& aTrack, const G4Step& aStep)
180 
181 // This routine simply calls the equivalent PostStepDoIt since all the
182 // necessary information resides in aStep.GetTotalEnergyDeposit()
183 
184 {
185  return G4Scintillation::PostStepDoIt(aTrack, aStep);
186 }
187 
188 // PostStepDoIt
189 // -------------
190 //
192 G4Scintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
193 
194 // This routine is called for each tracking step of a charged particle
195 // in a scintillator. A Poisson/Gauss-distributed number of photons is
196 // generated according to the scintillation yield formula, distributed
197 // evenly along the track segment and uniformly into 4pi.
198 
199 {
200  aParticleChange.Initialize(aTrack);
201 
202  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
203  const G4Material* aMaterial = aTrack.GetMaterial();
204 
205  G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
206  G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
207 
208  G4ThreeVector x0 = pPreStepPoint->GetPosition();
209  G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
210  G4double t0 = pPreStepPoint->GetGlobalTime();
211 
212  G4double TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();
213 
214  G4MaterialPropertiesTable* aMaterialPropertiesTable =
215  aMaterial->GetMaterialPropertiesTable();
216  if (!aMaterialPropertiesTable)
217  return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
218 
219  G4MaterialPropertyVector* Fast_Intensity =
220  aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
221  G4MaterialPropertyVector* Slow_Intensity =
222  aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
223 
224  if (!Fast_Intensity && !Slow_Intensity )
225  return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
226 
227  G4int nscnt = 1;
228  if (Fast_Intensity && Slow_Intensity) nscnt = 2;
229 
230  G4double ScintillationYield = 0.;
231 
232  // Scintillation depends on particle type, energy deposited
234 
235  ScintillationYield =
237 
238  // The default linear scintillation process
239  } else {
240 
241  ScintillationYield = aMaterialPropertiesTable->
242  GetConstProperty("SCINTILLATIONYIELD");
243 
244  // Units: [# scintillation photons / MeV]
245  ScintillationYield *= fYieldFactor;
246  }
247 
248  G4double ResolutionScale = aMaterialPropertiesTable->
249  GetConstProperty("RESOLUTIONSCALE");
250 
251  // Birks law saturation:
252 
253  //G4double constBirks = 0.0;
254 
255  //constBirks = aMaterial->GetIonisation()->GetBirksConstant();
256 
257  G4double MeanNumberOfPhotons;
258 
259  // Birk's correction via fEmSaturation and specifying scintillation by
260  // by particle type are physically mutually exclusive
261 
263  MeanNumberOfPhotons = ScintillationYield;
264  else if (fEmSaturation)
265  MeanNumberOfPhotons = ScintillationYield*
267  else
268  MeanNumberOfPhotons = ScintillationYield*TotalEnergyDeposit;
269 
270  if (MeanNumberOfPhotons > 10.)
271  {
272  G4double sigma = ResolutionScale * std::sqrt(MeanNumberOfPhotons);
273  fNumPhotons=G4int(G4RandGauss::shoot(MeanNumberOfPhotons,sigma)+0.5);
274  }
275  else
276  {
277  fNumPhotons = G4int(G4Poisson(MeanNumberOfPhotons));
278  }
279 
280  if ( fNumPhotons <= 0 || !fStackingFlag )
281  {
282  // return unchanged particle and no secondaries
283 
285 
286  return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
287  }
288 
290 
292 
294  if (aTrack.GetTrackStatus() == fAlive )
296  }
297 
299 
300  G4int materialIndex = aMaterial->GetIndex();
301 
302  // Retrieve the Scintillation Integral for this material
303  // new G4PhysicsOrderedFreeVector allocated to hold CII's
304 
305  G4int Num = fNumPhotons;
306 
307  for (G4int scnt = 1; scnt <= nscnt; scnt++) {
308 
309  G4double ScintillationTime = 0.*ns;
310  G4double ScintillationRiseTime = 0.*ns;
311  G4PhysicsOrderedFreeVector* ScintillationIntegral = nullptr;
312  G4ScintillationType ScintillationType = Slow;
313 
314  if (scnt == 1) {
315  if (nscnt == 1) {
316  if (Fast_Intensity) {
317  ScintillationTime = aMaterialPropertiesTable->
318  GetConstProperty("FASTTIMECONSTANT");
319  if (fFiniteRiseTime) {
320  ScintillationRiseTime = aMaterialPropertiesTable->
321  GetConstProperty("FASTSCINTILLATIONRISETIME");
322  }
323  ScintillationType = Fast;
324  ScintillationIntegral =
326  ((*fFastIntegralTable)(materialIndex));
327  }
328  if (Slow_Intensity) {
329  ScintillationTime = aMaterialPropertiesTable->
330  GetConstProperty("SLOWTIMECONSTANT");
331  if (fFiniteRiseTime) {
332  ScintillationRiseTime = aMaterialPropertiesTable->
333  GetConstProperty("SLOWSCINTILLATIONRISETIME");
334  }
335  ScintillationType = Slow;
336  ScintillationIntegral =
338  ((*fSlowIntegralTable)(materialIndex));
339  }
340  }
341  else {
342  G4double yieldRatio = aMaterialPropertiesTable->
343  GetConstProperty("YIELDRATIO");
344  if ( fExcitationRatio == 1.0 || fExcitationRatio == 0.0) {
345  Num = G4int (std::min(yieldRatio,1.0) * fNumPhotons);
346  }
347  else {
349  }
350  ScintillationTime = aMaterialPropertiesTable->
351  GetConstProperty("FASTTIMECONSTANT");
352  if (fFiniteRiseTime) {
353  ScintillationRiseTime = aMaterialPropertiesTable->
354  GetConstProperty("FASTSCINTILLATIONRISETIME");
355  }
356  ScintillationType = Fast;
357  ScintillationIntegral =
359  ((*fFastIntegralTable)(materialIndex));
360  }
361  }
362  else {
363  Num = fNumPhotons - Num;
364  ScintillationTime = aMaterialPropertiesTable->
365  GetConstProperty("SLOWTIMECONSTANT");
366  if (fFiniteRiseTime) {
367  ScintillationRiseTime = aMaterialPropertiesTable->
368  GetConstProperty("SLOWSCINTILLATIONRISETIME");
369  }
370  ScintillationType = Slow;
371  ScintillationIntegral =
373  ((*fSlowIntegralTable)(materialIndex));
374  }
375 
376  if (!ScintillationIntegral) continue;
377 
378  // Max Scintillation Integral
379 
380  G4double CIImax = ScintillationIntegral->GetMaxValue();
381 
382  for (G4int i = 0; i < Num; i++) {
383 
384  // Determine photon energy
385 
386  G4double CIIvalue = G4UniformRand()*CIImax;
387  G4double sampledEnergy =
388  ScintillationIntegral->GetEnergy(CIIvalue);
389 
390  if (verboseLevel>1) {
391  G4cout << "sampledEnergy = " << sampledEnergy << G4endl;
392  G4cout << "CIIvalue = " << CIIvalue << G4endl;
393  }
394 
395  // Generate random photon direction
396 
397  G4double cost = 1. - 2.*G4UniformRand();
398  G4double sint = std::sqrt((1.-cost)*(1.+cost));
399 
400  G4double phi = twopi*G4UniformRand();
401  G4double sinp = std::sin(phi);
402  G4double cosp = std::cos(phi);
403 
404  G4double px = sint*cosp;
405  G4double py = sint*sinp;
406  G4double pz = cost;
407 
408  // Create photon momentum direction vector
409 
410  G4ParticleMomentum photonMomentum(px, py, pz);
411 
412  // Determine polarization of new photon
413 
414  G4double sx = cost*cosp;
415  G4double sy = cost*sinp;
416  G4double sz = -sint;
417 
418  G4ThreeVector photonPolarization(sx, sy, sz);
419 
420  G4ThreeVector perp = photonMomentum.cross(photonPolarization);
421 
422  phi = twopi*G4UniformRand();
423  sinp = std::sin(phi);
424  cosp = std::cos(phi);
425 
426  photonPolarization = cosp * photonPolarization + sinp * perp;
427 
428  photonPolarization = photonPolarization.unit();
429 
430  // Generate a new photon:
431 
432  G4DynamicParticle* aScintillationPhoton =
434  photonMomentum);
435  aScintillationPhoton->SetPolarization
436  (photonPolarization.x(),
437  photonPolarization.y(),
438  photonPolarization.z());
439 
440  aScintillationPhoton->SetKineticEnergy(sampledEnergy);
441 
442  // Generate new G4Track object:
443 
444  G4double rand;
445 
446  if (aParticle->GetDefinition()->GetPDGCharge() != 0) {
447  rand = G4UniformRand();
448  } else {
449  rand = 1.0;
450  }
451 
452  G4double delta = rand * aStep.GetStepLength();
453  G4double deltaTime = delta / (pPreStepPoint->GetVelocity()+
454  rand*(pPostStepPoint->GetVelocity()-
455  pPreStepPoint->GetVelocity())/2.);
456 
457  // emission time distribution
458  if (ScintillationRiseTime==0.0) {
459  deltaTime = deltaTime -
460  ScintillationTime * std::log( G4UniformRand() );
461  } else {
462  deltaTime = deltaTime +
463  sample_time(ScintillationRiseTime, ScintillationTime);
464  }
465 
466  G4double aSecondaryTime = t0 + deltaTime;
467 
468  G4ThreeVector aSecondaryPosition =
469  x0 + rand * aStep.GetDeltaPosition();
470 
471  G4Track* aSecondaryTrack = new G4Track(aScintillationPhoton,
472  aSecondaryTime,
473  aSecondaryPosition);
474 
475  aSecondaryTrack->SetTouchableHandle(
477  // aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0);
478 
479  aSecondaryTrack->SetParentID(aTrack.GetTrackID());
480 
481  if (fScintillationTrackInfo) aSecondaryTrack->
482  SetUserInformation(new G4ScintillationTrackInformation(ScintillationType));
483 
484  aParticleChange.AddSecondary(aSecondaryTrack);
485 
486  }
487  }
488 
489  if (verboseLevel>0) {
490  G4cout << "\n Exiting from G4Scintillation::DoIt -- NumberOfSecondaries = "
492  }
493 
494  return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
495 }
496 
497 // BuildThePhysicsTable for the scintillation process
498 // --------------------------------------------------
499 //
500 
502 {
503  if (fFastIntegralTable && fSlowIntegralTable) return;
504 
505  const G4MaterialTable* theMaterialTable =
507  G4int numOfMaterials = G4Material::GetNumberOfMaterials();
508 
509  // create new physics table
510 
512  new G4PhysicsTable(numOfMaterials);
514  new G4PhysicsTable(numOfMaterials);
515 
516  // loop for materials
517 
518  for (G4int i=0 ; i < numOfMaterials; i++)
519  {
520  G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
522  G4PhysicsOrderedFreeVector* bPhysicsOrderedFreeVector =
524 
525  // Retrieve vector of scintillation wavelength intensity for
526  // the material from the material's optical properties table.
527 
528  G4Material* aMaterial = (*theMaterialTable)[i];
529 
530  G4MaterialPropertiesTable* aMaterialPropertiesTable =
531  aMaterial->GetMaterialPropertiesTable();
532 
533  if (aMaterialPropertiesTable) {
534 
535  G4MaterialPropertyVector* theFastLightVector =
536  aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
537 
538  if (theFastLightVector) {
539 
540  // Retrieve the first intensity point in vector
541  // of (photon energy, intensity) pairs
542 
543  G4double currentIN = (*theFastLightVector)[0];
544 
545  if (currentIN >= 0.0) {
546 
547  // Create first (photon energy, Scintillation
548  // Integral pair
549 
550  G4double currentPM = theFastLightVector->Energy(0);
551 
552  G4double currentCII = 0.0;
553 
554  aPhysicsOrderedFreeVector->
555  InsertValues(currentPM , currentCII);
556 
557  // Set previous values to current ones prior to loop
558 
559  G4double prevPM = currentPM;
560  G4double prevCII = currentCII;
561  G4double prevIN = currentIN;
562 
563  // loop over all (photon energy, intensity)
564  // pairs stored for this material
565 
566  for (size_t ii = 1;
567  ii < theFastLightVector->GetVectorLength();
568  ++ii)
569  {
570  currentPM = theFastLightVector->Energy(ii);
571  currentIN = (*theFastLightVector)[ii];
572 
573  currentCII = 0.5 * (prevIN + currentIN);
574 
575  currentCII = prevCII +
576  (currentPM - prevPM) * currentCII;
577 
578  aPhysicsOrderedFreeVector->
579  InsertValues(currentPM, currentCII);
580 
581  prevPM = currentPM;
582  prevCII = currentCII;
583  prevIN = currentIN;
584  }
585 
586  }
587  }
588 
589  G4MaterialPropertyVector* theSlowLightVector =
590  aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
591 
592  if (theSlowLightVector) {
593 
594  // Retrieve the first intensity point in vector
595  // of (photon energy, intensity) pairs
596 
597  G4double currentIN = (*theSlowLightVector)[0];
598 
599  if (currentIN >= 0.0) {
600 
601  // Create first (photon energy, Scintillation
602  // Integral pair
603 
604  G4double currentPM = theSlowLightVector->Energy(0);
605 
606  G4double currentCII = 0.0;
607 
608  bPhysicsOrderedFreeVector->
609  InsertValues(currentPM , currentCII);
610 
611  // Set previous values to current ones prior to loop
612 
613  G4double prevPM = currentPM;
614  G4double prevCII = currentCII;
615  G4double prevIN = currentIN;
616 
617  // loop over all (photon energy, intensity)
618  // pairs stored for this material
619 
620  for (size_t ii = 1;
621  ii < theSlowLightVector->GetVectorLength();
622  ++ii)
623  {
624  currentPM = theSlowLightVector->Energy(ii);
625  currentIN = (*theSlowLightVector)[ii];
626 
627  currentCII = 0.5 * (prevIN + currentIN);
628 
629  currentCII = prevCII +
630  (currentPM - prevPM) * currentCII;
631 
632  bPhysicsOrderedFreeVector->
633  InsertValues(currentPM, currentCII);
634 
635  prevPM = currentPM;
636  prevCII = currentCII;
637  prevIN = currentIN;
638  }
639 
640  }
641  }
642  }
643 
644  // The scintillation integral(s) for a given material
645  // will be inserted in the table(s) according to the
646  // position of the material in the material table.
647 
648  fFastIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
649  fSlowIntegralTable->insertAt(i,bPhysicsOrderedFreeVector);
650 
651  }
652 }
653 
655 {
656  if (fEmSaturation) {
657  G4Exception("G4Scintillation::SetScintillationByParticleType", "Scint02",
658  JustWarning, "Redefinition: Birks Saturation is replaced by ScintillationByParticleType!");
660  }
661  fScintillationByParticleType = scintType;
662 }
663 
664 // GetMeanFreePath
665 // ---------------
666 //
667 
669  G4double ,
671 {
672  *condition = StronglyForced;
673 
674  return DBL_MAX;
675 
676 }
677 
678 // GetMeanLifeTime
679 // ---------------
680 //
681 
684 {
685  *condition = Forced;
686 
687  return DBL_MAX;
688 
689 }
690 
692 {
693 // tau1: rise time and tau2: decay time
694 
695  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
696  while(1) {
697  // two random numbers
698  G4double ran1 = G4UniformRand();
699  G4double ran2 = G4UniformRand();
700  //
701  // exponential distribution as envelope function: very efficient
702  //
703  G4double d = (tau1+tau2)/tau2;
704  // make sure the envelope function is
705  // always larger than the bi-exponential
706  G4double t = -1.0*tau2*std::log(1-ran1);
707  G4double gg = d*single_exp(t,tau2);
708  if (ran2 <= bi_exp(t,tau1,tau2)/gg) return t;
709  }
710  return -1.0;
711 }
712 
715 {
717  // Get the scintillation yield vector //
719 
721 
722  G4MaterialPropertyVector *Scint_Yield_Vector = nullptr;
723 
724  G4MaterialPropertiesTable *aMaterialPropertiesTable
726 
727  // Get the G4MaterialPropertyVector containing the scintillation
728  // yield as a function of the energy deposited and particle type
729 
730  // Protons
731  if(pDef==G4Proton::ProtonDefinition())
732  Scint_Yield_Vector = aMaterialPropertiesTable->
733  GetProperty("PROTONSCINTILLATIONYIELD");
734 
735  // Deuterons
736  else if(pDef==G4Deuteron::DeuteronDefinition())
737  Scint_Yield_Vector = aMaterialPropertiesTable->
738  GetProperty("DEUTERONSCINTILLATIONYIELD");
739 
740  // Tritons
741  else if(pDef==G4Triton::TritonDefinition())
742  Scint_Yield_Vector = aMaterialPropertiesTable->
743  GetProperty("TRITONSCINTILLATIONYIELD");
744 
745  // Alphas
746  else if(pDef==G4Alpha::AlphaDefinition())
747  Scint_Yield_Vector = aMaterialPropertiesTable->
748  GetProperty("ALPHASCINTILLATIONYIELD");
749 
750  // Ions (particles derived from G4VIon and G4Ions) and recoil ions
751  // below the production cut from neutrons after hElastic
752  else if(pDef->GetParticleType()== "nucleus" ||
754  Scint_Yield_Vector = aMaterialPropertiesTable->
755  GetProperty("IONSCINTILLATIONYIELD");
756 
757  // Electrons (must also account for shell-binding energy
758  // attributed to gamma from standard photoelectric effect)
759  else if(pDef==G4Electron::ElectronDefinition() ||
760  pDef==G4Gamma::GammaDefinition())
761  Scint_Yield_Vector = aMaterialPropertiesTable->
762  GetProperty("ELECTRONSCINTILLATIONYIELD");
763 
764  // Default for particles not enumerated/listed above
765  else
766  Scint_Yield_Vector = aMaterialPropertiesTable->
767  GetProperty("ELECTRONSCINTILLATIONYIELD");
768 
769  // If the user has specified none of the above particles then the
770  // default is the electron scintillation yield
771  if(!Scint_Yield_Vector)
772  Scint_Yield_Vector = aMaterialPropertiesTable->
773  GetProperty("ELECTRONSCINTILLATIONYIELD");
774 
775  // Throw an exception if no scintillation yield vector is found
776  if (!Scint_Yield_Vector) {
778  ed << "\nG4Scintillation::PostStepDoIt(): "
779  << "Request for scintillation yield for energy deposit and particle\n"
780  << "type without correct entry in MaterialPropertiesTable.\n"
781  << "ScintillationByParticleType requires at minimum that \n"
782  << "ELECTRONSCINTILLATIONYIELD is set by the user\n"
783  << G4endl;
784  G4String comments = "Missing MaterialPropertiesTable entry - No correct entry in MaterialPropertiesTable";
785  G4Exception("G4Scintillation::PostStepDoIt","Scint01",
786  FatalException,ed,comments);
787  }
788 
790  // Calculate the scintillation light //
792  // To account for potential nonlinearity and scintillation photon
793  // density along the track, light (L) is produced according to:
794  //
795  // L_currentStep = L(PreStepKE) - L(PreStepKE - EDep)
796 
797  G4double ScintillationYield = 0.;
798 
799  G4double StepEnergyDeposit = aStep.GetTotalEnergyDeposit();
800  G4double PreStepKineticEnergy = aStep.GetPreStepPoint()->GetKineticEnergy();
801 
802  if(PreStepKineticEnergy <= Scint_Yield_Vector->GetMaxEnergy()){
803  G4double Yield1 = Scint_Yield_Vector->Value(PreStepKineticEnergy);
804  G4double Yield2 = Scint_Yield_Vector->
805  Value(PreStepKineticEnergy - StepEnergyDeposit);
806  ScintillationYield = Yield1 - Yield2;
807  } else {
809  ed << "\nG4Scintillation::GetScintillationYieldByParticleType(): Request\n"
810  << "for scintillation light yield above the available energy range\n"
811  << "specifed in G4MaterialPropertiesTable. A linear interpolation\n"
812  << "will be performed to compute the scintillation light yield using\n"
813  << "(L_max / E_max) as the photon yield per unit energy."
814  << G4endl;
815  G4String cmt = "\nScintillation yield may be unphysical!\n";
816  G4Exception("G4Scintillation::GetScintillationYieldByParticleType()",
817  "Scint03", JustWarning, ed, cmt);
818 
819  G4double LinearYield = Scint_Yield_Vector->GetMaxValue()
820  / Scint_Yield_Vector->GetMaxEnergy();
821 
822  // Units: [# scintillation photons]
823  ScintillationYield = LinearYield * StepEnergyDeposit;
824  }
825 
826 #ifdef G4DEBUG_SCINTILLATION
827 
828  // Increment track aggregators
829  ScintTrackYield += ScintillationYield;
830  ScintTrackEDep += StepEnergyDeposit;
831 
832  G4cout << "\n--- G4Scintillation::GetScintillationYieldByParticleType() ---\n"
833  << "--\n"
834  << "-- Name = " << aTrack.GetParticleDefinition()->GetParticleName() << "\n"
835  << "-- TrackID = " << aTrack.GetTrackID() << "\n"
836  << "-- ParentID = " << aTrack.GetParentID() << "\n"
837  << "-- Current KE = " << aTrack.GetKineticEnergy()/MeV << " MeV\n"
838  << "-- Step EDep = " << aStep.GetTotalEnergyDeposit()/MeV << " MeV\n"
839  << "-- Track EDep = " << ScintTrackEDep/MeV << " MeV\n"
840  << "-- Vertex KE = " << aTrack.GetVertexKineticEnergy()/MeV << " MeV\n"
841  << "-- Step yield = " << ScintillationYield << " photons\n"
842  << "-- Track yield = " << ScintTrackYield << " photons\n"
843  << G4endl;
844 
845  // The track has terminated within or has left the scintillator volume
846  if( (aTrack.GetTrackStatus() == fStopButAlive) or
847  (aStep.GetPostStepPoint()->GetStepStatus() == fGeomBoundary) ){
848 
849  // Reset aggregators for the next track
850  ScintTrackEDep = 0.;
851  ScintTrackYield = 0.;
852  }
853 
854 #endif
855 
856  return ScintillationYield;
857 }
void SetScintillationByParticleType(const G4bool)
G4double condition(const G4ErrorSymMatrix &m)
G4Scintillation(const G4String &processName="Scintillation", G4ProcessType type=fElectromagnetic)
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
static G4Triton * TritonDefinition()
Definition: G4Triton.cc:90
G4double bi_exp(G4double t, G4double tau1, G4double tau2)
G4int GetParentID() const
ThreeVector shoot(const G4int Ap, const G4int Af)
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4int GetNumberOfSecondaries() const
G4double GetMaxEnergy() const
G4int verboseLevel
Definition: G4VProcess.hh:368
G4bool fTrackSecondariesFirst
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
CLHEP::Hep3Vector G4ThreeVector
G4double GetStepLength() const
const G4DynamicParticle * GetDynamicParticle() const
size_t GetIndex() const
Definition: G4Material.hh:262
G4double single_exp(G4double t, G4double tau2)
G4PhysicsTable * fFastIntegralTable
G4StepStatus GetStepStatus() const
G4EmSaturation * fEmSaturation
G4TrackStatus GetTrackStatus() const
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
std::vector< G4Material * > G4MaterialTable
G4PhysicsTable * fSlowIntegralTable
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep) override
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4ParticleDefinition * GetDefinition() const
void BuildPhysicsTable(const G4ParticleDefinition &aParticleType) override
G4double GetVelocity() const
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
static constexpr double twopi
Definition: G4SIunits.hh:76
G4StepPoint * GetPreStepPoint() const
G4double GetKineticEnergy() const
G4double GetMeanLifeTime(const G4Track &aTrack, G4ForceCondition *) override
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4double GetVertexKineticEnergy() const
const G4ThreeVector & GetPosition() const
bool G4bool
Definition: G4Types.hh:79
const G4ParticleDefinition * GetParticleDefinition() const
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
const G4String & GetParticleType() const
Definition: G4Step.hh:76
G4int GetTrackID() const
void SetPolarization(G4double polX, G4double polY, G4double polZ)
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:594
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
G4double VisibleEnergyDepositionAtAStep(const G4Step *) const
void SetKineticEnergy(G4double aEnergy)
G4double GetEnergy(G4double aValue)
G4double GetTotalEnergyDeposit() const
G4Material * GetMaterial() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double fExcitationRatio
static G4OpticalPhoton * OpticalPhoton()
virtual void Initialize(const G4Track &)
G4bool fScintillationTrackInfo
G4MaterialPropertiesTable * GetMaterialPropertiesTable() const
Definition: G4Material.hh:252
void SetNumberOfSecondaries(G4int totSecondaries)
void SetParentID(const G4int aValue)
G4StepPoint * GetPostStepPoint() const
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void AddSecondary(G4Track *aSecondary)
void insertAt(size_t, G4PhysicsVector *)
#define G4endl
Definition: G4ios.hh:61
G4double GetGlobalTime() const
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *) override
G4double GetKineticEnergy() const
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4ForceCondition
G4double GetPDGCharge() const
G4ThreeVector G4ParticleMomentum
G4double sample_time(G4double tau1, G4double tau2)
G4ThreeVector GetDeltaPosition() const
#define DBL_MAX
Definition: templates.hh:83
G4double GetScintillationYieldByParticleType(const G4Track &aTrack, const G4Step &aStep)
static G4Deuteron * DeuteronDefinition()
Definition: G4Deuteron.cc:89
static G4Alpha * AlphaDefinition()
Definition: G4Alpha.cc:84
static G4Neutron * NeutronDefinition()
Definition: G4Neutron.cc:99
#define ns
Definition: xmlparse.cc:614
G4MaterialPropertyVector * GetProperty(const char *key)
const G4TouchableHandle & GetTouchableHandle() const
void clearAndDestroy()
G4VParticleChange * AtRestDoIt(const G4Track &aTrack, const G4Step &aStep) override
static G4Gamma * GammaDefinition()
Definition: G4Gamma.cc:81
G4ProcessType
G4bool fScintillationByParticleType