major upgrade: separated TG4ParticlesManager, TG4G3PhysicsManager (now aggregated)
[u/mrichter/AliRoot.git] / 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 "TG4ParticlesManager.h"
8 #include "TG4G3PhysicsManager.h"
9 #include "TG4PhysicsConstructorEM.h"
10 #include "TG4PhysicsConstructorOptical.h"
11 #include "TG4PhysicsConstructorHadron.h"
12 #include "TG4PhysicsConstructorSpecialCuts.h"
13 #include "TG4PhysicsConstructorSpecialControls.h"
14 #include "TG4G3Cut.h"
15 #include "TG4G3Control.h"
16
17 #include <G4ParticleDefinition.hh>
18 #include <G4VProcess.hh>
19 #include <G4VModularPhysicsList.hh>
20
21 #include <TDatabasePDG.h>
22
23 TG4PhysicsManager* TG4PhysicsManager::fgInstance = 0;
24
25 TG4PhysicsManager::TG4PhysicsManager(G4VModularPhysicsList* physicsList)
26   : fPhysicsList(physicsList),
27     fSetEMPhysics(true),
28     fSetOpticalPhysics(false),
29     fSetHadronPhysics(false),
30     fSetSpecialCutsPhysics(false),
31     fSetSpecialControlsPhysics(false)
32
33
34 //
35   if (fgInstance) {
36     TG4Globals::Exception(
37       "TG4PhysicsManager: attempt to create two instances of singleton.");
38   }
39       
40   fgInstance = this;  
41   
42   // create particles manager
43   fParticlesManager = new TG4ParticlesManager();
44       
45   // create G3 physics manager
46   fG3PhysicsManager = new TG4G3PhysicsManager();
47
48   // fill process name map
49   FillProcessMap();
50 }
51
52 TG4PhysicsManager::TG4PhysicsManager(){
53 //
54   delete fParticlesManager;
55   delete fG3PhysicsManager;
56 }
57
58 TG4PhysicsManager::TG4PhysicsManager(const TG4PhysicsManager& right) {
59 // 
60   TG4Globals::Exception(
61     "Attempt to copy TG4PhysicsManager singleton.");
62 }
63
64 TG4PhysicsManager::~TG4PhysicsManager() {
65 //
66 }
67
68 // operators
69
70 TG4PhysicsManager& 
71 TG4PhysicsManager::operator=(const TG4PhysicsManager& right)
72 {
73   // check assignement to self
74   if (this == &right) return *this;
75
76   TG4Globals::Exception(
77     "Attempt to assign TG4PhysicsManager singleton.");
78     
79   return *this;  
80 }    
81           
82 // private methods
83
84 void TG4PhysicsManager::FillProcessMap()
85 {
86 // Fills fProcessMap.
87 // The default G4 process names are used in the map.
88 // ---
89
90   // multiple scattering
91   fProcessMap.Add("msc",  kPMultipleScattering);
92   fProcessMap.Add("Imsc", kPMultipleScattering);
93     
94   // continuous energy loss
95   // !! including delta rays
96   fProcessMap.Add("eIoni",  kPEnergyLoss);
97   fProcessMap.Add("IeIoni", kPEnergyLoss);
98   fProcessMap.Add("LowEnergyIoni", kPEnergyLoss);
99   fProcessMap.Add("hIoni",  kPEnergyLoss);
100   fProcessMap.Add("IhIoni", kPEnergyLoss);
101   fProcessMap.Add("hLowEIoni", kPEnergyLoss);
102   fProcessMap.Add("MuIoni", kPEnergyLoss);
103   fProcessMap.Add("IMuIonisation", kPEnergyLoss);
104   fProcessMap.Add("ionIoni",  kPEnergyLoss);
105   fProcessMap.Add("ionLowEIoni",  kPEnergyLoss);
106   fProcessMap.Add("PAIonisation",  kPEnergyLoss);
107   
108   // bending in mag. field
109   // kPMagneticFieldL
110
111   // particle decay
112   fProcessMap.Add("Decay", kPDecay);
113   
114   // photon pair production or
115   // muon direct pair production
116   fProcessMap.Add("conv", kPPair);
117   fProcessMap.Add("LowEnConversion", kPPair);
118   fProcessMap.Add("MuPairProd", kPPair);
119   fProcessMap.Add("IMuPairProduction", kPPair);
120
121   // Compton scattering
122   fProcessMap.Add("compt", kPCompton);
123   fProcessMap.Add("LowEnCompton", kPCompton);
124   fProcessMap.Add("polarCompt", kPCompton);
125
126   // photoelectric effect
127   fProcessMap.Add("phot", kPPhotoelectric);
128   fProcessMap.Add("LowEnPhotoElec", kPPhotoelectric);
129
130   // bremsstrahlung
131   fProcessMap.Add("eBrem", kPBrem);
132   fProcessMap.Add("IeBrem", kPBrem);
133   fProcessMap.Add("MuBrem", kPBrem);
134   fProcessMap.Add("IMuBremsstrahlung", kPBrem);
135   fProcessMap.Add("LowEnBrem", kPBrem);
136
137   // delta-ray production
138   // kPDeltaRay
139   // has to be distinguished from kPEnergyLoss on flight
140   
141   // positron annihilation
142   fProcessMap.Add("annihil", kPAnnihilation);
143   fProcessMap.Add("Iannihil", kPAnnihilation);
144
145   // hadronic interaction
146   // kPHadronic
147
148   // nuclear evaporation
149   // kPEvaporation
150   
151   // nuclear fission
152   // kPNuclearFission
153
154   // nuclear absorption
155   fProcessMap.Add("PionMinusAbsorptionAtRest", kPNuclearAbsorption);
156   fProcessMap.Add("PiMinusAbsorptionAtRest", kPNuclearAbsorption);
157   fProcessMap.Add("KaonMinusAbsorption", kPNuclearAbsorption);         
158   fProcessMap.Add("KaonMinusAbsorptionAtRest", kPNuclearAbsorption);         
159   
160   // antiproton annihilation
161   fProcessMap.Add("AntiProtonAnnihilationAtRest", kPPbarAnnihilation);
162   // fProcessMap.Add("AntiNeutronAnnihilationAtRest", not defined);
163
164   // neutron capture    
165   fProcessMap.Add("NeutronCaptureAtRest", kPNCapture);
166   // fProcessMap.Add("LCapture", hadron capture not defined);
167
168   // hadronic elastic incoherent scattering
169   fProcessMap.Add("LElastic", kPHElastic);
170
171   // hadronic inelastic scattering
172   fProcessMap.Add("inelastic", kPHInhelastic);
173
174   // muon nuclear interaction
175   fProcessMap.Add("MuNucl", kPMuonNuclear);
176
177   // exceeded time of flight cut
178   // kPTOFlimit
179   
180   // nuclear photofission
181   // kPPhotoFission
182
183   // Rayleigh scattering
184   fProcessMap.Add("Rayleigh Scattering", kPRayleigh);
185
186   // no mechanism is active, usually at the entrance of a new volume
187   fProcessMap.Add("Transportation", kPNull);
188
189   // particle has fallen below energy threshold and tracking stops
190   // kPStop
191   
192   // Cerenkov photon absorption
193   fProcessMap.Add("Absorption", kPLightAbsorption);
194
195   // Cerenkov photon reflection/refraction
196   // kPLightScattering, kPLightReflection, kPLightRefraction
197   // has to be inquired from the G4OpBoundary process
198
199   // synchrotron radiation
200   fProcessMap.Add("SynchrotronRadiation", kPSynchrotron);
201 }  
202
203
204 // public methods
205
206 void TG4PhysicsManager::BuildPhysics()
207 {
208 // Empty function - not needed in G4.
209 // (Physics is built within /run/initialize.)
210 // ---
211
212   TG4Globals::Warning(
213     "TG4PhysicsManager::BuildPhysics: is empty function in G4 MC.");
214 }    
215
216 void TG4PhysicsManager::CreatePhysicsConstructors()
217 {
218 // Creates the selected physics constructors
219 // and registeres them in the modular physics list.
220 // ---
221
222   // electromagnetic physics
223   if (fSetEMPhysics) 
224     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorEM());
225
226   // optical physics
227   if (fSetOpticalPhysics) 
228     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorOptical());
229
230   // hadron physics
231   if (fSetHadronPhysics) 
232     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorHadron());
233
234   if (fSetSpecialCutsPhysics) 
235     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorSpecialCuts());
236
237   if (fSetSpecialControlsPhysics) 
238     fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorSpecialControls());
239
240          // all created physics constructors are deleted
241          // in the G4VModularPhysicsList destructor
242 }    
243
244 void TG4PhysicsManager::SetCut(const char* cutName, Float_t cutValue)
245 {
246 // Sets the specified cut.
247 // ---
248
249   fG3PhysicsManager->CheckLock();
250   TG4G3Cut g3Cut = fG3PhysicsManager->GetG3Cut(cutName);
251   if (g3Cut != kNoG3Cuts)
252     fG3PhysicsManager->SetCut(g3Cut, cutValue);
253   else {   
254     G4String text = "TG4PhysicsManager::SetCut:\n";
255     text = text + "    Parameter " + cutName;
256     text = text + " is not implemented.";
257     TG4Globals::Warning(text);
258   }  
259 }  
260   
261 void TG4PhysicsManager::SetProcess(const char* controlName, Int_t controlValue)
262 {
263 // Sets the specified process control.
264 // ---
265
266   fG3PhysicsManager->CheckLock();
267   TG4G3Control control = fG3PhysicsManager->GetG3Control(controlName);
268   if (control != kNoG3Controls)
269     fG3PhysicsManager->SetProcess(control, controlValue);
270   else {   
271     G4String text = "TG4PhysicsManager::SetProcess:\n";
272     text = text + "    Parameter " + controlName;
273     text = text + " is not implemented.";
274     TG4Globals::Warning(text);
275   }  
276 }  
277
278 Float_t TG4PhysicsManager::Xsec(char* ch, Float_t p1, Int_t i1, Int_t i2)
279 {
280 // Not yet implemented -> gives exception.
281 // ---
282
283   TG4Globals::Exception(
284     "TG4PhysicsManager::Xsec: not yet implemented.");
285
286   return 0.;
287 }    
288   
289 Int_t TG4PhysicsManager::IdFromPDG(Int_t pdgID) const
290 {
291 // G4 does not use the integer particle identifiers
292 // Id <-> PDG is identity.
293 // ---
294
295   return pdgID;
296 }  
297
298 Int_t TG4PhysicsManager::PDGFromId(Int_t mcID) const
299 {
300 // G4 does not use integer particle identifiers
301 // Id <-> PDG is identity.
302 // ---
303
304   return mcID;
305 }  
306
307 void  TG4PhysicsManager::DefineParticles()
308 {
309   // ======
310   // Taken from TGeant3
311   //
312   // Use ENDF-6 mapping for ions = 10000*z+10*a+iso
313   // and add 1 000 000
314   // and numbers above 5 000 000 for special applications
315   //
316
317   const Int_t kion=10000000;
318   const Int_t kspe=50000000;
319
320   const Double_t kGeV=0.9314943228;
321   const Double_t kHslash = 1.0545726663e-27;
322   const Double_t kErgGeV = 1/1.6021773349e-3;
323   const Double_t kHshGeV = kHslash*kErgGeV;
324   const Double_t kYearsToSec = 3600*24*365.25;
325
326   TDatabasePDG *pdgDB = TDatabasePDG::Instance();
327
328   pdgDB->AddParticle("Deuteron","Deuteron",2*kGeV+8.071e-3,kTRUE,
329                      0,1,"Ion",kion+10020);
330                      
331   pdgDB->AddParticle("Triton","Triton",3*kGeV+14.931e-3,kFALSE,
332                      kHshGeV/(12.33*kYearsToSec),1,"Ion",kion+10030);
333
334   pdgDB->AddParticle("Alpha","Alpha",4*kGeV+2.424e-3,kTRUE,
335                      kHshGeV/(12.33*kYearsToSec),2,"Ion",kion+20040);
336
337   pdgDB->AddParticle("HE3","HE3",3*kGeV+14.931e-3,kFALSE,
338                      0,2,"Ion",kion+20030);
339
340   pdgDB->AddParticle("Cherenkov","Cherenkov",0,kFALSE,
341                      0,0,"Special",kspe+50);
342
343   pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,
344                      0,0,"Special",kspe+51);
345
346
347   // To do: define the PDG database extension
348   // in a common part.
349   //
350   // AliMC::ExtendPDGDatabase(); 
351   //
352   // end of "common" implementation
353   // ======
354
355   fParticlesManager->MapParticles();
356 }    
357
358
359 void TG4PhysicsManager::SetProcessActivation()
360 {
361 // (In)Activates built processes according
362 // to the setup in fControlVector.
363 // ---
364
365   if (fPhysicsList) {
366     // temporarily excluded
367     // fPhysicsList->SetProcessActivation();
368   }  
369   else {
370     G4String text = "TG4PhysicsManager::SetProcessActivation:\n";
371     text = text +   "   There is no physics list set.";
372     TG4Globals::Exception(text);
373   }
374 }       
375
376
377 AliMCProcess TG4PhysicsManager::GetMCProcess(const G4VProcess* process)
378 {
379 // Returns the AliMCProcess code of the specified G4 process.
380 // ---
381  
382   if (!process) return kPNoProcess;
383
384   G4String name = process->GetProcessName();
385   G4int code = fProcessMap.GetSecond(name);
386   
387   if (code == 0) return kPNoProcess;
388   
389   return (AliMCProcess)code; 
390 }
391