ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
XrayFluoMercuryPrimaryGeneratorAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file XrayFluoMercuryPrimaryGeneratorAction.cc
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