TGeant4/TG4PhysicsManager.cxx

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