ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PrimaryGeneratorAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PrimaryGeneratorAction.cc
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