Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 "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 
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:
60  SetMinimumEta(0.);
61  SetMaximumEta(3.5);
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
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 
83  if (gunMessenger)
84  delete gunMessenger;
85  if (particleGun)
86  delete particleGun;
87 }
88 
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 #ifdef debug
107  if (verboseLevel >= 2 ) {
108  G4cout << "Energy " << particleEnergy/GeV << " GeV; Theta "
109  << theta/deg << " degree; Phi " << phi/deg << " degree" << G4endl;
110  G4cout << "Shooting in " << particleDir << " direction "<< G4endl;
111  }
112 #endif
113  } else if (generatorInput == singleScan) {
114  G4double scanEtaStep, scanPhiStep;
115  if (scanSteps == 0) {
116  scanPhiStep = (phiMax - phiMin) / phiSteps;
117  phiValue = phiMin - scanPhiStep; //first step is going to change scanCurrentPhiValue
118  etaValue = etaMin;
119  }
120 
121  scanEtaStep = (etaMax - etaMin) / etaSteps;
122  scanPhiStep = (phiMax - phiMin) / phiSteps;
123 #ifdef debug
124  if (verboseLevel > 2 ) {
125  G4cout << " scanEtaStep " << scanEtaStep << " # of Steps " << etaSteps
126  << G4endl;
127  G4cout << " scanPhiStep " << scanPhiStep << " # of Steps " << phiSteps
128  << G4endl;
129  }
130 #endif
131 
132  //----- First scan in phi, then in eta
133  if (phiMax - phiValue < 1.E-6 * scanPhiStep) { // !only <= 1.E6 steps allowed
134  if (etaMax - etaValue < 1.E-6 * scanEtaStep) { // !only <= 1.E6 steps allowed
135  G4cout << " Scan completed!" << G4endl;
136  return;
137  } else {
138  etaValue += scanEtaStep;
139  phiValue = phiMin;
140  }
141  } else {
142  phiValue += scanPhiStep;
143  }
144  G4double theta = std::atan(G4Exp(-etaValue))*2.;
145 
146  G4double scanX = std::sin(theta)*std::cos(phiValue);
147  G4double scanY = std::sin(theta)*std::sin(phiValue);
148  G4double scanZ = std::cos(theta);
149  if (verboseLevel >= 2 ) {
150  G4cout << "Scan eta " << etaValue << " Phi " << phiValue/deg
151  << " theta " << theta/deg << G4endl;
152  }
153  particleDir = G4ThreeVector(scanX,scanY,scanZ);
154  particleGun->SetParticleMomentumDirection(particleDir);
155 #ifdef debug
156  if (verboseLevel > 2 ) {
157  G4cout << "Shooting in " << particleDir << " direction "<< G4endl;
158  }
159 #endif
160  scanSteps++;
161  }
162 
163  // Generate GEANT4 primary vertex
164  particleGun->GeneratePrimaryVertex(anEvent);
165 }
166 
167 
169  verboseLevel = val;
170 }
171 
172 
174 
175  if (val=="on") {
176  generatorInput = singleRandom;
177  print (1);
178  } else {
179  generatorInput = singleFixed;
180  print (1);
181  }
182 }
183 
184 
186 
187  if (val=="on") {
188  generatorInput = singleScan;
189  scanSteps = 0;
190  print (1);
191  } else {
192  generatorInput = singleFixed;
193  print (1);
194  }
195 }
196 
197 
199 
200  if (p <= 0.) {
201  G4cerr << "CCalPrimaryGeneratorAction::SetMinimumEnergy: value " << p/GeV
202  << "GeV is out of bounds, it will not be used" << G4endl;
203  G4cerr << " Should be >0. Please check" << G4endl;
204  } else {
205  energyMin = p;
206 #ifdef debug
207  if (verboseLevel >= 1 ) {
208  G4cout << " CCalPrimaryGeneratorAction: setting min. value of energy to "
209  << energyMin/GeV << " GeV " << G4endl;
210  }
211 #endif
212  }
213 }
214 
215 
217 
218  if (p <= 0.) {
219  G4cerr << "CCalPrimaryGeneratorAction::SetMaximumEnergy: value " << p/GeV
220  << "GeV is out of bounds, it will not be used" << G4endl;
221  G4cerr << " Should be >0. Please check" << G4endl;
222  } else {
223  energyMax = p;
224 #ifdef debug
225  if (verboseLevel >= 1 ) {
226  G4cout << " CCalPrimaryGeneratorAction: setting max. value of energy to "
227  << energyMax/GeV << " GeV " << G4endl;
228  }
229 #endif
230  }
231 }
232 
233 
235 
236  if (std::fabs(p)>2.*pi) {
237  G4cerr << "CCalPrimaryGeneratorAction::SetMinimumPhi: setting value quite "
238  << "large" << G4endl;
239  G4cerr << " Should be given in radians - Please check" << G4endl;
240  } else {
241  phiMin = std::fabs(p);
242 #ifdef debug
243  if (verboseLevel >= 1 ) {
244  G4cout << " CCalPrimaryGeneratorAction: setting min. value of phi to "
245  << phiMin << G4endl;
246  }
247 #endif
248  }
249 }
250 
251 
253 
254  if (std::fabs(p)>2.*pi) {
255  G4cerr << "CCalPrimaryGeneratorAction::SetMaximumPhi: setting value quite "
256  << "large" << G4endl;
257  G4cerr << " Should be given in radians - Please check" << G4endl;
258  } else {
259  phiMax = std::fabs(p);
260 #ifdef debug
261  if (verboseLevel >= 1 ) {
262  G4cout << " CCalPrimaryGeneratorAction: setting max. value of phi to "
263  << phiMax << G4endl;
264  }
265 #endif
266  }
267 }
268 
269 
271 
272  if (val <= 0) {
273  G4cerr << "CCalPrimaryGeneratorAction::SetStepsPhi: value " << val
274  << " is out of bounds, it will not be used" << G4endl;
275  G4cerr << " Should be > 0 Please check" << G4endl;
276  } else {
277  phiSteps = val;
278 #ifdef debug
279  if (verboseLevel >= 1 ) {
280  G4cout << " CCalPrimaryGeneratorAction: setting no. of steps in phi to "
281  << phiSteps << G4endl;
282  }
283 #endif
284  }
285 }
286 
287 
289 
290  etaMin = p;
291 #ifdef debug
292  if (verboseLevel >= 1 ) {
293  G4cout << " CCalPrimaryGeneratorAction: setting min. value of eta to "
294  << etaMin << G4endl;
295  }
296 #endif
297 }
298 
299 
301 
302  etaMax = p;
303 #ifdef debug
304  if (verboseLevel >= 1 ) {
305  G4cout << " CCalPrimaryGeneratorAction: setting max. value of eta to "
306  << etaMax << G4endl;
307  }
308 #endif
309 }
310 
311 
313 
314  if (val <= 0) {
315  G4cerr<<"CCalPrimaryGeneratorAction::SetStepsEta: value " << val << " is out of bounds, it will not be used"<<G4endl;
316  G4cerr<<" Should be > 0 Please check"<<G4endl;
317  } else {
318  etaSteps = val;
319 #ifdef debug
320  if (verboseLevel >= 1 ) {
321  G4cout << " CCalPrimaryGeneratorAction: setting no. of steps in eta to "
322  << etaSteps << G4endl;
323  }
324 #endif
325  }
326 }
327 
329 
330  particleGun->SetParticlePosition(pos);
331 }
332 
335 }
336 
337 void CCalPrimaryGeneratorAction::initialize(){
338 
339  isInitialized = 1;
340 
341  print (1);
342 }
343 
344 
345 void CCalPrimaryGeneratorAction::print(G4int val){
346 
347 #ifdef debug
348  if (verboseLevel >= val) {
349 
350  if (generatorInput == singleRandom) {
351  G4cout << G4endl
352  << "**********************************************************************" << G4endl
353  << "* *" << G4endl
354  << "* CCalPrimaryGeneratorAction DEFAULT Random Energy/Direction setting:*" << G4endl
355  << "* *" << G4endl
356  << "* *" << G4endl
357  << "* Energy in [ "<< energyMin/GeV << " - " << energyMax/GeV << "] (GeV) "<< G4endl
358  << "* Phi angle in [ "<< phiMin << " - " << phiMax << "] (rad) "<< G4endl
359  << "* [ "<< phiMin/degree << " - " << phiMax/degree << "] (deg) "<< G4endl
360  << "* Eta in [ "<< etaMin << " - " << etaMax << "] "<< G4endl
361  << "* *" << G4endl
362  << "* *" << G4endl
363  << "**********************************************************************" << G4endl;
364  } else if (generatorInput == singleScan) {
365  G4cout << G4endl
366  << "**********************************************************************" << G4endl
367  << "* *" << G4endl
368  << "* CCalPrimaryGeneratorAction DEFAULT Scan Direction settings : *" << G4endl
369  << "* *" << G4endl
370  << "* *" << G4endl
371  << "* Phi angle in [ " << phiMin/degree << " - " << phiMax/degree << "] (deg) " << G4endl
372  << "* Eta in [ " << etaMin << " - " << etaMax << "] " << G4endl
373  << "* Steps along eta " << etaSteps << " and along phi " << phiSteps << G4endl
374  << "* *" << G4endl
375  << "* *" << G4endl
376  << "**********************************************************************" << G4endl;
377  } else if (generatorInput == singleFixed) {
378  G4cout << G4endl
379  << "*******************************************************************" << G4endl
380  << "* *" << G4endl
381  << "* CCalPrimaryGeneratorAction: Current settings : *" << G4endl
382  << "* *" << G4endl
383  << "* " << particleGun->GetNumberOfParticles()
384  << " " << particleGun->GetParticleDefinition()->GetParticleName()
385  << " of " << particleGun->GetParticleEnergy()/GeV << " GeV" << G4endl
386  << "* will be shot from " << particleGun->GetParticlePosition() << G4endl;
387  G4cout << "* in direction " << particleGun->GetParticleMomentumDirection() << G4endl;
388  G4cout << "* *" << G4endl
389  << "* *" << G4endl
390  << "*******************************************************************" << G4endl;
391  }
392  }
393 #endif
394 
395 }
396 
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
static constexpr double TeV
Definition: G4SIunits.hh:218
G4ParticleMomentum GetParticleMomentumDirection() const
static double shoot()
Definition: RandFlat.cc:59
void SetParticlePosition(G4ThreeVector aPosition)
G4GLOB_DLL std::ostream G4cout
static constexpr double degree
Definition: G4SIunits.hh:144
void SetGunPosition(const G4ThreeVector &pos) const
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4int GetNumberOfParticles() const
void SetParticleEnergy(G4double aKineticEnergy)
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:79
static G4ParticleTable * GetParticleTable()
void SetRunIDCounter(G4int i)
G4ParticleDefinition * GetParticleDefinition() const
static constexpr double GeV
Definition: G4SIunits.hh:217
#define G4endl
Definition: G4ios.hh:61
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
static constexpr double deg
Definition: G4SIunits.hh:152
G4bool isInitialized()
static const G4double pos
void SetParticleDefinition(G4ParticleDefinition *aParticleDefinition)
G4double GetParticleEnergy() const
G4GLOB_DLL std::ostream G4cerr