[u/mrichter/AliRoot.git] / TGeant4 / TG4PhysicsManager.cxx

// $Id$
// Category: physics
//
// Author: I. Hrivnacova
//
// Class TG4PhysicsManager
// -----------------------
// See the class description in the header file.

#include "TG4PhysicsManager.h"
#include "TG4ModularPhysicsList.h"
#include "TG4ParticlesManager.h"
#include "TG4G3PhysicsManager.h"
#include "TG4PhysicsConstructorGeneral.h"
#include "TG4PhysicsConstructorEM.h"
#include "TG4PhysicsConstructorMuon.h"
#include "TG4PhysicsConstructorHadron.h"
#include "TG4PhysicsConstructorIon.h"
#include "TG4PhysicsConstructorOptical.h"
#include "TG4PhysicsConstructorSpecialCuts.h"
#include "TG4PhysicsConstructorSpecialControls.h"
#include "TG4GeometryServices.h"
#include "TG4G3Cut.h"
#include "TG4G3Control.h"
#include "TG4G3Units.h"
#include "TG4Limits.h"
#include "AliDecayer.h"

#include <G4ParticleDefinition.hh>
#include <G4OpBoundaryProcess.hh>
#include <G4VProcess.hh>
#include <G3MedTable.hh>

#include <TDatabasePDG.h>

TG4PhysicsManager* TG4PhysicsManager::fgInstance = 0;

//_____________________________________________________________________________
TG4PhysicsManager::TG4PhysicsManager(TG4ModularPhysicsList* physicsList)
  : fPhysicsList(physicsList),
    fDecayer(0),
    fSetEMPhysics(true),
    fSetMuonPhysics(true),
    fSetHadronPhysics(false),
    fSetOpticalPhysics(false),
    fSetSpecialCutsPhysics(false),
    fSetSpecialControlsPhysics(false)

{ 
//
  if (fgInstance) {
    TG4Globals::Exception(
      "TG4PhysicsManager: attempt to create two instances of singleton.");
  }
      
  fgInstance = this;  
  
  // create particles manager
  fParticlesManager = new TG4ParticlesManager();
      
  // create G3 physics manager
  fG3PhysicsManager = new TG4G3PhysicsManager();

  // fill process name map
  // FillProcessMap();
}

//_____________________________________________________________________________
TG4PhysicsManager::TG4PhysicsManager(){
//
  delete fDecayer;
  delete fParticlesManager;
  delete fG3PhysicsManager;
}

//_____________________________________________________________________________
TG4PhysicsManager::TG4PhysicsManager(const TG4PhysicsManager& right) {
// 
  TG4Globals::Exception(
    "Attempt to copy TG4PhysicsManager singleton.");
}

//_____________________________________________________________________________
TG4PhysicsManager::~TG4PhysicsManager() {
//
}

// operators

//_____________________________________________________________________________
TG4PhysicsManager& 
TG4PhysicsManager::operator=(const TG4PhysicsManager& right)
{
  // check assignement to self
  if (this == &right) return *this;

  TG4Globals::Exception(
    "Attempt to assign TG4PhysicsManager singleton.");
    
  return *this;  
}    
          
// private methods

