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

// $Id$
// Category: physics
//
// See the class description in the header file.

#include "TG4PhysicsManager.h"
#include "TG4PhysicsList.h"
#include "TG4CutVector.h"
#include "TG4FlagVector.h"
#include "TG4G3Defaults.h"

#include <G4ParticleDefinition.hh>
#include <Randomize.hh>

#include <TDatabasePDG.h>

TG4PhysicsManager* TG4PhysicsManager::fgInstance = 0;

TG4PhysicsManager::TG4PhysicsManager()
  : fLock(false),
    fPhysicsList(0),
    fCutVector(0),
    fFlagVector(0) 
{ 
//
  if (fgInstance) {
    TG4Globals::Exception(
      "TG4PhysicsManager: attempt to create two instances of singleton.");
  }
      
  fgInstance = this;  
      
  // initialize fIsCutVector 
  fIsCutVector = new TG4boolVector;
  G4int i;
  //for (i=0; i<kNofParticlesWSP; i++) fIsCutVector->insert(false);  
  for (i=0; i<kNofParticlesWSP; i++) fIsCutVector->push_back(false);

  // initialize fIsFlagVector 
  fIsFlagVector = new TG4boolVector;
  //for (i=0; i<kNofParticlesWSP; i++) fIsFlagVector->insert(false);
  for (i=0; i<kNofParticlesWSP; i++) fIsFlagVector->push_back(false);

  // define fCutNameVector
  fG3CutNameVector.insert("CUTGAM");
  fG3CutNameVector.insert("CUTELE");
  fG3CutNameVector.insert("CUTNEU");
  fG3CutNameVector.insert("CUTHAD");
  fG3CutNameVector.insert("CUTMUO");
  fG3CutNameVector.insert("BCUTE");
  fG3CutNameVector.insert("BCUTM"); 
  fG3CutNameVector.insert("DCUTE");
  fG3CutNameVector.insert("DCUTM");
  fG3CutNameVector.insert("PPCUTM");

  // define fFlagNameVector
  fG3FlagNameVector.insert("PAIR");
  fG3FlagNameVector.insert("COMP");
  fG3FlagNameVector.insert("PHOT");
  fG3FlagNameVector.insert("PFIS");
  fG3FlagNameVector.insert("DRAY");
  fG3FlagNameVector.insert("ANNI");
  fG3FlagNameVector.insert("BREM");
  fG3FlagNameVector.insert("HADR");
  fG3FlagNameVector.insert("MUNU");
  fG3FlagNameVector.insert("DCAY");
  fG3FlagNameVector.insert("LOSS");
  fG3FlagNameVector.insert("MULS");
}

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

TG4PhysicsManager::~TG4PhysicsManager() {
//
  delete fIsCutVector;
  delete fIsFlagVector;
}

// 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::LockException() const
{
// Gives exception in case of attempt to modified physics
// setup after physics manager was locked.
// ---

  G4String text = "TG4PhysicsManager: \n";
  text = text + "    It is too late to change physics setup. \n";
  text = text + "    PhysicsManager has been already locked.";
  TG4Globals::Exception(text);
}

G4int TG4PhysicsManager::GetPDGEncoding(G4ParticleDefinition* particle)
{
// Returns the PDG code of particle;
// if standard PDG code is not defined the TDatabasePDG
// is used.
// ---

  // get PDG encoding from G4 particle definition
  G4int pdgEncoding = particle->GetPDGEncoding();

  if (pdgEncoding == 0) {
    // get PDG encoding from TDatabasePDG
  
    // get particle name from the name map
    G4String g4name = particle->GetParticleName();
    G4String tname = fParticleNameMap.GetSecond(g4name);
    if (tname == "Undefined") {
      G4String text = "TG4PhysicsManager::GetPDGEncoding: \n";
      text = text + "    Particle " + g4name;
      text = text + " was not found in the name map.";
      TG4Globals::Exception(text);
    }  
  
    // get particle from TDatabasePDG
    TDatabasePDG* pdgDB = TDatabasePDG::Instance();
    TParticlePDG* particle = pdgDB->GetParticle(tname);
    if (!particle) {
      G4String text = "TG4PhysicsManager::GetPDGEncoding: \n";
      text = text + "    Particle " + tname;
      text = text + " was not found in TDatabasePDG.";
      TG4Globals::Exception(text);
    }  
    
    // get PDG encoding
    pdgEncoding = particle->PdgCode();
  }
    
  return pdgEncoding;  
}  
     
G4int TG4PhysicsManager::GetPDGEncoding(G4String particleName)
{
// Returns the PDG code of particle sepcified by name.
// ---

  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

  G4ParticleDefinition* particle = 0;  
  particle = particleTable->FindParticle(particleName);
  if (!particle) {    
    G4String text = "TG4PhysicsManager::GetPDGEncoding:\n";
    text = text + "   G4ParticleTable::FindParticle() " + particleName;
    text = text + " failed.";
    TG4Globals::Exception(text);
  }	

  return GetPDGEncoding(particle);
}  
  
