/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

DsG4Scintillation.cc
Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 //  * DISCLAIMER                                                       *
00004 //  *                                                                  *
00005 //  * The following disclaimer summarizes all the specific disclaimers *
00006 //  * of contributors to this software. The specific disclaimers,which *
00007 //  * govern, are listed with their locations in:                      *
00008 //  *   http://cern.ch/geant4/license                                  *
00009 //  *                                                                  *
00010 //  * Neither the authors of this software system, nor their employing *
00011 //  * institutes,nor the agencies providing financial support for this *
00012 //  * work  make  any representation or  warranty, express or implied, *
00013 //  * regarding  this  software system or assume any liability for its *
00014 //  * use.                                                             *
00015 //  *                                                                  *
00016 //  * This  code  implementation is the  intellectual property  of the *
00017 //  * GEANT4 collaboration.                                            *
00018 //  * By copying,  distributing  or modifying the Program (or any work *
00019 //  * based  on  the Program)  you indicate  your  acceptance of  this *
00020 //  * statement, and all its terms.                                    *
00021 //  ********************************************************************
00022 // 
00023 // 
00024 // 
00025 // //////////////////////////////////////////////////////////////////////
00026 //  Scintillation Light Class Implementation
00027 // //////////////////////////////////////////////////////////////////////
00028 // 
00029 //  File:        G4Scintillation.cc 
00030 //  Description: RestDiscrete Process - Generation of Scintillation Photons
00031 //  Version:     1.0
00032 //  Created:     1998-11-07  
00033 //  Author:      Peter Gumplinger
00034 //  Updated:     2005-08-17 by Peter Gumplinger
00035 //               > change variable name MeanNumPhotons -> MeanNumberOfPhotons
00036 //               2005-07-28 by Peter Gumplinger
00037 //               > add G4ProcessType to constructor
00038 //               2004-08-05 by Peter Gumplinger
00039 //               > changed StronglyForced back to Forced in GetMeanLifeTime
00040 //               2002-11-21 by Peter Gumplinger
00041 //               > change to use G4Poisson for small MeanNumberOfPhotons
00042 //               2002-11-07 by Peter Gumplinger
00043 //               > now allow for fast and slow scintillation component
00044 //               2002-11-05 by Peter Gumplinger
00045 //               > now use scintillation constants from G4Material
00046 //               2002-05-09 by Peter Gumplinger
00047 //               > use only the PostStepPoint location for the origin of
00048 //                scintillation photons when energy is lost to the medium
00049 //                by a neutral particle
00050 //                2000-09-18 by Peter Gumplinger
00051 //               > change: aSecondaryPosition=x0+rand*aStep.GetDeltaPosition();
00052 //                aSecondaryTrack->SetTouchable(0);
00053 //                2001-09-17, migration of Materials to pure STL (mma) 
00054 //                2003-06-03, V.Ivanchenko fix compilation warnings
00055 //    
00056 //mail:        gum@triumf.ca
00057 //               
00059 
00060 //-------------------------------------------------------------------
00061 // DsG4Scintillation is a class modified from G4Scintillation
00062 // Birks' law is implemented 
00063 // Author: Liang Zhan, 2006/01/27
00064 // Added weighted photon track method based on GLG4Scint. Jianglai 09/05/2006
00065 // Modified: bv@bnl.gov, 2008/4/16 for DetSim
00066 //--------------------------------------------------------------------
00067 
00068 #include "DsG4Scintillation.h"
00069 #include "G4UnitsTable.hh"
00070 #include "G4LossTableManager.hh"
00071 #include "G4MaterialCutsCouple.hh"
00072 #include "G4Gamma.hh"
00073 #include "G4Electron.hh"
00074 #include "globals.hh"
00075 
00076 #include "DsPhotonTrackInfo.h"
00077 #include "G4DataHelpers/G4CompositeTrackInfo.h"
00078 
00080 
00081 using namespace std;
00082 
00084 // Class Implementation  
00086 
00088 // Operators
00090 
00091 // DsG4Scintillation::operator=(const DsG4Scintillation &right)
00092 // {
00093 // }
00094 
00096 // Constructors
00098 
00099 DsG4Scintillation::DsG4Scintillation(const G4String& processName,
00100                                      G4ProcessType type)
00101     : G4VRestDiscreteProcess(processName, type)
00102     , doReemission(true)
00103     , doBothProcess(true)
00104     , fPhotonWeight(1.0)
00105     , fApplyPreQE(false)
00106     , fPreQE(1.)
00107     , m_noop(false)
00108 {
00109     SetProcessSubType(fScintillation);
00110     fTrackSecondariesFirst = false;
00111 
00112     YieldFactor = 1.0;
00113     ExcitationRatio = 1.0;
00114 
00115     theFastIntegralTable = NULL;
00116     theSlowIntegralTable = NULL;
00117     theReemissionIntegralTable = NULL;
00118 
00119     //verboseLevel = 2;
00120     //G4cout << " DsG4Scintillation set verboseLevel by hand to " << verboseLevel << G4endl;
00121 
00122     if (verboseLevel > 0) {
00123         G4cout << GetProcessName() << " is created " << G4endl;
00124     }
00125 
00126     BuildThePhysicsTable();
00127 
00128 }
00129 
00131 // Destructors
00133 
00134 DsG4Scintillation::~DsG4Scintillation() 
00135 {
00136     if (theFastIntegralTable != NULL) {
00137         theFastIntegralTable->clearAndDestroy();
00138         delete theFastIntegralTable;
00139     }
00140     if (theSlowIntegralTable != NULL) {
00141         theSlowIntegralTable->clearAndDestroy();
00142         delete theSlowIntegralTable;
00143     }
00144     if (theReemissionIntegralTable != NULL) {
00145         theReemissionIntegralTable->clearAndDestroy();
00146         delete theReemissionIntegralTable;
00147     }
00148 }
00149 
00151 // Methods
00153 
00154 // AtRestDoIt
00155 // ----------
00156 //
00157 G4VParticleChange*
00158 DsG4Scintillation::AtRestDoIt(const G4Track& aTrack, const G4Step& aStep)
00159 
00160 // This routine simply calls the equivalent PostStepDoIt since all the
00161 // necessary information resides in aStep.GetTotalEnergyDeposit()
00162 
00163 {
00164     return DsG4Scintillation::PostStepDoIt(aTrack, aStep);
00165 }
00166 
00167 // PostStepDoIt
00168 // -------------
00169 //
00170 G4VParticleChange*
00171 DsG4Scintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
00172 
00173 // This routine is called for each tracking step of a charged particle
00174 // in a scintillator. A Poisson/Gauss-distributed number of photons is 
00175 // generated according to the scintillation yield formula, distributed 
00176 // evenly along the track segment and uniformly into 4pi.
00177 
00178 {
00179     aParticleChange.Initialize(aTrack);
00180 
00181     if (m_noop) {               // do nothing, bail
00182         aParticleChange.SetNumberOfSecondaries(0);
00183         return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
00184     }
00185 
00186 
00187     G4String pname="";
00188     G4ThreeVector vertpos;
00189     G4double vertenergy=0.0;
00190     G4double reem_d=0.0;
00191     G4bool flagReemission= false;
00192     DsPhotonTrackInfo* reemittedTI=0;
00193     if (aTrack.GetDefinition() == G4OpticalPhoton::OpticalPhoton()) {
00194         G4Track *track=aStep.GetTrack();
00195         G4CompositeTrackInfo* composite=dynamic_cast<G4CompositeTrackInfo*>(track->GetUserInformation());
00196         reemittedTI = composite?dynamic_cast<DsPhotonTrackInfo*>( composite->GetPhotonTrackInfo() ):0;
00197         
00198         const G4VProcess* process = track->GetCreatorProcess();
00199         if(process) pname = process->GetProcessName();
00200 
00201         if (verboseLevel > 0) { 
00202           G4cout<<"Optical photon. Process name is " << pname<<G4endl;
00203         } 
00204         if(doBothProcess) {
00205             flagReemission= doReemission
00206                 && aTrack.GetTrackStatus() == fStopAndKill
00207                 && aStep.GetPostStepPoint()->GetStepStatus() != fGeomBoundary;     
00208         }
00209         else{
00210             flagReemission= doReemission
00211                 && aTrack.GetTrackStatus() == fStopAndKill
00212                 && aStep.GetPostStepPoint()->GetStepStatus() != fGeomBoundary
00213                 && pname=="Cerenkov";
00214         }
00215         if(verboseLevel > 0) {
00216             G4cout<<"flag of Reemission is "<<flagReemission<<"!!"<<G4endl;
00217         }
00218         if (!flagReemission) {
00219             return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
00220         }
00221     }
00222 
00223     G4double TotalEnergyDeposit = aStep.GetTotalEnergyDeposit();
00224     if (verboseLevel > 0 ) { 
00225       G4cout << " TotalEnergyDeposit " << TotalEnergyDeposit 
00226              << " material " << aTrack.GetMaterial()->GetName() << G4endl;
00227     }
00228     if (TotalEnergyDeposit <= 0.0 && !flagReemission) {
00229         return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
00230     }
00231 
00232     const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00233     const G4String aParticleName = aParticle->GetDefinition()->GetParticleName();
00234     const G4Material* aMaterial = aTrack.GetMaterial();
00235 
00236     G4MaterialPropertiesTable* aMaterialPropertiesTable =
00237         aMaterial->GetMaterialPropertiesTable();
00238     if (!aMaterialPropertiesTable)
00239         return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
00240 
00241     G4String FastTimeConstant = "FASTTIMECONSTANT";
00242     G4String SlowTimeConstant = "SLOWTIMECONSTANT";
00243     G4String strYieldRatio = "YIELDRATIO";
00244 
00245     
00246     if (aParticleName == "opticalphoton") {
00247       FastTimeConstant = "ReemissionFASTTIMECONSTANT";
00248       SlowTimeConstant = "ReemissionSLOWTIMECONSTANT";
00249       strYieldRatio = "ReemissionYIELDRATIO";
00250     }
00251     else if(aParticleName == "gamma" || aParticleName == "e+" || aParticleName == "e-") {
00252       FastTimeConstant = "GammaFASTTIMECONSTANT";
00253       SlowTimeConstant = "GammaSLOWTIMECONSTANT";
00254       strYieldRatio = "GammaYIELDRATIO";
00255       slowerTimeConstant = gammaSlowerTime;
00256       slowerRatio = gammaSlowerRatio;
00257     }
00258     else if(aParticleName == "alpha") {
00259       FastTimeConstant = "AlphaFASTTIMECONSTANT";
00260       SlowTimeConstant = "AlphaSLOWTIMECONSTANT";
00261       strYieldRatio = "AlphaYIELDRATIO";
00262       slowerTimeConstant = alphaSlowerTime;
00263       slowerRatio = alphaSlowerRatio;
00264     }
00265     else {
00266       FastTimeConstant = "NeutronFASTTIMECONSTANT";
00267       SlowTimeConstant = "NeutronSLOWTIMECONSTANT";
00268       strYieldRatio = "NeutronYIELDRATIO";
00269       slowerTimeConstant = neutronSlowerTime;
00270       slowerRatio = neutronSlowerRatio;
00271     }
00272 
00273     const G4MaterialPropertyVector* Fast_Intensity = 
00274         aMaterialPropertiesTable->GetProperty("FASTCOMPONENT"); 
00275     const G4MaterialPropertyVector* Slow_Intensity =
00276         aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
00277     const G4MaterialPropertyVector* Reemission_Prob =
00278         aMaterialPropertiesTable->GetProperty("REEMISSIONPROB");
00279     if (verboseLevel > 0 ) {
00280       G4cout << " MaterialPropertyVectors: Fast_Intensity " << Fast_Intensity 
00281              << " Slow_Intensity " << Slow_Intensity << " Reemission_Prob " << Reemission_Prob << G4endl;
00282     }
00283     if (!Fast_Intensity && !Slow_Intensity )
00284         return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
00285 
00286     G4int nscnt = 1;
00287     if (Fast_Intensity && Slow_Intensity) nscnt = 2;
00288     if ( verboseLevel > 0) {
00289       G4cout << " Fast_Intensity " << Fast_Intensity << " Slow_Intensity " << Slow_Intensity << " nscnt " << nscnt << G4endl;
00290     }
00291     G4StepPoint* pPreStepPoint  = aStep.GetPreStepPoint();
00292     G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
00293 
00294     G4ThreeVector x0 = pPreStepPoint->GetPosition();
00295     G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
00296     G4double      t0 = pPreStepPoint->GetGlobalTime();
00297 
00298     //Replace NumPhotons by NumTracks
00299     G4int NumTracks=0;
00300     G4double weight=1.0;
00301     if (flagReemission) {   
00302         if(verboseLevel > 0){   
00303             G4cout<<"the process name is "<<pname<<"!!"<<G4endl;}
00304         
00305         if ( Reemission_Prob == 0)
00306             return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
00307         G4double p_reemission=
00308             Reemission_Prob->GetProperty(aTrack.GetKineticEnergy());
00309         if (G4UniformRand() >= p_reemission)
00310             return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
00311         NumTracks= 1;
00312         weight= aTrack.GetWeight();
00313         if (verboseLevel > 0 ) {
00314             G4cout << " flagReemission " << flagReemission << " weight " << weight << G4endl;}
00315     }
00316     else {
00318         // J.B.Birks. The theory and practice of Scintillation Counting. 
00319         // Pergamon Press, 1964.      
00320         // For particles with energy much smaller than minimum ionization 
00321         // energy, the scintillation response is non-linear because of quenching  
00322         // effect. The light output is reduced by a parametric factor: 
00323         // 1/(1 + birk1*delta + birk2* delta^2). 
00324         // Delta is the energy loss per unit mass thickness. birk1 and birk2 
00325         // were measured for several organic scintillators.         
00326         // Here we use birk1 = 0.0125*g/cm2/MeV and ignore birk2.               
00327         // R.L.Craun and D.L.Smith. Nucl. Inst. and Meth., 80:239-244, 1970.   
00328         // Liang Zhan  01/27/2006 
00329         // /////////////////////////////////////////////////////////////////////
00330         
00331         
00332         G4double ScintillationYield = 0;
00333         {// Yield.  Material must have this or we lack raisins dayetras
00334             const G4MaterialPropertyVector* ptable =
00335                 aMaterialPropertiesTable->GetProperty("SCINTILLATIONYIELD");
00336             if (!ptable) {
00337                 G4cout << "ConstProperty: failed to get SCINTILLATIONYIELD"
00338                        << G4endl;
00339                 return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
00340             }
00341             ScintillationYield = ptable->GetProperty(0);
00342         }
00343 
00344         G4double ResolutionScale    = 1;
00345         {// Resolution Scale
00346             const G4MaterialPropertyVector* ptable =
00347                 aMaterialPropertiesTable->GetProperty("RESOLUTIONSCALE");
00348             if (ptable)
00349                 ResolutionScale = ptable->GetProperty(0);
00350         }
00351 
00352         G4double dE = TotalEnergyDeposit;
00353         G4double dx = aStep.GetStepLength();
00354         G4double dE_dx = dE/dx;
00355         if(aTrack.GetDefinition() == G4Gamma::Gamma() && dE > 0)
00356         { 
00357           G4LossTableManager* manager = G4LossTableManager::Instance();
00358           dE_dx = dE/manager->GetRange(G4Electron::Electron(), dE, aTrack.GetMaterialCutsCouple());
00359           //G4cout<<"gamma dE_dx = "<<dE_dx/(MeV/mm)<<"MeV/mm"<<G4endl;
00360         }
00361         
00362         G4double delta = dE_dx/aMaterial->GetDensity();//get scintillator density 
00363         //G4double birk1 = 0.0125*g/cm2/MeV;
00364         G4double birk1 = birksConstant1;
00365         if(abs(aParticle->GetCharge())>1.5)//for particle charge greater than 1.
00366             birk1 = 0.57*birk1;
00367         
00368         G4double birk2 = 0;
00369         //birk2 = (0.0031*g/MeV/cm2)*(0.0031*g/MeV/cm2);
00370         birk2 = birksConstant2;
00371         
00372         G4double QuenchedTotalEnergyDeposit 
00373             = TotalEnergyDeposit/(1+birk1*delta+birk2*delta*delta);
00374 
00375        //Add 300ns trick for muon simuation, by Haoqi Jan 27, 2011  
00376        if(FastMu300nsTrick)  {
00377            // cout<<"GlobalTime ="<<aStep.GetTrack()->GetGlobalTime()/ns<<endl;
00378            if(aStep.GetTrack()->GetGlobalTime()/ns>300) {
00379                ScintillationYield = YieldFactor * ScintillationYield;
00380            }
00381            else{
00382             ScintillationYield=0.;
00383            }
00384         }
00385         else {    
00386             ScintillationYield = YieldFactor * ScintillationYield; 
00387         }
00388 
00389         G4double MeanNumberOfPhotons= ScintillationYield * QuenchedTotalEnergyDeposit;
00390    
00391         // Implemented the fast simulation method from GLG4Scint
00392         // Jianglai 09-05-2006
00393         
00394         // randomize number of TRACKS (not photons)
00395         // this gets statistics right for number of PE after applying
00396         // boolean random choice to final absorbed track (change from
00397         // old method of applying binomial random choice to final absorbed
00398         // track, which did want poissonian number of photons divided
00399         // as evenly as possible into tracks)
00400         // Note for weight=1, there's no difference between tracks and photons.
00401         G4double MeanNumberOfTracks= MeanNumberOfPhotons/fPhotonWeight; 
00402         if ( fApplyPreQE ) {
00403             MeanNumberOfTracks*=fPreQE;
00404         }
00405         if (MeanNumberOfTracks > 10.) {
00406             G4double sigma = ResolutionScale * sqrt(MeanNumberOfTracks);
00407             NumTracks = G4int(G4RandGauss::shoot(MeanNumberOfTracks,sigma)+0.5);
00408         }
00409         else {
00410             NumTracks = G4int(G4Poisson(MeanNumberOfTracks));
00411         }
00412         if ( verboseLevel > 0 ) {
00413           G4cout << " Generated " << NumTracks << " scint photons. mean(scint photons) = " << MeanNumberOfTracks << G4endl;
00414         }
00415     }
00416     weight*=fPhotonWeight;
00417     if ( verboseLevel > 0 ) {
00418       G4cout << " set scint photon weight to " << weight << " after multiplying original weight by fPhotonWeight " << fPhotonWeight 
00419              << " NumTracks = " << NumTracks
00420              << G4endl;
00421     }
00422     // G4cerr<<"Scint weight is "<<weight<<G4endl;
00423     if (NumTracks <= 0) {
00424         // return unchanged particle and no secondaries 
00425         aParticleChange.SetNumberOfSecondaries(0);
00426         return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
00427     }
00428 
00430 
00431     aParticleChange.SetNumberOfSecondaries(NumTracks);
00432 
00433     if (fTrackSecondariesFirst) {
00434         if (!flagReemission) 
00435             if (aTrack.GetTrackStatus() == fAlive )
00436                 aParticleChange.ProposeTrackStatus(fSuspend);
00437     }
00438         
00440 
00441     G4int materialIndex = aMaterial->GetIndex();
00442 
00443     G4PhysicsOrderedFreeVector* ReemissionIntegral = NULL;
00444     ReemissionIntegral =
00445         (G4PhysicsOrderedFreeVector*)((*theReemissionIntegralTable)(materialIndex));
00446 
00447     // Retrieve the Scintillation Integral for this material  
00448     // new G4PhysicsOrderedFreeVector allocated to hold CII's
00449 
00450     G4int Num = NumTracks; //# tracks is now the loop control
00451         
00452     G4double fastTimeConstant = 0.0;
00453     { // Fast Time Constant
00454         const G4MaterialPropertyVector* ptable =
00455         aMaterialPropertiesTable->GetProperty(FastTimeConstant.c_str());
00456         if (verboseLevel > 0) {
00457           G4cout << " MaterialPropertyVector table " << ptable << " for FASTTIMECONSTANT"<<G4endl;
00458         }
00459         if (!ptable) ptable = aMaterialPropertiesTable->GetProperty("FASTTIMECONSTANT");
00460         if (ptable) {
00461             fastTimeConstant = ptable->GetProperty(0);
00462           if (verboseLevel > 0) { 
00463             G4cout << " dump fast time constant table " << G4endl;
00464             const_cast <G4MaterialPropertyVector*>(ptable)->DumpVector();
00465           }
00466         }
00467     }
00468 
00469     G4double slowTimeConstant = 0.0;
00470     { // Slow Time Constant
00471         const G4MaterialPropertyVector* ptable =
00472         aMaterialPropertiesTable->GetProperty(SlowTimeConstant.c_str());
00473         if (verboseLevel > 0) {
00474           G4cout << " MaterialPropertyVector table " << ptable << " for SLOWTIMECONSTANT"<<G4endl;
00475         }
00476         if(!ptable) ptable = aMaterialPropertiesTable->GetProperty("SLOWTIMECONSTANT");
00477         if (ptable){
00478           slowTimeConstant = ptable->GetProperty(0);
00479           if (verboseLevel > 0) { 
00480             G4cout << " dump slow time constant table " << G4endl;
00481             const_cast <G4MaterialPropertyVector*>(ptable)->DumpVector();
00482           }
00483         }
00484     }
00485 
00486     G4double YieldRatio = 0.0;
00487     { // Slow Time Constant
00488         const G4MaterialPropertyVector* ptable =
00489             aMaterialPropertiesTable->GetProperty(strYieldRatio.c_str());
00490         if(!ptable) ptable = aMaterialPropertiesTable->GetProperty("YIELDRATIO");
00491         if (ptable)
00492             YieldRatio = ptable->GetProperty(0);
00493         if (verboseLevel > 0) {
00494             G4cout << " YieldRatio = "<< YieldRatio << " and dump yield ratio table (yield ratio = fast/(fast+slow): " << G4endl;
00495             const_cast <G4MaterialPropertyVector*>(ptable)->DumpVector();
00496           }
00497     }
00498 
00499 
00500     //loop over fast/slow scintillations
00501     for (G4int scnt = 1; scnt <= nscnt; scnt++) {
00502 
00503         G4double ScintillationTime = 0.*ns;
00504         G4PhysicsOrderedFreeVector* ScintillationIntegral = NULL;
00505          
00506         if (scnt == 1) {//fast
00507             if (nscnt == 1) {
00508                 if(Fast_Intensity){
00509                     ScintillationTime   = fastTimeConstant;
00510                     ScintillationIntegral =
00511                         (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
00512                 }
00513                 if(Slow_Intensity){
00514                     ScintillationTime   = slowTimeConstant;
00515                     ScintillationIntegral =
00516                         (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
00517                 }
00518             }
00519             else {
00520                 if ( ExcitationRatio == 1.0 ) {
00521                   Num = G4int( 0.5 +  (min(YieldRatio,1.0) * NumTracks) );  // round off, not truncation
00522                 }
00523                 else {
00524                   Num = G4int( 0.5 +  (min(ExcitationRatio,1.0) * NumTracks));
00525                 }
00526                 if ( verboseLevel>1 ){
00527                   G4cout << "Generate Num " << Num << " optical photons with fast component using NumTracks " 
00528                          << NumTracks << " YieldRatio " << YieldRatio << " ExcitationRatio " << ExcitationRatio 
00529                          << " min(YieldRatio,1.)*NumTracks = " <<  min(YieldRatio,1.)*NumTracks 
00530                          << " min(ExcitationRatio,1.)*NumTracks = " <<  min(ExcitationRatio,1.)*NumTracks 
00531                          << G4endl;
00532                 }
00533                 ScintillationTime   = fastTimeConstant;
00534                 ScintillationIntegral =
00535                     (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
00536             }
00537         }
00538         else {//slow
00539             Num = NumTracks - Num;
00540             ScintillationTime   =   slowTimeConstant;
00541             ScintillationIntegral =
00542                 (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
00543         }
00544         if (verboseLevel > 0) {
00545           G4cout << "generate " << Num << " optical photons with scintTime " << ScintillationTime 
00546                  << " slowTimeConstant " << slowTimeConstant << " fastTimeConstant " << fastTimeConstant << G4endl;
00547         }
00548 
00549         if (!ScintillationIntegral) continue;
00550         
00551         // Max Scintillation Integral
00552                 
00553         for (G4int i = 0; i < Num; i++) { //Num is # of 2ndary tracks now
00554             // Determine photon energy
00555 
00556         if(scnt == 2) {
00557             ScintillationTime   =   slowTimeConstant;
00558             if((G4UniformRand() < slowerRatio) && (!flagReemission)) { 
00559               ScintillationTime = slowerTimeConstant;
00560             }
00561         }
00562 
00563             G4double sampledEnergy;
00564             if ( !flagReemission ) {
00565                 // normal scintillation
00566                 G4double CIIvalue = G4UniformRand()*
00567                     ScintillationIntegral->GetMaxValue();
00568                 sampledEnergy=
00569                     ScintillationIntegral->GetEnergy(CIIvalue);
00570 
00571                 if (verboseLevel>1) 
00572                     {
00573                         G4cout << "sampledEnergy = " << sampledEnergy << G4endl;
00574                         G4cout << "CIIvalue =        " << CIIvalue << G4endl;
00575                     }
00576             }
00577             else {
00578                 // reemission, the sample method need modification
00579                 G4double CIIvalue = G4UniformRand()*
00580                     ScintillationIntegral->GetMaxValue();
00581                 if (CIIvalue == 0.0) {
00582                     // return unchanged particle and no secondaries  
00583                     aParticleChange.SetNumberOfSecondaries(0);
00584                     return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
00585                 }
00586                 sampledEnergy=
00587                     ScintillationIntegral->GetEnergy(CIIvalue);
00588                 if (verboseLevel>1) {
00589                     G4cout << "oldEnergy = " <<aTrack.GetKineticEnergy() << G4endl;
00590                     G4cout << "reemittedSampledEnergy = " << sampledEnergy
00591                            << "\nreemittedCIIvalue =        " << CIIvalue << G4endl;
00592                 }
00593             }
00594 
00595             // Generate random photon direction
00596 
00597             G4double cost = 1. - 2.*G4UniformRand();
00598             G4double sint = sqrt((1.-cost)*(1.+cost));
00599 
00600             G4double phi = twopi*G4UniformRand();
00601             G4double sinp = sin(phi);
00602             G4double cosp = cos(phi);
00603 
00604             G4double px = sint*cosp;
00605             G4double py = sint*sinp;
00606             G4double pz = cost;
00607 
00608             // Create photon momentum direction vector 
00609 
00610             G4ParticleMomentum photonMomentum(px, py, pz);
00611 
00612             // Determine polarization of new photon 
00613 
00614             G4double sx = cost*cosp;
00615             G4double sy = cost*sinp; 
00616             G4double sz = -sint;
00617 
00618             G4ThreeVector photonPolarization(sx, sy, sz);
00619 
00620             G4ThreeVector perp = photonMomentum.cross(photonPolarization);
00621 
00622             phi = twopi*G4UniformRand();
00623             sinp = sin(phi);
00624             cosp = cos(phi);
00625 
00626             photonPolarization = cosp * photonPolarization + sinp * perp;
00627 
00628             photonPolarization = photonPolarization.unit();
00629 
00630             // Generate a new photon:
00631 
00632             G4DynamicParticle* aScintillationPhoton =
00633                 new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(), 
00634                                       photonMomentum);
00635             aScintillationPhoton->SetPolarization
00636                 (photonPolarization.x(),
00637                  photonPolarization.y(),
00638                  photonPolarization.z());
00639 
00640             aScintillationPhoton->SetKineticEnergy(sampledEnergy);
00641 
00642             // Generate new G4Track object:
00643 
00644             G4double rand=0;
00645             G4ThreeVector aSecondaryPosition;
00646             G4double deltaTime;
00647             if (flagReemission) {
00648                 deltaTime= pPostStepPoint->GetGlobalTime() - t0
00649                            -ScintillationTime * log( G4UniformRand() );
00650                 aSecondaryPosition= pPostStepPoint->GetPosition();
00651                 vertpos = aTrack.GetVertexPosition();
00652                 vertenergy = aTrack.GetKineticEnergy();
00653                 reem_d = 
00654                     sqrt( pow( aSecondaryPosition.x()-vertpos.x(), 2)
00655                           + pow( aSecondaryPosition.y()-vertpos.y(), 2)
00656                           + pow( aSecondaryPosition.z()-vertpos.z(), 2) );
00657             }
00658             else {
00659                 if (aParticle->GetDefinition()->GetPDGCharge() != 0) 
00660                     {
00661                         rand = G4UniformRand();
00662                     }
00663                 else
00664                     {
00665                         rand = 1.0;
00666                     }
00667 
00668                 G4double delta = rand * aStep.GetStepLength();
00669                 deltaTime = delta /
00670                     ((pPreStepPoint->GetVelocity()+
00671                       pPostStepPoint->GetVelocity())/2.);
00672 
00673                 deltaTime = deltaTime - 
00674                     ScintillationTime * log( G4UniformRand() );
00675 
00676                 aSecondaryPosition =
00677                     x0 + rand * aStep.GetDeltaPosition();
00678             }
00679             G4double aSecondaryTime = t0 + deltaTime;
00680 
00681             if ( verboseLevel>1 ){
00682               G4cout << "Generate " << i << "th scintillation photon at relative time(ns) " << deltaTime 
00683                      << " with ScintillationTime " << ScintillationTime << " flagReemission " << flagReemission << G4endl;
00684             }
00685             G4Track* aSecondaryTrack = 
00686                 new G4Track(aScintillationPhoton,aSecondaryTime,aSecondaryPosition);
00687 
00688             G4CompositeTrackInfo* comp=new G4CompositeTrackInfo();
00689             DsPhotonTrackInfo* trackinf=new DsPhotonTrackInfo();
00690             if ( flagReemission ){
00691                 if ( reemittedTI ) *trackinf = *reemittedTI;
00692                 trackinf->SetReemitted();
00693             }
00694             else if ( fApplyPreQE ) {
00695                 trackinf->SetMode(DsPhotonTrackInfo::kQEPreScale);
00696                 trackinf->SetQE(fPreQE);
00697             }
00698             comp->SetPhotonTrackInfo(trackinf);
00699             aSecondaryTrack->SetUserInformation(comp);
00700                 
00701             aSecondaryTrack->SetWeight( weight );
00702             aSecondaryTrack->SetTouchableHandle(aStep.GetPreStepPoint()->GetTouchableHandle());
00703             // aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0);//this is wrong
00704                 
00705             aSecondaryTrack->SetParentID(aTrack.GetTrackID());
00706                 
00707             // add the secondary to the ParticleChange object
00708             aParticleChange.SetSecondaryWeightByProcess( true ); // recommended
00709             aParticleChange.AddSecondary(aSecondaryTrack);
00710                 
00711             aSecondaryTrack->SetWeight( weight );
00712             if ( verboseLevel > 0 ) {
00713               G4cout << " aSecondaryTrack->SetWeight( " << weight<< " ) ; aSecondaryTrack->GetWeight() = " << aSecondaryTrack->GetWeight() << G4endl;}
00714             // The above line is necessary because AddSecondary() 
00715             // overrides our setting of the secondary track weight, 
00716             // in Geant4.3.1 & earlier. (and also later, at least 
00717             // until Geant4.7 (and beyond?)
00718             //  -- maybe not required if SetWeightByProcess(true) called,
00719             //  but we do both, just to be sure)
00720         }
00721     } // end loop over fast/slow scints
00722 
00723     if (verboseLevel > 0) {
00724         G4cout << "\n Exiting from G4Scintillation::DoIt -- NumberOfSecondaries = " 
00725                << aParticleChange.GetNumberOfSecondaries() << G4endl;
00726     }
00727 
00728     return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
00729 }
00730 
00731 // BuildThePhysicsTable for the scintillation process
00732 // --------------------------------------------------
00733 //
00734 
00735 void DsG4Scintillation::BuildThePhysicsTable()
00736 {
00737     if (theFastIntegralTable && theSlowIntegralTable && theReemissionIntegralTable) return;
00738 
00739     const G4MaterialTable* theMaterialTable = 
00740         G4Material::GetMaterialTable();
00741     G4int numOfMaterials = G4Material::GetNumberOfMaterials();
00742 
00743     // create new physics table
00744     if (verboseLevel > 0) {
00745       G4cout << " theFastIntegralTable " << theFastIntegralTable 
00746              << " theSlowIntegralTable " << theSlowIntegralTable 
00747              << " theReemissionIntegralTable " << theReemissionIntegralTable << G4endl;
00748     }
00749     if(!theFastIntegralTable)theFastIntegralTable = new G4PhysicsTable(numOfMaterials);
00750     if(!theSlowIntegralTable)theSlowIntegralTable = new G4PhysicsTable(numOfMaterials);
00751     if(!theReemissionIntegralTable)theReemissionIntegralTable
00752                                        = new G4PhysicsTable(numOfMaterials);
00753     if (verboseLevel > 0) {
00754       G4cout << " building the physics tables for the scintillation process " << G4endl;
00755     }
00756     // loop for materials
00757 
00758     for (G4int i=0 ; i < numOfMaterials; i++) {
00759         G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
00760             new G4PhysicsOrderedFreeVector();
00761         G4PhysicsOrderedFreeVector* bPhysicsOrderedFreeVector =
00762             new G4PhysicsOrderedFreeVector();
00763         G4PhysicsOrderedFreeVector* cPhysicsOrderedFreeVector =
00764             new G4PhysicsOrderedFreeVector();
00765 
00766         // Retrieve vector of scintillation wavelength intensity for
00767         // the material from the material's optical properties table.
00768 
00769         G4Material* aMaterial = (*theMaterialTable)[i];
00770 
00771         G4MaterialPropertiesTable* aMaterialPropertiesTable =
00772             aMaterial->GetMaterialPropertiesTable();
00773 
00774         if (aMaterialPropertiesTable) {
00775 
00776             G4MaterialPropertyVector* theFastLightVector = 
00777                 aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
00778 
00779             if (theFastLightVector) {
00780               if (verboseLevel > 0) {
00781                 G4cout << " Building the material properties table for FASTCOMPONENT" << G4endl;
00782               }
00783                 // Retrieve the first intensity point in vector
00784                 // of (photon energy, intensity) pairs 
00785 
00786                 theFastLightVector->ResetIterator();
00787                 ++(*theFastLightVector);        // advance to 1st entry 
00788 
00789                 G4double currentIN = theFastLightVector->
00790                     GetProperty();
00791 
00792                 if (currentIN >= 0.0) {
00793 
00794                     // Create first (photon energy, Scintillation 
00795                     // Integral pair  
00796 
00797                     G4double currentPM = theFastLightVector->
00798                         GetPhotonEnergy();
00799 
00800                     G4double currentCII = 0.0;
00801 
00802                     aPhysicsOrderedFreeVector->
00803                         InsertValues(currentPM , currentCII);
00804 
00805                     // Set previous values to current ones prior to loop
00806 
00807                     G4double prevPM  = currentPM;
00808                     G4double prevCII = currentCII;
00809                     G4double prevIN  = currentIN;
00810 
00811                     // loop over all (photon energy, intensity)
00812                     // pairs stored for this material  
00813 
00814                     while(++(*theFastLightVector)) {
00815                         currentPM = theFastLightVector->
00816                             GetPhotonEnergy();
00817 
00818                         currentIN=theFastLightVector->  
00819                             GetProperty();
00820 
00821                         currentCII = 0.5 * (prevIN + currentIN);
00822 
00823                         currentCII = prevCII +
00824                             (currentPM - prevPM) * currentCII;
00825 
00826                         aPhysicsOrderedFreeVector->
00827                             InsertValues(currentPM, currentCII);
00828 
00829                         prevPM  = currentPM;
00830                         prevCII = currentCII;
00831                         prevIN  = currentIN;
00832                     }
00833 
00834                 }
00835             }
00836 
00837             G4MaterialPropertyVector* theSlowLightVector =
00838                 aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
00839 
00840             if (theSlowLightVector) {
00841               if (verboseLevel > 0) {
00842                 G4cout << " Building the material properties table for SLOWCOMPONENT" << G4endl;
00843               }
00844 
00845                 // Retrieve the first intensity point in vector
00846                 // of (photon energy, intensity) pairs
00847 
00848                 theSlowLightVector->ResetIterator();
00849                 ++(*theSlowLightVector);  // advance to 1st entry
00850 
00851                 G4double currentIN = theSlowLightVector->
00852                     GetProperty();
00853 
00854                 if (currentIN >= 0.0) {
00855 
00856                     // Create first (photon energy, Scintillation
00857                     // Integral pair
00858 
00859                     G4double currentPM = theSlowLightVector->
00860                         GetPhotonEnergy();
00861 
00862                     G4double currentCII = 0.0;
00863 
00864                     bPhysicsOrderedFreeVector->
00865                         InsertValues(currentPM , currentCII);
00866 
00867                     // Set previous values to current ones prior to loop
00868 
00869                     G4double prevPM  = currentPM;
00870                     G4double prevCII = currentCII;
00871                     G4double prevIN  = currentIN;
00872 
00873                     // loop over all (photon energy, intensity)
00874                     // pairs stored for this material
00875 
00876                     while(++(*theSlowLightVector)) {
00877                         currentPM = theSlowLightVector->
00878                             GetPhotonEnergy();
00879 
00880                         currentIN=theSlowLightVector->
00881                             GetProperty();
00882 
00883                         currentCII = 0.5 * (prevIN + currentIN);
00884 
00885                         currentCII = prevCII +
00886                             (currentPM - prevPM) * currentCII;
00887 
00888                         bPhysicsOrderedFreeVector->
00889                             InsertValues(currentPM, currentCII);
00890 
00891                         prevPM  = currentPM;
00892                         prevCII = currentCII;
00893                         prevIN  = currentIN;
00894                     }
00895 
00896                 }
00897             }
00898 
00899             G4MaterialPropertyVector* theReemissionVector =
00900                 aMaterialPropertiesTable->GetProperty("REEMISSIONPROB");
00901 
00902             if (theReemissionVector) {
00903               if (verboseLevel > 0) {
00904                 G4cout << " Building the material properties table for REEMISSIONPROB" << G4endl;
00905               }
00906 
00907                 // Retrieve the first intensity point in vector
00908                 // of (photon energy, intensity) pairs
00909 
00910                 theReemissionVector->ResetIterator();
00911                 ++(*theReemissionVector);  // advance to 1st entry
00912 
00913                 G4double currentIN = theReemissionVector->
00914                     GetProperty();
00915 
00916                 if (currentIN >= 0.0) {
00917 
00918                     // Create first (photon energy, Scintillation
00919                     // Integral pair
00920 
00921                     G4double currentPM = theReemissionVector->
00922                         GetPhotonEnergy();
00923 
00924                     G4double currentCII = 0.0;
00925 
00926                     cPhysicsOrderedFreeVector->
00927                         InsertValues(currentPM , currentCII);
00928 
00929                     // Set previous values to current ones prior to loop
00930 
00931                     G4double prevPM  = currentPM;
00932                     G4double prevCII = currentCII;
00933                     G4double prevIN  = currentIN;
00934 
00935                     // loop over all (photon energy, intensity)
00936                     // pairs stored for this material
00937 
00938                     while(++(*theReemissionVector)) {
00939                         currentPM = theReemissionVector->
00940                             GetPhotonEnergy();
00941 
00942                         currentIN=theReemissionVector->
00943                             GetProperty();
00944 
00945                         currentCII = 0.5 * (prevIN + currentIN);
00946 
00947                         currentCII = prevCII +
00948                             (currentPM - prevPM) * currentCII;
00949 
00950                         cPhysicsOrderedFreeVector->
00951                             InsertValues(currentPM, currentCII);
00952 
00953                         prevPM  = currentPM;
00954                         prevCII = currentCII;
00955                         prevIN  = currentIN;
00956                     }
00957 
00958                 }
00959             }
00960 
00961         }
00962 
00963         // The scintillation integral(s) for a given material
00964         // will be inserted in the table(s) according to the
00965         // position of the material in the material table.
00966 
00967         theFastIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
00968         theSlowIntegralTable->insertAt(i,bPhysicsOrderedFreeVector);
00969         theReemissionIntegralTable->insertAt(i,cPhysicsOrderedFreeVector);
00970     }
00971 }
00972 
00973 // GetMeanFreePath
00974 // ---------------
00975 //
00976 
00977 G4double DsG4Scintillation::GetMeanFreePath(const G4Track&,
00978                                             G4double ,
00979                                             G4ForceCondition* condition)
00980 {
00981     *condition = StronglyForced;
00982 
00983     return DBL_MAX;
00984 
00985 }
00986 
00987 // GetMeanLifeTime
00988 // ---------------
00989 //
00990 
00991 G4double DsG4Scintillation::GetMeanLifeTime(const G4Track&,
00992                                             G4ForceCondition* condition)
00993 {
00994     *condition = Forced;
00995 
00996     return DBL_MAX;
00997 
00998 }
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:17:58 for DetSim by doxygen 1.7.4