//_____________________________________________________________________________
void TG4PhysicsManager::FillProcessMap()
{
// Fills fProcessMCMap.
// The default G4 process names are used in the map.
// Not used - the map is filled in physics constructors
// only with processes that are really built.
// ---

  // multiple scattering
  fProcessMCMap.Add("msc",  kPMultipleScattering);
  fProcessMCMap.Add("Imsc", kPMultipleScattering);
    
  // continuous energy loss
  // !! including delta rays
  fProcessMCMap.Add("eIoni",  kPEnergyLoss);
  fProcessMCMap.Add("IeIoni", kPEnergyLoss);
  fProcessMCMap.Add("LowEnergyIoni", kPEnergyLoss);
  fProcessMCMap.Add("hIoni",  kPEnergyLoss);
  fProcessMCMap.Add("IhIoni", kPEnergyLoss);
  fProcessMCMap.Add("hLowEIoni", kPEnergyLoss);
  fProcessMCMap.Add("MuIoni", kPEnergyLoss);
  fProcessMCMap.Add("IMuIonisation", kPEnergyLoss);
  fProcessMCMap.Add("ionIoni",  kPEnergyLoss);
  fProcessMCMap.Add("ionLowEIoni",  kPEnergyLoss);
  fProcessMCMap.Add("PAIonisation",  kPEnergyLoss);
  
  // bending in mag. field
  // kPMagneticFieldL

  // particle decay
  fProcessMCMap.Add("Decay", kPDecay);
  
  // photon pair production or
  // muon direct pair production
  fProcessMCMap.Add("conv", kPPair);
  fProcessMCMap.Add("LowEnConversion", kPPair);
  fProcessMCMap.Add("MuPairProd", kPPair);
  fProcessMCMap.Add("IMuPairProduction", kPPair);

  // Compton scattering
  fProcessMCMap.Add("compt", kPCompton);
  fProcessMCMap.Add("LowEnCompton", kPCompton);
  fProcessMCMap.Add("polarCompt", kPCompton);

  // photoelectric effect
  fProcessMCMap.Add("phot", kPPhotoelectric);
  fProcessMCMap.Add("LowEnPhotoElec", kPPhotoelectric);

  // bremsstrahlung
  fProcessMCMap.Add("eBrem", kPBrem);
  fProcessMCMap.Add("IeBrem", kPBrem);
  fProcessMCMap.Add("MuBrems", kPBrem);
  fProcessMCMap.Add("IMuBremsstrahlung", kPBrem);
  fProcessMCMap.Add("LowEnBrem", kPBrem);

  // delta-ray production
  // kPDeltaRay
  // has to be distinguished from kPEnergyLoss on flight
  
  // positron annihilation
  fProcessMCMap.Add("annihil", kPAnnihilation);
  fProcessMCMap.Add("Iannihil", kPAnnihilation);

  // hadronic interaction
  // kPHadronic

  // nuclear evaporation
  // kPEvaporation
  
  // nuclear fission
  // kPNuclearFission

  // nuclear absorption
  fProcessMCMap.Add("PionMinusAbsorptionAtRest", kPNuclearAbsorption);
  fProcessMCMap.Add("PiMinusAbsorptionAtRest", kPNuclearAbsorption);
  fProcessMCMap.Add("KaonMinusAbsorption", kPNuclearAbsorption);         
  fProcessMCMap.Add("KaonMinusAbsorptionAtRest", kPNuclearAbsorption);         
  
  // antiproton annihilation
  fProcessMCMap.Add("AntiProtonAnnihilationAtRest", kPPbarAnnihilation);
  // fProcessMCMap.Add("AntiNeutronAnnihilationAtRest", not defined);

  // neutron capture    
  fProcessMCMap.Add("NeutronCaptureAtRest", kPNCapture);
  // fProcessMCMap.Add("LCapture", hadron capture not defined);

  // hadronic elastic incoherent scattering
  fProcessMCMap.Add("LElastic", kPHElastic);

  // hadronic inelastic scattering
  fProcessMCMap.Add("inelastic", kPHInhelastic);

  // muon nuclear interaction
  fProcessMCMap.Add("MuNucl", kPMuonNuclear);

  // exceeded time of flight cut
  // kPTOFlimit
  
  // nuclear photofission
  // kPPhotoFission

  // Rayleigh scattering
  fProcessMCMap.Add("Rayleigh Scattering", kPRayleigh);

  // no mechanism is active, usually at the entrance of a new volume
  fProcessMCMap.Add("Transportation", kPNull);

  // particle has fallen below energy threshold and tracking stops
  // kPStop
  
  // Cerenkov photon absorption
  fProcessMCMap.Add("Absorption", kPLightAbsorption);

  // Cerenkov photon reflection/refraction
  // kPLightScattering, kPLightReflection, kPLightRefraction
  // has to be inquired from the G4OpBoundary process

  // synchrotron radiation
  fProcessMCMap.Add("SynchrotronRadiation", kPSynchrotron);
}  

