| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Please cite the following paper if you use
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........oo
34
35 PrimaryGeneratorAction::PrimaryGeneratorAction
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()->FindP
48 fParticleGun->SetParticleDefinition(part
49
50 fParticleGun->SetParticleEnergy(3*MeV);
51
52 fParticleGun->SetParticleMomentumDirection(G
53
54 fGunMessenger = new PrimaryGeneratorMessenge
55
56 // Matrix
57
58 fBeamMatrix = CLHEP::HepMatrix(32,32);
59
60 }
61
62 //....oooOO0OOooo........oooOO0OOooo........oo
63
64 PrimaryGeneratorAction::~PrimaryGeneratorActio
65 {
66 delete fParticleGun;
67 delete fGunMessenger;
68 }
69
70 //....oooOO0OOooo........oooOO0OOooo........oo
71
72 void PrimaryGeneratorAction::GeneratePrimaries
73 {
74 G4int numEvent;
75
76 numEvent=anEvent->GetEventID()+1;
77
78 G4double x0,y0,z0,theta,phi,xMom0,yMom0,zMom
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
93 if (numEvent== 2) {theta = -0.005/3; phi
94 if (numEvent== 3) {theta = 0.005/3; phi
95 if (numEvent== 4) {theta = 0.00500; phi
96 if (numEvent== 5) {theta = -0.00500; phi
97 if (numEvent== 6) {theta = -0.005/3; phi
98 if (numEvent== 7) {theta = 0.005/3; phi
99 if (numEvent== 8) {theta = 0.00500; phi
100 if (numEvent== 9) {theta = -0.00500; phi
101 if (numEvent== 10) {theta = -0.005/3; phi
102 if (numEvent== 11) {theta = 0.005/3; phi
103 if (numEvent== 12) {theta = 0.00500; phi
104 if (numEvent== 13) {theta = -0.00500; phi
105 if (numEvent== 14) {theta = -0.005/3; phi
106 if (numEvent== 15) {theta = 0.005/3; phi
107 if (numEvent== 16) {theta = 0.00500; phi
108 if (numEvent== 17) {theta = -0.00500; phi
109 if (numEvent== 18) {theta = -0.005/3; phi
110 if (numEvent== 19) {theta = 0.005/3; phi
111 if (numEvent== 20) {theta = 0.00500; phi
112 if (numEvent== 21) {theta = -0.00500; phi
113 if (numEvent== 22) {theta = -0.005/3; phi
114 if (numEvent== 23) {theta = 0.005/3; phi
115 if (numEvent== 24) {theta = 0.00500; phi
116 if (numEvent== 25) {theta = -0.00500; phi
117 if (numEvent== 26) {theta = -0.005/3; phi
118 if (numEvent== 27) {theta = 0.005/3; phi
119 if (numEvent== 28) {theta = 0.00500; phi
120 if (numEvent== 29) {theta = -0.00500; phi
121 if (numEvent== 30) {theta = -0.005/3; phi
122 if (numEvent== 31) {theta = 0.005/3; phi
123 if (numEvent== 32) {theta = 0.00500; phi
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 * thet
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 * the
159 fBeamMatrix(numEvent,11) = theta * the
160 fBeamMatrix(numEvent,12) = theta * phi
161 fBeamMatrix(numEvent,13) = phi * phi
162 fBeamMatrix(numEvent,14) = theta * the
163 fBeamMatrix(numEvent,15) = theta * phi
164 fBeamMatrix(numEvent,16) = phi * phi
165
166 fBeamMatrix(numEvent,17) = theta * the
167 fBeamMatrix(numEvent,18) = theta * the
168 fBeamMatrix(numEvent,19) = theta * phi
169 fBeamMatrix(numEvent,20) = theta * the
170 fBeamMatrix(numEvent,21) = theta * the
171 fBeamMatrix(numEvent,22) = theta * phi
172 fBeamMatrix(numEvent,23) = phi * phi *
173
174 fBeamMatrix(numEvent,24) = theta * the
175 fBeamMatrix(numEvent,25) = theta * the
176 fBeamMatrix(numEvent,26) = theta * the
177 fBeamMatrix(numEvent,27) = theta * the
178 fBeamMatrix(numEvent,28) = theta * phi
179
180 fBeamMatrix(numEvent,29) = theta * the
181 fBeamMatrix(numEvent,30) = theta * the
182 fBeamMatrix(numEvent,31) = theta * the
183
184 fBeamMatrix(numEvent,32) = theta * the
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() << G4e
197 G4cout << fDetector->GetCoef() << G4en
198 G4cout << fDetector->GetGrid() << G4en
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);
214
215 while (aR < 0) aR = G4RandGauss::shoot(0.10e
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) *microm
230
231 theta = G4RandGauss::shoot(theta,0.03e-3/2.3
232
233 y0 = G4RandGauss::shoot(y0,220/2.35) *microm
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)
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(m
255 << " THETA(mrad)= " << theta/mrad <<
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(G
266
267 fParticleGun->SetParticlePosition(G4ThreeVec
268
269 G4ParticleDefinition* particle=
270 G4ParticleTable::GetParticleTable()->FindP
271
272 fParticleGun->SetParticleDefinition(particle
273
274 fParticleGun->GeneratePrimaryVertex(anEvent)
275 }
276
277 // end shoot
278
279 }
280
281 //....oooOO0OOooo........oooOO0OOooo........oo
282
283 void PrimaryGeneratorAction::SetEmission(G4int
284 {
285 fEmission = value;
286 }
287
288 //....oooOO0OOooo........oooOO0OOooo........oo
289
290 G4double PrimaryGeneratorAction::XYofAngle(G4d
291 {
292 return std::pow(20000*angle, 13);
293 }
294
295