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