Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 // -------------------------------------------------------------------
27 // $Id$
28 // -------------------------------------------------------------------
29 
30 #include "PrimaryGeneratorAction.hh"
31 #include "PrimaryGeneratorMessenger.hh"
32 #include "G4SystemOfUnits.hh"
33 
34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
35 
36 PrimaryGeneratorAction::PrimaryGeneratorAction(DetectorConstruction* DC)
37 :detector(DC)
38 {
39  angleMax = 0.09;
40  emission =0;
41 
42  G4int n_particle = 1;
43  particleGun = new G4ParticleGun(n_particle);
44 
45  G4ParticleDefinition* particle =
46  G4ParticleTable::GetParticleTable()->FindParticle("proton");
47  particleGun->SetParticleDefinition(particle);
48 
49  particleGun->SetParticleEnergy(3*MeV);
50 
51  particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
52 
53  gunMessenger = new PrimaryGeneratorMessenger(this);
54 
55  // Matrix
56 
57  beamMatrix = CLHEP::HepMatrix(32,32);
58 
59 }
60 
61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
62 
63 PrimaryGeneratorAction::~PrimaryGeneratorAction()
64 {
65  delete particleGun;
66  delete gunMessenger;
67 }
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70 
71 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
72 {
73  G4int numEvent;
74  numEvent=anEvent->GetEventID()+1;
75  G4double x0,y0,z0,theta,phi,xMom0,yMom0,zMom0,e0,de;
76  shoot=false;
77 
78  x0=0; y0=0; z0=0; theta=0; phi=0; e0=0;
79 
80 // Coefficient computation
81 
82 if (emission==1)
83 {
84  detector->SetCoef();
85  shoot=true;
86  de = 0;
87 
88  if (numEvent== 1) {theta = -0.00500; phi = -0.00500; de = -0.0050; }
89  if (numEvent== 2) {theta = -0.005/3; phi = -0.00500; de = -0.0050; }
90  if (numEvent== 3) {theta = 0.005/3; phi = -0.00500; de = -0.0050; }
91  if (numEvent== 4) {theta = 0.00500; phi = -0.00500; de = -0.0050; }
92  if (numEvent== 5) {theta = -0.00500; phi = -0.005/3; de = -0.0050; }
93  if (numEvent== 6) {theta = -0.005/3; phi = -0.005/3; de = -0.0050; }
94  if (numEvent== 7) {theta = 0.005/3; phi = -0.005/3; de = -0.0050; }
95  if (numEvent== 8) {theta = 0.00500; phi = -0.005/3; de = -0.0050; }
96  if (numEvent== 9) {theta = -0.00500; phi = 0.005/3; de = -0.0050; }
97  if (numEvent== 10) {theta = -0.005/3; phi = 0.005/3; de = -0.0050; }
98  if (numEvent== 11) {theta = 0.005/3; phi = 0.005/3; de = -0.0050; }
99  if (numEvent== 12) {theta = 0.00500; phi = 0.005/3; de = -0.0050; }
100  if (numEvent== 13) {theta = -0.00500; phi = 0.00500; de = -0.0050; }
101  if (numEvent== 14) {theta = -0.005/3; phi = 0.00500; de = -0.0050; }
102  if (numEvent== 15) {theta = 0.005/3; phi = 0.00500; de = -0.0050; }
103  if (numEvent== 16) {theta = 0.00500; phi = 0.00500; de = -0.0050; }
104  if (numEvent== 17) {theta = -0.00500; phi = -0.00500; de = 0.0050; }
105  if (numEvent== 18) {theta = -0.005/3; phi = -0.00500; de = 0.0050; }
106  if (numEvent== 19) {theta = 0.005/3; phi = -0.00500; de = 0.0050; }
107  if (numEvent== 20) {theta = 0.00500; phi = -0.00500; de = 0.0050; }
108  if (numEvent== 21) {theta = -0.00500; phi = -0.005/3; de = 0.0050; }
109  if (numEvent== 22) {theta = -0.005/3; phi = -0.005/3; de = 0.0050; }
110  if (numEvent== 23) {theta = 0.005/3; phi = -0.005/3; de = 0.0050; }
111  if (numEvent== 24) {theta = 0.00500; phi = -0.005/3; de = 0.0050; }
112  if (numEvent== 25) {theta = -0.00500; phi = 0.005/3; de = 0.0050; }
113  if (numEvent== 26) {theta = -0.005/3; phi = 0.005/3; de = 0.0050; }
114  if (numEvent== 27) {theta = 0.005/3; phi = 0.005/3; de = 0.0050; }
115  if (numEvent== 28) {theta = 0.00500; phi = 0.005/3; de = 0.0050; }
116  if (numEvent== 29) {theta = -0.00500; phi = 0.00500; de = 0.0050; }
117  if (numEvent== 30) {theta = -0.005/3; phi = 0.00500; de = 0.0050; }
118  if (numEvent== 31) {theta = 0.005/3; phi = 0.00500; de = 0.0050; }
119  if (numEvent== 32) {theta = 0.00500; phi = 0.00500; de = 0.0050; }
120 
121  theta=theta*200*angleMax*1e-3;
122  phi=phi*200*angleMax*1e-3;
123 
124  de=de/100;
125  e0 = 3*MeV + 2*de*3*MeV;
126 
127  x0 = 0;
128  y0 = 0;
129  z0 = -8870*mm;
130 
131  // triplet only
132  //z0 = -3230*mm;
133 
134  G4float cte = 1 ;
135  G4float DE = 100*de;
136 
137  phi = phi * 1000;
138  theta = theta * 1000;
139 
140  // Matrix filling up
141 
142  beamMatrix(numEvent,1) = cte;
143 
144  beamMatrix(numEvent,2) = theta;
145  beamMatrix(numEvent,3) = phi;
146  beamMatrix(numEvent,4) = DE;
147 
148  beamMatrix(numEvent,5) = theta * theta;
149  beamMatrix(numEvent,6) = theta * phi;
150  beamMatrix(numEvent,7) = phi * phi;
151  beamMatrix(numEvent,8) = theta * DE;
152  beamMatrix(numEvent,9) = phi * DE;
153 
154  beamMatrix(numEvent,10) = theta * theta * theta;
155  beamMatrix(numEvent,11) = theta * theta * phi;
156  beamMatrix(numEvent,12) = theta * phi * phi;
157  beamMatrix(numEvent,13) = phi * phi * phi;
158  beamMatrix(numEvent,14) = theta * theta * DE;
159  beamMatrix(numEvent,15) = theta * phi * DE;
160  beamMatrix(numEvent,16) = phi * phi * DE;
161 
162  beamMatrix(numEvent,17) = theta * theta * theta * phi;
163  beamMatrix(numEvent,18) = theta * theta * phi * phi;
164  beamMatrix(numEvent,19) = theta * phi * phi * phi;
165  beamMatrix(numEvent,20) = theta * theta * theta * DE;
166  beamMatrix(numEvent,21) = theta * theta * phi * DE;
167  beamMatrix(numEvent,22) = theta * phi * phi * DE;
168  beamMatrix(numEvent,23) = phi * phi * phi * DE;
169 
170  beamMatrix(numEvent,24) = theta * theta * theta * phi * phi;
171  beamMatrix(numEvent,25) = theta * theta * phi * phi * phi;
172  beamMatrix(numEvent,26) = theta * theta * theta * phi * DE;
173  beamMatrix(numEvent,27) = theta * theta * phi * phi * DE;
174  beamMatrix(numEvent,28) = theta * phi * phi * phi * DE;
175 
176  beamMatrix(numEvent,29) = theta * theta * theta * phi * phi * phi;
177  beamMatrix(numEvent,30) = theta * theta * theta * phi * phi * DE;
178  beamMatrix(numEvent,31) = theta * theta * phi * phi * phi * DE;
179 
180  beamMatrix(numEvent,32) = theta * theta * theta * phi * phi * phi * DE;
181 
182  //
183 
184  phi = phi / 1000;
185  theta = theta / 1000;
186 
187 } // end coefficient
188 
189 // Full beam
190 
191 if (emission==2)
192 {
193  shoot=false;
194  G4double aR, angle, rR;
195  aR = -1;
196  e0= G4RandGauss::shoot(3*MeV,5.0955e-5*MeV); // AIFIRA ENERGY RESOLUTION
197  while (aR < 0) aR = G4RandGauss::shoot(0.10e-3 , 0.06e-3/2.35) * rad; // old =0.08e-3 displacement
198  angle = G4UniformRand() * 2 * CLHEP::pi *rad;
199  theta = aR * std::cos(angle);
200  phi = aR * std::sin(angle);
201  rR = XYofAngle(aR);
202  x0 = rR*std::cos(angle);
203  y0 = rR*std::sin(angle);
204  x0 = G4RandGauss::shoot(x0,220/2.35) *micrometer;
205  theta = G4RandGauss::shoot(theta,0.03e-3/2.35);
206  y0 = G4RandGauss::shoot(y0,220/2.35) *micrometer;
207  phi = G4RandGauss::shoot(phi,0.03e-3/2.35);
208  z0 = (-9120+250)*mm; //position C0
209  if (std::sqrt(x0*x0+y0*y0)/micrometer<2000) shoot=true;
210 
211  /*
212  xMom0 = std::sin(theta);
213  yMom0 = std::sin(phi);
214  zMom0 = std::cos(theta);
215  */
216 
217 } // end full beam
218 
219 // shoot
220 
221 if (shoot)
222 {
223 
224  G4cout << "-> Event= " << numEvent<< ": X0(mm)= " << x0/mm << " X0(mm) = " << y0/mm << " Z0(m) = " << z0/m
225  << " THETA(mrad)= " << theta/mrad << " PHI(mrad)= " << phi/mrad << " E0(MeV)= " << e0/MeV << G4endl;
226 
227  xMom0 = std::sin(theta);
228  yMom0 = std::sin(phi);
229  zMom0 = std::sqrt(1.-xMom0*xMom0-yMom0*yMom0);
230 
231  particleGun->SetParticleEnergy(e0);
232 
233  particleGun->SetParticleMomentumDirection(G4ThreeVector(xMom0,yMom0,zMom0));
234 
235  particleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
236 
237  G4ParticleDefinition* particle=
238  G4ParticleTable::GetParticleTable()->FindParticle("proton");
239 
240  particleGun->SetParticleDefinition(particle);
241 
242  particleGun->GeneratePrimaryVertex(anEvent);
243 }
244 
245 // end shoot
246 
247 }
248 
249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
250 
251 void PrimaryGeneratorAction::SetEmission(G4int value)
252 {
253  emission = value;
254 }
255 
256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
257 
258 G4double PrimaryGeneratorAction::XYofAngle(G4double angle)// returns position in micrometers
259 {
260  return std::pow(20000*angle, 13);
261 }
262 
Generated on Sat May 25 2013 14:31:58 for Geant4 by   doxygen 1.8.4