void TG4PhysicsManager::SetCut(TG3Cut cut, G4double cutValue)
{  
// Sets kinetic energy cut (in a G3-like way).
// ---

  if (!fCutVector) {
    // create vector of kinetic energy cut values  
    fCutVector = new TG4CutVector();
  }  
  fCutVector->SetG3Cut(cut, cutValue);
  SwitchIsCutVector(cut);
}  

void TG4PhysicsManager::SetProcess(TG3Flag flag, G4int flagValue)
{
// Sets control process flag (in a G3-like way).
// ---

  if (!fFlagVector) {
    // create vector of control process flag values
    fFlagVector = new TG4FlagVector;
  }  
  fFlagVector->SetG3Flag(flag, flagValue);
} 


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
  // ======

  // map G4 particle names to TDatabasePDG names
  // (the map is built only for particles that have not
  //  defined standard PDG encoding)
  
  fParticleNameMap.Add("deuteron","Deuteron");
  fParticleNameMap.Add("triton",  "Triton");
  fParticleNameMap.Add("alpha",   "Alpha");
  fParticleNameMap.Add("He3",     "HE3");
  fParticleNameMap.Add("opticalphoton","Cherenkov");
  // fParticleNameMap.Add("???","FeedbackPhoton");
  fParticleNameMap.Add("geantino", "Rootino");
  
  // map G4 particle names to TDatabasePDG encodings
  fParticlePDGMap.Add("deuteron", GetPDGEncoding("deuteron"));
  fParticlePDGMap.Add("triton", GetPDGEncoding("triton"));
  fParticlePDGMap.Add("alpha", GetPDGEncoding("alpha"));
  fParticlePDGMap.Add("He3", GetPDGEncoding("He3") );
  fParticlePDGMap.Add("opticalphoton", GetPDGEncoding("opticalphoton"));
  // fParticlePDGMap.Add("???","FeedbackPhoton");
  fParticleNameMap.Add("geantino", GetPDGEncoding("geantino"));

  // add verbose
  G4cout << "Particle maps have been filled." << endl;
  //fParticleNameMap.PrintAll();
  //fParticlePDGMap.PrintAll();
}    

void TG4PhysicsManager::SwitchIsCutVector(TG3Cut cut)
{
// Updates the vector of booleans (fIsCutVector) for the specified cut.
// ---

  switch (cut) {
    case kCUTGAM: 
           (*fIsCutVector)[kGamma] = true; 
           break;
    case kBCUTE:
           (*fIsCutVector)[kGamma] = true; 
           break;
    case kBCUTM:
           (*fIsCutVector)[kGamma] = true; 
           break;
    case kCUTELE:
           (*fIsCutVector)[kElectron] = true; 
           break;
    case kDCUTE:
           (*fIsCutVector)[kElectron] = true; 
           break;
    case kDCUTM:
           (*fIsCutVector)[kElectron] = true; 
           break;
    case kCUTNEU:
           (*fIsCutVector)[kNeutralHadron] = true; 
           break;
    case kCUTHAD:
           (*fIsCutVector)[kChargedHadron] = true; 
           break;
    case kCUTMUO:
           (*fIsCutVector)[kMuon] = true; 
           break;
    default:
           break;
  }
}

void TG4PhysicsManager::SwitchIsFlagVector(TG3Flag flag)
{
// Updates the vector of booleans (fIsFlagVector) for the specified flag.
// ---

  switch (flag) {
    case kPAIR: 
           // gamma
           (*fIsFlagVector)[kGamma] = true; 
           break;
    case kCOMP:
           // gamma
           (*fIsFlagVector)[kGamma] = true; 
	   break;
    case kPHOT:
           // gamma
           (*fIsFlagVector)[kGamma] = true; 
	   break;
    case kPFIS:
           // gamma
           (*fIsFlagVector)[kGamma] = true; 
	   break;
    case kDRAY:	
           // all charged particles
           (*fIsFlagVector)[kElectron] = true; 
           (*fIsFlagVector)[kEplus] = true; 
           (*fIsFlagVector)[kChargedHadron] = true; 
           (*fIsFlagVector)[kMuon] = true; 
	   break;
    case kANNI:
	   // e+ only
           (*fIsFlagVector)[kEplus] = true; 
	   break;
    case kBREM:
	   // e-/e+, muons
           (*fIsFlagVector)[kElectron] = true; 
           (*fIsFlagVector)[kEplus] = true; 
           (*fIsFlagVector)[kMuon] = true; 
	   break;
    case kHADR:
	   // hadrons
           (*fIsFlagVector)[kNeutralHadron] = true; 
           (*fIsFlagVector)[kChargedHadron] = true; 
	   break;
    case kMUNU:
           // muons
           (*fIsFlagVector)[kMuon] = true; 
	   break;
    case kDCAY:
           // any
           (*fIsFlagVector)[kAny] = true; 
	   break;
    case kLOSS:
           // all charged particles
           (*fIsFlagVector)[kElectron] = true; 
           (*fIsFlagVector)[kEplus] = true; 
           (*fIsFlagVector)[kChargedHadron] = true; 
           (*fIsFlagVector)[kMuon] = true; 
	   break;
    case kMULS:
           // all charged particles
           (*fIsFlagVector)[kElectron] = true; 
           (*fIsFlagVector)[kEplus] = true; 
           (*fIsFlagVector)[kChargedHadron] = true; 
           (*fIsFlagVector)[kMuon] = true; 
	   break;
    default:
          break;
  }
}

