diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction
index 9e39063fe67f67b66cdf9882d2fdb7ae761af4a4..280f56a4f04ca961f55b02431f62cbead8e13941 100644
--- a/Inputs/EventGenerator/10He.reaction
+++ b/Inputs/EventGenerator/10He.reaction
@@ -17,7 +17,7 @@ Beam
%XThetaXProfilePath=
%YPhiYProfilePath=
-%%Beam energy given in MeV ; Excitation in MeV
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TwoBodyReaction
Beam= 11Li
@@ -26,29 +26,34 @@ TwoBodyReaction
Heavy= 10He
ExcitationEnergyLight= 0.0
ExcitationEnergyHeavy= 10.0
- CrossSectionPath= 11Li(d,3He)10He.txt
+ CrossSectionPath= 11Li(d,3He)10He.txt CS10He
ShootLight= 1
ShootHeavy= 1
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ParticleDecay 10He
- Daughter= 9He n
- ExcitationEnergy= 0.5 0
- DifferentialCrossSection= flat.txt
- shoot= 1 1
+ Daughter= 9He n
+ ExcitationEnergy= 0.5 0
+ DifferentialCrossSection= flat.txt
+ shoot= 1 1
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
GammaDecay 9He
- Cascade
- BranchingRatio= 30
- Energies= 0.1
- DifferentialCrossSection= flat.txt
- Cascade
- BranchingRatio= 70
- Energies= 0.2 0.3
+ Cascade
+ BranchingRatio= 30
+ Energies= 0.1
+ DifferentialCrossSection= flat.txt
+ Cascade
+ BranchingRatio= 70
+ Energies= 0.2 0.3
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ParticleDecay 9He
- Daughter= 8He n
- DifferentialCrossSection= flat.txt
- shoot= 1 1
+ Daughter= 8He n
+ DifferentialCrossSection= flat.txt
+ shoot= 1 1
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/NPLib/Physics/NPBeam.cxx b/NPLib/Physics/NPBeam.cxx
index f0eb5953182a3ca87e4a61a1efd00d5857c903a7..85ae58a0dc03d84bdc2f8610ba36d178cec5e502 100644
--- a/NPLib/Physics/NPBeam.cxx
+++ b/NPLib/Physics/NPBeam.cxx
@@ -33,6 +33,7 @@
// NPL header
#include "NPBeam.h"
#include "NPFunction.h"
+#include "NPOptionManager.h"
// Use CLHEP System of unit and Physical Constant
#include "CLHEP/Units/GlobalSystemOfUnits.h"
@@ -60,7 +61,7 @@ Beam::Beam()
fEffectiveTargetThickness = 0 ;
fTargetAngle = 0 ;
fTargetZ = 0 ;
-
+ fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
// case of user given distribution
fEnergyHist = new TH1F("EnergyHist","EnergyHist",1,0,1);
fXThetaXHist = new TH2F("XThetaXHis","XThetaXHis",1,0,1,1,0,1);
@@ -115,7 +116,7 @@ void Beam::ReadConfigurationFile(string Path){
getline(ReactionFile, LineBuffer);
if (LineBuffer.compare(0, 4, "Beam") == 0) {
- cout << "Beam Found" << endl ;
+ if(fVerboseLevel==1) cout << "Beam Found" << endl ;
ReadingStatus = true ;
}
@@ -131,83 +132,84 @@ void Beam::ReadConfigurationFile(string Path){
ReactionFile >> DataBuffer;
delete fBeamNucleus;
fBeamNucleus = new Nucleus(DataBuffer);
- cout << "Beam Particle: " << fBeamNucleus->GetName() << endl;
+ if(fVerboseLevel==1) cout << "Beam Particle: " << fBeamNucleus->GetName() << endl;
}
else if (DataBuffer == "Energy=") {
check_Energy = true ;
ReactionFile >> DataBuffer;
fEnergy = atof(DataBuffer.c_str()) * MeV;
- cout << "Beam Energy: " << fEnergy / MeV << " MeV" << endl;
+ if(fVerboseLevel==1) cout << "Beam Energy: " << fEnergy / MeV << " MeV" << endl;
}
else if (DataBuffer == "SigmaEnergy=") {
check_SigmaEnergy = true ;
ReactionFile >> DataBuffer;
fSigmaEnergy= atof(DataBuffer.c_str()) * MeV;
- cout << "Beam Energy Sigma: " << fSigmaEnergy / MeV << " MeV" << endl;
+ if(fVerboseLevel==1) cout << "Beam Energy Sigma: " << fSigmaEnergy / MeV << " MeV" << endl;
}
else if (DataBuffer=="MeanX=") {
check_MeanX = true ;
ReactionFile >> DataBuffer;
fMeanX = atof(DataBuffer.c_str()) * mm;
- cout << "Mean X: " << fMeanX / mm << " mm" << endl;
+ if(fVerboseLevel==1) cout << "Mean X: " << fMeanX / mm << " mm" << endl;
}
else if (DataBuffer=="MeanY=") {
check_MeanY = true ;
ReactionFile >> DataBuffer;
fMeanY = atof(DataBuffer.c_str()) * mm;
- cout << "Mean Y: " << fMeanY / mm << " mm" << endl;
+ if(fVerboseLevel==1) cout << "Mean Y: " << fMeanY / mm << " mm" << endl;
}
else if (DataBuffer=="SigmaX=") {
check_SigmaX = true ;
ReactionFile >> DataBuffer;
fSigmaX = atof(DataBuffer.c_str()) * mm;
- cout << "Sigma X: " << fSigmaX / mm << " mm" << endl;
+ if(fVerboseLevel==1) cout << "Sigma X: " << fSigmaX / mm << " mm" << endl;
}
else if (DataBuffer=="SigmaY=") {
check_SigmaY = true ;
ReactionFile >> DataBuffer;
fSigmaY = atof(DataBuffer.c_str()) * mm;
- cout << "Sigma Y: " << fSigmaY / mm << " mm" << endl;
+ if(fVerboseLevel==1) cout << "Sigma Y: " << fSigmaY / mm << " mm" << endl;
}
else if (DataBuffer == "MeanThetaX=" ) {
check_MeanThetaX = true ;
ReactionFile >> DataBuffer;
fMeanThetaX = atof(DataBuffer.c_str()) * deg;
- cout << "Mean Theta X: " << fMeanThetaX / deg << " deg" << endl;
+ if(fVerboseLevel==1) cout << "Mean Theta X: " << fMeanThetaX / deg << " deg" << endl;
}
else if (DataBuffer == "MeanPhiY=") {
check_MeanPhiY = true ;
ReactionFile >> DataBuffer;
fMeanPhiY = atof(DataBuffer.c_str()) * deg;
- cout << "Mean Phi Y: " << fMeanPhiY / deg << " deg" << endl;
+ if(fVerboseLevel==1) cout << "Mean Phi Y: " << fMeanPhiY / deg << " deg" << endl;
}
else if (DataBuffer == "SigmaThetaX=" ) {
check_SigmaThetaX = true ;
ReactionFile >> DataBuffer;
fSigmaThetaX = atof(DataBuffer.c_str()) * deg;
- cout << "Sigma Theta X: " << fSigmaThetaX / deg << " deg" << endl;
+ if(fVerboseLevel==1) cout << "Sigma Theta X: " << fSigmaThetaX / deg << " deg" << endl;
}
else if (DataBuffer == "SigmaPhiY=") {
check_SigmaPhiY = true ;
ReactionFile >> DataBuffer;
fSigmaPhiY = atof(DataBuffer.c_str()) * deg;
- cout << "Sigma Phi Y: " << fSigmaPhiY / deg << " deg" << endl;
+ if(fVerboseLevel==1) cout << "Sigma Phi Y: " << fSigmaPhiY / deg << " deg" << endl;
}
else if (DataBuffer == "EnergyProfilePath=") {
check_EnergyProfilePath = true ;
string FileName,HistName;
ReactionFile >> FileName >> HistName;
+ if(fVerboseLevel==1) cout << "Reading Energy profile file: " << FileName << endl;
fEnergyHist = Read1DProfile(FileName, HistName );
}
@@ -215,6 +217,7 @@ void Beam::ReadConfigurationFile(string Path){
check_XThetaXPath = true ;
string FileName,HistName;
ReactionFile >> FileName >> HistName;
+ if(fVerboseLevel==1) cout << "Reading X-ThetaX profile file: " << FileName << endl;
fXThetaXHist = Read2DProfile(FileName, HistName );
}
@@ -222,6 +225,7 @@ void Beam::ReadConfigurationFile(string Path){
check_YPhiYPath = true ;
string FileName,HistName;
ReactionFile >> FileName >> HistName;
+ if(fVerboseLevel==1) cout << "Reading Y-ThetaY profile file: " << FileName << endl;
fYPhiYHist = Read2DProfile(FileName, HistName );
}
diff --git a/NPLib/Physics/NPBeam.h b/NPLib/Physics/NPBeam.h
index 622c4978343241d020c79768286455b96640046a..e3180879923152c79d1f205da0491d7a9807a457 100755
--- a/NPLib/Physics/NPBeam.h
+++ b/NPLib/Physics/NPBeam.h
@@ -48,7 +48,10 @@ namespace NPL{
public: // Various Method
void ReadConfigurationFile(string Path);
-
+
+ private:
+ int fVerboseLevel;
+
private:
Nucleus* fBeamNucleus;
double fEnergy;
@@ -84,7 +87,8 @@ namespace NPL{
void SetEnergyHist (TH1F* EnergyHist) {delete fEnergyHist; fEnergyHist = new TH1F(*EnergyHist);}
void SetXThetaXHist (TH2F* XThetaXHist) {delete fXThetaXHist; fXThetaXHist = new TH2F(*XThetaXHist);}
void SetYPhiYHist (TH2F* YPhiYHist) {delete fYPhiYHist; fYPhiYHist = new TH2F(*YPhiYHist);}
-
+ void SetVerboseLevel(int verbose) {fVerboseLevel = verbose;}
+
// Get
Nucleus* GetNucleus () const {return fBeamNucleus;}
double GetEnergy () const {return fEnergy;}
@@ -100,7 +104,8 @@ namespace NPL{
TH1F* GetEnergyHist () const {return fEnergyHist;}
TH2F* GetXThetaXHist () const {return fXThetaXHist;}
TH2F* GetYPhiYHist () const {return fYPhiYHist;}
-
+ int GetVerboseLevel() const {return fVerboseLevel;}
+
private: // Event Generation private variable
double fTargetSize;
double fEffectiveTargetSize; // target size has seen from the beam axis
diff --git a/NPLib/Physics/NPFunction.cxx b/NPLib/Physics/NPFunction.cxx
index 72ca6e7a12333d8245a6bcb5eba0db194a2275b8..9e3fc528aa87680c89590cbb2c50318faa578072 100644
--- a/NPLib/Physics/NPFunction.cxx
+++ b/NPLib/Physics/NPFunction.cxx
@@ -35,7 +35,6 @@ TH1F* Read1DProfile(string filename,string HistName){
bool type = OpenASCIIorROOTFile(filename, ASCII , ROOT);
-
// ASCII File case
if(type){
string LineBuffer;
@@ -57,7 +56,7 @@ TH1F* Read1DProfile(string filename,string HistName){
w.push_back(wb);
}
}
-
+
// Look for the step size, min and max of the distribution
double min = 0 ;
double max = 0 ;
@@ -70,7 +69,7 @@ TH1F* Read1DProfile(string filename,string HistName){
max = x[0] ;
previous = x[0] ;
}
-
+
for(unsigned int i = 0 ; i < size ; i++){
if(x[i] > max) max = x[i] ;
if(x[i] < min) min = x[i] ;
@@ -79,12 +78,12 @@ TH1F* Read1DProfile(string filename,string HistName){
previous = x[i] ;
}
-
h = new TH1F(HistName.c_str(),HistName.c_str(),step,min,max);
-
+
for(unsigned int i = 0 ; i < size ; i++){
h->Fill(x[i],w[i]);
}
+
}
// ROOT File case
@@ -95,6 +94,7 @@ TH1F* Read1DProfile(string filename,string HistName){
exit(1);
}
}
+
return h;
}
@@ -199,15 +199,21 @@ bool OpenASCIIorROOTFile(string filename, ifstream &ASCII , TFile &ROOT){
// look for .root extension
size_t pos = filename.find(".root");
+ string GlobalPath = getenv("NPTOOL");
+ string StandardPath = GlobalPath + "/Inputs/CrossSection/" + filename;
+
// extension not found, file is assume to be a ASCII file
if(pos > filename.size()) {
ASCII.open(filename.c_str());
if(!ASCII.is_open()){
- cout << "Error, profile file " << filename << " not found " << endl ;
- exit(1);
-
+ ASCII.open(StandardPath.c_str());
+ if(!ASCII.is_open()){
+ cout << "Error, file " << filename << " not found " << endl ;
+ exit(1);
+ }
}
+
return true;
}
diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index e1308fdda082ebe1c32d13d8a49ca2a5f0d6fab7..b7336b51bfbefafdb4a804a8f83591b6d07953ea 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -18,9 +18,9 @@
* (2003 nuclear table of isotopes mass). *
* *
* KineRelativistic: Used in NPSimulation *
- * A relativistic calculation is made to compute Light and Heavy nuclei *
+ * A relativistic calculation is made to compute Light and Heavy nuclei *
* angle given the Theta CM angle. *
- * *
+ * *
* ReconstructRelativistic: Used in NPAnalysis *
* A relativistic calculation is made to compute Excitation energy given the*
* light angle and energy in Lab frame. *
@@ -32,17 +32,19 @@
* Based on previous work by N.de Sereville *
* *
*****************************************************************************/
-
+
#include <iostream>
#include <fstream>
#include <string>
#include <cstdlib>
#include <sstream>
-//#include <stdlib.h>
#include <cmath>
#include <vector>
#include "NPReaction.h"
+#include "NPBeam.h"
+#include "NPOptionManager.h"
+#include "NPFunction.h"
// Use CLHEP System of unit and Physical Constant
#include "CLHEP/Units/GlobalSystemOfUnits.h"
@@ -52,125 +54,38 @@
using namespace NPL;
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-Reaction::Reaction()
-{
- //------------- Default Constructor -------------
-
- fNuclei1 = new Nucleus();
- fNuclei2 = new Nucleus();
- fNuclei3 = new Nucleus();
- fNuclei4 = new Nucleus();
- fBeamEnergy = 0;
- fThetaCM = 0;
- fExcitation3 = 0;
- fExcitation4 = 0;
- fQValue = 0;
- fCrossSectionAngleMin = 0;
- fCrossSectionAngleMax = 180;
- initializePrecomputeVariable();
- fCrossSectionHist = new TH1F("Reaction_CS","Reaction_CS",180,0,180);
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-Reaction::Reaction(string name1, string name2, string name3, string name4, double BeamEnergy, double ExcitationEnergyLight, double ExcitationEnergyHeavy ,string Path, double CSThetaMin, double CSThetaMax)
-{
- SetEveryThing( name1, name2, name3, name4, BeamEnergy, ExcitationEnergyLight, ExcitationEnergyHeavy, Path, CSThetaMin, CSThetaMax);
-}
-
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void Reaction::SetEveryThing(string name1, string name2, string name3, string name4, double BeamEnergy, double ExcitationEnergyLight, double ExcitationEnergyHeavy, string Path, double CSThetaMin, double CSThetaMax)
-{
- //------------- Constructor with nuclei names and beam energy ------------
-
- fNuclei1 = new Nucleus(name1);
- fNuclei2 = new Nucleus(name2);
- fNuclei3 = new Nucleus(name3);
- fNuclei4 = new Nucleus(name4);
- fBeamEnergy = BeamEnergy;
- fThetaCM = 0;
- fExcitation3 = ExcitationEnergyLight;
- fExcitation4 = ExcitationEnergyHeavy;
- fQValue = (fNuclei1->GetMassExcess() + fNuclei2->GetMassExcess()
- - fNuclei3->GetMassExcess() - fNuclei4->GetMassExcess()) / 1000;
- fCrossSectionHist = new TH1F("Reaction_CS","Reaction_CS",180,0,180);
-
-
- int masse = fNuclei1->GetA() + fNuclei2->GetA() - fNuclei3->GetA() - fNuclei4->GetA();
- int charge = fNuclei1->GetZ() + fNuclei2->GetZ() - fNuclei3->GetZ() - fNuclei4->GetZ();
- if (masse || charge) {
- cout << endl;
- cout << "********************************** Error **********************************" << endl;
- cout << "* NPReaction: charge of mass not conserved. Check you event generator file *" << endl;
- cout << "***************************************************************************************" << endl;
- exit(1);
- }
-
- ///Read the differential cross section
- string GlobalPath = getenv("NPTOOL");
- string StandardPath = GlobalPath + "/Inputs/CrossSection/" + Path;
- ifstream CSFile;
- CSFile.open( StandardPath.c_str() );
-
- if(CSFile.is_open()) cout << "Reading Cross Section File " << Path << endl;
-
- // In case the file is not found in the standard path, the programm try to interpret the file name as an absolute or relative file path.
- else
- {
- CSFile.open( Path.c_str() );
- if(CSFile.is_open()) { cout << "Reading Cross Section File " << Path << endl;}
-
- else {cout << "Cross Section File " << Path << " not found" << endl;return;}
- }
-
- double CSBuffer,AngleBuffer;
- vector<double> CrossSectionBuffer;
- double thetamin = 200;
- double thetamax = -10;
- while(!CSFile.eof()) {
- CSFile >> AngleBuffer;
- CSFile >> CSBuffer;
- // only treat angular range defined by user
- if (AngleBuffer < CSThetaMin || AngleBuffer > CSThetaMax) continue;
- double CSFinal = CSBuffer*sin(AngleBuffer*deg);
- CrossSectionBuffer.push_back(CSFinal);
- fCrossSectionHist->Fill(AngleBuffer,CSFinal);
- // determine theta min and max
- if (AngleBuffer < thetamin) thetamin = AngleBuffer;
- if (AngleBuffer > thetamax) thetamax = AngleBuffer;
- }
-
- // set theta min and max values
- fCrossSectionAngleMin = thetamin;
- fCrossSectionAngleMax = thetamax;
-
- CSFile.close();
- fCrossSectionSize = CrossSectionBuffer.size();
- fCrossSection = new double[fCrossSectionSize] ;
- for(int i = 0 ; i < fCrossSectionSize ; i++ ){
- fCrossSection[i] = CrossSectionBuffer[i];
- }
- initializePrecomputeVariable();
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+Reaction::Reaction(){
+ //------------- Default Constructor -------------
+
+ fNuclei1 = new Nucleus();
+ fNuclei2 = new Nucleus();
+ fNuclei3 = new Nucleus();
+ fNuclei4 = new Nucleus();
+ fBeamEnergy = 0;
+ fThetaCM = 0;
+ fExcitation3 = 0;
+ fExcitation4 = 0;
+ fQValue = 0;
+ fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
+ initializePrecomputeVariable();
+ fCrossSectionHist = new TH1F("Reaction_CS","Reaction_CS",180,0,180);
+ fshoot3=true;
+ fshoot4=true;
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-Reaction::~Reaction()
-{
- //------------- Default Destructor ------------
-
- delete fNuclei1;
- delete fNuclei2;
- delete fNuclei3;
- delete fNuclei4;
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+Reaction::~Reaction(){
+ //------------- Default Destructor ------------
+
+ delete fNuclei1;
+ delete fNuclei2;
+ delete fNuclei3;
+ delete fNuclei4;
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-bool Reaction::CheckKinematic()
-{
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+bool Reaction::CheckKinematic(){
double theta = fThetaCM;
if (m1 > m2) theta = M_PI - fThetaCM;
@@ -186,7 +101,7 @@ bool Reaction::CheckKinematic()
cout << "Problem with energy conservation" << endl;
return false;
}
-
+
else{
//cout << "Kinematic OK" << endl;
return true;
@@ -199,18 +114,17 @@ double Reaction::ShootRandomThetaCM(){
SetThetaCM( theta=fCrossSectionHist->GetRandom()*deg );
return theta;
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Reaction::KineRelativistic(double &ThetaLab3, double &KineticEnergyLab3,
- double &ThetaLab4, double &KineticEnergyLab4)
-{
+ double &ThetaLab4, double &KineticEnergyLab4){
// 2-body relativistic kinematics: direct + inverse
// EnergieLab3,4 : lab energy in MeV of the 2 ejectiles
// ThetaLab3,4 : angles in rad
-
+
// case of inverse kinematics
double theta = fThetaCM;
if (m1 > m2) theta = M_PI - fThetaCM;
-
+
fEnergyImpulsionCM_3 = TLorentzVector(pCM_3*sin(theta),0,pCM_3*cos(theta),ECM_3);
fEnergyImpulsionCM_4 = fTotalEnergyImpulsionCM - fEnergyImpulsionCM_3;
@@ -222,7 +136,7 @@ void Reaction::KineRelativistic(double &ThetaLab3, double &KineticEnergyLab3,
// Angle in the lab frame
ThetaLab3 = fEnergyImpulsionLab_3.Angle(fEnergyImpulsionLab_1.Vect());
if (ThetaLab3 < 0) ThetaLab3 += M_PI;
-
+
ThetaLab4 = fEnergyImpulsionLab_4.Angle(fEnergyImpulsionLab_1.Vect());
if (fabs(ThetaLab3) < 1e-6) ThetaLab3 = 0;
ThetaLab4 = fabs(ThetaLab4);
@@ -233,16 +147,15 @@ void Reaction::KineRelativistic(double &ThetaLab3, double &KineticEnergyLab3,
KineticEnergyLab4 = fEnergyImpulsionLab_4.E() - m4;
// test for total energy conversion
- if (fabs(fTotalEnergyImpulsionLab.E() - (fEnergyImpulsionLab_3.E()+fEnergyImpulsionLab_4.E())) > 1e-6)
+ if (fabs(fTotalEnergyImpulsionLab.E() - (fEnergyImpulsionLab_3.E()+fEnergyImpulsionLab_4.E())) > 1e-6)
cout << "Problem for energy conservation" << endl;
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-double Reaction::ReconstructRelativistic(double EnergyLab, double ThetaLab)
-{
- // EnergyLab in MeV
- // ThetaLab in rad
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+double Reaction::ReconstructRelativistic(double EnergyLab, double ThetaLab){
+ // EnergyLab in MeV
+ // ThetaLab in rad
double E3 = m3 + EnergyLab;
double p_Lab_3 = sqrt(E3*E3 - m3*m3);
@@ -250,227 +163,242 @@ double Reaction::ReconstructRelativistic(double EnergyLab, double ThetaLab)
fEnergyImpulsionLab_4 = fTotalEnergyImpulsionLab - fEnergyImpulsionLab_3;
double Eex = fEnergyImpulsionLab_4.Mag() - fNuclei4->Mass();
-
+
return Eex;
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//Return ThetaCM
-double Reaction::EnergyLabToThetaCM(double EnergyLab, double ThetaLab)
- {
- double E3 = m3 + EnergyLab;
- double p_Lab_3 = sqrt(E3*E3 - m3*m3);
-
- fEnergyImpulsionLab_3 = TLorentzVector(p_Lab_3*sin(ThetaLab),0,p_Lab_3*cos(ThetaLab),E3);
- fEnergyImpulsionCM_3 = fEnergyImpulsionLab_3;
- fEnergyImpulsionCM_3.Boost(0,0,-BetaCM);
-
- double ThetaCM = CLHEP::pi - fEnergyImpulsionCM_1.Angle(fEnergyImpulsionCM_3.Vect());
-
- return(ThetaCM);
- }
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void Reaction::Print() const
- {
- // Print informations concerning the reaction
-
- cout << "Reaction : " << fNuclei2->GetName() << "(" << fNuclei1->GetName()
- << "," << fNuclei3->GetName() << ")" << fNuclei4->GetName() << " @ "
- << fBeamEnergy << " MeV"
- << endl ;
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//Return ThetaCM
+double Reaction::EnergyLabToThetaCM(double EnergyLab, double ThetaLab){
+ double E3 = m3 + EnergyLab;
+ double p_Lab_3 = sqrt(E3*E3 - m3*m3);
+
+ fEnergyImpulsionLab_3 = TLorentzVector(p_Lab_3*sin(ThetaLab),0,p_Lab_3*cos(ThetaLab),E3);
+ fEnergyImpulsionCM_3 = fEnergyImpulsionLab_3;
+ fEnergyImpulsionCM_3.Boost(0,0,-BetaCM);
+
+ double ThetaCM = CLHEP::pi - fEnergyImpulsionCM_1.Angle(fEnergyImpulsionCM_3.Vect());
+
+ return(ThetaCM);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Reaction::Print() const{
+ // Print informations concerning the reaction
+
+ cout << "Reaction : " << fNuclei2->GetName() << "(" << fNuclei1->GetName()
+ << "," << fNuclei3->GetName() << ")" << fNuclei4->GetName() << " @ "
+ << fBeamEnergy << " MeV"
+ << endl ;
+
+ cout << "Exc Nuclei 3 = " << fExcitation3 << " MeV" << endl;
+ cout << "Exc Nuclei 4 = " << fExcitation4 << " MeV" << endl;
+ cout << "Qgg = " << fQValue << " MeV" << endl;
+}
+
+
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////
+void Reaction::ReadConfigurationFile(string Path){
+ ////////General Reading needs////////
+ string LineBuffer;
+ string DataBuffer;
+
+ ////////Reaction Setting needs///////
+ string Beam, Target, Heavy, Light, CrossSectionPath ;
+ double ExcitationEnergy3 , ExcitationEnergy4 ;
+ double CSHalfOpenAngleMin = 0, CSHalfOpenAngleMax = 0;
+ bool ReadingStatus = false ;
+ bool check_Beam = false ;
+ bool check_Target = false ;
+ bool check_Light = false ;
+ bool check_Heavy = false ;
+ bool check_ExcitationEnergy3 = false ;
+ bool check_ExcitationEnergy4 = false ;
+ bool check_CrossSectionPath = false ;
+ bool check_shoot3 = false ;
+ bool check_shoot4 = false;
+
+ //////////////////////////////////////////////////////////////////////////////////////////
+ ifstream ReactionFile;
+ string GlobalPath = getenv("NPTOOL");
+ string StandardPath = GlobalPath + "/Inputs/EventGenerator/" + Path;
+ ReactionFile.open(StandardPath.c_str());
+ if (ReactionFile.is_open()) {cout << "Reading Reaction File " << Path << endl ;}
+
+ // In case the file is not found in the standard path, the programm try to interpret the file name as an absolute or relative file path.
+ else{
+ ReactionFile.open( Path.c_str() );
+ if(ReactionFile.is_open()) { if(fVerboseLevel==1) cout << "Reading Reaction File " << Path << endl;}
+
+ else {cout << "Reaction File " << Path << " not found" << endl;exit(1);}
+ }
+
+ while (!ReactionFile.eof()) {
+ //Pick-up next line
+ getline(ReactionFile, LineBuffer);
+
+ if (LineBuffer.compare(0, 15, "TwoBodyReaction") == 0) { ReadingStatus = true ;}
+
+
+ while(ReadingStatus){
+
+ ReactionFile >> DataBuffer;
+
+ //Search for comment Symbol %
+ if (LineBuffer.compare(0, 1, "%") == 0) {/* Do Nothing */;}
+
+ else if (DataBuffer=="Beam=") {
+ check_Beam = true ;
+ ReactionFile >> DataBuffer;
+ fNuclei1 = new Nucleus(DataBuffer);
+ Beam = DataBuffer;
+ if(fVerboseLevel==1) cout << "Beam " << fNuclei1->GetName() << endl;
+ }
+
+ else if (DataBuffer=="Target=") {
+ check_Target = true ;
+ ReactionFile >> DataBuffer;
+ fNuclei2 = new Nucleus(DataBuffer);
+ if(fVerboseLevel==1) cout << "Target " << fNuclei2->GetName() << endl;
+ }
+
+ else if (DataBuffer=="Light=" || DataBuffer=="Nuclei3=") {
+ check_Light = true ;
+ ReactionFile >> DataBuffer;
+ fNuclei3 = new Nucleus(DataBuffer);
+ if(fVerboseLevel==1) cout << "Light " << fNuclei3->GetName() << endl;
+ }
+
+ else if (DataBuffer== "Heavy="|| DataBuffer=="Nuclei4=") {
+ check_Heavy = true ;
+ ReactionFile >> DataBuffer;
+ fNuclei4 = new Nucleus(DataBuffer);
+ if(fVerboseLevel==1) cout << "Heavy " << fNuclei4->GetName() << endl;
+ }
+
+ else if (DataBuffer=="ExcitationEnergy3=" || DataBuffer=="ExcitationEnergyLight=") {
+ check_ExcitationEnergy3 = true ;
+ ReactionFile >> DataBuffer;
+ fExcitation3 = atof(DataBuffer.c_str()) * MeV;
+ if(fVerboseLevel==1) cout << "Excitation Energy Nuclei 3: " << ExcitationEnergy3 / MeV << " MeV" << endl;
+ }
+
+ else if (DataBuffer=="ExcitationEnergy4=" || DataBuffer=="ExcitationEnergyHeavy=") {
+ check_ExcitationEnergy4 = true ;
+ ReactionFile >> DataBuffer;
+ fExcitation4 = atof(DataBuffer.c_str()) * MeV;
+ if(fVerboseLevel==1) cout << "Excitation Energy Nuclei 4: " << ExcitationEnergy4 / MeV << " MeV" << endl;
+ }
+
+ else if (DataBuffer== "CrossSectionPath=") {
+ check_CrossSectionPath = true ;
+ string FileName,HistName;
+ ReactionFile >> FileName >> HistName;
+ if(fVerboseLevel==1) cout << "Reading Cross Section file: " << FileName << endl;
+ fCrossSectionHist = Read1DProfile(FileName, HistName );
+ }
+
+ else if (DataBuffer.compare(0, 17, "HalfOpenAngleMin=") == 0) {
+ ReactionFile >> DataBuffer;
+ CSHalfOpenAngleMin = atof(DataBuffer.c_str()) * deg;
+ if(fVerboseLevel==1) cout << "HalfOpenAngleMin " << CSHalfOpenAngleMin / deg << " degree" << endl;
+ }
+
+ else if (DataBuffer.compare(0, 17, "HalfOpenAngleMax=") == 0) {
+ ReactionFile >> DataBuffer;
+ CSHalfOpenAngleMax = atof(DataBuffer.c_str()) * deg;
+ if(fVerboseLevel==1) cout << "HalfOpenAngleMax " << CSHalfOpenAngleMax / deg << " degree" << endl;
+ }
- cout << "Exc Nuclei 3 = " << fExcitation3 << " MeV" << endl;
- cout << "Exc Nuclei 4 = " << fExcitation4 << " MeV" << endl;
- cout << "Qgg = " << fQValue << " MeV" << endl;
- }
-
-
-
-//////////////////////////////////////////////////////////////////////////////////////////////////////////
-void Reaction::ReadConfigurationFile(string Path)
- {
- ////////General Reading needs////////
- string LineBuffer;
- string DataBuffer;
-
- ////////Reaction Setting needs///////
- string Beam, Target, Heavy, Light, CrossSectionPath ;
- double BeamEnergy = 0 , ExcitationEnergy3 = 0, ExcitationEnergy4 = 0;
- double CSHalfOpenAngleMin = 0, CSHalfOpenAngleMax = 0;
- bool ReadingStatus = false ;
- bool check_Beam = false ;
- bool check_Target = false ;
- bool check_Light = false ;
- bool check_Heavy = false ;
- bool check_ExcitationEnergy3 = false ;
- bool check_ExcitationEnergy4 = false ;
- bool check_BeamEnergy = false ;
- bool check_CrossSectionPath = false ;
-
-
- //////////////////////////////////////////////////////////////////////////////////////////
- ifstream ReactionFile;
- string GlobalPath = getenv("NPTOOL");
- string StandardPath = GlobalPath + "/Inputs/EventGenerator/" + Path;
- ReactionFile.open(StandardPath.c_str());
-
- if (ReactionFile.is_open()) {cout << "Reading Reaction File " << Path << endl ;}
-
- // In case the file is not found in the standard path, the programm try to interpret the file name as an absolute or relative file path.
- else
- {
- ReactionFile.open( Path.c_str() );
- if(ReactionFile.is_open()) { cout << "Reading Reaction File " << Path << endl;}
-
- else {cout << "Reaction File " << Path << " not found" << endl;return;}
- }
-
- while (!ReactionFile.eof()) {
- //Pick-up next line
- getline(ReactionFile, LineBuffer);
-
-
-
- if (LineBuffer.compare(0, 9, "Transfert") == 0) { ReadingStatus = true ;}
-
-
- while(ReadingStatus){
-
- ReactionFile >> DataBuffer;
-
- //Search for comment Symbol %
- if (LineBuffer.compare(0, 1, "%") == 0) {/* Do Nothing */;}
-
- else if (DataBuffer=="Beam=") {
- check_Beam = true ;
- ReactionFile >> DataBuffer;
- Beam = DataBuffer;
- cout << "Beam " << Beam << endl;
- }
-
- else if (DataBuffer=="Target=") {
- check_Target = true ;
- ReactionFile >> DataBuffer;
- Target = DataBuffer;
- cout << "Target " << Target << endl;
- }
-
- else if (DataBuffer=="Light=") {
- check_Light = true ;
- ReactionFile >> DataBuffer;
- Light = DataBuffer;
- cout << "Light " << Light << endl;
- }
-
- else if (DataBuffer== "Heavy=") {
- check_Heavy = true ;
- ReactionFile >> DataBuffer;
- Heavy = DataBuffer;
- cout << "Heavy " << Heavy << endl;
- }
-
- else if (DataBuffer=="ExcitationEnergy3=" || DataBuffer=="ExcitationEnergyLight=") {
- check_ExcitationEnergy3 = true ;
- ReactionFile >> DataBuffer;
- ExcitationEnergy3 = atof(DataBuffer.c_str()) * MeV;
- cout << "Excitation Energy Nuclei 3: " << ExcitationEnergy3 / MeV << " MeV" << endl;
- }
-
- else if (DataBuffer=="ExcitationEnergy4=" || DataBuffer=="ExcitationEnergyHeavy=") {
- check_ExcitationEnergy4 = true ;
- ReactionFile >> DataBuffer;
- ExcitationEnergy4 = atof(DataBuffer.c_str()) * MeV;
- cout << "Excitation Energy Nuclei 4: " << ExcitationEnergy4 / MeV << " MeV" << endl;
- }
-
- else if (DataBuffer=="BeamEnergy=") {
- check_BeamEnergy = true ;
- ReactionFile >> DataBuffer;
- BeamEnergy = atof(DataBuffer.c_str()) * MeV;
- cout << "Beam Energy " << BeamEnergy / MeV << " MeV" << endl;
- }
-
- else if (DataBuffer== "CrossSectionPath=") {
- check_CrossSectionPath = true ;
- ReactionFile >> CrossSectionPath;
- cout << "Cross Section File: " << CrossSectionPath << endl ;
- }
-
- else if (DataBuffer.compare(0, 17, "HalfOpenAngleMin=") == 0) {
- ReactionFile >> DataBuffer;
- CSHalfOpenAngleMin = atof(DataBuffer.c_str()) * deg;
- cout << "HalfOpenAngleMin " << CSHalfOpenAngleMin / deg << " degree" << endl;
- }
-
- else if (DataBuffer.compare(0, 17, "HalfOpenAngleMax=") == 0) {
- ReactionFile >> DataBuffer;
- CSHalfOpenAngleMax = atof(DataBuffer.c_str()) * deg;
- cout << "HalfOpenAngleMax " << CSHalfOpenAngleMax / deg << " degree" << endl;
- }
-
-
- ///////////////////////////////////////////////////
- // If no Transfert Token and no comment, toggle out
- else
- {/*Ignore Token used by G4*/}
-
- ///////////////////////////////////////////////////
- // If all Token found toggle out
- if (check_Beam && check_Target && check_Light && check_Heavy && check_ExcitationEnergy3
- && check_ExcitationEnergy4 && check_BeamEnergy && check_CrossSectionPath)
- ReadingStatus = false;
- }
- }
-
- SetEveryThing(Beam, Target, Light, Heavy,BeamEnergy,ExcitationEnergy3, ExcitationEnergy4,CrossSectionPath, CSHalfOpenAngleMin, CSHalfOpenAngleMax);
-
- ReactionFile.close();
- }
-
-
-////////////////////////////////////////////////////////////////////////////////////////////
-void Reaction::initializePrecomputeVariable()
- {
- m1 = fNuclei1->Mass();
- m2 = fNuclei2->Mass();
- m3 = fNuclei3->Mass() + fExcitation3;
- m4 = fNuclei4->Mass() + fExcitation4;
-
- s = m1*m1 + m2*m2 + 2*m2*(fBeamEnergy + m1);
-
- fTotalEnergyImpulsionCM = TLorentzVector(0,0,0,sqrt(s));
-
- ECM_1 = (s + m1*m1 - m2*m2)/(2*sqrt(s));
- ECM_2 = (s + m2*m2 - m1*m1)/(2*sqrt(s));
- ECM_3 = (s + m3*m3 - m4*m4)/(2*sqrt(s));
- ECM_4 = (s + m4*m4 - m3*m3)/(2*sqrt(s));
-
- pCM_3 = sqrt(ECM_3*ECM_3 - m3*m3);
- pCM_4 = sqrt(ECM_4*ECM_4 - m4*m4);
-
- fImpulsionLab_1 = TVector3(0,0,sqrt(fBeamEnergy*fBeamEnergy + 2*fBeamEnergy*m1));
- fImpulsionLab_2 = TVector3(0,0,0);
-
- fEnergyImpulsionLab_1 = TLorentzVector(fImpulsionLab_1,m1+fBeamEnergy);
- fEnergyImpulsionLab_2 = TLorentzVector(fImpulsionLab_2,m2);
-
- fTotalEnergyImpulsionLab = fEnergyImpulsionLab_1 + fEnergyImpulsionLab_2;
-
- BetaCM = fTotalEnergyImpulsionLab.Beta();
-
- fEnergyImpulsionCM_1 = fEnergyImpulsionLab_1;
- fEnergyImpulsionCM_1.Boost(0,0,-BetaCM);
-
- fEnergyImpulsionCM_2 = fEnergyImpulsionLab_2;
- fEnergyImpulsionCM_2.Boost(0,0,-BetaCM);
+ else if (DataBuffer== "Shoot3=" || DataBuffer== "ShootLight=") {
+ check_shoot3 = true ;
+ ReactionFile >> DataBuffer;
+
+ if(atoi(DataBuffer.c_str()) == 0 )
+ fshoot3 = true;
+
+ if(fVerboseLevel==1 && fshoot3) cout << "Shoot 3 : Yes " << endl;
+ else if (fVerboseLevel==1 ) cout << "Shoot 3 : No " << endl;
+ }
+
+ else if (DataBuffer== "Shoot4=" || DataBuffer== "ShootHeavy=") {
+ check_shoot4 = true ;
+ ReactionFile >> DataBuffer;
+
+ if(atoi(DataBuffer.c_str()) == 0 )
+ fshoot4 = true;
+
+ if(fVerboseLevel==1 && fshoot4) cout << "Shoot 4 : Yes " << endl;
+ else if (fVerboseLevel==1 ) cout << "Shoot 4 : No " << endl;
+ }
+
+ ///////////////////////////////////////////////////
+ // If no Transfert Token and no comment, toggle out
+ else
+ {/*Ignore Token used by G4*/}
+
+ ///////////////////////////////////////////////////
+ // If all Token found toggle out
+ if (check_Beam && check_Target && check_Light && check_Heavy && check_ExcitationEnergy3
+ && check_ExcitationEnergy4 && check_CrossSectionPath && check_shoot4 && check_shoot3)
+ ReadingStatus = false;
+ }
+ }
+
+ // Pick up the beam energy from the Beam event generator
+ NPL::Beam* localBeam= new NPL::Beam();
+ // localBeam->SetVerboseLevel(0);
+ localBeam->ReadConfigurationFile(Path);
+
+ // Modifiy the CS to shoot only within ]HalfOpenAngleMin,HalfOpenAngleMax[
+ SetCSAngle(CSHalfOpenAngleMin,CSHalfOpenAngleMax);
+ delete localBeam;
+ ReactionFile.close();
+ initializePrecomputeVariable();
}
-////////////////////////////////////////////////////////////////////////////////////////////
-void Reaction::SetNuclei3(double EnergyLab, double ThetaLab)
-{
+
+////////////////////////////////////////////////////////////////////////////////////////////
+void Reaction::initializePrecomputeVariable(){
+ m1 = fNuclei1->Mass();
+ m2 = fNuclei2->Mass();
+ m3 = fNuclei3->Mass() + fExcitation3;
+ m4 = fNuclei4->Mass() + fExcitation4;
+
+ s = m1*m1 + m2*m2 + 2*m2*(fBeamEnergy + m1);
+
+ fTotalEnergyImpulsionCM = TLorentzVector(0,0,0,sqrt(s));
+
+ ECM_1 = (s + m1*m1 - m2*m2)/(2*sqrt(s));
+ ECM_2 = (s + m2*m2 - m1*m1)/(2*sqrt(s));
+ ECM_3 = (s + m3*m3 - m4*m4)/(2*sqrt(s));
+ ECM_4 = (s + m4*m4 - m3*m3)/(2*sqrt(s));
+
+ pCM_3 = sqrt(ECM_3*ECM_3 - m3*m3);
+ pCM_4 = sqrt(ECM_4*ECM_4 - m4*m4);
+
+ fImpulsionLab_1 = TVector3(0,0,sqrt(fBeamEnergy*fBeamEnergy + 2*fBeamEnergy*m1));
+ fImpulsionLab_2 = TVector3(0,0,0);
+
+ fEnergyImpulsionLab_1 = TLorentzVector(fImpulsionLab_1,m1+fBeamEnergy);
+ fEnergyImpulsionLab_2 = TLorentzVector(fImpulsionLab_2,m2);
+
+ fTotalEnergyImpulsionLab = fEnergyImpulsionLab_1 + fEnergyImpulsionLab_2;
+
+ BetaCM = fTotalEnergyImpulsionLab.Beta();
+
+ fEnergyImpulsionCM_1 = fEnergyImpulsionLab_1;
+ fEnergyImpulsionCM_1.Boost(0,0,-BetaCM);
+
+ fEnergyImpulsionCM_2 = fEnergyImpulsionLab_2;
+ fEnergyImpulsionCM_2.Boost(0,0,-BetaCM);
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+void Reaction::SetNuclei3(double EnergyLab, double ThetaLab){
double p3 = sqrt(pow(EnergyLab,2) + 2*m3*EnergyLab);
fEnergyImpulsionLab_3 = TLorentzVector(p3*sin(ThetaLab),0,p3*cos(ThetaLab),EnergyLab+m3);
@@ -486,23 +414,22 @@ void Reaction::SetNuclei3(double EnergyLab, double ThetaLab)
-////////////////////////////////////////////////////////////////////////////////////////////
-TGraph* Reaction::GetKinematicLine3()
-{
- int size = 360;
+////////////////////////////////////////////////////////////////////////////////////////////
+TGraph* Reaction::GetKinematicLine3(){
+ int size = 360;
double x[size];
double y[size];
double theta3,E3,theta4,E4;
- for (int i = 0; i < size; ++i)
- {
- SetThetaCM(((double)i)/2*deg);
+ for (int i = 0; i < size; ++i){
+ SetThetaCM(((double)i)/2*deg);
KineRelativistic(theta3, E3, theta4, E4);
fNuclei3->SetKineticEnergy(E3);
- x[i] = theta3/deg;
- y[i] = E3;
- }
+ x[i] = theta3/deg;
+ y[i] = E3;
+ }
+
TGraph* KineLine3 = new TGraph(size,x,y);
KineLine3->SetTitle("Kinematic Line of particule 3");
return(KineLine3);
@@ -510,85 +437,80 @@ TGraph* Reaction::GetKinematicLine3()
-////////////////////////////////////////////////////////////////////////////////////////////
-TGraph* Reaction::GetKinematicLine4()
-{
- int size = 360;
+////////////////////////////////////////////////////////////////////////////////////////////
+TGraph* Reaction::GetKinematicLine4(){
+ int size = 360;
double x[size];
double y[size];
double theta3,E3,theta4,E4;
-
+
for (int i = 0; i < size; ++i)
- {
- SetThetaCM(((double)i)/2*deg);
+ {
+ SetThetaCM(((double)i)/2*deg);
KineRelativistic(theta3, E3, theta4, E4);
fNuclei4->SetKineticEnergy(E4);
- x[i] = theta4/deg;
- y[i] = E4;
- }
+ x[i] = theta4/deg;
+ y[i] = E4;
+ }
TGraph* KineLine4= new TGraph(size,x,y);
KineLine4->SetTitle("Kinematic Line of particule 4");
return(KineLine4);
}
-////////////////////////////////////////////////////////////////////////////////////////////
-TGraph* Reaction::GetBrhoLine3()
-{
- int size = 360;
+////////////////////////////////////////////////////////////////////////////////////////////
+TGraph* Reaction::GetBrhoLine3(){
+ int size = 360;
double x[size];
double y[size];
double theta3,E3,theta4,E4;
double Brho;
for (int i = 0; i < size; ++i)
- {
- SetThetaCM(((double)i)/2*deg);
+ {
+ SetThetaCM(((double)i)/2*deg);
KineRelativistic(theta3, E3, theta4, E4);
fNuclei3->SetKineticEnergy(E3);
Brho = fNuclei3->GetBrho();
- x[i] = theta3/deg;
- y[i] = Brho;
- }
+ x[i] = theta3/deg;
+ y[i] = Brho;
+ }
TGraph* LineBrho3= new TGraph(size,x,y);
return(LineBrho3);
}
-////////////////////////////////////////////////////////////////////////////////////////////
-TGraph* Reaction::GetThetaLabVersusThetaCM()
-{
- int size = 360;
+////////////////////////////////////////////////////////////////////////////////////////////
+TGraph* Reaction::GetThetaLabVersusThetaCM(){
+ int size = 360;
double x[size];
double y[size];
double theta3,E3,theta4,E4;
- for (int i = 0; i < size; ++i)
- {
- SetThetaCM(((double)i)/2*deg);
+ for (int i = 0; i < size; ++i){
+ SetThetaCM(((double)i)/2*deg);
KineRelativistic(theta3, E3, theta4, E4);
- x[i] = fThetaCM/deg;
- y[i] = theta3/deg;
- }
+ x[i] = fThetaCM/deg;
+ y[i] = theta3/deg;
+ }
+
TGraph* AngleLine= new TGraph(size,x,y);
return(AngleLine);
}
-////////////////////////////////////////////////////////////////////////////////////////////
-void Reaction::PrintKinematic()
-{
- int size = 360;
+////////////////////////////////////////////////////////////////////////////////////////////
+void Reaction::PrintKinematic(){
+ int size = 360;
double theta3,E3,theta4,E4,Brho3,Brho4;
cout << endl;
cout << "*********************** Print Kinematic ***********************" << endl;
cout << "ThetaCM" << " " << "ThetaLab" << " " << "EnergyLab3" << " " << "Brho3" << " " << "EnergyLab4" << " " << "Brho4" << endl;
- for (int i = 0; i < size; ++i)
- {
- SetThetaCM(((double)i)/2*deg);
+ for (int i = 0; i < size; ++i){
+ SetThetaCM(((double)i)/2*deg);
KineRelativistic(theta3, E3, theta4, E4);
fNuclei3->SetKineticEnergy(E3);
@@ -598,11 +520,18 @@ void Reaction::PrintKinematic()
Brho4 = fNuclei4->GetBrho();
cout << (double)i/2 << " " << theta3/deg << " " << E3 << " " << Brho3 << " " << E4 << " " << Brho4 << endl;
- }
+ }
}
-
+////////////////////////////////////////////////////////////////////////////////////////////
+void Reaction::SetCSAngle(double CSHalfOpenAngleMin,double CSHalfOpenAngleMax){
+
+ for (int i = 0 ; i< fCrossSectionHist->GetNbinsX(); i++)
+ if( fCrossSectionHist->GetBinCenter(i) > CSHalfOpenAngleMax && fCrossSectionHist->GetBinCenter(i) < CSHalfOpenAngleMin)
+ fCrossSectionHist->SetBinContent(i,0);
+
+}
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index 576ef726bdfc905c73312aac6ed65482b616550d..7e0e2537f2602044ccfce65ec40c1cd7d6ad7e6a 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -34,12 +34,12 @@
* *
* *
*****************************************************************************/
-// C++ header
+ // C++ header
#include <string>
#include "NPNucleus.h"
-// ROOT header
+ // ROOT header
#include "TLorentzVector.h"
#include "TLorentzRotation.h"
#include "TVector3.h"
@@ -54,13 +54,18 @@ namespace NPL{
public: // Constructors and Destructors
Reaction();
- Reaction(string name1, string name2, string name3, string name4, double BeamEnergy , double ExcitationEnergyLight, double ExcitationEnergyHeavy, string Path, double CSThetaMin = 0, double CSThetaMax = 180);
~Reaction();
public: // Various Method
- void SetEveryThing(string name1, string name2, string name3, string name4, double BeamEnergy, double ExcitationEnergyLight, double ExcitationEnergyHeavy, string Path, double CSThetaMin, double CSThetaMax);
void ReadConfigurationFile(string Path);
+ private:
+ int fVerboseLevel;
+
+ private: // use for Monte Carlo simulation
+ bool fshoot3;
+ bool fshoot4;
+
private:
Nucleus *fNuclei1; // Beam
Nucleus *fNuclei2; // Target
@@ -72,17 +77,18 @@ namespace NPL{
double fExcitation3; // Excitation energy in MeV
double fExcitation4; // Excitation energy in MeV
TH1F* fCrossSectionHist;
- double* fCrossSection; // Differential CrossSection
- int fCrossSectionSize; // Size of array containing Differention CrossSection
- double fCrossSectionAngleMin; // Minimum angle of the differential cross-section given by the user
- double fCrossSectionAngleMax; // Maximum angle of the differential cross-section given by the user
-
+
public:
- // Getters and Setters
+ // Getters and Setters
void SetBeamEnergy(double eBeam) {fBeamEnergy = eBeam; initializePrecomputeVariable();}
void SetThetaCM(double angle) {fThetaCM = angle; initializePrecomputeVariable();}
void SetExcitation3(double exci) {fExcitation3 = exci; initializePrecomputeVariable();}
void SetExcitation4(double exci) {fExcitation4 = exci; initializePrecomputeVariable();}
+ // For retro compatibility
+ void SetExcitationLight(double exci) {fExcitation3 = exci; initializePrecomputeVariable();}
+ void SetExcitationHeavy(double exci) {fExcitation4 = exci; initializePrecomputeVariable();}
+ void SetVerboseLevel(int verbose) {fVerboseLevel = verbose;}
+
double GetBeamEnergy() const {return fBeamEnergy;}
double GetThetaCM() const {return fThetaCM;}
double GetExcitation3() const {return fExcitation3;}
@@ -93,11 +99,13 @@ namespace NPL{
Nucleus* GetNucleus3() const {return fNuclei3;}
Nucleus* GetNucleus4() const {return fNuclei4;}
TH1F* GetCrossSectionHist() const {return fCrossSectionHist;}
- double* GetCrossSection() const {return fCrossSection;}
- int GetCrossSectionSize() const {return fCrossSectionSize;}
- double GetCrossSectionAngleMin() const {return fCrossSectionAngleMin;}
- double GetCrossSectionAngleMax() const {return fCrossSectionAngleMax;}
-
+ int GetVerboseLevel() const {return fVerboseLevel;}
+ bool GetShoot3() const {return fshoot3;}
+ bool GetShoot4() const {return fshoot4;}
+
+ public:
+ // Modify the CS histo to so cross section shoot is within ]HalfOpenAngleMin,HalfOpenAngleMax[
+ void SetCSAngle(double CSHalfOpenAngleMin,double CSHalfOpenAngleMax);
private: // intern precompute variable
void initializePrecomputeVariable();
@@ -106,7 +114,7 @@ namespace NPL{
double m3;
double m4;
- // Lorents Vector
+ // Lorents Vector
TLorentzVector fEnergyImpulsionLab_1;
TLorentzVector fEnergyImpulsionLab_2;
TLorentzVector fEnergyImpulsionLab_3;
@@ -119,7 +127,7 @@ namespace NPL{
TLorentzVector fEnergyImpulsionCM_4;
TLorentzVector fTotalEnergyImpulsionCM;
- // Impulsion Vector3
+ // Impulsion Vector3
TVector3 fImpulsionLab_1;
TVector3 fImpulsionLab_2;
TVector3 fImpulsionLab_3;
@@ -130,7 +138,7 @@ namespace NPL{
TVector3 fImpulsionCM_3;
TVector3 fImpulsionCM_4;
- // CM Energy composante & CM impulsion norme
+ // CM Energy composante & CM impulsion norme
Double_t ECM_1;
Double_t ECM_2;
Double_t ECM_3;
@@ -138,43 +146,41 @@ namespace NPL{
Double_t pCM_3;
Double_t pCM_4;
- // Mandelstam variable
+ // Mandelstam variable
Double_t s;
- // Center of Mass Kinematic
+ // Center of Mass Kinematic
Double_t BetaCM;
public: // Kinematics
- // Check that the reaction is alowed
+ // Check that the reaction is alowed
bool CheckKinematic();
- // Use fCrossSectionHist to shoot a Random ThetaCM and set fThetaCM to this value
+ // Use fCrossSectionHist to shoot a Random ThetaCM and set fThetaCM to this value
double ShootRandomThetaCM();
- // Compute ThetaLab and EnergyLab for product of reaction
+ // Compute ThetaLab and EnergyLab for product of reaction
void KineRelativistic(double &ThetaLab3, double &KineticEnergyLab3,
double &ThetaLab4, double &KineticEnergyLab4);
- // Return Excitation Energy
+ // Return Excitation Energy
double ReconstructRelativistic(double EnergyLab, double ThetaLab);
-
- // Return ThetaCM
- // EnergyLab: energy measured in the laboratory frame
- // ExcitationEnergy: excitation energy previously calculated.
+ // Return ThetaCM
+ // EnergyLab: energy measured in the laboratory frame
+ // ExcitationEnergy: excitation energy previously calculated.
double EnergyLabToThetaCM(double EnergyLab, double ThetaLab);
void SetNuclei3(double EnergyLab, double ThetaLab);
-
TGraph* GetKinematicLine3();
TGraph* GetKinematicLine4();
TGraph* GetBrhoLine3();
TGraph* GetThetaLabVersusThetaCM();
void PrintKinematic();
- // Print private paremeter
+ // Print private paremeter
void Print() const;
};
}
diff --git a/NPLib/Tools/NPOptionManager.cxx b/NPLib/Tools/NPOptionManager.cxx
index 37d595360dd2b2d485ebd9be519933571f143f7c..b359dd5f969768c657e239f4825d9717fa90ba67 100644
--- a/NPLib/Tools/NPOptionManager.cxx
+++ b/NPLib/Tools/NPOptionManager.cxx
@@ -51,6 +51,7 @@ NPOptionManager::NPOptionManager(int argc, char** argv)
fOutputFileName = fDefaultOutputFileName;
fRunToReadFileName = fDefaultRunToReadFileName;
fCalibrationFileName = fDefaultCalibrationFileName;
+ fVerboseLevel = 1;
fDisableAllBranchOption = false;
fInputPhysicalTreeOption = false;
@@ -61,29 +62,33 @@ NPOptionManager::NPOptionManager(int argc, char** argv)
else if (argument == "--event-generator" && argc >= i + 1) fReactionFileName = argv[i+1] ;
- else if (argument == "-E" && argc >= i + 1) fReactionFileName = argv[i+1] ;
+ else if (argument == "-E" && argc >= i + 1) fReactionFileName = argv[i+1] ;
- else if (argument == "--detector" && argc >= i + 1) fDetectorFileName = argv[i+1] ;
+ else if (argument == "--detector" && argc >= i + 1) fDetectorFileName = argv[i+1] ;
- else if (argument == "-D" && argc >= i + 1) fDetectorFileName = argv[i+1] ;
+ else if (argument == "-D" && argc >= i + 1) fDetectorFileName = argv[i+1] ;
+
+ else if (argument == "--output" && argc >= i + 1) fOutputFileName = argv[i+1] ;
- else if (argument == "--output" && argc >= i + 1) fOutputFileName = argv[i+1] ;
+ else if (argument == "-O" && argc >= i + 1) fOutputFileName = argv[i+1] ;
- else if (argument == "-O" && argc >= i + 1) fOutputFileName = argv[i+1] ;
+ else if (argument == "--run" && argc >= i + 1) fRunToReadFileName = argv[i+1] ;
- else if (argument == "--run" && argc >= i + 1) fRunToReadFileName = argv[i+1] ;
+ else if (argument == "-R" && argc >= i + 1) fRunToReadFileName = argv[i+1] ;
- else if (argument == "-R" && argc >= i + 1) fRunToReadFileName = argv[i+1] ;
+ else if (argument == "--cal" && argc >= i + 1) fCalibrationFileName = argv[i+1] ;
- else if (argument == "--cal" && argc >= i + 1) fCalibrationFileName = argv[i+1] ;
+ else if (argument == "-C" && argc >= i + 1) fCalibrationFileName = argv[i+1] ;
- else if (argument == "-C" && argc >= i + 1) fCalibrationFileName = argv[i+1] ;
+ else if (argument == "-V" && argc >= i + 1) fVerboseLevel = atoi(argv[i+1]) ;
- else if (argument == "--disable-branch") fDisableAllBranchOption = true ;
+ else if (argument == "--verbose" && argc >= i + 1) fVerboseLevel = atoi(argv[i+1]) ;
- else if (argument == "--input-physical") fInputPhysicalTreeOption = true ;
+ else if (argument == "--disable-branch") fDisableAllBranchOption = true ;
- else if (argument == "-IP") fInputPhysicalTreeOption = true ;
+ else if (argument == "--input-physical") fInputPhysicalTreeOption = true ;
+
+ else if (argument == "-IP") fInputPhysicalTreeOption = true ;
//else ;
}
@@ -231,6 +236,7 @@ void NPOptionManager::DisplayHelp()
cout << "\t --detector　-D <arg>\t \t \t \t \t \t　Set arg as the detector configuration file" << endl ;
cout << "\t --event-generator　-E <arg>\t \t \t \t \t　Set arg as the event generator file" << endl ;
cout << "\t --output　-O <arg>\t \t \t \t \t \t　Set arg as the Output File Name (output tree)" << endl ;
+ cout << "\t --verbose -V <arg>\t \t \t \t \t \t \t　Set the verbose level of some of the object, 0 for nothing, 1 for normal printout. Error and warning are not affected" << endl ;
cout << endl << "NPAnalysis only:"<<endl;
cout << "\t --run -R <arg>\t \t \t \t \t \t \t　Set arg as the run to read file list" << endl ;
cout << "\t --cal -C <arg>\t \t \t \t \t \t \t　Set arg as the calibration file list" << endl ;
diff --git a/NPLib/Tools/NPOptionManager.h b/NPLib/Tools/NPOptionManager.h
index 410ffd87d85dcb89ea9efbcbf8aecb201f38e4dc..440d8530e12e54ae272d4517bd17657d567799eb 100644
--- a/NPLib/Tools/NPOptionManager.h
+++ b/NPLib/Tools/NPOptionManager.h
@@ -84,12 +84,14 @@ class NPOptionManager
string GetOutputFile() {return fOutputFileName;}
bool GetDisableAllBranchOption() {return fDisableAllBranchOption;}
bool GetInputPhysicalTreeOption() {return fInputPhysicalTreeOption;}
+ int GetVerboseLevel() {return fVerboseLevel;}
// Setters
- void SetReactionFile(string name) {fReactionFileName = name;}
- void SetDetectorFile(string name) {fDetectorFileName = name;}
- void SetRunToReadFile(string name) {fRunToReadFileName = name;}
-
+ void SetReactionFile(string name) {fReactionFileName = name;}
+ void SetDetectorFile(string name) {fDetectorFileName = name;}
+ void SetRunToReadFile(string name) {fRunToReadFileName = name;}
+ void GetVerboseLevel(int VerboseLevel) {fVerboseLevel = VerboseLevel;}
+
private:
// default values
string fDefaultReactionFileName;
@@ -105,6 +107,7 @@ class NPOptionManager
string fOutputFileName;
bool fDisableAllBranchOption;
bool fInputPhysicalTreeOption;
+ int fVerboseLevel; // 0 for not talk, 1 for talking
};
#endif
diff --git a/NPSimulation/include/EventGeneratorTwoBodyReaction.hh b/NPSimulation/include/EventGeneratorTwoBodyReaction.hh
index b7895e025f719576193371e71e5b009394a13160..b4143e2463abf0cbb9851a149e5148fcc6ebee14 100644
--- a/NPSimulation/include/EventGeneratorTwoBodyReaction.hh
+++ b/NPSimulation/include/EventGeneratorTwoBodyReaction.hh
@@ -11,7 +11,7 @@
* Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
* *
* Creation Date : January 2009 *
- * Last update : January 2011 *
+ * Last update : January 2013 *
*---------------------------------------------------------------------------*
* Decription: *
* This event Generator is used to simulated two body two body reaction. *
@@ -49,19 +49,6 @@ class EventGeneratorTwoBodyReaction : public VEventGenerator
// Default constructor used to allocate memory
EventGeneratorTwoBodyReaction();
- // This is the constructor to be used
- EventGeneratorTwoBodyReaction(string name1 , // Beam nuclei
- string name2 , // Target nuclei
- string name3 , // Product of reaction
- string name4 , // Product of reaction
- double ExcitationEnergyLight , // Excitation of Light Nuclei
- double ExcitationEnergyHeavy , // Excitation of Heavy Nuclei
- bool ShootLight ,
- bool ShootHeavy ,
- string Path,
- double CSThetaMin,
- double CSThetaMax); // Path of the differentiel Cross Section
-
// Default Destructor
virtual ~EventGeneratorTwoBodyReaction();
@@ -84,8 +71,6 @@ class EventGeneratorTwoBodyReaction : public VEventGenerator
private: // Reaction and CrossSection Shoot
Reaction* m_Reaction;
- G4double m_HalfOpenAngleMin; // Min Half open angle of the source
- G4double m_HalfOpenAngleMax; // Max Half open angle of the source
private: // Beam Parameters
@@ -94,18 +79,6 @@ class EventGeneratorTwoBodyReaction : public VEventGenerator
// Other methods
void Print() const;
void InitializeRootOutput();
-
- void SetEverything(string name1 , // Beam nuclei
- string name2 , // Target nuclei
- string name3 , // Product of reaction
- string name4 , // Product of reaction
- double ExcitationEnergyLight , // Excitation of Light Nuclei
- double ExcitationEnergyHeavy , // Excitation of Heavy Nuclei
- bool ShootLight ,
- bool ShootHeavy ,
- string Path,
- double CSThetaMin,
- double CSThetaMax); // Path of the differentiel Cross Section
};
#endif
diff --git a/NPSimulation/src/EventGeneratorTwoBodyReaction.cc b/NPSimulation/src/EventGeneratorTwoBodyReaction.cc
index 4d11ea980dcfecc3516769e10fc64a65a751673d..c71dd490707ebcd5aeb18a0ec1d7c36ac40394eb 100644
--- a/NPSimulation/src/EventGeneratorTwoBodyReaction.cc
+++ b/NPSimulation/src/EventGeneratorTwoBodyReaction.cc
@@ -9,7 +9,7 @@
* Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
* *
* Creation Date : January 2009 *
- * Last update : January 2011 *
+ * Last update : January 2013 *
*---------------------------------------------------------------------------*
* Decription: *
* This event Generator is used to simulated two body two body reaction. *
@@ -23,7 +23,6 @@
* + 23/01/2013: Class change name (ild name EventGeneratorTransfert) *
* (A. MATTA) *
* *
- * *
*****************************************************************************/
// C++ headers
#include <iostream>
@@ -52,9 +51,7 @@ EventGeneratorTwoBodyReaction::EventGeneratorTwoBodyReaction()
: m_ShootLight(0),
m_ShootHeavy(0),
m_Target(0),
-m_Reaction(new Reaction),
-m_HalfOpenAngleMin(0),
-m_HalfOpenAngleMax(180)
+m_Reaction(new Reaction)
{
//------------- Default Constructor -------------
m_ParticleStack= ParticleStack::getInstance();
@@ -75,32 +72,6 @@ EventGeneratorTwoBodyReaction::~EventGeneratorTwoBodyReaction(){
delete m_Reaction;
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-EventGeneratorTwoBodyReaction::EventGeneratorTwoBodyReaction( string name1 , // Beam nuclei
- string name2 , // Target nuclei
- string name3 , // Product of reaction
- string name4 , // Product of reaction
- double ExcitationEnergyLight , // Excitation of Light Nuclei
- double ExcitationEnergyHeavy , // Excitation of Heavy Nuclei
- bool ShootLight ,
- bool ShootHeavy ,
- string Path ,
- double CSThetaMin ,
- double CSThetaMax) // Path of the differentiel Cross Section
-{
- SetEverything( name1,
- name2,
- name3,
- name4,
- ExcitationEnergyLight,
- ExcitationEnergyHeavy,
- ShootLight,
- ShootHeavy,
- Path,
- CSThetaMin,
- CSThetaMax);
-
-}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorTwoBodyReaction::InitializeRootOutput(){
@@ -121,162 +92,11 @@ void EventGeneratorTwoBodyReaction::Print() const{
void EventGeneratorTwoBodyReaction::ReadConfiguration(string Path, int dump){
dump = 0 ;
- ////////General Reading needs////////
- string LineBuffer;
- string DataBuffer;
-
- ////////Reaction Setting needs///////
- string Beam, Target, Heavy, Light, CrossSectionPath ;
- G4double ExcitationEnergyLight = 0, ExcitationEnergyHeavy = 0;
- G4double CSHalfOpenAngleMin = 0, CSHalfOpenAngleMax = 180;
-
- bool ShootLight = false ;
- bool ShootHeavy = false ;
-
- bool ReadingStatus = false ;
- bool check_Beam = false ;
- bool check_Target = false ;
- bool check_Light = false ;
- bool check_Heavy = false ;
- bool check_ExcitationEnergyLight = false ;
- bool check_ExcitationEnergyHeavy = false ;
- bool check_CrossSectionPath = false ;
- bool check_ShootLight = false ;
- bool check_ShootHeavy = false ;
-
- //////////////////////////////////////////////////////////////////////////////////////////
- ifstream ReactionFile;
- ReactionFile.open(Path.c_str());
-
- if (ReactionFile.is_open()) {} else {
- return;
- }
+ m_Reaction->ReadConfigurationFile(Path);
- while (!ReactionFile.eof()) {
- //Pick-up next line
- getline(ReactionFile, LineBuffer);
- if (LineBuffer.compare(0, 15, "TwoBodyReaction") == 0) {
- ReadingStatus = true ;
-
- G4cout << "////////// Two Body Reaction found ///////////" << G4endl;
- }
-
-
- while(ReadingStatus){
-
- ReactionFile >> DataBuffer;
-
- //Search for comment Symbol %
- if (DataBuffer.compare(0, 1, "%") == 0) { ReactionFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
-
- else if (DataBuffer=="Beam=") {
- check_Beam = true ;
- ReactionFile >> DataBuffer;
- Beam = DataBuffer;
- m_BeamName = m_ParticleStack->ChangeNameToG4Standard(m_BeamName);
- G4cout << "Beam " << Beam << G4endl;
- }
-
- else if (DataBuffer=="Target=") {
- check_Target = true ;
- ReactionFile >> DataBuffer;
- Target = DataBuffer;
- G4cout << "Target " << Target << G4endl;
- }
-
- else if (DataBuffer=="Light=") {
- check_Light = true ;
- ReactionFile >> DataBuffer;
- Light = DataBuffer;
- G4cout << "Light " << Light << G4endl;
- }
-
- else if (DataBuffer=="Heavy=") {
- check_Heavy = true ;
- ReactionFile >> DataBuffer;
- Heavy = DataBuffer;
- G4cout << "Heavy " << Heavy << G4endl;
- }
-
- else if (DataBuffer=="ExcitationEnergy3=" || DataBuffer=="ExcitationEnergyLight=") {
- check_ExcitationEnergyLight = true ;
- ReactionFile >> DataBuffer;
- ExcitationEnergyLight = atof(DataBuffer.c_str()) * MeV;
- G4cout << "Excitation Energy Nuclei 3: " << ExcitationEnergyLight / MeV << " MeV" << G4endl;
- }
-
- else if (DataBuffer=="ExcitationEnergy4=" || DataBuffer=="ExcitationEnergyHeavy=") {
- check_ExcitationEnergyHeavy = true ;
- ReactionFile >> DataBuffer;
- ExcitationEnergyHeavy = atof(DataBuffer.c_str()) * MeV;
- G4cout << "Excitation Energy Nuclei 4: " << ExcitationEnergyHeavy / MeV << " MeV" << G4endl;
- }
+ m_ShootLight = m_Reaction->GetShoot3();
+ m_ShootHeavy = m_Reaction->GetShoot4();;
- else if (DataBuffer=="CrossSectionPath=") {
- check_CrossSectionPath = true ;
- ReactionFile >> CrossSectionPath;
- G4cout << "Cross Section File: " << CrossSectionPath << G4endl ;
- }
-
- else if (DataBuffer=="HalfOpenAngleMin=") {
- ReactionFile >> DataBuffer;
- CSHalfOpenAngleMin = atof(DataBuffer.c_str()) * deg;
- G4cout << "HalfOpenAngleMin " << CSHalfOpenAngleMin / deg << " degree" << G4endl;
- }
-
- else if (DataBuffer=="HalfOpenAngleMax=") {
- ReactionFile >> DataBuffer;
- CSHalfOpenAngleMax = atof(DataBuffer.c_str()) * deg;
- G4cout << "HalfOpenAngleMax " << CSHalfOpenAngleMax / deg << " degree" << G4endl;
- }
-
- else if (DataBuffer=="ShootLight=") {
- check_ShootLight = true ;
- ReactionFile >> DataBuffer;
- if (atof(DataBuffer.c_str()) == 1) ShootLight = true ;
- if (ShootLight) G4cout << "Shoot Light particle : yes" << G4endl;
- else G4cout << "Shoot Light particle : no" << G4endl;
- }
-
- else if (DataBuffer=="ShootHeavy=") {
- check_ShootHeavy = true ;
- ReactionFile >> DataBuffer;
- if (atof(DataBuffer.c_str()) == 1) ShootHeavy = true ;
- if (ShootHeavy) G4cout << "Shoot Heavy particle : yes" << G4endl;
- else G4cout << "Shoot Heavy particle : no" << G4endl;
- }
-
-
- ///////////////////////////////////////////////////
- // If no Transfert Token and no comment, toggle out
- else
- {ReadingStatus = false; G4cout << "WARNING : Wrong Token Sequence: Getting out " << G4endl ;}
-
- ///////////////////////////////////////////////////
- // If all Token found toggle out
- if(check_Beam && check_Target && check_Light && check_Heavy && check_ExcitationEnergyLight && check_ExcitationEnergyHeavy
- && check_CrossSectionPath && check_ShootLight && check_ShootHeavy){
- ReadingStatus = false ;
- G4cout << "////////////////////////////////" << G4endl;
- }
- }
-
-
- }
-
- SetEverything( Beam,
- Target,
- Light,
- Heavy,
- ExcitationEnergyLight,
- ExcitationEnergyHeavy,
- ShootLight,
- ShootHeavy,
- CrossSectionPath,
- CSHalfOpenAngleMin/deg,
- CSHalfOpenAngleMax/deg);
-
- ReactionFile.close();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -289,7 +109,7 @@ void EventGeneratorTwoBodyReaction::GenerateEvent(G4Event* anEvent){
G4int LightAx = m_Reaction->GetNucleus3()->GetA();
if (m_Target != 0) {
- //m_Target->WriteDEDXTable(G4ParticleTable::GetParticleTable()->GetIon(LightZx,LightAx, 0.) ,0, 1*GeV);
+ m_Target->WriteDEDXTable(G4ParticleTable::GetParticleTable()->GetIon(LightZx,LightAx, 0.) ,0, 2*m_Reaction->GetBeamEnergy());
}
}
@@ -407,25 +227,3 @@ void EventGeneratorTwoBodyReaction::GenerateEvent(G4Event* anEvent){
//Add the particle to the particle stack
m_ParticleStack->AddParticleToStack(HeavyParticle);
}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorTwoBodyReaction::SetEverything(string name1, // Beam nuclei
- string name2, // Target nuclei
- string name3, // Product of reaction
- string name4, // Product of reaction
- double ExcitationEnergyLight,// Excitation of Light Nuclei
- double ExcitationEnergyHeavy,// Excitation of Heavy Nuclei
- bool ShootLight,
- bool ShootHeavy,
- string Path,
- double CSThetaMin,
- double CSThetaMax)
-{
- m_Reaction = new Reaction(name1, name2, name3, name4, 0,ExcitationEnergyLight, ExcitationEnergyHeavy, Path, CSThetaMin, CSThetaMax);
- m_ShootLight = ShootLight;
- m_ShootHeavy = ShootHeavy;
- m_HalfOpenAngleMin = CSThetaMin;
- m_HalfOpenAngleMax = CSThetaMax;
-}
-