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