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