TG3Cut TG4PhysicsManager::GetG3Cut(G4String cutName)
{
// Retrieves corresponding TG3Cut constant from the cutName.
// ---

  if      (cutName == fG3CutNameVector[kCUTGAM]) return kCUTGAM; 
  else if (cutName == fG3CutNameVector[kBCUTE])  return kBCUTE; 
  else if (cutName == fG3CutNameVector[kBCUTM])  return kBCUTM; 
  else if (cutName == fG3CutNameVector[kCUTELE]) return kCUTELE;
  else if (cutName == fG3CutNameVector[kDCUTE])  return kDCUTE;
  else if (cutName == fG3CutNameVector[kDCUTM])  return kDCUTM;
  else if (cutName == fG3CutNameVector[kCUTNEU]) return kCUTNEU;
  else if (cutName == fG3CutNameVector[kCUTHAD]) return kCUTHAD;
  else if (cutName == fG3CutNameVector[kCUTMUO]) return kCUTMUO;
  else return kNoG3Cuts;
}

TG3Flag TG4PhysicsManager::GetG3Flag(G4String flagName)
{
// Retrieves corresponding TG3Flag constant from the flagName.
// ---

  if      (flagName == fG3FlagNameVector[kPAIR]) return kPAIR;
  else if (flagName == fG3FlagNameVector[kCOMP]) return kCOMP;
  else if (flagName == fG3FlagNameVector[kPHOT]) return kPHOT;
  else if (flagName == fG3FlagNameVector[kPFIS]) return kPFIS;
  else if (flagName == fG3FlagNameVector[kDRAY]) return kDRAY;
  else if (flagName == fG3FlagNameVector[kANNI]) return kANNI;
  else if (flagName == fG3FlagNameVector[kBREM]) return kBREM;
  else if (flagName == fG3FlagNameVector[kHADR]) return kHADR;
  else if (flagName == fG3FlagNameVector[kMUNU]) return kMUNU;
  else if (flagName == fG3FlagNameVector[kDCAY]) return kDCAY;
  else if (flagName == fG3FlagNameVector[kLOSS]) return kLOSS;
  else if (flagName == fG3FlagNameVector[kMULS]) return kMULS;
  else return kNoG3Flags;
}

// public methods

void TG4PhysicsManager::BuildPhysics()
{
// Empty function - not needed in G4.
// (Physics is built within /run/initialize.)

  TG4Globals::Warning(
    "TG4PhysicsManager::BuildPhysics: is empty function in G4 MC.");
}    

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

  if (fLock) LockException();
  TG3Cut g3Cut = GetG3Cut(cutName);
  if (g3Cut != kNoG3Cuts)
    SetCut(g3Cut, cutValue);
  else {   
    G4String text = "TG4PhysicsManager::SetCut:\n";
    text = text + "    Parameter " + cutName;
    text = text + " is not implemented.";
    TG4Globals::Warning(text);
  }  
}  
  
void TG4PhysicsManager::SetProcess(const char* flagName, Int_t flagValue)
{
// Sets the specified process control.
// ---

  if (fLock) LockException();
  TG3Flag g3Flag = GetG3Flag(flagName);
  if (g3Flag != kNoG3Flags)
    SetProcess(g3Flag, flagValue);
  else {   
    G4String text = "TG4PhysicsManager::SetProcess:\n";
    text = text + "    Parameter " + flagName;
    text = text + " is not implemented.";
    TG4Globals::Warning(text);
  }  
}  

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

  if (fPhysicsList) {
    // temporarily excluded
    // fPhysicsList->SetProcessActivation();
  }  
  else {
    G4String text = "TG4PhysicsManager::SetProcessActivation:\n";
    text = text +   "   There is no physics list set.";
    TG4Globals::Exception(text);
  }
}       

G4int TG4PhysicsManager::GetPDGEncodingFast(G4ParticleDefinition* particle)
{
// Returns the PDG code of particle;
// if standard PDG code is not defined the preregistred
// fParticlePDGMap is used.
// ---

  // get PDG encoding from G4 particle definition
  G4int pdgEncoding = particle->GetPDGEncoding();

  if (pdgEncoding == 0) {
    // use FParticlePDGMap if standard PDG code is not defined
    G4String name = particle->GetParticleName();
    pdgEncoding = fParticlePDGMap.GetSecond(name);
  }
    
  return pdgEncoding;  
}  
     
