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