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