00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00030
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00045
00046 #include "G4OpWLS.hh"
00047
00048 #include "G4ios.hh"
00049 #include "G4PhysicalConstants.hh"
00050 #include "G4SystemOfUnits.hh"
00051 #include "G4OpProcessSubType.hh"
00052
00053 #include "G4WLSTimeGeneratorProfileDelta.hh"
00054 #include "G4WLSTimeGeneratorProfileExponential.hh"
00055
00057
00059
00061
00063
00064 G4OpWLS::G4OpWLS(const G4String& processName, G4ProcessType type)
00065 : G4VDiscreteProcess(processName, type)
00066 {
00067 SetProcessSubType(fOpWLS);
00068
00069 theIntegralTable = 0;
00070
00071 if (verboseLevel>0) {
00072 G4cout << GetProcessName() << " is created " << G4endl;
00073 }
00074
00075 WLSTimeGeneratorProfile =
00076 new G4WLSTimeGeneratorProfileDelta("WLSTimeGeneratorProfileDelta");
00077
00078 BuildThePhysicsTable();
00079 }
00080
00082
00084
00085 G4OpWLS::~G4OpWLS()
00086 {
00087 if (theIntegralTable != 0) {
00088 theIntegralTable->clearAndDestroy();
00089 delete theIntegralTable;
00090 }
00091 delete WLSTimeGeneratorProfile;
00092 }
00093
00095
00097
00098
00099
00100
00101 G4VParticleChange*
00102 G4OpWLS::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
00103 {
00104 aParticleChange.Initialize(aTrack);
00105
00106 aParticleChange.ProposeTrackStatus(fStopAndKill);
00107
00108 if (verboseLevel>0) {
00109 G4cout << "\n** Photon absorbed! **" << G4endl;
00110 }
00111
00112 const G4Material* aMaterial = aTrack.GetMaterial();
00113
00114 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
00115
00116 G4MaterialPropertiesTable* aMaterialPropertiesTable =
00117 aMaterial->GetMaterialPropertiesTable();
00118 if (!aMaterialPropertiesTable)
00119 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00120
00121 const G4MaterialPropertyVector* WLS_Intensity =
00122 aMaterialPropertiesTable->GetProperty("WLSCOMPONENT");
00123
00124 if (!WLS_Intensity)
00125 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00126
00127 G4int NumPhotons = 1;
00128
00129 if (aMaterialPropertiesTable->ConstPropertyExists("WLSMEANNUMBERPHOTONS")) {
00130
00131 G4double MeanNumberOfPhotons = aMaterialPropertiesTable->
00132 GetConstProperty("WLSMEANNUMBERPHOTONS");
00133
00134 NumPhotons = G4int(G4Poisson(MeanNumberOfPhotons));
00135
00136 if (NumPhotons <= 0) {
00137
00138
00139
00140 aParticleChange.SetNumberOfSecondaries(0);
00141
00142 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00143
00144 }
00145
00146 }
00147
00148 aParticleChange.SetNumberOfSecondaries(NumPhotons);
00149
00150 G4int materialIndex = aMaterial->GetIndex();
00151
00152
00153
00154
00155 G4double WLSTime = 0.*ns;
00156 G4PhysicsOrderedFreeVector* WLSIntegral = 0;
00157
00158 WLSTime = aMaterialPropertiesTable->
00159 GetConstProperty("WLSTIMECONSTANT");
00160 WLSIntegral =
00161 (G4PhysicsOrderedFreeVector*)((*theIntegralTable)(materialIndex));
00162
00163
00164
00165 G4double CIImax = WLSIntegral->GetMaxValue();
00166
00167 for (G4int i = 0; i < NumPhotons; i++) {
00168
00169
00170
00171 G4double CIIvalue = G4UniformRand()*CIImax;
00172 G4double sampledEnergy =
00173 WLSIntegral->GetEnergy(CIIvalue);
00174
00175 if (verboseLevel>1) {
00176 G4cout << "sampledEnergy = " << sampledEnergy << G4endl;
00177 G4cout << "CIIvalue = " << CIIvalue << G4endl;
00178 }
00179
00180
00181
00182 G4double cost = 1. - 2.*G4UniformRand();
00183 G4double sint = std::sqrt((1.-cost)*(1.+cost));
00184
00185 G4double phi = twopi*G4UniformRand();
00186 G4double sinp = std::sin(phi);
00187 G4double cosp = std::cos(phi);
00188
00189 G4double px = sint*cosp;
00190 G4double py = sint*sinp;
00191 G4double pz = cost;
00192
00193
00194
00195 G4ParticleMomentum photonMomentum(px, py, pz);
00196
00197
00198
00199 G4double sx = cost*cosp;
00200 G4double sy = cost*sinp;
00201 G4double sz = -sint;
00202
00203 G4ThreeVector photonPolarization(sx, sy, sz);
00204
00205 G4ThreeVector perp = photonMomentum.cross(photonPolarization);
00206
00207 phi = twopi*G4UniformRand();
00208 sinp = std::sin(phi);
00209 cosp = std::cos(phi);
00210
00211 photonPolarization = cosp * photonPolarization + sinp * perp;
00212
00213 photonPolarization = photonPolarization.unit();
00214
00215
00216
00217 G4DynamicParticle* aWLSPhoton =
00218 new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(),
00219 photonMomentum);
00220 aWLSPhoton->SetPolarization
00221 (photonPolarization.x(),
00222 photonPolarization.y(),
00223 photonPolarization.z());
00224
00225 aWLSPhoton->SetKineticEnergy(sampledEnergy);
00226
00227
00228
00229
00230
00231 G4double TimeDelay = WLSTimeGeneratorProfile->GenerateTime(WLSTime);
00232 G4double aSecondaryTime = (pPostStepPoint->GetGlobalTime()) + TimeDelay;
00233
00234 G4ThreeVector aSecondaryPosition = pPostStepPoint->GetPosition();
00235
00236 G4Track* aSecondaryTrack =
00237 new G4Track(aWLSPhoton,aSecondaryTime,aSecondaryPosition);
00238
00239 aSecondaryTrack->SetTouchableHandle(aTrack.GetTouchableHandle());
00240
00241
00242 aSecondaryTrack->SetParentID(aTrack.GetTrackID());
00243
00244 aParticleChange.AddSecondary(aSecondaryTrack);
00245 }
00246
00247 if (verboseLevel>0) {
00248 G4cout << "\n Exiting from G4OpWLS::DoIt -- NumberOfSecondaries = "
00249 << aParticleChange.GetNumberOfSecondaries() << G4endl;
00250 }
00251
00252 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
00253 }
00254
00255
00256
00257
00258
00259 void G4OpWLS::BuildThePhysicsTable()
00260 {
00261 if (theIntegralTable) return;
00262
00263 const G4MaterialTable* theMaterialTable =
00264 G4Material::GetMaterialTable();
00265 G4int numOfMaterials = G4Material::GetNumberOfMaterials();
00266
00267
00268
00269 if(!theIntegralTable)theIntegralTable = new G4PhysicsTable(numOfMaterials);
00270
00271
00272
00273 for (G4int i=0 ; i < numOfMaterials; i++)
00274 {
00275 G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
00276 new G4PhysicsOrderedFreeVector();
00277
00278
00279
00280
00281 G4Material* aMaterial = (*theMaterialTable)[i];
00282
00283 G4MaterialPropertiesTable* aMaterialPropertiesTable =
00284 aMaterial->GetMaterialPropertiesTable();
00285
00286 if (aMaterialPropertiesTable) {
00287
00288 G4MaterialPropertyVector* theWLSVector =
00289 aMaterialPropertiesTable->GetProperty("WLSCOMPONENT");
00290
00291 if (theWLSVector) {
00292
00293
00294
00295
00296 G4double currentIN = (*theWLSVector)[0];
00297
00298 if (currentIN >= 0.0) {
00299
00300
00301
00302 G4double currentPM = theWLSVector->Energy(0);
00303
00304 G4double currentCII = 0.0;
00305
00306 aPhysicsOrderedFreeVector->
00307 InsertValues(currentPM , currentCII);
00308
00309
00310
00311 G4double prevPM = currentPM;
00312 G4double prevCII = currentCII;
00313 G4double prevIN = currentIN;
00314
00315
00316
00317
00318 for (size_t j = 1;
00319 j < theWLSVector->GetVectorLength();
00320 j++)
00321 {
00322 currentPM = theWLSVector->Energy(j);
00323 currentIN = (*theWLSVector)[j];
00324
00325 currentCII = 0.5 * (prevIN + currentIN);
00326
00327 currentCII = prevCII +
00328 (currentPM - prevPM) * currentCII;
00329
00330 aPhysicsOrderedFreeVector->
00331 InsertValues(currentPM, currentCII);
00332
00333 prevPM = currentPM;
00334 prevCII = currentCII;
00335 prevIN = currentIN;
00336 }
00337 }
00338 }
00339 }
00340
00341
00342
00343
00344 theIntegralTable->insertAt(i,aPhysicsOrderedFreeVector);
00345 }
00346 }
00347
00348
00349
00350
00351 G4double G4OpWLS::GetMeanFreePath(const G4Track& aTrack,
00352 G4double ,
00353 G4ForceCondition* )
00354 {
00355 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00356 const G4Material* aMaterial = aTrack.GetMaterial();
00357
00358 G4double thePhotonEnergy = aParticle->GetTotalEnergy();
00359
00360 G4MaterialPropertiesTable* aMaterialPropertyTable;
00361 G4MaterialPropertyVector* AttenuationLengthVector;
00362
00363 G4double AttenuationLength = DBL_MAX;
00364
00365 aMaterialPropertyTable = aMaterial->GetMaterialPropertiesTable();
00366
00367 if ( aMaterialPropertyTable ) {
00368 AttenuationLengthVector = aMaterialPropertyTable->
00369 GetProperty("WLSABSLENGTH");
00370 if ( AttenuationLengthVector ){
00371 AttenuationLength = AttenuationLengthVector->
00372 Value(thePhotonEnergy);
00373 }
00374 else {
00375
00376 }
00377 }
00378 else {
00379
00380 }
00381
00382 return AttenuationLength;
00383 }
00384
00385 void G4OpWLS::UseTimeProfile(const G4String name)
00386 {
00387 if (name == "delta")
00388 {
00389 delete WLSTimeGeneratorProfile;
00390 WLSTimeGeneratorProfile =
00391 new G4WLSTimeGeneratorProfileDelta("delta");
00392 }
00393 else if (name == "exponential")
00394 {
00395 delete WLSTimeGeneratorProfile;
00396 WLSTimeGeneratorProfile =
00397 new G4WLSTimeGeneratorProfileExponential("exponential");
00398 }
00399 else
00400 {
00401 G4Exception("G4OpWLS::UseTimeProfile", "em0202",
00402 FatalException,
00403 "generator does not exist");
00404 }
00405 }