Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
examples/advanced/nanobeam/src/PrimaryGeneratorAction.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 // Please cite the following paper if you use this software
27 // Nucl.Instrum.Meth.B260:20-27, 2007
28 
29 #include "PrimaryGeneratorAction.hh"
30 #include "PrimaryGeneratorMessenger.hh"
31 #include "G4SystemOfUnits.hh"
32 
33 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
34 
35 PrimaryGeneratorAction::PrimaryGeneratorAction(DetectorConstruction* DC)
36 :fDetector(DC)
37 {
38  fAngleMax = 0.09;
39 
40  // Default
41  fEmission = 1;
42 
43  G4int n_particle = 1;
44  fParticleGun = new G4ParticleGun(n_particle);
45 
46  G4ParticleDefinition* particle =
47  G4ParticleTable::GetParticleTable()->FindParticle("proton");
48  fParticleGun->SetParticleDefinition(particle);
49 
50  fParticleGun->SetParticleEnergy(3*MeV);
51 
52  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
53 
54  fGunMessenger = new PrimaryGeneratorMessenger(this);
55 
56  // Matrix
57 
58  fBeamMatrix = CLHEP::HepMatrix(32,32);
59 
60 }
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 
64 PrimaryGeneratorAction::~PrimaryGeneratorAction()
65 {
66  delete fParticleGun;
67  delete fGunMessenger;
68 }
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 
72 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
73 {
74  G4int numEvent;
75 
76  numEvent=anEvent->GetEventID()+1;
77 
78  G4double x0,y0,z0,theta,phi,xMom0,yMom0,zMom0,e0,de;
79 
80  fShoot=false;
81 
82  x0=0; y0=0; z0=0; theta=0; phi=0; e0=0;
83 
84 // Coefficient computation
85 
86 if (fEmission==1)
87 {
88  fDetector->SetCoef(1);
89  fShoot=true;
90  de = 0;
91 
92  if (numEvent== 1) {theta = -0.00500; phi = -0.00500; de = -0.0050; }
93  if (numEvent== 2) {theta = -0.005/3; phi = -0.00500; de = -0.0050; }
94  if (numEvent== 3) {theta = 0.005/3; phi = -0.00500; de = -0.0050; }
95  if (numEvent== 4) {theta = 0.00500; phi = -0.00500; de = -0.0050; }
96  if (numEvent== 5) {theta = -0.00500; phi = -0.005/3; de = -0.0050; }
97  if (numEvent== 6) {theta = -0.005/3; phi = -0.005/3; de = -0.0050; }
98  if (numEvent== 7) {theta = 0.005/3; phi = -0.005/3; de = -0.0050; }
99  if (numEvent== 8) {theta = 0.00500; phi = -0.005/3; de = -0.0050; }
100  if (numEvent== 9) {theta = -0.00500; phi = 0.005/3; de = -0.0050; }
101  if (numEvent== 10) {theta = -0.005/3; phi = 0.005/3; de = -0.0050; }
102  if (numEvent== 11) {theta = 0.005/3; phi = 0.005/3; de = -0.0050; }
103  if (numEvent== 12) {theta = 0.00500; phi = 0.005/3; de = -0.0050; }
104  if (numEvent== 13) {theta = -0.00500; phi = 0.00500; de = -0.0050; }
105  if (numEvent== 14) {theta = -0.005/3; phi = 0.00500; de = -0.0050; }
106  if (numEvent== 15) {theta = 0.005/3; phi = 0.00500; de = -0.0050; }
107  if (numEvent== 16) {theta = 0.00500; phi = 0.00500; de = -0.0050; }
108  if (numEvent== 17) {theta = -0.00500; phi = -0.00500; de = 0.0050; }
109  if (numEvent== 18) {theta = -0.005/3; phi = -0.00500; de = 0.0050; }
110  if (numEvent== 19) {theta = 0.005/3; phi = -0.00500; de = 0.0050; }
111  if (numEvent== 20) {theta = 0.00500; phi = -0.00500; de = 0.0050; }
112  if (numEvent== 21) {theta = -0.00500; phi = -0.005/3; de = 0.0050; }
113  if (numEvent== 22) {theta = -0.005/3; phi = -0.005/3; de = 0.0050; }
114  if (numEvent== 23) {theta = 0.005/3; phi = -0.005/3; de = 0.0050; }
115  if (numEvent== 24) {theta = 0.00500; phi = -0.005/3; de = 0.0050; }
116  if (numEvent== 25) {theta = -0.00500; phi = 0.005/3; de = 0.0050; }
117  if (numEvent== 26) {theta = -0.005/3; phi = 0.005/3; de = 0.0050; }
118  if (numEvent== 27) {theta = 0.005/3; phi = 0.005/3; de = 0.0050; }
119  if (numEvent== 28) {theta = 0.00500; phi = 0.005/3; de = 0.0050; }
120  if (numEvent== 29) {theta = -0.00500; phi = 0.00500; de = 0.0050; }
121  if (numEvent== 30) {theta = -0.005/3; phi = 0.00500; de = 0.0050; }
122  if (numEvent== 31) {theta = 0.005/3; phi = 0.00500; de = 0.0050; }
123  if (numEvent== 32) {theta = 0.00500; phi = 0.00500; de = 0.0050; }
124 
125  theta=theta*200*fAngleMax*1e-3;
126  phi=phi*200*fAngleMax*1e-3;
127 
128  de=de/100;
129  e0 = 3*MeV + 2*de*3*MeV;
130 
131  x0 = 0;
132  y0 = 0;
133  z0 = -8870*mm;
134 
135  // triplet only
136  //z0 = -3230*mm;
137 
138  G4float cte = 1 ;
139  G4float DE = 100*de;
140 
141  phi = phi * 1000;
142  theta = theta * 1000;
143 
144  // Matrix filling up
145 
146  fBeamMatrix(numEvent,1) = cte;
147 
148  fBeamMatrix(numEvent,2) = theta;
149  fBeamMatrix(numEvent,3) = phi;
150  fBeamMatrix(numEvent,4) = DE;
151 
152  fBeamMatrix(numEvent,5) = theta * theta;
153  fBeamMatrix(numEvent,6) = theta * phi;
154  fBeamMatrix(numEvent,7) = phi * phi;
155  fBeamMatrix(numEvent,8) = theta * DE;
156  fBeamMatrix(numEvent,9) = phi * DE;
157 
158  fBeamMatrix(numEvent,10) = theta * theta * theta;
159  fBeamMatrix(numEvent,11) = theta * theta * phi;
160  fBeamMatrix(numEvent,12) = theta * phi * phi;
161  fBeamMatrix(numEvent,13) = phi * phi * phi;
162  fBeamMatrix(numEvent,14) = theta * theta * DE;
163  fBeamMatrix(numEvent,15) = theta * phi * DE;
164  fBeamMatrix(numEvent,16) = phi * phi * DE;
165 
166  fBeamMatrix(numEvent,17) = theta * theta * theta * phi;
167  fBeamMatrix(numEvent,18) = theta * theta * phi * phi;
168  fBeamMatrix(numEvent,19) = theta * phi * phi * phi;
169  fBeamMatrix(numEvent,20) = theta * theta * theta * DE;
170  fBeamMatrix(numEvent,21) = theta * theta * phi * DE;
171  fBeamMatrix(numEvent,22) = theta * phi * phi * DE;
172  fBeamMatrix(numEvent,23) = phi * phi * phi * DE;
173 
174  fBeamMatrix(numEvent,24) = theta * theta * theta * phi * phi;
175  fBeamMatrix(numEvent,25) = theta * theta * phi * phi * phi;
176  fBeamMatrix(numEvent,26) = theta * theta * theta * phi * DE;
177  fBeamMatrix(numEvent,27) = theta * theta * phi * phi * DE;
178  fBeamMatrix(numEvent,28) = theta * phi * phi * phi * DE;
179 
180  fBeamMatrix(numEvent,29) = theta * theta * theta * phi * phi * phi;
181  fBeamMatrix(numEvent,30) = theta * theta * theta * phi * phi * DE;
182  fBeamMatrix(numEvent,31) = theta * theta * phi * phi * phi * DE;
183 
184  fBeamMatrix(numEvent,32) = theta * theta * theta * phi * phi * phi * DE;
185 
186  //
187 
188  phi = phi / 1000;
189  theta = theta / 1000;
190 
191  /*
192  G4cout << fDetector->GetG1() << G4endl;
193  G4cout << fDetector->GetG2() << G4endl;
194  G4cout << fDetector->GetG3() << G4endl;
195  G4cout << fDetector->GetG4() << G4endl;
196  G4cout << fDetector->GetModel() << G4endl;
197  G4cout << fDetector->GetCoef() << G4endl;
198  G4cout << fDetector->GetGrid() << G4endl;
199  */
200 
201 } // end coefficient
202 
203 // Full beam
204 
205 if (fEmission==2)
206 {
207  fDetector->SetCoef(0);
208  fShoot=false;
209 
210  G4double aR, angle, rR;
211  aR = -1;
212 
213  e0= G4RandGauss::shoot(3*MeV,5.0955e-5*MeV); // AIFIRA ENERGY RESOLUTION
214 
215  while (aR < 0) aR = G4RandGauss::shoot(0.10e-3 , 0.06e-3/2.35) * rad; // old =0.08e-3 displacement
216 
217  angle = G4UniformRand() * 2 * CLHEP::pi *rad;
218 
219  theta = aR * std::cos(angle);
220 
221  phi = aR * std::sin(angle);
222 
223  rR = XYofAngle(aR);
224 
225  x0 = rR*std::cos(angle);
226 
227  y0 = rR*std::sin(angle);
228 
229  x0 = G4RandGauss::shoot(x0,220/2.35) *micrometer;
230 
231  theta = G4RandGauss::shoot(theta,0.03e-3/2.35);
232 
233  y0 = G4RandGauss::shoot(y0,220/2.35) *micrometer;
234 
235  phi = G4RandGauss::shoot(phi,0.03e-3/2.35);
236 
237  z0 = (-9120+250)*mm; //position C0
238 
239  if (std::sqrt(x0*x0+y0*y0)/micrometer<2000) fShoot=true;
240 
241  /*
242  xMom0 = std::sin(theta);
243  yMom0 = std::sin(phi);
244  zMom0 = std::cos(theta);
245  */
246 
247 } // end full beam
248 
249 // shoot
250 
251 if (fShoot)
252 {
253 
254  G4cout << "-> Event= " << numEvent<< ": X0(mm)= " << x0/mm << " X0(mm) = " << y0/mm << " Z0(m) = " << z0/m
255  << " THETA(mrad)= " << theta/mrad << " PHI(mrad)= " << phi/mrad << " E0(MeV)= " << e0/MeV << G4endl;
256 
257  xMom0 = std::sin(theta);
258 
259  yMom0 = std::sin(phi);
260 
261  zMom0 = std::sqrt(1.-xMom0*xMom0-yMom0*yMom0);
262 
263  fParticleGun->SetParticleEnergy(e0);
264 
265  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(xMom0,yMom0,zMom0));
266 
267  fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
268 
269  G4ParticleDefinition* particle=
270  G4ParticleTable::GetParticleTable()->FindParticle("proton");
271 
272  fParticleGun->SetParticleDefinition(particle);
273 
274  fParticleGun->GeneratePrimaryVertex(anEvent);
275 }
276 
277 // end shoot
278 
279 }
280 
281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
282 
283 void PrimaryGeneratorAction::SetEmission(G4int value)
284 {
285  fEmission = value;
286 }
287 
288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
289 
290 G4double PrimaryGeneratorAction::XYofAngle(G4double angle)// returns position in micrometers
291 {
292  return std::pow(20000*angle, 13);
293 }
294 
G4INCL::DeJongSpin::shoot
ThreeVector shoot(const G4int Ap, const G4int Af)
Definition: G4INCLDeJongSpin.cc:67
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition: G4ParticleTable.cc:620
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4float
float G4float
Definition: G4Types.hh:77
G4ParticleGun::SetParticleMomentumDirection
void SetParticleMomentumDirection(G4ParticleMomentum aMomentumDirection)
Definition: G4ParticleGun.hh:96
G4int
int G4int
Definition: G4Types.hh:78
G4ParticleGun::GeneratePrimaryVertex
virtual void GeneratePrimaryVertex(G4Event *evt)
Definition: G4ParticleGun.cc:212
G4Event::GetEventID
G4int GetEventID() const
Definition: G4Event.hh:140
G4VPrimaryGenerator::SetParticlePosition
void SetParticlePosition(G4ThreeVector aPosition)
Definition: G4VPrimaryGenerator.hh:69
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4cout
G4GLOB_DLL std::ostream G4cout
G4ParticleGun::SetParticleEnergy
void SetParticleEnergy(G4double aKineticEnergy)
Definition: G4ParticleGun.cc:132
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:95
value
const XML_Char int const XML_Char * value
Definition: include/expat.h:331
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
python.hepunit.micrometer
int micrometer
Definition: hepunit.py:34
G4ParticleGun::SetParticleDefinition
void SetParticleDefinition(G4ParticleDefinition *aParticleDefinition)
Definition: G4ParticleGun.cc:102
Generated on Wed Apr 30 2014 15:55:22 for Geant4.10 by   doxygen 1.8.7