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