Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 // 28 // Author: Alfonso Mantero (Alfonso.Mantero@ge.infn.it) 29 // 30 // History: 31 // ----------- 32 // 02 Sep 2003 Alfonso Mantero created 33 // 34 // ------------------------------------------------------------------- 35 36 #include "XrayFluoMercuryPrimaryGeneratorAction.hh" 37 #include "XrayFluoMercuryDetectorConstruction.hh" 38 #include "XrayFluoMercuryPrimaryGeneratorMessenger.hh" 39 #include "XrayFluoRunAction.hh" 40 #include "XrayFluoAnalysisManager.hh" 41 #include "XrayFluoDataSet.hh" 42 #include "G4PhysicalConstants.hh" 43 #include "G4SystemOfUnits.hh" 44 #include "G4DataVector.hh" 45 #include "G4Event.hh" 46 #include "G4ParticleGun.hh" 47 #include "G4ParticleTable.hh" 48 #include "G4ParticleDefinition.hh" 49 #include "Randomize.hh" 50 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 52 53 XrayFluoMercuryPrimaryGeneratorAction::XrayFluoMercuryPrimaryGeneratorAction(const XrayFluoMercuryDetectorConstruction* XrayFluoDC) 54 :globalFlag(false),spectrum("off") 55 { 56 57 XrayFluoDetector = XrayFluoDC; 58 59 G4int n_particle = 1; 60 particleGun = new G4ParticleGun(n_particle); 61 62 //create a messenger for this class 63 gunMessenger = new XrayFluoMercuryPrimaryGeneratorMessenger(this); 64 runManager = new XrayFluoRunAction(); 65 66 // default particle kinematic 67 68 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable(); 69 G4String particleName; 70 G4ParticleDefinition* particle 71 = particleTable->FindParticle(particleName="gamma"); 72 particleGun->SetParticleDefinition(particle); 73 particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,-1.)); 74 75 76 particleGun->SetParticleEnergy(10.*keV); 77 G4double position = -0.5*(XrayFluoDetector->GetWorldSizeZ()); 78 particleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,position)); 79 80 G4cout << "XrayFluoMercuryPrimaryGeneratorAction created" << G4endl; 81 82 } 83 84 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 86 87 XrayFluoMercuryPrimaryGeneratorAction::~XrayFluoMercuryPrimaryGeneratorAction() 88 { 89 delete particleGun; 90 delete gunMessenger; 91 delete runManager; 92 93 G4cout << "XrayFluoMercuryPrimaryGeneratorAction deleted" << G4endl; 94 95 } 96 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 98 99 void XrayFluoMercuryPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) 100 { 101 //this function is called at the begining of event 102 // 103 104 // Conidering the sunas a Poin-like source. 105 106 G4double z0 = -0.5*(XrayFluoDetector->GetWorldSizeZ()); 107 G4double y0 = 0.*m, x0 = 0.*m; 108 109 110 // Let's try to illuminate only the prtion of Mercury surface that can be seen by the detector. 111 112 G4double spacecraftLatitude = XrayFluoDetector->GetOrbitInclination(); 113 G4double mercuryDia = XrayFluoDetector->GetMercuryDia(); 114 G4double sunDia = XrayFluoDetector->GetSunDia(); 115 G4double opticField = XrayFluoDetector->GetOpticAperture(); 116 117 118 G4double a = 2*std::tan(opticField/2); 119 120 // if (!pointLikeFlag) { 121 122 // let's decide from wich point of the sun surface the particle is coming: 123 124 G4double theta = std::acos(2.*G4UniformRand() - 1.0); 125 G4double phi = 2. * pi * G4UniformRand(); 126 G4double rho = sunDia/2; 127 128 G4double sunPosX = x0 + rho * std::sin(theta) * std::cos(phi); 129 G4double sunPosY = y0 + rho * std::sin(theta) * std::sin(phi); 130 G4double sunPosZ = z0 + rho * std::cos(theta); 131 132 particleGun->SetParticlePosition(G4ThreeVector(sunPosX,sunPosY,sunPosZ)); 133 134 // the angle at the center of Mercury subtending the area seen by the optics: 135 G4double alpha = 2 * a/mercuryDia; 136 137 if(!globalFlag){ 138 theta = alpha * G4UniformRand() + (180.*deg - spacecraftLatitude)-alpha/2.; 139 phi = alpha * G4UniformRand() + 90. * deg - alpha/2.; 140 } 141 142 else if(globalFlag){ 143 theta = pi/2. * rad * G4UniformRand() + 90.*deg ; //was 900., probably an error 144 phi = 2*pi*rad * G4UniformRand() ; 145 } 146 147 rho = mercuryDia/2.; 148 149 G4double mercuryPosX = rho * std::sin(theta) * std::cos(phi); 150 G4double mercuryPosY = rho * std::sin(theta) * std::sin(phi); 151 G4double mercuryPosZ = rho * std::cos(theta); 152 153 particleGun->SetParticleMomentumDirection( 154 G4ThreeVector(mercuryPosX-sunPosX ,mercuryPosY-sunPosY,mercuryPosZ-sunPosZ)); 155 156 // } 157 // if (pointLikeFlag) { 158 159 // // theta is the angle that the mean direction of the incident light (on the desired 160 // // point of the surface of Mercury) makes with the Z-axis 161 // G4double theta = std::asin( mercuryDia/2. * std::sin(spacecraftLatitude) / 162 // std::sqrt(std::pow(z0,2)+std::pow(mercuryDia/2.,2)-2*mercuryDia/2.*z0*std::cos(spacecraftLatitude)) ); 163 164 // // on the y axis, the light emitted from the Sun must be in [theta-phi;theta+phi] 165 // G4double phi = std::asin( mercuryDia/2.*std::sin(spacecraftLatitude) + a*std::cos(spacecraftLatitude) / 166 // std::sqrt( std::pow(mercuryDia/2.*std::sin(spacecraftLatitude) + a*std::cos(spacecraftLatitude) , 2) + 167 // std::pow(z0 - mercuryDia/2.*std::cos(spacecraftLatitude) - a*std::sin(spacecraftLatitude) , 2)) ) 168 // - theta; 169 170 // // on the x axis, the light emitted from the Sun must be in [-zeta;zeta] 171 // G4double zeta = std::atan( a/std::sqrt(std::pow(z0,2)+std::pow(mercuryDia,2)-2*mercuryDia*z0*std::cos(spacecraftLatitude)) ); 172 173 174 175 // //alpha in [-zeta;zeta] 176 // G4double alpha = (2*zeta)*G4UniformRand() - zeta; 177 // //beta in [theta-phi;theta+phi] 178 // G4double beta = (G4UniformRand()*2*phi) - phi + theta; 179 180 // G4double dirY = std::sin(beta); 181 // G4double dirX = std::sin(alpha); 182 183 // particleGun->SetParticleMomentumDirection(G4ThreeVector(dirX.,dirY,1.)); 184 185 // particleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0)); 186 187 // } 188 189 190 191 //shoot particles according to a certain spectrum 192 if (spectrum =="on") 193 { 194 G4String particle = particleGun->GetParticleDefinition() 195 ->GetParticleName(); 196 if(particle == "proton"|| particle == "alpha") 197 { 198 G4DataVector* energies = runManager->GetEnergies(); 199 G4DataVector* data = runManager->GetData(); 200 201 G4double sum = runManager->GetDataSum(); 202 G4double partSum = 0; 203 G4int j = 0; 204 G4double random= sum*G4UniformRand(); 205 while (partSum<random) 206 { 207 partSum += (*data)[j]; 208 j++; 209 } 210 211 particleGun->SetParticleEnergy((*energies)[j]); 212 213 } 214 else if (particle == "gamma") 215 { 216 const XrayFluoDataSet* dataSet = runManager->GetGammaSet(); 217 218 G4int i = 0; 219 G4int id = 0; 220 G4double minEnergy = 0. * keV; 221 G4double particleEnergy= 0.; 222 G4double maxEnergy = 10. * keV; 223 G4double energyRange = maxEnergy - minEnergy; 224 225 while ( i == 0) 226 { 227 G4double random = G4UniformRand(); 228 229 G4double randomNum = G4UniformRand(); //*5.0E6; 230 231 particleEnergy = (random*energyRange) + minEnergy; 232 233 if ((dataSet->FindValue(particleEnergy,id)) > randomNum) 234 { 235 i = 1; 236 237 } 238 } 239 particleGun->SetParticleEnergy(particleEnergy); 240 } 241 } 242 243 244 #ifdef G4ANALYSIS_USE 245 246 G4double partEnergy = particleGun->GetParticleEnergy(); 247 XrayFluoAnalysisManager* analysis = XrayFluoAnalysisManager::getInstance(); 248 analysis->analysePrimaryGenerator(partEnergy/keV); 249 250 #endif 251 252 particleGun->GeneratePrimaryVertex(anEvent); 253 } 254 255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 256 257 258 259 260 261 262 263 264 265