ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
XrayFluoPrimaryGeneratorAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file XrayFluoPrimaryGeneratorAction.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: Elena Guardincerri (Elena.Guardincerri@ge.infn.it)
29 //
30 // History:
31 // -----------
32 // 28 Nov 2001 Elena Guardincerri Created
33 //
34 // -------------------------------------------------------------------
35 
36 #include "XrayFluoPrimaryGeneratorAction.hh"
37 #include "XrayFluoDetectorConstruction.hh"
38 #include "XrayFluoPrimaryGeneratorMessenger.hh"
39 #include "XrayFluoRunAction.hh"
40 #include "G4Event.hh"
41 #include "G4Gamma.hh"
42 #include "G4ParticleGun.hh"
43 #include "G4ParticleTable.hh"
44 #include "G4ParticleDefinition.hh"
45 #include "G4MTRunManager.hh"
46 #include "Randomize.hh"
47 #include "XrayFluoAnalysisManager.hh"
48 #include "XrayFluoDataSet.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53 
54 XrayFluoPrimaryGeneratorAction::XrayFluoPrimaryGeneratorAction(const
55  XrayFluoDetectorConstruction* XrayFluoDC)
56  :rndmFlag("off"),beam("off"),spectrum("off"),isoVert("off"),phaseSpaceGunFlag(false),
57  rayleighFlag(true), detectorPosition(0)
58 {
59  runAction = 0;
60  XrayFluoDetector = XrayFluoDC;
61 
62  G4int n_particle = 1;
63  particleGun = new G4ParticleGun(n_particle);
64 
65  //create a messenger for this class
66  gunMessenger = new XrayFluoPrimaryGeneratorMessenger(this);
67 
68  // default particle kinematic
69  G4ParticleDefinition* particle
70  = G4Gamma::Definition();
71  particleGun->SetParticleDefinition(particle);
72  particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
73  particleGun->SetParticleEnergy(10. * keV);
74 
75  G4double position = -0.5*(XrayFluoDetector->GetWorldSizeZ());
76  particleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,position));
77 
78  G4cout << "XrayFluoPrimaryGeneratorAction created" << G4endl;
79 
80 }
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 
84 void XrayFluoPrimaryGeneratorAction::ActivatePhaseSpace(G4String fileName) {
85 
86  // load phase-space
87  phaseSpaceGunFlag = true;
88 
89  // reads the data stored on disk form previous runs
90  // and get these data to data members
91 
92  XrayFluoAnalysisManager* analysis = XrayFluoAnalysisManager::getInstance();
93  analysis->LoadGunData(fileName, rayleighFlag);
94  detectorPosition = XrayFluoDetector->GetDetectorPosition();
95  detectorPosition.setR(detectorPosition.r()-(5.*cm)); // 5 cm before the detector, so in front of it.
96 
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 
101 void XrayFluoPrimaryGeneratorAction::SetRayleighFlag (G4bool value)
102 {
103  rayleighFlag = value;
104 }
105 
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
108 
109 XrayFluoPrimaryGeneratorAction::~XrayFluoPrimaryGeneratorAction()
110 {
111  delete particleGun;
112  delete gunMessenger;
113 }
114 
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116 
117 void XrayFluoPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
118 {
119  //retrieve runAction, if not done
120  if (!runAction)
121  {
122  //Sequential runaction
123  if (G4RunManager::GetRunManager()->GetRunManagerType() ==
124  G4RunManager::sequentialRM)
125  runAction = static_cast<const XrayFluoRunAction*>
126  (G4RunManager::GetRunManager()->GetUserRunAction());
127  else //MT master runaction
128  runAction = static_cast<const XrayFluoRunAction*>
129  (G4MTRunManager::GetMasterRunManager()->GetUserRunAction());
130  if (!runAction)
131  G4cout << "Something wrong here!" << G4endl;
132  }
133 
134  //this function is called at the begining of event
135  //
136  G4double z0 = -0.5*(XrayFluoDetector->GetWorldSizeZ());
137  G4double y0 = 0.*cm, x0 = 0.*cm;
138  if (rndmFlag == "on")
139  {
140  y0 = (XrayFluoDetector->GetDia3SizeXY())/std::sqrt(2.)*(G4UniformRand()-0.5); // it was GetSampleSizeXY(),
141  x0 = (XrayFluoDetector->GetDia3SizeXY())/std::sqrt(2.)*(G4UniformRand()-0.5); // not divided by std::sqrt(2.)
142  }
143  particleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
144 
145  //randomize starting point
146  if (beam == "on")
147  {
148  G4double radius = 0.5 * mm;
149  G4double rho = radius*std::sqrt(G4UniformRand());
150  G4double theta = 2*pi*G4UniformRand()*rad;
151  G4double position = -0.5*(XrayFluoDetector->GetWorldSizeZ());
152 
153  G4double y = rho * std::sin(theta);
154  G4double x = rho * std::cos(theta);
155 
156  particleGun->SetParticlePosition(G4ThreeVector(x,y,position));
157  }
158  //shoot particles according to a certain spectrum
159  if (spectrum =="on")
160  {
161  G4String particle = particleGun->GetParticleDefinition()
162  ->GetParticleName();
163  if(particle == "proton"|| particle == "alpha")
164  {
165  G4DataVector* energies = runAction->GetEnergies();
166  G4DataVector* data = runAction->GetData();
167 
168  G4double sum = runAction->GetDataSum();
169  G4double partSum = 0;
170  G4int j = 0;
171  G4double random= sum*G4UniformRand();
172  while (partSum<random)
173  {
174  partSum += (*data)[j];
175  j++;
176  }
177 
178  particleGun->SetParticleEnergy((*energies)[j]);
179 
180  }
181  else if (particle == "gamma")
182  {
183  const XrayFluoDataSet* dataSet = runAction->GetGammaSet();
184 
185  G4int i = 0;
186  G4int id = 0;
187  G4double minEnergy = 0. * keV;
188  G4double particleEnergy= 0.;
189  G4double maxEnergy = 10. * keV;
190  G4double energyRange = maxEnergy - minEnergy;
191 
192  while ( i == 0)
193  {
194  G4double random = G4UniformRand();
195 
196  G4double randomNum = G4UniformRand(); //*5.0E6;
197 
198  particleEnergy = (random*energyRange) + minEnergy;
199 
200  if ((dataSet->FindValue(particleEnergy,id)) > randomNum)
201  {
202  i = 1;
203 
204  }
205  }
206  particleGun->SetParticleEnergy(particleEnergy);
207  }
208  }
209 
210  // Randomize starting point and direction
211 
212  if (isoVert == "on")
213  {
214  G4double rho = 1. *m;
215  //theta in [0;pi/2]
216  G4double theta = (pi/2)*G4UniformRand();
217  //phi in [-pi;pi]
218  G4double phi = (G4UniformRand()*2*pi)- pi;
219  G4double x = rho*std::sin(theta)*std::sin(phi);
220  G4double y = rho*std::sin(theta)*std::cos(phi);
221  G4double z = -(rho*std::cos(theta));
222  particleGun->SetParticlePosition(G4ThreeVector(x,y,z));
223 
224  G4double Xdim = XrayFluoDetector->GetSampleSizeXY();
225  G4double Ydim = XrayFluoDetector->GetSampleSizeXY();
226 
227  G4double Dx = Xdim*(G4UniformRand()-0.5);
228 
229  G4double Dy = Ydim*(G4UniformRand()-0.5);
230 
231  particleGun->SetParticleMomentumDirection(G4ThreeVector(-x+Dx,-y+Dy,-z));
232 
233  }
234 
235  // using prevoiously genereated emissions from sample.....
236 
237  if (phaseSpaceGunFlag){
238 
239  particleGun->SetParticlePosition(detectorPosition);
240  particleGun->SetParticleMomentumDirection(detectorPosition);
241 
242  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
243 
244  const std::pair<G4double,G4String> kine =
245  XrayFluoAnalysisManager::getInstance()->GetEmittedParticleEnergyAndType();
246 
247  G4double energy = kine.first;
248  G4ParticleDefinition* particle = particleTable->FindParticle(kine.second);
249 
250  particleGun->SetParticleEnergy(energy);
251  particleGun->SetParticleDefinition(particle);
252 
253 
254  }
255 
256  G4double partEnergy = particleGun->GetParticleEnergy();
257  XrayFluoAnalysisManager* analysis = XrayFluoAnalysisManager::getInstance();
258  analysis->analysePrimaryGenerator(partEnergy/keV);
259 
260 
261  particleGun->GeneratePrimaryVertex(anEvent);
262 }
263 
264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
265 
266 
267 
268 
269 
270 
271 
272 
273