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 // 27 // 28 // Author: Alfonso Mantero (Alfonso.Mantero@ge 28 // Author: Alfonso Mantero (Alfonso.Mantero@ge.infn.it) 29 // 29 // 30 // History: 30 // History: 31 // ----------- 31 // ----------- 32 // 02 Sep 2003 Alfonso Mantero created 32 // 02 Sep 2003 Alfonso Mantero created 33 // 33 // 34 // ------------------------------------------- 34 // ------------------------------------------------------------------- 35 35 36 #include "XrayFluoMercuryPrimaryGeneratorActio 36 #include "XrayFluoMercuryPrimaryGeneratorAction.hh" 37 #include "XrayFluoMercuryDetectorConstruction. 37 #include "XrayFluoMercuryDetectorConstruction.hh" 38 #include "XrayFluoMercuryPrimaryGeneratorMesse 38 #include "XrayFluoMercuryPrimaryGeneratorMessenger.hh" 39 #include "XrayFluoRunAction.hh" 39 #include "XrayFluoRunAction.hh" 40 #include "XrayFluoAnalysisManager.hh" 40 #include "XrayFluoAnalysisManager.hh" 41 #include "XrayFluoDataSet.hh" 41 #include "XrayFluoDataSet.hh" 42 #include "G4PhysicalConstants.hh" 42 #include "G4PhysicalConstants.hh" 43 #include "G4SystemOfUnits.hh" 43 #include "G4SystemOfUnits.hh" 44 #include "G4DataVector.hh" 44 #include "G4DataVector.hh" 45 #include "G4Event.hh" 45 #include "G4Event.hh" 46 #include "G4ParticleGun.hh" 46 #include "G4ParticleGun.hh" 47 #include "G4ParticleTable.hh" 47 #include "G4ParticleTable.hh" 48 #include "G4ParticleDefinition.hh" 48 #include "G4ParticleDefinition.hh" 49 #include "Randomize.hh" 49 #include "Randomize.hh" 50 50 51 //....oooOO0OOooo........oooOO0OOooo........oo 51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 52 52 53 XrayFluoMercuryPrimaryGeneratorAction::XrayFlu 53 XrayFluoMercuryPrimaryGeneratorAction::XrayFluoMercuryPrimaryGeneratorAction(const XrayFluoMercuryDetectorConstruction* XrayFluoDC) 54 :globalFlag(false),spectrum("off") 54 :globalFlag(false),spectrum("off") 55 { 55 { 56 56 57 XrayFluoDetector = XrayFluoDC; 57 XrayFluoDetector = XrayFluoDC; 58 58 59 G4int n_particle = 1; 59 G4int n_particle = 1; 60 particleGun = new G4ParticleGun(n_particle) 60 particleGun = new G4ParticleGun(n_particle); 61 61 62 //create a messenger for this class 62 //create a messenger for this class 63 gunMessenger = new XrayFluoMercuryPrimaryGen 63 gunMessenger = new XrayFluoMercuryPrimaryGeneratorMessenger(this); 64 runManager = new XrayFluoRunAction(); 64 runManager = new XrayFluoRunAction(); 65 65 66 // default particle kinematic 66 // default particle kinematic 67 67 68 G4ParticleTable* particleTable = G4ParticleT 68 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable(); 69 G4String particleName; 69 G4String particleName; 70 G4ParticleDefinition* particle 70 G4ParticleDefinition* particle 71 = particleTable->FindParticle(particleName 71 = particleTable->FindParticle(particleName="gamma"); 72 particleGun->SetParticleDefinition(particle) 72 particleGun->SetParticleDefinition(particle); 73 particleGun->SetParticleMomentumDirection(G4 73 particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,-1.)); 74 74 75 75 76 particleGun->SetParticleEnergy(10.*keV); 76 particleGun->SetParticleEnergy(10.*keV); 77 G4double position = -0.5*(XrayFluoDetector-> 77 G4double position = -0.5*(XrayFluoDetector->GetWorldSizeZ()); 78 particleGun->SetParticlePosition(G4ThreeVect 78 particleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,position)); 79 79 80 G4cout << "XrayFluoMercuryPrimaryGeneratorAc 80 G4cout << "XrayFluoMercuryPrimaryGeneratorAction created" << G4endl; 81 81 82 } 82 } 83 83 84 84 85 //....oooOO0OOooo........oooOO0OOooo........oo 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 86 86 87 XrayFluoMercuryPrimaryGeneratorAction::~XrayFl 87 XrayFluoMercuryPrimaryGeneratorAction::~XrayFluoMercuryPrimaryGeneratorAction() 88 { 88 { 89 delete particleGun; 89 delete particleGun; 90 delete gunMessenger; 90 delete gunMessenger; 91 delete runManager; 91 delete runManager; 92 92 93 G4cout << "XrayFluoMercuryPrimaryGeneratorAc 93 G4cout << "XrayFluoMercuryPrimaryGeneratorAction deleted" << G4endl; 94 94 95 } 95 } 96 96 97 //....oooOO0OOooo........oooOO0OOooo........oo 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 98 98 99 void XrayFluoMercuryPrimaryGeneratorAction::Ge 99 void XrayFluoMercuryPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) 100 { 100 { 101 //this function is called at the begining of 101 //this function is called at the begining of event 102 // 102 // 103 103 104 // Conidering the sunas a Poin-like source. 104 // Conidering the sunas a Poin-like source. 105 105 106 G4double z0 = -0.5*(XrayFluoDetector->GetWor 106 G4double z0 = -0.5*(XrayFluoDetector->GetWorldSizeZ()); 107 G4double y0 = 0.*m, x0 = 0.*m; 107 G4double y0 = 0.*m, x0 = 0.*m; 108 108 109 109 110 // Let's try to illuminate only the prtion o 110 // Let's try to illuminate only the prtion of Mercury surface that can be seen by the detector. 111 111 112 G4double spacecraftLatitude = XrayFluoDetect 112 G4double spacecraftLatitude = XrayFluoDetector->GetOrbitInclination(); 113 G4double mercuryDia = XrayFluoDetector->GetM 113 G4double mercuryDia = XrayFluoDetector->GetMercuryDia(); 114 G4double sunDia = XrayFluoDetector->GetSunDi 114 G4double sunDia = XrayFluoDetector->GetSunDia(); 115 G4double opticField = XrayFluoDetector->GetO 115 G4double opticField = XrayFluoDetector->GetOpticAperture(); 116 116 117 117 118 G4double a = 2*std::tan(opticField/2); 118 G4double a = 2*std::tan(opticField/2); 119 119 120 // if (!pointLikeFlag) { 120 // if (!pointLikeFlag) { 121 121 122 // let's decide from wich point of the sun 122 // let's decide from wich point of the sun surface the particle is coming: 123 123 124 G4double theta = std::acos(2.*G4UniformRand( 124 G4double theta = std::acos(2.*G4UniformRand() - 1.0); 125 G4double phi = 2. * pi * G4UniformRand(); 125 G4double phi = 2. * pi * G4UniformRand(); 126 G4double rho = sunDia/2; 126 G4double rho = sunDia/2; 127 127 128 G4double sunPosX = x0 + rho * std::sin(the 128 G4double sunPosX = x0 + rho * std::sin(theta) * std::cos(phi); 129 G4double sunPosY = y0 + rho * std::sin(the 129 G4double sunPosY = y0 + rho * std::sin(theta) * std::sin(phi); 130 G4double sunPosZ = z0 + rho * std::cos(the 130 G4double sunPosZ = z0 + rho * std::cos(theta); 131 131 132 particleGun->SetParticlePosition(G4ThreeVe 132 particleGun->SetParticlePosition(G4ThreeVector(sunPosX,sunPosY,sunPosZ)); 133 133 134 // the angle at the center of Mercury subt 134 // the angle at the center of Mercury subtending the area seen by the optics: 135 G4double alpha = 2 * a/mercuryDia; 135 G4double alpha = 2 * a/mercuryDia; 136 136 137 if(!globalFlag){ 137 if(!globalFlag){ 138 theta = alpha * G4UniformRand() + (180.* 138 theta = alpha * G4UniformRand() + (180.*deg - spacecraftLatitude)-alpha/2.; 139 phi = alpha * G4UniformRand() + 90. * de 139 phi = alpha * G4UniformRand() + 90. * deg - alpha/2.; 140 } 140 } 141 141 142 else if(globalFlag){ 142 else if(globalFlag){ 143 theta = pi/2. * rad * G4UniformRand() + 143 theta = pi/2. * rad * G4UniformRand() + 90.*deg ; //was 900., probably an error 144 phi = 2*pi*rad * G4UniformRand() ; 144 phi = 2*pi*rad * G4UniformRand() ; 145 } 145 } 146 146 147 rho = mercuryDia/2.; 147 rho = mercuryDia/2.; 148 148 149 G4double mercuryPosX = rho * std::sin(thet 149 G4double mercuryPosX = rho * std::sin(theta) * std::cos(phi); 150 G4double mercuryPosY = rho * std::sin(thet 150 G4double mercuryPosY = rho * std::sin(theta) * std::sin(phi); 151 G4double mercuryPosZ = rho * std::cos(thet 151 G4double mercuryPosZ = rho * std::cos(theta); 152 152 153 particleGun->SetParticleMomentumDirection( 153 particleGun->SetParticleMomentumDirection( 154 G4ThreeVector(mercuryPosX-sunPosX ,m 154 G4ThreeVector(mercuryPosX-sunPosX ,mercuryPosY-sunPosY,mercuryPosZ-sunPosZ)); 155 155 156 // } 156 // } 157 // if (pointLikeFlag) { 157 // if (pointLikeFlag) { 158 158 159 // // theta is the angle that the mean direc 159 // // theta is the angle that the mean direction of the incident light (on the desired 160 // // point of the surface of Mercury) makes 160 // // point of the surface of Mercury) makes with the Z-axis 161 // G4double theta = std::asin( mercuryDia/2. 161 // G4double theta = std::asin( mercuryDia/2. * std::sin(spacecraftLatitude) / 162 // std::sqrt(std::pow(z0,2)+std::pow(mer 162 // std::sqrt(std::pow(z0,2)+std::pow(mercuryDia/2.,2)-2*mercuryDia/2.*z0*std::cos(spacecraftLatitude)) ); 163 163 164 // // on the y axis, the light emitted from 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.*s 165 // G4double phi = std::asin( mercuryDia/2.*std::sin(spacecraftLatitude) + a*std::cos(spacecraftLatitude) / 166 // std::sqrt( std::pow(mercuryDia/2. 166 // std::sqrt( std::pow(mercuryDia/2.*std::sin(spacecraftLatitude) + a*std::cos(spacecraftLatitude) , 2) + 167 // std::pow(z0 - mercuryDia/2.*std:: 167 // std::pow(z0 - mercuryDia/2.*std::cos(spacecraftLatitude) - a*std::sin(spacecraftLatitude) , 2)) ) 168 // - theta; 168 // - theta; 169 169 170 // // on the x axis, the light emitted from 170 // // on the x axis, the light emitted from the Sun must be in [-zeta;zeta] 171 // G4double zeta = std::atan( a/std::sqrt(st 171 // G4double zeta = std::atan( a/std::sqrt(std::pow(z0,2)+std::pow(mercuryDia,2)-2*mercuryDia*z0*std::cos(spacecraftLatitude)) ); 172 172 173 173 174 174 175 // //alpha in [-zeta;zeta] 175 // //alpha in [-zeta;zeta] 176 // G4double alpha = (2*zeta)*G4UniformRand() 176 // G4double alpha = (2*zeta)*G4UniformRand() - zeta; 177 // //beta in [theta-phi;theta+phi] 177 // //beta in [theta-phi;theta+phi] 178 // G4double beta = (G4UniformRand()*2*phi) - 178 // G4double beta = (G4UniformRand()*2*phi) - phi + theta; 179 179 180 // G4double dirY = std::sin(beta); 180 // G4double dirY = std::sin(beta); 181 // G4double dirX = std::sin(alpha); 181 // G4double dirX = std::sin(alpha); 182 182 183 // particleGun->SetParticleMomentumDirection 183 // particleGun->SetParticleMomentumDirection(G4ThreeVector(dirX.,dirY,1.)); 184 184 185 // particleGun->SetParticlePosition(G4ThreeV 185 // particleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0)); 186 186 187 // } 187 // } 188 188 189 189 190 190 191 //shoot particles according to a certain spe 191 //shoot particles according to a certain spectrum 192 if (spectrum =="on") 192 if (spectrum =="on") 193 { 193 { 194 G4String particle = particleGun->GetPar 194 G4String particle = particleGun->GetParticleDefinition() 195 ->GetParticleName(); 195 ->GetParticleName(); 196 if(particle == "proton"|| particle == "a 196 if(particle == "proton"|| particle == "alpha") 197 { 197 { 198 G4DataVector* energies = runManager->GetE 198 G4DataVector* energies = runManager->GetEnergies(); 199 G4DataVector* data = runManager->GetData( 199 G4DataVector* data = runManager->GetData(); 200 200 201 G4double sum = runManager->GetDataSum(); 201 G4double sum = runManager->GetDataSum(); 202 G4double partSum = 0; 202 G4double partSum = 0; 203 G4int j = 0; 203 G4int j = 0; 204 G4double random= sum*G4UniformRand(); 204 G4double random= sum*G4UniformRand(); 205 while (partSum<random) 205 while (partSum<random) 206 { 206 { 207 partSum += (*data)[j]; 207 partSum += (*data)[j]; 208 j++; 208 j++; 209 } 209 } 210 210 211 particleGun->SetParticleEnergy((*energies) 211 particleGun->SetParticleEnergy((*energies)[j]); 212 212 213 } 213 } 214 else if (particle == "gamma") 214 else if (particle == "gamma") 215 { 215 { 216 const XrayFluoDataSet* dataSet = runManage 216 const XrayFluoDataSet* dataSet = runManager->GetGammaSet(); 217 217 218 G4int i = 0; 218 G4int i = 0; 219 G4int id = 0; 219 G4int id = 0; 220 G4double minEnergy = 0. * keV; 220 G4double minEnergy = 0. * keV; 221 G4double particleEnergy= 0.; 221 G4double particleEnergy= 0.; 222 G4double maxEnergy = 10. * keV; 222 G4double maxEnergy = 10. * keV; 223 G4double energyRange = maxEnergy - minEner 223 G4double energyRange = maxEnergy - minEnergy; 224 224 225 while ( i == 0) 225 while ( i == 0) 226 { 226 { 227 G4double random = G4UniformRand(); 227 G4double random = G4UniformRand(); 228 228 229 G4double randomNum = G4UniformRand(); 229 G4double randomNum = G4UniformRand(); //*5.0E6; 230 230 231 particleEnergy = (random*energyRange) 231 particleEnergy = (random*energyRange) + minEnergy; 232 232 233 if ((dataSet->FindValue(particleEnergy 233 if ((dataSet->FindValue(particleEnergy,id)) > randomNum) 234 { 234 { 235 i = 1; 235 i = 1; 236 236 237 } 237 } 238 } 238 } 239 particleGun->SetParticleEnergy(particleEn 239 particleGun->SetParticleEnergy(particleEnergy); 240 } 240 } 241 } 241 } 242 242 243 243 244 #ifdef G4ANALYSIS_USE 244 #ifdef G4ANALYSIS_USE 245 245 246 G4double partEnergy = particleGun->GetPartic 246 G4double partEnergy = particleGun->GetParticleEnergy(); 247 XrayFluoAnalysisManager* analysis = XrayFlu 247 XrayFluoAnalysisManager* analysis = XrayFluoAnalysisManager::getInstance(); 248 analysis->analysePrimaryGenerator(partEnergy 248 analysis->analysePrimaryGenerator(partEnergy/keV); 249 249 250 #endif 250 #endif 251 251 252 particleGun->GeneratePrimaryVertex(anEvent); 252 particleGun->GeneratePrimaryVertex(anEvent); 253 } 253 } 254 254 255 //....oooOO0OOooo........oooOO0OOooo........oo 255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 256 256 257 257 258 258 259 259 260 260 261 261 262 262 263 263 264 264 265 265