diff --git a/NPLib/Physics/CMakeLists.txt b/NPLib/Physics/CMakeLists.txt
index ed76e00bfe5f56332048acb0e91d6fb441c178dc..60c974730c14561ed36ad99e3a31031f0d34e78d 100644
--- a/NPLib/Physics/CMakeLists.txt
+++ b/NPLib/Physics/CMakeLists.txt
@@ -26,5 +26,4 @@ target_link_libraries(NPInteractionCoordinates ${ROOT_LIBRARIES} )
add_library(NPReactionConditions SHARED TReactionConditions.cxx TReactionConditionsDict.cxx)
target_link_libraries(NPReactionConditions ${ROOT_LIBRARIES} )
-
install(FILES NPDecay.h NPBeam.h NPEnergyLoss.h NPFunction.h NPNucleus.h NPReaction.h NPQFS.h TInitialConditions.h TInteractionCoordinates.h TReactionConditions.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/Physics/NPBeam.cxx b/NPLib/Physics/NPBeam.cxx
index e7644c7e4076621954b17203ce309dcb608e451a..f7c7551413e47cd51ee2955f29fe38101eca41a9 100644
--- a/NPLib/Physics/NPBeam.cxx
+++ b/NPLib/Physics/NPBeam.cxx
@@ -81,9 +81,10 @@ Beam::Beam(){
fXThetaXHist = new TH2F(Form("XThetaXHis_%i",offset),"XThetaXHis",1,0,1,1,0,1);
fYPhiYHist = new TH2F(Form("YPhiYHist_%i",offset),"YPhiYHist",1,0,1,1,0,1);
}
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-Beam::Beam(string isotope){
- SetUp(isotope);
+Beam::Beam(string name){
+ SetUp(name);
fEnergy = 0;
fExcitationEnergy = 0;
fSigmaEnergy = -1 ;
@@ -116,7 +117,6 @@ Beam::Beam(string isotope){
fYPhiYHist = new TH2F(Form("YPhiYHist_%i",offset),"YPhiYHist",1,0,1,1,0,1);
}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Beam::~Beam(){
}
@@ -126,6 +126,7 @@ void Beam::ReadConfigurationFile(string Path){
NPL::InputParser parser(Path);
ReadConfigurationFile(parser);
}
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Beam::ReadConfigurationFile(NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Beam");
@@ -253,20 +254,20 @@ void Beam::Print() const {
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Beam::SetTargetSize(double TargetSize){
fTargetSize = TargetSize;
fEffectiveTargetSize = fTargetSize*cos(fTargetAngle);
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Beam::SetTargetThickness(double TargetThickness){
fTargetThickness = TargetThickness;
fEffectiveTargetThickness = fTargetThickness/cos(fTargetAngle);
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Beam::SetTargetAngle(double TargetAngle){
fTargetAngle = TargetAngle;
fEffectiveTargetSize = fTargetSize*cos(fTargetAngle);
diff --git a/NPLib/Physics/NPBeam.h b/NPLib/Physics/NPBeam.h
index a82e28e80179b722d49eb0aa5a44c14466c6c49b..c20c3f519a6a18c5e283551c9ffc653a9c9c1761 100644
--- a/NPLib/Physics/NPBeam.h
+++ b/NPLib/Physics/NPBeam.h
@@ -111,7 +111,7 @@ namespace NPL{
TH1F* GetEnergyHist () const {return fEnergyHist;}
TH2F* GetXThetaXHist () const {return fXThetaXHist;}
TH2F* GetYPhiYHist () const {return fYPhiYHist;}
- int GetVerboseLevel() const {return fVerboseLevel;}
+ int GetVerboseLevel() const {return fVerboseLevel;}
private: // Event Generation private variable
double fTargetSize;
diff --git a/NPLib/Physics/NPNucleus.cxx b/NPLib/Physics/NPNucleus.cxx
index cf7029c70c3f8840f21435c3e3fbeb1c973eeb19..59c02dc918177e6277a20e18953b149782e9f2cc 100644
--- a/NPLib/Physics/NPNucleus.cxx
+++ b/NPLib/Physics/NPNucleus.cxx
@@ -52,6 +52,7 @@ Nucleus::Nucleus(){
fCharge= 0;
fAtomicWeight= 0;
fMassExcess= 0;
+ fMass=0;
fExcitationEnergy=0;
fSpinParity= "";
fSpin= 0;
@@ -60,14 +61,14 @@ Nucleus::Nucleus(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-Nucleus::Nucleus(string isotope){
- SetUp(isotope);
+Nucleus::Nucleus(string name){
+ SetUp(name);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-Nucleus::Nucleus(string isotope, const string& path){
- SetUp(isotope);
- LoadENSDF(isotope, path);
+Nucleus::Nucleus(string name, const string& path){
+ SetUp(name);
+ LoadENSDF(name, path);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -92,6 +93,7 @@ Nucleus::Nucleus(string name, vector<string> subpart, double binding,double Ex,
fSpin= Spin;
fParity= Parity;
fLifeTime = LifeTime;
+ fMass=Mass;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -177,21 +179,50 @@ void Nucleus::LoadENSDF(const string& isotope, const string& pathENSDF)
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void Nucleus::SetUp(string isotope){
+void Nucleus::SetUp(string name){
+ if(name=="electron"){
+ fName = "electron";
+ fCharge = -1;
+ fAtomicWeight = 0;
+ fMassExcess = 0;
+ fMass = electron_mass_c2;
+ fExcitationEnergy =0;
+ fSpinParity = "0.5";
+ fSpin = 0.5;
+ fParity = "";
+ fLifeTime = -1;
+ return;
+ }
+
+ else if(name=="gamma"){
+ fName = "gamma";
+ fCharge = 0;
+ fAtomicWeight = 0;
+ fMassExcess = 0;
+ fMass = 0;
+ fExcitationEnergy =0;
+ fSpinParity = "0";
+ fSpin = 0;
+ fParity = "";
+ fLifeTime = -1;
+ return;
+ }
+
+
//----------- Constructor Using nubtab12.asc by name----------
// open the file to read and check if it is open
fExcitationEnergy=0;
// Replace the n,p,d,t,a by there standard name:
- if(isotope=="p") isotope="1H";
- else if(isotope=="n") isotope="1n";
- else if(isotope=="d") isotope="2H";
- else if(isotope=="t") isotope="3H";
- else if(isotope=="a") isotope="4He";
- else if(isotope=="n") isotope="1n";
- else if(isotope=="neutron") isotope="1n";
- else if(isotope=="g") isotope="gamma";
- else if(isotope=="gamma") isotope="gamma";
- else if(isotope=="4n"){ // tetraneutron @Eres = 0
+ if(name=="p") name="1H";
+ else if(name=="n") name="1n";
+ else if(name=="d") name="2H";
+ else if(name=="t") name="3H";
+ else if(name=="a") name="4He";
+ else if(name=="n") name="1n";
+ else if(name=="neutron") name="1n";
+ else if(name=="g") name="gamma";
+ else if(name=="gamma") name="gamma";
+ else if(name=="4n"){ // tetraneutron @Eres = 0
string line = "004 0000 4n 32285.268 0.0005 219.4 ys 0.6 1/2+ 00 02PaDGt B-=100";
Extract(line.data());
return;
@@ -213,7 +244,7 @@ void Nucleus::SetUp(string isotope){
space = s_name.find_first_of(" ");
s_name.resize(space);
- if (s_name.find(isotope) != string::npos && s_name.length() == isotope.length()) break;
+ if (s_name.find(name) != string::npos && s_name.length() == name.length()) break;
}
Extract(line.data());
}
@@ -410,7 +441,7 @@ void Nucleus::Extract(string line){
if (s_spinparity.find("19/2") != string::npos) fSpin = 9.5 ;
if (s_spinparity.find("21/2") != string::npos) fSpin = 10.5 ;
GetNucleusName();
-
+ fMass=Mass();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -636,26 +667,3 @@ void Nucleus::PrintThreshold() const {
cout << " Sa : " << GetSa() << " MeV" << endl;
}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/NPLib/Physics/NPNucleus.h b/NPLib/Physics/NPNucleus.h
index 789a75d571766bac09c0afdec7e3ac7c507ac5a3..a44acc39f6a1a1cf499b2e927723440d0e49e512 100644
--- a/NPLib/Physics/NPNucleus.h
+++ b/NPLib/Physics/NPNucleus.h
@@ -64,6 +64,7 @@ namespace NPL {
int fCharge; // Nucleus charge
int fAtomicWeight; // Nucleus atomic weight
double fMassExcess; // Nucleus mass excess in keV
+ double fMass;
string fSpinParity; // Nucleus spin and parity
double fSpin; // Nucleus spin
string fParity; // Nucleus parity
@@ -85,14 +86,14 @@ namespace NPL {
public:
void EnergyToBrho(double Q=-1000);
void EnergyToTof();
- void BetaToVelocity();
+ void BetaToVelocity();
void BrhoToEnergy(double Q=-1000);
void BrhoToTof() {BrhoToEnergy(); EnergyToTof();}
- void TofToEnergy();
- void TofToBrho() {TofToEnergy(); EnergyToBrho();}
- void EnergyToBeta();
- void BetaToEnergy();
- void BetaToGamma();
+ void TofToEnergy();
+ void TofToBrho() {TofToEnergy(); EnergyToBrho();}
+ void EnergyToBeta();
+ void BetaToEnergy();
+ void BetaToGamma();
double DopplerCorrection(double EnergyLabGamma, double ThetaLabGamma);
@@ -101,23 +102,23 @@ namespace NPL {
public :
void GetNucleusName();
- string GetName() const {return fNucleusName;}
+ string GetName() const {return fName;}
int GetZ() const {return fCharge;}
int GetA() const {return fAtomicWeight;}
double GetMassExcess() const {return fMassExcess;}
string GetSpinParity() const {return fSpinParity;}
double GetSpin() const {return fSpin;}
string GetParity() const {return fParity;}
- double GetLifeTime() const {return fLifeTime;}
- double GetLifeTimeError() const {return fLifeTimeErr;}
+ double GetLifeTime() const {return fLifeTime;}
+ double GetLifeTimeError() const {return fLifeTimeErr;}
double GetEnergy() const {return fKineticEnergy;}
double GetBrho() const {return fBrho;}
double GetTimeOfFlight() const {return fTimeOfFlight;}
double GetBeta() const {return fBeta;}
double GetGamma() const {return fGamma;}
double GetVelocity() const {return fVelocity;}
- TLorentzVector GetEnergyImpulsion() const {return fEnergyImpulsion;}
- double GetExcitationEnergy() const {return fExcitationEnergy;}
+ TLorentzVector GetEnergyImpulsion() const {return fEnergyImpulsion;}
+ double GetExcitationEnergy() const {return fExcitationEnergy;}
void SetName(const char* name) {fName = name;}
void SetZ(int charge) {fCharge = charge;}
void SetA(int mass) {fAtomicWeight = mass;}
diff --git a/NPLib/Physics/nubtab16.asc b/NPLib/Physics/nubtab16.asc
index 414448a3b97010b9440ef5aa89987ed66d681c0c..63da4518a2a376db085b36bbe80915e9b04401f6 100644
--- a/NPLib/Physics/nubtab16.asc
+++ b/NPLib/Physics/nubtab16.asc
@@ -1,4 +1,3 @@
-000 0000 gamma 0.0 0.0000 stbl 0+ 06 1932 IS=0.000000 00
001 0000 1n 8071.3171 0.0005 613.9 s 0.6 1/2+ 06 1932 B-=100
001 0010 1H 7288.9706 0.0001 stbl 1/2+ 06 11Be53d 1920 IS=99.9885 70
002 0010 2H 13135.7217 0.0001 stbl 1+ 03 1932 IS=0.0115 70
diff --git a/NPSimulation/EventGenerator/EventGeneratorBeam.cc b/NPSimulation/EventGenerator/EventGeneratorBeam.cc
index 748fa663713a2de716f2c4edf0685262abf71866..d169d083067a2a4f8f3c1217c4e8a81957db65c4 100644
--- a/NPSimulation/EventGenerator/EventGeneratorBeam.cc
+++ b/NPSimulation/EventGenerator/EventGeneratorBeam.cc
@@ -76,8 +76,12 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent){
// Define the particle to be shoot
if(m_Beam->GetZ()==0 && m_Beam->GetA()==1)
m_particle = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
-
- else
+ else if(m_Beam->GetName()=="electron")
+ m_particle = G4ParticleTable::GetParticleTable()->FindParticle("e-");
+ else if(m_Beam->GetName()=="gamma")
+ m_particle = G4ParticleTable::GetParticleTable()->FindParticle("gamma");
+
+ else
m_particle =
G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(m_Beam->GetZ(), m_Beam->GetA() ,m_Beam->GetExcitationEnergy());
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index f10c3859c91d038df1d7ff0e8757921553aac9c2..9d9d3795ba1f64c0d2a14cd5eab9392e4ceb49ad 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -101,6 +101,7 @@ G4bool NPS::BeamReaction::IsApplicable(const G4ParticleDefinition& particleType)
static std::string particleName;
particleName = particleType.GetParticleName();
if(particleName=="neutron") particleName="n1";
+ else if(particleName=="e-") particleName="electron";
if (particleName.find(m_BeamName) != std::string::npos) {
return true;
}
@@ -127,12 +128,12 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
bool is_first = (to_entrance==0);
if(is_first && m_shoot){
- /* std::cout << "Something went wrong in beam reaction, m_shoot cannot be true at beginning of event" << std::endl;
- std::cout << "rand: " << m_rand << std::endl;
- std::cout << "length: " << m_length << std::endl;
- std::cout << "step: " << m_StepSize << std::endl;
- std::cout << "Z: " << m_Z << std::endl;
- std::cout << "S: " << m_S << std::endl;*/
+ /* std::cout << "Something went wrong in beam reaction, m_shoot cannot be true at beginning of event" << std::endl;
+ std::cout << "rand: " << m_rand << std::endl;
+ std::cout << "length: " << m_length << std::endl;
+ std::cout << "step: " << m_StepSize << std::endl;
+ std::cout << "Z: " << m_Z << std::endl;
+ std::cout << "S: " << m_S << std::endl;*/
m_shoot = false;
}
@@ -144,7 +145,7 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
m_ReactionConditions->Clear();
m_shoot=true;
}
-
+
// curviligne coordinate along beam path
m_S = to_entrance - 0.5*(to_exit+to_entrance);
m_length = m_Z-m_S;
@@ -194,18 +195,18 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
->GetVolume()
->GetLogicalVolume()
->GetMaterial();
-
+
double energy = PrimaryTrack->GetKineticEnergy();
double speed = PrimaryTrack->GetVelocity();
double time = PrimaryTrack->GetGlobalTime()+m_length/speed;
-
+
double reac_energy = SlowDownBeam (
- PrimaryTrack->GetParticleDefinition(),
- energy,
- m_length,
- material);
-
+ PrimaryTrack->GetParticleDefinition(),
+ energy,
+ m_length,
+ material);
+
G4ThreeVector ldir = pdirection;
ldir *= m_length;
@@ -225,27 +226,32 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
///////////////////////////////
// Two-Body Reaction Case /////
///////////////////////////////
-
- //////Define the kind of particle to shoot////////
- // Nucleus 3
- int LightZ = m_Reaction.GetNucleus3()->GetZ();
- int LightA = m_Reaction.GetNucleus3()->GetA();
static G4IonTable* IonTable
= G4ParticleTable::GetParticleTable()->GetIonTable();
+ //////Define the kind of particle to shoot////////
+ // Particle 3
G4ParticleDefinition* LightName;
+ if(m_Reaction.GetNucleus3()->GetName()=="electron"){
+ LightName=G4Electron::Definition();
+ }
+ else{
+ int LightZ = m_Reaction.GetNucleus3()->GetZ();
+ int LightA = m_Reaction.GetNucleus3()->GetA();
+ if (LightZ == 0 && LightA == 1){
+ LightName = G4Neutron::Definition();
+ }
+ else {
+ if (m_Reaction.GetUseExInGeant4())
+ LightName
+ = IonTable->GetIon(LightZ, LightA, m_Reaction.GetExcitation3() * MeV);
+ else
+ LightName = IonTable->GetIon(LightZ, LightA);
+ }
- if (LightZ == 0 && LightA == 1) // neutron is special case
- {
- LightName = G4Neutron::Definition();
- } else {
- if (m_Reaction.GetUseExInGeant4())
- LightName
- = IonTable->GetIon(LightZ, LightA, m_Reaction.GetExcitation3() * MeV);
- else
- LightName = IonTable->GetIon(LightZ, LightA);
}
+
// Nucleus 4
G4int HeavyZ = m_Reaction.GetNucleus4()->GetZ();
G4int HeavyA = m_Reaction.GetNucleus4()->GetA();
@@ -490,7 +496,7 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
TLorentzVector* P_A = m_QFS.GetEnergyImpulsionLab_A();
TLorentzVector* P_B = m_QFS.GetEnergyImpulsionLab_B();
-
+
G4ThreeVector momentum_kineB_beam( P_B->Px(), P_B->Py(), P_B->Pz() );
momentum_kineB_beam = momentum_kineB_beam.unit();
TKEB = P_B->Energy() - m_QFS.GetNucleusB()->Mass();
diff --git a/Projects/Strasse/PhysicsListOption.txt b/Projects/Strasse/PhysicsListOption.txt
index ff38d11c7d6361c44a8d00159a7a5543a22c25cd..56772eef01e5db197b2641b663872b7d1506c70f 100644
--- a/Projects/Strasse/PhysicsListOption.txt
+++ b/Projects/Strasse/PhysicsListOption.txt
@@ -8,4 +8,4 @@ StoppingPhysics 0
OpticalPhysics 0
HadronPhysicsINCLXX 0
HadronPhysicsQGSP_BIC_HP 0
-Decay 1
+Decay 0