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