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 "TG4PhysicsList.h"
   8 #include "TG4CutVector.h"
   9 #include "TG4FlagVector.h"
  10 #include "TG4G3Defaults.h"
  11
  12 #include <G4ParticleDefinition.hh>
  13
  14 #include <TDatabasePDG.h>
  15
  16 TG4PhysicsManager* TG4PhysicsManager::fgInstance = 0;
  17
  18 TG4PhysicsManager::TG4PhysicsManager()
  19   : fLock(false),
  20     fPhysicsList(0),
  21     fCutVector(0),
  22     fFlagVector(0)
  23 {
  24 //
  25   if (fgInstance) {
  26     TG4Globals::Exception(
  27       "TG4PhysicsManager: attempt to create two instances of singleton.");
  28   }
  29
  30   fgInstance = this;
  31
  32   // initialize fIsCutVector
  33   fIsCutVector = new TG4boolVector;
  34   G4int i;
  35   //for (i=0; i<kNofParticlesWSP; i++) fIsCutVector->insert(false);
  36   for (i=0; i<kNofParticlesWSP; i++) fIsCutVector->push_back(false);
  37
  38   // initialize fIsFlagVector
  39   fIsFlagVector = new TG4boolVector;
  40   //for (i=0; i<kNofParticlesWSP; i++) fIsFlagVector->insert(false);
  41   for (i=0; i<kNofParticlesWSP; i++) fIsFlagVector->push_back(false);
  42
  43   // define fCutNameVector, fFlagNameVector
  44   FillG3CutNameVector();
  45   FillG3FlagNameVector();
  46
  47   // fill process name map
  48   FillProcessMap();
  49 }
  50
  51 TG4PhysicsManager::TG4PhysicsManager(const TG4PhysicsManager& right) {
  52 //
  53   TG4Globals::Exception(
  54     "Attempt to copy TG4PhysicsManager singleton.");
  55 }
  56
  57 TG4PhysicsManager::~TG4PhysicsManager() {
  58 //
  59   delete fIsCutVector;
  60   delete fIsFlagVector;
  61 }
  62
  63 // operators
  64
  65 TG4PhysicsManager&
  66 TG4PhysicsManager::operator=(const TG4PhysicsManager& right)
  67 {
  68   // check assignement to self
  69   if (this == &right) return *this;
  70
  71   TG4Globals::Exception(
  72     "Attempt to assign TG4PhysicsManager singleton.");
  73
  74   return *this;
  75 }
  76
  77 // private methods
  78
  79 void TG4PhysicsManager::LockException() const
  80 {
  81 // Gives exception in case of attempt to modified physics
  82 // setup after physics manager was locked.
  83 // ---
  84
  85   G4String text = "TG4PhysicsManager: \n";
  86   text = text + "    It is too late to change physics setup. \n";
  87   text = text + "    PhysicsManager has been already locked.";
  88   TG4Globals::Exception(text);
  89 }
  90
  91 void TG4PhysicsManager::FillG3CutNameVector()
  92 {
  93 // Defines fCutNameVector.
  94 // ---
  95
  96   fG3CutNameVector.insert("CUTGAM");
  97   fG3CutNameVector.insert("CUTELE");
  98   fG3CutNameVector.insert("CUTNEU");
  99   fG3CutNameVector.insert("CUTHAD");
 100   fG3CutNameVector.insert("CUTMUO");
 101   fG3CutNameVector.insert("BCUTE");
 102   fG3CutNameVector.insert("BCUTM");
 103   fG3CutNameVector.insert("DCUTE");
 104   fG3CutNameVector.insert("DCUTM");
 105   fG3CutNameVector.insert("PPCUTM");
 106 }
 107
 108 void TG4PhysicsManager::FillG3FlagNameVector()
 109 {
 110 // Defines fFlagNameVector.
 111 // ---
 112
 113   fG3FlagNameVector.insert("PAIR");
 114   fG3FlagNameVector.insert("COMP");
 115   fG3FlagNameVector.insert("PHOT");
 116   fG3FlagNameVector.insert("PFIS");
 117   fG3FlagNameVector.insert("DRAY");
 118   fG3FlagNameVector.insert("ANNI");
 119   fG3FlagNameVector.insert("BREM");
 120   fG3FlagNameVector.insert("HADR");
 121   fG3FlagNameVector.insert("MUNU");
 122   fG3FlagNameVector.insert("DCAY");
 123   fG3FlagNameVector.insert("LOSS");
 124   fG3FlagNameVector.insert("MULS");
 125 }
 126
 127 void TG4PhysicsManager::FillProcessMap()
 128 {
 129 // Fills fProcessMap.
 130 // The default G4 process names are used in the map.
 131 // ---
 132
 133   // multiple scattering
 134   fProcessMap.Add("msc",  kPMultipleScattering);
 135   fProcessMap.Add("Imsc", kPMultipleScattering);
 136
 137   // continuous energy loss
 138   // !! including delta rays
 139   fProcessMap.Add("eIoni",  kPEnergyLoss);
 140   fProcessMap.Add("IeIoni", kPEnergyLoss);
 141   fProcessMap.Add("LowEnergyIoni", kPEnergyLoss);
 142   fProcessMap.Add("hIoni",  kPEnergyLoss);
 143   fProcessMap.Add("IhIoni", kPEnergyLoss);
 144   fProcessMap.Add("hLowEIoni", kPEnergyLoss);
 145   fProcessMap.Add("MuIoni", kPEnergyLoss);
 146   fProcessMap.Add("IMuIonisation", kPEnergyLoss);
 147   fProcessMap.Add("ionIoni",  kPEnergyLoss);
 148   fProcessMap.Add("ionLowEIoni",  kPEnergyLoss);
 149   fProcessMap.Add("PAIonisation",  kPEnergyLoss);
 150
 151   // bending in mag. field
 152   // kPMagneticFieldL
 153
 154   // particle decay
 155   fProcessMap.Add("Decay", kPDecay);
 156
 157   // photon pair production or
 158   // muon direct pair production
 159   fProcessMap.Add("conv", kPPair);
 160   fProcessMap.Add("LowEnConversion", kPPair);
 161   fProcessMap.Add("MuPairProd", kPPair);
 162   fProcessMap.Add("IMuPairProduction", kPPair);
 163
 164   // Compton scattering
 165   fProcessMap.Add("compt", kPCompton);
 166   fProcessMap.Add("LowEnCompton", kPCompton);
 167   fProcessMap.Add("polarCompt", kPCompton);
 168
 169   // photoelectric effect
 170   fProcessMap.Add("phot", kPPhotoelectric);
 171   fProcessMap.Add("LowEnPhotoElec", kPPhotoelectric);
 172
 173   // bremsstrahlung
 174   fProcessMap.Add("eBrem", kPBrem);
 175   fProcessMap.Add("IeBrem", kPBrem);
 176   fProcessMap.Add("MuBrem", kPBrem);
 177   fProcessMap.Add("IMuBremsstrahlung", kPBrem);
 178   fProcessMap.Add("LowEnBrem", kPBrem);
 179
 180   // delta-ray production
 181   // kPDeltaRay
 182   // has to be distinguished from kPEnergyLoss on flight
 183
 184   // positron annihilation
 185   fProcessMap.Add("annihil", kPAnnihilation);
 186   fProcessMap.Add("Iannihil", kPAnnihilation);
 187
 188   // hadronic interaction
 189   // kPHadronic
 190
 191   // nuclear evaporation
 192   // kPEvaporation
 193
 194   // nuclear fission
 195   // kPNuclearFission
 196
 197   // nuclear absorption
 198   fProcessMap.Add("PionMinusAbsorptionAtRest", kPNuclearAbsorption);
 199   fProcessMap.Add("PiMinusAbsorptionAtRest", kPNuclearAbsorption);
 200   fProcessMap.Add("KaonMinusAbsorption", kPNuclearAbsorption);
 201   fProcessMap.Add("KaonMinusAbsorptionAtRest", kPNuclearAbsorption);
 202
 203   // antiproton annihilation
 204   fProcessMap.Add("AntiProtonAnnihilationAtRest", kPPbarAnnihilation);
 205   // fProcessMap.Add("AntiNeutronAnnihilationAtRest", not defined);
 206
 207   // neutron capture
 208   fProcessMap.Add("NeutronCaptureAtRest", kPNCapture);
 209   // fProcessMap.Add("LCapture", hadron capture not defined);
 210
 211   // hadronic elastic incoherent scattering
 212   fProcessMap.Add("LElastic", kPHElastic);
 213
 214   // hadronic inelastic scattering
 215   fProcessMap.Add("inelastic", kPHInhelastic);
 216
 217   // muon nuclear interaction
 218   fProcessMap.Add("MuNucl", kPMuonNuclear);
 219
 220   // exceeded time of flight cut
 221   // kPTOFlimit
 222
 223   // nuclear photofission
 224   // kPPhotoFission
 225
 226   // Rayleigh scattering
 227   fProcessMap.Add("Rayleigh Scattering", kPRayleigh);
 228
 229   // no mechanism is active, usually at the entrance of a new volume
 230   // kPNull
 231
 232   // particle has fallen below energy threshold and tracking stops
 233   // kPStop
 234
 235   // Cerenkov photon absorption
 236   fProcessMap.Add("Absorption", kPLightAbsorption);
 237
 238   // Cerenkov photon reflection/refraction
 239   // kPLightScattering, kPLightReflection, kPLightRefraction
 240   // has to be inquired from the G4OpBoundary process
 241
 242   // synchrotron radiation
 243   fProcessMap.Add("SynchrotronRadiation", kPSynchrotron);
 244 }
 245
 246
 247 G4int TG4PhysicsManager::GetPDGEncoding(G4ParticleDefinition* particle)
 248 {
 249 // Returns the PDG code of particle;
 250 // if standard PDG code is not defined the TDatabasePDG
 251 // is used.
 252 // ---
 253
 254   // get PDG encoding from G4 particle definition
 255   G4int pdgEncoding = particle->GetPDGEncoding();
 256
 257   if (pdgEncoding == 0) {
 258     // get PDG encoding from TDatabasePDG
 259
 260     // get particle name from the name map
 261     G4String g4name = particle->GetParticleName();
 262     G4String tname = fParticleNameMap.GetSecond(g4name);
 263     if (tname == "Undefined") {
 264       G4String text = "TG4PhysicsManager::GetPDGEncoding: \n";
 265       text = text + "    Particle " + g4name;
 266       text = text + " was not found in the name map.";
 267       TG4Globals::Exception(text);
 268     }
 269
 270     // get particle from TDatabasePDG
 271     TDatabasePDG* pdgDB = TDatabasePDG::Instance();
 272     TParticlePDG* particle = pdgDB->GetParticle(tname);
 273     if (!particle) {
 274       G4String text = "TG4PhysicsManager::GetPDGEncoding: \n";
 275       text = text + "    Particle " + tname;
 276       text = text + " was not found in TDatabasePDG.";
 277       TG4Globals::Exception(text);
 278     }
 279
 280     // get PDG encoding
 281     pdgEncoding = particle->PdgCode();
 282   }
 283
 284   return pdgEncoding;
 285 }
 286
 287 AliMCProcess TG4PhysicsManager::GetMCProcess(const G4String& g4ProcessName)
 288 {
 289 // Returns the AliMCProcess code of process specified by name.
 290 // ---
 291
 292   G4int processCode = fProcessMap.GetSecond(g4ProcessName);
 293
 294   if (processCode == 0) return kPNoProcess;
 295
 296   return (AliMCProcess)processCode;
 297 }
 298
 299 G4int TG4PhysicsManager::GetPDGEncoding(G4String particleName)
 300 {
 301 // Returns the PDG code of particle specified by name.
 302 // ---
 303
 304   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
 305
 306   G4ParticleDefinition* particle = 0;
 307   particle = particleTable->FindParticle(particleName);
 308   if (!particle) {
 309     G4String text = "TG4PhysicsManager::GetPDGEncoding:\n";
 310     text = text + "   G4ParticleTable::FindParticle() " + particleName;
 311     text = text + " failed.";
 312     TG4Globals::Exception(text);
 313   }
 314
 315   return GetPDGEncoding(particle);
 316 }
 317
 318 void TG4PhysicsManager::SetCut(TG3Cut cut, G4double cutValue)
 319 {
 320 // Sets kinetic energy cut (in a G3-like way).
 321 // ---
 322
 323   if (!fCutVector) {
 324     // create vector of kinetic energy cut values
 325     fCutVector = new TG4CutVector();
 326   }
 327   fCutVector->SetG3Cut(cut, cutValue);
 328   SwitchIsCutVector(cut);
 329 }
 330
 331 void TG4PhysicsManager::SetProcess(TG3Flag flag, G4int flagValue)
 332 {
 333 // Sets control process flag (in a G3-like way).
 334 // ---
 335
 336   if (!fFlagVector) {
 337     // create vector of control process flag values
 338     fFlagVector = new TG4FlagVector;
 339   }
 340   fFlagVector->SetG3Flag(flag, flagValue);
 341 }
 342
 343
 344 Float_t TG4PhysicsManager::Xsec(char* ch, Float_t p1, Int_t i1, Int_t i2)
 345 {
 346 // Not yet implemented -> gives exception.
 347 // ---
 348
 349   TG4Globals::Exception(
 350     "TG4PhysicsManager::Xsec: not yet implemented.");
 351
 352   return 0.;
 353 }
 354
 355 Int_t TG4PhysicsManager::IdFromPDG(Int_t pdgID) const
 356 {
 357 // G4 does not use the integer particle identifiers
 358 // Id <-> PDG is identity.
 359 // ---
 360
 361   return pdgID;
 362 }
 363
 364 Int_t TG4PhysicsManager::PDGFromId(Int_t mcID) const
 365 {
 366 // G4 does not use integer particle identifiers
 367 // Id <-> PDG is identity.
 368 // ---
 369
 370   return mcID;
 371 }
 372
 373 void  TG4PhysicsManager::DefineParticles()
 374 {
 375   // ======
 376   // Taken from TGeant3
 377   //
 378   // Use ENDF-6 mapping for ions = 10000*z+10*a+iso
 379   // and add 1 000 000
 380   // and numbers above 5 000 000 for special applications
 381   //
 382
 383   const Int_t kion=10000000;
 384   const Int_t kspe=50000000;
 385
 386   const Double_t kGeV=0.9314943228;
 387   const Double_t kHslash = 1.0545726663e-27;
 388   const Double_t kErgGeV = 1/1.6021773349e-3;
 389   const Double_t kHshGeV = kHslash*kErgGeV;
 390   const Double_t kYearsToSec = 3600*24*365.25;
 391
 392   TDatabasePDG *pdgDB = TDatabasePDG::Instance();
 393
 394   pdgDB->AddParticle("Deuteron","Deuteron",2*kGeV+8.071e-3,kTRUE,
 395                      0,1,"Ion",kion+10020);
 396
 397   pdgDB->AddParticle("Triton","Triton",3*kGeV+14.931e-3,kFALSE,
 398                      kHshGeV/(12.33*kYearsToSec),1,"Ion",kion+10030);
 399
 400   pdgDB->AddParticle("Alpha","Alpha",4*kGeV+2.424e-3,kTRUE,
 401                      kHshGeV/(12.33*kYearsToSec),2,"Ion",kion+20040);
 402
 403   pdgDB->AddParticle("HE3","HE3",3*kGeV+14.931e-3,kFALSE,
 404                      0,2,"Ion",kion+20030);
 405
 406   pdgDB->AddParticle("Cherenkov","Cherenkov",0,kFALSE,
 407                      0,0,"Special",kspe+50);
 408
 409   pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,
 410                      0,0,"Special",kspe+51);
 411
 412
 413   // To do: define the PDG database extension
 414   // in a common part.
 415   //
 416   // AliMC::ExtendPDGDatabase();
 417   //
 418   // end of "common" implementation
 419   // ======
 420
 421   // map G4 particle names to TDatabasePDG names
 422   // (the map is built only for particles that have not
 423   //  defined standard PDG encoding)
 424
 425   fParticleNameMap.Add("deuteron","Deuteron");
 426   fParticleNameMap.Add("triton",  "Triton");
 427   fParticleNameMap.Add("alpha",   "Alpha");
 428   fParticleNameMap.Add("He3",     "HE3");
 429   fParticleNameMap.Add("opticalphoton","Cherenkov");
 430   // fParticleNameMap.Add("???","FeedbackPhoton");
 431   fParticleNameMap.Add("geantino", "Rootino");
 432
 433   // map G4 particle names to TDatabasePDG encodings
 434   fParticlePDGMap.Add("deuteron", GetPDGEncoding("deuteron"));
 435   fParticlePDGMap.Add("triton", GetPDGEncoding("triton"));
 436   fParticlePDGMap.Add("alpha", GetPDGEncoding("alpha"));
 437   fParticlePDGMap.Add("He3", GetPDGEncoding("He3") );
 438   fParticlePDGMap.Add("opticalphoton", GetPDGEncoding("opticalphoton"));
 439   // fParticlePDGMap.Add("???","FeedbackPhoton");
 440   fParticleNameMap.Add("geantino", GetPDGEncoding("geantino"));
 441
 442   // add verbose
 443   G4cout << "Particle maps have been filled." << G4endl;
 444   //fParticleNameMap.PrintAll();
 445   //fParticlePDGMap.PrintAll();
 446 }
 447
 448 void TG4PhysicsManager::SwitchIsCutVector(TG3Cut cut)
 449 {
 450 // Updates the vector of booleans (fIsCutVector) for the specified cut.
 451 // ---
 452
 453   switch (cut) {
 454     case kCUTGAM:
 455            (*fIsCutVector)[kGamma] = true;
 456            break;
 457     case kBCUTE:
 458            (*fIsCutVector)[kGamma] = true;
 459            break;
 460     case kBCUTM:
 461            (*fIsCutVector)[kGamma] = true;
 462            break;
 463     case kCUTELE:
 464            (*fIsCutVector)[kElectron] = true;
 465            break;
 466     case kDCUTE:
 467            (*fIsCutVector)[kElectron] = true;
 468            break;
 469     case kDCUTM:
 470            (*fIsCutVector)[kElectron] = true;
 471            break;
 472     case kCUTNEU:
 473            (*fIsCutVector)[kNeutralHadron] = true;
 474            break;
 475     case kCUTHAD:
 476            (*fIsCutVector)[kChargedHadron] = true;
 477            break;
 478     case kCUTMUO:
 479            (*fIsCutVector)[kMuon] = true;
 480            break;
 481     default:
 482            break;
 483   }
 484 }
 485
 486 void TG4PhysicsManager::SwitchIsFlagVector(TG3Flag flag)
 487 {
 488 // Updates the vector of booleans (fIsFlagVector) for the specified flag.
 489 // ---
 490
 491   switch (flag) {
 492     case kPAIR:
 493            // gamma
 494            (*fIsFlagVector)[kGamma] = true;
 495            break;
 496     case kCOMP:
 497            // gamma
 498            (*fIsFlagVector)[kGamma] = true;
 499            break;
 500     case kPHOT:
 501            // gamma
 502            (*fIsFlagVector)[kGamma] = true;
 503            break;
 504     case kPFIS:
 505            // gamma
 506            (*fIsFlagVector)[kGamma] = true;
 507            break;
 508     case kDRAY:
 509            // all charged particles
 510            (*fIsFlagVector)[kElectron] = true;
 511            (*fIsFlagVector)[kEplus] = true;
 512            (*fIsFlagVector)[kChargedHadron] = true;
 513            (*fIsFlagVector)[kMuon] = true;
 514            break;
 515     case kANNI:
 516            // e+ only
 517            (*fIsFlagVector)[kEplus] = true;
 518            break;
 519     case kBREM:
 520            // e-/e+, muons
 521            (*fIsFlagVector)[kElectron] = true;
 522            (*fIsFlagVector)[kEplus] = true;
 523            (*fIsFlagVector)[kMuon] = true;
 524            break;
 525     case kHADR:
 526            // hadrons
 527            (*fIsFlagVector)[kNeutralHadron] = true;
 528            (*fIsFlagVector)[kChargedHadron] = true;
 529            break;
 530     case kMUNU:
 531            // muons
 532            (*fIsFlagVector)[kMuon] = true;
 533            break;
 534     case kDCAY:
 535            // any
 536            (*fIsFlagVector)[kAny] = true;
 537            break;
 538     case kLOSS:
 539            // all charged particles
 540            (*fIsFlagVector)[kElectron] = true;
 541            (*fIsFlagVector)[kEplus] = true;
 542            (*fIsFlagVector)[kChargedHadron] = true;
 543            (*fIsFlagVector)[kMuon] = true;
 544            break;
 545     case kMULS:
 546            // all charged particles
 547            (*fIsFlagVector)[kElectron] = true;
 548            (*fIsFlagVector)[kEplus] = true;
 549            (*fIsFlagVector)[kChargedHadron] = true;
 550            (*fIsFlagVector)[kMuon] = true;
 551            break;
 552     default:
 553           break;
 554   }
 555 }
 556
 557 TG3Cut TG4PhysicsManager::GetG3Cut(G4String cutName)
 558 {
 559 // Retrieves corresponding TG3Cut constant from the cutName.
 560 // ---
 561
 562   if      (cutName == fG3CutNameVector[kCUTGAM]) return kCUTGAM;
 563   else if (cutName == fG3CutNameVector[kBCUTE])  return kBCUTE;
 564   else if (cutName == fG3CutNameVector[kBCUTM])  return kBCUTM;
 565   else if (cutName == fG3CutNameVector[kCUTELE]) return kCUTELE;
 566   else if (cutName == fG3CutNameVector[kDCUTE])  return kDCUTE;
 567   else if (cutName == fG3CutNameVector[kDCUTM])  return kDCUTM;
 568   else if (cutName == fG3CutNameVector[kCUTNEU]) return kCUTNEU;
 569   else if (cutName == fG3CutNameVector[kCUTHAD]) return kCUTHAD;
 570   else if (cutName == fG3CutNameVector[kCUTMUO]) return kCUTMUO;
 571   else return kNoG3Cuts;
 572 }
 573
 574 TG3Flag TG4PhysicsManager::GetG3Flag(G4String flagName)
 575 {
 576 // Retrieves corresponding TG3Flag constant from the flagName.
 577 // ---
 578
 579   if      (flagName == fG3FlagNameVector[kPAIR]) return kPAIR;
 580   else if (flagName == fG3FlagNameVector[kCOMP]) return kCOMP;
 581   else if (flagName == fG3FlagNameVector[kPHOT]) return kPHOT;
 582   else if (flagName == fG3FlagNameVector[kPFIS]) return kPFIS;
 583   else if (flagName == fG3FlagNameVector[kDRAY]) return kDRAY;
 584   else if (flagName == fG3FlagNameVector[kANNI]) return kANNI;
 585   else if (flagName == fG3FlagNameVector[kBREM]) return kBREM;
 586   else if (flagName == fG3FlagNameVector[kHADR]) return kHADR;
 587   else if (flagName == fG3FlagNameVector[kMUNU]) return kMUNU;
 588   else if (flagName == fG3FlagNameVector[kDCAY]) return kDCAY;
 589   else if (flagName == fG3FlagNameVector[kLOSS]) return kLOSS;
 590   else if (flagName == fG3FlagNameVector[kMULS]) return kMULS;
 591   else return kNoG3Flags;
 592 }
 593
 594 // public methods
 595
 596 void TG4PhysicsManager::BuildPhysics()
 597 {
 598 // Empty function - not needed in G4.
 599 // (Physics is built within /run/initialize.)
 600
 601   TG4Globals::Warning(
 602     "TG4PhysicsManager::BuildPhysics: is empty function in G4 MC.");
 603 }
 604
 605 void TG4PhysicsManager::SetCut(const char* cutName, Float_t cutValue)
 606 {
 607 // Sets the specified cut.
 608 // ---
 609
 610   if (fLock) LockException();
 611   TG3Cut g3Cut = GetG3Cut(cutName);
 612   if (g3Cut != kNoG3Cuts)
 613     SetCut(g3Cut, cutValue);
 614   else {
 615     G4String text = "TG4PhysicsManager::SetCut:\n";
 616     text = text + "    Parameter " + cutName;
 617     text = text + " is not implemented.";
 618     TG4Globals::Warning(text);
 619   }
 620 }
 621
 622 void TG4PhysicsManager::SetProcess(const char* flagName, Int_t flagValue)
 623 {
 624 // Sets the specified process control.
 625 // ---
 626
 627   if (fLock) LockException();
 628   TG3Flag g3Flag = GetG3Flag(flagName);
 629   if (g3Flag != kNoG3Flags)
 630     SetProcess(g3Flag, flagValue);
 631   else {
 632     G4String text = "TG4PhysicsManager::SetProcess:\n";
 633     text = text + "    Parameter " + flagName;
 634     text = text + " is not implemented.";
 635     TG4Globals::Warning(text);
 636   }
 637 }
 638
 639 void TG4PhysicsManager::SetProcessActivation()
 640 {
 641 // (In)Activates built processes according
 642 // to the setup in fFlagVector.
 643 // ---
 644
 645   if (fPhysicsList) {
 646     // temporarily excluded
 647     // fPhysicsList->SetProcessActivation();
 648   }
 649   else {
 650     G4String text = "TG4PhysicsManager::SetProcessActivation:\n";
 651     text = text +   "   There is no physics list set.";
 652     TG4Globals::Exception(text);
 653   }
 654 }
 655
 656 G4int TG4PhysicsManager::GetPDGEncodingFast(G4ParticleDefinition* particle)
 657 {
 658 // Returns the PDG code of particle;
 659 // if standard PDG code is not defined the preregistred
 660 // fParticlePDGMap is used.
 661 // ---
 662
 663   // get PDG encoding from G4 particle definition
 664   G4int pdgEncoding = particle->GetPDGEncoding();
 665
 666   if (pdgEncoding == 0) {
 667     // use FParticlePDGMap if standard PDG code is not defined
 668     G4String name = particle->GetParticleName();
 669     pdgEncoding = fParticlePDGMap.GetSecond(name);
 670   }
 671
 672   return pdgEncoding;
 673 }
 674
 675 G4bool TG4PhysicsManager::CheckCutWithCutVector(G4String name,
 676                              G4double value, TG3Cut& cut)
 677 {
 678 // Retrieves corresponding TG3Cut from the name and
 679 // in case the value is different from the value in cutVector
 680 // sets true the value of the fIsCutVector element
 681 // corresponding to this cut and returns true;
 682 // returns false otherwise.
 683 // ---
 684
 685   // convert cut name -> TG3Cut
 686   cut = GetG3Cut(name);
 687
 688   // set switch vector element only if the value
 689   // is different from the value in cutVector
 690   if (cut !=kNoG3Cuts) {
 691     // get tolerance from TG4G3Defaults in GeV
 692     G4double tolerance = TG4G3Defaults::CutTolerance()/GeV;
 693     if (!(fCutVector) || (abs(value - (*fCutVector)[cut]) > tolerance)) {
 694       SwitchIsCutVector(cut);
 695       return true;
 696     }
 697     else return false;
 698   }
 699   return false;
 700 }
 701
 702 G4bool TG4PhysicsManager::CheckFlagWithFlagVector(G4String name,
 703                               G4double value, TG3Flag& flag)
 704 {
 705 // Retrieves corresponding TG3Flag from the name and
 706 // in case the value is different from the value in flagVector
 707 // sets true the value of the fIsFlagVector element
 708 // corresponding to this flag and returns true;
 709 // returns false otherwise.
 710 // ---
 711
 712   // convert flag name -> TG3Flag
 713   flag = GetG3Flag(name);
 714
 715   // set switch vector element only if the value
 716   // is different from the value in flagVector
 717   if (flag !=kNoG3Flags) {
 718     if (!(fFlagVector) || (abs(value - (*fFlagVector)[flag]) > 0.01)) {
 719       SwitchIsFlagVector(flag);
 720       return true;
 721     }
 722     else return false;
 723   }
 724   return false;
 725 }
 726
 727 G4bool TG4PhysicsManager::CheckCutWithG3Defaults(G4String name,
 728                               G4double value, TG3Cut& cut)
 729 {
 730 // Retrieves corresponding TG3Cut from the name and
 731 // in case the value is different from the G3 default value
 732 // sets true the value of the SwitchCutVector element
 733 // corresponding to this cut and returns true;
 734 // returns false otherwise.
 735 // ---
 736
 737   // convert cut name -> TG3Cut
 738   cut = GetG3Cut(name);
 739
 740   // set switch vector element only if the value
 741   // is different from G3 default
 742   if (cut !=kNoG3Cuts) {
 743     if (!TG4G3Defaults::IsDefaultCut(cut, value)) {
 744       SwitchIsCutVector(cut);
 745       return true;
 746     }
 747     else return false;
 748   }
 749   return false;
 750 }
 751
 752 G4bool TG4PhysicsManager::CheckFlagWithG3Defaults(G4String name,
 753                               G4double value, TG3Flag& flag)
 754 {
 755 // Retrieves corresponding TG3Flag from the name and
 756 // in case the value is different from the G3 default value
 757 // sets true the value of the SwitchFlagVector element
 758 // corresponding to this flag and returns true;
 759 // returns false otherwise.
 760 // ---
 761
 762   // convert flag name -> TG3Flag
 763   flag = GetG3Flag(name);
 764
 765   // set switch vector element only if the value
 766   // is different from G3 default
 767   if (flag !=kNoG3Flags) {
 768     if (!TG4G3Defaults::IsDefaultFlag(flag, value)) {
 769       SwitchIsFlagVector(flag);
 770       return true;
 771     }
 772     else return false;
 773   }
 774   return false;
 775 }
 776
 777 void TG4PhysicsManager::SetG3DefaultCuts()
 778 {
 779 // Sets G3 default values of kinetic energy cuts.
 780 // ---
 781
 782   if (fLock) LockException();
 783   if (!fCutVector) {
 784     // create vector of kinetic energy cut values
 785     fCutVector = new TG4CutVector();
 786   }
 787   fCutVector->SetG3Defaults();
 788 }
 789
 790 void TG4PhysicsManager::SetG3DefaultProcesses()
 791 {
 792 // Sets G3 default values of control process flags.
 793 // ---
 794
 795   if (fLock) LockException();
 796   if (!fFlagVector) {
 797     // create vector of control process flag values
 798     fFlagVector = new TG4FlagVector;
 799   }
 800   fFlagVector->SetG3Defaults();
 801 }
 802
 803 G4bool TG4PhysicsManager::IsSpecialCuts() const
 804 {
 805 // Returns true if any special cut value is set.
 806 // ---
 807
 808   for (G4int i=0; i<kNofParticlesWSP; i++)
 809   {  if ((*fIsCutVector)[i]) return true; }
 810
 811   return false;
 812 }
 813
 814 G4bool TG4PhysicsManager::IsSpecialFlags() const
 815 {
 816 // Returns true if any special flag value is set.
 817 // ---
 818
 819   for (G4int i=0; i<kNofParticlesWSP; i++)
 820   {  if ((*fIsFlagVector)[i]) return true; }
 821
 822   return false;
 823 }
 824
 825 TG3ParticleWSP TG4PhysicsManager::GetG3ParticleWSP(
 826                                   G4ParticleDefinition* particle) const
 827 {
 828 // Returns TG3ParticleWSP constant for the specified particle.
 829 // (See TG3ParticleWSP.h, too.)
 830 // ---
 831
 832   G4String name = particle->GetParticleName();
 833   G4String pType = particle->GetParticleType();
 834
 835   if (name == "gamma") {
 836     return kGamma;
 837   }
 838   else if (name == "e-") {
 839     return kElectron;
 840   }
 841   else if (name == "e+") {
 842     return kEplus;
 843   }
 844   else if (( pType == "baryon" || pType == "meson" || pType == "nucleus" )) {
 845     if (particle->GetPDGCharge() == 0) {
 846       return kNeutralHadron;
 847     }
 848     else
 849       return kChargedHadron;
 850   }
 851   else if ( name == "mu-" || name == "mu+" ) {
 852     return kMuon;
 853   }
 854   else {
 855     return kNofParticlesWSP;
 856   }
 857 }
 858
 859 void TG4PhysicsManager::GetG3ParticleWSPName(G4int particleWSP,
 860                                              G4String& name) const
 861 {
 862 // Fills the passed name with the name/type of particle specified
 863 // by TG3ParticleWSP constant.
 864 // (See TG3ParticleWSP.h, too.)
 865 // ---
 866
 867   switch (particleWSP) {
 868     case kGamma:
 869       name = "Gamma";
 870       break;
 871     case kElectron:
 872       name = "Electron";
 873       break;
 874     case kEplus:
 875       name = "Eplus";
 876       break;
 877     case kNeutralHadron:
 878       name = "NeutralHadron";
 879       break;
 880     case kChargedHadron:
 881       name = "ChargedHadron";
 882       break;
 883     case kMuon:
 884       name = "Muon";
 885       break;
 886     case kAny:
 887       name = "Any";
 888       break;
 889     case kNofParticlesWSP:
 890       name = "NoSP";
 891       break;
 892     default:
 893       G4String text = "TG4PhysicsList::GetG3ParticleWSPName:\n";
 894       text = text + "   Wrong particleWSP.";
 895       TG4Globals::Exception(text);
 896       break;
 897   }
 898 }
 899