source: trunk/examples/advanced/composite_calorimeter/src/CCalPrimaryGeneratorAction.cc@ 1218

Last change on this file since 1218 was 807, checked in by garnier, 17 years ago

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