Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/optical/src/G4OpWLS.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/optical/src/G4OpWLS.cc (Version 11.3.0) and /processes/optical/src/G4OpWLS.cc (Version 10.0.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
                                                   >>  27 // $Id: G4OpWLS.cc 79221 2014-02-20 14:58:02Z gcosmo $
 27 //                                                 28 //
 28 //////////////////////////////////////////////     29 ////////////////////////////////////////////////////////////////////////
 29 // Optical Photon WaveLength Shifting (WLS) Cl     30 // Optical Photon WaveLength Shifting (WLS) Class Implementation
 30 //////////////////////////////////////////////     31 ////////////////////////////////////////////////////////////////////////
 31 //                                                 32 //
 32 // File:        G4OpWLS.cc                         33 // File:        G4OpWLS.cc
 33 // Description: Discrete Process -- Wavelength     34 // Description: Discrete Process -- Wavelength Shifting of Optical Photons
 34 // Version:     1.0                                35 // Version:     1.0
 35 // Created:     2003-05-13                         36 // Created:     2003-05-13
 36 // Author:      John Paul Archambault              37 // Author:      John Paul Archambault
 37 //              (Adaptation of G4Scintillation     38 //              (Adaptation of G4Scintillation and G4OpAbsorption)
 38 // Updated:     2005-07-28 - add G4ProcessType     39 // Updated:     2005-07-28 - add G4ProcessType to constructor
 39 //              2006-05-07 - add G4VWLSTimeGen     40 //              2006-05-07 - add G4VWLSTimeGeneratorProfile
                                                   >>  41 // mail:        gum@triumf.ca
                                                   >>  42 //              jparcham@phys.ualberta.ca
 40 //                                                 43 //
 41 //////////////////////////////////////////////     44 ////////////////////////////////////////////////////////////////////////
 42                                                    45 
 43 #include "G4OpWLS.hh"                              46 #include "G4OpWLS.hh"
                                                   >>  47 
 44 #include "G4ios.hh"                                48 #include "G4ios.hh"
 45 #include "G4PhysicalConstants.hh"                  49 #include "G4PhysicalConstants.hh"
 46 #include "G4SystemOfUnits.hh"                      50 #include "G4SystemOfUnits.hh"
 47 #include "G4OpProcessSubType.hh"                   51 #include "G4OpProcessSubType.hh"
 48 #include "G4Poisson.hh"                        <<  52 
 49 #include "G4OpticalParameters.hh"              << 
 50 #include "G4WLSTimeGeneratorProfileDelta.hh"       53 #include "G4WLSTimeGeneratorProfileDelta.hh"
 51 #include "G4WLSTimeGeneratorProfileExponential     54 #include "G4WLSTimeGeneratorProfileExponential.hh"
 52                                                    55 
 53 //....oooOO0OOooo........oooOO0OOooo........oo <<  56 /////////////////////////
                                                   >>  57 // Class Implementation
                                                   >>  58 /////////////////////////
                                                   >>  59 
                                                   >>  60         //////////////////////
                                                   >>  61         // static data members
                                                   >>  62         //////////////////////
                                                   >>  63 
                                                   >>  64 G4VWLSTimeGeneratorProfile* G4OpWLS::WLSTimeGeneratorProfile = 0;
                                                   >>  65 
                                                   >>  66 /////////////////
                                                   >>  67 // Constructors
                                                   >>  68 /////////////////
                                                   >>  69 
 54 G4OpWLS::G4OpWLS(const G4String& processName,      70 G4OpWLS::G4OpWLS(const G4String& processName, G4ProcessType type)
 55   : G4VDiscreteProcess(processName, type)          71   : G4VDiscreteProcess(processName, type)
 56 {                                                  72 {
 57   WLSTimeGeneratorProfile = nullptr;           << 
 58   Initialise();                                << 
 59   SetProcessSubType(fOpWLS);                       73   SetProcessSubType(fOpWLS);
 60   theIntegralTable = nullptr;                  << 
 61                                                    74 
 62   if(verboseLevel > 0)                         <<  75   theIntegralTable = NULL;
                                                   >>  76   if(!WLSTimeGeneratorProfile)
                                                   >>  77   { WLSTimeGeneratorProfile = new G4WLSTimeGeneratorProfileDelta("WLSTimeGeneratorProfileDelta"); }
                                                   >>  78  
                                                   >>  79   if (verboseLevel>0) {
 63     G4cout << GetProcessName() << " is created     80     G4cout << GetProcessName() << " is created " << G4endl;
                                                   >>  81   }
 64 }                                                  82 }
 65                                                    83 
 66 //....oooOO0OOooo........oooOO0OOooo........oo <<  84 ////////////////
                                                   >>  85 // Destructors
                                                   >>  86 ////////////////
                                                   >>  87 
 67 G4OpWLS::~G4OpWLS()                                88 G4OpWLS::~G4OpWLS()
 68 {                                                  89 {
 69   if(theIntegralTable)                         <<  90   if (theIntegralTable != 0) {
 70   {                                            << 
 71     theIntegralTable->clearAndDestroy();           91     theIntegralTable->clearAndDestroy();
 72     delete theIntegralTable;                       92     delete theIntegralTable;
 73   }                                                93   }
 74   delete WLSTimeGeneratorProfile;              << 
 75 }                                                  94 }
 76                                                    95 
 77 //....oooOO0OOooo........oooOO0OOooo........oo <<  96 ////////////
 78 void G4OpWLS::PreparePhysicsTable(const G4Part <<  97 // Methods
                                                   >>  98 ////////////
 79                                                    99 
 80 //....oooOO0OOooo........oooOO0OOooo........oo << 100 void G4OpWLS::BuildPhysicsTable(const G4ParticleDefinition&)
 81 void G4OpWLS::Initialise()                     << 
 82 {                                                 101 {
 83   G4OpticalParameters* params = G4OpticalParam << 102     if (!theIntegralTable) BuildThePhysicsTable();
 84   SetVerboseLevel(params->GetWLSVerboseLevel() << 
 85   UseTimeProfile(params->GetWLSTimeProfile()); << 
 86 }                                                 103 }
 87                                                   104 
 88 //....oooOO0OOooo........oooOO0OOooo........oo << 105 // PostStepDoIt
 89 G4VParticleChange* G4OpWLS::PostStepDoIt(const << 106 // -------------
 90                                          const << 107 //
                                                   >> 108 G4VParticleChange*
                                                   >> 109 G4OpWLS::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
 91 {                                                 110 {
 92   std::vector<G4Track*> proposedSecondaries;   << 
 93   aParticleChange.Initialize(aTrack);             111   aParticleChange.Initialize(aTrack);
                                                   >> 112   
 94   aParticleChange.ProposeTrackStatus(fStopAndK    113   aParticleChange.ProposeTrackStatus(fStopAndKill);
 95                                                   114 
 96   if(verboseLevel > 1)                         << 115   if (verboseLevel>0) {
 97   {                                            << 116     G4cout << "\n** Photon absorbed! **" << G4endl;
 98     G4cout << "\n** G4OpWLS: Photon absorbed!  << 
 99   }                                               117   }
                                                   >> 118   
                                                   >> 119   const G4Material* aMaterial = aTrack.GetMaterial();
100                                                   120 
101   G4StepPoint* pPostStepPoint = aStep.GetPostS    121   G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
102   G4MaterialPropertiesTable* MPT =             << 122     
103     aTrack.GetMaterial()->GetMaterialPropertie << 123   G4MaterialPropertiesTable* aMaterialPropertiesTable =
104   if(!MPT)                                     << 124     aMaterial->GetMaterialPropertiesTable();
105   {                                            << 125   if (!aMaterialPropertiesTable)
106     return G4VDiscreteProcess::PostStepDoIt(aT    126     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
107   }                                            << 127 
108   if(!MPT->GetProperty(kWLSCOMPONENT))         << 128   const G4MaterialPropertyVector* WLS_Intensity = 
109   {                                            << 129     aMaterialPropertiesTable->GetProperty("WLSCOMPONENT"); 
                                                   >> 130 
                                                   >> 131   if (!WLS_Intensity)
110     return G4VDiscreteProcess::PostStepDoIt(aT    132     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
111   }                                            << 
112                                                   133 
113   G4int NumPhotons = 1;                           134   G4int NumPhotons = 1;
114   if(MPT->ConstPropertyExists(kWLSMEANNUMBERPH << 135 
115   {                                            << 136   if (aMaterialPropertiesTable->ConstPropertyExists("WLSMEANNUMBERPHOTONS")) {
116     G4double MeanNumberOfPhotons = MPT->GetCon << 137 
117     NumPhotons                   = G4int(G4Poi << 138      G4double MeanNumberOfPhotons = aMaterialPropertiesTable->
118     if(NumPhotons <= 0)                        << 139                                     GetConstProperty("WLSMEANNUMBERPHOTONS");
119     {                                          << 140 
120       // return unchanged particle and no seco << 141      NumPhotons = G4int(G4Poisson(MeanNumberOfPhotons));
121       aParticleChange.SetNumberOfSecondaries(0 << 142 
122       return G4VDiscreteProcess::PostStepDoIt( << 143      if (NumPhotons <= 0) {
123     }                                          << 144 
                                                   >> 145         // return unchanged particle and no secondaries
                                                   >> 146 
                                                   >> 147         aParticleChange.SetNumberOfSecondaries(0);
                                                   >> 148 
                                                   >> 149         return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
                                                   >> 150 
                                                   >> 151      }
                                                   >> 152 
124   }                                               153   }
125                                                   154 
126   // Retrieve the WLS Integral for this materi << 155   aParticleChange.SetNumberOfSecondaries(NumPhotons);
127   // new G4PhysicsFreeVector allocated to hold << 156 
128   G4double primaryEnergy = aTrack.GetDynamicPa    157   G4double primaryEnergy = aTrack.GetDynamicParticle()->GetKineticEnergy();
129   G4double WLSTime       = 0.;                 << 
130   G4PhysicsFreeVector* WLSIntegral = nullptr;  << 
131                                                   158 
132   WLSTime     = MPT->GetConstProperty(kWLSTIME << 159   G4int materialIndex = aMaterial->GetIndex();
133   WLSIntegral = (G4PhysicsFreeVector*) ((*theI << 160 
134     aTrack.GetMaterial()->GetIndex()));        << 161   // Retrieve the WLS Integral for this material
                                                   >> 162   // new G4PhysicsOrderedFreeVector allocated to hold CII's
135                                                   163 
                                                   >> 164   G4double WLSTime = 0.*ns;
                                                   >> 165   G4PhysicsOrderedFreeVector* WLSIntegral = 0;
                                                   >> 166 
                                                   >> 167   WLSTime   = aMaterialPropertiesTable->
                                                   >> 168     GetConstProperty("WLSTIMECONSTANT");
                                                   >> 169   WLSIntegral =
                                                   >> 170     (G4PhysicsOrderedFreeVector*)((*theIntegralTable)(materialIndex));
                                                   >> 171    
136   // Max WLS Integral                             172   // Max WLS Integral
137   G4double CIImax       = WLSIntegral->GetMaxV << 173   
                                                   >> 174   G4double CIImax = WLSIntegral->GetMaxValue();
                                                   >> 175  
138   G4int NumberOfPhotons = NumPhotons;             176   G4int NumberOfPhotons = NumPhotons;
                                                   >> 177  
                                                   >> 178   for (G4int i = 0; i < NumPhotons; i++) {
139                                                   179 
140   for(G4int i = 0; i < NumPhotons; ++i)        << 
141   {                                            << 
142     G4double sampledEnergy;                       180     G4double sampledEnergy;
                                                   >> 181     
143     // Make sure the energy of the secondary i    182     // Make sure the energy of the secondary is less than that of the primary
144     for(G4int j = 1; j <= 100; ++j)            << 183 
145     {                                          << 184     for (G4int j = 1; j <= 100; j++) {
146       // Determine photon energy               << 185 
147       G4double CIIvalue = G4UniformRand() * CI << 186         // Determine photon energy
148       sampledEnergy     = WLSIntegral->GetEner << 187 
149       if(sampledEnergy <= primaryEnergy)       << 188         G4double CIIvalue = G4UniformRand()*CIImax;
150         break;                                 << 189         sampledEnergy = WLSIntegral->GetEnergy(CIIvalue);
151     }                                          << 190 
152     // If no such energy can be sampled, retur << 191         if (verboseLevel>1) {
153     if(sampledEnergy > primaryEnergy)          << 192            G4cout << "sampledEnergy = " << sampledEnergy << G4endl;
154     {                                          << 193            G4cout << "CIIvalue =      " << CIIvalue << G4endl;
155       if(verboseLevel > 1)                     << 
156       {                                        << 
157         G4cout << " *** G4OpWLS: One less WLS  << 
158                << G4endl;                      << 
159       }                                        << 
160       NumberOfPhotons--;                       << 
161       if(NumberOfPhotons == 0)                 << 
162       {                                        << 
163         if(verboseLevel > 1)                   << 
164         {                                      << 
165           G4cout                               << 
166             << " *** G4OpWLS: No WLS photon ca << 
167             << G4endl;                         << 
168         }                                         194         }
169         // return unchanged particle and no se << 195 
170         aParticleChange.SetNumberOfSecondaries << 196         if (sampledEnergy <= primaryEnergy) break;
171         return G4VDiscreteProcess::PostStepDoI << 
172       }                                        << 
173       continue;                                << 
174     }                                             197     }
175     else if(verboseLevel > 1)                  << 198 
176     {                                          << 199     // If no such energy can be sampled, return one less secondary, or none
177       G4cout << "G4OpWLS: Created photon with  << 200 
178              << G4endl;                        << 201     if (sampledEnergy > primaryEnergy) {
                                                   >> 202        if (verboseLevel>1)
                                                   >> 203        G4cout << " *** One less WLS photon will be returned ***" << G4endl;
                                                   >> 204        NumberOfPhotons--;
                                                   >> 205        aParticleChange.SetNumberOfSecondaries(NumberOfPhotons);
                                                   >> 206        if (NumberOfPhotons == 0) {
                                                   >> 207           if (verboseLevel>1)
                                                   >> 208           G4cout << " *** No WLS photon can be sampled for this primary ***"
                                                   >> 209                  << G4endl;
                                                   >> 210           // return unchanged particle and no secondaries
                                                   >> 211           return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
                                                   >> 212        }
                                                   >> 213        continue;
179     }                                             214     }
180                                                   215 
181     // Generate random photon direction           216     // Generate random photon direction
182     G4double cost = 1. - 2. * G4UniformRand(); << 217     
183     G4double sint = std::sqrt((1. - cost) * (1 << 218     G4double cost = 1. - 2.*G4UniformRand();
184     G4double phi  = twopi * G4UniformRand();   << 219     G4double sint = std::sqrt((1.-cost)*(1.+cost));
                                                   >> 220 
                                                   >> 221     G4double phi = twopi*G4UniformRand();
185     G4double sinp = std::sin(phi);                222     G4double sinp = std::sin(phi);
186     G4double cosp = std::cos(phi);                223     G4double cosp = std::cos(phi);
187     G4ParticleMomentum photonMomentum(sint * c << 224     
188                                                << 225     G4double px = sint*cosp;
189     G4ThreeVector photonPolarization(cost * co << 226     G4double py = sint*sinp;
                                                   >> 227     G4double pz = cost;
                                                   >> 228     
                                                   >> 229     // Create photon momentum direction vector
                                                   >> 230     
                                                   >> 231     G4ParticleMomentum photonMomentum(px, py, pz);
                                                   >> 232     
                                                   >> 233     // Determine polarization of new photon
                                                   >> 234     
                                                   >> 235     G4double sx = cost*cosp;
                                                   >> 236     G4double sy = cost*sinp;
                                                   >> 237     G4double sz = -sint;
                                                   >> 238     
                                                   >> 239     G4ThreeVector photonPolarization(sx, sy, sz);
                                                   >> 240     
190     G4ThreeVector perp = photonMomentum.cross(    241     G4ThreeVector perp = photonMomentum.cross(photonPolarization);
191                                                << 242     
192     phi                = twopi * G4UniformRand << 243     phi = twopi*G4UniformRand();
193     sinp               = std::sin(phi);        << 244     sinp = std::sin(phi);
194     cosp               = std::cos(phi);        << 245     cosp = std::cos(phi);
195     photonPolarization = (cosp * photonPolariz << 246     
196                                                << 247     photonPolarization = cosp * photonPolarization + sinp * perp;
                                                   >> 248     
                                                   >> 249     photonPolarization = photonPolarization.unit();
                                                   >> 250     
197     // Generate a new photon:                     251     // Generate a new photon:
198     auto sec_dp =                              << 252     
199       new G4DynamicParticle(G4OpticalPhoton::O << 253     G4DynamicParticle* aWLSPhoton =
200     sec_dp->SetPolarization(photonPolarization << 254       new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(),
201     sec_dp->SetKineticEnergy(sampledEnergy);   << 255           photonMomentum);
202                                                << 256     aWLSPhoton->SetPolarization
203     G4double secTime = pPostStepPoint->GetGlob << 257       (photonPolarization.x(),
204                        WLSTimeGeneratorProfile << 258        photonPolarization.y(),
205     G4ThreeVector secPos = pPostStepPoint->Get << 259        photonPolarization.z());
206     G4Track* secTrack    = new G4Track(sec_dp, << 260     
207                                                << 261     aWLSPhoton->SetKineticEnergy(sampledEnergy);
208     secTrack->SetTouchableHandle(aTrack.GetTou << 262     
209     secTrack->SetParentID(aTrack.GetTrackID()) << 263     // Generate new G4Track object:
210                                                << 264     
211     proposedSecondaries.push_back(secTrack);   << 265     // Must give position of WLS optical photon
212   }                                            << 266 
213                                                << 267     G4double TimeDelay = WLSTimeGeneratorProfile->GenerateTime(WLSTime);
214   aParticleChange.SetNumberOfSecondaries((G4in << 268     G4double aSecondaryTime = (pPostStepPoint->GetGlobalTime()) + TimeDelay;
215   for(auto sec : proposedSecondaries)          << 269 
216   {                                            << 270     G4ThreeVector aSecondaryPosition = pPostStepPoint->GetPosition();
217     aParticleChange.AddSecondary(sec);         << 271 
218   }                                            << 272     G4Track* aSecondaryTrack = 
219   if(verboseLevel > 1)                         << 273       new G4Track(aWLSPhoton,aSecondaryTime,aSecondaryPosition);
220   {                                            << 274    
221     G4cout << "\n Exiting from G4OpWLS::DoIt - << 275     aSecondaryTrack->SetTouchableHandle(aTrack.GetTouchableHandle()); 
222            << aParticleChange.GetNumberOfSecon << 276     // aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0);
                                                   >> 277     
                                                   >> 278     aSecondaryTrack->SetParentID(aTrack.GetTrackID());
                                                   >> 279     
                                                   >> 280     aParticleChange.AddSecondary(aSecondaryTrack);
                                                   >> 281   }
                                                   >> 282 
                                                   >> 283   if (verboseLevel>0) {
                                                   >> 284     G4cout << "\n Exiting from G4OpWLS::DoIt -- NumberOfSecondaries = " 
                                                   >> 285      << aParticleChange.GetNumberOfSecondaries() << G4endl;  
223   }                                               286   }
224                                                << 287   
225   return G4VDiscreteProcess::PostStepDoIt(aTra    288   return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
226 }                                                 289 }
227                                                   290 
228 //....oooOO0OOooo........oooOO0OOooo........oo << 291 // BuildThePhysicsTable for the wavelength shifting process
229 void G4OpWLS::BuildPhysicsTable(const G4Partic << 292 // --------------------------------------------------
230 {                                              << 293 //
231   if(theIntegralTable)                         << 
232   {                                            << 
233     theIntegralTable->clearAndDestroy();       << 
234     delete theIntegralTable;                   << 
235     theIntegralTable = nullptr;                << 
236   }                                            << 
237                                                << 
238   const G4MaterialTable* materialTable = G4Mat << 
239   std::size_t numOfMaterials           = G4Mat << 
240   theIntegralTable                     = new G << 
241                                                   294 
                                                   >> 295 void G4OpWLS::BuildThePhysicsTable()
                                                   >> 296 {
                                                   >> 297   if (theIntegralTable) return;
                                                   >> 298   
                                                   >> 299   const G4MaterialTable* theMaterialTable = 
                                                   >> 300     G4Material::GetMaterialTable();
                                                   >> 301   G4int numOfMaterials = G4Material::GetNumberOfMaterials();
                                                   >> 302   
                                                   >> 303   // create new physics table
                                                   >> 304   
                                                   >> 305   if(!theIntegralTable)theIntegralTable = new G4PhysicsTable(numOfMaterials);
                                                   >> 306   
242   // loop for materials                           307   // loop for materials
243   for(std::size_t i = 0; i < numOfMaterials; + << 308   
244   {                                            << 309   for (G4int i=0 ; i < numOfMaterials; i++)
245     auto physVector = new G4PhysicsFreeVector( << 
246                                                << 
247     // Retrieve vector of WLS wavelength inten << 
248     // the material from the material's optica << 
249     G4MaterialPropertiesTable* MPT =           << 
250       (*materialTable)[i]->GetMaterialProperti << 
251     if(MPT)                                    << 
252     {                                             310     {
253       G4MaterialPropertyVector* wlsVector = MP << 311       G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
254       if(wlsVector)                            << 312   new G4PhysicsOrderedFreeVector();
255       {                                        << 313       
256         // Retrieve the first intensity point  << 314       // Retrieve vector of WLS wavelength intensity for
257         // of (photon energy, intensity) pairs << 315       // the material from the material's optical properties table.
258         G4double currentIN = (*wlsVector)[0];  << 316       
259         if(currentIN >= 0.0)                   << 317       G4Material* aMaterial = (*theMaterialTable)[i];
260         {                                      << 318 
261           // Create first (photon energy)      << 319       G4MaterialPropertiesTable* aMaterialPropertiesTable =
262           G4double currentPM  = wlsVector->Ene << 320   aMaterial->GetMaterialPropertiesTable();
263           G4double currentCII = 0.0;           << 321 
264           physVector->InsertValues(currentPM,  << 322       if (aMaterialPropertiesTable) {
265                                                << 323 
266           // Set previous values to current on << 324   G4MaterialPropertyVector* theWLSVector = 
267           G4double prevPM  = currentPM;        << 325     aMaterialPropertiesTable->GetProperty("WLSCOMPONENT");
268           G4double prevCII = currentCII;       << 326 
269           G4double prevIN  = currentIN;        << 327   if (theWLSVector) {
270                                                << 328     
271           // loop over all (photon energy, int << 329     // Retrieve the first intensity point in vector
272           // pairs stored for this material    << 330     // of (photon energy, intensity) pairs
273           for(std::size_t j = 1; j < wlsVector << 331     
274           {                                    << 332     G4double currentIN = (*theWLSVector)[0];
275             currentPM = wlsVector->Energy(j);  << 333     
276             currentIN = (*wlsVector)[j];       << 334     if (currentIN >= 0.0) {
277             currentCII =                       << 335 
278               prevCII + 0.5 * (currentPM - pre << 336       // Create first (photon energy) 
279                                                << 337      
280             physVector->InsertValues(currentPM << 338       G4double currentPM = theWLSVector->Energy(0);
281                                                << 339       
282             prevPM  = currentPM;               << 340       G4double currentCII = 0.0;
283             prevCII = currentCII;              << 341       
284             prevIN  = currentIN;               << 342       aPhysicsOrderedFreeVector->
285           }                                    << 343         InsertValues(currentPM , currentCII);
286         }                                      << 344       
                                                   >> 345       // Set previous values to current ones prior to loop
                                                   >> 346       
                                                   >> 347       G4double prevPM  = currentPM;
                                                   >> 348       G4double prevCII = currentCII;
                                                   >> 349       G4double prevIN  = currentIN;
                                                   >> 350       
                                                   >> 351       // loop over all (photon energy, intensity)
                                                   >> 352       // pairs stored for this material
                                                   >> 353 
                                                   >> 354             for (size_t j = 1;
                                                   >> 355                  j < theWLSVector->GetVectorLength();
                                                   >> 356                  j++)     
                                                   >> 357         {
                                                   >> 358     currentPM = theWLSVector->Energy(j);
                                                   >> 359     currentIN = (*theWLSVector)[j];
                                                   >> 360     
                                                   >> 361     currentCII = 0.5 * (prevIN + currentIN);
                                                   >> 362     
                                                   >> 363     currentCII = prevCII +
                                                   >> 364       (currentPM - prevPM) * currentCII;
                                                   >> 365     
                                                   >> 366     aPhysicsOrderedFreeVector->
                                                   >> 367       InsertValues(currentPM, currentCII);
                                                   >> 368     
                                                   >> 369     prevPM  = currentPM;
                                                   >> 370     prevCII = currentCII;
                                                   >> 371     prevIN  = currentIN;
                                                   >> 372         }
                                                   >> 373     }
                                                   >> 374   }
287       }                                           375       }
                                                   >> 376   // The WLS integral for a given material
                                                   >> 377   // will be inserted in the table according to the
                                                   >> 378   // position of the material in the material table.
                                                   >> 379 
                                                   >> 380   theIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
288     }                                             381     }
289     theIntegralTable->insertAt(i, physVector); << 
290   }                                            << 
291 }                                                 382 }
292                                                   383 
293 //....oooOO0OOooo........oooOO0OOooo........oo << 384 // GetMeanFreePath
294 G4double G4OpWLS::GetMeanFreePath(const G4Trac << 385 // ---------------
295                                   G4ForceCondi << 386 //
                                                   >> 387 G4double G4OpWLS::GetMeanFreePath(const G4Track& aTrack,
                                                   >> 388                  G4double ,
                                                   >> 389                  G4ForceCondition* )
296 {                                                 390 {
297   G4double thePhotonEnergy = aTrack.GetDynamic << 391   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
298   G4double attLength       = DBL_MAX;          << 392   const G4Material* aMaterial = aTrack.GetMaterial();
299   G4MaterialPropertiesTable* MPT =             << 393 
300     aTrack.GetMaterial()->GetMaterialPropertie << 394   G4double thePhotonEnergy = aParticle->GetTotalEnergy();
301                                                << 395 
302   if(MPT)                                      << 396   G4MaterialPropertiesTable* aMaterialPropertyTable;
303   {                                            << 397   G4MaterialPropertyVector* AttenuationLengthVector;
304     G4MaterialPropertyVector* attVector = MPT- << 398   
305     if(attVector)                              << 399   G4double AttenuationLength = DBL_MAX;
306     {                                          << 400 
307       attLength = attVector->Value(thePhotonEn << 401   aMaterialPropertyTable = aMaterial->GetMaterialPropertiesTable();
                                                   >> 402 
                                                   >> 403   if ( aMaterialPropertyTable ) {
                                                   >> 404     AttenuationLengthVector = aMaterialPropertyTable->
                                                   >> 405       GetProperty("WLSABSLENGTH");
                                                   >> 406     if ( AttenuationLengthVector ){
                                                   >> 407       AttenuationLength = AttenuationLengthVector->
                                                   >> 408   Value(thePhotonEnergy);
                                                   >> 409     }
                                                   >> 410     else {
                                                   >> 411       //             G4cout << "No WLS absorption length specified" << G4endl;
308     }                                             412     }
309   }                                               413   }
310   return attLength;                            << 414   else {
                                                   >> 415     //           G4cout << "No WLS absortion length specified" << G4endl;
                                                   >> 416   }
                                                   >> 417   
                                                   >> 418   return AttenuationLength;
311 }                                                 419 }
312                                                   420 
313 //....oooOO0OOooo........oooOO0OOooo........oo << 421 #include "G4AutoLock.hh"
314 void G4OpWLS::UseTimeProfile(const G4String na << 422 namespace {
315 {                                              << 423  G4Mutex wlsCmdHandlingMutex = G4MUTEX_INITIALIZER;
316   if(WLSTimeGeneratorProfile)                  << 
317   {                                            << 
318     delete WLSTimeGeneratorProfile;            << 
319     WLSTimeGeneratorProfile = nullptr;         << 
320   }                                            << 
321   if(name == "delta")                          << 
322   {                                            << 
323     WLSTimeGeneratorProfile = new G4WLSTimeGen << 
324   }                                            << 
325   else if(name == "exponential")               << 
326   {                                            << 
327     WLSTimeGeneratorProfile =                  << 
328       new G4WLSTimeGeneratorProfileExponential << 
329   }                                            << 
330   else                                         << 
331   {                                            << 
332     G4Exception("G4OpWLS::UseTimeProfile", "em << 
333                 "generator does not exist");   << 
334   }                                            << 
335   G4OpticalParameters::Instance()->SetWLSTimeP << 
336 }                                                 424 }
337                                                   425 
338 //....oooOO0OOooo........oooOO0OOooo........oo << 426 void G4OpWLS::UseTimeProfile(const G4String name)
339 void G4OpWLS::SetVerboseLevel(G4int verbose)   << 
340 {                                                 427 {
341   verboseLevel = verbose;                      << 428   G4AutoLock l(&wlsCmdHandlingMutex);
342   G4OpticalParameters::Instance()->SetWLSVerbo << 429   if (name == "delta")
                                                   >> 430     {
                                                   >> 431       delete WLSTimeGeneratorProfile;
                                                   >> 432       WLSTimeGeneratorProfile = 
                                                   >> 433              new G4WLSTimeGeneratorProfileDelta("delta");
                                                   >> 434     }
                                                   >> 435   else if (name == "exponential")
                                                   >> 436     {
                                                   >> 437       delete WLSTimeGeneratorProfile;
                                                   >> 438       WLSTimeGeneratorProfile =
                                                   >> 439              new G4WLSTimeGeneratorProfileExponential("exponential");
                                                   >> 440     }
                                                   >> 441   else
                                                   >> 442     {
                                                   >> 443       G4Exception("G4OpWLS::UseTimeProfile", "em0202",
                                                   >> 444                   FatalException,
                                                   >> 445                   "generator does not exist");
                                                   >> 446     }
343 }                                                 447 }
344                                                   448