//_____________________________________________________________________________
void TG4PhysicsManager::GstparCut(G4int itmed, TG4G3Cut par, G4double parval)
{
// Sets special tracking medium parameter. 
// It is applied to all logical volumes that use the specified 
// tracking medium.
// ---

  // get medium from table
  G3MedTableEntry* medium = G3Med.get(itmed);
  if (!medium) {
    G4String text = "TG4PhysicsManager::GstparCut: \n";
    text = text + "    Medium not found."; 
    G4Exception(text);
  }  

  // get/create user limits
  TG4Limits* limits 
    = TG4GeometryServices::Instance()->GetLimits(medium->GetLimits());
    
  if (!limits) {
    limits = new TG4Limits(*fG3PhysicsManager->GetCutVector(),
                           *fG3PhysicsManager->GetControlVector());
    medium->SetLimits(limits);

    // add verbose 
    G4cout << "TG4PhysicsManager::GstparCut: new TG4Limits() for medium " 
           << itmed << " has been created." << G4endl;  
  }	   

  // add units
  if (par == kTOFMAX) parval *= TG4G3Units::Time();
  else                parval *= TG4G3Units::Energy();

  // set parameter
  limits->SetG3Cut(par, parval);
}


//_____________________________________________________________________________
void TG4PhysicsManager::GstparControl(G4int itmed, TG4G3Control par, 
                                      TG4G3ControlValue parval)
{
// Sets special tracking medium parameter. 
// It is applied to all logical volumes that use the specified 
// tracking medium.
// ---

  // get medium from table
  G3MedTableEntry* medium = G3Med.get(itmed);
  if (!medium) {
    G4String text = "TG4PhysicsManager::GstparControl: \n";
    text = text + "    Medium not found."; 
    G4Exception(text);
  }  

  // get/create user limits
  TG4Limits* limits 
    = TG4GeometryServices::Instance()->GetLimits(medium->GetLimits());

  if (!limits) {
    limits = new TG4Limits(*fG3PhysicsManager->GetCutVector(),
                           *fG3PhysicsManager->GetControlVector());
    medium->SetLimits(limits);

    // add verbose 
    G4cout << "TG4PhysicsManager::GstparControl: new TG4Limits() for medium" 
           << itmed << " has been created." << G4endl;  
  }	   
  
  // set parameter
  limits->SetG3Control(par, parval);
}

// public methods

//_____________________________________________________________________________
void  TG4PhysicsManager::Gstpar(Int_t itmed, const char *param, Float_t parval) 
{ 
// Passes the tracking medium parameter to TG4Limits.
// The tracking medium parameter is set only in case
// its value is different from the "global" physics setup.
// (If: CheckCut/ControlWithG3Defaults is used checking
//  is performed with respect to G3 default values.)
// When any cut/control parameter is set in limits
// the physics manager is locked and the physics setup
// cannot be changed.
// Applying this TG4Limits to particles and physical
// processes is still in development.
//
//  Geant3 desription:
//  ==================
//  To change the value of cut  or mechanism "CHPAR"
//  to a new value PARVAL  for tracking medium ITMED
//  The  data   structure  JTMED   contains  the   standard  tracking
//  parameters (CUTS and flags to control the physics processes)  which
//  are used  by default  for all  tracking media.   It is  possible to
//  redefine individually  with GSTPAR  any of  these parameters  for a
//  given tracking medium. 
//  ITMED     tracking medium number 
//  CHPAR     is a character string (variable name) 
//  PARVAL    must be given as a floating point.
// ---

  G4String name = TG4GeometryServices::Instance()->CutName(param); 
  TG4G3Cut cut;
  if (fG3PhysicsManager->CheckCutWithTheVector(name, parval, cut)) {
      GstparCut(itmed, cut, parval);
      fG3PhysicsManager->Lock();
  }  
  else {
    TG4G3Control control;
    TG4G3ControlValue controlValue; 
    if (fG3PhysicsManager
         ->CheckControlWithTheVector(name, parval, control, controlValue)) {
      GstparControl(itmed, control, controlValue);
      fG3PhysicsManager->Lock();
    } 
    else if (cut==kNoG3Cuts && control==kNoG3Controls) { 
      G4String text = "TG4PhysicsManager::Gstpar:";
      text = text + name;
      text = text + " parameter is not yet implemented.";
      TG4Globals::Warning(text);
    }	
  }   
} 
 
