Commit 32378f47 authored by Adrien Matta's avatar Adrien Matta
Browse files

* Adjusting NPBeam and NPNucleus to deal with ee reaction

        - NPNucleus will be renamed subsequently to NPParticle
parent ecd97054
Pipeline #83602 passed with stages
in 15 minutes and 59 seconds
......@@ -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})
......@@ -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);
......
......@@ -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;
......
......@@ -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;
}
......@@ -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;}
......
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
......
......@@ -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());
......
......@@ -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();
......
......@@ -8,4 +8,4 @@ StoppingPhysics 0
OpticalPhysics 0
HadronPhysicsINCLXX 0
HadronPhysicsQGSP_BIC_HP 0
Decay 1
Decay 0
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment