TGeant4/TG4G3PhysicsManager.cxx

   1 // $Id$
   2 // Category: physics
   3 //
   4 // See the class description in the header file.
   5
   6 #include "TG4G3PhysicsManager.h"
   7 #include "TG4G3CutVector.h"
   8 #include "TG4G3ControlVector.h"
   9 #include "TG4G3Defaults.h"
  10
  11 #include <G4ParticleDefinition.hh>
  12 #include <G4VProcess.hh>
  13 #include <G4UImessenger.hh>
  14 #include <G4ProcessTable.hh>
  15
  16 TG4G3PhysicsManager* TG4G3PhysicsManager::fgInstance = 0;
  17
  18 TG4G3PhysicsManager::TG4G3PhysicsManager()
  19   : fLock(false),
  20     fCutVector(0),
  21     fControlVector(0)
  22 {
  23   if (fgInstance) {
  24     TG4Globals::Exception(
  25       "TG4G3PhysicsManager: attempt to create two instances of singleton.");
  26   }
  27
  28   fgInstance = this;
  29
  30   // initialize fIsCutVector
  31   fIsCutVector = new TG4boolVector;
  32   G4int i;
  33   //for (i=0; i<kNofParticlesWSP; i++) fIsCutVector->insert(false);
  34   for (i=0; i<kNofParticlesWSP; i++) fIsCutVector->push_back(false);
  35
  36   // initialize fIsControlVector
  37   fIsControlVector = new TG4boolVector;
  38   //for (i=0; i<kNofParticlesWSP; i++) fIsControlVector->insert(false);
  39   for (i=0; i<kNofParticlesWSP; i++) fIsControlVector->push_back(false);
  40
  41   // define fCutNameVector, fControlNameVector
  42   FillG3CutNameVector();
  43   FillG3ControlNameVector();
  44 }
  45
  46 TG4G3PhysicsManager::TG4G3PhysicsManager(const TG4G3PhysicsManager& right) {
  47 //
  48   TG4Globals::Exception(
  49     "Attempt to copy TG4G3PhysicsManager singleton.");
  50 }
  51
  52 TG4G3PhysicsManager::~TG4G3PhysicsManager() {
  53 //
  54   delete fIsCutVector;
  55   delete fIsControlVector;
  56 }
  57
  58 // operators
  59
  60 TG4G3PhysicsManager&
  61 TG4G3PhysicsManager::operator=(const TG4G3PhysicsManager& right)
  62 {
  63   // check assignement to self
  64   if (this == &right) return *this;
  65
  66   TG4Globals::Exception(
  67     "Attempt to assign TG4G3PhysicsManager singleton.");
  68
  69   return *this;
  70 }
  71
  72 // private methods
  73
  74 void TG4G3PhysicsManager::FillG3CutNameVector()
  75 {
  76 // Defines fCutNameVector.
  77 // ---
  78
  79   fG3CutNameVector.insert("CUTGAM");
  80   fG3CutNameVector.insert("CUTELE");
  81   fG3CutNameVector.insert("CUTNEU");
  82   fG3CutNameVector.insert("CUTHAD");
  83   fG3CutNameVector.insert("CUTMUO");
  84   fG3CutNameVector.insert("BCUTE");
  85   fG3CutNameVector.insert("BCUTM");
  86   fG3CutNameVector.insert("DCUTE");
  87   fG3CutNameVector.insert("DCUTM");
  88   fG3CutNameVector.insert("PPCUTM");
  89 }
  90
  91 void TG4G3PhysicsManager::FillG3ControlNameVector()
  92 {
  93 // Defines fControlNameVector.
  94 // ---
  95
  96   fG3ControlNameVector.insert("PAIR");
  97   fG3ControlNameVector.insert("COMP");
  98   fG3ControlNameVector.insert("PHOT");
  99   fG3ControlNameVector.insert("PFIS");
 100   fG3ControlNameVector.insert("DRAY");
 101   fG3ControlNameVector.insert("ANNI");
 102   fG3ControlNameVector.insert("BREM");
 103   fG3ControlNameVector.insert("HADR");
 104   fG3ControlNameVector.insert("MUNU");
 105   fG3ControlNameVector.insert("DCAY");
 106   fG3ControlNameVector.insert("LOSS");
 107   fG3ControlNameVector.insert("MULS");
 108 }
 109
 110 void TG4G3PhysicsManager::SetCut(TG4G3Cut cut, G4double cutValue)
 111 {
 112 // Sets kinetic energy cut (in a G3-like way).
 113 // ---
 114
 115   if (!fCutVector) {
 116     // create vector of kinetic energy cut values
 117     fCutVector = new TG4G3CutVector();
 118   }
 119   fCutVector->SetG3Cut(cut, cutValue);
 120   SwitchIsCutVector(cut);
 121 }
 122
 123 void TG4G3PhysicsManager::SetProcess(TG4G3Control control, G4int controlValue)
 124 {
 125 // Sets control process control (in a G3-like way).
 126 // ---
 127
 128   if (!fControlVector) {
 129     // create vector of control process control values
 130     fControlVector = new TG4G3ControlVector;
 131   }
 132   fControlVector->SetG3Control(control, controlValue);
 133 }
 134
 135
 136 void TG4G3PhysicsManager::SwitchIsCutVector(TG4G3Cut cut)
 137 {
 138 // Updates the vector of booleans (fIsCutVector) for the specified cut.
 139 // ---
 140
 141   switch (cut) {
 142     case kCUTGAM:
 143            (*fIsCutVector)[kGamma] = true;
 144            break;
 145     case kBCUTE:
 146            (*fIsCutVector)[kGamma] = true;
 147            break;
 148     case kBCUTM:
 149            (*fIsCutVector)[kGamma] = true;
 150            break;
 151     case kCUTELE:
 152            (*fIsCutVector)[kElectron] = true;
 153            break;
 154     case kDCUTE:
 155            (*fIsCutVector)[kElectron] = true;
 156            break;
 157     case kDCUTM:
 158            (*fIsCutVector)[kElectron] = true;
 159            break;
 160     case kCUTNEU:
 161            (*fIsCutVector)[kNeutralHadron] = true;
 162            break;
 163     case kCUTHAD:
 164            (*fIsCutVector)[kChargedHadron] = true;
 165            break;
 166     case kCUTMUO:
 167            (*fIsCutVector)[kMuon] = true;
 168            break;
 169     default:
 170            break;
 171   }
 172 }
 173
 174 void TG4G3PhysicsManager::SwitchIsControlVector(TG4G3Control control)
 175 {
 176 // Updates the vector of booleans (fIsControlVector) for the specified control.
 177 // ---
 178
 179   switch (control) {
 180     case kPAIR:
 181            // gamma
 182            (*fIsControlVector)[kGamma] = true;
 183            break;
 184     case kCOMP:
 185            // gamma
 186            (*fIsControlVector)[kGamma] = true;
 187            break;
 188     case kPHOT:
 189            // gamma
 190            (*fIsControlVector)[kGamma] = true;
 191            break;
 192     case kPFIS:
 193            // gamma
 194            (*fIsControlVector)[kGamma] = true;
 195            break;
 196     case kDRAY:
 197            // all charged particles
 198            (*fIsControlVector)[kElectron] = true;
 199            (*fIsControlVector)[kEplus] = true;
 200            (*fIsControlVector)[kChargedHadron] = true;
 201            (*fIsControlVector)[kMuon] = true;
 202            break;
 203     case kANNI:
 204            // e+ only
 205            (*fIsControlVector)[kEplus] = true;
 206            break;
 207     case kBREM:
 208            // e-/e+, muons
 209            (*fIsControlVector)[kElectron] = true;
 210            (*fIsControlVector)[kEplus] = true;
 211            (*fIsControlVector)[kMuon] = true;
 212            break;
 213     case kHADR:
 214            // hadrons
 215            (*fIsControlVector)[kNeutralHadron] = true;
 216            (*fIsControlVector)[kChargedHadron] = true;
 217            break;
 218     case kMUNU:
 219            // muons
 220            (*fIsControlVector)[kMuon] = true;
 221            break;
 222     case kDCAY:
 223            // any
 224            (*fIsControlVector)[kAny] = true;
 225            break;
 226     case kLOSS:
 227            // all charged particles
 228            (*fIsControlVector)[kElectron] = true;
 229            (*fIsControlVector)[kEplus] = true;
 230            (*fIsControlVector)[kChargedHadron] = true;
 231            (*fIsControlVector)[kMuon] = true;
 232            break;
 233     case kMULS:
 234            // all charged particles
 235            (*fIsControlVector)[kElectron] = true;
 236            (*fIsControlVector)[kEplus] = true;
 237            (*fIsControlVector)[kChargedHadron] = true;
 238            (*fIsControlVector)[kMuon] = true;
 239            break;
 240     default:
 241           break;
 242   }
 243 }
 244
 245 TG4G3Cut TG4G3PhysicsManager::GetG3Cut(G4String cutName)
 246 {
 247 // Retrieves corresponding TG4G3Cut constant from the cutName.
 248 // ---
 249
 250   if      (cutName == fG3CutNameVector[kCUTGAM]) return kCUTGAM;
 251   else if (cutName == fG3CutNameVector[kBCUTE])  return kBCUTE;
 252   else if (cutName == fG3CutNameVector[kBCUTM])  return kBCUTM;
 253   else if (cutName == fG3CutNameVector[kCUTELE]) return kCUTELE;
 254   else if (cutName == fG3CutNameVector[kDCUTE])  return kDCUTE;
 255   else if (cutName == fG3CutNameVector[kDCUTM])  return kDCUTM;
 256   else if (cutName == fG3CutNameVector[kCUTNEU]) return kCUTNEU;
 257   else if (cutName == fG3CutNameVector[kCUTHAD]) return kCUTHAD;
 258   else if (cutName == fG3CutNameVector[kCUTMUO]) return kCUTMUO;
 259   else return kNoG3Cuts;
 260 }
 261
 262 TG4G3Control TG4G3PhysicsManager::GetG3Control(G4String controlName)
 263 {
 264 // Retrieves corresponding TG4G3Control constant from the controlName.
 265 // ---
 266
 267   if      (controlName == fG3ControlNameVector[kPAIR]) return kPAIR;
 268   else if (controlName == fG3ControlNameVector[kCOMP]) return kCOMP;
 269   else if (controlName == fG3ControlNameVector[kPHOT]) return kPHOT;
 270   else if (controlName == fG3ControlNameVector[kPFIS]) return kPFIS;
 271   else if (controlName == fG3ControlNameVector[kDRAY]) return kDRAY;
 272   else if (controlName == fG3ControlNameVector[kANNI]) return kANNI;
 273   else if (controlName == fG3ControlNameVector[kBREM]) return kBREM;
 274   else if (controlName == fG3ControlNameVector[kHADR]) return kHADR;
 275   else if (controlName == fG3ControlNameVector[kMUNU]) return kMUNU;
 276   else if (controlName == fG3ControlNameVector[kDCAY]) return kDCAY;
 277   else if (controlName == fG3ControlNameVector[kLOSS]) return kLOSS;
 278   else if (controlName == fG3ControlNameVector[kMULS]) return kMULS;
 279   else return kNoG3Controls;
 280 }
 281
 282 // public methods
 283
 284
 285 void TG4G3PhysicsManager::CheckLock()
 286 {
 287 // Gives exception in case the physics manager is locked.
 288 // Prevents from modifying physics setup after the physics manager is locked.
 289 // ---
 290
 291   if (fLock) {
 292     G4String text = "TG4PhysicsManager: \n";
 293     text = text + "    It is too late to change physics setup. \n";
 294     text = text + "    PhysicsManager has been already locked.";
 295     TG4Globals::Exception(text);
 296   }
 297 }
 298
 299
 300 G4VProcess* TG4G3PhysicsManager::FindProcess(G4String processName) const
 301 {
 302 // Finds G4VProcess with specified name.
 303 // ---
 304
 305   G4ProcessTable* processTable = G4ProcessTable::GetProcessTable();
 306
 307   G4ProcessVector* processVector
 308     = processTable->FindProcesses(processName);
 309   G4VProcess* firstFoundProcess = 0;
 310   if (processVector->entries()>0) firstFoundProcess= (*processVector)[0];
 311
 312   processVector->clear();
 313   delete processVector;
 314
 315   return firstFoundProcess;
 316 }
 317
 318
 319 G4bool TG4G3PhysicsManager::CheckCutWithTheVector(G4String name,
 320                                  G4double value, TG4G3Cut& cut)
 321 {
 322 // Retrieves corresponding TG4G3Cut from the name and
 323 // in case the value is different from the value in cutVector
 324 // sets true the value of the fIsCutVector element
 325 // corresponding to this cut and returns true;
 326 // returns false otherwise.
 327 // ---
 328
 329   // convert cut name -> TG4G3Cut
 330   cut = GetG3Cut(name);
 331
 332   // set switch vector element only if the value
 333   // is different from the value in cutVector
 334   if (cut !=kNoG3Cuts) {
 335     // get tolerance from TG4G3Defaults in GeV
 336     G4double tolerance = TG4G3Defaults::CutTolerance()/GeV;
 337     if (!(fCutVector) || (abs(value - (*fCutVector)[cut]) > tolerance)) {
 338       SwitchIsCutVector(cut);
 339       return true;
 340     }
 341     else return false;
 342   }
 343   return false;
 344 }
 345
 346 G4bool TG4G3PhysicsManager::CheckControlWithTheVector(G4String name,
 347                                  G4double value, TG4G3Control& control)
 348 {
 349 // Retrieves corresponding TG4G3Control from the name and
 350 // in case the value is different from the value in controlVector
 351 // sets true the value of the fIsControlVector element
 352 // corresponding to this control and returns true;
 353 // returns false otherwise.
 354 // ---
 355
 356   // convert control name -> TG4G3Control
 357   control = GetG3Control(name);
 358
 359   // set switch vector element only if the value
 360   // is different from the value in controlVector
 361   if (control !=kNoG3Controls) {
 362     if (!(fControlVector) || (abs(value - (*fControlVector)[control]) > 0.01)) {
 363       SwitchIsControlVector(control);
 364       return true;
 365     }
 366     else return false;
 367   }
 368   return false;
 369 }
 370
 371 G4bool TG4G3PhysicsManager::CheckCutWithG3Defaults(G4String name,
 372                                  G4double value, TG4G3Cut& cut)
 373 {
 374 // Retrieves corresponding TG4G3Cut from the name and
 375 // in case the value is different from the G3 default value
 376 // sets true the value of the SwitchCutVector element
 377 // corresponding to this cut and returns true;
 378 // returns false otherwise.
 379 // ---
 380
 381   // convert cut name -> TG4G3Cut
 382   cut = GetG3Cut(name);
 383
 384   // set switch vector element only if the value
 385   // is different from G3 default
 386   if (cut !=kNoG3Cuts) {
 387     if (!TG4G3Defaults::IsDefaultCut(cut, value)) {
 388       SwitchIsCutVector(cut);
 389       return true;
 390     }
 391     else return false;
 392   }
 393   return false;
 394 }
 395
 396 G4bool TG4G3PhysicsManager::CheckControlWithG3Defaults(G4String name,
 397                                  G4double value, TG4G3Control& control)
 398 {
 399 // Retrieves corresponding TG4G3Control from the name and
 400 // in case the value is different from the G3 default value
 401 // sets true the value of the SwitchControlVector element
 402 // corresponding to this control and returns true;
 403 // returns false otherwise.
 404 // ---
 405
 406   // convert control name -> TG4G3Control
 407   control = GetG3Control(name);
 408
 409   // set switch vector element only if the value
 410   // is different from G3 default
 411   if (control !=kNoG3Controls) {
 412     if (!TG4G3Defaults::IsDefaultControl(control, value)) {
 413       SwitchIsControlVector(control);
 414       return true;
 415     }
 416     else return false;
 417   }
 418   return false;
 419 }
 420
 421 void TG4G3PhysicsManager::SetG3DefaultCuts()
 422 {
 423 // Sets G3 default values of kinetic energy cuts.
 424 // ---
 425
 426   CheckLock();
 427   if (!fCutVector) {
 428     // create vector of kinetic energy cut values
 429     fCutVector = new TG4G3CutVector();
 430   }
 431   fCutVector->SetG3Defaults();
 432 }
 433
 434 void TG4G3PhysicsManager::SetG3DefaultControls()
 435 {
 436 // Sets G3 default values of control process controls.
 437 // ---
 438
 439   CheckLock();
 440   if (!fControlVector) {
 441     // create vector of control process control values
 442     fControlVector = new TG4G3ControlVector;
 443   }
 444   fControlVector->SetG3Defaults();
 445 }
 446
 447 G4bool TG4G3PhysicsManager::IsSpecialCuts() const
 448 {
 449 // Returns true if any special cut value is set.
 450 // ---
 451
 452   for (G4int i=0; i<kNofParticlesWSP; i++)
 453   {  if ((*fIsCutVector)[i]) return true; }
 454
 455   return false;
 456 }
 457
 458 G4bool TG4G3PhysicsManager::IsSpecialControls() const
 459 {
 460 // Returns true if any special control value is set.
 461 // ---
 462
 463   for (G4int i=0; i<kNofParticlesWSP; i++)
 464   {  if ((*fIsControlVector)[i]) return true; }
 465
 466   return false;
 467 }
 468
 469 TG4G3ParticleWSP TG4G3PhysicsManager::GetG3ParticleWSP(
 470                                       G4ParticleDefinition* particle) const
 471 {
 472 // Returns TG4G3ParticleWSP constant for the specified particle.
 473 // (See TG4G3ParticleWSP.h, too.)
 474 // ---
 475
 476   G4String name = particle->GetParticleName();
 477   G4String pType = particle->GetParticleType();
 478
 479   if (name == "gamma") {
 480     return kGamma;
 481   }
 482   else if (name == "e-") {
 483     return kElectron;
 484   }
 485   else if (name == "e+") {
 486     return kEplus;
 487   }
 488   else if (( pType == "baryon" || pType == "meson" || pType == "nucleus" )) {
 489     if (particle->GetPDGCharge() == 0) {
 490       return kNeutralHadron;
 491     }
 492     else
 493       return kChargedHadron;
 494   }
 495   else if ( name == "mu-" || name == "mu+" ) {
 496     return kMuon;
 497   }
 498   else {
 499     return kNofParticlesWSP;
 500   }
 501 }
 502
 503 void TG4G3PhysicsManager::GetG3ParticleWSPName(G4int particleWSP,
 504                                                G4String& name) const
 505 {
 506 // Fills the passed name with the name/type of particle specified
 507 // by TG4G3ParticleWSP constant.
 508 // (See TG4G3ParticleWSP.h, too.)
 509 // ---
 510
 511   switch (particleWSP) {
 512     case kGamma:
 513       name = "Gamma";
 514       break;
 515     case kElectron:
 516       name = "Electron";
 517       break;
 518     case kEplus:
 519       name = "Eplus";
 520       break;
 521     case kNeutralHadron:
 522       name = "NeutralHadron";
 523       break;
 524     case kChargedHadron:
 525       name = "ChargedHadron";
 526       break;
 527     case kMuon:
 528       name = "Muon";
 529       break;
 530     case kAny:
 531       name = "Any";
 532       break;
 533     case kNofParticlesWSP:
 534       name = "NoSP";
 535       break;
 536     default:
 537       G4String text = "TG4G3PhysicsManager::GetG3ParticleWSPName:\n";
 538       text = text + "   Wrong particleWSP.";
 539       TG4Globals::Exception(text);
 540       break;
 541   }
 542 }
 543