//_____________________________________________________________________________
void TG4PhysicsManager::CreatePhysicsConstructors()
{
// Creates the selected physics constructors
// and registeres them in the modular physics list.
// ---

  // general physics
  fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorGeneral());

  // electromagnetic physics
  if (fSetEMPhysics) 
    fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorEM());

  // muon physics
  if (fSetMuonPhysics) 
    fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorMuon());

  // hadron physics
  if (fSetHadronPhysics) {
    fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorIon());
    fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorHadron());
  }  

  // optical physics
  if (fSetOpticalPhysics) 
    fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorOptical());

  if (fSetSpecialCutsPhysics) 
    fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorSpecialCuts());

  if (fSetSpecialControlsPhysics) 
    fPhysicsList->RegisterPhysics(new TG4PhysicsConstructorSpecialControls());

         // all created physics constructors are deleted
	 // in the TG4ModularPhysicsList destructor
}    

//_____________________________________________________________________________
void TG4PhysicsManager::SetCut(const char* cutName, Float_t cutValue)
{
// Sets the specified cut.
// ---

  fG3PhysicsManager->CheckLock();
  TG4G3Cut g3Cut = TG4G3CutVector::GetCut(cutName);

  if (g3Cut == kNoG3Cuts) {
    G4String text = "TG4PhysicsManager::SetCut:\n";
    text = text + "    Parameter " + cutName;
    text = text + " is not implemented.";
    TG4Globals::Warning(text);
    return;
  }  
  
  // add units
  if (g3Cut == kTOFMAX) cutValue *= TG4G3Units::Time();
  else                  cutValue *= TG4G3Units::Energy();

  fG3PhysicsManager->SetCut(g3Cut, cutValue);    
}  
  
//_____________________________________________________________________________
void TG4PhysicsManager::SetProcess(const char* controlName, Int_t value)
{
// Sets the specified process control.
// ---

  fG3PhysicsManager->CheckLock();
  TG4G3Control control = TG4G3ControlVector::GetControl(controlName);
  
  if (control != kNoG3Controls) {
    TG4G3ControlValue controlValue 
      = TG4G3ControlVector::GetControlValue(value, control);
    fG3PhysicsManager->SetProcess(control, controlValue);
  }  
  else {   
    G4String text = "TG4PhysicsManager::SetProcess:\n";
    text = text + "    Parameter " + controlName;
    text = text + " is not implemented.";
    TG4Globals::Warning(text);
  }  
}  

//_____________________________________________________________________________
Float_t TG4PhysicsManager::Xsec(char* ch, Float_t p1, Int_t i1, Int_t i2)
{
// Not yet implemented -> gives exception.
// ---

  TG4Globals::Exception(
    "TG4PhysicsManager::Xsec: not yet implemented.");

  return 0.;
}    
  
//_____________________________________________________________________________
Int_t TG4PhysicsManager::IdFromPDG(Int_t pdgID) const
{
// G4 does not use the integer particle identifiers
// Id <-> PDG is identity.
// ---

  return pdgID;
}  