G4bool TG4PhysicsManager::CheckCutWithCutVector(G4String name, 
                             G4double value, TG3Cut& cut)
{
// Retrieves corresponding TG3Cut from the name and 
// in case the value is different from the value in cutVector
// sets true the value of the fIsCutVector element 
// corresponding to this cut and returns true; 
// returns false otherwise.
// ---

  // convert cut name -> TG3Cut
  cut = GetG3Cut(name);

  // set switch vector element only if the value
  // is different from the value in cutVector
  if (cut !=kNoG3Cuts) {
    // get tolerance from TG4G3Defaults in GeV
    G4double tolerance = TG4G3Defaults::CutTolerance()/GeV;
    if (!(fCutVector) || (abs(value - (*fCutVector)[cut]) > tolerance)) {
      SwitchIsCutVector(cut);      
      return true;
    }  
    else return false;  
  }   	   	         
  return false;
}

G4bool TG4PhysicsManager::CheckFlagWithFlagVector(G4String name, 
                              G4double value, TG3Flag& flag)
{
// Retrieves corresponding TG3Flag from the name and 
// in case the value is different from the value in flagVector
// sets true the value of the fIsFlagVector element 
// corresponding to this flag and returns true; 
// returns false otherwise.
// ---

  // convert flag name -> TG3Flag
  flag = GetG3Flag(name);

  // set switch vector element only if the value
  // is different from the value in flagVector
  if (flag !=kNoG3Flags) {
    if (!(fFlagVector) || (abs(value - (*fFlagVector)[flag]) > 0.01)) {
      SwitchIsFlagVector(flag);      
      return true;
    }  
    else return false;  
  }   	   	         
  return false;
}

G4bool TG4PhysicsManager::CheckCutWithG3Defaults(G4String name, 
                              G4double value, TG3Cut& cut)
{
// Retrieves corresponding TG3Cut from the name and 
// in case the value is different from the G3 default value
// sets true the value of the SwitchCutVector element 
// corresponding to this cut and returns true; 
// returns false otherwise.
// ---

  // convert cut name -> TG3Cut
  cut = GetG3Cut(name);

  // set switch vector element only if the value
  // is different from G3 default
  if (cut !=kNoG3Cuts) {
    if (!TG4G3Defaults::IsDefaultCut(cut, value)) {
      SwitchIsCutVector(cut);      
      return true;
    }  
    else return false;  
  }   	   	         
  return false;
}

G4bool TG4PhysicsManager::CheckFlagWithG3Defaults(G4String name, 
                              G4double value, TG3Flag& flag)
{
// Retrieves corresponding TG3Flag from the name and 
// in case the value is different from the G3 default value
// sets true the value of the SwitchFlagVector element 
// corresponding to this flag and returns true; 
// returns false otherwise.
// ---

  // convert flag name -> TG3Flag
  flag = GetG3Flag(name);

  // set switch vector element only if the value
  // is different from G3 default
  if (flag !=kNoG3Flags) {
    if (!TG4G3Defaults::IsDefaultFlag(flag, value)) {
      SwitchIsFlagVector(flag);      
      return true;
    }  
    else return false;  
  }   	   	         
  return false;
}

void TG4PhysicsManager::SetG3DefaultCuts() 
{
// Sets G3 default values of kinetic energy cuts.
// ---

  if (fLock) LockException();
  if (!fCutVector) {
    // create vector of kinetic energy cut values  
    fCutVector = new TG4CutVector();
  }  
  fCutVector->SetG3Defaults();
}

void TG4PhysicsManager::SetG3DefaultProcesses()
{
// Sets G3 default values of control process flags.
// ---

  if (fLock) LockException();
  if (!fFlagVector) {
    // create vector of control process flag values
    fFlagVector = new TG4FlagVector;
  }  
  fFlagVector->SetG3Defaults();
}  

G4bool TG4PhysicsManager::IsSpecialCuts() const
{
// Returns true if any special cut value is set.
// ---

  for (G4int i=0; i<kNofParticlesWSP; i++)
  {  if ((*fIsCutVector)[i]) return true; }

  return false;
}

G4bool TG4PhysicsManager::IsSpecialFlags() const
{
// Returns true if any special flag value is set.
// ---

  for (G4int i=0; i<kNofParticlesWSP; i++)
  {  if ((*fIsFlagVector)[i]) return true; }

  return false;
}

TG3ParticleWSP TG4PhysicsManager::GetG3ParticleWSP(
                                  G4ParticleDefinition* particle) const 
{
// Returns TG3ParticleWSP constant for the specified particle.
// (See TG3ParticleWSP.h, too.)
// ---

  G4String name = particle->GetParticleName();     
  G4String pType = particle->GetParticleType();
    
  if (name == "gamma") {
    return kGamma;
  }  
  else if (name == "e-") {    
    return kElectron;
  }  
  else if (name == "e+") {   
    return kEplus;
  }  
  else if (( pType == "baryon" || pType == "meson" || pType == "nucleus" )) {
    if (particle->GetPDGCharge() == 0) { 
      return kNeutralHadron;
    }
    else  
      return kChargedHadron;
  }    
  else if ( name == "mu-" || name == "mu+" ) {
    return kMuon;
  }  
  else {
    return kNofParticlesWSP;
  }    
}  

