Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ML2PrimaryGenerationAction.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 // The code was written by :
27 // ^Claudio Andenna claudio.andenna@ispesl.it, claudio.andenna@iss.infn.it
28 // *Barbara Caccia barbara.caccia@iss.it
29 // with the support of Pablo Cirrone (LNS, INFN Catania Italy)
30 // with the contribute of Alessandro Occhigrossi*
31 //
32 // ^INAIL DIPIA - ex ISPESL and INFN Roma, gruppo collegato Sanità, Italy
33 // *Istituto Superiore di Sanità and INFN Roma, gruppo collegato Sanità, Italy
34 // Viale Regina Elena 299, 00161 Roma (Italy)
35 // tel (39) 06 49902246
36 // fax (39) 06 49387075
37 //
38 // more information:
39 // http://g4advancedexamples.lngs.infn.it/Examples/medical-linac
40 //
41 //*******************************************************//
42 
44 
45 using namespace CLHEP;
46 
48 :particleGun(0),gamma(0),electron(0),positron(0),primaryParticleData(0),particles(0)
49 {
50 }
51 CML2PrimaryGenerationAction* CML2PrimaryGenerationAction::instance = 0;
52 
54 {
55  if (instance == 0)
56  {
57  instance = new CML2PrimaryGenerationAction();
58  }
59  return instance;
60 }
62 {
63  rm=new G4RotationMatrix();
64  PrimaryGenerationActionMessenger=new CML2PrimaryGenerationActionMessenger(this);
65  particle=new Sparticle;
66  nParticle=nPhSpParticles=nRandomParticles=0;
67 
69 
70  gamma=particleTable->FindParticle("gamma");
71  electron=particleTable->FindParticle("e-");
72  positron=particleTable->FindParticle("e+");
73  particleGun=new G4ParticleGun();
74 
75  primaryParticleData=pData;
76  primaryParticleData->nPrimaryParticle=0;
77  primaryParticleData->partPDGE=0;
78 }
79 
81 {
82  accTargetZPosition=aTZ;
83  switch (idParticleSource)
84  {
85  case id_randomTarget:
86  setGunRandom();
87  break;
88  case id_phaseSpace:
89  setGunCalculatedPhaseSpace();
90  break;
91  }
92 }
93 
94 void CML2PrimaryGenerationAction::setGunRandom()
95 {
96  particleGun->SetParticleDefinition(electron);
97  particleGun->SetNumberOfParticles(1);
98  idCurrentParticleSource=idParticleSource;
99 }
100 
101 void CML2PrimaryGenerationAction::setGunCalculatedPhaseSpace()
102 {
103  particles=new Sparticle[nMaxParticlesInRamPhaseSpace];
104  particleGun->SetNumberOfParticles(1);
105  idCurrentParticleSource=idParticleSource;
106 }
107 
109 {
110  delete particleGun;
111  delete [] particles;
112  delete particles;
113 }
115 {
116  static int currentRecycle=nRecycling;
117  static G4ThreeVector pos0, dir0;
118  if (currentRecycle==nRecycling)
119  {
120  currentRecycle=0;
121  switch (idCurrentParticleSource)
122  {
123  case id_randomTarget:
124  GenerateFromRandom();
125  break;
126  case id_phaseSpace:
127  GenerateFromCalculatedPhaseSpace();
128  break;
129  }
130  pos0=pos;
131  dir0=dir;
132  }
133  currentRecycle++;
134  pos=pos0;
135  dir=dir0;
136  applySourceRotation(); // to follow the accelerator rotation
137 
138  primaryParticleData->partPDGE=particleGun->GetParticleDefinition()->GetPDGEncoding();
139  primaryParticleData->nPrimaryParticle++;
140 
141  particleGun->SetParticleEnergy(ek*MeV);
142  particleGun->SetParticlePosition(pos*mm);
144  particleGun->GeneratePrimaryVertex(anEvent);
145 }
146 void CML2PrimaryGenerationAction::GenerateFromRandom()
147 {
148  sinTheta=RandGauss::shoot(0., 0.003);
149  cosTheta=std::sqrt(1 - sinTheta*sinTheta);
150  phi=twopi*G4UniformRand();
151  dir.set(sinTheta*std::cos(phi), sinTheta*std::sin(phi), cosTheta);
152 
153  ro=G4UniformRand()*GunRadious;
154  alfa=G4UniformRand()*twopi;
155 
156  pos.setX(ro*std::sin(alfa));
157  pos.setY(ro*std::cos(alfa));
158  pos.setZ(-(accTargetZPosition +5.)*mm); // the primary electrons are generated 5 mm before the target
159  ek=RandGauss::shoot(GunMeanEnegy, GunStdEnegy);
160  nRandomParticles++;
161 }
162 void CML2PrimaryGenerationAction::GenerateFromCalculatedPhaseSpace()
163 {
164  static bool bFirstTime=true;
165  if (bFirstTime)
166  {bFirstTime=false;fillParticlesContainer();}
167  if (nParticle==nMaxParticlesInRamPhaseSpace) // once all the particles stored in RAM hae been processed a new set is loaded
168  {
169  fillParticlesContainer();
170  nParticle=0;
171  }
172 
173  pos=particles[nParticle].pos;
174  dir=particles[nParticle].dir;
175  ek=particles[nParticle].kinEnergy;
176  switch (particles[nParticle].partPDGE)
177  {
178  case -11:
179  particleGun->SetParticleDefinition(positron);
180  break;
181  case 11:
182  particleGun->SetParticleDefinition(electron);
183  break;
184  case 22:
185  particleGun->SetParticleDefinition(gamma);
186  break;
187  }
188  nPhSpParticles++;
189  nParticle++;
190 }
191 void CML2PrimaryGenerationAction::applySourceRotation()
192 {
193  pos=*rm*pos;
194  dir=*rm*dir;
195 }
196 void CML2PrimaryGenerationAction::fillParticlesContainer()
197 {
198  static int currentFilePosition=0;
199  static int currentFileSize=0;
200  int startDataFilePosition;
201  std::ifstream in;
202  in.open(calculatedPhaseSpaceFileIN, std::ios::in);
203  if (in)
204  {
205  std::cout <<"ERROR phase space file: "<< calculatedPhaseSpaceFileIN << " NOT found. Run abort "<< G4endl;
207  }
208 
209  static bool bFirstTime=true;
210  if (bFirstTime)
211  {
212  in.seekg(-1,std::ios::end);
213  currentFileSize=in.tellg();
214  in.seekg(0,std::ios::beg);
215  bFirstTime=false;
216  }
217 
218  char a[1000];
219  in.getline(a,1000,'\n');
220  in.getline(a,1000,'\n');
221  startDataFilePosition=in.tellg();
222  if (currentFilePosition>0)
223  {in.seekg(currentFilePosition, std::ios::beg);}
224  int i;
225  G4double x,y,z;
226  G4int d;
227 
228  static bool checkFileRewind=false;
229  static bool bRewindTheFile=false;
230  static int nPhSpFileRewind=0;
231 
232  for (i=0;i<nMaxParticlesInRamPhaseSpace;i++)
233  {
234  if (bRewindTheFile) // to read the phase space file again to fill the container
235  {
236  in.close();
237  in.open(calculatedPhaseSpaceFileIN, std::ios::in);
238  in.seekg(startDataFilePosition, std::ios::beg);
239  checkFileRewind=true;
240  bRewindTheFile=false;
241  std::cout<<"\n################\nI have reached the end of the phase space file "<<++nPhSpFileRewind <<" times, I rewind the file\n" << G4endl;
242  std::cout <<"loaded " <<i <<"/"<< nMaxParticlesInRamPhaseSpace<<" particles" << G4endl;
243  }
244  in >> d;
245  in >> x; in >>y; in >> z;
246 /* std::cout <<"x:" <<x << G4endl;
247  std::cout <<"y:" <<y << G4endl;
248  std::cout <<"z:" <<z << G4endl;*/
249  particles[i].pos.set(x,y,z-accTargetZPosition);
250  in >> x; in >>y; in >> z;
251  particles[i].dir.set(x,y,z);
252  in >> x;
253  particles[i].kinEnergy=x;
254  in >> d;
255  particles[i].partPDGE=d;
256  in >> d; in >> d;
257  if (in.eof()) {bRewindTheFile=true;}
258  if (checkFileRewind) {checkFileRewind=false;}
259  }
260  std::cout <<"loaded " <<i <<"/"<< nMaxParticlesInRamPhaseSpace<<" particles" << G4endl;
261  currentFilePosition=in.tellg(); // to remind the actual position in the phase space file
262  if (currentFilePosition>=currentFileSize) // to read the phase space file again
263  {currentFilePosition=startDataFilePosition;}
264  in.close();
265 }
266 
void set(double x, double y, double z)
virtual void AbortRun(G4bool softAbort=false)
ThreeVector shoot(const G4int Ap, const G4int Af)
idParticleSource
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static constexpr double mm
Definition: G4SIunits.hh:115
CLHEP::HepRotation G4RotationMatrix
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4int partPDGE
tuple x
Definition: test.py:50
void SetParticleMomentumDirection(G4ParticleMomentum aMomentumDirection)
G4ThreeVector dir
int G4int
Definition: G4Types.hh:78
void setY(double)
virtual void GeneratePrimaryVertex(G4Event *evt)
void setZ(double)
void setX(double)
static constexpr double twopi
Definition: G4SIunits.hh:76
void design(G4double accTargetZPosition)
void inizialize(SPrimaryParticle *primaryParticleData)
void SetParticlePosition(G4ThreeVector aPosition)
#define G4UniformRand()
Definition: Randomize.hh:97
void SetNumberOfParticles(G4int i)
void SetParticleEnergy(G4double aKineticEnergy)
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:79
G4ThreeVector pos
static G4ParticleTable * GetParticleTable()
G4double kinEnergy
G4ParticleDefinition * GetParticleDefinition() const
tuple z
Definition: test.py:28
static CML2PrimaryGenerationAction * GetInstance(void)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76
void SetParticleDefinition(G4ParticleDefinition *aParticleDefinition)