//_____________________________________________________________________________
Int_t TG4PhysicsManager::PDGFromId(Int_t mcID) const
{
// G4 does not use integer particle identifiers
// Id <-> PDG is identity.
// ---

  return mcID;
}  

//_____________________________________________________________________________
void  TG4PhysicsManager::DefineParticles()
{
  // ======
  // Taken from TGeant3
  //
  // Use ENDF-6 mapping for ions = 10000*z+10*a+iso
  // and add 1 000 000
  // and numbers above 5 000 000 for special applications
  //

  const Int_t kion=10000000;
  const Int_t kspe=50000000;

  const Double_t kGeV=0.9314943228;
  const Double_t kHslash = 1.0545726663e-27;
  const Double_t kErgGeV = 1/1.6021773349e-3;
  const Double_t kHshGeV = kHslash*kErgGeV;
  const Double_t kYearsToSec = 3600*24*365.25;

  TDatabasePDG *pdgDB = TDatabasePDG::Instance();

  pdgDB->AddParticle("Deuteron","Deuteron",2*kGeV+8.071e-3,kTRUE,
		     0,1,"Ion",kion+10020);
		     
  pdgDB->AddParticle("Triton","Triton",3*kGeV+14.931e-3,kFALSE,
		     kHshGeV/(12.33*kYearsToSec),1,"Ion",kion+10030);

  pdgDB->AddParticle("Alpha","Alpha",4*kGeV+2.424e-3,kTRUE,
		     kHshGeV/(12.33*kYearsToSec),2,"Ion",kion+20040);

  pdgDB->AddParticle("HE3","HE3",3*kGeV+14.931e-3,kFALSE,
		     0,2,"Ion",kion+20030);

  pdgDB->AddParticle("Cherenkov","Cherenkov",0,kFALSE,
		     0,0,"Special",kspe+50);

  pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,
		     0,0,"Special",kspe+51);


  // To do: define the PDG database extension
  // in a common part.
  //
  // AliMC::ExtendPDGDatabase(); 
  //
  // end of "common" implementation
  // ======

  fParticlesManager->MapParticles();
}    

//_____________________________________________________________________________
void TG4PhysicsManager::SetProcessActivation()
{
// (In)Activates built processes according
// to the setup in TG4G3PhysicsManager::fControlVector.
// ---

  fPhysicsList->SetProcessActivation();
}       


//_____________________________________________________________________________
AliMCProcess TG4PhysicsManager::GetMCProcess(const G4VProcess* process)
{
// Returns the AliMCProcess code of the specified G4 process.
// ---
 
  if (!process) return kPNoProcess;

  return fProcessMCMap.GetMCProcess(process);
}

//_____________________________________________________________________________
AliMCProcess TG4PhysicsManager::GetOpBoundaryStatus(const G4VProcess* process)
{
// Returns the AliMCProcess code according to the OpBoundary process
// status.
// ---
 
  if (!process) return kPNoProcess;

#ifdef TGEANT4_DEBUG
  G4OpBoundaryProcess* opBoundary
    = dynamic_cast<G4OpBoundaryProcess*>(process);
    
  if (!opBoundary) 
    TG4Globals::Exception(
      "TG4PhysicsManager::GetOpBoundaryStatus: Wrong process type.");
    return kPNoProcess;
  }
  
  return opBoundary;  
#else
  G4OpBoundaryProcess* opBoundary = (G4OpBoundaryProcess*)process;
#endif  

  switch (opBoundary->GetStatus()) {
    // reflection
    case FresnelReflection: 
    case TotalInternalReflection:
    case LambertianReflection: 
    case LobeReflection:
    case SpikeReflection: 
    case BackScattering:
       return kPLightReflection;
       ;;

    // refraction
    case FresnelRefraction: 
       return kPLightReflection;
       ;;

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