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

Last change on this file since 1327 was 807, checked in by garnier, 16 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.