void TG4PhysicsManager::GetG3ParticleWSPName(G4int particleWSP,
                                             G4String& name) const 
{
// Fills the passed name with the name/type of particle specified
// by TG3ParticleWSP constant.
// (See TG3ParticleWSP.h, too.)
// ---

  switch (particleWSP) {
    case kGamma:
      name = "Gamma";
      break;
    case kElectron:
      name = "Electron";
      break;
    case kEplus:
      name = "Eplus";
      break;
    case kNeutralHadron:
      name = "NeutralHadron";
      break;
    case kChargedHadron:
      name = "ChargedHadron";
      break;
    case kMuon:
      name = "Muon";
      break;
    case kAny:
      name = "Any";
      break;
    case kNofParticlesWSP:
      name = "NoSP";
      break;
    default:
      G4String text = "TG4PhysicsList::GetG3ParticleWSPName:\n";
      text = text + "   Wrong particleWSP."; 
      TG4Globals::Exception(text);
      break;      
  }
}
Commit	Line	Data
2817d3e2	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	#include <Randomize.hh>
	14
	15	#include <TDatabasePDG.h>
	16
	17	TG4PhysicsManager* TG4PhysicsManager::fgInstance = 0;
	18
	19	TG4PhysicsManager::TG4PhysicsManager()
	20	: fLock(false),
	21	fPhysicsList(0),
	22	fCutVector(0),
	23	fFlagVector(0)
	24	{
	25	//
	26	if (fgInstance) {
	27	TG4Globals::Exception(
	28	"TG4PhysicsManager: attempt to create two instances of singleton.");
	29	}
	30
	31	fgInstance = this;
	32
	33	// initialize fIsCutVector
	34	fIsCutVector = new TG4boolVector;
	35	G4int i;
	36	//for (i=0; i<kNofParticlesWSP; i++) fIsCutVector->insert(false);
	37	for (i=0; i<kNofParticlesWSP; i++) fIsCutVector->push_back(false);
	38
	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);
	43
	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");
	55
	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");
