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 "TG4GeometryServices.h"
  15 #include "TG4G3Cut.h"
  16 #include "TG4G3Control.h"
  17 #include "TG4Limits.h"
  18 #include "AliDecayer.h"
  19
  20 #include <G4ParticleDefinition.hh>
  21 #include <G4VProcess.hh>
  22 #include <G4VModularPhysicsList.hh>
  23 #include <G3MedTable.hh>
  24
  25 #include <TDatabasePDG.h>
  26
  27 TG4PhysicsManager* TG4PhysicsManager::fgInstance = 0;
  28
  29 //_____________________________________________________________________________
  30 TG4PhysicsManager::TG4PhysicsManager(G4VModularPhysicsList* physicsList)
  31   : fPhysicsList(physicsList),
  32     fDecayer(0),
  33     fSetEMPhysics(true),
  34     fSetOpticalPhysics(false),
  35     fSetHadronPhysics(false),
  36     fSetSpecialCutsPhysics(false),
  37     fSetSpecialControlsPhysics(false)
  38
  39 {
  40 //
  41   if (fgInstance) {
  42     TG4Globals::Exception(
  43       "TG4PhysicsManager: attempt to create two instances of singleton.");
  44   }
  45
  46   fgInstance = this;
  47
  48   // create particles manager
  49   fParticlesManager = new TG4ParticlesManager();
  50
  51   // create G3 physics manager
  52   fG3PhysicsManager = new TG4G3PhysicsManager();
  53
  54   // fill process name map
  55   FillProcessMap();
  56 }
  57
  58 //_____________________________________________________________________________
  59 TG4PhysicsManager::TG4PhysicsManager(){
  60 //
  61   delete fDecayer;
  62   delete fParticlesManager;
  63   delete fG3PhysicsManager;
  64 }
  65
  66 //_____________________________________________________________________________
  67 TG4PhysicsManager::TG4PhysicsManager(const TG4PhysicsManager& right) {
  68 //
  69   TG4Globals::Exception(
  70     "Attempt to copy TG4PhysicsManager singleton.");
  71 }
  72
  73 //_____________________________________________________________________________
  74 TG4PhysicsManager::~TG4PhysicsManager() {
  75 //
  76 }
  77
  78 // operators
  79
  80 //_____________________________________________________________________________
  81 TG4PhysicsManager&
  82 TG4PhysicsManager::operator=(const TG4PhysicsManager& right)
  83 {
  84   // check assignement to self
  85   if (this == &right) return *this;
  86
  87   TG4Globals::Exception(
  88     "Attempt to assign TG4PhysicsManager singleton.");
  89
  90   return *this;
  91 }
  92
  93 // private methods
  94
  95 //_____________________________________________________________________________
  96 void TG4PhysicsManager::FillProcessMap()
  97 {
  98 // Fills fProcessMap.
  99 // The default G4 process names are used in the map.
 100 // ---
 101
 102   // multiple scattering
 103   fProcessMap.Add("msc",  kPMultipleScattering);
 104   fProcessMap.Add("Imsc", kPMultipleScattering);
 105
 106   // continuous energy loss
 107   // !! including delta rays
 108   fProcessMap.Add("eIoni",  kPEnergyLoss);
 109   fProcessMap.Add("IeIoni", kPEnergyLoss);
 110   fProcessMap.Add("LowEnergyIoni", kPEnergyLoss);
 111   fProcessMap.Add("hIoni",  kPEnergyLoss);
 112   fProcessMap.Add("IhIoni", kPEnergyLoss);
 113   fProcessMap.Add("hLowEIoni", kPEnergyLoss);
 114   fProcessMap.Add("MuIoni", kPEnergyLoss);
 115   fProcessMap.Add("IMuIonisation", kPEnergyLoss);
 116   fProcessMap.Add("ionIoni",  kPEnergyLoss);
 117   fProcessMap.Add("ionLowEIoni",  kPEnergyLoss);
 118   fProcessMap.Add("PAIonisation",  kPEnergyLoss);
 119
 120   // bending in mag. field
 121   // kPMagneticFieldL
 122
 123   // particle decay
 124   fProcessMap.Add("Decay", kPDecay);
 125
 126   // photon pair production or
 127   // muon direct pair production
 128   fProcessMap.Add("conv", kPPair);
 129   fProcessMap.Add("LowEnConversion", kPPair);
 130   fProcessMap.Add("MuPairProd", kPPair);
 131   fProcessMap.Add("IMuPairProduction", kPPair);
 132
 133   // Compton scattering
 134   fProcessMap.Add("compt", kPCompton);
 135   fProcessMap.Add("LowEnCompton", kPCompton);
 136   fProcessMap.Add("polarCompt", kPCompton);
 137
 138   // photoelectric effect
 139   fProcessMap.Add("phot", kPPhotoelectric);
 140   fProcessMap.Add("LowEnPhotoElec", kPPhotoelectric);
 141
 142   // bremsstrahlung
 143   fProcessMap.Add("eBrem", kPBrem);
 144   fProcessMap.Add("IeBrem", kPBrem);
 145   fProcessMap.Add("MuBrems", kPBrem);
 146   fProcessMap.Add("IMuBremsstrahlung", kPBrem);
 147   fProcessMap.Add("LowEnBrem", kPBrem);
 148
 149   // delta-ray production
 150   // kPDeltaRay
 151   // has to be distinguished from kPEnergyLoss on flight
 152
 153   // positron annihilation
 154   fProcessMap.Add("annihil", kPAnnihilation);
 155   fProcessMap.Add("Iannihil", kPAnnihilation);
 156
 157   // hadronic interaction
 158   // kPHadronic
 159
 160   // nuclear evaporation
 161   // kPEvaporation
 162
 163   // nuclear fission
 164   // kPNuclearFission
 165
 166   // nuclear absorption
 167   fProcessMap.Add("PionMinusAbsorptionAtRest", kPNuclearAbsorption);
 168   fProcessMap.Add("PiMinusAbsorptionAtRest", kPNuclearAbsorption);
 169   fProcessMap.Add("KaonMinusAbsorption", kPNuclearAbsorption);
 170   fProcessMap.Add("KaonMinusAbsorptionAtRest", kPNuclearAbsorption);
 171
 172   // antiproton annihilation
 173   fProcessMap.Add("AntiProtonAnnihilationAtRest", kPPbarAnnihilation);
 174   // fProcessMap.Add("AntiNeutronAnnihilationAtRest", not defined);
 175
 176   // neutron capture
 177   fProcessMap.Add("NeutronCaptureAtRest", kPNCapture);
 178   // fProcessMap.Add("LCapture", hadron capture not defined);
 179
 180   // hadronic elastic incoherent scattering
 181   fProcessMap.Add("LElastic", kPHElastic);
 182
 183   // hadronic inelastic scattering
 184   fProcessMap.Add("inelastic", kPHInhelastic);
 185
 186   // muon nuclear interaction
 187   fProcessMap.Add("MuNucl", kPMuonNuclear);
 188
 189   // exceeded time of flight cut
 190   // kPTOFlimit
 191
 192   // nuclear photofission
 193   // kPPhotoFission
 194
 195   // Rayleigh scattering
 196   fProcessMap.Add("Rayleigh Scattering", kPRayleigh);
 197
 198   // no mechanism is active, usually at the entrance of a new volume
 199   fProcessMap.Add("Transportation", kPNull);
 200
 201   // particle has fallen below energy threshold and tracking stops
 202   // kPStop
 203
 204   // Cerenkov photon absorption
 205   fProcessMap.Add("Absorption", kPLightAbsorption);
 206
 207   // Cerenkov photon reflection/refraction
 208   // kPLightScattering, kPLightReflection, kPLightRefraction
 209   // has to be inquired from the G4OpBoundary process
 210
 211   // synchrotron radiation
 212   fProcessMap.Add("SynchrotronRadiation", kPSynchrotron);
 213 }
 214
 215 //_____________________________________________________________________________
 216 void TG4PhysicsManager::GstparCut(G4int itmed, TG4G3Cut par, G4double parval)
 217 {
 218 // Sets special tracking medium parameter.
 219 // It is applied to all logical volumes that use the specified
 220 // tracking medium.
 221 // ---
 222
 223   // get medium from table
 224   G3MedTableEntry* medium = G3Med.get(itmed);
 225   if (!medium) {
 226     G4String text = "TG4PhysicsManager::GstparCut: \n";
 227     text = text + "    Medium not found.";
 228     G4Exception(text);
 229   }
 230
 231   // get/create user limits
 232   G4UserLimits* limits = medium->GetLimits();
 233   TG4Limits* tg4Limits;
 234   if (limits) {
 235     tg4Limits = dynamic_cast<TG4Limits*> (limits);
 236     if (!tg4Limits)
 237       G4Exception("TG4PhysicsManager::GstparCut: Wrong limits type.");
 238   }
 239   else {
 240     tg4Limits = new TG4Limits();
 241     medium->SetLimits(tg4Limits);
 242
 243     // add verbose
 244     G4cout << "TG4PhysicsManager::GstparCut: new TG4Limits() for medium "
 245            << itmed << " has been created." << G4endl;
 246   }
 247   // set parameter
 248   tg4Limits->SetG3Cut(par, parval*GeV);
 249 }
 250
 251
 252 //_____________________________________________________________________________
 253 void TG4PhysicsManager::GstparControl(G4int itmed, TG4G3Control par,
 254                                       G4double parval)
 255 {
 256 // Sets special tracking medium parameter.
 257 // It is applied to all logical volumes that use the specified
 258 // tracking medium.
 259 // ---
 260
 261   // get medium from table
 262   G3MedTableEntry* medium = G3Med.get(itmed);
 263   if (!medium) {
 264     G4String text = "TG4PhysicsManager::GstparControl: \n";
 265     text = text + "    Medium not found.";
 266     G4Exception(text);
 267   }
 268
 269   // get/create user limits
 270   G4UserLimits* limits = medium->GetLimits();
 271   TG4Limits* tg4Limits;
 272   if (limits) {
 273     tg4Limits = dynamic_cast<TG4Limits*> (limits);
 274     if (!tg4Limits)
 275       G4Exception("TG4PhysicsManager::GstparControl: Wrong limits type.");
 276   }
 277   else {
 278     tg4Limits = new TG4Limits();
 279     medium->SetLimits(tg4Limits);
 280
 281     // add verbose
 282     G4cout << "TG4PhysicsManager::GstparControl: new TG4Limits() for medium"
 283            << itmed << " has been created." << G4endl;
 284   }
 285   // set parameter
 286   tg4Limits->SetG3Control(par, parval);
 287 }
 288
 289 // public methods
 290
 291 //_____________________________________________________________________________
 292 void TG4PhysicsManager::BuildPhysics()
 293 {
 294 // Empty function - not needed in G4.
 295 // (Physics is built within /run/initialize.)
 296 // ---
 297
 298   TG4Globals::Warning(
 299     "TG4PhysicsManager::BuildPhysics: is empty function in G4 MC.");
 300 }
 301
 302 //_____________________________________________________________________________
 303 void  TG4PhysicsManager::Gstpar(Int_t itmed, const char *param,
 304            Float_t parval)
 305 {
 306 // Passes the tracking medium parameter to TG4Limits.
 307 // The tracking medium parameter is set only in case
 308 // its value is different from the "global" physics setup.
 309 // (If: CheckCut/ControlWithG3Defaults is used checking
 310 //  is performed with respect to G3 default values.)
 311 // When any cut/control parameter is set in limits
 312 // the physics manager is locked and the physics setup
 313 // cannot be changed.
 314 // Applying this TG4Limits to particles and physical
 315 // processes is still in development.
 316 //
 317 //  Geant3 desription:
 318 //  ==================
 319 //  To change the value of cut  or mechanism "CHPAR"
 320 //  to a new value PARVAL  for tracking medium ITMED
 321 //  The  data   structure  JTMED   contains  the   standard  tracking
 322 //  parameters (CUTS and flags to control the physics processes)  which
 323 //  are used  by default  for all  tracking media.   It is  possible to
 324 //  redefine individually  with GSTPAR  any of  these parameters  for a
 325 //  given tracking medium.
 326 //  ITMED     tracking medium number
 327 //  CHPAR     is a character string (variable name)
 328 //  PARVAL    must be given as a floating point.
 329 // ---
 330
 331   G4String name = TG4GeometryServices::Instance()->CutName(param);
 332   TG4G3Cut cut;
 333   if (fG3PhysicsManager->CheckCutWithTheVector(name, parval, cut)) {
 334       GstparCut(itmed, cut, parval);
 335       fG3PhysicsManager->Lock();
 336   }
 337   else {
 338     TG4G3Control control;
 339     if (fG3PhysicsManager->CheckControlWithTheVector(name, parval, control)) {
 340       GstparControl(itmed, control, parval);
 341       fG3PhysicsManager->Lock();
 342     }
 343     else if (cut==kNoG3Cuts && control==kNoG3Controls) {
 344       G4String text = "TG4PhysicsManager::Gstpar:";
 345       text = text + name;
 346       text = text + " parameter is not yet implemented.";
 347       TG4Globals::Warning(text);
 348     }
 349   }
 350 }
 351
 352 //_____________________________________________________________________________
 353 void TG4PhysicsManager::CreatePhysicsConstructors()
 354 {
 355 // Creates the selected physics constructors
 356 // and registeres them in the modular physics list.
 357 // ---
 358
 359   // electromagnetic physics
 360   if (fSetEMPhysics)
 361     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorEM());
 362
 363   // optical physics
 364   if (fSetOpticalPhysics)
 365     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorOptical());
 366
 367   // hadron physics
 368   if (fSetHadronPhysics)
 369     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorHadron());
 370
 371   if (fSetSpecialCutsPhysics)
 372     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorSpecialCuts());
 373
 374   if (fSetSpecialControlsPhysics)
 375     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorSpecialControls());
 376
 377          // all created physics constructors are deleted
 378          // in the G4VModularPhysicsList destructor
 379 }
 380
 381 //_____________________________________________________________________________
 382 void TG4PhysicsManager::SetCut(const char* cutName, Float_t cutValue)
 383 {
 384 // Sets the specified cut.
 385 // ---
 386
 387   fG3PhysicsManager->CheckLock();
 388   TG4G3Cut g3Cut = fG3PhysicsManager->GetG3Cut(cutName);
 389   if (g3Cut != kNoG3Cuts)
 390     fG3PhysicsManager->SetCut(g3Cut, cutValue);
 391   else {
 392     G4String text = "TG4PhysicsManager::SetCut:\n";
 393     text = text + "    Parameter " + cutName;
 394     text = text + " is not implemented.";
 395     TG4Globals::Warning(text);
 396   }
 397 }
 398
 399 //_____________________________________________________________________________
 400 void TG4PhysicsManager::SetProcess(const char* controlName, Int_t controlValue)
 401 {
 402 // Sets the specified process control.
 403 // ---
 404
 405   fG3PhysicsManager->CheckLock();
 406   TG4G3Control control = fG3PhysicsManager->GetG3Control(controlName);
 407   if (control != kNoG3Controls)
 408     fG3PhysicsManager->SetProcess(control, controlValue);
 409   else {
 410     G4String text = "TG4PhysicsManager::SetProcess:\n";
 411     text = text + "    Parameter " + controlName;
 412     text = text + " is not implemented.";
 413     TG4Globals::Warning(text);
 414   }
 415 }
 416
 417 //_____________________________________________________________________________
 418 Float_t TG4PhysicsManager::Xsec(char* ch, Float_t p1, Int_t i1, Int_t i2)
 419 {
 420 // Not yet implemented -> gives exception.
 421 // ---
 422
 423   TG4Globals::Exception(
 424     "TG4PhysicsManager::Xsec: not yet implemented.");
 425
 426   return 0.;
 427 }
 428
 429 //_____________________________________________________________________________
 430 Int_t TG4PhysicsManager::IdFromPDG(Int_t pdgID) const
 431 {
 432 // G4 does not use the integer particle identifiers
 433 // Id <-> PDG is identity.
 434 // ---
 435
 436   return pdgID;
 437 }
 438
 439 //_____________________________________________________________________________
 440 Int_t TG4PhysicsManager::PDGFromId(Int_t mcID) const
 441 {
 442 // G4 does not use integer particle identifiers
 443 // Id <-> PDG is identity.
 444 // ---
 445
 446   return mcID;
 447 }
 448
 449 //_____________________________________________________________________________
 450 void  TG4PhysicsManager::DefineParticles()
 451 {
 452   // ======
 453   // Taken from TGeant3
 454   //
 455   // Use ENDF-6 mapping for ions = 10000*z+10*a+iso
 456   // and add 1 000 000
 457   // and numbers above 5 000 000 for special applications
 458   //
 459
 460   const Int_t kion=10000000;
 461   const Int_t kspe=50000000;
 462
 463   const Double_t kGeV=0.9314943228;
 464   const Double_t kHslash = 1.0545726663e-27;
 465   const Double_t kErgGeV = 1/1.6021773349e-3;
 466   const Double_t kHshGeV = kHslash*kErgGeV;
 467   const Double_t kYearsToSec = 3600*24*365.25;
 468
 469   TDatabasePDG *pdgDB = TDatabasePDG::Instance();
 470
 471   pdgDB->AddParticle("Deuteron","Deuteron",2*kGeV+8.071e-3,kTRUE,
 472                      0,1,"Ion",kion+10020);
 473
 474   pdgDB->AddParticle("Triton","Triton",3*kGeV+14.931e-3,kFALSE,
 475                      kHshGeV/(12.33*kYearsToSec),1,"Ion",kion+10030);
 476
 477   pdgDB->AddParticle("Alpha","Alpha",4*kGeV+2.424e-3,kTRUE,
 478                      kHshGeV/(12.33*kYearsToSec),2,"Ion",kion+20040);
 479
 480   pdgDB->AddParticle("HE3","HE3",3*kGeV+14.931e-3,kFALSE,
 481                      0,2,"Ion",kion+20030);
 482
 483   pdgDB->AddParticle("Cherenkov","Cherenkov",0,kFALSE,
 484                      0,0,"Special",kspe+50);
 485
 486   pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,
 487                      0,0,"Special",kspe+51);
 488
 489
 490   // To do: define the PDG database extension
 491   // in a common part.
 492   //
 493   // AliMC::ExtendPDGDatabase();
 494   //
 495   // end of "common" implementation
 496   // ======
 497
 498   fParticlesManager->MapParticles();
 499 }
 500
 501 //_____________________________________________________________________________
 502 void TG4PhysicsManager::SetProcessActivation()
 503 {
 504 // (In)Activates built processes according
 505 // to the setup in fControlVector.
 506 // ---
 507
 508   if (fPhysicsList) {
 509     // temporarily excluded
 510     // fPhysicsList->SetProcessActivation();
 511   }
 512   else {
 513     G4String text = "TG4PhysicsManager::SetProcessActivation:\n";
 514     text = text +   "   There is no physics list set.";
 515     TG4Globals::Exception(text);
 516   }
 517 }
 518
 519
 520 //_____________________________________________________________________________
 521 AliMCProcess TG4PhysicsManager::GetMCProcess(const G4VProcess* process)
 522 {
 523 // Returns the AliMCProcess code of the specified G4 process.
 524 // ---
 525
 526   if (!process) return kPNoProcess;
 527
 528   G4String name = process->GetProcessName();
 529   G4int code = fProcessMap.GetSecond(name);
 530
 531   if (code == 0) return kPNoProcess;
 532
 533   return (AliMCProcess)code;
 534 }
 535