ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CCalPrimaryGeneratorAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CCalPrimaryGeneratorAction.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 //
27 // File: CCalPrimaryGeneratorAction.cc
28 // Description: CCalPrimaryGeneratorAction Sets up particle beam
30 
31 #include <CLHEP/Random/RandFlat.h>
32 
33 #include "CCalPrimaryGeneratorAction.hh"
34 #include "CCalPrimaryGeneratorMessenger.hh"
35 #include "G4HEPEvtInterface.hh"
36 
37 #include "G4Exp.hh"
38 #include "G4PhysicalConstants.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4Event.hh"
41 #include "G4ParticleGun.hh"
42 #include "G4ParticleTable.hh"
43 #include "G4ParticleDefinition.hh"
44 #include "G4HEPEvtInterface.hh"
45 #include "G4RunManager.hh"
46 
47 //#define debug
48 
49 CCalPrimaryGeneratorAction::CCalPrimaryGeneratorAction(): particleGun(0),
50  generatorInput(singleFixed), verboseLevel(0), n_particle(1),
51  particleName("pi-"), particleEnergy(100*GeV), particlePosition(0.,0.,0.),
52  particleDir(1.,1.,0.1), isInitialized(0), scanSteps(0) {
53 
54  //Initialise the messenger
55  gunMessenger = new CCalPrimaryGeneratorMessenger(this);
56 
57  //Default settings:
58  SetMinimumEnergy(1.*GeV);
59  SetMaximumEnergy(1.*TeV);
60  SetMinimumEta(0.);
61  SetMaximumEta(3.5);
62  SetMinimumPhi(0.*degree);
63  SetMaximumPhi(360.*degree);
64  SetStepsPhi(1);
65  SetStepsEta(1);
66 
67  // Number of particles per gun
68  particleGun = new G4ParticleGun(n_particle);
69 
70  // Getting particle definition
71  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
72  G4ParticleDefinition* particle = particleTable->FindParticle(particleName);
73  particleGun->SetParticleDefinition(particle);
74 
75  // Setting particle properties
76  particleGun->SetParticleEnergy(particleEnergy);
77  particleGun->SetParticleMomentumDirection(particleDir);
78  particleGun->SetParticlePosition(particlePosition);
79  print(0);
80 }
81 
82 CCalPrimaryGeneratorAction::~CCalPrimaryGeneratorAction() {
83  if (gunMessenger)
84  delete gunMessenger;
85  if (particleGun)
86  delete particleGun;
87 }
88 
89 void CCalPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) {
90 
91  if (isInitialized == 0) initialize();
92 
93  if (generatorInput == singleRandom) {
94  particleEnergy = CLHEP::RandFlat::shoot(energyMin,energyMax);
95  particleGun->SetParticleEnergy(particleEnergy);
96 
97  G4double eta = CLHEP::RandFlat::shoot(etaMin,etaMax);
98  G4double phi = CLHEP::RandFlat::shoot(phiMin,phiMax);
99  G4double theta = std::atan(G4Exp(-eta))*2.;
100  G4double randomX = std::sin(theta)*std::cos(phi);
101  G4double randomY = std::sin(theta)*std::sin(phi);
102  G4double randomZ = std::cos(theta);
103 
104  particleDir = G4ThreeVector(randomX,randomY,randomZ);
105  particleGun->SetParticleMomentumDirection(particleDir);
106  if (verboseLevel >= 2 ) {
107  G4cout << "Energy " << particleEnergy/GeV << " GeV; Theta "
108  << theta/deg << " degree; Phi " << phi/deg << " degree" << G4endl;
109  G4cout << "Shooting in " << particleDir << " direction "<< G4endl;
110  }
111  } else if (generatorInput == singleScan) {
112  G4double scanEtaStep, scanPhiStep;
113  if (scanSteps == 0) {
114  scanPhiStep = (phiMax - phiMin) / phiSteps;
115  phiValue = phiMin - scanPhiStep; //first step is going to change scanCurrentPhiValue
116  etaValue = etaMin;
117  }
118 
119  scanEtaStep = (etaMax - etaMin) / etaSteps;
120  scanPhiStep = (phiMax - phiMin) / phiSteps;
121 #ifdef debug
122  if (verboseLevel > 2 ) {
123  G4cout << " scanEtaStep " << scanEtaStep << " # of Steps " << etaSteps
124  << G4endl;
125  G4cout << " scanPhiStep " << scanPhiStep << " # of Steps " << phiSteps
126  << G4endl;
127  }
128 #endif
129 
130  //----- First scan in phi, then in eta
131  if (phiMax - phiValue < 1.E-6 * scanPhiStep) { // !only <= 1.E6 steps allowed
132  if (etaMax - etaValue < 1.E-6 * scanEtaStep) { // !only <= 1.E6 steps allowed
133  G4cout << " Scan completed!" << G4endl;
134  return;
135  } else {
136  etaValue += scanEtaStep;
137  phiValue = phiMin;
138  }
139  } else {
140  phiValue += scanPhiStep;
141  }
142  G4double theta = std::atan(G4Exp(-etaValue))*2.;
143 
144  G4double scanX = std::sin(theta)*std::cos(phiValue);
145  G4double scanY = std::sin(theta)*std::sin(phiValue);
146  G4double scanZ = std::cos(theta);
147  if (verboseLevel >= 2 ) {
148  G4cout << "Scan eta " << etaValue << " Phi " << phiValue/deg
149  << " theta " << theta/deg << G4endl;
150  }
151  particleDir = G4ThreeVector(scanX,scanY,scanZ);
152  particleGun->SetParticleMomentumDirection(particleDir);
153 #ifdef debug
154  if (verboseLevel > 2 ) {
155  G4cout << "Shooting in " << particleDir << " direction "<< G4endl;
156  }
157 #endif
158  scanSteps++;
159  }
160 
161  // Generate GEANT4 primary vertex
162  particleGun->GeneratePrimaryVertex(anEvent);
163 }
164 
165 
166 void CCalPrimaryGeneratorAction::SetVerboseLevel(G4int val){
167  verboseLevel = val;
168 }
169 
170 
171 void CCalPrimaryGeneratorAction::SetRandom(G4String val) {
172 
173  if (val=="on") {
174  generatorInput = singleRandom;
175  print (1);
176  } else {
177  generatorInput = singleFixed;
178  print (1);
179  }
180 }
181 
182 
183 void CCalPrimaryGeneratorAction::SetScan(G4String val) {
184 
185  if (val=="on") {
186  generatorInput = singleScan;
187  scanSteps = 0;
188  print (1);
189  } else {
190  generatorInput = singleFixed;
191  print (1);
192  }
193 }
194 
195 
196 void CCalPrimaryGeneratorAction::SetMinimumEnergy(G4double p){
197 
198  if (p <= 0.) {
199  G4cerr << "CCalPrimaryGeneratorAction::SetMinimumEnergy: value " << p/GeV
200  << "GeV is out of bounds, it will not be used" << G4endl;
201  G4cerr << " Should be >0. Please check" << G4endl;
202  } else {
203  energyMin = p;
204 #ifdef debug
205  if (verboseLevel >= 1 ) {
206  G4cout << " CCalPrimaryGeneratorAction: setting min. value of energy to "
207  << energyMin/GeV << " GeV " << G4endl;
208  }
209 #endif
210  }
211 }
212 
213 
214 void CCalPrimaryGeneratorAction::SetMaximumEnergy(G4double p){
215 
216  if (p <= 0.) {
217  G4cerr << "CCalPrimaryGeneratorAction::SetMaximumEnergy: value " << p/GeV
218  << "GeV is out of bounds, it will not be used" << G4endl;
219  G4cerr << " Should be >0. Please check" << G4endl;
220  } else {
221  energyMax = p;
222 #ifdef debug
223  if (verboseLevel >= 1 ) {
224  G4cout << " CCalPrimaryGeneratorAction: setting max. value of energy to "
225  << energyMax/GeV << " GeV " << G4endl;
226  }
227 #endif
228  }
229 }
230 
231 
232 void CCalPrimaryGeneratorAction::SetMinimumPhi(G4double p){
233 
234  if (std::fabs(p)>2.*pi) {
235  G4cerr << "CCalPrimaryGeneratorAction::SetMinimumPhi: setting value quite "
236  << "large" << G4endl;
237  G4cerr << " Should be given in radians - Please check" << G4endl;
238  } else {
239  phiMin = std::fabs(p);
240 #ifdef debug
241  if (verboseLevel >= 1 ) {
242  G4cout << " CCalPrimaryGeneratorAction: setting min. value of phi to "
243  << phiMin << G4endl;
244  }
245 #endif
246  }
247 }
248 
249 
250 void CCalPrimaryGeneratorAction::SetMaximumPhi(G4double p){
251 
252  if (std::fabs(p)>2.*pi) {
253  G4cerr << "CCalPrimaryGeneratorAction::SetMaximumPhi: setting value quite "
254  << "large" << G4endl;
255  G4cerr << " Should be given in radians - Please check" << G4endl;
256  } else {
257  phiMax = std::fabs(p);
258 #ifdef debug
259  if (verboseLevel >= 1 ) {
260  G4cout << " CCalPrimaryGeneratorAction: setting max. value of phi to "
261  << phiMax << G4endl;
262  }
263 #endif
264  }
265 }
266 
267 
268 void CCalPrimaryGeneratorAction::SetStepsPhi(G4int val){
269 
270  if (val <= 0) {
271  G4cerr << "CCalPrimaryGeneratorAction::SetStepsPhi: value " << val
272  << " is out of bounds, it will not be used" << G4endl;
273  G4cerr << " Should be > 0 Please check" << G4endl;
274  } else {
275  phiSteps = val;
276 #ifdef debug
277  if (verboseLevel >= 1 ) {
278  G4cout << " CCalPrimaryGeneratorAction: setting no. of steps in phi to "
279  << phiSteps << G4endl;
280  }
281 #endif
282  }
283 }
284 
285 
286 void CCalPrimaryGeneratorAction::SetMinimumEta(G4double p){
287 
288  etaMin = p;
289 #ifdef debug
290  if (verboseLevel >= 1 ) {
291  G4cout << " CCalPrimaryGeneratorAction: setting min. value of eta to "
292  << etaMin << G4endl;
293  }
294 #endif
295 }
296 
297 
298 void CCalPrimaryGeneratorAction::SetMaximumEta(G4double p){
299 
300  etaMax = p;
301 #ifdef debug
302  if (verboseLevel >= 1 ) {
303  G4cout << " CCalPrimaryGeneratorAction: setting max. value of eta to "
304  << etaMax << G4endl;
305  }
306 #endif
307 }
308 
309 
310 void CCalPrimaryGeneratorAction::SetStepsEta(G4int val){
311 
312  if (val <= 0) {
313  G4cerr<<"CCalPrimaryGeneratorAction::SetStepsEta: value " << val << " is out of bounds, it will not be used"<<G4endl;
314  G4cerr<<" Should be > 0 Please check"<<G4endl;
315  } else {
316  etaSteps = val;
317 #ifdef debug
318  if (verboseLevel >= 1 ) {
319  G4cout << " CCalPrimaryGeneratorAction: setting no. of steps in eta to "
320  << etaSteps << G4endl;
321  }
322 #endif
323  }
324 }
325 
326 void CCalPrimaryGeneratorAction::SetGunPosition(const G4ThreeVector & pos) const {
327 
328  particleGun->SetParticlePosition(pos);
329 }
330 
331 void CCalPrimaryGeneratorAction::SetRunNo(G4int val){
332  G4RunManager::GetRunManager()->SetRunIDCounter( val );
333 }
334 
335 void CCalPrimaryGeneratorAction::initialize(){
336 
337  isInitialized = 1;
338 
339  print (1);
340 }
341 
342 
343 void CCalPrimaryGeneratorAction::print(G4int val){
344 
345  if (verboseLevel >= val) {
346 
347  if (generatorInput == singleRandom) {
348  G4cout << G4endl
349  << "**********************************************************************" << G4endl
350  << "* *" << G4endl
351  << "* CCalPrimaryGeneratorAction DEFAULT Random Energy/Direction setting:*" << G4endl
352  << "* *" << G4endl
353  << "* *" << G4endl
354  << "* Energy in [ "<< energyMin/GeV << " - " << energyMax/GeV << "] (GeV) "<< G4endl
355  << "* Phi angle in [ "<< phiMin << " - " << phiMax << "] (rad) "<< G4endl
356  << "* [ "<< phiMin/degree << " - " << phiMax/degree << "] (deg) "<< G4endl
357  << "* Eta in [ "<< etaMin << " - " << etaMax << "] "<< G4endl
358  << "* *" << G4endl
359  << "* *" << G4endl
360  << "**********************************************************************" << G4endl;
361  } else if (generatorInput == singleScan) {
362  G4cout << G4endl
363  << "**********************************************************************" << G4endl
364  << "* *" << G4endl
365  << "* CCalPrimaryGeneratorAction DEFAULT Scan Direction settings : *" << G4endl
366  << "* *" << G4endl
367  << "* *" << G4endl
368  << "* Phi angle in [ " << phiMin/degree << " - " << phiMax/degree << "] (deg) " << G4endl
369  << "* Eta in [ " << etaMin << " - " << etaMax << "] " << G4endl
370  << "* Steps along eta " << etaSteps << " and along phi " << phiSteps << G4endl
371  << "* *" << G4endl
372  << "* *" << G4endl
373  << "**********************************************************************" << G4endl;
374  } else if (generatorInput == singleFixed) {
375  G4cout << G4endl
376  << "*******************************************************************" << G4endl
377  << "* *" << G4endl
378  << "* CCalPrimaryGeneratorAction: Current settings : *" << G4endl
379  << "* *" << G4endl
380  << "* " << particleGun->GetNumberOfParticles()
381  << " " << particleGun->GetParticleDefinition()->GetParticleName()
382  << " of " << particleGun->GetParticleEnergy()/GeV << " GeV" << G4endl
383  << "* will be shot from " << particleGun->GetParticlePosition() << G4endl;
384  G4cout << "* in direction " << particleGun->GetParticleMomentumDirection() << G4endl;
385  G4cout << "* *" << G4endl
386  << "* *" << G4endl
387  << "*******************************************************************" << G4endl;
388  }
389  }
390 }