69	}
70
71	TG4PhysicsManager::TG4PhysicsManager(const TG4PhysicsManager& right) {
72	//
73	TG4Globals::Exception(
74	"Attempt to copy TG4PhysicsManager singleton.");
75	}
76
77	TG4PhysicsManager::~TG4PhysicsManager() {
78	//
79	delete fIsCutVector;
80	delete fIsFlagVector;
81	}
82
83	// operators
84
85	TG4PhysicsManager&
86	TG4PhysicsManager::operator=(const TG4PhysicsManager& right)
87	{
88	// check assignement to self
89	if (this == &right) return *this;
90
91	TG4Globals::Exception(
92	"Attempt to assign TG4PhysicsManager singleton.");
93
94	return *this;
95	}
96
97	// private methods
98
99	void TG4PhysicsManager::LockException() const
100	{
101	// Gives exception in case of attempt to modified physics
102	// setup after physics manager was locked.
103	// ---
104
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);
109	}
110
111	G4int TG4PhysicsManager::GetPDGEncoding(G4ParticleDefinition* particle)
112	{
113	// Returns the PDG code of particle;
114	// if standard PDG code is not defined the TDatabasePDG
115	// is used.
116	// ---
117
118	// get PDG encoding from G4 particle definition
119	G4int pdgEncoding = particle->GetPDGEncoding();
120
121	if (pdgEncoding == 0) {
122	// get PDG encoding from TDatabasePDG
123
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);
132	}
133
134	// get particle from TDatabasePDG
135	TDatabasePDG* pdgDB = TDatabasePDG::Instance();
136	TParticlePDG* particle = pdgDB->GetParticle(tname);
137	if (!particle) {
138	G4String text = "TG4PhysicsManager::GetPDGEncoding: \n";
139	text = text + " Particle " + tname;
140	text = text + " was not found in TDatabasePDG.";
141	TG4Globals::Exception(text);
142	}
143
144	// get PDG encoding
145	pdgEncoding = particle->PdgCode();
146	}
147
148	return pdgEncoding;
149	}
150
151	G4int TG4PhysicsManager::GetPDGEncoding(G4String particleName)
152	{
153	// Returns the PDG code of particle sepcified by name.
154	// ---
155
156	G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
157
158	G4ParticleDefinition* particle = 0;
159	particle = particleTable->FindParticle(particleName);
160	if (!particle) {
161	G4String text = "TG4PhysicsManager::GetPDGEncoding:\n";
162	text = text + " G4ParticleTable::FindParticle() " + particleName;
163	text = text + " failed.";
164	TG4Globals::Exception(text);
165	}
166
167	return GetPDGEncoding(particle);
168	}
169
170	void TG4PhysicsManager::SetCut(TG3Cut cut, G4double cutValue)
171	{
172	// Sets kinetic energy cut (in a G3-like way).
173	// ---
174
175	if (!fCutVector) {
176	// create vector of kinetic energy cut values
177	fCutVector = new TG4CutVector();
178	}
179	fCutVector->SetG3Cut(cut, cutValue);
180	SwitchIsCutVector(cut);
181	}
182
183	void TG4PhysicsManager::SetProcess(TG3Flag flag, G4int flagValue)
184	{
185	// Sets control process flag (in a G3-like way).
186	// ---
187
188	if (!fFlagVector) {
189	// create vector of control process flag values
190	fFlagVector = new TG4FlagVector;
191	}
192	fFlagVector->SetG3Flag(flag, flagValue);
193	}
194
195
196	Float_t TG4PhysicsManager::Xsec(char* ch, Float_t p1, Int_t i1, Int_t i2)
197	{
198	// Not yet implemented -> gives exception.
199	// ---
200
201	TG4Globals::Exception(
202	"TG4PhysicsManager::Xsec: not yet implemented.");
203
204	return 0.;
205	}
206
207	Int_t TG4PhysicsManager::IdFromPDG(Int_t pdgID) const
208	{
209	// G4 does not use the integer particle identifiers
210	// Id <-> PDG is identity.
211	// ---
212
213	return pdgID;
214	}
215
216	Int_t TG4PhysicsManager::PDGFromId(Int_t mcID) const
217	{
218	// G4 does not use integer particle identifiers
219	// Id <-> PDG is identity.
220	// ---
221
222	return mcID;
223	}
224
225	void TG4PhysicsManager::DefineParticles()
226	{
227	// ======
228	// Taken from TGeant3
229	//
230	// Use ENDF-6 mapping for ions = 10000z+10a+iso
231	// and add 1 000 000
232	// and numbers above 5 000 000 for special applications
233	//
234
235	const Int_t kion=10000000;
236	const Int_t kspe=50000000;
237
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 = 360024365.25;
243
244	TDatabasePDG *pdgDB = TDatabasePDG::Instance();
245
246	pdgDB->AddParticle("Deuteron","Deuteron",2*kGeV+8.071e-3,kTRUE,
247	0,1,"Ion",kion+10020);
248
249	pdgDB->AddParticle("Triton","Triton",3*kGeV+14.931e-3,kFALSE,
250	kHshGeV/(12.33*kYearsToSec),1,"Ion",kion+10030);
251
252	pdgDB->AddParticle("Alpha","Alpha",4*kGeV+2.424e-3,kTRUE,
253	kHshGeV/(12.33*kYearsToSec),2,"Ion",kion+20040);
254
255	pdgDB->AddParticle("HE3","HE3",3*kGeV+14.931e-3,kFALSE,
256	0,2,"Ion",kion+20030);
257
258	pdgDB->AddParticle("Cherenkov","Cherenkov",0,kFALSE,
259	0,0,"Special",kspe+50);
260
261	pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,
262	0,0,"Special",kspe+51);
263
264
265	// To do: define the PDG database extension
266	// in a common part.
267	//
268	// AliMC::ExtendPDGDatabase();
269	//
270	// end of "common" implementation
271	// ======
272
273	// map G4 particle names to TDatabasePDG names
274	// (the map is built only for particles that have not
275	// defined standard PDG encoding)
276
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");
284
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"));
293
294	// add verbose
295	G4cout << "Particle maps have been filled." << endl;
296	//fParticleNameMap.PrintAll();
297	//fParticlePDGMap.PrintAll();
298	}
299
300	void TG4PhysicsManager::SwitchIsCutVector(TG3Cut cut)
301	{
302	// Updates the vector of booleans (fIsCutVector) for the specified cut.
303	// ---
304
305	switch (cut) {
306	case kCUTGAM:
307	(*fIsCutVector)[kGamma] = true;
308	break;
309	case kBCUTE:
310	(*fIsCutVector)[kGamma] = true;
311	break;
312	case kBCUTM:
313	(*fIsCutVector)[kGamma] = true;
314	break;
315	case kCUTELE:
316	(*fIsCutVector)[kElectron] = true;
317	break;
318	case kDCUTE:
319	(*fIsCutVector)[kElectron] = true;
320	break;
321	case kDCUTM:
322	(*fIsCutVector)[kElectron] = true;
323	break;
324	case kCUTNEU:
325	(*fIsCutVector)[kNeutralHadron] = true;
326	break;
327	case kCUTHAD:
328	(*fIsCutVector)[kChargedHadron] = true;
329	break;
330	case kCUTMUO:
331	(*fIsCutVector)[kMuon] = true;
332	break;
333	default:
334	break;
335	}
336	}
337
338	void TG4PhysicsManager::SwitchIsFlagVector(TG3Flag flag)
339	{
340	// Updates the vector of booleans (fIsFlagVector) for the specified flag.
341	// ---
342
343	switch (flag) {
344	case kPAIR:
345	// gamma
346	(*fIsFlagVector)[kGamma] = true;
347	break;
348	case kCOMP:
349	// gamma
350	(*fIsFlagVector)[kGamma] = true;
351	break;
352	case kPHOT:
353	// gamma
354	(*fIsFlagVector)[kGamma] = true;
355	break;
356	case kPFIS:
357	// gamma
358	(*fIsFlagVector)[kGamma] = true;
359	break;
360	case kDRAY:
361	// all charged particles
362	(*fIsFlagVector)[kElectron] = true;
363	(*fIsFlagVector)[kEplus] = true;
364	(*fIsFlagVector)[kChargedHadron] = true;
365	(*fIsFlagVector)[kMuon] = true;
366	break;
367	case kANNI:
368	// e+ only
369	(*fIsFlagVector)[kEplus] = true;
370	break;
371	case kBREM:
372	// e-/e+, muons
373	(*fIsFlagVector)[kElectron] = true;
374	(*fIsFlagVector)[kEplus] = true;
375	(*fIsFlagVector)[kMuon] = true;
376	break;
377	case kHADR:
378	// hadrons
379	(*fIsFlagVector)[kNeutralHadron] = true;
380	(*fIsFlagVector)[kChargedHadron] = true;
381	break;
382	case kMUNU:
383	// muons
384	(*fIsFlagVector)[kMuon] = true;
385	break;
386	case kDCAY:
387	// any
388	(*fIsFlagVector)[kAny] = true;
389	break;
390	case kLOSS:
391	// all charged particles
392	(*fIsFlagVector)[kElectron] = true;
393	(*fIsFlagVector)[kEplus] = true;
394	(*fIsFlagVector)[kChargedHadron] = true;
395	(*fIsFlagVector)[kMuon] = true;
396	break;
397	case kMULS:
398	// all charged particles
399	(*fIsFlagVector)[kElectron] = true;
400	(*fIsFlagVector)[kEplus] = true;
401	(*fIsFlagVector)[kChargedHadron] = true;
402	(*fIsFlagVector)[kMuon] = true;
403	break;
404	default:
405	break;
406	}
407	}
408
409	TG3Cut TG4PhysicsManager::GetG3Cut(G4String cutName)
410	{
411	// Retrieves corresponding TG3Cut constant from the cutName.
412	// ---
413
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;
424	}
425
426	TG3Flag TG4PhysicsManager::GetG3Flag(G4String flagName)
427	{
428	// Retrieves corresponding TG3Flag constant from the flagName.
429	// ---
430
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;
444	}
445
446	// public methods
447
448	void TG4PhysicsManager::BuildPhysics()
449	{
450	// Empty function - not needed in G4.
451	// (Physics is built within /run/initialize.)
452
453	TG4Globals::Warning(
454	"TG4PhysicsManager::BuildPhysics: is empty function in G4 MC.");
455	}
456
457	void TG4PhysicsManager::SetCut(const char* cutName, Float_t cutValue)
458	{
459	// Sets the specified cut.
460	// ---
461
462	if (fLock) LockException();
463	TG3Cut g3Cut = GetG3Cut(cutName);
464	if (g3Cut != kNoG3Cuts)
465	SetCut(g3Cut, cutValue);
466	else {
467	G4String text = "TG4PhysicsManager::SetCut:\n";
468	text = text + " Parameter " + cutName;
469	text = text + " is not implemented.";
470	TG4Globals::Warning(text);
471	}
472	}
473
474	void TG4PhysicsManager::SetProcess(const char* flagName, Int_t flagValue)
475	{
476	// Sets the specified process control.
477	// ---
478
479	if (fLock) LockException();
480	TG3Flag g3Flag = GetG3Flag(flagName);
481	if (g3Flag != kNoG3Flags)
482	SetProcess(g3Flag, flagValue);
483	else {
484	G4String text = "TG4PhysicsManager::SetProcess:\n";
485	text = text + " Parameter " + flagName;
486	text = text + " is not implemented.";
487	TG4Globals::Warning(text);
488	}
489	}
490
491	void TG4PhysicsManager::SetProcessActivation()
492	{
493	// (In)Activates built processes according
494	// to the setup in fFlagVector.
495	// ---
496
497	if (fPhysicsList) {
498	// temporarily excluded
499	// fPhysicsList->SetProcessActivation();
500	}
501	else {
502	G4String text = "TG4PhysicsManager::SetProcessActivation:\n";
503	text = text + " There is no physics list set.";
504	TG4Globals::Exception(text);
505	}
506	}
507
508	G4int TG4PhysicsManager::GetPDGEncodingFast(G4ParticleDefinition* particle)
509	{
510	// Returns the PDG code of particle;
511	// if standard PDG code is not defined the preregistred
512	// fParticlePDGMap is used.
513	// ---
514
515	// get PDG encoding from G4 particle definition
516	G4int pdgEncoding = particle->GetPDGEncoding();
517
518	if (pdgEncoding == 0) {
519	// use FParticlePDGMap if standard PDG code is not defined
520	G4String name = particle->GetParticleName();
521	pdgEncoding = fParticlePDGMap.GetSecond(name);
522	}
523
524	return pdgEncoding;
525	}
526
527	G4bool TG4PhysicsManager::CheckCutWithCutVector(G4String name,
528	G4double value, TG3Cut& cut)
529	{
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.
535	// ---
536
537	// convert cut name -> TG3Cut
538	cut = GetG3Cut(name);
539
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);
547	return true;
548	}
549	else return false;
550	}
551	return false;
552	}
553
554	G4bool TG4PhysicsManager::CheckFlagWithFlagVector(G4String name,
555	G4double value, TG3Flag& flag)
556	{
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.
562	// ---
563
564	// convert flag name -> TG3Flag
565	flag = GetG3Flag(name);
566
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);
572	return true;
573	}
574	else return false;
575	}
576	return false;
577	}
578
579	G4bool TG4PhysicsManager::CheckCutWithG3Defaults(G4String name,
580	G4double value, TG3Cut& cut)
581	{
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.
587	// ---
588
589	// convert cut name -> TG3Cut
590	cut = GetG3Cut(name);
591
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);
597	return true;
598	}
599	else return false;
600	}
601	return false;
602	}
603
604	G4bool TG4PhysicsManager::CheckFlagWithG3Defaults(G4String name,
605	G4double value, TG3Flag& flag)
606	{
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.
612	// ---
613
614	// convert flag name -> TG3Flag
615	flag = GetG3Flag(name);
616
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);
622	return true;
623	}
624	else return false;
625	}
626	return false;
627	}
628
629	void TG4PhysicsManager::SetG3DefaultCuts()
630	{
631	// Sets G3 default values of kinetic energy cuts.
632	// ---
633
634	if (fLock) LockException();
635	if (!fCutVector) {
636	// create vector of kinetic energy cut values
637	fCutVector = new TG4CutVector();
638	}
639	fCutVector->SetG3Defaults();
640	}
641
642	void TG4PhysicsManager::SetG3DefaultProcesses()
643	{
644	// Sets G3 default values of control process flags.
645	// ---
646
647	if (fLock) LockException();
648	if (!fFlagVector) {
649	// create vector of control process flag values
650	fFlagVector = new TG4FlagVector;
651	}
652	fFlagVector->SetG3Defaults();
653	}
654
655	G4bool TG4PhysicsManager::IsSpecialCuts() const
656	{
657	// Returns true if any special cut value is set.
658	// ---
659
660	for (G4int i=0; i<kNofParticlesWSP; i++)
661	{ if ((*fIsCutVector)[i]) return true; }
662
663	return false;
664	}
665
666	G4bool TG4PhysicsManager::IsSpecialFlags() const
667	{
668	// Returns true if any special flag value is set.
669	// ---
670
671	for (G4int i=0; i<kNofParticlesWSP; i++)
672	{ if ((*fIsFlagVector)[i]) return true; }
673
674	return false;
675	}
676
677	TG3ParticleWSP TG4PhysicsManager::GetG3ParticleWSP(
678	G4ParticleDefinition* particle) const
679	{
680	// Returns TG3ParticleWSP constant for the specified particle.
681	// (See TG3ParticleWSP.h, too.)
682	// ---
683
684	G4String name = particle->GetParticleName();
685	G4String pType = particle->GetParticleType();
686
687	if (name == "gamma") {
688	return kGamma;
689	}
690	else if (name == "e-") {
691	return kElectron;
692	}
693	else if (name == "e+") {
694	return kEplus;
695	}
696	else if (( pType == "baryon" \|\| pType == "meson" \|\| pType == "nucleus" )) {
697	if (particle->GetPDGCharge() == 0) {
698	return kNeutralHadron;
699	}
700	else
701	return kChargedHadron;
702	}
703	else if ( name == "mu-" \|\| name == "mu+" ) {
704	return kMuon;
705	}
706	else {
707	return kNofParticlesWSP;
708	}
709	}
710
711	void TG4PhysicsManager::GetG3ParticleWSPName(G4int particleWSP,
712	G4String& name) const
713	{
714	// Fills the passed name with the name/type of particle specified
715	// by TG3ParticleWSP constant.
716	// (See TG3ParticleWSP.h, too.)
717	// ---
718
719	switch (particleWSP) {
720	case kGamma:
721	name = "Gamma";
722	break;
723	case kElectron:
724	name = "Electron";
725	break;
726	case kEplus:
727	name = "Eplus";
728	break;
729	case kNeutralHadron:
730	name = "NeutralHadron";
731	break;
732	case kChargedHadron:
733	name = "ChargedHadron";
734	break;
735	case kMuon:
736	name = "Muon";
737	break;
738	case kAny:
739	name = "Any";
740	break;
741	case kNofParticlesWSP:
742	name = "NoSP";
743	break;
744	default:
745	G4String text = "TG4PhysicsList::GetG3ParticleWSPName:\n";
746	text = text + " Wrong particleWSP.";
747	TG4Globals::Exception(text);
748	break;
749	}
750	}
751