Skip to content
Snippets Groups Projects
Commit 091c5e99 authored by adrien-matta's avatar adrien-matta
Browse files

* Adding support for double differential cross section 

* The code compile and run but still need confirmation that it work as
* expected
* Modification does not affect the usual single differential way of
* working
parent 768addcd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
NPLib/Physics/NPReaction.cxx
+ 42
− 32
View file @ 091c5e99
......@@ -27,9 +27,10 @@
* *
*---------------------------------------------------------------------------*
* Comment: *
* + Redesign using LorentzVector by Pierre Morfouace *
* + 20/01/2011: Add support for excitation energy for light ejectile *
* (N. de Sereville) *
* Based on previous work by N.de Sereville *
* + Based on previous work by N.de Sereville *
* *
*****************************************************************************/
......@@ -55,7 +56,6 @@ using namespace NPUNITS;
// ROOT
#include"TF1.h"
ClassImp(Reaction)
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......@@ -89,7 +89,8 @@ Reaction::Reaction(){
fCrossSectionHist = new TH1F(Form("EnergyHist_%i",offset),"Reaction_CS",1,0,180);
fExcitationEnergyHist = NULL;
fDoubleDifferentialCrossSectionHist = NULL ;
fshoot3=true;
fshoot4=true;
......@@ -154,8 +155,9 @@ Reaction::Reaction(string reaction){
fCrossSectionHist = new TH1F(Form("EnergyHist_%i",offset),"Reaction_CS",1,0,180);
fCrossSectionHist = NULL;
fshoot3=true;
fDoubleDifferentialCrossSectionHist = NULL ;
fshoot3=true;
fshoot4=true;
initializePrecomputeVariable();
......@@ -205,7 +207,15 @@ bool Reaction::CheckKinematic(){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Reaction::ShootRandomThetaCM(){
double theta;
SetThetaCM( theta=fCrossSectionHist->GetRandom()*deg );
if(fDoubleDifferentialCrossSectionHist){
// Take a slice in energy
TAxis* Y = fDoubleDifferentialCrossSectionHist->GetYaxis();
TH1D* Proj = fDoubleDifferentialCrossSectionHist
->ProjectionY("proj",Y->FindBin(fBeamEnergy),Y->FindBin(fBeamEnergy));
SetThetaCM( theta=Proj->GetRandom()*deg );
}
else
SetThetaCM( theta=fCrossSectionHist->GetRandom()*deg );
return theta;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......@@ -252,7 +262,6 @@ void Reaction::KineRelativistic(double &ThetaLab3, double &KineticEnergyLab3,
cout << "Problem for energy conservation" << endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Reaction::ReconstructRelativistic(double EnergyLab, double ThetaLab){
// EnergyLab in MeV
......@@ -268,8 +277,6 @@ double Reaction::ReconstructRelativistic(double EnergyLab, double ThetaLab){
return Eex;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//Return ThetaCM
double Reaction::EnergyLabToThetaCM(double EnergyLab, double ThetaLab){
......@@ -316,6 +323,7 @@ void Reaction::ReadConfigurationFile(string Path){
bool check_ExcitationEnergy3 = false ;
bool check_ExcitationEnergy4 = false ;
bool check_ExcitationEnergyDistribution = false;
bool check_DoubleDifferentialCrossSectionPath = false;
bool check_CrossSectionPath = false ;
bool check_shoot3 = false ;
bool check_shoot4 = false;
......@@ -408,7 +416,7 @@ void Reaction::ReadConfigurationFile(string Path){
string FileName,HistName;
ReactionFile >> FileName >> HistName;
if(fVerboseLevel==1) cout << "Reading Cross Section file: " << FileName << endl;
TH1F* CStemp = Read1DProfile(FileName, HistName );
TH1F* CStemp = Read1DProfile(FileName, HistName);
// multiply CStemp by sin(theta)
TF1* fsin = new TF1("sin",Form("1/(sin(x*%f/180.))",M_PI),0,180);
......@@ -417,6 +425,27 @@ void Reaction::ReadConfigurationFile(string Path){
delete fsin;
}
// Use for multi energy beam simulation
// This distribution is the differential Cross section for each Beam energy
// The Beam energy distribution itself should not be included originally
else if (DataBuffer== "DoubleDifferentialCrossSectionPath=") {
check_DoubleDifferentialCrossSectionPath = true ;
string FileName,HistName;
ReactionFile >> FileName >> HistName;
if(fVerboseLevel==1) cout << "Reading Double Differential Cross Section file: " << FileName << endl;
TH2F* CStemp = Read2DProfile(FileName, HistName);
// multiply CStemp by sin(theta)
// X axis is theta CM
// Y axis is beam energy
// Division affect only X axis
TF1* fsin = new TF1("sin",Form("1/(sin(x*%f/180.))",M_PI),0,180);
CStemp->Divide(fsin,1);
SetDoubleDifferentialCrossSectionHist(CStemp);
delete fsin;
}
else if (DataBuffer.compare(0, 17, "HalfOpenAngleMin=") == 0) {
ReactionFile >> DataBuffer;
CSHalfOpenAngleMin = atof(DataBuffer.c_str()) * deg;
......@@ -459,7 +488,9 @@ void Reaction::ReadConfigurationFile(string Path){
///////////////////////////////////////////////////
// If all Token found toggle out
if (check_Beam && check_Target && check_Light && check_Heavy && check_ExcitationEnergy3
&& (check_ExcitationEnergy4 || check_ExcitationEnergyDistribution ) && check_CrossSectionPath && check_shoot4 && check_shoot3)
&& (check_ExcitationEnergy4 || check_ExcitationEnergyDistribution )
&& (check_CrossSectionPath || check_DoubleDifferentialCrossSectionPath)
&& check_shoot4 && check_shoot3)
ReadingStatus = false;
}
}
......@@ -470,7 +501,6 @@ void Reaction::ReadConfigurationFile(string Path){
initializePrecomputeVariable();
}
////////////////////////////////////////////////////////////////////////////////////////////
void Reaction::initializePrecomputeVariable(){
......@@ -524,9 +554,6 @@ void Reaction::SetNuclei3(double EnergyLab, double ThetaLab){
fExcitation4 = ReconstructRelativistic(EnergyLab, ThetaLab);
}
////////////////////////////////////////////////////////////////////////////////////////////
TGraph* Reaction::GetKinematicLine3(double AngleStep_CM){
......@@ -547,7 +574,6 @@ TGraph* Reaction::GetKinematicLine3(double AngleStep_CM){
return(fKineLine3);
}
////////////////////////////////////////////////////////////////////////////////////////////
TGraph* Reaction::GetKinematicLine4(double AngleStep_CM){
......@@ -627,8 +653,6 @@ TGraph* Reaction::GetThetaLabVersusThetaCM(double AngleStep_CM){
return(fAngleLine);
}
////////////////////////////////////////////////////////////////////////////////////////////
void Reaction::PrintKinematic(){
int size = 360;
......@@ -651,7 +675,6 @@ void Reaction::PrintKinematic(){
}
}
////////////////////////////////////////////////////////////////////////////////////////////
void Reaction::SetCSAngle(double CSHalfOpenAngleMin,double CSHalfOpenAngleMax){
......@@ -661,16 +684,3 @@ void Reaction::SetCSAngle(double CSHalfOpenAngleMin,double CSHalfOpenAngleMax){
}
NPLib/Physics/NPReaction.h
+ 5
− 1
View file @ 091c5e99
......@@ -90,6 +90,7 @@ namespace NPL{
double fExcitation3; // Excitation energy in MeV
double fExcitation4; // Excitation energy in MeV
TH1F* fCrossSectionHist; // Differential cross section in CM frame
TH2F* fDoubleDifferentialCrossSectionHist; // Diff. CS CM frame vs Beam E
TH1F* fExcitationEnergyHist; // Distribution of Excitation energy
public:
// Getters and Setters
......@@ -101,8 +102,11 @@ namespace NPL{
void SetExcitationLight(double exci) {fExcitation3 = exci; initializePrecomputeVariable();}
void SetExcitationHeavy(double exci) {fExcitation4 = exci; initializePrecomputeVariable();}
void SetVerboseLevel(int verbose) {fVerboseLevel = verbose;}
void SetCrossSectionHist (TH1F* CrossSectionHist) {delete fCrossSectionHist; fCrossSectionHist = CrossSectionHist;}
void SetCrossSectionHist (TH1F* CrossSectionHist)
{delete fCrossSectionHist; fCrossSectionHist = CrossSectionHist;}
void SetDoubleDifferentialCrossSectionHist(TH2F* CrossSectionHist)
{fDoubleDifferentialCrossSectionHist = CrossSectionHist;}
double GetBeamEnergy() const {return fBeamEnergy;}
double GetThetaCM() const {return fThetaCM;}
double GetExcitation3() const {return fExcitation3;}
......
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