Geant4  10.01.p01
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 85355 2014-10-28 09:58:59Z 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"
77 
78 #include "G4Scintillation.hh"
79 
81 // Class Implementation
83 
85  // static data members
87 
88 
89 //G4bool G4Scintillation::fTrackSecondariesFirst = false;
90 //G4bool G4Scintillation::fFiniteRiseTime = false;
91 //G4double G4Scintillation::fYieldFactor = 1.0;
92 //G4double G4Scintillation::fExcitationRatio = 1.0;
93 //G4bool G4Scintillation::fScintillationByParticleType = false;
94 //G4EmSaturation* G4Scintillation::fEmSaturation = NULL;
95 
97  // Operators
99 
100 // G4Scintillation::operator=(const G4Scintillation &right)
101 // {
102 // }
103 
105  // Constructors
107 
109  G4ProcessType type)
110  : G4VRestDiscreteProcess(processName, type) ,
111  fTrackSecondariesFirst(false),
112  fFiniteRiseTime(false),
113  fYieldFactor(1.0),
114  fExcitationRatio(1.0),
115  fScintillationByParticleType(false),
116  fEmSaturation(0)
117 {
119 
120 #ifdef G4DEBUG_SCINTILLATION
121  ScintTrackEDep = 0.;
122  ScintTrackYield = 0.;
123 #endif
124 
125  fFastIntegralTable = NULL;
126  fSlowIntegralTable = NULL;
127 
128  if (verboseLevel>0) {
129  G4cout << GetProcessName() << " is created " << G4endl;
130  }
131 }
132 
134  // Destructors
136 
138 {
139  if (fFastIntegralTable != NULL) {
141  delete fFastIntegralTable;
142  }
143  if (fSlowIntegralTable != NULL) {
145  delete fSlowIntegralTable;
146  }
147 }
148 
150  // Methods
152 
154 {
155  if (fFastIntegralTable != NULL) {
157  delete fFastIntegralTable;
158  fFastIntegralTable = NULL;
159  }
160  if (fSlowIntegralTable != NULL) {
162  delete fSlowIntegralTable;
163  fSlowIntegralTable = NULL;
164  }
166 }
167 
168 
169 // AtRestDoIt
170 // ----------
171 //
173 G4Scintillation::AtRestDoIt(const G4Track& aTrack, const G4Step& aStep)
174 
175 // This routine simply calls the equivalent PostStepDoIt since all the
176 // necessary information resides in aStep.GetTotalEnergyDeposit()
177 
178 {
179  return G4Scintillation::PostStepDoIt(aTrack, aStep);
180 }
181 
182 // PostStepDoIt
183 // -------------
184 //
186 G4Scintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
187 
188 // This routine is called for each tracking step of a charged particle
189 // in a scintillator. A Poisson/Gauss-distributed number of photons is
190 // generated according to the scintillation yield formula, distributed
191 // evenly along the track segment and uniformly into 4pi.
192 
193 {
194  aParticleChange.Initialize(aTrack);
195 
196  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
197  const G4Material* aMaterial = aTrack.GetMaterial();
198 
199  G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
200  G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
201 
202  G4ThreeVector x0 = pPreStepPoint->GetPosition();
203  G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
204  G4double t0 = pPreStepPoint->GetGlobalTime();
205 
206  G4double TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();
207 
208  G4MaterialPropertiesTable* aMaterialPropertiesTable =
209  aMaterial->GetMaterialPropertiesTable();
210  if (!aMaterialPropertiesTable)
211  return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
212 
213  G4MaterialPropertyVector* Fast_Intensity =
214  aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
215  G4MaterialPropertyVector* Slow_Intensity =
216  aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
217 
218  if (!Fast_Intensity && !Slow_Intensity )
219  return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
220 
221  G4int nscnt = 1;
222  if (Fast_Intensity && Slow_Intensity) nscnt = 2;
223 
224  G4double ScintillationYield = 0.;
225 
226  // Scintillation depends on particle type, energy deposited
228 
229  ScintillationYield =
231 
232  // The default linear scintillation process
233  } else {
234 
235  ScintillationYield = aMaterialPropertiesTable->
236  GetConstProperty("SCINTILLATIONYIELD");
237 
238  // Units: [# scintillation photons / MeV]
239  ScintillationYield *= fYieldFactor;
240  }
241 
242  G4double ResolutionScale = aMaterialPropertiesTable->
243  GetConstProperty("RESOLUTIONSCALE");
244 
245  // Birks law saturation:
246 
247  //G4double constBirks = 0.0;
248 
249  //constBirks = aMaterial->GetIonisation()->GetBirksConstant();
250 
251  G4double MeanNumberOfPhotons;
252 
253  // Birk's correction via fEmSaturation and specifying scintillation by
254  // by particle type are physically mutually exclusive
255 
257  MeanNumberOfPhotons = ScintillationYield;
258  else if (fEmSaturation)
259  MeanNumberOfPhotons = ScintillationYield*
261  else
262  MeanNumberOfPhotons = ScintillationYield*TotalEnergyDeposit;
263 
264  G4int NumPhotons;
265 
266  if (MeanNumberOfPhotons > 10.)
267  {
268  G4double sigma = ResolutionScale * std::sqrt(MeanNumberOfPhotons);
269  NumPhotons = G4int(G4RandGauss::shoot(MeanNumberOfPhotons,sigma)+0.5);
270  }
271  else
272  {
273  NumPhotons = G4int(G4Poisson(MeanNumberOfPhotons));
274  }
275 
276  if (NumPhotons <= 0)
277  {
278  // return unchanged particle and no secondaries
279 
281 
282  return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
283  }
284 
286 
288 
290  if (aTrack.GetTrackStatus() == fAlive )
292  }
293 
295 
296  G4int materialIndex = aMaterial->GetIndex();
297 
298  // Retrieve the Scintillation Integral for this material
299  // new G4PhysicsOrderedFreeVector allocated to hold CII's
300 
301  G4int Num = NumPhotons;
302 
303  for (G4int scnt = 1; scnt <= nscnt; scnt++) {
304 
305  G4double ScintillationTime = 0.*ns;
306  G4double ScintillationRiseTime = 0.*ns;
307  G4PhysicsOrderedFreeVector* ScintillationIntegral = NULL;
308 
309  if (scnt == 1) {
310  if (nscnt == 1) {
311  if (Fast_Intensity) {
312  ScintillationTime = aMaterialPropertiesTable->
313  GetConstProperty("FASTTIMECONSTANT");
314  if (fFiniteRiseTime) {
315  ScintillationRiseTime = aMaterialPropertiesTable->
316  GetConstProperty("FASTSCINTILLATIONRISETIME");
317  }
318  ScintillationIntegral =
320  ((*fFastIntegralTable)(materialIndex));
321  }
322  if (Slow_Intensity) {
323  ScintillationTime = aMaterialPropertiesTable->
324  GetConstProperty("SLOWTIMECONSTANT");
325  if (fFiniteRiseTime) {
326  ScintillationRiseTime = aMaterialPropertiesTable->
327  GetConstProperty("SLOWSCINTILLATIONRISETIME");
328  }
329  ScintillationIntegral =
331  ((*fSlowIntegralTable)(materialIndex));
332  }
333  }
334  else {
335  G4double yieldRatio = aMaterialPropertiesTable->
336  GetConstProperty("YIELDRATIO");
337  if ( fExcitationRatio == 1.0 || fExcitationRatio == 0.0) {
338  Num = G4int (std::min(yieldRatio,1.0) * NumPhotons);
339  }
340  else {
341  Num = G4int (std::min(fExcitationRatio,1.0) * NumPhotons);
342  }
343  ScintillationTime = aMaterialPropertiesTable->
344  GetConstProperty("FASTTIMECONSTANT");
345  if (fFiniteRiseTime) {
346  ScintillationRiseTime = aMaterialPropertiesTable->
347  GetConstProperty("FASTSCINTILLATIONRISETIME");
348  }
349  ScintillationIntegral =
351  ((*fFastIntegralTable)(materialIndex));
352  }
353  }
354  else {
355  Num = NumPhotons - Num;
356  ScintillationTime = aMaterialPropertiesTable->
357  GetConstProperty("SLOWTIMECONSTANT");
358  if (fFiniteRiseTime) {
359  ScintillationRiseTime = aMaterialPropertiesTable->
360  GetConstProperty("SLOWSCINTILLATIONRISETIME");
361  }
362  ScintillationIntegral =
364  ((*fSlowIntegralTable)(materialIndex));
365  }
366 
367  if (!ScintillationIntegral) continue;
368 
369  // Max Scintillation Integral
370 
371  G4double CIImax = ScintillationIntegral->GetMaxValue();
372 
373  for (G4int i = 0; i < Num; i++) {
374 
375  // Determine photon energy
376 
377  G4double CIIvalue = G4UniformRand()*CIImax;
378  G4double sampledEnergy =
379  ScintillationIntegral->GetEnergy(CIIvalue);
380 
381  if (verboseLevel>1) {
382  G4cout << "sampledEnergy = " << sampledEnergy << G4endl;
383  G4cout << "CIIvalue = " << CIIvalue << G4endl;
384  }
385 
386  // Generate random photon direction
387 
388  G4double cost = 1. - 2.*G4UniformRand();
389  G4double sint = std::sqrt((1.-cost)*(1.+cost));
390 
391  G4double phi = twopi*G4UniformRand();
392  G4double sinp = std::sin(phi);
393  G4double cosp = std::cos(phi);
394 
395  G4double px = sint*cosp;
396  G4double py = sint*sinp;
397  G4double pz = cost;
398 
399  // Create photon momentum direction vector
400 
401  G4ParticleMomentum photonMomentum(px, py, pz);
402 
403  // Determine polarization of new photon
404 
405  G4double sx = cost*cosp;
406  G4double sy = cost*sinp;
407  G4double sz = -sint;
408 
409  G4ThreeVector photonPolarization(sx, sy, sz);
410 
411  G4ThreeVector perp = photonMomentum.cross(photonPolarization);
412 
413  phi = twopi*G4UniformRand();
414  sinp = std::sin(phi);
415  cosp = std::cos(phi);
416 
417  photonPolarization = cosp * photonPolarization + sinp * perp;
418 
419  photonPolarization = photonPolarization.unit();
420 
421  // Generate a new photon:
422 
423  G4DynamicParticle* aScintillationPhoton =
425  photonMomentum);
426  aScintillationPhoton->SetPolarization
427  (photonPolarization.x(),
428  photonPolarization.y(),
429  photonPolarization.z());
430 
431  aScintillationPhoton->SetKineticEnergy(sampledEnergy);
432 
433  // Generate new G4Track object:
434 
435  G4double rand;
436 
437  if (aParticle->GetDefinition()->GetPDGCharge() != 0) {
438  rand = G4UniformRand();
439  } else {
440  rand = 1.0;
441  }
442 
443  G4double delta = rand * aStep.GetStepLength();
444  G4double deltaTime = delta /
445  ((pPreStepPoint->GetVelocity()+
446  pPostStepPoint->GetVelocity())/2.);
447 
448  // emission time distribution
449  if (ScintillationRiseTime==0.0) {
450  deltaTime = deltaTime -
451  ScintillationTime * std::log( G4UniformRand() );
452  } else {
453  deltaTime = deltaTime +
454  sample_time(ScintillationRiseTime, ScintillationTime);
455  }
456 
457  G4double aSecondaryTime = t0 + deltaTime;
458 
459  G4ThreeVector aSecondaryPosition =
460  x0 + rand * aStep.GetDeltaPosition();
461 
462  G4Track* aSecondaryTrack = new G4Track(aScintillationPhoton,
463  aSecondaryTime,
464  aSecondaryPosition);
465 
466  aSecondaryTrack->SetTouchableHandle(
468  // aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0);
469 
470  aSecondaryTrack->SetParentID(aTrack.GetTrackID());
471 
472  aParticleChange.AddSecondary(aSecondaryTrack);
473 
474  }
475  }
476 
477  if (verboseLevel>0) {
478  G4cout << "\n Exiting from G4Scintillation::DoIt -- NumberOfSecondaries = "
480  }
481 
482  return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
483 }
484 
485 // BuildThePhysicsTable for the scintillation process
486 // --------------------------------------------------
487 //
488 
490 {
491  if (fFastIntegralTable && fSlowIntegralTable) return;
492 
493  const G4MaterialTable* theMaterialTable =
495  G4int numOfMaterials = G4Material::GetNumberOfMaterials();
496 
497  // create new physics table
498 
500  new G4PhysicsTable(numOfMaterials);
502  new G4PhysicsTable(numOfMaterials);
503 
504  // loop for materials
505 
506  for (G4int i=0 ; i < numOfMaterials; i++)
507  {
508  G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
510  G4PhysicsOrderedFreeVector* bPhysicsOrderedFreeVector =
512 
513  // Retrieve vector of scintillation wavelength intensity for
514  // the material from the material's optical properties table.
515 
516  G4Material* aMaterial = (*theMaterialTable)[i];
517 
518  G4MaterialPropertiesTable* aMaterialPropertiesTable =
519  aMaterial->GetMaterialPropertiesTable();
520 
521  if (aMaterialPropertiesTable) {
522 
523  G4MaterialPropertyVector* theFastLightVector =
524  aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
525 
526  if (theFastLightVector) {
527 
528  // Retrieve the first intensity point in vector
529  // of (photon energy, intensity) pairs
530 
531  G4double currentIN = (*theFastLightVector)[0];
532 
533  if (currentIN >= 0.0) {
534 
535  // Create first (photon energy, Scintillation
536  // Integral pair
537 
538  G4double currentPM = theFastLightVector->Energy(0);
539 
540  G4double currentCII = 0.0;
541 
542  aPhysicsOrderedFreeVector->
543  InsertValues(currentPM , currentCII);
544 
545  // Set previous values to current ones prior to loop
546 
547  G4double prevPM = currentPM;
548  G4double prevCII = currentCII;
549  G4double prevIN = currentIN;
550 
551  // loop over all (photon energy, intensity)
552  // pairs stored for this material
553 
554  for (size_t ii = 1;
555  ii < theFastLightVector->GetVectorLength();
556  ++ii)
557  {
558  currentPM = theFastLightVector->Energy(ii);
559  currentIN = (*theFastLightVector)[ii];
560 
561  currentCII = 0.5 * (prevIN + currentIN);
562 
563  currentCII = prevCII +
564  (currentPM - prevPM) * currentCII;
565 
566  aPhysicsOrderedFreeVector->
567  InsertValues(currentPM, currentCII);
568 
569  prevPM = currentPM;
570  prevCII = currentCII;
571  prevIN = currentIN;
572  }
573 
574  }
575  }
576 
577  G4MaterialPropertyVector* theSlowLightVector =
578  aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
579 
580  if (theSlowLightVector) {
581 
582  // Retrieve the first intensity point in vector
583  // of (photon energy, intensity) pairs
584 
585  G4double currentIN = (*theSlowLightVector)[0];
586 
587  if (currentIN >= 0.0) {
588 
589  // Create first (photon energy, Scintillation
590  // Integral pair
591 
592  G4double currentPM = theSlowLightVector->Energy(0);
593 
594  G4double currentCII = 0.0;
595 
596  bPhysicsOrderedFreeVector->
597  InsertValues(currentPM , currentCII);
598 
599  // Set previous values to current ones prior to loop
600 
601  G4double prevPM = currentPM;
602  G4double prevCII = currentCII;
603  G4double prevIN = currentIN;
604 
605  // loop over all (photon energy, intensity)
606  // pairs stored for this material
607 
608  for (size_t ii = 1;
609  ii < theSlowLightVector->GetVectorLength();
610  ++ii)
611  {
612  currentPM = theSlowLightVector->Energy(ii);
613  currentIN = (*theSlowLightVector)[ii];
614 
615  currentCII = 0.5 * (prevIN + currentIN);
616 
617  currentCII = prevCII +
618  (currentPM - prevPM) * currentCII;
619 
620  bPhysicsOrderedFreeVector->
621  InsertValues(currentPM, currentCII);
622 
623  prevPM = currentPM;
624  prevCII = currentCII;
625  prevIN = currentIN;
626  }
627 
628  }
629  }
630  }
631 
632  // The scintillation integral(s) for a given material
633  // will be inserted in the table(s) according to the
634  // position of the material in the material table.
635 
636  fFastIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
637  fSlowIntegralTable->insertAt(i,bPhysicsOrderedFreeVector);
638 
639  }
640 }
641 
643 {
644  fTrackSecondariesFirst = state;
645 }
646 
648 {
649  fFiniteRiseTime = state;
650 }
651 
653 {
654  fYieldFactor = yieldfactor;
655 }
656 
658  fExcitationRatio = excitationratio;
659 }
660 
661 // Called by the user to set the scintillation yield as a function
662 // of energy deposited by particle type
663 
665 {
666  if (fEmSaturation) {
667  G4Exception("G4Scintillation::SetScintillationByParticleType", "Scint02",
668  JustWarning, "Redefinition: Birks Saturation is replaced by ScintillationByParticleType!");
670  }
671  fScintillationByParticleType = scintType;
672 }
673 
675 {
676  fEmSaturation = sat;
677 }
678 
680 {
681  fEmSaturation = NULL;
682 }
683 
684 // GetMeanFreePath
685 // ---------------
686 //
687 
689  G4double ,
691 {
692  *condition = StronglyForced;
693 
694  return DBL_MAX;
695 
696 }
697 
698 // GetMeanLifeTime
699 // ---------------
700 //
701 
704 {
705  *condition = Forced;
706 
707  return DBL_MAX;
708 
709 }
710 
712 {
713 // tau1: rise time and tau2: decay time
714 
715  while(1) {
716  // two random numbers
717  G4double ran1 = G4UniformRand();
718  G4double ran2 = G4UniformRand();
719  //
720  // exponential distribution as envelope function: very efficient
721  //
722  G4double d = (tau1+tau2)/tau2;
723  // make sure the envelope function is
724  // always larger than the bi-exponential
725  G4double t = -1.0*tau2*std::log(1-ran1);
726  G4double gg = d*single_exp(t,tau2);
727  if (ran2 <= bi_exp(t,tau1,tau2)/gg) return t;
728  }
729  return -1.0;
730 }
731 
734 {
736  // Get the scintillation yield vector //
738 
740 
741  G4MaterialPropertyVector *Scint_Yield_Vector = NULL;
742 
743  G4MaterialPropertiesTable *aMaterialPropertiesTable
745 
746  // Get the G4MaterialPropertyVector containing the scintillation
747  // yield as a function of the energy deposited and particle type
748 
749  // Protons
750  if(pDef==G4Proton::ProtonDefinition())
751  Scint_Yield_Vector = aMaterialPropertiesTable->
752  GetProperty("PROTONSCINTILLATIONYIELD");
753 
754  // Deuterons
755  else if(pDef==G4Deuteron::DeuteronDefinition())
756  Scint_Yield_Vector = aMaterialPropertiesTable->
757  GetProperty("DEUTERONSCINTILLATIONYIELD");
758 
759  // Tritons
760  else if(pDef==G4Triton::TritonDefinition())
761  Scint_Yield_Vector = aMaterialPropertiesTable->
762  GetProperty("TRITONSCINTILLATIONYIELD");
763 
764  // Alphas
765  else if(pDef==G4Alpha::AlphaDefinition())
766  Scint_Yield_Vector = aMaterialPropertiesTable->
767  GetProperty("ALPHASCINTILLATIONYIELD");
768 
769  // Ions (particles derived from G4VIon and G4Ions) and recoil ions
770  // below the production cut from neutrons after hElastic
771  else if(pDef->GetParticleType()== "nucleus" ||
773  Scint_Yield_Vector = aMaterialPropertiesTable->
774  GetProperty("IONSCINTILLATIONYIELD");
775 
776  // Electrons (must also account for shell-binding energy
777  // attributed to gamma from standard photoelectric effect)
778  else if(pDef==G4Electron::ElectronDefinition() ||
779  pDef==G4Gamma::GammaDefinition())
780  Scint_Yield_Vector = aMaterialPropertiesTable->
781  GetProperty("ELECTRONSCINTILLATIONYIELD");
782 
783  // Default for particles not enumerated/listed above
784  else
785  Scint_Yield_Vector = aMaterialPropertiesTable->
786  GetProperty("ELECTRONSCINTILLATIONYIELD");
787 
788  // If the user has specified none of the above particles then the
789  // default is the electron scintillation yield
790  if(!Scint_Yield_Vector)
791  Scint_Yield_Vector = aMaterialPropertiesTable->
792  GetProperty("ELECTRONSCINTILLATIONYIELD");
793 
794  // Throw an exception if no scintillation yield vector is found
795  if (!Scint_Yield_Vector) {
797  ed << "\nG4Scintillation::PostStepDoIt(): "
798  << "Request for scintillation yield for energy deposit and particle\n"
799  << "type without correct entry in MaterialPropertiesTable.\n"
800  << "ScintillationByParticleType requires at minimum that \n"
801  << "ELECTRONSCINTILLATIONYIELD is set by the user\n"
802  << G4endl;
803  G4String comments = "Missing MaterialPropertiesTable entry - No correct entry in MaterialPropertiesTable";
804  G4Exception("G4Scintillation::PostStepDoIt","Scint01",
805  FatalException,ed,comments);
806  }
807 
809  // Calculate the scintillation light //
811  // To account for potential nonlinearity and scintillation photon
812  // density along the track, light (L) is produced according to:
813  //
814  // L_currentStep = L(PreStepKE) - L(PreStepKE - EDep)
815 
816  G4double ScintillationYield = 0.;
817 
818  G4double StepEnergyDeposit = aStep.GetTotalEnergyDeposit();
819  G4double PreStepKineticEnergy = aStep.GetPreStepPoint()->GetKineticEnergy();
820 
821  if(PreStepKineticEnergy <= Scint_Yield_Vector->GetMaxEnergy()){
822  G4double Yield1 = Scint_Yield_Vector->Value(PreStepKineticEnergy);
823  G4double Yield2 = Scint_Yield_Vector->
824  Value(PreStepKineticEnergy - StepEnergyDeposit);
825  ScintillationYield = Yield1 - Yield2;
826  } else {
828  ed << "\nG4Scintillation::GetScintillationYieldByParticleType(): Request\n"
829  << "for scintillation light yield above the available energy range\n"
830  << "specifed in G4MaterialPropertiesTable. A linear interpolation\n"
831  << "will be performed to compute the scintillation light yield using\n"
832  << "(L_max / E_max) as the photon yield per unit energy."
833  << G4endl;
834  G4String cmt = "\nScintillation yield may be unphysical!\n";
835  G4Exception("G4Scintillation::GetScintillationYieldByParticleType()",
836  "Scint03", JustWarning, ed, cmt);
837 
838  G4double LinearYield = Scint_Yield_Vector->GetMaxValue()
839  / Scint_Yield_Vector->GetMaxEnergy();
840 
841  // Units: [# scintillation photons]
842  ScintillationYield = LinearYield * StepEnergyDeposit;
843  }
844 
845 #ifdef G4DEBUG_SCINTILLATION
846 
847  // Increment track aggregators
848  ScintTrackYield += ScintillationYield;
849  ScintTrackEDep += StepEnergyDeposit;
850 
851  G4cout << "\n--- G4Scintillation::GetScintillationYieldByParticleType() ---\n"
852  << "--\n"
853  << "-- Name = " << aTrack.GetParticleDefinition()->GetParticleName() << "\n"
854  << "-- TrackID = " << aTrack.GetTrackID() << "\n"
855  << "-- ParentID = " << aTrack.GetParentID() << "\n"
856  << "-- Current KE = " << aTrack.GetKineticEnergy()/MeV << " MeV\n"
857  << "-- Step EDep = " << aStep.GetTotalEnergyDeposit()/MeV << " MeV\n"
858  << "-- Track EDep = " << ScintTrackEDep/MeV << " MeV\n"
859  << "-- Vertex KE = " << aTrack.GetVertexKineticEnergy()/MeV << " MeV\n"
860  << "-- Step yield = " << ScintillationYield << " photons\n"
861  << "-- Track yield = " << ScintTrackYield << " photons\n"
862  << G4endl;
863 
864  // The track has terminated within or has left the scintillator volume
865  if( (aTrack.GetTrackStatus() == fStopButAlive) or
866  (aStep.GetPostStepPoint()->GetStepStatus() == fGeomBoundary) ){
867 
868  // Reset aggregators for the next track
869  ScintTrackEDep = 0.;
870  ScintTrackYield = 0.;
871  }
872 
873 #endif
874 
875  return ScintillationYield;
876 }
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)
static const double MeV
Definition: G4SIunits.hh:193
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4int GetNumberOfSecondaries() const
G4MaterialPropertyVector * GetProperty(const char *key)
G4double GetMaxEnergy() const
void SetFiniteRiseTime(const G4bool state)
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
G4double GetMeanLifeTime(const G4Track &aTrack, G4ForceCondition *)
size_t GetIndex() const
Definition: G4Material.hh:262
G4double single_exp(G4double t, G4double tau2)
G4PhysicsTable * fFastIntegralTable
G4StepStatus GetStepStatus() const
virtual G4double VisibleEnergyDeposition(const G4ParticleDefinition *, const G4MaterialCutsCouple *, G4double length, G4double edepTotal, G4double edepNIEL=0.0)
G4EmSaturation * fEmSaturation
G4TrackStatus GetTrackStatus() const
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:588
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
std::vector< G4Material * > G4MaterialTable
G4PhysicsTable * fSlowIntegralTable
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4VParticleChange * AtRestDoIt(const G4Track &aTrack, const G4Step &aStep)
G4ParticleDefinition * GetDefinition() const
G4double GetVelocity() const
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
G4StepPoint * GetPreStepPoint() const
void SetScintillationYieldFactor(const G4double yieldfactor)
G4double GetKineticEnergy() const
#define G4UniformRand()
Definition: Randomize.hh:95
G4GLOB_DLL std::ostream G4cout
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
G4double GetVertexKineticEnergy() const
const G4ThreeVector & GetPosition() const
void AddSaturation(G4EmSaturation *)
bool G4bool
Definition: G4Types.hh:79
G4double GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *)
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:595
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
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
void SetScintillationExcitationRatio(const G4double ratio)
G4double fExcitationRatio
static G4OpticalPhoton * OpticalPhoton()
virtual void Initialize(const G4Track &)
void SetTrackSecondariesFirst(const G4bool state)
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
G4double GetKineticEnergy() const
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4ForceCondition
void BuildPhysicsTable(const G4ParticleDefinition &aParticleType)
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:597
const G4TouchableHandle & GetTouchableHandle() const
void clearAndDestroy()
static G4Gamma * GammaDefinition()
Definition: G4Gamma.cc:81
G4ProcessType
G4bool fScintillationByParticleType