| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // Please cite the following paper if you use 26 // Please cite the following paper if you use this software
27 // Nucl.Instrum.Meth.B260:20-27, 2007 27 // Nucl.Instrum.Meth.B260:20-27, 2007
28 28
29 #include "PrimaryGeneratorAction.hh" 29 #include "PrimaryGeneratorAction.hh"
30 #include "PrimaryGeneratorMessenger.hh" 30 #include "PrimaryGeneratorMessenger.hh"
31 #include "G4SystemOfUnits.hh" 31 #include "G4SystemOfUnits.hh"
32 32
33 //....oooOO0OOooo........oooOO0OOooo........oo 33 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
34 34
35 PrimaryGeneratorAction::PrimaryGeneratorAction 35 PrimaryGeneratorAction::PrimaryGeneratorAction(DetectorConstruction* DC)
36 :fDetector(DC) 36 :fDetector(DC)
37 { 37 {
38 fAngleMax = 0.09; 38 fAngleMax = 0.09;
39 39
40 // Default 40 // Default
41 fEmission = 1; 41 fEmission = 1;
42 42
43 G4int n_particle = 1; 43 G4int n_particle = 1;
44 fParticleGun = new G4ParticleGun(n_particle 44 fParticleGun = new G4ParticleGun(n_particle);
45 45
46 G4ParticleDefinition* particle = 46 G4ParticleDefinition* particle =
47 G4ParticleTable::GetParticleTable()->FindP 47 G4ParticleTable::GetParticleTable()->FindParticle("proton");
48 fParticleGun->SetParticleDefinition(part 48 fParticleGun->SetParticleDefinition(particle);
49 49
50 fParticleGun->SetParticleEnergy(3*MeV); 50 fParticleGun->SetParticleEnergy(3*MeV);
51 51
52 fParticleGun->SetParticleMomentumDirection(G 52 fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
53 53
54 fGunMessenger = new PrimaryGeneratorMessenge 54 fGunMessenger = new PrimaryGeneratorMessenger(this);
55 55
56 // Matrix 56 // Matrix
57 57
58 fBeamMatrix = CLHEP::HepMatrix(32,32); 58 fBeamMatrix = CLHEP::HepMatrix(32,32);
59 59
60 } 60 }
61 61
62 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 63
64 PrimaryGeneratorAction::~PrimaryGeneratorActio 64 PrimaryGeneratorAction::~PrimaryGeneratorAction()
65 { 65 {
66 delete fParticleGun; 66 delete fParticleGun;
67 delete fGunMessenger; 67 delete fGunMessenger;
68 } 68 }
69 69
70 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 71
72 void PrimaryGeneratorAction::GeneratePrimaries 72 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
73 { 73 {
74 G4int numEvent; 74 G4int numEvent;
75 75
76 numEvent=anEvent->GetEventID()+1; 76 numEvent=anEvent->GetEventID()+1;
77 77
78 G4double x0,y0,z0,theta,phi,xMom0,yMom0,zMom 78 G4double x0,y0,z0,theta,phi,xMom0,yMom0,zMom0,e0,de;
79 79
80 fShoot=false; 80 fShoot=false;
81 81
82 x0=0; y0=0; z0=0; theta=0; phi=0; e0=0; 82 x0=0; y0=0; z0=0; theta=0; phi=0; e0=0;
83 83
84 // Coefficient computation 84 // Coefficient computation
85 85
86 if (fEmission==1) 86 if (fEmission==1)
87 { 87 {
88 fDetector->SetCoef(1); 88 fDetector->SetCoef(1);
89 fShoot=true; 89 fShoot=true;
90 de = 0; 90 de = 0;
91 91
92 if (numEvent== 1) {theta = -0.00500; phi 92 if (numEvent== 1) {theta = -0.00500; phi = -0.00500; de = -0.0050; }
93 if (numEvent== 2) {theta = -0.005/3; phi 93 if (numEvent== 2) {theta = -0.005/3; phi = -0.00500; de = -0.0050; }
94 if (numEvent== 3) {theta = 0.005/3; phi 94 if (numEvent== 3) {theta = 0.005/3; phi = -0.00500; de = -0.0050; }
95 if (numEvent== 4) {theta = 0.00500; phi 95 if (numEvent== 4) {theta = 0.00500; phi = -0.00500; de = -0.0050; }
96 if (numEvent== 5) {theta = -0.00500; phi 96 if (numEvent== 5) {theta = -0.00500; phi = -0.005/3; de = -0.0050; }
97 if (numEvent== 6) {theta = -0.005/3; phi 97 if (numEvent== 6) {theta = -0.005/3; phi = -0.005/3; de = -0.0050; }
98 if (numEvent== 7) {theta = 0.005/3; phi 98 if (numEvent== 7) {theta = 0.005/3; phi = -0.005/3; de = -0.0050; }
99 if (numEvent== 8) {theta = 0.00500; phi 99 if (numEvent== 8) {theta = 0.00500; phi = -0.005/3; de = -0.0050; }
100 if (numEvent== 9) {theta = -0.00500; phi 100 if (numEvent== 9) {theta = -0.00500; phi = 0.005/3; de = -0.0050; }
101 if (numEvent== 10) {theta = -0.005/3; phi 101 if (numEvent== 10) {theta = -0.005/3; phi = 0.005/3; de = -0.0050; }
102 if (numEvent== 11) {theta = 0.005/3; phi 102 if (numEvent== 11) {theta = 0.005/3; phi = 0.005/3; de = -0.0050; }
103 if (numEvent== 12) {theta = 0.00500; phi 103 if (numEvent== 12) {theta = 0.00500; phi = 0.005/3; de = -0.0050; }
104 if (numEvent== 13) {theta = -0.00500; phi 104 if (numEvent== 13) {theta = -0.00500; phi = 0.00500; de = -0.0050; }
105 if (numEvent== 14) {theta = -0.005/3; phi 105 if (numEvent== 14) {theta = -0.005/3; phi = 0.00500; de = -0.0050; }
106 if (numEvent== 15) {theta = 0.005/3; phi 106 if (numEvent== 15) {theta = 0.005/3; phi = 0.00500; de = -0.0050; }
107 if (numEvent== 16) {theta = 0.00500; phi 107 if (numEvent== 16) {theta = 0.00500; phi = 0.00500; de = -0.0050; }
108 if (numEvent== 17) {theta = -0.00500; phi 108 if (numEvent== 17) {theta = -0.00500; phi = -0.00500; de = 0.0050; }
109 if (numEvent== 18) {theta = -0.005/3; phi 109 if (numEvent== 18) {theta = -0.005/3; phi = -0.00500; de = 0.0050; }
110 if (numEvent== 19) {theta = 0.005/3; phi 110 if (numEvent== 19) {theta = 0.005/3; phi = -0.00500; de = 0.0050; }
111 if (numEvent== 20) {theta = 0.00500; phi 111 if (numEvent== 20) {theta = 0.00500; phi = -0.00500; de = 0.0050; }
112 if (numEvent== 21) {theta = -0.00500; phi 112 if (numEvent== 21) {theta = -0.00500; phi = -0.005/3; de = 0.0050; }
113 if (numEvent== 22) {theta = -0.005/3; phi 113 if (numEvent== 22) {theta = -0.005/3; phi = -0.005/3; de = 0.0050; }
114 if (numEvent== 23) {theta = 0.005/3; phi 114 if (numEvent== 23) {theta = 0.005/3; phi = -0.005/3; de = 0.0050; }
115 if (numEvent== 24) {theta = 0.00500; phi 115 if (numEvent== 24) {theta = 0.00500; phi = -0.005/3; de = 0.0050; }
116 if (numEvent== 25) {theta = -0.00500; phi 116 if (numEvent== 25) {theta = -0.00500; phi = 0.005/3; de = 0.0050; }
117 if (numEvent== 26) {theta = -0.005/3; phi 117 if (numEvent== 26) {theta = -0.005/3; phi = 0.005/3; de = 0.0050; }
118 if (numEvent== 27) {theta = 0.005/3; phi 118 if (numEvent== 27) {theta = 0.005/3; phi = 0.005/3; de = 0.0050; }
119 if (numEvent== 28) {theta = 0.00500; phi 119 if (numEvent== 28) {theta = 0.00500; phi = 0.005/3; de = 0.0050; }
120 if (numEvent== 29) {theta = -0.00500; phi 120 if (numEvent== 29) {theta = -0.00500; phi = 0.00500; de = 0.0050; }
121 if (numEvent== 30) {theta = -0.005/3; phi 121 if (numEvent== 30) {theta = -0.005/3; phi = 0.00500; de = 0.0050; }
122 if (numEvent== 31) {theta = 0.005/3; phi 122 if (numEvent== 31) {theta = 0.005/3; phi = 0.00500; de = 0.0050; }
123 if (numEvent== 32) {theta = 0.00500; phi 123 if (numEvent== 32) {theta = 0.00500; phi = 0.00500; de = 0.0050; }
124 124
125 theta=theta*200*fAngleMax*1e-3; 125 theta=theta*200*fAngleMax*1e-3;
126 phi=phi*200*fAngleMax*1e-3; 126 phi=phi*200*fAngleMax*1e-3;
127 127
128 de=de/100; 128 de=de/100;
129 e0 = 3*MeV + 2*de*3*MeV; 129 e0 = 3*MeV + 2*de*3*MeV;
130 130
131 x0 = 0; 131 x0 = 0;
132 y0 = 0; 132 y0 = 0;
133 z0 = -8870*mm; 133 z0 = -8870*mm;
134 134
135 // triplet only 135 // triplet only
136 //z0 = -3230*mm; 136 //z0 = -3230*mm;
137 137
138 G4float cte = 1 ; 138 G4float cte = 1 ;
139 G4float DE = 100*de; 139 G4float DE = 100*de;
140 140
141 phi = phi * 1000; 141 phi = phi * 1000;
142 theta = theta * 1000; 142 theta = theta * 1000;
143 143
144 // Matrix filling up 144 // Matrix filling up
145 145
146 fBeamMatrix(numEvent,1) = cte; 146 fBeamMatrix(numEvent,1) = cte;
147 147
148 fBeamMatrix(numEvent,2) = theta; 148 fBeamMatrix(numEvent,2) = theta;
149 fBeamMatrix(numEvent,3) = phi; 149 fBeamMatrix(numEvent,3) = phi;
150 fBeamMatrix(numEvent,4) = DE; 150 fBeamMatrix(numEvent,4) = DE;
151 151
152 fBeamMatrix(numEvent,5) = theta * thet 152 fBeamMatrix(numEvent,5) = theta * theta;
153 fBeamMatrix(numEvent,6) = theta * phi; 153 fBeamMatrix(numEvent,6) = theta * phi;
154 fBeamMatrix(numEvent,7) = phi * phi; 154 fBeamMatrix(numEvent,7) = phi * phi;
155 fBeamMatrix(numEvent,8) = theta * DE; 155 fBeamMatrix(numEvent,8) = theta * DE;
156 fBeamMatrix(numEvent,9) = phi * DE; 156 fBeamMatrix(numEvent,9) = phi * DE;
157 157
158 fBeamMatrix(numEvent,10) = theta * the 158 fBeamMatrix(numEvent,10) = theta * theta * theta;
159 fBeamMatrix(numEvent,11) = theta * the 159 fBeamMatrix(numEvent,11) = theta * theta * phi;
160 fBeamMatrix(numEvent,12) = theta * phi 160 fBeamMatrix(numEvent,12) = theta * phi * phi;
161 fBeamMatrix(numEvent,13) = phi * phi 161 fBeamMatrix(numEvent,13) = phi * phi * phi;
162 fBeamMatrix(numEvent,14) = theta * the 162 fBeamMatrix(numEvent,14) = theta * theta * DE;
163 fBeamMatrix(numEvent,15) = theta * phi 163 fBeamMatrix(numEvent,15) = theta * phi * DE;
164 fBeamMatrix(numEvent,16) = phi * phi 164 fBeamMatrix(numEvent,16) = phi * phi * DE;
165 165
166 fBeamMatrix(numEvent,17) = theta * the 166 fBeamMatrix(numEvent,17) = theta * theta * theta * phi;
167 fBeamMatrix(numEvent,18) = theta * the 167 fBeamMatrix(numEvent,18) = theta * theta * phi * phi;
168 fBeamMatrix(numEvent,19) = theta * phi 168 fBeamMatrix(numEvent,19) = theta * phi * phi * phi;
169 fBeamMatrix(numEvent,20) = theta * the 169 fBeamMatrix(numEvent,20) = theta * theta * theta * DE;
170 fBeamMatrix(numEvent,21) = theta * the 170 fBeamMatrix(numEvent,21) = theta * theta * phi * DE;
171 fBeamMatrix(numEvent,22) = theta * phi 171 fBeamMatrix(numEvent,22) = theta * phi * phi * DE;
172 fBeamMatrix(numEvent,23) = phi * phi * 172 fBeamMatrix(numEvent,23) = phi * phi * phi * DE;
173 173
174 fBeamMatrix(numEvent,24) = theta * the 174 fBeamMatrix(numEvent,24) = theta * theta * theta * phi * phi;
175 fBeamMatrix(numEvent,25) = theta * the 175 fBeamMatrix(numEvent,25) = theta * theta * phi * phi * phi;
176 fBeamMatrix(numEvent,26) = theta * the 176 fBeamMatrix(numEvent,26) = theta * theta * theta * phi * DE;
177 fBeamMatrix(numEvent,27) = theta * the 177 fBeamMatrix(numEvent,27) = theta * theta * phi * phi * DE;
178 fBeamMatrix(numEvent,28) = theta * phi 178 fBeamMatrix(numEvent,28) = theta * phi * phi * phi * DE;
179 179
180 fBeamMatrix(numEvent,29) = theta * the 180 fBeamMatrix(numEvent,29) = theta * theta * theta * phi * phi * phi;
181 fBeamMatrix(numEvent,30) = theta * the 181 fBeamMatrix(numEvent,30) = theta * theta * theta * phi * phi * DE;
182 fBeamMatrix(numEvent,31) = theta * the 182 fBeamMatrix(numEvent,31) = theta * theta * phi * phi * phi * DE;
183 183
184 fBeamMatrix(numEvent,32) = theta * the 184 fBeamMatrix(numEvent,32) = theta * theta * theta * phi * phi * phi * DE;
185 185
186 // 186 //
187 187
188 phi = phi / 1000; 188 phi = phi / 1000;
189 theta = theta / 1000; 189 theta = theta / 1000;
190 190
191 /* 191 /*
192 G4cout << fDetector->GetG1() << G4endl 192 G4cout << fDetector->GetG1() << G4endl;
193 G4cout << fDetector->GetG2() << G4endl 193 G4cout << fDetector->GetG2() << G4endl;
194 G4cout << fDetector->GetG3() << G4endl 194 G4cout << fDetector->GetG3() << G4endl;
195 G4cout << fDetector->GetG4() << G4endl 195 G4cout << fDetector->GetG4() << G4endl;
196 G4cout << fDetector->GetModel() << G4e 196 G4cout << fDetector->GetModel() << G4endl;
197 G4cout << fDetector->GetCoef() << G4en 197 G4cout << fDetector->GetCoef() << G4endl;
198 G4cout << fDetector->GetGrid() << G4en 198 G4cout << fDetector->GetGrid() << G4endl;
199 */ 199 */
200 200
201 } // end coefficient 201 } // end coefficient
202 202
203 // Full beam 203 // Full beam
204 204
205 if (fEmission==2) 205 if (fEmission==2)
206 { 206 {
207 fDetector->SetCoef(0); 207 fDetector->SetCoef(0);
208 fShoot=false; 208 fShoot=false;
209 209
210 G4double aR, angle, rR; 210 G4double aR, angle, rR;
211 aR = -1; 211 aR = -1;
212 212
213 e0= G4RandGauss::shoot(3*MeV,5.0955e-5*MeV); 213 e0= G4RandGauss::shoot(3*MeV,5.0955e-5*MeV); // AIFIRA ENERGY RESOLUTION
214 214
215 while (aR < 0) aR = G4RandGauss::shoot(0.10e 215 while (aR < 0) aR = G4RandGauss::shoot(0.10e-3 , 0.06e-3/2.35) * rad; // old =0.08e-3 displacement
216 216
217 angle = G4UniformRand() * 2 * CLHEP::pi *rad 217 angle = G4UniformRand() * 2 * CLHEP::pi *rad;
218 218
219 theta = aR * std::cos(angle); 219 theta = aR * std::cos(angle);
220 220
221 phi = aR * std::sin(angle); 221 phi = aR * std::sin(angle);
222 222
223 rR = XYofAngle(aR); 223 rR = XYofAngle(aR);
224 224
225 x0 = rR*std::cos(angle); 225 x0 = rR*std::cos(angle);
226 226
227 y0 = rR*std::sin(angle); 227 y0 = rR*std::sin(angle);
228 228
229 x0 = G4RandGauss::shoot(x0,220/2.35) *microm 229 x0 = G4RandGauss::shoot(x0,220/2.35) *micrometer;
230 230
231 theta = G4RandGauss::shoot(theta,0.03e-3/2.3 231 theta = G4RandGauss::shoot(theta,0.03e-3/2.35);
232 232
233 y0 = G4RandGauss::shoot(y0,220/2.35) *microm 233 y0 = G4RandGauss::shoot(y0,220/2.35) *micrometer;
234 234
235 phi = G4RandGauss::shoot(phi,0.03e-3/2.35); 235 phi = G4RandGauss::shoot(phi,0.03e-3/2.35);
236 236
237 z0 = (-9120+250)*mm; //position C0 237 z0 = (-9120+250)*mm; //position C0
238 238
239 if (std::sqrt(x0*x0+y0*y0)/micrometer<2000) 239 if (std::sqrt(x0*x0+y0*y0)/micrometer<2000) fShoot=true;
240 240
241 /* 241 /*
242 xMom0 = std::sin(theta); 242 xMom0 = std::sin(theta);
243 yMom0 = std::sin(phi); 243 yMom0 = std::sin(phi);
244 zMom0 = std::cos(theta); 244 zMom0 = std::cos(theta);
245 */ 245 */
246 246
247 } // end full beam 247 } // end full beam
248 248
249 // shoot 249 // shoot
250 250
251 if (fShoot) 251 if (fShoot)
252 { 252 {
253 253
254 G4cout << "-> Event= " << numEvent<< ": X0(m 254 G4cout << "-> Event= " << numEvent<< ": X0(mm)= " << x0/mm << " X0(mm) = " << y0/mm << " Z0(m) = " << z0/m
255 << " THETA(mrad)= " << theta/mrad << 255 << " THETA(mrad)= " << theta/mrad << " PHI(mrad)= " << phi/mrad << " E0(MeV)= " << e0/MeV << G4endl;
256 256
257 xMom0 = std::sin(theta); 257 xMom0 = std::sin(theta);
258 258
259 yMom0 = std::sin(phi); 259 yMom0 = std::sin(phi);
260 260
261 zMom0 = std::sqrt(1.-xMom0*xMom0-yMom0*yMom0 261 zMom0 = std::sqrt(1.-xMom0*xMom0-yMom0*yMom0);
262 262
263 fParticleGun->SetParticleEnergy(e0); 263 fParticleGun->SetParticleEnergy(e0);
264 264
265 fParticleGun->SetParticleMomentumDirection(G 265 fParticleGun->SetParticleMomentumDirection(G4ThreeVector(xMom0,yMom0,zMom0));
266 266
267 fParticleGun->SetParticlePosition(G4ThreeVec 267 fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
268 268
269 G4ParticleDefinition* particle= 269 G4ParticleDefinition* particle=
270 G4ParticleTable::GetParticleTable()->FindP 270 G4ParticleTable::GetParticleTable()->FindParticle("proton");
271 271
272 fParticleGun->SetParticleDefinition(particle 272 fParticleGun->SetParticleDefinition(particle);
273 273
274 fParticleGun->GeneratePrimaryVertex(anEvent) 274 fParticleGun->GeneratePrimaryVertex(anEvent);
275 } 275 }
276 276
277 // end shoot 277 // end shoot
278 278
279 } 279 }
280 280
281 //....oooOO0OOooo........oooOO0OOooo........oo 281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
282 282
283 void PrimaryGeneratorAction::SetEmission(G4int 283 void PrimaryGeneratorAction::SetEmission(G4int value)
284 { 284 {
285 fEmission = value; 285 fEmission = value;
286 } 286 }
287 287
288 //....oooOO0OOooo........oooOO0OOooo........oo 288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
289 289
290 G4double PrimaryGeneratorAction::XYofAngle(G4d 290 G4double PrimaryGeneratorAction::XYofAngle(G4double angle)// returns position in micrometers
291 { 291 {
292 return std::pow(20000*angle, 13); 292 return std::pow(20000*angle, 13);
293 } 293 }
294 294
295 295