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