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