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 // $Id: GammaRayTelPrimaryGeneratorAction.cc 82630 2014-07-01 09:43:00Z gcosmo $ 27 // ------------------------------------------- 28 // ------------------------------------------------------------ 28 // GEANT 4 class implementation file 29 // GEANT 4 class implementation file 29 // CERN Geneva Switzerland 30 // CERN Geneva Switzerland 30 // 31 // 31 // 32 // 32 // ------------ GammaRayTelPrimaryGenerat 33 // ------------ GammaRayTelPrimaryGeneratorAction ------ 33 // by G.Santin, F.Longo & R.Giannit 34 // by G.Santin, F.Longo & R.Giannitrapani (13 nov 2000) 34 // 35 // 35 // ******************************************* 36 // ************************************************************ 36 37 >> 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 40 37 #include "G4RunManager.hh" 41 #include "G4RunManager.hh" 38 #include "GammaRayTelPrimaryGeneratorAction.hh 42 #include "GammaRayTelPrimaryGeneratorAction.hh" 39 43 40 #include "GammaRayTelDetectorConstruction.hh" 44 #include "GammaRayTelDetectorConstruction.hh" 41 #include "GammaRayTelPrimaryGeneratorMessenger 45 #include "GammaRayTelPrimaryGeneratorMessenger.hh" 42 46 43 #include "G4PhysicalConstants.hh" 47 #include "G4PhysicalConstants.hh" 44 #include "G4SystemOfUnits.hh" 48 #include "G4SystemOfUnits.hh" 45 #include "G4Event.hh" 49 #include "G4Event.hh" 46 #include "G4ParticleGun.hh" 50 #include "G4ParticleGun.hh" 47 #include "G4GeneralParticleSource.hh" 51 #include "G4GeneralParticleSource.hh" 48 #include "G4ParticleTable.hh" 52 #include "G4ParticleTable.hh" 49 #include "G4ParticleDefinition.hh" 53 #include "G4ParticleDefinition.hh" 50 #include "Randomize.hh" 54 #include "Randomize.hh" 51 55 52 //....oooOO0OOooo........oooOO0OOooo........oo 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 53 57 54 GammaRayTelPrimaryGeneratorAction::GammaRayTel << 58 GammaRayTelPrimaryGeneratorAction::GammaRayTelPrimaryGeneratorAction() 55 detector = static_cast<const GammaRayTelDe << 59 :rndmFlag("off"),nSourceType(0),nSpectrumType(0),sourceGun(false) 56 << 60 { 57 // create a messenger for this class << 61 GammaRayTelDetector = static_cast<const GammaRayTelDetectorConstruction*> 58 << 62 (G4RunManager::GetRunManager()->GetUserDetectorConstruction()); 59 gunMessenger = new GammaRayTelPrimaryGener << 63 60 << 64 //create a messenger for this class 61 constexpr auto NUMBER_OF_PARTICLES{1}; << 65 62 particleGun = new G4ParticleGun(NUMBER_OF_ << 66 gunMessenger = new GammaRayTelPrimaryGeneratorMessenger(this); 63 << 67 64 // default particle kinematic << 68 G4int n_particle = 1; >> 69 >> 70 particleGun = new G4ParticleGun(n_particle); >> 71 // default particle kinematic >> 72 >> 73 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable(); >> 74 G4String particleName; >> 75 G4ParticleDefinition* particle >> 76 = particleTable->FindParticle(particleName="e-"); >> 77 particleGun->SetParticleDefinition(particle); >> 78 particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,-1.)); >> 79 particleGun->SetParticleEnergy(30.*MeV); >> 80 G4double position = 0.5*(GammaRayTelDetector->GetWorldSizeZ()); >> 81 particleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,position)); >> 82 particleSource = new G4GeneralParticleSource(); >> 83 >> 84 } 65 85 66 auto *particleTable = G4ParticleTable::Get << 86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 67 auto *particle = particleTable->FindPartic << 68 particleGun->SetParticleDefinition(particl << 69 particleGun->SetParticleMomentumDirection( << 70 87 71 constexpr auto PARTICLE_ENERGY{30. * MeV}; << 88 GammaRayTelPrimaryGeneratorAction::~GammaRayTelPrimaryGeneratorAction() 72 particleGun->SetParticleEnergy(PARTICLE_EN << 89 { >> 90 >> 91 delete particleGun; >> 92 delete particleSource; 73 93 74 auto position = 0.5 * (detector->GetWorldS << 94 delete gunMessenger; 75 particleGun->SetParticlePosition(G4ThreeVe << 76 particleSource = new G4GeneralParticleSour << 77 } 95 } 78 96 79 //....oooOO0OOooo........oooOO0OOooo........oo 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 80 98 81 GammaRayTelPrimaryGeneratorAction::~GammaRayTe << 99 void GammaRayTelPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) 82 delete particleGun; << 100 { 83 delete particleSource; << 101 if (sourceGun) 84 delete gunMessenger; << 102 { >> 103 >> 104 //this function is called at the begining of event >> 105 // >> 106 G4double z0 = 0.5*(GammaRayTelDetector->GetWorldSizeZ()); >> 107 G4double x0 = 0.*cm, y0 = 0.*cm; >> 108 >> 109 G4ThreeVector pos0; >> 110 G4ThreeVector dir0; >> 111 G4ThreeVector vertex0 = G4ThreeVector(x0,y0,z0); >> 112 >> 113 dir0 = G4ThreeVector(0.,0.,-1.); >> 114 >> 115 G4double theta, phi, y, f; >> 116 G4double theta0=0.; >> 117 G4double phi0=0.; >> 118 >> 119 switch(nSourceType) { >> 120 case 0: >> 121 particleGun->SetParticlePosition(vertex0); >> 122 particleGun->SetParticleMomentumDirection(dir0); >> 123 break; >> 124 case 1: >> 125 // GS: Generate random position on the 4PIsphere to create a unif. distrib. >> 126 // GS: on the sphere >> 127 phi = G4UniformRand() * twopi; >> 128 do { >> 129 y = G4UniformRand()*1.0; >> 130 theta = G4UniformRand() * pi; >> 131 f = std::sin(theta); >> 132 } while (y > f); >> 133 vertex0 = G4ThreeVector(1.,0.,0.); >> 134 vertex0.setMag(dVertexRadius); >> 135 vertex0.setTheta(theta); >> 136 vertex0.setPhi(phi); >> 137 particleGun->SetParticlePosition(vertex0); >> 138 >> 139 dir0 = G4ThreeVector(1.,0.,0.); >> 140 do { >> 141 phi = G4UniformRand() * twopi; >> 142 do { >> 143 y = G4UniformRand()*1.0; >> 144 theta = G4UniformRand() * pi; >> 145 f = std::sin(theta); >> 146 } while (y > f); >> 147 dir0.setPhi(phi); >> 148 dir0.setTheta(theta); >> 149 } while (vertex0.dot(dir0) >= -0.7 * vertex0.mag()); >> 150 particleGun->SetParticleMomentumDirection((G4ParticleMomentum)dir0); >> 151 >> 152 break; >> 153 case 2: >> 154 // GS: Generate random position on the upper semi-sphere z>0 to create a unif. distrib. >> 155 // GS: on a plane >> 156 phi = G4UniformRand() * twopi; >> 157 do { >> 158 y = G4UniformRand()*1.0; >> 159 theta = G4UniformRand() * halfpi; >> 160 f = std::sin(theta) * std::cos(theta); >> 161 } while (y > f); >> 162 vertex0 = G4ThreeVector(1.,0.,0.); >> 163 >> 164 G4double xy = GammaRayTelDetector->GetWorldSizeXY(); >> 165 G4double z = GammaRayTelDetector->GetWorldSizeZ(); >> 166 >> 167 if (dVertexRadius > xy*0.5) >> 168 { >> 169 G4cout << "vertexRadius too big " << G4endl; >> 170 G4cout << "vertexRadius setted to " << xy*0.45 << G4endl; >> 171 dVertexRadius = xy*0.45; >> 172 } >> 173 >> 174 if (dVertexRadius > z*0.5) >> 175 { >> 176 G4cout << "vertexRadius too high " << G4endl; >> 177 G4cout << "vertexRadius setted to " << z*0.45 << G4endl; >> 178 dVertexRadius = z*0.45; >> 179 } >> 180 >> 181 >> 182 vertex0.setMag(dVertexRadius); >> 183 vertex0.setTheta(theta); >> 184 vertex0.setPhi(phi); >> 185 >> 186 // GS: Get the user defined direction for the primaries and >> 187 // GS: Rotate the random position according to the user defined direction for the particle >> 188 >> 189 dir0 = particleGun->GetParticleMomentumDirection(); >> 190 if (dir0.mag() > 0.001) >> 191 { >> 192 theta0 = dir0.theta(); >> 193 phi0 = dir0.phi(); >> 194 } >> 195 >> 196 if (theta0!=0.) >> 197 { >> 198 G4ThreeVector rotationAxis(1.,0.,0.); >> 199 rotationAxis.setPhi(phi0+halfpi); >> 200 vertex0.rotate(theta0+pi,rotationAxis); >> 201 } >> 202 particleGun->SetParticlePosition(vertex0); >> 203 break; >> 204 } >> 205 >> 206 >> 207 G4double pEnergy; >> 208 >> 209 switch(nSpectrumType) { >> 210 case 0: >> 211 break; >> 212 case 1: >> 213 break; >> 214 case 2: >> 215 do { >> 216 y = G4UniformRand()*100000.0; >> 217 pEnergy = G4UniformRand() * 10. * GeV; >> 218 f = std::pow(pEnergy * (1/GeV), -4.); >> 219 } while (y > f); >> 220 >> 221 particleGun->SetParticleEnergy(pEnergy); >> 222 >> 223 break; >> 224 case 3: >> 225 break; >> 226 } >> 227 >> 228 particleGun->GeneratePrimaryVertex(anEvent); >> 229 } >> 230 else >> 231 { >> 232 particleSource->GeneratePrimaryVertex(anEvent); >> 233 } >> 234 85 } 235 } 86 236 87 //....oooOO0OOooo........oooOO0OOooo........oo 237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 88 238 89 void GammaRayTelPrimaryGeneratorAction::Genera << 239 90 if (sourceGun) { << 240 91 G4cout << "Using G4ParticleGun... " << << 241 92 << 242 93 // this function is called at the begi << 243 94 // << 244 95 G4double x0 = 0. * cm; << 245 96 G4double y0 = 0. * cm; << 97 G4double z0 = 0.5 * (detector->GetWorl << 98 << 99 G4ThreeVector pos0; << 100 auto vertex0 = G4ThreeVector(x0, y0, z << 101 auto momentumDirection0 = G4ThreeVecto << 102 << 103 G4double theta; << 104 G4double phi; << 105 G4double y = 0.; << 106 G4double f = 0.; << 107 G4double theta0 = 0.; << 108 G4double phi0 = 0.; << 109 << 110 switch (sourceType) { << 111 case 0: << 112 particleGun->SetParticlePosition(v << 113 particleGun->SetParticleMomentumDi << 114 break; << 115 case 1: << 116 // GS: Generate random position on << 117 // GS: on the sphere << 118 phi = G4UniformRand() * twopi; << 119 do { << 120 y = G4UniformRand() * 1.0; << 121 theta = G4UniformRand() * pi; << 122 f = std::sin(theta); << 123 } while (y > f); << 124 vertex0 = G4ThreeVector(1., 0., 0. << 125 vertex0.setMag(vertexRadius); << 126 vertex0.setTheta(theta); << 127 vertex0.setPhi(phi); << 128 particleGun->SetParticlePosition(v << 129 << 130 momentumDirection0 = G4ThreeVector << 131 << 132 do { << 133 phi = G4UniformRand() * twopi; << 134 do { << 135 y = G4UniformRand() * 1.0; << 136 theta = G4UniformRand() * << 137 f = std::sin(theta); << 138 } while (y > f); << 139 momentumDirection0.setPhi(phi) << 140 momentumDirection0.setTheta(th << 141 } while (vertex0.dot(momentumDirec << 142 << 143 particleGun->SetParticleMomentumDi << 144 << 145 break; << 146 case 2: << 147 // GS: Generate random position on << 148 // GS: on a plane << 149 phi = G4UniformRand() * twopi; << 150 << 151 do { << 152 y = G4UniformRand() * 1.0; << 153 theta = G4UniformRand() * half << 154 f = std::sin(theta) * std::cos << 155 } while (y > f); << 156 << 157 vertex0 = G4ThreeVector(1., 0., 0. << 158 << 159 auto xy = detector->GetWorldSizeXY << 160 auto z = detector->GetWorldSizeZ() << 161 << 162 if (vertexRadius > xy * 0.5) { << 163 G4cout << "vertexRadius too bi << 164 G4cout << "vertexRadius set to << 165 vertexRadius = xy * 0.45; << 166 } << 167 << 168 if (vertexRadius > z * 0.5) { << 169 G4cout << "vertexRadius too hi << 170 G4cout << "vertexRadius set to << 171 vertexRadius = z * 0.45; << 172 } << 173 << 174 vertex0.setMag(vertexRadius); << 175 vertex0.setTheta(theta); << 176 vertex0.setPhi(phi); << 177 << 178 // GS: Get the user defined direct << 179 // GS: Rotate the random position << 180 << 181 momentumDirection0 = particleGun-> << 182 if (momentumDirection0.mag() > 0.0 << 183 theta0 = momentumDirection0.th << 184 phi0 = momentumDirection0.phi( << 185 } << 186 << 187 if (theta0 != 0.) { << 188 G4ThreeVector rotationAxis(1., << 189 rotationAxis.setPhi(phi0 + hal << 190 vertex0.rotate(theta0 + pi, ro << 191 } << 192 particleGun->SetParticlePosition(v << 193 break; << 194 } << 195 << 196 constexpr auto INITIAL_PARTICLE_ENERGY << 197 G4double particleEnergy = INITIAL_PART << 198 << 199 switch (spectrumType) { << 200 case 0: // Uniform energy (1 GeV - 10 << 201 y = G4UniformRand(); << 202 particleEnergy = y * 9.0 * GeV + 1 << 203 G4cout << "Particle energy: " << p << 204 break; << 205 case 1: // Logarithmic energy << 206 y = G4UniformRand(); << 207 particleEnergy = std::pow(10, y) * << 208 G4cout << "Particle energy: " << p << 209 break; << 210 case 2: // Power law (-4) << 211 do { << 212 y = G4UniformRand() * 100000.0 << 213 particleEnergy = G4UniformRand << 214 f = std::pow(particleEnergy * << 215 } while (y > f); << 216 // particleGun->SetParticleEnergy( << 217 break; << 218 case 3: // Monochromatic << 219 particleEnergy = particleGun->GetP << 220 // 100 MeV; << 221 G4cout << "Particle energy: " << p << 222 break; << 223 } << 224 particleGun->SetParticleEnergy(particl << 225 G4cout << "Particle: " << particleGun- << 226 particleGun->GeneratePrimaryVertex(eve << 227 } else { << 228 particleSource->GeneratePrimaryVertex( << 229 } << 230 } << 231 246