composite_calorimeter/src/CCalPrimaryGeneratorAction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 ///////////////////////////////////////////////////////////////////////////////
0027 // File: CCalPrimaryGeneratorAction.cc
0028 // Description: CCalPrimaryGeneratorAction Sets up particle beam
0029 ///////////////////////////////////////////////////////////////////////////////
0030
0031 #include <CLHEP/Random/RandFlat.h>
0032
0033 #include "CCalPrimaryGeneratorAction.hh"
0034 #include "CCalPrimaryGeneratorMessenger.hh"
0035 #include "G4HEPEvtInterface.hh"
0036
0037 #include "G4Exp.hh"
0038 #include "G4PhysicalConstants.hh"
0039 #include "G4SystemOfUnits.hh"
0040 #include "G4Event.hh"
0041 #include "G4ParticleGun.hh"
0042 #include "G4ParticleTable.hh"
0043 #include "G4ParticleDefinition.hh"
0044 #include "G4HEPEvtInterface.hh"
0045 #include "G4RunManager.hh"
0046
0047 //#define debug
0048
0049 CCalPrimaryGeneratorAction::CCalPrimaryGeneratorAction(): particleGun(0),
0050   generatorInput(singleFixed),  verboseLevel(0), n_particle(1),
0051   particleName("pi-"), particleEnergy(100*GeV), particlePosition(0.,0.,0.),
0052   particleDir(1.,1.,0.1), isInitialized(0), scanSteps(0) {
0053
0054   //Initialise the messenger
0055   gunMessenger = new CCalPrimaryGeneratorMessenger(this);
0056
0057   //Default settings:
0058   SetMinimumEnergy(1.*GeV);
0059   SetMaximumEnergy(1.*TeV);
0060   SetMinimumEta(0.);
0061   SetMaximumEta(3.5);
0062   SetMinimumPhi(0.*degree);
0063   SetMaximumPhi(360.*degree);
0064   SetStepsPhi(1);
0065   SetStepsEta(1);
0066
0067   // Number of particles per gun
0068   particleGun = new G4ParticleGun(n_particle);
0069
0070   // Getting particle definition
0071   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
0072   G4ParticleDefinition* particle = particleTable->FindParticle(particleName);
0073   particleGun->SetParticleDefinition(particle);
0074
0075   // Setting particle properties
0076   particleGun->SetParticleEnergy(particleEnergy);
0077   particleGun->SetParticleMomentumDirection(particleDir);
0078   particleGun->SetParticlePosition(particlePosition);
0079   print(0);
0080 }
0081
0082 CCalPrimaryGeneratorAction::~CCalPrimaryGeneratorAction() {
0083   if (gunMessenger)
0084     delete gunMessenger;
0085   if (particleGun)
0086     delete particleGun;
0087 }
0088
0089 void CCalPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) {
0090
0091   if (isInitialized == 0) initialize();
0092
0093   if (generatorInput == singleRandom) {
0094     particleEnergy = CLHEP::RandFlat::shoot(energyMin,energyMax);
0095     particleGun->SetParticleEnergy(particleEnergy);
0096
0097     G4double eta = CLHEP::RandFlat::shoot(etaMin,etaMax);
0098     G4double phi = CLHEP::RandFlat::shoot(phiMin,phiMax);
0099     G4double theta = std::atan(G4Exp(-eta))*2.;
0100     G4double randomX = std::sin(theta)*std::cos(phi);
0101     G4double randomY = std::sin(theta)*std::sin(phi);
0102     G4double randomZ = std::cos(theta);
0103
0104     particleDir = G4ThreeVector(randomX,randomY,randomZ);
0105     particleGun->SetParticleMomentumDirection(particleDir);
0106     if (verboseLevel >= 2 ) {
0107       G4cout << "Energy " << particleEnergy/GeV << " GeV; Theta "
0108              << theta/deg << " degree; Phi " << phi/deg << " degree" << G4endl;
0109       G4cout << "Shooting in " << particleDir << " direction "<< G4endl;
0110     }
0111   } else if (generatorInput == singleScan) {
0112     G4double scanEtaStep, scanPhiStep;
0113     if (scanSteps == 0) {
0114       scanPhiStep = (phiMax - phiMin) / phiSteps;
0115       phiValue = phiMin - scanPhiStep; //first step is going to change scanCurrentPhiValue
0116       etaValue = etaMin;
0117     }
0118
0119     scanEtaStep = (etaMax - etaMin) / etaSteps;
0120     scanPhiStep = (phiMax - phiMin) / phiSteps;
0121 #ifdef debug
0122     if (verboseLevel > 2 ) {
0123       G4cout << " scanEtaStep " << scanEtaStep << " # of Steps " << etaSteps
0124              << G4endl;
0125       G4cout << " scanPhiStep " << scanPhiStep << " # of Steps " << phiSteps
0126              << G4endl;
0127     }
0128 #endif
0129
0130     //----- First scan in phi, then in eta
0131     if (phiMax - phiValue < 1.E-6 * scanPhiStep) { // !only <= 1.E6 steps allowed
0132       if (etaMax - etaValue < 1.E-6 * scanEtaStep) { // !only <= 1.E6 steps allowed
0133         G4cout << " Scan completed!" << G4endl;
0134         return;
0135       } else {
0136         etaValue += scanEtaStep;
0137         phiValue  = phiMin;
0138       }
0139     } else {
0140       phiValue += scanPhiStep;
0141     }
0142     G4double theta = std::atan(G4Exp(-etaValue))*2.;
0143
0144     G4double scanX = std::sin(theta)*std::cos(phiValue);
0145     G4double scanY = std::sin(theta)*std::sin(phiValue);
0146     G4double scanZ = std::cos(theta);
0147     if (verboseLevel >= 2 ) {
0148       G4cout << "Scan eta " << etaValue << " Phi " << phiValue/deg
0149              << " theta " << theta/deg << G4endl;
0150     }
0151     particleDir = G4ThreeVector(scanX,scanY,scanZ);
0152     particleGun->SetParticleMomentumDirection(particleDir);
0153 #ifdef debug
0154     if (verboseLevel > 2 ) {
0155       G4cout  << "Shooting in " << particleDir << " direction "<< G4endl;
0156     }
0157 #endif
0158     scanSteps++;
0159   }
0160
0161   // Generate GEANT4 primary vertex
0162   particleGun->GeneratePrimaryVertex(anEvent);
0163 }
0164
0165
0166 void CCalPrimaryGeneratorAction::SetVerboseLevel(G4int val){
0167   verboseLevel = val;
0168 }
0169
0170
0171 void CCalPrimaryGeneratorAction::SetRandom(G4String val) {
0172
0173   if (val=="on") {
0174     generatorInput = singleRandom;
0175     print (1);
0176   } else {
0177     generatorInput = singleFixed;
0178     print (1);
0179   }
0180 }
0181
0182
0183 void CCalPrimaryGeneratorAction::SetScan(G4String val) {
0184
0185   if (val=="on") {
0186     generatorInput = singleScan;
0187     scanSteps = 0;
0188     print (1);
0189   } else {
0190     generatorInput = singleFixed;
0191     print (1);
0192   }
0193 }
0194
0195
0196 void CCalPrimaryGeneratorAction::SetMinimumEnergy(G4double p){
0197
0198   if (p <= 0.) {
0199     G4cerr << "CCalPrimaryGeneratorAction::SetMinimumEnergy: value " << p/GeV
0200            << "GeV is out of bounds, it will not be used" << G4endl;
0201     G4cerr << " Should be  >0.  Please check" << G4endl;
0202   } else {
0203     energyMin = p;
0204 #ifdef debug
0205     if (verboseLevel >= 1 ) {
0206       G4cout << " CCalPrimaryGeneratorAction: setting min. value of energy to "
0207              << energyMin/GeV << " GeV " << G4endl;
0208     }
0209 #endif
0210   }
0211 }
0212
0213
0214 void CCalPrimaryGeneratorAction::SetMaximumEnergy(G4double p){
0215
0216   if (p <= 0.) {
0217     G4cerr << "CCalPrimaryGeneratorAction::SetMaximumEnergy: value " << p/GeV
0218            << "GeV is out of bounds, it will not be used" << G4endl;
0219     G4cerr << " Should be  >0.  Please check" << G4endl;
0220   } else {
0221     energyMax = p;
0222 #ifdef debug
0223     if (verboseLevel >= 1 ) {
0224       G4cout << " CCalPrimaryGeneratorAction: setting max. value of energy to "
0225              << energyMax/GeV << " GeV " << G4endl;
0226     }
0227 #endif
0228   }
0229 }
0230
0231
0232 void CCalPrimaryGeneratorAction::SetMinimumPhi(G4double p){
0233
0234   if (std::fabs(p)>2.*pi) {
0235     G4cerr << "CCalPrimaryGeneratorAction::SetMinimumPhi: setting value quite "
0236            << "large" << G4endl;
0237     G4cerr << " Should be given in radians - Please check" << G4endl;
0238   } else {
0239     phiMin = std::fabs(p);
0240 #ifdef debug
0241     if (verboseLevel >= 1 ) {
0242       G4cout << " CCalPrimaryGeneratorAction: setting min. value of phi to "
0243              << phiMin << G4endl;
0244     }
0245 #endif
0246   }
0247 }
0248
0249
0250 void CCalPrimaryGeneratorAction::SetMaximumPhi(G4double p){
0251
0252   if (std::fabs(p)>2.*pi) {
0253     G4cerr << "CCalPrimaryGeneratorAction::SetMaximumPhi: setting value quite "
0254            << "large" << G4endl;
0255     G4cerr << " Should be given in radians - Please check" << G4endl;
0256   } else {
0257     phiMax = std::fabs(p);
0258 #ifdef debug
0259     if (verboseLevel >= 1 ) {
0260       G4cout << " CCalPrimaryGeneratorAction: setting max. value of phi to "
0261              << phiMax << G4endl;
0262     }
0263 #endif
0264   }
0265 }
0266
0267
0268 void CCalPrimaryGeneratorAction::SetStepsPhi(G4int val){
0269
0270   if (val <= 0) {
0271     G4cerr << "CCalPrimaryGeneratorAction::SetStepsPhi: value " << val
0272            << " is out of bounds, it will not be used" << G4endl;
0273     G4cerr << " Should be  > 0  Please check" << G4endl;
0274   } else {
0275     phiSteps = val;
0276 #ifdef debug
0277     if (verboseLevel >= 1 ) {
0278       G4cout << " CCalPrimaryGeneratorAction: setting no. of steps in phi to "
0279              << phiSteps << G4endl;
0280     }
0281 #endif
0282   }
0283 }
0284
0285
0286 void CCalPrimaryGeneratorAction::SetMinimumEta(G4double p){
0287
0288   etaMin = p;
0289 #ifdef debug
0290   if (verboseLevel >= 1 ) {
0291     G4cout << " CCalPrimaryGeneratorAction: setting min. value of eta to "
0292            << etaMin << G4endl;
0293   }
0294 #endif
0295 }
0296
0297
0298 void CCalPrimaryGeneratorAction::SetMaximumEta(G4double p){
0299
0300   etaMax = p;
0301 #ifdef debug
0302   if (verboseLevel >= 1 ) {
0303     G4cout << " CCalPrimaryGeneratorAction: setting max. value of eta to "
0304            << etaMax << G4endl;
0305   }
0306 #endif
0307 }
0308
0309
0310 void CCalPrimaryGeneratorAction::SetStepsEta(G4int val){
0311
0312   if (val <= 0) {
0313     G4cerr<<"CCalPrimaryGeneratorAction::SetStepsEta: value " << val << " is out of bounds, it will not be used"<<G4endl;
0314     G4cerr<<" Should be  > 0  Please check"<<G4endl;
0315   } else {
0316     etaSteps = val;
0317 #ifdef debug
0318     if (verboseLevel >= 1 ) {
0319       G4cout << " CCalPrimaryGeneratorAction: setting no. of steps in eta to "
0320              << etaSteps << G4endl;
0321     }
0322 #endif
0323   }
0324 }
0325
0326 void CCalPrimaryGeneratorAction::SetGunPosition(const G4ThreeVector & pos) const {
0327
0328   particleGun->SetParticlePosition(pos);
0329 }
0330
0331 void CCalPrimaryGeneratorAction::SetRunNo(G4int val){
0332   G4RunManager::GetRunManager()->SetRunIDCounter( val );
0333 }
0334
0335 void CCalPrimaryGeneratorAction::initialize(){
0336
0337   isInitialized = 1;
0338
0339   print (1);
0340 }
0341
0342
0343 void CCalPrimaryGeneratorAction::print(G4int val){
0344
0345   if (verboseLevel >= val) {
0346
0347     if (generatorInput == singleRandom) {
0348       G4cout << G4endl
0349              << "**********************************************************************" << G4endl
0350              << "*                                                                    *" << G4endl
0351              << "* CCalPrimaryGeneratorAction DEFAULT Random Energy/Direction setting:*" << G4endl
0352              << "*                                                                    *" << G4endl
0353              << "*                                                                    *" << G4endl
0354              << "*   Energy in    [ "<< energyMin/GeV   << " - " << energyMax/GeV   << "] (GeV) "<< G4endl
0355              << "*   Phi angle in [ "<< phiMin          << " - " << phiMax          << "] (rad) "<< G4endl
0356              << "*                [ "<< phiMin/degree   << " - " << phiMax/degree   << "] (deg) "<< G4endl
0357              << "*   Eta in       [ "<< etaMin          << " - " << etaMax          << "]       "<< G4endl
0358              << "*                                                                    *" << G4endl
0359              << "*                                                                    *" << G4endl
0360                    << "**********************************************************************" << G4endl;
0361     } else if (generatorInput == singleScan) {
0362       G4cout << G4endl
0363              << "**********************************************************************" << G4endl
0364              << "*                                                                    *" << G4endl
0365              << "* CCalPrimaryGeneratorAction DEFAULT Scan Direction settings :       *" << G4endl
0366              << "*                                                                    *" << G4endl
0367              << "*                                                                    *" << G4endl
0368              << "*   Phi angle in [ " << phiMin/degree   << " - " << phiMax/degree << "] (deg) " << G4endl
0369              << "*   Eta in       [ " << etaMin          << " - " << etaMax        << "]       " << G4endl
0370              << "*   Steps along eta " << etaSteps << " and along phi " << phiSteps << G4endl
0371              << "*                                                                    *" << G4endl
0372              << "*                                                                    *" << G4endl
0373                    << "**********************************************************************" << G4endl;
0374     } else  if (generatorInput == singleFixed) {
0375       G4cout << G4endl
0376              << "*******************************************************************" << G4endl
0377              << "*                                                                 *" << G4endl
0378              << "* CCalPrimaryGeneratorAction: Current settings :                  *" << G4endl
0379              << "*                                                                 *" << G4endl
0380              << "* " << particleGun->GetNumberOfParticles()
0381              << "  " << particleGun->GetParticleDefinition()->GetParticleName()
0382              << " of " << particleGun->GetParticleEnergy()/GeV << " GeV" << G4endl
0383              << "* will be shot from " << particleGun->GetParticlePosition() << G4endl;
0384       G4cout << "* in direction " << particleGun->GetParticleMomentumDirection() << G4endl;
0385       G4cout << "*                                                                 *" << G4endl
0386              << "*                                                                 *" << G4endl
0387              << "*******************************************************************" << G4endl;
0388     }
0389   }
0390 }