eic-opticks/u4/Local_DsG4Scintillation.cc

0001 //
0002 // ********************************************************************
0003 //  * DISCLAIMER                                                       *
0004 //  *                                                                  *
0005 //  * The following disclaimer summarizes all the specific disclaimers *
0006 //  * of contributors to this software. The specific disclaimers,which *
0007 //  * govern, are listed with their locations in:                      *
0008 //  *   http://cern.ch/geant4/license                                  *
0009 //  *                                                                  *
0010 //  * Neither the authors of this software system, nor their employing *
0011 //  * institutes,nor the agencies providing financial support for this *
0012 //  * work  make  any representation or  warranty, express or implied, *
0013 //  * regarding  this  software system or assume any liability for its *
0014 //  * use.                                                             *
0015 //  *                                                                  *
0016 //  * This  code  implementation is the  intellectual property  of the *
0017 //  * GEANT4 collaboration.                                            *
0018 //  * By copying,  distributing  or modifying the Program (or any work *
0019 //  * based  on  the Program)  you indicate  your  acceptance of  this *
0020 //  * statement, and all its terms.                                    *
0021 //  ********************************************************************
0022 //
0023 //
0024 //
0025 // //////////////////////////////////////////////////////////////////////
0026 //  Scintillation Light Class Implementation
0027 // //////////////////////////////////////////////////////////////////////
0028 //
0029 //  File:        G4Scintillation.cc
0030 //  Description: RestDiscrete Process - Generation of Scintillation Photons
0031 //  Version:     1.0
0032 //  Created:     1998-11-07
0033 //  Author:      Peter Gumplinger
0034 //  Updated:     2005-08-17 by Peter Gumplinger
0035 //               > change variable name MeanNumPhotons -> MeanNumberOfPhotons
0036 //               2005-07-28 by Peter Gumplinger
0037 //               > add G4ProcessType to constructor
0038 //               2004-08-05 by Peter Gumplinger
0039 //               > changed StronglyForced back to Forced in GetMeanLifeTime
0040 //               2002-11-21 by Peter Gumplinger
0041 //               > change to use G4Poisson for small MeanNumberOfPhotons
0042 //               2002-11-07 by Peter Gumplinger
0043 //               > now allow for fast and slow scintillation component
0044 //               2002-11-05 by Peter Gumplinger
0045 //               > now use scintillation constants from G4Material
0046 //               2002-05-09 by Peter Gumplinger
0047 //               > use only the PostStepPoint location for the origin of
0048 //                scintillation photons when energy is lost to the medium
0049 //                by a neutral particle
0050 //                2000-09-18 by Peter Gumplinger
0051 //               > change: aSecondaryPosition=x0+rand*aStep.GetDeltaPosition();
0052 //                aSecondaryTrack->SetTouchable(0);
0053 //                2001-09-17, migration of Materials to pure STL (mma)
0054 //                2003-06-03, V.Ivanchenko fix compilation warnings
0055 //
0056 //mail:        gum@triumf.ca
0057 //
0058 //////////////////////////////////////////////////////////////////////////
0059
0060 //-------------------------------------------------------------------
0061 // Local_DsG4Scintillation is a class modified from G4Scintillation
0062 // Birks' law is implemented
0063 // Author: Liang Zhan, 2006/01/27
0064 // Added weighted photon track method based on GLG4Scint. Jianglai 09/05/2006
0065 // Modified: bv@bnl.gov, 2008/4/16 for DetSim
0066 //--------------------------------------------------------------------
0067 //
0068 #ifdef STANDALONE
0069
0070 #include "SLOG.hh"
0071 #include "U4.hh"
0072 #else
0073 #include <boost/python.hpp>
0074 #endif
0075
0076 #include "Local_DsG4Scintillation.hh"
0077 #include "globals.hh"
0078 #include "G4PhysicalConstants.hh"
0079 #include "G4SystemOfUnits.hh"
0080 #include "G4UnitsTable.hh"
0081 #include "G4LossTableManager.hh"
0082 #include "G4MaterialCutsCouple.hh"
0083 #include "G4Gamma.hh"
0084 #include "G4Electron.hh"
0085 #include "globals.hh"
0086
0087 #ifdef WITH_G4OPTICKS
0088 #include "G4Opticks.hh"
0089 #include "CGenstep.hh"
0090 #include "CTrack.hh"
0091 #include "CPhotonInfo.hh"
0092 #include "SLOG.hh"
0093 #endif
0094
0095
0096
0097 #ifdef WITH_G4OPTICKS
0098 //const plog::Severity Local_DsG4Scintillation::LEVEL = SLOG::EnvLevel("Local_DsG4Scintillation", "DEBUG") ;
0099 const plog::Severity Local_DsG4Scintillation::LEVEL = error ;
0100 #endif
0101
0102
0103 #include "U4Stack.h"
0104 #include "SEvt.hh"
0105
0106 const bool Local_DsG4Scintillation::FLOAT = getenv("Local_DsG4Scintillation_FLOAT") != nullptr ;
0107 const int  Local_DsG4Scintillation::PIDX  = std::atoi( getenv("PIDX") ? getenv("PIDX") : "-1" );
0108
0109 //void Local_DsG4Scintillation::ResetNumberOfInteractionLengthLeft(){ G4VProcess::ResetNumberOfInteractionLengthLeft(); }
0110 void Local_DsG4Scintillation::ResetNumberOfInteractionLengthLeft()
0111 {
0112     //std::cout << "Local_DsG4Scintillation::FLOAT " << FLOAT << std::endl ;
0113     G4double u = G4UniformRand() ;
0114
0115 #ifndef PRODUCTION
0116 #ifdef DEBUG_TAG
0117     SEvt::AddTag( 1, U4Stack_ScintDiscreteReset, u );
0118 #endif
0119 #endif
0120
0121
0122
0123     if(FLOAT)
0124     {
0125         float f = -1.f*std::log( float(u) ) ;
0126         theNumberOfInteractionLengthLeft = f ;
0127     }
0128     else
0129     {
0130         theNumberOfInteractionLengthLeft = -1.*G4Log(u) ;
0131     }
0132     theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
0133 }
0134
0135
0136
0137
0138
0139
0140
0141 //#include "DsPhotonTrackInfo.h"
0142 //#include "G4DataHelpers/G4CompositeTrackInfo.h"
0143
0144 ///////////////////////////////////////////////////////////////////
0145
0146 using namespace std;
0147
0148 /////////////////////////
0149 // Class Implementation
0150 /////////////////////////
0151
0152 //////////////
0153 // Operators
0154 //////////////
0155
0156 // Local_DsG4Scintillation::operator=(const Local_DsG4Scintillation &right)
0157 // {
0158 // }
0159
0160 /////////////////
0161 // Constructors
0162 /////////////////
0163
0164 Local_DsG4Scintillation::Local_DsG4Scintillation(G4int opticksMode, const G4String& processName,
0165                                      G4ProcessType type)
0166     : G4VRestDiscreteProcess(processName, type)
0167     , doReemission(true)   // SCB set false to simplify debug
0168     , doBothProcess(true)
0169     , doReemissionOnly(false)
0170     , fEnableQuenching(true)
0171     , slowerTimeConstant(0) , slowerRatio(0)
0172     , gammaSlowerTime(0) , gammaSlowerRatio(0)
0173     , neutronSlowerTime(0) , neutronSlowerRatio(0)
0174     , alphaSlowerTime(0) , alphaSlowerRatio(0)
0175     , flagDecayTimeFast(true), flagDecayTimeSlow(true)
0176     , fPhotonWeight(1.0)
0177     , fApplyPreQE(false)
0178     , fPreQE(1.)
0179     , m_noop(false)
0180     , m_opticksMode(opticksMode)
0181 {
0182     SetProcessSubType(fScintillation);
0183     fTrackSecondariesFirst = false;
0184
0185     YieldFactor = 1.0;
0186     ExcitationRatio = 1.0;
0187
0188     theFastIntegralTable = NULL;
0189     theSlowIntegralTable = NULL;
0190     theReemissionIntegralTable = NULL;
0191
0192     //verboseLevel = 2;
0193     //G4cout << " Local_DsG4Scintillation set verboseLevel by hand to " << verboseLevel << G4endl;
0194
0195 #ifdef STANDALONE
0196     {
0197         const char* level_ = getenv("Local_DsG4Scintillation_verboseLevel") ;
0198         const char* fallback = "0" ;
0199         int level =  std::atoi(level_ ? level_ : fallback) ;
0200         SetVerboseLevel(level);
0201         std::cout << "Local_DsG4Scintillation::Local_DsG4Scintillation level " << level << " verboseLevel " << verboseLevel << std::endl ;
0202     }
0203 #endif
0204
0205     if (verboseLevel > 0) {
0206         G4cout << GetProcessName() << " is created " << G4endl;
0207     }
0208
0209     BuildThePhysicsTable();
0210
0211     // FORCE reemission only
0212     //doReemissionOnly = true;   // SCB commented
0213 }
0214
0215 ////////////////
0216 // Destructors
0217 ////////////////
0218
0219 Local_DsG4Scintillation::~Local_DsG4Scintillation()
0220 {
0221     if (theFastIntegralTable != NULL) {
0222         theFastIntegralTable->clearAndDestroy();
0223         delete theFastIntegralTable;
0224     }
0225     if (theSlowIntegralTable != NULL) {
0226         theSlowIntegralTable->clearAndDestroy();
0227         delete theSlowIntegralTable;
0228     }
0229     if (theReemissionIntegralTable != NULL) {
0230         theReemissionIntegralTable->clearAndDestroy();
0231         delete theReemissionIntegralTable;
0232     }
0233 }
0234
0235 ////////////
0236 // Methods
0237 ////////////
0238
0239 // AtRestDoIt
0240 // ----------
0241 //
0242 G4VParticleChange*
0243 Local_DsG4Scintillation::AtRestDoIt(const G4Track& aTrack, const G4Step& aStep)
0244
0245 // This routine simply calls the equivalent PostStepDoIt since all the
0246 // necessary information resides in aStep.GetTotalEnergyDeposit()
0247
0248 {
0249     return Local_DsG4Scintillation::PostStepDoIt(aTrack, aStep);
0250 }
0251
0252 // PostStepDoIt
0253 // -------------
0254 //
0255 G4VParticleChange*
0256 Local_DsG4Scintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
0257
0258 // This routine is called for each tracking step of a charged particle
0259 // in a scintillator. A Poisson/Gauss-distributed number of photons is
0260 // generated according to the scintillation yield formula, distributed
0261 // evenly along the track segment and uniformly into 4pi.
0262
0263 {
0264     aParticleChange.Initialize(aTrack);
0265     // SCB CHECK SEE notes/issues/Geant4_UseGivenVelocity_after_refraction_is_there_a_better_way_than_the_kludge_fix.rst
0266
0267     if (m_noop) {               // do nothing, bail
0268         aParticleChange.SetNumberOfSecondaries(0);
0269         return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
0270     }
0271
0272
0273     G4String pname="";
0274     G4ThreeVector vertpos;
0275     //G4double vertenergy=0.0;
0276     //G4double reem_d=0.0;
0277     G4bool flagReemission= false;
0278     //DsPhotonTrackInfo* reemittedTI=0;
0279     if (aTrack.GetDefinition() == G4OpticalPhoton::OpticalPhoton()) {
0280         G4Track *track=aStep.GetTrack();
0281         //G4CompositeTrackInfo* composite=dynamic_cast<G4CompositeTrackInfo*>(track->GetUserInformation());
0282         //reemittedTI = composite?dynamic_cast<DsPhotonTrackInfo*>( composite->GetPhotonTrackInfo() ):0;
0283
0284         const G4VProcess* process = track->GetCreatorProcess();
0285         if(process) pname = process->GetProcessName();
0286
0287         if (verboseLevel > 0) {
0288           G4cout<<"Optical photon. Process name is " << pname<<G4endl;
0289         }
0290         if(doBothProcess) {
0291             flagReemission= doReemission
0292                 && aTrack.GetTrackStatus() == fStopAndKill
0293                 && aStep.GetPostStepPoint()->GetStepStatus() != fGeomBoundary;
0294         }
0295         else{
0296             flagReemission= doReemission
0297                 && aTrack.GetTrackStatus() == fStopAndKill
0298                 && aStep.GetPostStepPoint()->GetStepStatus() != fGeomBoundary
0299                 && pname=="Cerenkov";
0300         }
0301         if(verboseLevel > 0) {
0302             G4cout<<"flag of Reemission is "<<flagReemission<<"!!"<<G4endl;
0303         }
0304         if (!flagReemission) {
0305             return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
0306         }
0307     }
0308
0309     G4double TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();
0310 #ifdef STANDALONE
0311     //LOG(info) << "TotalEnergyDeposit " << TotalEnergyDeposit ;
0312 #endif
0313
0314
0315     if (verboseLevel > 0 ) {
0316       G4cout << " TotalEnergyDeposit " << TotalEnergyDeposit
0317              << " material " << aTrack.GetMaterial()->GetName() << G4endl;
0318     }
0319     if (TotalEnergyDeposit <= 0.0 && !flagReemission) {
0320         return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
0321     }
0322
0323     const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
0324     const G4String aParticleName = aParticle->GetDefinition()->GetParticleName();
0325     const G4Material* aMaterial = aTrack.GetMaterial();
0326
0327     G4MaterialPropertiesTable* aMaterialPropertiesTable =
0328         aMaterial->GetMaterialPropertiesTable();
0329
0330     //aMaterialPropertiesTable-> DumpTable();
0331
0332     if (!aMaterialPropertiesTable)
0333         return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
0334
0335    // G4String FastTimeConstant = "FASTTIMECONSTANT";
0336    // G4String SlowTimeConstant = "SLOWTIMECONSTANT";
0337    // G4String strYieldRatio = "YIELDRATIO";
0338
0339     // reset the slower time constant and ratio
0340
0341     const G4MaterialPropertyVector* Fast_Intensity =
0342         aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
0343     const G4MaterialPropertyVector* Slow_Intensity =
0344         aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
0345     const G4MaterialPropertyVector* Reemission_Prob =
0346         aMaterialPropertiesTable->GetProperty("REEMISSIONPROB");
0347     if (verboseLevel > 0 ) {
0348       G4cout << " MaterialPropertyVectors: Fast_Intensity " << Fast_Intensity
0349              << " Slow_Intensity " << Slow_Intensity << " Reemission_Prob " << Reemission_Prob << G4endl;
0350     }
0351     if (!Fast_Intensity && !Slow_Intensity )
0352         return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
0353
0354     //-------------find the type of particle------------------------------//
0355     /*
0356         Find the particle type and register the scintillation time constant corresponding.
0357         We save the yield ratio and time constant in the form of G4PhysicVector. In this kind of G4PhysicVector, we interprete Energy as scintillation time and interprete Value as the yield ratio.
0358
0359     */
0360     G4MaterialPropertyVector* Ratio_timeconstant = 0 ;
0361     if (aParticleName == "opticalphoton") {
0362       Ratio_timeconstant = aMaterialPropertiesTable->GetProperty("OpticalCONSTANT");
0363     }
0364     else if(aParticleName == "gamma" || aParticleName == "e+" || aParticleName == "e-") {
0365       Ratio_timeconstant = aMaterialPropertiesTable->GetProperty("GammaCONSTANT");
0366     }
0367     else if(aParticleName == "alpha") {
0368       Ratio_timeconstant = aMaterialPropertiesTable->GetProperty("AlphaCONSTANT");
0369     }
0370     else {
0371       Ratio_timeconstant = aMaterialPropertiesTable->GetProperty("NeutronCONSTANT");
0372     }
0373
0374   //-----------------------------------------------------//
0375
0376     G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
0377     G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
0378
0379     G4ThreeVector x0 = pPreStepPoint->GetPosition();
0380     G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
0381     G4double      t0 = pPreStepPoint->GetGlobalTime();
0382
0383
0384     //Replace NumPhotons by NumTracks
0385     G4int NumTracks=0;
0386     G4double weight=1.0;
0387     if (flagReemission) {
0388         if(verboseLevel > 0){
0389             G4cout<<"the process name is "<<pname<<"!!"<<G4endl;}
0390
0391         if ( Reemission_Prob == 0)
0392             return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
0393         G4double p_reemission=
0394             Reemission_Prob->Value(aTrack.GetKineticEnergy());
0395
0396         G4double u_reemission = G4UniformRand() ;
0397
0398 #ifndef PRODUCTION
0399 #ifdef DEBUG_TAG
0400         SEvt::AddTag(1, U4Stack_Reemission, u_reemission);
0401 #endif
0402 #endif
0403
0404         if (u_reemission >= p_reemission)
0405             return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
0406         NumTracks= 1;
0407         weight= aTrack.GetWeight();
0408         if (verboseLevel > 0 ) {
0409             G4cout << " flagReemission " << flagReemission << " weight " << weight << G4endl;}
0410     }
0411     else {
0412         //////////////////////////////////// Birks' law ////////////////////////
0413         // J.B.Birks. The theory and practice of Scintillation Counting.
0414         // Pergamon Press, 1964.
0415         // For particles with energy much smaller than minimum ionization
0416         // energy, the scintillation response is non-linear because of quenching
0417         // effect. The light output is reduced by a parametric factor:
0418         // 1/(1 + birk1*delta + birk2* delta^2).
0419         // Delta is the energy loss per unit mass thickness. birk1 and birk2
0420         // were measured for several organic scintillators.
0421         // Here we use birk1 = 0.0125*g/cm2/MeV and ignore birk2.
0422         // R.L.Craun and D.L.Smith. Nucl. Inst. and Meth., 80:239-244, 1970.
0423         // Liang Zhan  01/27/2006
0424         // /////////////////////////////////////////////////////////////////////
0425
0426
0427         G4double ScintillationYield = 0;
0428         {// Yield.  Material must have this or we lack raisins dayetras
0429            /* const G4MaterialPropertyVector* ptable =
0430                 aMaterialPropertiesTable->GetProperty("SCINTILLATIONYIELD");
0431             if (!ptable) {
0432                 G4cout << "ConstProperty: failed to get SCINTILLATIONYIELD"
0433                        << G4endl;
0434                 return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
0435             }
0436             ScintillationYield = ptable->Value(0);
0437             std::cout<<"sci ScintillationYield = "<<ScintillationYield<<std::endl;*/
0438             ScintillationYield = aMaterialPropertiesTable->GetConstProperty("SCINTILLATIONYIELD");
0439            // std::cout<<"sci const ScintillationYield = "<<ScintillationYield<<std::endl;
0440         }
0441
0442         G4double ResolutionScale    = 1;
0443         {// Resolution Scale
0444             const G4MaterialPropertyVector* ptable =
0445                 aMaterialPropertiesTable->GetProperty("RESOLUTIONSCALE");
0446             if (ptable)
0447                 ResolutionScale = ptable->Value(0);
0448         }
0449
0450         G4double dE = TotalEnergyDeposit;
0451         G4double dx = aStep.GetStepLength();
0452         G4double dE_dx = dE/dx;
0453         if(aTrack.GetDefinition() == G4Gamma::Gamma() && dE > 0)
0454         {
0455           G4LossTableManager* manager = G4LossTableManager::Instance();
0456           dE_dx = dE/manager->GetRange(G4Electron::Electron(), dE, aTrack.GetMaterialCutsCouple());
0457           //G4cout<<"gamma dE_dx = "<<dE_dx/(MeV/mm)<<"MeV/mm"<<G4endl;
0458         }
0459
0460         G4double delta = dE_dx/aMaterial->GetDensity();//get scintillator density
0461         //G4double birk1 = 0.0125*g/cm2/MeV;
0462         G4double birk1 = birksConstant1;
0463         if(abs(aParticle->GetCharge())>1.5)//for particle charge greater than 1.
0464             birk1 = 0.57*birk1;
0465
0466         G4double birk2 = 0;
0467         //birk2 = (0.0031*g/MeV/cm2)*(0.0031*g/MeV/cm2);
0468         birk2 = birksConstant2;
0469
0470         G4double QuenchedTotalEnergyDeposit = TotalEnergyDeposit;
0471         // if quenching is enabled, apply the birks law
0472         if (fEnableQuenching) {
0473             QuenchedTotalEnergyDeposit
0474             = TotalEnergyDeposit/(1+birk1*delta+birk2*delta*delta);
0475         }
0476
0477        //Add 300ns trick for muon simuation, by Haoqi Jan 27, 2011
0478        if(FastMu300nsTrick)  {
0479            // cout<<"GlobalTime ="<<aStep.GetTrack()->GetGlobalTime()/ns<<endl;
0480            if(aStep.GetTrack()->GetGlobalTime()/ns>300) {
0481                ScintillationYield = YieldFactor * ScintillationYield;
0482            }
0483            else{
0484             ScintillationYield=0.;
0485            }
0486         }
0487         else {
0488             ScintillationYield = YieldFactor * ScintillationYield;
0489         }
0490
0491         G4double MeanNumberOfPhotons= ScintillationYield * QuenchedTotalEnergyDeposit;
0492
0493         // Implemented the fast simulation method from GLG4Scint
0494         // Jianglai 09-05-2006
0495
0496         // randomize number of TRACKS (not photons)
0497         // this gets statistics right for number of PE after applying
0498         // boolean random choice to final absorbed track (change from
0499         // old method of applying binomial random choice to final absorbed
0500         // track, which did want poissonian number of photons divided
0501         // as evenly as possible into tracks)
0502         // Note for weight=1, there's no difference between tracks and photons.
0503         G4double MeanNumberOfTracks= MeanNumberOfPhotons/fPhotonWeight;
0504         if ( fApplyPreQE ) {
0505             MeanNumberOfTracks*=fPreQE;
0506         }
0507         if (MeanNumberOfTracks > 10.) {
0508             G4double sigma = ResolutionScale * sqrt(MeanNumberOfTracks);
0509             NumTracks = G4int(G4RandGauss::shoot(MeanNumberOfTracks,sigma)+0.5);
0510         }
0511         else {
0512             NumTracks = G4int(G4Poisson(MeanNumberOfTracks));
0513         }
0514         if ( verboseLevel > 0 ) {
0515           G4cout << " Generated " << NumTracks << " scint photons. mean(scint photons) = " << MeanNumberOfTracks << G4endl;
0516         }
0517     }
0518
0519     weight*=fPhotonWeight;
0520     if ( verboseLevel > 0 ) {
0521       G4cout << " set scint photon weight to " << weight << " after multiplying original weight by fPhotonWeight " << fPhotonWeight
0522              << " NumTracks = " << NumTracks
0523              << G4endl;
0524     }
0525     // G4cerr<<"Scint weight is "<<weight<<G4endl;
0526     if (NumTracks <= 0) {
0527         // return unchanged particle and no secondaries
0528         aParticleChange.SetNumberOfSecondaries(0);
0529         return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
0530     }
0531
0532     ////////////////////////////////////////////////////////////////
0533
0534     aParticleChange.SetNumberOfSecondaries(NumTracks);
0535
0536     if (fTrackSecondariesFirst) {
0537         if (!flagReemission)
0538             if (aTrack.GetTrackStatus() == fAlive )
0539                 aParticleChange.ProposeTrackStatus(fSuspend);
0540     }
0541
0542     ////////////////////////////////////////////////////////////////
0543
0544     G4int materialIndex = aMaterial->GetIndex();
0545
0546     //G4PhysicsOrderedFreeVector* ReemissionIntegral = NULL;
0547     //ReemissionIntegral =
0548     //    (G4PhysicsOrderedFreeVector*)((*theReemissionIntegralTable)(materialIndex));
0549
0550     // Retrieve the Scintillation Integral for this material
0551     // new G4PhysicsOrderedFreeVector allocated to hold CII's
0552
0553    // G4int Num = NumTracks; //# tracks is now the loop control
0554
0555    /*
0556       Determine the photon number of echo component at first.
0557       we determine the photon number by Monte Carlo sampling.
0558       reason:
0559       It may lose the little photons if we just use total number times yield ratio when total number is small
0560     */
0561     int nscnt = Ratio_timeconstant->GetVectorLength();
0562     std::vector<G4int> NumVec(nscnt);
0563     NumVec.clear();
0564     for(G4int i = 0 ; i < NumTracks ; i++){
0565        G4double p = G4UniformRand();
0566
0567 #ifndef PRODUCTION
0568 #ifdef DEBUG_TAG
0569        SEvt::AddTag(1, U4Stack_Reemission, p );
0570 #endif
0571 #endif
0572
0573        G4double p_count = 0;
0574        for(G4int j = 0 ; j < nscnt; j++)
0575        {
0576             p_count += (*Ratio_timeconstant)[j] ;
0577             if( p < p_count ){
0578                NumVec[j]++ ;
0579                break;
0580             }
0581         }
0582
0583      }
0584
0585 //-------------------------------------------------//
0586
0587 #ifdef WITH_G4OPTICKS
0588
0589     CPho ancestor = CPhotonInfo::Get(&aTrack);
0590     int ancestor_id = ancestor.get_id() ;
0591
0592
0593     /**
0594     ancestor_id:-1
0595         normal case, meaning that aTrack was not a photon
0596         so the generation loop will yield "primary" photons
0597         with id : gs.offset + i
0598
0599     ancestor_id>-1
0600         aTrack is a photon that may be about to reemit (Num=0 or 1)
0601         ancestor_id is the absolute id of the primary parent photon,
0602         this id is retained thru any subsequent remission secondary generations
0603     **/
0604
0605 #endif
0606
0607 #ifdef STANDALONE
0608     U4::GenPhotonAncestor(&aTrack);
0609 #endif
0610
0611
0612 //-------------------------------------------------//
0613
0614     for(G4int scnt = 0 ; scnt < nscnt ; scnt++){
0615
0616          G4double ScintillationTime = 0.*ns;
0617          G4PhysicsOrderedFreeVector* ScintillationIntegral = NULL;
0618
0619          if ( scnt == 0 ){
0620               ScintillationIntegral =
0621                     (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
0622          }
0623          else{
0624               ScintillationIntegral =
0625                     (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
0626          }
0627        // std::cout<<"scnt == "<<scnt <<"  Num =="<<NumVec[scnt]<<" time =="<<Ratio_timeconstant->Energy(scnt)<<std::endl;
0628
0629        //  G4int m_Num =G4int(NumTracks * (*Ratio_timeconstant)[scnt]);
0630          ScintillationTime = Ratio_timeconstant->Energy(scnt);
0631          if (!flagDecayTimeFast && scnt == 0){
0632                ScintillationTime = 0.*ns  ;
0633          }
0634
0635          if (!flagDecayTimeSlow && scnt != 0){
0636
0637                ScintillationTime = 0.*ns  ;
0638          }
0639
0640         G4int NumPhoton =  NumVec[scnt] ;
0641
0642
0643 #ifdef WITH_G4OPTICKS
0644         if(flagReemission) assert( NumPhoton == 0 || NumPhoton == 1);   // expecting only 0 or 1 remission photons
0645         bool is_opticks_genstep = NumPhoton > 0 && !flagReemission ;
0646
0647         CGenstep gs ;
0648         if(is_opticks_genstep && (m_opticksMode & 1))
0649         {
0650             gs = G4Opticks::Get()->collectGenstep_Local_DsG4Scintillation_r4695( &aTrack, &aStep, NumPhoton, scnt, ScintillationTime);
0651         }
0652 #endif
0653
0654 #ifdef STANDALONE
0655         if(flagReemission) assert( NumPhoton == 0 || NumPhoton == 1);   // expecting only 0 or 1 remission photons
0656         bool is_opticks_genstep = NumPhoton > 0 && !flagReemission ;
0657         if(is_opticks_genstep && (m_opticksMode & 1))
0658         {
0659             NumPhoton = std::min( NumPhoton, 3 );  // for debugging purposes it helps to have less photons
0660             U4::CollectGenstep_DsG4Scintillation_r4695( &aTrack, &aStep, NumPhoton, scnt, ScintillationTime);
0661         }
0662 #endif
0663
0664          if( m_opticksMode == 0 || (m_opticksMode & 2) )
0665          {
0666
0667          for(G4int i = 0 ; i < NumPhoton ; i++) {
0668 #ifdef WITH_G4OPTICKS
0669            G4Opticks::Get()->setAlignIndex( ancestor_id > -1 ? ancestor_id : gs.offset + i );  // aka photon_id
0670 #endif
0671 #ifdef STANDALONE
0672            U4::GenPhotonBegin(i);
0673 #endif
0674
0675            G4double sampledEnergy;
0676            if ( !flagReemission ) {
0677                 // normal scintillation
0678
0679                G4double u_sampledEnergy = G4UniformRand() ;
0680
0681                G4double CIIvalue = u_sampledEnergy*
0682                     ScintillationIntegral->GetMaxValue();
0683                sampledEnergy=
0684                     ScintillationIntegral->GetEnergy(CIIvalue);
0685
0686                if (verboseLevel>1)
0687                     {
0688                         G4cout << "sampledEnergy = " << sampledEnergy << G4endl;
0689                         G4cout << "CIIvalue =        " << CIIvalue << G4endl;
0690                     }
0691             }
0692          else {
0693                 // reemission, the sample method need modification
0694
0695                 G4double u_sampledEnergy = G4UniformRand() ;
0696 #ifndef PRODUCTION
0697 #ifdef DEBUG_TAG
0698                 SEvt::AddTag(1, U4Stack_Reemission, u_sampledEnergy );
0699 #endif
0700 #endif
0701
0702                 G4double CIIvalue = u_sampledEnergy*
0703                     ScintillationIntegral->GetMaxValue();
0704                 if (CIIvalue == 0.0) {
0705                     // return unchanged particle and no secondaries
0706                     aParticleChange.SetNumberOfSecondaries(0);
0707                     return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
0708                    }
0709                 sampledEnergy=
0710                     ScintillationIntegral->GetEnergy(CIIvalue);
0711                 if (verboseLevel>1) {
0712                     G4cout << "oldEnergy = " <<aTrack.GetKineticEnergy() << G4endl;
0713                     G4cout << "reemittedSampledEnergy = " << sampledEnergy
0714                            << "\nreemittedCIIvalue =        " << CIIvalue << G4endl;
0715                    }
0716              }
0717
0718            // Generate random photon direction
0719
0720             G4double u_cost = G4UniformRand();
0721 #ifndef PRODUCTION
0722 #ifdef DEBUG_TAG
0723             SEvt::AddTag(1, U4Stack_Reemission, u_cost );
0724 #endif
0725 #endif
0726
0727             G4double cost = 1. - 2.*u_cost ;
0728             G4double sint = sqrt((1.-cost)*(1.+cost));
0729
0730             G4double u_phi = G4UniformRand() ;
0731 #ifndef PRODUCTION
0732 #ifdef DEBUG_TAG
0733             SEvt::AddTag(1, U4Stack_Reemission, u_phi );
0734 #endif
0735 #endif
0736
0737             G4double phi = twopi*u_phi ;
0738             G4double sinp = sin(phi);
0739             G4double cosp = cos(phi);
0740
0741             G4double px = sint*cosp;
0742             G4double py = sint*sinp;
0743             G4double pz = cost;
0744
0745             // Create photon momentum direction vector
0746             G4ParticleMomentum photonMomentum(px, py, pz);
0747
0748             // Determine polarization of new photon
0749             G4double sx = cost*cosp;
0750             G4double sy = cost*sinp;
0751             G4double sz = -sint;
0752
0753             G4ThreeVector photonPolarization(sx, sy, sz);
0754
0755             G4ThreeVector perp = photonMomentum.cross(photonPolarization);
0756
0757             u_phi = G4UniformRand() ;
0758 #ifndef PRODUCTION
0759 #ifdef DEBUG_TAG
0760             SEvt::AddTag(1, U4Stack_Reemission, u_phi );
0761 #endif
0762 #endif
0763
0764             phi = twopi*u_phi ;
0765             sinp = sin(phi);
0766             cosp = cos(phi);
0767
0768             photonPolarization = cosp * photonPolarization + sinp * perp;
0769
0770             photonPolarization = photonPolarization.unit();
0771
0772             // Generate a new photon:
0773
0774             G4DynamicParticle* aScintillationPhoton =
0775                 new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(),
0776                                       photonMomentum);
0777             aScintillationPhoton->SetPolarization
0778                 (photonPolarization.x(),
0779                  photonPolarization.y(),
0780                  photonPolarization.z());
0781
0782             aScintillationPhoton->SetKineticEnergy(sampledEnergy);
0783
0784             // Generate new G4Track object:
0785             G4double rand=0;
0786             G4ThreeVector aSecondaryPosition;
0787             G4double deltaTime;
0788             if (flagReemission) {
0789
0790
0791                 G4double u_deltaTime = G4UniformRand() ;
0792 #ifndef PRODUCTION
0793 #ifdef DEBUG_TAG
0794                 SEvt::AddTag(1, U4Stack_Reemission, u_deltaTime );
0795 #endif
0796 #endif
0797
0798                 deltaTime= pPostStepPoint->GetGlobalTime() - t0
0799                            -ScintillationTime * log( u_deltaTime );
0800                 aSecondaryPosition= pPostStepPoint->GetPosition();
0801                 vertpos = aTrack.GetVertexPosition();
0802                 //vertenergy = aTrack.GetKineticEnergy();
0803                 //reem_d =
0804                 //    sqrt( pow( aSecondaryPosition.x()-vertpos.x(), 2)
0805                 //          + pow( aSecondaryPosition.y()-vertpos.y(), 2)
0806                 //          + pow( aSecondaryPosition.z()-vertpos.z(), 2) );
0807             }
0808             else {
0809                 if (aParticle->GetDefinition()->GetPDGCharge() != 0)
0810                     {
0811                         rand = G4UniformRand();
0812                     }
0813                 else
0814                     {
0815                         rand = 1.0;
0816                     }
0817
0818                 G4double delta = rand * aStep.GetStepLength();
0819                 deltaTime = delta /
0820                     ((pPreStepPoint->GetVelocity()+
0821                       pPostStepPoint->GetVelocity())/2.);
0822
0823                 deltaTime = deltaTime -
0824                     ScintillationTime * log( G4UniformRand() );
0825
0826                 aSecondaryPosition =
0827                     x0 + rand * aStep.GetDeltaPosition();
0828             }
0829             G4double aSecondaryTime = t0 + deltaTime;
0830             if ( verboseLevel>1 ){
0831               G4cout << "Generate " << i << "th scintillation photon at relative time(ns) " << deltaTime
0832                      << " with ScintillationTime " << ScintillationTime << " flagReemission " << flagReemission << G4endl;
0833             }
0834             G4Track* aSecondaryTrack =
0835                 new G4Track(aScintillationPhoton,aSecondaryTime,aSecondaryPosition);
0836
0837             aSecondaryTrack->SetWeight( weight );
0838             aSecondaryTrack->SetTouchableHandle(aStep.GetPreStepPoint()->GetTouchableHandle());
0839             aSecondaryTrack->SetParentID(aTrack.GetTrackID());
0840             // add the secondary to the ParticleChange object
0841             aParticleChange.SetSecondaryWeightByProcess( true ); // recommended
0842             aParticleChange.AddSecondary(aSecondaryTrack);
0843
0844             aSecondaryTrack->SetWeight( weight );
0845             if ( verboseLevel > 0 ) {
0846               G4cout << " aSecondaryTrack->SetWeight( " << weight<< " ) ; aSecondaryTrack->GetWeight() = " << aSecondaryTrack->GetWeight() << G4endl;}
0847
0848 #ifdef WITH_G4OPTICKS
0849            aSecondaryTrack->SetUserInformation(CPhotonInfo::MakeScintillation(gs, i, ancestor ));
0850            G4Opticks::Get()->setAlignIndex(-1);
0851 #endif
0852
0853 #ifdef STANDALONE
0854            U4::GenPhotonEnd(i, aSecondaryTrack);
0855 #endif
0856
0857          }    // i:genloop over NumPhoton
0858
0859
0860          }   //  (opticksMode == 0) || (opticksMode & 2 )   : opticks not enabled, or opticks enabled and doing Geant4 comparison
0861
0862
0863    }         // scntloop
0864
0865
0866
0867     G4VParticleChange* change = G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
0868
0869
0870     if (verboseLevel > 0) {
0871         G4cout << "\n Exiting from G4Scintillation::DoIt -- NumberOfSecondaries = "
0872                << aParticleChange.GetNumberOfSecondaries() << G4endl;
0873     }
0874
0875     return change ;
0876 }
0877
0878
0879
0880 #ifdef WITH_G4OPTICKS
0881 G4MaterialPropertyVector* Local_DsG4Scintillation::getMaterialProperty(const char* name, G4int materialIndex)
0882 {
0883     const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
0884     G4int numOfMaterials = G4Material::GetNumberOfMaterials();
0885     assert( materialIndex < numOfMaterials );
0886
0887     G4Material* aMaterial = (*theMaterialTable)[materialIndex];
0888     G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
0889     G4MaterialPropertyVector* prop = aMaterialPropertiesTable ? aMaterialPropertiesTable->GetProperty(name) : nullptr ;
0890     return prop ;
0891 }
0892
0893 G4PhysicsOrderedFreeVector* Local_DsG4Scintillation::getScintillationIntegral(G4int scnt, G4int materialIndex) const
0894 {
0895     G4PhysicsOrderedFreeVector* ScintillationIntegral = NULL;
0896
0897     if ( scnt == 0 ){
0898           ScintillationIntegral =
0899                 (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
0900     }
0901     else{
0902           ScintillationIntegral =
0903                 (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
0904     }
0905
0906     return ScintillationIntegral ;
0907 }
0908
0909
0910 G4double Local_DsG4Scintillation::getSampledEnergy(G4int scnt, G4int materialIndex) const
0911 {
0912     G4PhysicsOrderedFreeVector* ScintillationIntegral = getScintillationIntegral(scnt, materialIndex);
0913     G4double CIIvalue = G4UniformRand()*ScintillationIntegral->GetMaxValue();
0914     G4double sampledEnergy = ScintillationIntegral->GetEnergy(CIIvalue);
0915     return sampledEnergy ;
0916 }
0917
0918 G4double Local_DsG4Scintillation::getSampledWavelength(G4int scnt, G4int materialIndex) const
0919 {
0920     G4double sampledEnergy = getSampledEnergy(scnt, materialIndex );
0921     G4double wavelength = h_Planck*c_light/sampledEnergy ;
0922     G4double wavelength_nm = wavelength/nm ;
0923     return wavelength_nm ;
0924 }
0925 #endif
0926
0927
0928 // BuildThePhysicsTable for the scintillation process
0929 // --------------------------------------------------
0930 //
0931
0932 void Local_DsG4Scintillation::BuildThePhysicsTable()
0933 {
0934     if (theFastIntegralTable && theSlowIntegralTable && theReemissionIntegralTable) return;
0935
0936     const G4MaterialTable* theMaterialTable =
0937         G4Material::GetMaterialTable();
0938     G4int numOfMaterials = G4Material::GetNumberOfMaterials();
0939
0940     // create new physics table
0941     if (verboseLevel > 0) {
0942       G4cout << " theFastIntegralTable " << theFastIntegralTable
0943              << " theSlowIntegralTable " << theSlowIntegralTable
0944              << " theReemissionIntegralTable " << theReemissionIntegralTable << G4endl;
0945     }
0946     if(!theFastIntegralTable)theFastIntegralTable = new G4PhysicsTable(numOfMaterials);
0947     if(!theSlowIntegralTable)theSlowIntegralTable = new G4PhysicsTable(numOfMaterials);
0948     if(!theReemissionIntegralTable)theReemissionIntegralTable
0949                                        = new G4PhysicsTable(numOfMaterials);
0950     if (verboseLevel > 0) {
0951       G4cout << " building the physics tables for the scintillation process " << G4endl;
0952     }
0953     // loop for materials
0954
0955     for (G4int i=0 ; i < numOfMaterials; i++) {
0956         G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
0957             new G4PhysicsOrderedFreeVector();
0958         G4PhysicsOrderedFreeVector* bPhysicsOrderedFreeVector =
0959             new G4PhysicsOrderedFreeVector();
0960         G4PhysicsOrderedFreeVector* cPhysicsOrderedFreeVector =
0961             new G4PhysicsOrderedFreeVector();
0962
0963         // Retrieve vector of scintillation wavelength intensity for
0964         // the material from the material's optical properties table.
0965
0966         G4Material* aMaterial = (*theMaterialTable)[i];
0967
0968         G4MaterialPropertiesTable* aMaterialPropertiesTable =
0969             aMaterial->GetMaterialPropertiesTable();
0970
0971         if (aMaterialPropertiesTable) {
0972
0973             G4MaterialPropertyVector* theFastLightVector =
0974                 aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
0975
0976             if (theFastLightVector) {
0977               if (verboseLevel > 0) {
0978                 G4cout << " Building the material properties table for FASTCOMPONENT" << G4endl;
0979               }
0980                 // Retrieve the first intensity point in vector
0981                 // of (photon energy, intensity) pairs
0982
0983                 G4double currentIN = (*theFastLightVector)[0];
0984
0985                 if (currentIN >= 0.0) {
0986
0987                     // Create first (photon energy, Scintillation
0988                     // Integral pair
0989
0990                     G4double currentPM = theFastLightVector->
0991                         Energy(0);
0992
0993                     G4double currentCII = 0.0;
0994
0995                     aPhysicsOrderedFreeVector->
0996                         InsertValues(currentPM , currentCII);
0997
0998                     // Set previous values to current ones prior to loop
0999
1000                     G4double prevPM  = currentPM;
1001                     G4double prevCII = currentCII;
1002                     G4double prevIN  = currentIN;
1003
1004                     // loop over all (photon energy, intensity)
1005                     // pairs stored for this material
1006
1007                     for(size_t ii = 1;
1008                               ii < theFastLightVector->GetVectorLength();
1009                               ++ii)
1010                     {
1011                         currentPM = theFastLightVector->Energy(ii);
1012
1013                         currentIN= (*theFastLightVector)[ii];
1014
1015                         currentCII = 0.5 * (prevIN + currentIN);
1016
1017                         currentCII = prevCII +
1018                             (currentPM - prevPM) * currentCII;
1019
1020                         aPhysicsOrderedFreeVector->
1021                             InsertValues(currentPM, currentCII);
1022
1023                         prevPM  = currentPM;
1024                         prevCII = currentCII;
1025                         prevIN  = currentIN;
1026                     }
1027
1028                 }
1029             }
1030
1031             G4MaterialPropertyVector* theSlowLightVector =
1032                 aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
1033
1034             if (theSlowLightVector) {
1035                 if (verboseLevel > 0) {
1036                   G4cout << " Building the material properties table for SLOWCOMPONENT" << G4endl;
1037                 }
1038
1039                 // Retrieve the first intensity point in vector
1040                 // of (photon energy, intensity) pairs
1041
1042                 G4double currentIN = (*theSlowLightVector)[0];
1043
1044                 if (currentIN >= 0.0) {
1045
1046                     // Create first (photon energy, Scintillation
1047                     // Integral pair
1048
1049                     G4double currentPM = theSlowLightVector->Energy(0);
1050
1051                     G4double currentCII = 0.0;
1052
1053                     bPhysicsOrderedFreeVector->
1054                         InsertValues(currentPM , currentCII);
1055
1056                     // Set previous values to current ones prior to loop
1057
1058                     G4double prevPM  = currentPM;
1059                     G4double prevCII = currentCII;
1060                     G4double prevIN  = currentIN;
1061
1062                     // loop over all (photon energy, intensity)
1063                     // pairs stored for this material
1064
1065                     for (size_t ii = 1;
1066                          ii < theSlowLightVector->GetVectorLength();
1067                          ++ii)
1068                     {
1069                         currentPM = theSlowLightVector->Energy(ii);
1070
1071                         currentIN = (*theSlowLightVector)[ii];
1072
1073                         currentCII = 0.5 * (prevIN + currentIN);
1074
1075                         currentCII = prevCII +
1076                             (currentPM - prevPM) * currentCII;
1077
1078                         bPhysicsOrderedFreeVector->
1079                             InsertValues(currentPM, currentCII);
1080
1081                         prevPM  = currentPM;
1082                         prevCII = currentCII;
1083                         prevIN  = currentIN;
1084                     }
1085
1086                 }
1087             }
1088
1089             G4MaterialPropertyVector* theReemissionVector =
1090                 aMaterialPropertiesTable->GetProperty("REEMISSIONPROB");
1091
1092             if (theReemissionVector) {
1093               if (verboseLevel > 0) {
1094                 G4cout << " Building the material properties table for REEMISSIONPROB" << G4endl;
1095               }
1096
1097                 // Retrieve the first intensity point in vector
1098                 // of (photon energy, intensity) pairs
1099
1100               G4double currentIN = (*theReemissionVector)[0];
1101
1102                 if (currentIN >= 0.0) {
1103
1104                     // Create first (photon energy, Scintillation
1105                     // Integral pair
1106
1107                     G4double currentPM = theReemissionVector->Energy(0);
1108
1109                     G4double currentCII = 0.0;
1110
1111                     cPhysicsOrderedFreeVector->
1112                         InsertValues(currentPM , currentCII);
1113
1114                     // Set previous values to current ones prior to loop
1115
1116                     G4double prevPM  = currentPM;
1117                     G4double prevCII = currentCII;
1118                     G4double prevIN  = currentIN;
1119
1120                     // loop over all (photon energy, intensity)
1121                     // pairs stored for this material
1122
1123                     for (size_t ii = 1;
1124                          ii < theReemissionVector->GetVectorLength();
1125                          ++ii)
1126                     {
1127
1128                         currentPM = theReemissionVector->Energy(ii);
1129
1130                         currentIN = (*theReemissionVector)[ii];
1131
1132                         currentCII = 0.5 * (prevIN + currentIN);
1133
1134                         currentCII = prevCII +
1135                             (currentPM - prevPM) * currentCII;
1136
1137                         cPhysicsOrderedFreeVector->
1138                             InsertValues(currentPM, currentCII);
1139
1140                         prevPM  = currentPM;
1141                         prevCII = currentCII;
1142                         prevIN  = currentIN;
1143                     }
1144
1145                 }
1146             }
1147
1148         }
1149
1150         // The scintillation integral(s) for a given material
1151         // will be inserted in the table(s) according to the
1152         // position of the material in the material table.
1153
1154         theFastIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
1155         theSlowIntegralTable->insertAt(i,bPhysicsOrderedFreeVector);
1156         theReemissionIntegralTable->insertAt(i,cPhysicsOrderedFreeVector);
1157     }
1158
1159     if (verboseLevel > 0) {
1160       G4cout << " theFastIntegralTable " << theFastIntegralTable
1161              << " theSlowIntegralTable " << theSlowIntegralTable
1162              << " theReemissionIntegralTable " << theReemissionIntegralTable << G4endl;
1163     }
1164
1165
1166 }
1167
1168 // GetMeanFreePath
1169 // ---------------
1170 //
1171
1172 G4double Local_DsG4Scintillation::GetMeanFreePath(const G4Track&,
1173                                             G4double ,
1174                                             G4ForceCondition* condition)
1175 {
1176     *condition = StronglyForced;
1177
1178     return DBL_MAX;
1179
1180 }
1181
1182 // GetMeanLifeTime
1183 // ---------------
1184 //
1185
1186 G4double Local_DsG4Scintillation::GetMeanLifeTime(const G4Track&,
1187                                             G4ForceCondition* condition)
1188 {
1189     *condition = Forced;
1190
1191     return DBL_MAX;
1192
1193 }
1194
1195 // OP simulator
1196 G4PhysicsTable* Local_DsG4Scintillation::getSlowIntegralTable() {
1197     return theSlowIntegralTable;
1198 }
1199 G4PhysicsTable* Local_DsG4Scintillation::getFastIntegralTable() {
1200     return theFastIntegralTable;
1201 }
1202 G4PhysicsTable* Local_DsG4Scintillation::getReemissionIntegralTable() {
1203     return theReemissionIntegralTable;
1204 }