4 // See the class description in the header file.
6 #include "TG4PhysicsManager.h"
7 #include "TG4PhysicsList.h"
8 #include "TG4CutVector.h"
9 #include "TG4FlagVector.h"
10 #include "TG4G3Defaults.h"
12 #include <G4ParticleDefinition.hh>
13 #include <Randomize.hh>
15 #include <TDatabasePDG.h>
17 TG4PhysicsManager* TG4PhysicsManager::fgInstance = 0;
19 TG4PhysicsManager::TG4PhysicsManager()
27 TG4Globals::Exception(
28 "TG4PhysicsManager: attempt to create two instances of singleton.");
33 // initialize fIsCutVector
34 fIsCutVector = new TG4boolVector;
36 //for (i=0; i<kNofParticlesWSP; i++) fIsCutVector->insert(false);
37 for (i=0; i<kNofParticlesWSP; i++) fIsCutVector->push_back(false);
39 // initialize fIsFlagVector
40 fIsFlagVector = new TG4boolVector;
41 //for (i=0; i<kNofParticlesWSP; i++) fIsFlagVector->insert(false);
42 for (i=0; i<kNofParticlesWSP; i++) fIsFlagVector->push_back(false);
44 // define fCutNameVector
45 fG3CutNameVector.insert("CUTGAM");
46 fG3CutNameVector.insert("CUTELE");
47 fG3CutNameVector.insert("CUTNEU");
48 fG3CutNameVector.insert("CUTHAD");
49 fG3CutNameVector.insert("CUTMUO");
50 fG3CutNameVector.insert("BCUTE");
51 fG3CutNameVector.insert("BCUTM");
52 fG3CutNameVector.insert("DCUTE");
53 fG3CutNameVector.insert("DCUTM");
54 fG3CutNameVector.insert("PPCUTM");
56 // define fFlagNameVector
57 fG3FlagNameVector.insert("PAIR");
58 fG3FlagNameVector.insert("COMP");
59 fG3FlagNameVector.insert("PHOT");
60 fG3FlagNameVector.insert("PFIS");
61 fG3FlagNameVector.insert("DRAY");
62 fG3FlagNameVector.insert("ANNI");
63 fG3FlagNameVector.insert("BREM");
64 fG3FlagNameVector.insert("HADR");
65 fG3FlagNameVector.insert("MUNU");
66 fG3FlagNameVector.insert("DCAY");
67 fG3FlagNameVector.insert("LOSS");
68 fG3FlagNameVector.insert("MULS");
71 TG4PhysicsManager::TG4PhysicsManager(const TG4PhysicsManager& right) {
73 TG4Globals::Exception(
74 "Attempt to copy TG4PhysicsManager singleton.");
77 TG4PhysicsManager::~TG4PhysicsManager() {
86 TG4PhysicsManager::operator=(const TG4PhysicsManager& right)
88 // check assignement to self
89 if (this == &right) return *this;
91 TG4Globals::Exception(
92 "Attempt to assign TG4PhysicsManager singleton.");
99 void TG4PhysicsManager::LockException() const
101 // Gives exception in case of attempt to modified physics
102 // setup after physics manager was locked.
105 G4String text = "TG4PhysicsManager: \n";
106 text = text + " It is too late to change physics setup. \n";
107 text = text + " PhysicsManager has been already locked.";
108 TG4Globals::Exception(text);
111 G4int TG4PhysicsManager::GetPDGEncoding(G4ParticleDefinition* particle)
113 // Returns the PDG code of particle;
114 // if standard PDG code is not defined the TDatabasePDG
118 // get PDG encoding from G4 particle definition
119 G4int pdgEncoding = particle->GetPDGEncoding();
121 if (pdgEncoding == 0) {
122 // get PDG encoding from TDatabasePDG
124 // get particle name from the name map
125 G4String g4name = particle->GetParticleName();
126 G4String tname = fParticleNameMap.GetSecond(g4name);
127 if (tname == "Undefined") {
128 G4String text = "TG4PhysicsManager::GetPDGEncoding: \n";
129 text = text + " Particle " + g4name;
130 text = text + " was not found in the name map.";
131 TG4Globals::Exception(text);
134 // get particle from TDatabasePDG
135 TDatabasePDG* pdgDB = TDatabasePDG::Instance();
136 TParticlePDG* particle = pdgDB->GetParticle(tname);
138 G4String text = "TG4PhysicsManager::GetPDGEncoding: \n";
139 text = text + " Particle " + tname;
140 text = text + " was not found in TDatabasePDG.";
141 TG4Globals::Exception(text);
145 pdgEncoding = particle->PdgCode();
151 G4int TG4PhysicsManager::GetPDGEncoding(G4String particleName)
153 // Returns the PDG code of particle sepcified by name.
156 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
158 G4ParticleDefinition* particle = 0;
159 particle = particleTable->FindParticle(particleName);
161 G4String text = "TG4PhysicsManager::GetPDGEncoding:\n";
162 text = text + " G4ParticleTable::FindParticle() " + particleName;
163 text = text + " failed.";
164 TG4Globals::Exception(text);
167 return GetPDGEncoding(particle);
170 void TG4PhysicsManager::SetCut(TG3Cut cut, G4double cutValue)
172 // Sets kinetic energy cut (in a G3-like way).
176 // create vector of kinetic energy cut values
177 fCutVector = new TG4CutVector();
179 fCutVector->SetG3Cut(cut, cutValue);
180 SwitchIsCutVector(cut);
183 void TG4PhysicsManager::SetProcess(TG3Flag flag, G4int flagValue)
185 // Sets control process flag (in a G3-like way).
189 // create vector of control process flag values
190 fFlagVector = new TG4FlagVector;
192 fFlagVector->SetG3Flag(flag, flagValue);
196 Float_t TG4PhysicsManager::Xsec(char* ch, Float_t p1, Int_t i1, Int_t i2)
198 // Not yet implemented -> gives exception.
201 TG4Globals::Exception(
202 "TG4PhysicsManager::Xsec: not yet implemented.");
207 Int_t TG4PhysicsManager::IdFromPDG(Int_t pdgID) const
209 // G4 does not use the integer particle identifiers
210 // Id <-> PDG is identity.
216 Int_t TG4PhysicsManager::PDGFromId(Int_t mcID) const
218 // G4 does not use integer particle identifiers
219 // Id <-> PDG is identity.
225 void TG4PhysicsManager::DefineParticles()
228 // Taken from TGeant3
230 // Use ENDF-6 mapping for ions = 10000*z+10*a+iso
232 // and numbers above 5 000 000 for special applications
235 const Int_t kion=10000000;
236 const Int_t kspe=50000000;
238 const Double_t kGeV=0.9314943228;
239 const Double_t kHslash = 1.0545726663e-27;
240 const Double_t kErgGeV = 1/1.6021773349e-3;
241 const Double_t kHshGeV = kHslash*kErgGeV;
242 const Double_t kYearsToSec = 3600*24*365.25;
244 TDatabasePDG *pdgDB = TDatabasePDG::Instance();
246 pdgDB->AddParticle("Deuteron","Deuteron",2*kGeV+8.071e-3,kTRUE,
247 0,1,"Ion",kion+10020);
249 pdgDB->AddParticle("Triton","Triton",3*kGeV+14.931e-3,kFALSE,
250 kHshGeV/(12.33*kYearsToSec),1,"Ion",kion+10030);
252 pdgDB->AddParticle("Alpha","Alpha",4*kGeV+2.424e-3,kTRUE,
253 kHshGeV/(12.33*kYearsToSec),2,"Ion",kion+20040);
255 pdgDB->AddParticle("HE3","HE3",3*kGeV+14.931e-3,kFALSE,
256 0,2,"Ion",kion+20030);
258 pdgDB->AddParticle("Cherenkov","Cherenkov",0,kFALSE,
259 0,0,"Special",kspe+50);
261 pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,
262 0,0,"Special",kspe+51);
265 // To do: define the PDG database extension
268 // AliMC::ExtendPDGDatabase();
270 // end of "common" implementation
273 // map G4 particle names to TDatabasePDG names
274 // (the map is built only for particles that have not
275 // defined standard PDG encoding)
277 fParticleNameMap.Add("deuteron","Deuteron");
278 fParticleNameMap.Add("triton", "Triton");
279 fParticleNameMap.Add("alpha", "Alpha");
280 fParticleNameMap.Add("He3", "HE3");
281 fParticleNameMap.Add("opticalphoton","Cherenkov");
282 // fParticleNameMap.Add("???","FeedbackPhoton");
283 fParticleNameMap.Add("geantino", "Rootino");
285 // map G4 particle names to TDatabasePDG encodings
286 fParticlePDGMap.Add("deuteron", GetPDGEncoding("deuteron"));
287 fParticlePDGMap.Add("triton", GetPDGEncoding("triton"));
288 fParticlePDGMap.Add("alpha", GetPDGEncoding("alpha"));
289 fParticlePDGMap.Add("He3", GetPDGEncoding("He3") );
290 fParticlePDGMap.Add("opticalphoton", GetPDGEncoding("opticalphoton"));
291 // fParticlePDGMap.Add("???","FeedbackPhoton");
292 fParticleNameMap.Add("geantino", GetPDGEncoding("geantino"));
295 G4cout << "Particle maps have been filled." << endl;
296 //fParticleNameMap.PrintAll();
297 //fParticlePDGMap.PrintAll();
300 void TG4PhysicsManager::SwitchIsCutVector(TG3Cut cut)
302 // Updates the vector of booleans (fIsCutVector) for the specified cut.
307 (*fIsCutVector)[kGamma] = true;
310 (*fIsCutVector)[kGamma] = true;
313 (*fIsCutVector)[kGamma] = true;
316 (*fIsCutVector)[kElectron] = true;
319 (*fIsCutVector)[kElectron] = true;
322 (*fIsCutVector)[kElectron] = true;
325 (*fIsCutVector)[kNeutralHadron] = true;
328 (*fIsCutVector)[kChargedHadron] = true;
331 (*fIsCutVector)[kMuon] = true;
338 void TG4PhysicsManager::SwitchIsFlagVector(TG3Flag flag)
340 // Updates the vector of booleans (fIsFlagVector) for the specified flag.
346 (*fIsFlagVector)[kGamma] = true;
350 (*fIsFlagVector)[kGamma] = true;
354 (*fIsFlagVector)[kGamma] = true;
358 (*fIsFlagVector)[kGamma] = true;
361 // all charged particles
362 (*fIsFlagVector)[kElectron] = true;
363 (*fIsFlagVector)[kEplus] = true;
364 (*fIsFlagVector)[kChargedHadron] = true;
365 (*fIsFlagVector)[kMuon] = true;
369 (*fIsFlagVector)[kEplus] = true;
373 (*fIsFlagVector)[kElectron] = true;
374 (*fIsFlagVector)[kEplus] = true;
375 (*fIsFlagVector)[kMuon] = true;
379 (*fIsFlagVector)[kNeutralHadron] = true;
380 (*fIsFlagVector)[kChargedHadron] = true;
384 (*fIsFlagVector)[kMuon] = true;
388 (*fIsFlagVector)[kAny] = true;
391 // all charged particles
392 (*fIsFlagVector)[kElectron] = true;
393 (*fIsFlagVector)[kEplus] = true;
394 (*fIsFlagVector)[kChargedHadron] = true;
395 (*fIsFlagVector)[kMuon] = true;
398 // all charged particles
399 (*fIsFlagVector)[kElectron] = true;
400 (*fIsFlagVector)[kEplus] = true;
401 (*fIsFlagVector)[kChargedHadron] = true;
402 (*fIsFlagVector)[kMuon] = true;
409 TG3Cut TG4PhysicsManager::GetG3Cut(G4String cutName)
411 // Retrieves corresponding TG3Cut constant from the cutName.
414 if (cutName == fG3CutNameVector[kCUTGAM]) return kCUTGAM;
415 else if (cutName == fG3CutNameVector[kBCUTE]) return kBCUTE;
416 else if (cutName == fG3CutNameVector[kBCUTM]) return kBCUTM;
417 else if (cutName == fG3CutNameVector[kCUTELE]) return kCUTELE;
418 else if (cutName == fG3CutNameVector[kDCUTE]) return kDCUTE;
419 else if (cutName == fG3CutNameVector[kDCUTM]) return kDCUTM;
420 else if (cutName == fG3CutNameVector[kCUTNEU]) return kCUTNEU;
421 else if (cutName == fG3CutNameVector[kCUTHAD]) return kCUTHAD;
422 else if (cutName == fG3CutNameVector[kCUTMUO]) return kCUTMUO;
423 else return kNoG3Cuts;
426 TG3Flag TG4PhysicsManager::GetG3Flag(G4String flagName)
428 // Retrieves corresponding TG3Flag constant from the flagName.
431 if (flagName == fG3FlagNameVector[kPAIR]) return kPAIR;
432 else if (flagName == fG3FlagNameVector[kCOMP]) return kCOMP;
433 else if (flagName == fG3FlagNameVector[kPHOT]) return kPHOT;
434 else if (flagName == fG3FlagNameVector[kPFIS]) return kPFIS;
435 else if (flagName == fG3FlagNameVector[kDRAY]) return kDRAY;
436 else if (flagName == fG3FlagNameVector[kANNI]) return kANNI;
437 else if (flagName == fG3FlagNameVector[kBREM]) return kBREM;
438 else if (flagName == fG3FlagNameVector[kHADR]) return kHADR;
439 else if (flagName == fG3FlagNameVector[kMUNU]) return kMUNU;
440 else if (flagName == fG3FlagNameVector[kDCAY]) return kDCAY;
441 else if (flagName == fG3FlagNameVector[kLOSS]) return kLOSS;
442 else if (flagName == fG3FlagNameVector[kMULS]) return kMULS;
443 else return kNoG3Flags;
448 void TG4PhysicsManager::BuildPhysics()
450 // Empty function - not needed in G4.
451 // (Physics is built within /run/initialize.)
454 "TG4PhysicsManager::BuildPhysics: is empty function in G4 MC.");
457 void TG4PhysicsManager::SetCut(const char* cutName, Float_t cutValue)
459 // Sets the specified cut.
462 if (fLock) LockException();
463 TG3Cut g3Cut = GetG3Cut(cutName);
464 if (g3Cut != kNoG3Cuts)
465 SetCut(g3Cut, cutValue);
467 G4String text = "TG4PhysicsManager::SetCut:\n";
468 text = text + " Parameter " + cutName;
469 text = text + " is not implemented.";
470 TG4Globals::Warning(text);
474 void TG4PhysicsManager::SetProcess(const char* flagName, Int_t flagValue)
476 // Sets the specified process control.
479 if (fLock) LockException();
480 TG3Flag g3Flag = GetG3Flag(flagName);
481 if (g3Flag != kNoG3Flags)
482 SetProcess(g3Flag, flagValue);
484 G4String text = "TG4PhysicsManager::SetProcess:\n";
485 text = text + " Parameter " + flagName;
486 text = text + " is not implemented.";
487 TG4Globals::Warning(text);
491 void TG4PhysicsManager::SetProcessActivation()
493 // (In)Activates built processes according
494 // to the setup in fFlagVector.
498 // temporarily excluded
499 // fPhysicsList->SetProcessActivation();
502 G4String text = "TG4PhysicsManager::SetProcessActivation:\n";
503 text = text + " There is no physics list set.";
504 TG4Globals::Exception(text);
508 G4int TG4PhysicsManager::GetPDGEncodingFast(G4ParticleDefinition* particle)
510 // Returns the PDG code of particle;
511 // if standard PDG code is not defined the preregistred
512 // fParticlePDGMap is used.
515 // get PDG encoding from G4 particle definition
516 G4int pdgEncoding = particle->GetPDGEncoding();
518 if (pdgEncoding == 0) {
519 // use FParticlePDGMap if standard PDG code is not defined
520 G4String name = particle->GetParticleName();
521 pdgEncoding = fParticlePDGMap.GetSecond(name);
527 G4bool TG4PhysicsManager::CheckCutWithCutVector(G4String name,
528 G4double value, TG3Cut& cut)
530 // Retrieves corresponding TG3Cut from the name and
531 // in case the value is different from the value in cutVector
532 // sets true the value of the fIsCutVector element
533 // corresponding to this cut and returns true;
534 // returns false otherwise.
537 // convert cut name -> TG3Cut
538 cut = GetG3Cut(name);
540 // set switch vector element only if the value
541 // is different from the value in cutVector
542 if (cut !=kNoG3Cuts) {
543 // get tolerance from TG4G3Defaults in GeV
544 G4double tolerance = TG4G3Defaults::CutTolerance()/GeV;
545 if (!(fCutVector) || (abs(value - (*fCutVector)[cut]) > tolerance)) {
546 SwitchIsCutVector(cut);
554 G4bool TG4PhysicsManager::CheckFlagWithFlagVector(G4String name,
555 G4double value, TG3Flag& flag)
557 // Retrieves corresponding TG3Flag from the name and
558 // in case the value is different from the value in flagVector
559 // sets true the value of the fIsFlagVector element
560 // corresponding to this flag and returns true;
561 // returns false otherwise.
564 // convert flag name -> TG3Flag
565 flag = GetG3Flag(name);
567 // set switch vector element only if the value
568 // is different from the value in flagVector
569 if (flag !=kNoG3Flags) {
570 if (!(fFlagVector) || (abs(value - (*fFlagVector)[flag]) > 0.01)) {
571 SwitchIsFlagVector(flag);
579 G4bool TG4PhysicsManager::CheckCutWithG3Defaults(G4String name,
580 G4double value, TG3Cut& cut)
582 // Retrieves corresponding TG3Cut from the name and
583 // in case the value is different from the G3 default value
584 // sets true the value of the SwitchCutVector element
585 // corresponding to this cut and returns true;
586 // returns false otherwise.
589 // convert cut name -> TG3Cut
590 cut = GetG3Cut(name);
592 // set switch vector element only if the value
593 // is different from G3 default
594 if (cut !=kNoG3Cuts) {
595 if (!TG4G3Defaults::IsDefaultCut(cut, value)) {
596 SwitchIsCutVector(cut);
604 G4bool TG4PhysicsManager::CheckFlagWithG3Defaults(G4String name,
605 G4double value, TG3Flag& flag)
607 // Retrieves corresponding TG3Flag from the name and
608 // in case the value is different from the G3 default value
609 // sets true the value of the SwitchFlagVector element
610 // corresponding to this flag and returns true;
611 // returns false otherwise.
614 // convert flag name -> TG3Flag
615 flag = GetG3Flag(name);
617 // set switch vector element only if the value
618 // is different from G3 default
619 if (flag !=kNoG3Flags) {
620 if (!TG4G3Defaults::IsDefaultFlag(flag, value)) {
621 SwitchIsFlagVector(flag);
629 void TG4PhysicsManager::SetG3DefaultCuts()
631 // Sets G3 default values of kinetic energy cuts.
634 if (fLock) LockException();
636 // create vector of kinetic energy cut values
637 fCutVector = new TG4CutVector();
639 fCutVector->SetG3Defaults();
642 void TG4PhysicsManager::SetG3DefaultProcesses()
644 // Sets G3 default values of control process flags.
647 if (fLock) LockException();
649 // create vector of control process flag values
650 fFlagVector = new TG4FlagVector;
652 fFlagVector->SetG3Defaults();
655 G4bool TG4PhysicsManager::IsSpecialCuts() const
657 // Returns true if any special cut value is set.
660 for (G4int i=0; i<kNofParticlesWSP; i++)
661 { if ((*fIsCutVector)[i]) return true; }
666 G4bool TG4PhysicsManager::IsSpecialFlags() const
668 // Returns true if any special flag value is set.
671 for (G4int i=0; i<kNofParticlesWSP; i++)
672 { if ((*fIsFlagVector)[i]) return true; }
677 TG3ParticleWSP TG4PhysicsManager::GetG3ParticleWSP(
678 G4ParticleDefinition* particle) const
680 // Returns TG3ParticleWSP constant for the specified particle.
681 // (See TG3ParticleWSP.h, too.)
684 G4String name = particle->GetParticleName();
685 G4String pType = particle->GetParticleType();
687 if (name == "gamma") {
690 else if (name == "e-") {
693 else if (name == "e+") {
696 else if (( pType == "baryon" || pType == "meson" || pType == "nucleus" )) {
697 if (particle->GetPDGCharge() == 0) {
698 return kNeutralHadron;
701 return kChargedHadron;
703 else if ( name == "mu-" || name == "mu+" ) {
707 return kNofParticlesWSP;
711 void TG4PhysicsManager::GetG3ParticleWSPName(G4int particleWSP,
712 G4String& name) const
714 // Fills the passed name with the name/type of particle specified
715 // by TG3ParticleWSP constant.
716 // (See TG3ParticleWSP.h, too.)
719 switch (particleWSP) {
730 name = "NeutralHadron";
733 name = "ChargedHadron";
741 case kNofParticlesWSP:
745 G4String text = "TG4PhysicsList::GetG3ParticleWSPName:\n";
746 text = text + " Wrong particleWSP.";
747 TG4Globals::Exception(text);