diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index 93e30380c1d118b8cb038ef8806da790f2c5e4d7..0d5dd385cf1b736ecb32492b02b640ecc9746434 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -72,6 +72,7 @@ ClassImp(Reaction)
//
fBeamEnergy = 0;
fThetaCM = 0;
+ fExcitation1 = 0;
fExcitation3 = 0;
fExcitation4 = 0;
fQValue = 0;
@@ -84,6 +85,9 @@ ClassImp(Reaction)
fshoot3=true;
fshoot4=true;
+ RandGen=new TRandom3();
+
+ fLabCrossSection=false; // flag if the provided cross-section is in the lab or not
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -133,6 +137,7 @@ Reaction::Reaction(string reaction){
fNuclei4 = Nucleus(D);
fBeamEnergy = atof(E.c_str());
fThetaCM = 0;
+ fExcitation1 = 0;
fExcitation3 = 0;
fExcitation4 = 0;
fQValue = 0;
@@ -146,9 +151,12 @@ Reaction::Reaction(string reaction){
fCrossSectionHist = new TH1F(Form("EnergyHist_%i",offset),"Reaction_CS",1,0,180);
fDoubleDifferentialCrossSectionHist = NULL ;
-
+
fshoot3=true;
fshoot4=true;
+ RandGen=new TRandom3();
+
+ fLabCrossSection=false;
initializePrecomputeVariable();
}
@@ -182,7 +190,7 @@ bool Reaction::CheckKinematic(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Reaction::ShootRandomThetaCM(){
- double theta;
+ double theta; // CM
if(fDoubleDifferentialCrossSectionHist){
// Take a slice in energy
TAxis* Y = fDoubleDifferentialCrossSectionHist->GetYaxis();
@@ -203,6 +211,12 @@ double Reaction::ShootRandomThetaCM(){
TH1D* Proj = fDoubleDifferentialCrossSectionHist->ProjectionX("proj",binY,binY);
SetThetaCM( theta=Proj->GetRandom()*deg );
}
+ else if (fLabCrossSection){
+ double thetalab=fCrossSectionHist->GetRandom()*deg;//shot in lab
+ double energylab=EnergyLabFromThetaLab(thetalab); //get corresponding energy
+ theta = EnergyLabToThetaCM(energylab, thetalab); //transform to theta CM
+ SetThetaCM( theta );
+ }
else
SetThetaCM( theta=fCrossSectionHist->GetRandom()*deg );
return theta;
@@ -286,6 +300,53 @@ double Reaction::EnergyLabToThetaCM(double EnergyLab, double ThetaLab, double P
return(ThetaCM);
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//Return EnergyLab
+double Reaction::EnergyLabFromThetaLab(double ThetaLab){
+ //ThetaLab in rad
+
+ // NB
+ // Tis member function is in progress
+ // Classic treatment need to be changed to relativistic!!
+ //
+
+ //Treat Exeptions
+ if(fBeamEnergy==0) return m4*fQValue/(m3+m4);
+
+ //Solve the equation
+ double m4e = m4 + fExcitation4;
+ //Possible problem with this calculation, decide on m4e or m4 in the following formula
+ double r = sqrt(m1*m3*fBeamEnergy)/(m3+m4e);
+ double s = ((m1-m4e)*fBeamEnergy-m4e*fQValue)/(m3+m4e);
+
+ //Solve quadratic euation of T3 : T3-(2*r*cos(thetaLab)*sqrt(T3)+s = 0
+ //
+ //Solution: sqrt(T3) = r * (cos +/- sqrt(Delta)) ;
+ //where Delta =cos*cos - s/(r*r)
+
+ //Select solutions according to delta values
+ double Delta = cos(ThetaLab)*cos(ThetaLab)-s/(r*r);
+ double EnergyLab=0;
+
+ if(Delta<0) return 0;
+ else
+ if(Delta==0) EnergyLab =r*cos(ThetaLab);
+ else {
+ double sol1 = r*(cos(ThetaLab) + sqrt(Delta));
+ double sol2 = r*(cos(ThetaLab) - sqrt(Delta));
+
+ if(sol1<0 && sol2<0) return 0;
+ else if(sol1>0 && sol2<0) EnergyLab =sol1*sol1;
+ else if(sol1<0 && sol2>0) EnergyLab = sol2*sol2;
+ else {
+ if(RandGen->Rndm()<0.5) EnergyLab = sol1*sol1; // choose either of the solutions
+ else EnergyLab = sol2*sol2;
+ }
+ }
+
+ return(EnergyLab);
+}
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Reaction::Print() const{
@@ -296,6 +357,7 @@ void Reaction::Print() const{
<< fBeamEnergy << " MeV"
<< endl ;
+ cout << "Exc Nuclei 1 = " << fExcitation1 << " MeV" << endl;
cout << "Exc Nuclei 3 = " << fExcitation3 << " MeV" << endl;
cout << "Exc Nuclei 4 = " << fExcitation4 << " MeV" << endl;
cout << "Qgg = " << fQValue << " MeV" << endl;
@@ -355,6 +417,13 @@ void Reaction::ReadConfigurationFile(NPL::InputParser parser){
exit(1);
}
+ if(blocks[i]->HasToken("ExcitationEnergyBeam")){
+ fExcitation1 = blocks[i]->GetDouble("ExcitationEnergyBeam","MeV");
+ }
+ else if(blocks[i]->HasToken("ExcitationEnergy1")){
+ fExcitation1 = blocks[i]->GetDouble("ExcitationEnergy1","MeV");
+ }
+
if(blocks[i]->HasToken("ExcitationEnergyLight"))
fExcitation3 = blocks[i]->GetDouble("ExcitationEnergyLight","MeV");
else if(blocks[i]->HasToken("ExcitationEnergy3"))
@@ -382,6 +451,20 @@ void Reaction::ReadConfigurationFile(NPL::InputParser parser){
delete fsin;
}
+ if(blocks[i]->HasToken("LabCrossSectionPath")){
+ fLabCrossSection=true;
+
+ vector<string> file = blocks[i]->GetVectorString("LabCrossSectionPath");
+ TH1F* CStemp = Read1DProfile(file[0], file[1]);
+
+ // multiply CStemp by sin(theta)
+ TF1* fsin = new TF1("sin",Form("1/(sin(x*%f/180.))",M_PI),0,180);
+ CStemp->Divide(fsin,1);
+ SetCrossSectionHist(CStemp);
+ delete fsin;
+ }
+
+
if(blocks[i]->HasToken("DoubleDifferentialCrossSectionPath")){
vector<string> file = blocks[i]->GetVectorString("DoubleDifferentialCrossSectionPath");
TH2F* CStemp = Read2DProfile(file[0],file[1]);
@@ -427,7 +510,10 @@ void Reaction::initializePrecomputeVariable(){
if(fBeamEnergy < 0)
fBeamEnergy = 0 ;
- m1 = fNuclei1.Mass();
+ if(fExcitation1>=0) fNuclei1.SetExcitationEnergy(fExcitation1); // Write over the beam excitation energy
+
+ //fNuclei1.GetExcitationEnergy() is a copy of fExcitation1
+ m1 = fNuclei1.Mass() + fNuclei1.GetExcitationEnergy();// in case of isomeric state,
m2 = fNuclei2.Mass();
m3 = fNuclei3.Mass() + fExcitation3;
m4 = fNuclei4.Mass() + fExcitation4;
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index 2d70cb1dec26a35df7dea3597b319e7361cf0e56..20d9c2e0e4348b0a693eb15a5e209fd7bcc621f7 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -44,6 +44,7 @@ using namespace NPL;
#include "TGraph.h"
#include "TCanvas.h"
#include "TH1F.h"
+#include "TRandom3.h"
using namespace std;
@@ -62,13 +63,16 @@ namespace NPL{
void ReadConfigurationFile(string Path);
void ReadConfigurationFile(NPL::InputParser);
-
+ private:
+ TRandom3* RandGen;
+
private:
int fVerboseLevel;
private: // use for Monte Carlo simulation
bool fshoot3;
bool fshoot4;
+ bool fLabCrossSection;
private: // use to display the kinematical line
TGraph* fKineLine3 ;
@@ -84,6 +88,7 @@ namespace NPL{
double fQValue; // Q-value in MeV
double fBeamEnergy; // Beam energy in MeV
double fThetaCM; // Center-of-mass angle in radian
+ double fExcitation1; // Excitation energy in MeV of the beam, useful for isomers
double fExcitation3; // Excitation energy in MeV
double fExcitation4; // Excitation energy in MeV
TH1F* fCrossSectionHist; // Differential cross section in CM frame
@@ -93,9 +98,11 @@ namespace NPL{
// Getters and Setters
void SetBeamEnergy(double eBeam) {fBeamEnergy = eBeam; initializePrecomputeVariable();}
void SetThetaCM(double angle) {fThetaCM = angle; initializePrecomputeVariable();}
+ void SetExcitation1(double exci) {fExcitation1 = exci; initializePrecomputeVariable();}
void SetExcitation3(double exci) {fExcitation3 = exci; initializePrecomputeVariable();}
void SetExcitation4(double exci) {fExcitation4 = exci; initializePrecomputeVariable();}
// For retro compatibility
+ void SetExcitationBeam(double exci) {fExcitation1 = exci; initializePrecomputeVariable();}
void SetExcitationLight(double exci) {fExcitation3 = exci; initializePrecomputeVariable();}
void SetExcitationHeavy(double exci) {fExcitation4 = exci; initializePrecomputeVariable();}
void SetVerboseLevel(int verbose) {fVerboseLevel = verbose;}
@@ -106,6 +113,7 @@ namespace NPL{
{fDoubleDifferentialCrossSectionHist = CrossSectionHist;}
double GetBeamEnergy() const {return fBeamEnergy;}
double GetThetaCM() const {return fThetaCM;}
+ double GetExcitation1() const {return fExcitation1;}
double GetExcitation3() const {return fExcitation3;}
double GetExcitation4() const {return fExcitation4;}
double GetQValue() const {return fQValue;}
@@ -120,7 +128,7 @@ namespace NPL{
public:
- // Modify the CS histo to so cross section shoot is within ]HalfOpenAngleMin,HalfOpenAngleMax[
+ // Modify the CS histo so cross section shoot is within ]HalfOpenAngleMin,HalfOpenAngleMax[
void SetCSAngle(double CSHalfOpenAngleMin,double CSHalfOpenAngleMax);
private: // intern precompute variable
@@ -180,6 +188,7 @@ namespace NPL{
public: // Kinematics
// Check that the reaction is alowed
bool CheckKinematic();
+ bool IsLabCrossSection(){return fLabCrossSection;};
// Use fCrossSectionHist to shoot a Random ThetaCM and set fThetaCM to this value
double ShootRandomThetaCM();
@@ -198,6 +207,10 @@ namespace NPL{
// EnergyLab: energy measured in the laboratory frame
// ExcitationEnergy: excitation energy previously calculated.
double EnergyLabToThetaCM(double EnergyLab, double ThetaLab, double PhiLab=0);
+ // Return theoretical EnergyLab, useful for a random distribution in the lab frame
+ // ThetaLab: angle measured in the laboratory frame
+ // This uses the fExcitation4 and fQValue both set previously
+ double EnergyLabFromThetaLab(double ThetaLab);
void SetNuclei3(double EnergyLab, double ThetaLab, double PhiLab=0);