Skip to content
Snippets Groups Projects
Commit 309fd4a8 authored by Morfouace's avatar Morfouace
Browse files

Adding GetEcm function to NPReaction / usefull for resonant scattering 

experiment
parent 7f9d7212
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
NPLib/Physics/NPReaction.cxx
+ 28
− 39
View file @ 309fd4a8
...@@ -80,7 +80,7 @@ ClassImp(Reaction) ...@@ -80,7 +80,7 @@ ClassImp(Reaction)
fCrossSectionHist = NULL; fCrossSectionHist = NULL;
fExcitationEnergyHist = NULL; fExcitationEnergyHist = NULL;
fDoubleDifferentialCrossSectionHist = NULL ; fDoubleDifferentialCrossSectionHist = NULL ;
fshoot3=true; fshoot3=true;
fshoot4=true; fshoot4=true;
...@@ -90,7 +90,7 @@ ClassImp(Reaction) ...@@ -90,7 +90,7 @@ ClassImp(Reaction)
// This constructor aim to provide a fast way to instantiate a reaction without input file // This constructor aim to provide a fast way to instantiate a reaction without input file
// The string should be of the form "A(b,c)D@E" with E the ernegy of the beam in MeV // The string should be of the form "A(b,c)D@E" with E the ernegy of the beam in MeV
Reaction::Reaction(string reaction){ Reaction::Reaction(string reaction){
// Instantiate the parameter to default // Instantiate the parameter to default
// Analyse the reaction and extract revelant information // Analyse the reaction and extract revelant information
string A,b,c,D,E; string A,b,c,D,E;
unsigned int i=0; unsigned int i=0;
...@@ -189,7 +189,7 @@ double Reaction::ShootRandomThetaCM(){ ...@@ -189,7 +189,7 @@ double Reaction::ShootRandomThetaCM(){
int binY; int binY;
// Those test are there for the tail event of the energy distribution // Those test are there for the tail event of the energy distribution
// In case the energy is outside the range of the 2D histo we take the // In case the energy is outside the range of the 2D histo we take the
// closest available CS // closest available CS
if(Y->FindBin(fBeamEnergy) > Y->GetLast()) if(Y->FindBin(fBeamEnergy) > Y->GetLast())
binY = Y->GetLast()-1; binY = Y->GetLast()-1;
...@@ -252,33 +252,27 @@ void Reaction::KineRelativistic(double& ThetaLab3, double& KineticEnergyLab3, ...@@ -252,33 +252,27 @@ void Reaction::KineRelativistic(double& ThetaLab3, double& KineticEnergyLab3,
} }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Reaction::ReconstructRelativistic(double EnergyLab, double ThetaLab, double PhiLab){ double Reaction::ReconstructRelativistic(double EnergyLab, double ThetaLab){
// EnergyLab in MeV // EnergyLab in MeV
// ThetaLab in rad // ThetaLab in rad
double E3 = m3 + EnergyLab; double E3 = m3 + EnergyLab;
double p_Lab_3 = sqrt(E3*E3 - m3*m3); double p_Lab_3 = sqrt(E3*E3 - m3*m3);
TVector3 p_Lab_3_vec = TVector3(0,0,1); fEnergyImpulsionLab_3 = TLorentzVector(p_Lab_3*sin(ThetaLab),0,p_Lab_3*cos(ThetaLab),E3);
p_Lab_3_vec.SetMagThetaPhi(p_Lab_3, ThetaLab, PhiLab);
fEnergyImpulsionLab_3 = TLorentzVector(p_Lab_3_vec, E3);
fEnergyImpulsionLab_4 = fTotalEnergyImpulsionLab - fEnergyImpulsionLab_3; fEnergyImpulsionLab_4 = fTotalEnergyImpulsionLab - fEnergyImpulsionLab_3;
double Eex = fEnergyImpulsionLab_4.Mag() - fNuclei4.Mass(); double Eex = fEnergyImpulsionLab_4.Mag() - fNuclei4.Mass();
//by Shuya 171005
//cout << m1 << " " << m2 << " " << fNuclei4.Mass() << " " << m3 << " " << EnergyLab << " " << ThetaLab << " " << ThetaLab/deg << " " << Eex << endl;
return Eex; return Eex;
} }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//Return ThetaCM //Return ThetaCM
double Reaction::EnergyLabToThetaCM(double EnergyLab, double ThetaLab, double PhiLab){ double Reaction::EnergyLabToThetaCM(double EnergyLab, double ThetaLab){
double E3 = m3 + EnergyLab; double E3 = m3 + EnergyLab;
double p_Lab_3 = sqrt(E3*E3 - m3*m3); double p_Lab_3 = sqrt(E3*E3 - m3*m3);
TVector3 p_Lab_3_vec = TVector3(0,0,1);
p_Lab_3_vec.SetMagThetaPhi(p_Lab_3, ThetaLab, PhiLab);
fEnergyImpulsionLab_3 = TLorentzVector(p_Lab_3_vec, E3); fEnergyImpulsionLab_3 = TLorentzVector(p_Lab_3*sin(ThetaLab),0,p_Lab_3*cos(ThetaLab),E3);
fEnergyImpulsionCM_3 = fEnergyImpulsionLab_3; fEnergyImpulsionCM_3 = fEnergyImpulsionLab_3;
fEnergyImpulsionCM_3.Boost(0,0,-BetaCM); fEnergyImpulsionCM_3.Boost(0,0,-BetaCM);
...@@ -322,7 +316,7 @@ void Reaction::ReadConfigurationFile(NPL::InputParser parser){ ...@@ -322,7 +316,7 @@ void Reaction::ReadConfigurationFile(NPL::InputParser parser){
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("TwoBodyReaction"); vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("TwoBodyReaction");
if(blocks.size()>0 && NPOptionManager::getInstance()->GetVerboseLevel()) if(blocks.size()>0 && NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "\033[1;35m//// Two body reaction found " << endl; cout << endl << "\033[1;35m//// Two body reaction found " << endl;
vector<string> token1 = {"Beam","Target","Light","Heavy"}; vector<string> token1 = {"Beam","Target","Light","Heavy"};
vector<string> token2 = {"Beam","Target","Nuclei3","Nuclei4"}; vector<string> token2 = {"Beam","Target","Nuclei3","Nuclei4"};
...@@ -349,21 +343,21 @@ void Reaction::ReadConfigurationFile(NPL::InputParser parser){ ...@@ -349,21 +343,21 @@ void Reaction::ReadConfigurationFile(NPL::InputParser parser){
fNuclei2 = Nucleus(blocks[i]->GetString("Target")); fNuclei2 = Nucleus(blocks[i]->GetString("Target"));
fNuclei3 = Nucleus(blocks[i]->GetString("Nuclei3")); fNuclei3 = Nucleus(blocks[i]->GetString("Nuclei3"));
fNuclei4 = Nucleus(blocks[i]->GetString("Nuclei4")); fNuclei4 = Nucleus(blocks[i]->GetString("Nuclei4"));
} }
else{ else{
cout << "ERROR: check your input file formatting \033[0m" << endl; cout << "ERROR: check your input file formatting \033[0m" << endl;
exit(1); exit(1);
} }
if(blocks[i]->HasToken("ExcitationEnergyLight")) if(blocks[i]->HasToken("ExcitationEnergyLight"))
fExcitation3 = blocks[i]->GetDouble("ExcitationEnergyLight","MeV"); fExcitation3 = blocks[i]->GetDouble("ExcitationEnergyLight","MeV");
else if(blocks[i]->HasToken("ExcitationEnergy3")) else if(blocks[i]->HasToken("ExcitationEnergy3"))
fExcitation3 = blocks[i]->GetDouble("ExcitationEnergy3","MeV"); fExcitation3 = blocks[i]->GetDouble("ExcitationEnergy3","MeV");
if(blocks[i]->HasToken("ExcitationEnergyHeavy")) if(blocks[i]->HasToken("ExcitationEnergyHeavy"))
fExcitation4 = blocks[i]->GetDouble("ExcitationEnergyHeavy","MeV"); fExcitation4 = blocks[i]->GetDouble("ExcitationEnergyHeavy","MeV");
else if(blocks[i]->HasToken("ExcitationEnergy4")) else if(blocks[i]->HasToken("ExcitationEnergy4"))
fExcitation4 = blocks[i]->GetDouble("ExcitationEnergy4","MeV"); fExcitation4 = blocks[i]->GetDouble("ExcitationEnergy4","MeV");
if(blocks[i]->HasToken("ExcitationEnergyDistribution")){ if(blocks[i]->HasToken("ExcitationEnergyDistribution")){
vector<string> file = blocks[i]->GetVectorString("ExcitationEnergyDistribution"); vector<string> file = blocks[i]->GetVectorString("ExcitationEnergyDistribution");
...@@ -405,13 +399,13 @@ void Reaction::ReadConfigurationFile(NPL::InputParser parser){ ...@@ -405,13 +399,13 @@ void Reaction::ReadConfigurationFile(NPL::InputParser parser){
} }
if(blocks[i]->HasToken("Shoot3")){ if(blocks[i]->HasToken("Shoot3")){
fshoot3 = blocks[i]->GetInt("Shoot3"); fshoot3 = blocks[i]->GetInt("Shoot3");
} }
if(blocks[i]->HasToken("Shoot4")){ if(blocks[i]->HasToken("Shoot4")){
fshoot4 = blocks[i]->GetInt("Shoot4"); fshoot4 = blocks[i]->GetInt("Shoot4");
} }
if(blocks[i]->HasToken("ShootHeavy")){ if(blocks[i]->HasToken("ShootHeavy")){
fshoot4 = blocks[i]->GetInt("ShootHeavy"); fshoot4 = blocks[i]->GetInt("ShootHeavy");
} }
if(blocks[i]->HasToken("ShootLight")){ if(blocks[i]->HasToken("ShootLight")){
fshoot3 = blocks[i]->GetInt("ShootLight"); fshoot3 = blocks[i]->GetInt("ShootLight");
} }
...@@ -435,6 +429,7 @@ void Reaction::initializePrecomputeVariable(){ ...@@ -435,6 +429,7 @@ void Reaction::initializePrecomputeVariable(){
s = m1*m1 + m2*m2 + 2*m2*(fBeamEnergy + m1); s = m1*m1 + m2*m2 + 2*m2*(fBeamEnergy + m1);
fTotalEnergyImpulsionCM = TLorentzVector(0,0,0,sqrt(s)); fTotalEnergyImpulsionCM = TLorentzVector(0,0,0,sqrt(s));
fEcm = sqrt(s) - m1 -m2;
ECM_1 = (s + m1*m1 - m2*m2)/(2*sqrt(s)); ECM_1 = (s + m1*m1 - m2*m2)/(2*sqrt(s));
ECM_2 = (s + m2*m2 - m1*m1)/(2*sqrt(s)); ECM_2 = (s + m2*m2 - m1*m1)/(2*sqrt(s));
...@@ -464,20 +459,15 @@ void Reaction::initializePrecomputeVariable(){ ...@@ -464,20 +459,15 @@ void Reaction::initializePrecomputeVariable(){
} }
//////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////
void Reaction::SetNuclei3(double EnergyLab, double ThetaLab, double PhiLab){ void Reaction::SetNuclei3(double EnergyLab, double ThetaLab){
double p3 = sqrt(pow(EnergyLab,2) + 2*m3*EnergyLab);
double E3 = m3 + EnergyLab;
double p_Lab_3 = sqrt(E3*E3 - m3*m3);
TVector3 p_Lab_3_vec = TVector3(0,0,1);
p_Lab_3_vec.SetMagThetaPhi(p_Lab_3, ThetaLab, PhiLab);
fEnergyImpulsionLab_3 = TLorentzVector(p_Lab_3_vec, E3); fEnergyImpulsionLab_3 = TLorentzVector(p3*sin(ThetaLab),0,p3*cos(ThetaLab),EnergyLab+m3);
fEnergyImpulsionLab_4 = fTotalEnergyImpulsionLab - fEnergyImpulsionLab_3; fEnergyImpulsionLab_4 = fTotalEnergyImpulsionLab - fEnergyImpulsionLab_3;
fNuclei3.SetEnergyImpulsion(fEnergyImpulsionLab_3); fNuclei3.SetEnergyImpulsion(fEnergyImpulsionLab_3);
fNuclei4.SetEnergyImpulsion(fEnergyImpulsionLab_4); fNuclei4.SetEnergyImpulsion(fEnergyImpulsionLab_4);
//ThetaCM and Energy do not depend on PhiLab
fThetaCM = EnergyLabToThetaCM(EnergyLab, ThetaLab); fThetaCM = EnergyLabToThetaCM(EnergyLab, ThetaLab);
fExcitation4 = ReconstructRelativistic(EnergyLab, ThetaLab); fExcitation4 = ReconstructRelativistic(EnergyLab, ThetaLab);
} }
...@@ -525,20 +515,20 @@ TGraph* Reaction::GetKinematicLine4(double AngleStep_CM){ ...@@ -525,20 +515,20 @@ TGraph* Reaction::GetKinematicLine4(double AngleStep_CM){
return(fKineLine4); return(fKineLine4);
} }
//////////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////
TGraph* Reaction::GetTheta3VsTheta4(double AngleStep_CM) TGraph* Reaction::GetTheta3VsTheta4(double AngleStep_CM)
{ {
vector<double> vx; vector<double> vx;
vector<double> vy; vector<double> vy;
double theta3,E3,theta4,E4; double theta3,E3,theta4,E4;
for (double angle=0 ; angle < 360 ; angle+=AngleStep_CM){ for (double angle=0 ; angle < 360 ; angle+=AngleStep_CM){
SetThetaCM(angle*deg); SetThetaCM(angle*deg);
KineRelativistic(theta3, E3, theta4, E4); KineRelativistic(theta3, E3, theta4, E4);
vx.push_back(theta3/deg); vx.push_back(theta3/deg);
vy.push_back(theta4/deg); vy.push_back(theta4/deg);
} }
fTheta3VsTheta4= new TGraph(vx.size(),&vx[0],&vy[0]); fTheta3VsTheta4= new TGraph(vx.size(),&vx[0],&vy[0]);
return(fTheta3VsTheta4); return(fTheta3VsTheta4);
...@@ -644,4 +634,3 @@ void Reaction::SetCSAngle(double CSHalfOpenAngleMin,double CSHalfOpenAngleMax){ ...@@ -644,4 +634,3 @@ void Reaction::SetCSAngle(double CSHalfOpenAngleMin,double CSHalfOpenAngleMax){
fCrossSectionHist->SetBinContent(i,0); fCrossSectionHist->SetBinContent(i,0);
} }
} }
NPLib/Physics/NPReaction.h
+ 8
− 5
View file @ 309fd4a8
...@@ -82,6 +82,7 @@ namespace NPL{ ...@@ -82,6 +82,7 @@ namespace NPL{
Nucleus fNuclei3; // Light ejectile Nucleus fNuclei3; // Light ejectile
Nucleus fNuclei4; // Heavy ejectile Nucleus fNuclei4; // Heavy ejectile
double fQValue; // Q-value in MeV double fQValue; // Q-value in MeV
double fEcm; // Ecm in MeV
double fBeamEnergy; // Beam energy in MeV double fBeamEnergy; // Beam energy in MeV
double fThetaCM; // Center-of-mass angle in radian double fThetaCM; // Center-of-mass angle in radian
double fExcitation3; // Excitation energy in MeV double fExcitation3; // Excitation energy in MeV
...@@ -92,6 +93,7 @@ namespace NPL{ ...@@ -92,6 +93,7 @@ namespace NPL{
public: public:
// Getters and Setters // Getters and Setters
void SetBeamEnergy(double eBeam) {fBeamEnergy = eBeam; initializePrecomputeVariable();} void SetBeamEnergy(double eBeam) {fBeamEnergy = eBeam; initializePrecomputeVariable();}
void SetEcm(double Ecm) {fEcm= Ecm; s=pow(Ecm+m1+m2,2); fBeamEnergy=(s-m1*m1-m2*m2)/(2*m2)-m1; initializePrecomputeVariable();}
void SetThetaCM(double angle) {fThetaCM = angle; initializePrecomputeVariable();} void SetThetaCM(double angle) {fThetaCM = angle; initializePrecomputeVariable();}
void SetExcitation3(double exci) {fExcitation3 = exci; initializePrecomputeVariable();} void SetExcitation3(double exci) {fExcitation3 = exci; initializePrecomputeVariable();}
void SetExcitation4(double exci) {fExcitation4 = exci; initializePrecomputeVariable();} void SetExcitation4(double exci) {fExcitation4 = exci; initializePrecomputeVariable();}
...@@ -99,16 +101,17 @@ namespace NPL{ ...@@ -99,16 +101,17 @@ namespace NPL{
void SetExcitationLight(double exci) {fExcitation3 = exci; initializePrecomputeVariable();} void SetExcitationLight(double exci) {fExcitation3 = exci; initializePrecomputeVariable();}
void SetExcitationHeavy(double exci) {fExcitation4 = exci; initializePrecomputeVariable();} void SetExcitationHeavy(double exci) {fExcitation4 = exci; initializePrecomputeVariable();}
void SetVerboseLevel(int verbose) {fVerboseLevel = verbose;} void SetVerboseLevel(int verbose) {fVerboseLevel = verbose;}
void SetCrossSectionHist (TH1F* CrossSectionHist) void SetCrossSectionHist (TH1F* CrossSectionHist)
{delete fCrossSectionHist; fCrossSectionHist = CrossSectionHist;} {delete fCrossSectionHist; fCrossSectionHist = CrossSectionHist;}
void SetDoubleDifferentialCrossSectionHist(TH2F* CrossSectionHist) void SetDoubleDifferentialCrossSectionHist(TH2F* CrossSectionHist)
{fDoubleDifferentialCrossSectionHist = CrossSectionHist;} {fDoubleDifferentialCrossSectionHist = CrossSectionHist;}
double GetBeamEnergy() const {return fBeamEnergy;} double GetBeamEnergy() const {return fBeamEnergy;}
double GetThetaCM() const {return fThetaCM;} double GetThetaCM() const {return fThetaCM;}
double GetExcitation3() const {return fExcitation3;} double GetExcitation3() const {return fExcitation3;}
double GetExcitation4() const {return fExcitation4;} double GetExcitation4() const {return fExcitation4;}
double GetQValue() const {return fQValue;} double GetQValue() const {return fQValue;}
double GetEcm() const {return fEcm;}
Nucleus GetNucleus1() const {return fNuclei1;} Nucleus GetNucleus1() const {return fNuclei1;}
Nucleus GetNucleus2() const {return fNuclei2;} Nucleus GetNucleus2() const {return fNuclei2;}
Nucleus GetNucleus3() const {return fNuclei3;} Nucleus GetNucleus3() const {return fNuclei3;}
...@@ -192,14 +195,14 @@ namespace NPL{ ...@@ -192,14 +195,14 @@ namespace NPL{
double &ThetaLab4, double &KineticEnergyLab4); double &ThetaLab4, double &KineticEnergyLab4);
// Return Excitation Energy // Return Excitation Energy
double ReconstructRelativistic(double EnergyLab, double ThetaLab, double PhiLab=0); double ReconstructRelativistic(double EnergyLab, double ThetaLab);
// Return ThetaCM // Return ThetaCM
// EnergyLab: energy measured in the laboratory frame // EnergyLab: energy measured in the laboratory frame
// ExcitationEnergy: excitation energy previously calculated. // ExcitationEnergy: excitation energy previously calculated.
double EnergyLabToThetaCM(double EnergyLab, double ThetaLab, double PhiLab=0); double EnergyLabToThetaCM(double EnergyLab, double ThetaLab);
void SetNuclei3(double EnergyLab, double ThetaLab, double PhiLab=0); void SetNuclei3(double EnergyLab, double ThetaLab);
TGraph* GetKinematicLine3(double AngleStep_CM=1); TGraph* GetKinematicLine3(double AngleStep_CM=1);
TGraph* GetKinematicLine4(double AngleStep_CM=1); TGraph* GetKinematicLine4(double AngleStep_CM=1);
......
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