TGeant4/TG4PhysicsManager.cxx

   1 // $Id$
   2 // Category: physics
   3 //
   4 // See the class description in the header file.
   5
   6 #include "TG4PhysicsManager.h"
   7 #include "TG4ParticlesManager.h"
   8 #include "TG4G3PhysicsManager.h"
   9 #include "TG4PhysicsConstructorEM.h"
  10 #include "TG4PhysicsConstructorOptical.h"
  11 #include "TG4PhysicsConstructorHadron.h"
  12 #include "TG4PhysicsConstructorSpecialCuts.h"
  13 #include "TG4PhysicsConstructorSpecialControls.h"
  14 #include "TG4G3Cut.h"
  15 #include "TG4G3Control.h"
  16
  17 #include <G4ParticleDefinition.hh>
  18 #include <G4VProcess.hh>
  19 #include <G4VModularPhysicsList.hh>
  20
  21 #include <TDatabasePDG.h>
  22
  23 TG4PhysicsManager* TG4PhysicsManager::fgInstance = 0;
  24
  25 TG4PhysicsManager::TG4PhysicsManager(G4VModularPhysicsList* physicsList)
  26   : fPhysicsList(physicsList),
  27     fSetEMPhysics(true),
  28     fSetOpticalPhysics(false),
  29     fSetHadronPhysics(false),
  30     fSetSpecialCutsPhysics(false),
  31     fSetSpecialControlsPhysics(false)
  32
  33 {
  34 //
  35   if (fgInstance) {
  36     TG4Globals::Exception(
  37       "TG4PhysicsManager: attempt to create two instances of singleton.");
  38   }
  39
  40   fgInstance = this;
  41
  42   // create particles manager
  43   fParticlesManager = new TG4ParticlesManager();
  44
  45   // create G3 physics manager
  46   fG3PhysicsManager = new TG4G3PhysicsManager();
  47
  48   // fill process name map
  49   FillProcessMap();
  50 }
  51
  52 TG4PhysicsManager::TG4PhysicsManager(){
  53 //
  54   delete fParticlesManager;
  55   delete fG3PhysicsManager;
  56 }
  57
  58 TG4PhysicsManager::TG4PhysicsManager(const TG4PhysicsManager& right) {
  59 //
  60   TG4Globals::Exception(
  61     "Attempt to copy TG4PhysicsManager singleton.");
  62 }
  63
  64 TG4PhysicsManager::~TG4PhysicsManager() {
  65 //
  66 }
  67
  68 // operators
  69
  70 TG4PhysicsManager&
  71 TG4PhysicsManager::operator=(const TG4PhysicsManager& right)
  72 {
  73   // check assignement to self
  74   if (this == &right) return *this;
  75
  76   TG4Globals::Exception(
  77     "Attempt to assign TG4PhysicsManager singleton.");
  78
  79   return *this;
  80 }
  81
  82 // private methods
  83
  84 void TG4PhysicsManager::FillProcessMap()
  85 {
  86 // Fills fProcessMap.
  87 // The default G4 process names are used in the map.
  88 // ---
  89
  90   // multiple scattering
  91   fProcessMap.Add("msc",  kPMultipleScattering);
  92   fProcessMap.Add("Imsc", kPMultipleScattering);
  93
  94   // continuous energy loss
  95   // !! including delta rays
  96   fProcessMap.Add("eIoni",  kPEnergyLoss);
  97   fProcessMap.Add("IeIoni", kPEnergyLoss);
  98   fProcessMap.Add("LowEnergyIoni", kPEnergyLoss);
  99   fProcessMap.Add("hIoni",  kPEnergyLoss);
 100   fProcessMap.Add("IhIoni", kPEnergyLoss);
 101   fProcessMap.Add("hLowEIoni", kPEnergyLoss);
 102   fProcessMap.Add("MuIoni", kPEnergyLoss);
 103   fProcessMap.Add("IMuIonisation", kPEnergyLoss);
 104   fProcessMap.Add("ionIoni",  kPEnergyLoss);
 105   fProcessMap.Add("ionLowEIoni",  kPEnergyLoss);
 106   fProcessMap.Add("PAIonisation",  kPEnergyLoss);
 107
 108   // bending in mag. field
 109   // kPMagneticFieldL
 110
 111   // particle decay
 112   fProcessMap.Add("Decay", kPDecay);
 113
 114   // photon pair production or
 115   // muon direct pair production
 116   fProcessMap.Add("conv", kPPair);
 117   fProcessMap.Add("LowEnConversion", kPPair);
 118   fProcessMap.Add("MuPairProd", kPPair);
 119   fProcessMap.Add("IMuPairProduction", kPPair);
 120
 121   // Compton scattering
 122   fProcessMap.Add("compt", kPCompton);
 123   fProcessMap.Add("LowEnCompton", kPCompton);
 124   fProcessMap.Add("polarCompt", kPCompton);
 125
 126   // photoelectric effect
 127   fProcessMap.Add("phot", kPPhotoelectric);
 128   fProcessMap.Add("LowEnPhotoElec", kPPhotoelectric);
 129
 130   // bremsstrahlung
 131   fProcessMap.Add("eBrem", kPBrem);
 132   fProcessMap.Add("IeBrem", kPBrem);
 133   fProcessMap.Add("MuBrem", kPBrem);
 134   fProcessMap.Add("IMuBremsstrahlung", kPBrem);
 135   fProcessMap.Add("LowEnBrem", kPBrem);
 136
 137   // delta-ray production
 138   // kPDeltaRay
 139   // has to be distinguished from kPEnergyLoss on flight
 140
 141   // positron annihilation
 142   fProcessMap.Add("annihil", kPAnnihilation);
 143   fProcessMap.Add("Iannihil", kPAnnihilation);
 144
 145   // hadronic interaction
 146   // kPHadronic
 147
 148   // nuclear evaporation
 149   // kPEvaporation
 150
 151   // nuclear fission
 152   // kPNuclearFission
 153
 154   // nuclear absorption
 155   fProcessMap.Add("PionMinusAbsorptionAtRest", kPNuclearAbsorption);
 156   fProcessMap.Add("PiMinusAbsorptionAtRest", kPNuclearAbsorption);
 157   fProcessMap.Add("KaonMinusAbsorption", kPNuclearAbsorption);
 158   fProcessMap.Add("KaonMinusAbsorptionAtRest", kPNuclearAbsorption);
 159
 160   // antiproton annihilation
 161   fProcessMap.Add("AntiProtonAnnihilationAtRest", kPPbarAnnihilation);
 162   // fProcessMap.Add("AntiNeutronAnnihilationAtRest", not defined);
 163
 164   // neutron capture
 165   fProcessMap.Add("NeutronCaptureAtRest", kPNCapture);
 166   // fProcessMap.Add("LCapture", hadron capture not defined);
 167
 168   // hadronic elastic incoherent scattering
 169   fProcessMap.Add("LElastic", kPHElastic);
 170
 171   // hadronic inelastic scattering
 172   fProcessMap.Add("inelastic", kPHInhelastic);
 173
 174   // muon nuclear interaction
 175   fProcessMap.Add("MuNucl", kPMuonNuclear);
 176
 177   // exceeded time of flight cut
 178   // kPTOFlimit
 179
 180   // nuclear photofission
 181   // kPPhotoFission
 182
 183   // Rayleigh scattering
 184   fProcessMap.Add("Rayleigh Scattering", kPRayleigh);
 185
 186   // no mechanism is active, usually at the entrance of a new volume
 187   fProcessMap.Add("Transportation", kPNull);
 188
 189   // particle has fallen below energy threshold and tracking stops
 190   // kPStop
 191
 192   // Cerenkov photon absorption
 193   fProcessMap.Add("Absorption", kPLightAbsorption);
 194
 195   // Cerenkov photon reflection/refraction
 196   // kPLightScattering, kPLightReflection, kPLightRefraction
 197   // has to be inquired from the G4OpBoundary process
 198
 199   // synchrotron radiation
 200   fProcessMap.Add("SynchrotronRadiation", kPSynchrotron);
 201 }
 202
 203
 204 // public methods
 205
 206 void TG4PhysicsManager::BuildPhysics()
 207 {
 208 // Empty function - not needed in G4.
 209 // (Physics is built within /run/initialize.)
 210 // ---
 211
 212   TG4Globals::Warning(
 213     "TG4PhysicsManager::BuildPhysics: is empty function in G4 MC.");
 214 }
 215
 216 void TG4PhysicsManager::CreatePhysicsConstructors()
 217 {
 218 // Creates the selected physics constructors
 219 // and registeres them in the modular physics list.
 220 // ---
 221
 222   // electromagnetic physics
 223   if (fSetEMPhysics)
 224     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorEM());
 225
 226   // optical physics
 227   if (fSetOpticalPhysics)
 228     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorOptical());
 229
 230   // hadron physics
 231   if (fSetHadronPhysics)
 232     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorHadron());
 233
 234   if (fSetSpecialCutsPhysics)
 235     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorSpecialCuts());
 236
 237   if (fSetSpecialControlsPhysics)
 238     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorSpecialControls());
 239
 240          // all created physics constructors are deleted
 241          // in the G4VModularPhysicsList destructor
 242 }
 243
 244 void TG4PhysicsManager::SetCut(const char* cutName, Float_t cutValue)
 245 {
 246 // Sets the specified cut.
 247 // ---
 248
 249   fG3PhysicsManager->CheckLock();
 250   TG4G3Cut g3Cut = fG3PhysicsManager->GetG3Cut(cutName);
 251   if (g3Cut != kNoG3Cuts)
 252     fG3PhysicsManager->SetCut(g3Cut, cutValue);
 253   else {
 254     G4String text = "TG4PhysicsManager::SetCut:\n";
 255     text = text + "    Parameter " + cutName;
 256     text = text + " is not implemented.";
 257     TG4Globals::Warning(text);
 258   }
 259 }
 260
 261 void TG4PhysicsManager::SetProcess(const char* controlName, Int_t controlValue)
 262 {
 263 // Sets the specified process control.
 264 // ---
 265
 266   fG3PhysicsManager->CheckLock();
 267   TG4G3Control control = fG3PhysicsManager->GetG3Control(controlName);
 268   if (control != kNoG3Controls)
 269     fG3PhysicsManager->SetProcess(control, controlValue);
 270   else {
 271     G4String text = "TG4PhysicsManager::SetProcess:\n";
 272     text = text + "    Parameter " + controlName;
 273     text = text + " is not implemented.";
 274     TG4Globals::Warning(text);
 275   }
 276 }
 277
 278 Float_t TG4PhysicsManager::Xsec(char* ch, Float_t p1, Int_t i1, Int_t i2)
 279 {
 280 // Not yet implemented -> gives exception.
 281 // ---
 282
 283   TG4Globals::Exception(
 284     "TG4PhysicsManager::Xsec: not yet implemented.");
 285
 286   return 0.;
 287 }
 288
 289 Int_t TG4PhysicsManager::IdFromPDG(Int_t pdgID) const
 290 {
 291 // G4 does not use the integer particle identifiers
 292 // Id <-> PDG is identity.
 293 // ---
 294
 295   return pdgID;
 296 }
 297
 298 Int_t TG4PhysicsManager::PDGFromId(Int_t mcID) const
 299 {
 300 // G4 does not use integer particle identifiers
 301 // Id <-> PDG is identity.
 302 // ---
 303
 304   return mcID;
 305 }
 306
 307 void  TG4PhysicsManager::DefineParticles()
 308 {
 309   // ======
 310   // Taken from TGeant3
 311   //
 312   // Use ENDF-6 mapping for ions = 10000*z+10*a+iso
 313   // and add 1 000 000
 314   // and numbers above 5 000 000 for special applications
 315   //
 316
 317   const Int_t kion=10000000;
 318   const Int_t kspe=50000000;
 319
 320   const Double_t kGeV=0.9314943228;
 321   const Double_t kHslash = 1.0545726663e-27;
 322   const Double_t kErgGeV = 1/1.6021773349e-3;
 323   const Double_t kHshGeV = kHslash*kErgGeV;
 324   const Double_t kYearsToSec = 3600*24*365.25;
 325
 326   TDatabasePDG *pdgDB = TDatabasePDG::Instance();
 327
 328   pdgDB->AddParticle("Deuteron","Deuteron",2*kGeV+8.071e-3,kTRUE,
 329                      0,1,"Ion",kion+10020);
 330
 331   pdgDB->AddParticle("Triton","Triton",3*kGeV+14.931e-3,kFALSE,
 332                      kHshGeV/(12.33*kYearsToSec),1,"Ion",kion+10030);
 333
 334   pdgDB->AddParticle("Alpha","Alpha",4*kGeV+2.424e-3,kTRUE,
 335                      kHshGeV/(12.33*kYearsToSec),2,"Ion",kion+20040);
 336
 337   pdgDB->AddParticle("HE3","HE3",3*kGeV+14.931e-3,kFALSE,
 338                      0,2,"Ion",kion+20030);
 339
 340   pdgDB->AddParticle("Cherenkov","Cherenkov",0,kFALSE,
 341                      0,0,"Special",kspe+50);
 342
 343   pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,
 344                      0,0,"Special",kspe+51);
 345
 346
 347   // To do: define the PDG database extension
 348   // in a common part.
 349   //
 350   // AliMC::ExtendPDGDatabase();
 351   //
 352   // end of "common" implementation
 353   // ======
 354
 355   fParticlesManager->MapParticles();
 356 }
 357
 358
 359 void TG4PhysicsManager::SetProcessActivation()
 360 {
 361 // (In)Activates built processes according
 362 // to the setup in fControlVector.
 363 // ---
 364
 365   if (fPhysicsList) {
 366     // temporarily excluded
 367     // fPhysicsList->SetProcessActivation();
 368   }
 369   else {
 370     G4String text = "TG4PhysicsManager::SetProcessActivation:\n";
 371     text = text +   "   There is no physics list set.";
 372     TG4Globals::Exception(text);
 373   }
 374 }
 375
 376
 377 AliMCProcess TG4PhysicsManager::GetMCProcess(const G4VProcess* process)
 378 {
 379 // Returns the AliMCProcess code of the specified G4 process.
 380 // ---
 381
 382   if (!process) return kPNoProcess;
 383
 384   G4String name = process->GetProcessName();
 385   G4int code = fProcessMap.GetSecond(name);
 386
 387   if (code == 0) return kPNoProcess;
 388
 389   return (AliMCProcess)code;
 390 }
 391