| Geant4 Cross Reference |
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 /// \file optical/wls/src/WLSPrimaryGeneratorA 27 /// \file optical/wls/src/WLSPrimaryGeneratorAction.cc
28 /// \brief Implementation of the WLSPrimaryGen 28 /// \brief Implementation of the WLSPrimaryGeneratorAction class
29 // 29 //
30 // 30 //
31 31
32 #include "WLSPrimaryGeneratorAction.hh" 32 #include "WLSPrimaryGeneratorAction.hh"
33 33
34 #include "WLSDetectorConstruction.hh" 34 #include "WLSDetectorConstruction.hh"
35 #include "WLSPrimaryGeneratorMessenger.hh" 35 #include "WLSPrimaryGeneratorMessenger.hh"
36 36
37 #include "G4AutoLock.hh" 37 #include "G4AutoLock.hh"
38 #include "G4Event.hh" 38 #include "G4Event.hh"
39 #include "G4GeneralParticleSource.hh" 39 #include "G4GeneralParticleSource.hh"
40 #include "G4Material.hh" 40 #include "G4Material.hh"
41 #include "G4MaterialPropertiesTable.hh" 41 #include "G4MaterialPropertiesTable.hh"
42 #include "G4OpticalPhoton.hh" 42 #include "G4OpticalPhoton.hh"
43 #include "G4PhysicsTable.hh" 43 #include "G4PhysicsTable.hh"
44 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh"
45 #include "G4UImanager.hh" 45 #include "G4UImanager.hh"
46 #include "Randomize.hh" 46 #include "Randomize.hh"
47 47
48 namespace 48 namespace
49 { 49 {
50 G4Mutex gen_mutex = G4MUTEX_INITIALIZER; << 50 G4Mutex gen_mutex = G4MUTEX_INITIALIZER;
51 } 51 }
52 52
53 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 54
55 G4bool WLSPrimaryGeneratorAction::fFirst = fal 55 G4bool WLSPrimaryGeneratorAction::fFirst = false;
56 56
57 WLSPrimaryGeneratorAction::WLSPrimaryGenerator << 57 WLSPrimaryGeneratorAction::WLSPrimaryGeneratorAction(
>> 58 WLSDetectorConstruction* dc)
58 { 59 {
59 fDetector = dc; << 60 fDetector = dc;
>> 61 fIntegralTable = nullptr;
60 62
61 fParticleGun = new G4GeneralParticleSource() 63 fParticleGun = new G4GeneralParticleSource();
>> 64
62 fGunMessenger = new WLSPrimaryGeneratorMesse 65 fGunMessenger = new WLSPrimaryGeneratorMessenger(this);
>> 66
>> 67 fTimeConstant = 0.;
>> 68 fUseSampledEnergy = false;
63 } 69 }
64 70
65 //....oooOO0OOooo........oooOO0OOooo........oo 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
66 72
67 WLSPrimaryGeneratorAction::~WLSPrimaryGenerato 73 WLSPrimaryGeneratorAction::~WLSPrimaryGeneratorAction()
68 { 74 {
69 delete fParticleGun; 75 delete fParticleGun;
70 delete fGunMessenger; 76 delete fGunMessenger;
71 if (fIntegralTable) { << 77 if(fIntegralTable)
>> 78 {
72 fIntegralTable->clearAndDestroy(); 79 fIntegralTable->clearAndDestroy();
73 delete fIntegralTable; 80 delete fIntegralTable;
74 } 81 }
75 } 82 }
76 83
77 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 85
79 void WLSPrimaryGeneratorAction::SetDecayTimeCo 86 void WLSPrimaryGeneratorAction::SetDecayTimeConstant(G4double time)
80 { 87 {
81 fTimeConstant = time; 88 fTimeConstant = time;
82 } 89 }
83 90
84 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85 92
86 void WLSPrimaryGeneratorAction::BuildEmissionS 93 void WLSPrimaryGeneratorAction::BuildEmissionSpectrum()
87 { 94 {
88 if (fIntegralTable) return; << 95 if(fIntegralTable)
>> 96 return;
89 97
90 const G4MaterialTable* theMaterialTable = G4 98 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
91 G4int numOfMaterials = G4Material::GetNumber << 99 G4int numOfMaterials = G4Material::GetNumberOfMaterials();
92 100
93 if (!fIntegralTable) fIntegralTable = new G4 << 101 if(!fIntegralTable)
>> 102 fIntegralTable = new G4PhysicsTable(numOfMaterials);
94 103
95 for (G4int i = 0; i < numOfMaterials; ++i) { << 104 for(G4int i = 0; i < numOfMaterials; ++i)
96 auto vec = new G4PhysicsFreeVector(); << 105 {
>> 106 G4PhysicsFreeVector* vec = new G4PhysicsFreeVector();
97 107
98 G4MaterialPropertiesTable* MPT = (*theMate << 108 G4MaterialPropertiesTable* MPT =
>> 109 (*theMaterialTable)[i]->GetMaterialPropertiesTable();
99 110
100 if (MPT) { << 111 if(MPT)
>> 112 {
101 G4MaterialPropertyVector* theWLSVector = 113 G4MaterialPropertyVector* theWLSVector = MPT->GetProperty("WLSCOMPONENT");
102 114
103 if (theWLSVector) { << 115 if(theWLSVector)
>> 116 {
104 G4double currentIN = (*theWLSVector)[0 117 G4double currentIN = (*theWLSVector)[0];
105 if (currentIN >= 0.0) { << 118 if(currentIN >= 0.0)
106 G4double currentPM = theWLSVector->E << 119 {
>> 120 G4double currentPM = theWLSVector->Energy(0);
107 G4double currentCII = 0.0; 121 G4double currentCII = 0.0;
108 vec->InsertValues(currentPM, current 122 vec->InsertValues(currentPM, currentCII);
109 G4double prevPM = currentPM; << 123 G4double prevPM = currentPM;
110 G4double prevCII = currentCII; 124 G4double prevCII = currentCII;
111 G4double prevIN = currentIN; << 125 G4double prevIN = currentIN;
112 126
113 for (size_t j = 1; j < theWLSVector- << 127 for(size_t j = 1; j < theWLSVector->GetVectorLength(); ++j)
114 currentPM = theWLSVector->Energy(j << 128 {
115 currentIN = (*theWLSVector)[j]; << 129 currentPM = theWLSVector->Energy(j);
>> 130 currentIN = (*theWLSVector)[j];
116 currentCII = 0.5 * (prevIN + curre 131 currentCII = 0.5 * (prevIN + currentIN);
117 currentCII = prevCII + (currentPM 132 currentCII = prevCII + (currentPM - prevPM) * currentCII;
118 vec->InsertValues(currentPM, curre 133 vec->InsertValues(currentPM, currentCII);
119 prevPM = currentPM; << 134 prevPM = currentPM;
120 prevCII = currentCII; 135 prevCII = currentCII;
121 prevIN = currentIN; << 136 prevIN = currentIN;
122 } 137 }
123 } 138 }
124 } 139 }
125 } 140 }
126 fIntegralTable->insertAt(i, vec); 141 fIntegralTable->insertAt(i, vec);
127 } 142 }
128 } 143 }
129 144
130 //....oooOO0OOooo........oooOO0OOooo........oo 145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131 146
132 void WLSPrimaryGeneratorAction::GeneratePrimar 147 void WLSPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
133 { 148 {
134 if (!fFirst) { << 149 if(!fFirst)
>> 150 {
135 fFirst = true; 151 fFirst = true;
136 BuildEmissionSpectrum(); 152 BuildEmissionSpectrum();
137 } 153 }
138 154
139 if (fUseSampledEnergy) { << 155 if(fUseSampledEnergy)
>> 156 {
140 const G4MaterialTable* theMaterialTable = 157 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
141 158
142 G4double sampledEnergy = 3. * eV; 159 G4double sampledEnergy = 3. * eV;
143 160
144 for (size_t j = 0; j < theMaterialTable->s << 161 for(size_t j = 0; j < theMaterialTable->size(); ++j)
>> 162 {
145 G4Material* fMaterial = (*theMaterialTab 163 G4Material* fMaterial = (*theMaterialTable)[j];
146 if (fMaterial->GetName() == "PMMA") { << 164 if(fMaterial->GetName() == "PMMA")
147 auto WLSIntensity = fMaterial->GetMate << 165 {
148 << 166 auto WLSIntensity =
149 if (WLSIntensity) { << 167 fMaterial->GetMaterialPropertiesTable()->GetProperty("WLSCOMPONENT");
150 auto WLSIntegral = (G4PhysicsFreeVec << 168
>> 169 if(WLSIntensity)
>> 170 {
>> 171 G4PhysicsFreeVector* WLSIntegral =
>> 172 (G4PhysicsFreeVector*) ((*fIntegralTable)(
>> 173 fMaterial->GetIndex()));
151 174
152 G4double CIImax = WLSIntegral->GetMa << 175 G4double CIImax = WLSIntegral->GetMaxValue();
153 G4double CIIvalue = G4UniformRand() 176 G4double CIIvalue = G4UniformRand() * CIImax;
154 177
155 sampledEnergy = WLSIntegral->GetEner 178 sampledEnergy = WLSIntegral->GetEnergy(CIIvalue);
156 } 179 }
157 } 180 }
158 } 181 }
159 182
160 // this does not work. 183 // this does not work.
161 G4String cmd = "/gun/energy " + G4UIcomman 184 G4String cmd = "/gun/energy " + G4UIcommand::ConvertToString(sampledEnergy / eV) + " eV";
162 G4UImanager::GetUIpointer()->ApplyCommand( 185 G4UImanager::GetUIpointer()->ApplyCommand(cmd);
163 } 186 }
164 187
165 // The code behind this line is not thread s 188 // The code behind this line is not thread safe because polarization
166 // and time are randomly selected and GPS pr 189 // and time are randomly selected and GPS properties are global
167 190
168 G4AutoLock l(&gen_mutex); 191 G4AutoLock l(&gen_mutex);
169 if (fParticleGun->GetParticleDefinition() == << 192 if(fParticleGun->GetParticleDefinition() == G4OpticalPhoton::Definition())
>> 193 {
170 SetOptPhotonPolar(); 194 SetOptPhotonPolar();
171 SetOptPhotonTime(); 195 SetOptPhotonTime();
172 } 196 }
173 197
174 fParticleGun->GeneratePrimaryVertex(anEvent) 198 fParticleGun->GeneratePrimaryVertex(anEvent);
175 } 199 }
176 200
177 //....oooOO0OOooo........oooOO0OOooo........oo 201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
178 202
179 void WLSPrimaryGeneratorAction::SetOptPhotonPo 203 void WLSPrimaryGeneratorAction::SetOptPhotonPolar()
180 { 204 {
181 G4double angle = G4UniformRand() * 360.0 * d 205 G4double angle = G4UniformRand() * 360.0 * deg;
182 SetOptPhotonPolar(angle); 206 SetOptPhotonPolar(angle);
183 } 207 }
184 208
185 //....oooOO0OOooo........oooOO0OOooo........oo 209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
186 210
187 void WLSPrimaryGeneratorAction::SetOptPhotonPo 211 void WLSPrimaryGeneratorAction::SetOptPhotonPolar(G4double angle)
188 { 212 {
189 if (fParticleGun->GetParticleDefinition()->G << 213 if(fParticleGun->GetParticleDefinition()->GetParticleName() !=
>> 214 "opticalphoton")
>> 215 {
190 G4cout << "-> warning from WLSPrimaryGener 216 G4cout << "-> warning from WLSPrimaryGeneratorAction::SetOptPhotonPolar()"
191 << ": the ParticleGun is not an op 217 << ": the ParticleGun is not an opticalphoton" << G4endl;
192 return; 218 return;
193 } 219 }
194 220
195 G4ThreeVector normal(1., 0., 0.); 221 G4ThreeVector normal(1., 0., 0.);
196 G4ThreeVector kphoton = fParticleGun->GetPar 222 G4ThreeVector kphoton = fParticleGun->GetParticleMomentumDirection();
197 G4ThreeVector product = normal.cross(kphoton 223 G4ThreeVector product = normal.cross(kphoton);
198 G4double modul2 = product * product; << 224 G4double modul2 = product * product;
199 225
200 G4ThreeVector e_perpend(0., 0., 1.); 226 G4ThreeVector e_perpend(0., 0., 1.);
201 if (modul2 > 0.) e_perpend = (1. / std::sqrt << 227 if(modul2 > 0.)
>> 228 e_perpend = (1. / std::sqrt(modul2)) * product;
202 G4ThreeVector e_paralle = e_perpend.cross(kp 229 G4ThreeVector e_paralle = e_perpend.cross(kphoton);
203 230
204 G4ThreeVector polar = std::cos(angle) * e_pa << 231 G4ThreeVector polar =
>> 232 std::cos(angle) * e_paralle + std::sin(angle) * e_perpend;
205 fParticleGun->SetParticlePolarization(polar) 233 fParticleGun->SetParticlePolarization(polar);
206 } 234 }
207 235
208 //....oooOO0OOooo........oooOO0OOooo........oo 236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
209 237
210 void WLSPrimaryGeneratorAction::SetOptPhotonTi 238 void WLSPrimaryGeneratorAction::SetOptPhotonTime()
211 { 239 {
212 G4double time = -std::log(G4UniformRand()) * 240 G4double time = -std::log(G4UniformRand()) * fTimeConstant;
213 fParticleGun->SetParticleTime(time); 241 fParticleGun->SetParticleTime(time);
214 } 242 }
215 243