diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction
index 95cc7b9f6dd4c32d16aa7c6b74c9cec49afa7306..2a8c6d3a37d9a5888036750491f39f73b8980713 100644
--- a/Inputs/EventGenerator/10He.reaction
+++ b/Inputs/EventGenerator/10He.reaction
@@ -20,14 +20,14 @@ Transfert
ShootHeavy= 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-ParticleDecay 10He
+%ParticleDecay 10He
Daughter= 9He n
ExcitationEnergy= 0.5 0
DifferentialCrossSection= flat.txt
shoot= 1 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-GammaDecay 9He
+%GammaDecay 9He
Cascade
BranchingRatio= 30
Energies= 0.1
@@ -37,7 +37,7 @@ GammaDecay 9He
Energies= 0.2 0.3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-ParticleDecay 9He
+%ParticleDecay 9He
Daughter= 8He n
DifferentialCrossSection= flat.txt
shoot= 1 1
\ No newline at end of file
diff --git a/NPLib/Chio/TChio_anPhysics.h b/NPLib/Chio/TChio_anPhysics.h
index 7497d5c68a1bae28fb1d90967745d601a6535eb8..504888cc9136f4a1296ea56539f46e8c734a4e20 100644
--- a/NPLib/Chio/TChio_anPhysics.h
+++ b/NPLib/Chio/TChio_anPhysics.h
@@ -72,7 +72,7 @@ class TChio_anPhysics : public TObject, public NPA::VDetector
vector<double> rawAmplitude;
vector<double> calAmplitude;
- TF1 *fdecay;
+ TF1 *fdecay; //!
private: // Root Input and Output tree classes
diff --git a/NPLib/ComptonTelescope/TComptonTelescopeDataDict.cxx b/NPLib/ComptonTelescope/TComptonTelescopeDataDict.cxx
index 0b23325edd01454050ef86f2f6551126d4bc70d1..ce11ff788a132cca4164e4099d6373502099a09a 100644
--- a/NPLib/ComptonTelescope/TComptonTelescopeDataDict.cxx
+++ b/NPLib/ComptonTelescope/TComptonTelescopeDataDict.cxx
@@ -1,5 +1,5 @@
//
-// File generated by rootcint at Wed Jul 18 14:41:32 2012
+// File generated by rootcint at Thu Nov 1 10:31:35 2012
// Do NOT change. Changes will be lost next time file is generated
//
@@ -1643,12 +1643,16 @@ static void G__cpp_setup_global0() {
}
static void G__cpp_setup_global1() {
+}
+
+static void G__cpp_setup_global2() {
G__resetglobalenv();
}
extern "C" void G__cpp_setup_globalTComptonTelescopeDataDict() {
G__cpp_setup_global0();
G__cpp_setup_global1();
+ G__cpp_setup_global2();
}
/*********************************************************
diff --git a/NPLib/ComptonTelescope/TComptonTelescopeProcessDataDict.cxx b/NPLib/ComptonTelescope/TComptonTelescopeProcessDataDict.cxx
index 7e2d1a123a3830d98bb7e0042e2ca36cc0d450a5..51d3d8d45305e8f9ea57b414f4fc996e689ab8c8 100644
--- a/NPLib/ComptonTelescope/TComptonTelescopeProcessDataDict.cxx
+++ b/NPLib/ComptonTelescope/TComptonTelescopeProcessDataDict.cxx
@@ -1,5 +1,5 @@
//
-// File generated by rootcint at Wed Jul 18 14:41:34 2012
+// File generated by rootcint at Thu Nov 1 10:31:36 2012
// Do NOT change. Changes will be lost next time file is generated
//
@@ -747,12 +747,16 @@ static void G__cpp_setup_global0() {
}
static void G__cpp_setup_global1() {
+}
+
+static void G__cpp_setup_global2() {
G__resetglobalenv();
}
extern "C" void G__cpp_setup_globalTComptonTelescopeProcessDataDict() {
G__cpp_setup_global0();
G__cpp_setup_global1();
+ G__cpp_setup_global2();
}
/*********************************************************
diff --git a/NPLib/Eurogam/TEurogamDataDict.cxx b/NPLib/Eurogam/TEurogamDataDict.cxx
index 670ac9d7564023d2a0e3e590482a924c57d7986a..61d9df7e77ee90e36faaca510504b37481a56396 100644
--- a/NPLib/Eurogam/TEurogamDataDict.cxx
+++ b/NPLib/Eurogam/TEurogamDataDict.cxx
@@ -1,5 +1,5 @@
//
-// File generated by rootcint at Wed Jul 18 14:41:36 2012
+// File generated by rootcint at Thu Nov 1 10:31:39 2012
// Do NOT change. Changes will be lost next time file is generated
//
@@ -753,12 +753,16 @@ static void G__cpp_setup_global0() {
}
static void G__cpp_setup_global1() {
+}
+
+static void G__cpp_setup_global2() {
G__resetglobalenv();
}
extern "C" void G__cpp_setup_globalTEurogamDataDict() {
G__cpp_setup_global0();
G__cpp_setup_global1();
+ G__cpp_setup_global2();
}
/*********************************************************
diff --git a/NPLib/InitialConditions/TInitialConditions.h b/NPLib/InitialConditions/TInitialConditions.h
index 3b11f4a1ec66537f5b82053857932c29b7a74168..7a0001e05d36eeaa187a07df0512925b7bec0de5 100644
--- a/NPLib/InitialConditions/TInitialConditions.h
+++ b/NPLib/InitialConditions/TInitialConditions.h
@@ -34,6 +34,30 @@ using namespace std ;
class TInitialConditions : public TObject
{
private:
+ /*
+ // Beam
+ TLorentzVector fIC_Beam_LV;
+ double fIC_Beam_ThetaX;
+ double fIC_Beam_PhiY;
+
+ // Particle
+ // Emitted particle properties (after interactions in the target)
+ vector<Double_t> fIC_Emitted_Angle_ThetaCM;
+ // Emitted particle angles in the incident frame
+ vector<Double_t> fIC_Emitted_Angle_Theta_LocalFrame;
+ vector<Double_t> fIC_Emitted_Angle_Phi_LocalFrame;
+ // Emitted particle angles in the world frame
+ vector<Double_t> fIC_Emitted_Angle_Theta_LabFrame;
+ vector<Double_t> fIC_Emitted_Angle_Phi_LabFrame;
+ // Emittedparticle energy
+ vector<Double_t> fIC_Emitted_Energy;
+ vector<int> fIC_Process;
+ vector<string> fIC_Name;
+ */
+
+
+
+
// Incident particle properties (before interactions in the target)
// Vertex of interaction
vector<Double_t> fIC_Position_X;
diff --git a/NPLib/Tools/NPOptionManager.cxx b/NPLib/Tools/NPOptionManager.cxx
index 9ffec525f9d324592327e9800dddca83cae85227..37d595360dd2b2d485ebd9be519933571f143f7c 100644
--- a/NPLib/Tools/NPOptionManager.cxx
+++ b/NPLib/Tools/NPOptionManager.cxx
@@ -6,7 +6,7 @@
*****************************************************************************/
/*****************************************************************************
- * Original Author: A. MATTA contact address: matta@ipno.in2p3.fr *
+ * Original Author: A. MATTA contact address: matta@ipno.in2p3.fr *
* *
* Creation Date : 21/07/09 *
* Last update : *
@@ -30,137 +30,137 @@ NPOptionManager* NPOptionManager::instance = 0 ;
NPOptionManager* NPOptionManager::getInstance(int argc, char** argv)
{
- if (instance == 0) instance = new NPOptionManager(argc, argv);
-
- return instance ;
+ if (instance == 0) instance = new NPOptionManager(argc, argv);
+
+ return instance ;
}
NPOptionManager::NPOptionManager(int argc, char** argv)
{
- // Default Setting
- fDefaultReactionFileName = "defaultReaction.reaction";
- fDefaultDetectorFileName = "defaultDetector.detector";
- fDefaultOutputFileName = "myResult.root";
- fDefaultRunToReadFileName = "defaultRunToTreat.txt";
- fDefaultCalibrationFileName = "defaultCalibration.txt";
- // Assigned values
- fReactionFileName = fDefaultReactionFileName;
- fDetectorFileName = fDefaultDetectorFileName;
- fOutputFileName = fDefaultOutputFileName;
- fRunToReadFileName = fDefaultRunToReadFileName;
- fCalibrationFileName = fDefaultCalibrationFileName;
- fDisableAllBranchOption = false;
- fInputPhysicalTreeOption = false;
-
- for (int i = 0; i < argc; i++) {
- string argument = argv[i];
-
- if (argument == "-H" || argument == "-h" || argument == "--help") DisplayHelp();
-
- 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 == "--detector" && 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 == "-O" && argc >= i + 1) fOutputFileName = 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 == "--cal" && argc >= i + 1) fCalibrationFileName = argv[i+1] ;
-
- else if (argument == "-C" && argc >= i + 1) fCalibrationFileName = argv[i+1] ;
-
- else if (argument == "--disable-branch") fDisableAllBranchOption = true ;
-
- else if (argument == "--input-physical") fInputPhysicalTreeOption = true ;
-
- else if (argument == "-IP") fInputPhysicalTreeOption = true ;
-
- else ;
- }
-
- CheckArguments();
+ // Default Setting
+ fDefaultReactionFileName = "defaultReaction.reaction";
+ fDefaultDetectorFileName = "defaultDetector.detector";
+ fDefaultOutputFileName = "myResult.root";
+ fDefaultRunToReadFileName = "defaultRunToTreat.txt";
+ fDefaultCalibrationFileName = "defaultCalibration.txt";
+ // Assigned values
+ fReactionFileName = fDefaultReactionFileName;
+ fDetectorFileName = fDefaultDetectorFileName;
+ fOutputFileName = fDefaultOutputFileName;
+ fRunToReadFileName = fDefaultRunToReadFileName;
+ fCalibrationFileName = fDefaultCalibrationFileName;
+ fDisableAllBranchOption = false;
+ fInputPhysicalTreeOption = false;
+
+ for (int i = 0; i < argc; i++) {
+ string argument = argv[i];
+
+ if (argument == "-H" || argument == "-h" || argument == "--help") DisplayHelp();
+
+ 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 == "--detector" && 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 == "-O" && argc >= i + 1) fOutputFileName = 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 == "--cal" && argc >= i + 1) fCalibrationFileName = argv[i+1] ;
+
+ else if (argument == "-C" && argc >= i + 1) fCalibrationFileName = argv[i+1] ;
+
+ else if (argument == "--disable-branch") fDisableAllBranchOption = true ;
+
+ else if (argument == "--input-physical") fInputPhysicalTreeOption = true ;
+
+ else if (argument == "-IP") fInputPhysicalTreeOption = true ;
+
+ //else ;
+ }
+
+ CheckArguments();
}
void NPOptionManager::CheckArguments()
{
- CheckEventGenerator();
- CheckDetectorConfiguration();
+ CheckEventGenerator();
+ CheckDetectorConfiguration();
}
void NPOptionManager::CheckEventGenerator()
-{
- bool checkFile = true;
-
- // NPTool path
- string GlobalPath = getenv("NPTOOL");
- string StandardPath = GlobalPath + "/Inputs/EventGenerator/" + fReactionFileName;
-
- // ifstream to configfile
- ifstream ConfigFile;
-
- // test if config file is in local path
- ConfigFile.open(fReactionFileName.c_str());
- if (!ConfigFile.is_open()) {
- ConfigFile.open(StandardPath.c_str());
- if (!ConfigFile.is_open()) { // if not, assign standard path
- checkFile = false;
- }
- else {
- fReactionFileName = StandardPath;
- }
- }
- if (!checkFile && fReactionFileName != fDefaultReactionFileName) { // if file does not exist
- SendErrorAndExit("EventGenerator");
- }
-
- // close ConfigFile
- ConfigFile.close();
+{
+ bool checkFile = true;
+
+ // NPTool path
+ string GlobalPath = getenv("NPTOOL");
+ string StandardPath = GlobalPath + "/Inputs/EventGenerator/" + fReactionFileName;
+
+ // ifstream to configfile
+ ifstream ConfigFile;
+
+ // test if config file is in local path
+ ConfigFile.open(fReactionFileName.c_str());
+ if (!ConfigFile.is_open()) {
+ ConfigFile.open(StandardPath.c_str());
+ if (!ConfigFile.is_open()) { // if not, assign standard path
+ checkFile = false;
+ }
+ else {
+ fReactionFileName = StandardPath;
+ }
+ }
+ if (!checkFile && fReactionFileName != fDefaultReactionFileName) { // if file does not exist
+ SendErrorAndExit("EventGenerator");
+ }
+
+ // close ConfigFile
+ ConfigFile.close();
}
void NPOptionManager::CheckDetectorConfiguration()
{
- bool checkFile = true;
-
- // NPTool path
- string GlobalPath = getenv("NPTOOL");
- string StandardPath = GlobalPath + "/Inputs/DetectorConfiguration/" + fDetectorFileName;
-
- // ifstream to configfile
- ifstream ConfigFile;
-
- // test if config file is in local path
- ConfigFile.open(fDetectorFileName.c_str());
- if (!ConfigFile.is_open()) {
- ConfigFile.open(StandardPath.c_str());
- if (!ConfigFile.is_open()) { // if not, assign standard path
- checkFile = false;
- }
- else {
- fDetectorFileName = StandardPath;
- }
- }
- if (!checkFile && fDetectorFileName != fDefaultDetectorFileName) { // if file does not exist
- SendErrorAndExit("DetectorConfiguration");
- }
-
- // close ConfigFile
- ConfigFile.close();
+ bool checkFile = true;
+
+ // NPTool path
+ string GlobalPath = getenv("NPTOOL");
+ string StandardPath = GlobalPath + "/Inputs/DetectorConfiguration/" + fDetectorFileName;
+
+ // ifstream to configfile
+ ifstream ConfigFile;
+
+ // test if config file is in local path
+ ConfigFile.open(fDetectorFileName.c_str());
+ if (!ConfigFile.is_open()) {
+ ConfigFile.open(StandardPath.c_str());
+ if (!ConfigFile.is_open()) { // if not, assign standard path
+ checkFile = false;
+ }
+ else {
+ fDetectorFileName = StandardPath;
+ }
+ }
+ if (!checkFile && fDetectorFileName != fDefaultDetectorFileName) { // if file does not exist
+ SendErrorAndExit("DetectorConfiguration");
+ }
+
+ // close ConfigFile
+ ConfigFile.close();
}
@@ -168,26 +168,26 @@ void NPOptionManager::CheckDetectorConfiguration()
// This method tests if the input files are the default ones
bool NPOptionManager::IsDefault(const char* type) const
{
- bool result = false;
-
- string stype = type;
- if (stype == "EventGenerator") {
- if (fReactionFileName == fDefaultReactionFileName) result = true;
- }
- else if (stype == "DetectorConfiguration") {
- if (fDetectorFileName == fDefaultDetectorFileName) result = true;
- }
- else if (stype == "Calibration") {
- if (fCalibrationFileName == fDefaultCalibrationFileName) result = true;
- }
- else if (stype == "RunToTreat") {
- if (fRunToReadFileName == fDefaultRunToReadFileName) result = true;
- }
- else {
- cout << "NPOptionManager::IsDefault() unkwown keyword" << endl;
- }
-
- return result;
+ bool result = false;
+
+ string stype = type;
+ if (stype == "EventGenerator") {
+ if (fReactionFileName == fDefaultReactionFileName) result = true;
+ }
+ else if (stype == "DetectorConfiguration") {
+ if (fDetectorFileName == fDefaultDetectorFileName) result = true;
+ }
+ else if (stype == "Calibration") {
+ if (fCalibrationFileName == fDefaultCalibrationFileName) result = true;
+ }
+ else if (stype == "RunToTreat") {
+ if (fRunToReadFileName == fDefaultRunToReadFileName) result = true;
+ }
+ else {
+ cout << "NPOptionManager::IsDefault() unkwown keyword" << endl;
+ }
+
+ return result;
}
@@ -195,61 +195,61 @@ bool NPOptionManager::IsDefault(const char* type) const
// This method tests if the input files are the default ones
void NPOptionManager::SendErrorAndExit(const char* type) const
{
- string stype = type;
- if (stype == "EventGenerator") {
- cout << endl;
- cout << "********************************** Error **********************************" << endl;
- cout << "* No event generator file found in $NPTool/Inputs/EventGenerator or local directories *" << endl;
- cout << "***************************************************************************************" << endl;
- cout << endl;
- exit(1);
- }
- else if (stype == "DetectorConfiguration") {
- cout << endl;
- cout << "*********************************** Error ***********************************" << endl;
- cout << "* No detector geometry file found in $NPTool/Inputs/EventGenerator or local directories *" << endl;
- cout << "*****************************************************************************************" << endl;
- cout << endl;
- exit(1);
- }
- else if (stype == "Calibration") {
- }
- else if (stype == "RunToTreat") {
- }
- else {
- cout << "NPOptionManager::SendErrorAndExit() unkwown keyword" << endl;
- }
+ string stype = type;
+ if (stype == "EventGenerator") {
+ cout << endl;
+ cout << "********************************** Error **********************************" << endl;
+ cout << "* No event generator file found in $NPTool/Inputs/EventGenerator or local directories *" << endl;
+ cout << "***************************************************************************************" << endl;
+ cout << endl;
+ exit(1);
+ }
+ else if (stype == "DetectorConfiguration") {
+ cout << endl;
+ cout << "*********************************** Error ***********************************" << endl;
+ cout << "* No detector geometry file found in $NPTool/Inputs/EventGenerator or local directories *" << endl;
+ cout << "*****************************************************************************************" << endl;
+ cout << endl;
+ exit(1);
+ }
+ else if (stype == "Calibration") {
+ }
+ else if (stype == "RunToTreat") {
+ }
+ else {
+ cout << "NPOptionManager::SendErrorAndExit() unkwown keyword" << endl;
+ }
}
void NPOptionManager::DisplayHelp()
{
- cout << endl << "----NPOptionManager Help----" << endl << endl ;
- cout << "List of Option " << endl ;
- cout << "\t --help　-H -h\t \t \t \t \t \t \t　Display this help message" << endl ;
- 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 << 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 ;
- cout << "\t --disable-branch\t \t \t \t \t \t　Disable of branch of Input tree except the one of the detector (faster)" << endl ;
- cout << "\t --input-physical -IP\t \t \t \t \t \t　Consider the Input file is containing Physics Class instead of Data Class. Output branches associate to the detector are not activated" << endl ;
- cout << endl << endl ;
-
- // exit current program
- exit(1);
+ cout << endl << "----NPOptionManager Help----" << endl << endl ;
+ cout << "List of Option " << endl ;
+ cout << "\t --help　-H -h\t \t \t \t \t \t \t　Display this help message" << endl ;
+ 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 << 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 ;
+ cout << "\t --disable-branch\t \t \t \t \t \t　Disable of branch of Input tree except the one of the detector (faster)" << endl ;
+ cout << "\t --input-physical -IP\t \t \t \t \t \t　Consider the Input file is containing Physics Class instead of Data Class. Output branches associate to the detector are not activated" << endl ;
+ cout << endl << endl ;
+
+ // exit current program
+ exit(1);
}
void NPOptionManager::Destroy()
{
- if (instance != 0) {
- delete instance;
- instance = 0;
- }
+ if (instance != 0) {
+ delete instance;
+ instance = 0;
+ }
}
diff --git a/NPLib/Vamos/TVamosCHIODataDict.cxx b/NPLib/Vamos/TVamosCHIODataDict.cxx
index f99e9be062d532b4286918ce33443ae61badac5e..cf3cb80b115a57abd3cd95ef792dfa4505caab8b 100644
--- a/NPLib/Vamos/TVamosCHIODataDict.cxx
+++ b/NPLib/Vamos/TVamosCHIODataDict.cxx
@@ -1,5 +1,5 @@
//
-// File generated by rootcint at Wed Jul 18 14:42:27 2012
+// File generated by rootcint at Thu Nov 1 10:32:32 2012
// Do NOT change. Changes will be lost next time file is generated
//
@@ -821,12 +821,16 @@ static void G__cpp_setup_global0() {
}
static void G__cpp_setup_global1() {
+}
+
+static void G__cpp_setup_global2() {
G__resetglobalenv();
}
extern "C" void G__cpp_setup_globalTVamosCHIODataDict() {
G__cpp_setup_global0();
G__cpp_setup_global1();
+ G__cpp_setup_global2();
}
/*********************************************************
diff --git a/NPLib/Vamos/TVamosDCDataDict.cxx b/NPLib/Vamos/TVamosDCDataDict.cxx
index 350ba7911ba2dd9f5746ba14a30f923661403177..ed866b517f7dddb457870688d822d07e1da30c33 100644
--- a/NPLib/Vamos/TVamosDCDataDict.cxx
+++ b/NPLib/Vamos/TVamosDCDataDict.cxx
@@ -1,5 +1,5 @@
//
-// File generated by rootcint at Wed Jul 18 14:42:28 2012
+// File generated by rootcint at Thu Nov 1 10:32:33 2012
// Do NOT change. Changes will be lost next time file is generated
//
@@ -1081,12 +1081,16 @@ static void G__cpp_setup_global0() {
}
static void G__cpp_setup_global1() {
+}
+
+static void G__cpp_setup_global2() {
G__resetglobalenv();
}
extern "C" void G__cpp_setup_globalTVamosDCDataDict() {
G__cpp_setup_global0();
G__cpp_setup_global1();
+ G__cpp_setup_global2();
}
/*********************************************************
diff --git a/NPLib/Vamos/TVamosFingerDataDict.cxx b/NPLib/Vamos/TVamosFingerDataDict.cxx
index 8edb1bab5ece566ee163c488b0b601e9b1bb7b21..5c60db6f5b3d95cc13636fa80b607939dbcb14aa 100644
--- a/NPLib/Vamos/TVamosFingerDataDict.cxx
+++ b/NPLib/Vamos/TVamosFingerDataDict.cxx
@@ -1,5 +1,5 @@
//
-// File generated by rootcint at Wed Jul 18 14:42:27 2012
+// File generated by rootcint at Thu Nov 1 10:32:31 2012
// Do NOT change. Changes will be lost next time file is generated
//
@@ -510,12 +510,16 @@ static void G__cpp_setup_global0() {
}
static void G__cpp_setup_global1() {
+}
+
+static void G__cpp_setup_global2() {
G__resetglobalenv();
}
extern "C" void G__cpp_setup_globalTVamosFingerDataDict() {
G__cpp_setup_global0();
G__cpp_setup_global1();
+ G__cpp_setup_global2();
}
/*********************************************************
diff --git a/NPLib/Vamos/TVamosPlasticDataDict.cxx b/NPLib/Vamos/TVamosPlasticDataDict.cxx
index 1de81741e1576e8fd03cbbba3a4531552178e89a..195c100caf2a6c8530b9cf08f907ea0dd04e8631 100644
--- a/NPLib/Vamos/TVamosPlasticDataDict.cxx
+++ b/NPLib/Vamos/TVamosPlasticDataDict.cxx
@@ -1,5 +1,5 @@
//
-// File generated by rootcint at Wed Jul 18 14:42:26 2012
+// File generated by rootcint at Thu Nov 1 10:32:30 2012
// Do NOT change. Changes will be lost next time file is generated
//
@@ -548,12 +548,16 @@ static void G__cpp_setup_global0() {
}
static void G__cpp_setup_global1() {
+}
+
+static void G__cpp_setup_global2() {
G__resetglobalenv();
}
extern "C" void G__cpp_setup_globalTVamosPlasticDataDict() {
G__cpp_setup_global0();
G__cpp_setup_global1();
+ G__cpp_setup_global2();
}
/*********************************************************
diff --git a/NPSimulation/include/EventGeneratorPhaseSpace.hh b/NPSimulation/include/EventGeneratorPhaseSpace.hh
deleted file mode 100644
index a2bd958a46f3b099eeffc544e558b72fab4d55c5..0000000000000000000000000000000000000000
--- a/NPSimulation/include/EventGeneratorPhaseSpace.hh
+++ /dev/null
@@ -1,106 +0,0 @@
-#ifndef EventGeneratorPhaseSpace_H
-#define EventGeneratorPhaseSpace_H
-/*****************************************************************************
- * Copyright (C) 2009-2010 this file is part of the NPTool Project *
- * *
- * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
- * For the list of contributors see $NPTOOL/Licence/Contributors *
- *****************************************************************************/
-
-/*****************************************************************************
- * Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
- * *
- * Creation Date : Feb 2010 *
- * Last update : *
- *---------------------------------------------------------------------------*
- * Decription: *
- * This event Generator is used to simulated pure phase space event and *
- * especially evaluate the background contribution of different phase space *
- * channel during an experiment. Manage up to 18 body phase space *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-// C++ header
-#include <string>
-
-// G4 headers
-#include "G4ParticleDefinition.hh"
-
-// NPSimulation
-#include "VEventGenerator.hh"
-#include "Target.hh"
-
-// NPLib
-#include "TInitialConditions.h"
-
-using namespace std;
-
-class EventGeneratorPhaseSpace : public VEventGenerator
-{
- public: // Constructors and Destructors
-
- // Default constructor used to allocate memory
- EventGeneratorPhaseSpace();
-
- // This is the constructor to be used
- EventGeneratorPhaseSpace( G4ParticleDefinition* BeamParticle,
- G4ParticleDefinition* TargetParticle,
- vector<G4ParticleDefinition*> ReactionProducts,
- double BeamEnergy,
- double BeamEnergySpread,
- double SigmaX,
- double SigmaY,
- double SigmaThetaX,
- double SigmaPhiY );
-
- // Default Destructor
- virtual ~EventGeneratorPhaseSpace();
-
- public: // Inherit from VEventGenerator class
- void ReadConfiguration(string);
- void GenerateEvent(G4Event*, G4ParticleGun*);
- void SetTarget(Target* Target);
-
- private: // Particle Shoot Option
- bool m_ShootLight;
- bool m_ShootHeavy;
- bool m_ShootDecayProduct;
-
- private: // Value to store in the output tree : initial value of beam and reaction and phase space event weight
- TInitialConditions* m_InitConditions;
- // When the Phase Space Generator is called, the weight of the current configuration is return and stored in this variable
- // Spectrum then need to be weighted by this paramater to be realistic
- // NB: This procedure avoid long calculation time of the rejection methods previously instantiate and therefore allow simulation of manybody phase space decay
- double m_EventWeight;
-
- private: // Target Parameter
- Target* m_Target;
-
- private: // Experimental Parameter
- G4ParticleDefinition* m_BeamParticle ;
- G4ParticleDefinition* m_TargetParticle ;
- double m_BeamEnergy ;
- double m_BeamEnergySpread ;
- double m_SigmaX ;
- double m_SigmaY ;
- double m_SigmaThetaX ;
- double m_SigmaPhiY ;
- vector<G4ParticleDefinition*> m_ReactionProducts ;
-
- private: //Other
- void InitializeRootOutput() ;
-
- void SetEverything( G4ParticleDefinition* BeamParticle ,
- G4ParticleDefinition* TargetParticle ,
- vector<G4ParticleDefinition*> ReactionProducts ,
- double BeamEnergy ,
- double BeamEnergySpread ,
- double SigmaX ,
- double SigmaY ,
- double SigmaThetaX ,
- double SigmaPhiY );
-};
-
-#endif
diff --git a/NPSimulation/include/ParticleStack.hh b/NPSimulation/include/ParticleStack.hh
index 62c432e67e4f8476d079581cf5962a818c75ddf1..5dabb7492bad1f55984802ba53434a00e301cbe1 100644
--- a/NPSimulation/include/ParticleStack.hh
+++ b/NPSimulation/include/ParticleStack.hh
@@ -23,11 +23,14 @@
*****************************************************************************/
// NPS
-#include"Particle.hh"
+#include "Particle.hh"
+
+// NPL
+#include "TInitialConditions.h"
// STL
-#include<vector>
-#include<string>
+#include <vector>
+#include <string>
using namespace std;
// G4
@@ -64,7 +67,8 @@ private:
static ParticleStack* instance;
// The particle gun
G4ParticleGun* m_particleGun;
-
+ // Host the Initial conditions TObject
+ TInitialConditions* m_InitialConditions;
private: // Private Member
vector<Particle> m_ParticleStack;
diff --git a/NPSimulation/src/EventGeneratorParticleDecay.cc b/NPSimulation/src/EventGeneratorParticleDecay.cc
index 813d002113ac5d5b545cd92deac2187d15ba2aea..e454ca14d0f8d804d4574db2dbec8759536f671f 100644
--- a/NPSimulation/src/EventGeneratorParticleDecay.cc
+++ b/NPSimulation/src/EventGeneratorParticleDecay.cc
@@ -40,395 +40,395 @@ using namespace NPL;
#include "TLorentzVector.h"
#include "TVector3.h"
EventGeneratorParticleDecay::EventGeneratorParticleDecay(){
- m_ParticleStack = ParticleStack::getInstance();
+ m_ParticleStack = ParticleStack::getInstance();
}
EventGeneratorParticleDecay::~EventGeneratorParticleDecay(){
-
+
}
void EventGeneratorParticleDecay::ReadConfiguration(string Path,int Occurence){
- ////////General Reading needs////////
- string LineBuffer;
- string DataBuffer;
- istringstream LineStream;
- int TokkenOccurence = 0;
- //////// Setting needs///////
- bool ReadingStatusParticleDecay = false ;
-
- bool check_Daughter = false ;
- bool check_CrossSection = false ;
- bool check_ExcitationEnergy = false ;
- bool check_shoot = false ;
- bool check_created = false ;
-
- // Instantiate new variable for the up coming Particle
- vector<string> DaughterName;
- vector<double> ExcitationEnergy;
- vector<bool> shoot;
- string CSPath = "TGenPhaseSpace";
-
- //////////////////////////////////////////////////////////////////////////////////////////
- ifstream InputFile;
- InputFile.open(Path.c_str());
+ ////////General Reading needs////////
+ string LineBuffer;
+ string DataBuffer;
+ istringstream LineStream;
+ int TokkenOccurence = 0;
+ //////// Setting needs///////
+ bool ReadingStatusParticleDecay = false ;
+
+ bool check_Daughter = false ;
+ bool check_CrossSection = false ;
+ bool check_ExcitationEnergy = false ;
+ bool check_shoot = false ;
+ bool check_created = false ;
+
+ // Instantiate new variable for the up coming Particle
+ vector<string> DaughterName;
+ vector<double> ExcitationEnergy;
+ vector<bool> shoot;
+ string CSPath = "TGenPhaseSpace";
+
+ //////////////////////////////////////////////////////////////////////////////////////////
+ ifstream InputFile;
+ InputFile.open(Path.c_str());
+
+ if (InputFile.is_open()) {}
+
+ else {
+ return;
+ }
+
+ while (!InputFile.eof()&& !check_created) {
+ //Pick-up next line
+ getline(InputFile, LineBuffer);
- if (InputFile.is_open()) {}
-
- else {
- return;
+ if (LineBuffer.compare(0, 13, "ParticleDecay") == 0) {
+ TokkenOccurence++;
+ if(TokkenOccurence==Occurence){
+ ReadingStatusParticleDecay = true ;
+ G4cout << "///////////////////////////////////////// " << G4endl;
+ // Get the nuclei name
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ LineStream >> DataBuffer;
+ DataBuffer.erase();
+ LineStream >> DataBuffer;
+ m_MotherNucleiName = DataBuffer ;
+ G4cout << "Particle Decay for " << m_MotherNucleiName << G4endl;
+ }
}
- while (!InputFile.eof()&& !check_created) {
- //Pick-up next line
+ ///////////////////////////////
+ /// Gamma Decay case
+ while(ReadingStatusParticleDecay){
+
+ InputFile >> DataBuffer;
+ //Search for comment Symbol %
+ if (DataBuffer.compare(0, 1, "%") == 0) {
+ InputFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );
+ }
+
+ else if (DataBuffer == "Daughter=") {
+ check_Daughter = true ;
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+
getline(InputFile, LineBuffer);
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ G4cout << " Daughter: " ;
+ while(LineStream >> DataBuffer){
+ DaughterName.push_back(DataBuffer);
+ G4cout << DataBuffer << " " ;
+ }
- if (LineBuffer.compare(0, 13, "ParticleDecay") == 0) {
- TokkenOccurence++;
- if(TokkenOccurence==Occurence){
- ReadingStatusParticleDecay = true ;
- G4cout << "///////////////////////////////////////// " << G4endl;
- // Get the nuclei name
- LineStream.clear();
- LineStream.str(LineBuffer);
- LineStream >> DataBuffer;
- DataBuffer.erase();
- LineStream >> DataBuffer;
- m_MotherNucleiName = DataBuffer ;
- G4cout << "Particle Decay for " << m_MotherNucleiName << G4endl;
- }
+ G4cout << G4endl;
+
+ }
+
+ else if(DataBuffer == "ExcitationEnergy=") {
+ check_ExcitationEnergy = true;
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+
+ getline(InputFile, LineBuffer);
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ G4cout << " Excitation Energy: " ;
+ while(LineStream >> DataBuffer){
+ ExcitationEnergy.push_back( atof(DataBuffer.c_str()) );
+ G4cout << DataBuffer << " " ;
}
- ///////////////////////////////
- /// Gamma Decay case
- while(ReadingStatusParticleDecay){
-
- InputFile >> DataBuffer;
- //Search for comment Symbol %
- if (DataBuffer.compare(0, 1, "%") == 0) {
- InputFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );
- }
-
- else if (DataBuffer == "Daughter=") {
- check_Daughter = true ;
- LineStream.clear();
- LineStream.str(LineBuffer);
-
- getline(InputFile, LineBuffer);
- LineStream.clear();
- LineStream.str(LineBuffer);
- G4cout << " Daughter: " ;
- while(LineStream >> DataBuffer){
- DaughterName.push_back(DataBuffer);
- G4cout << DataBuffer << " " ;
- }
+ G4cout << G4endl;
+ }
+
+ else if(DataBuffer == "DifferentialCrossSection=") {
+ check_CrossSection = true;
+ InputFile >> CSPath ;
+ G4cout << " Cross Section: " << CSPath << G4endl;
+ }
+
+ else if(DataBuffer == "shoot=") {
+ check_shoot = true;
+ LineStream.clear();
+ LineStream.str(LineBuffer);
- G4cout << G4endl;
-
- }
-
- else if(DataBuffer == "ExcitationEnergy=") {
- check_ExcitationEnergy = true;
- LineStream.clear();
- LineStream.str(LineBuffer);
-
- getline(InputFile, LineBuffer);
- LineStream.clear();
- LineStream.str(LineBuffer);
- G4cout << " Excitation Energy: " ;
- while(LineStream >> DataBuffer){
- ExcitationEnergy.push_back( atof(DataBuffer.c_str()) );
- G4cout << DataBuffer << " " ;
- }
-
- G4cout << G4endl;
- }
-
- else if(DataBuffer == "DifferentialCrossSection=") {
- check_CrossSection = true;
- InputFile >> CSPath ;
- G4cout << " Cross Section: " << CSPath << G4endl;
- }
-
- else if(DataBuffer == "shoot=") {
- check_shoot = true;
- LineStream.clear();
- LineStream.str(LineBuffer);
-
- getline(InputFile, LineBuffer);
- LineStream.clear();
- LineStream.str(LineBuffer);
- G4cout << " Shoot Particle: " ;
- while(LineStream >> DataBuffer){
- shoot.push_back( atof(DataBuffer.c_str()) );
- G4cout << DataBuffer << " " ;
- }
-
- G4cout << G4endl;
- }
-
-
-
- //////////////////////////////////////////////////////
- // If no Token and no comment, toggle out //
- else
- {ReadingStatusParticleDecay = false; G4cout << "ERROR : Wrong Token Sequence: Getting out " << G4endl ;
- exit(1);
- }
-
- // Decay ended
- if(check_Daughter && check_shoot){
- SetDecay(DaughterName,shoot,ExcitationEnergy,CSPath);
- ReadingStatusParticleDecay = false;
- check_created=true;
- }
+ getline(InputFile, LineBuffer);
+ LineStream.clear();
+ LineStream.str(LineBuffer);
+ G4cout << " Shoot Particle: " ;
+ while(LineStream >> DataBuffer){
+ shoot.push_back( atof(DataBuffer.c_str()) );
+ G4cout << DataBuffer << " " ;
}
+
+ G4cout << G4endl;
+ }
+
+
+
+ //////////////////////////////////////////////////////
+ // If no Token and no comment, toggle out //
+ else
+ {ReadingStatusParticleDecay = false; G4cout << "ERROR : Wrong Token Sequence: Getting out " << G4endl ;
+ exit(1);
+ }
+
+ // Decay ended
+ if(check_Daughter && check_shoot){
+ SetDecay(DaughterName,shoot,ExcitationEnergy,CSPath);
+ ReadingStatusParticleDecay = false;
+ check_created=true;
+ }
}
-
- G4cout << "///////////////////////////////////////// " << G4endl;
- InputFile.close();
+ }
+
+ G4cout << "///////////////////////////////////////// " << G4endl;
+ InputFile.close();
}
void EventGeneratorParticleDecay::GenerateEvent(G4Event* anEvent){
- // Look for the decaying nucleus
- Particle decayingParticle = m_ParticleStack->SearchAndRemoveParticle(m_MotherNucleiName);
- if(decayingParticle.GetParticleDefinition()==NULL){
- G4cout << "Gamma Decay Warning: The decaying particle " << m_MotherNucleiName
- << " was not found in the particle stack " << G4endl;
- return ;
- }
-
- G4ParticleDefinition* decayingParticleDefinition = decayingParticle.GetParticleDefinition();
- // Build the decaying particle four momenta vector:
-
- double NucleiEnergy= decayingParticle.GetParticleKineticEnergy()+decayingParticleDefinition->GetPDGMass();
+ // Look for the decaying nucleus
+ Particle decayingParticle = m_ParticleStack->SearchAndRemoveParticle(m_MotherNucleiName);
+ if(decayingParticle.GetParticleDefinition()==NULL){
+ G4cout << "Gamma Decay Warning: The decaying particle " << m_MotherNucleiName
+ << " was not found in the particle stack " << G4endl;
+ return ;
+ }
+
+ G4ParticleDefinition* decayingParticleDefinition = decayingParticle.GetParticleDefinition();
+ // Build the decaying particle four momenta vector:
+
+ double NucleiEnergy= decayingParticle.GetParticleKineticEnergy()+decayingParticleDefinition->GetPDGMass();
+
+ double NucleiMomentum=sqrt(NucleiEnergy*NucleiEnergy -
+ decayingParticleDefinition->GetPDGMass()*decayingParticleDefinition->GetPDGMass());
+
+ G4ThreeVector Momentum = decayingParticle.GetParticleMomentumDirection().unit();
+
+ // Case of one particle decaying with a user given cross section
+ if(m_DifferentialCrossSection!="TGenPhaseSpace"){
+ TLorentzVector NucleiLV( NucleiMomentum*Momentum.x(),
+ NucleiMomentum*Momentum.y(),
+ NucleiMomentum*Momentum.z(),
+ NucleiEnergy);
+ // Shoot the angle in Center of Mass (CM) frame
+ G4double ThetaCM = (m_CrossSectionThetaMin + m_CrossSectionShoot->shoot() * (m_CrossSectionThetaMax - m_CrossSectionThetaMin)) * deg;
+ G4double phi = RandFlat::shoot()*2.*pi;
- double NucleiMomentum=sqrt(NucleiEnergy*NucleiEnergy -
- decayingParticleDefinition->GetPDGMass()*decayingParticleDefinition->GetPDGMass());
+ // Build daughter particule CM LV
+ // Pre compute variable for the decay
+ double M = decayingParticleDefinition->GetPDGMass();
+ double m1 = m_DaughterNuclei[0]->GetPDGMass();
+ double m2 = m_DaughterNuclei[1]->GetPDGMass();
- G4ThreeVector Momentum = decayingParticle.GetParticleMomentumDirection().unit();
+ if(M<(m1+m2))
+ cout << "Warning: Particle Decay forbiden by kinematic, no particle emitted "<<endl;
- // Case of one particle decaying with a user given cross section
- if(m_DifferentialCrossSection!="TGenPhaseSpace"){
- TLorentzVector NucleiLV( NucleiMomentum*Momentum.x(),
- NucleiMomentum*Momentum.y(),
- NucleiMomentum*Momentum.z(),
- NucleiEnergy);
- // Shoot the angle in Center of Mass (CM) frame
- G4double ThetaCM = (m_CrossSectionThetaMin + m_CrossSectionShoot->shoot() * (m_CrossSectionThetaMax - m_CrossSectionThetaMin)) * deg;
- G4double phi = RandFlat::shoot()*2.*pi;
-
- // Build daughter particule CM LV
- // Pre compute variable for the decay
- double M = decayingParticleDefinition->GetPDGMass();
- double m1 = m_DaughterNuclei[0]->GetPDGMass();
- double m2 = m_DaughterNuclei[1]->GetPDGMass();
-
- if(M<(m1+m2))
- cout << "Warning: Particle Decay forbiden by kinematic, no particle emitted "<<endl;
-
- else {
- double Energy = ( 1./(2.*M) )*( M*M + m1*m1 - m2*m2);
- double Momentum = sqrt(Energy*Energy - m1*m1);
-
- TVector3 FirstDaughterMomentum = Momentum * TVector3( sin(ThetaCM) * cos(phi),
- sin(ThetaCM) * sin(phi),
- cos(ThetaCM));
-
- TLorentzVector FirstDaughterLV(FirstDaughterMomentum,Energy);
-
- FirstDaughterLV.Boost( NucleiLV.BoostVector() );
- TLorentzVector SecondDaughterLV = NucleiLV - FirstDaughterLV;
-
- G4ThreeVector DaughterDirection = G4ThreeVector( FirstDaughterLV.X() ,
- FirstDaughterLV.Y() ,
- FirstDaughterLV.Z() );
-
- Particle FirstDaughterParticle( m_DaughterNuclei[0],
- FirstDaughterLV.E()-m_DaughterNuclei[0]->GetPDGMass(),
- DaughterDirection.unit(),
- decayingParticle.GetParticlePosition(),
- m_shoot[0]);
-
- DaughterDirection = G4ThreeVector( SecondDaughterLV.X() ,
- SecondDaughterLV.Y() ,
- SecondDaughterLV.Z() );
-
- Particle SecondDaughterParticle( m_DaughterNuclei[1],
- SecondDaughterLV.E()-m_DaughterNuclei[1]->GetPDGMass(),
- DaughterDirection.unit(),
- decayingParticle.GetParticlePosition(),
- m_shoot[1]);
-
- ParticleStack::getInstance()->AddParticleToStack(FirstDaughterParticle);
- ParticleStack::getInstance()->AddParticleToStack(SecondDaughterParticle);
- }
+ else {
+ double Energy = ( 1./(2.*M) )*( M*M + m1*m1 - m2*m2);
+ double Momentum = sqrt(Energy*Energy - m1*m1);
+
+ TVector3 FirstDaughterMomentum = Momentum * TVector3( sin(ThetaCM) * cos(phi),
+ sin(ThetaCM) * sin(phi),
+ cos(ThetaCM));
+
+ TLorentzVector FirstDaughterLV(FirstDaughterMomentum,Energy);
+
+ FirstDaughterLV.Boost( NucleiLV.BoostVector() );
+ TLorentzVector SecondDaughterLV = NucleiLV - FirstDaughterLV;
+
+ G4ThreeVector DaughterDirection = G4ThreeVector( FirstDaughterLV.X() ,
+ FirstDaughterLV.Y() ,
+ FirstDaughterLV.Z() );
+
+ Particle FirstDaughterParticle( m_DaughterNuclei[0],
+ FirstDaughterLV.E()-m_DaughterNuclei[0]->GetPDGMass(),
+ DaughterDirection.unit(),
+ decayingParticle.GetParticlePosition(),
+ m_shoot[0]);
+
+ DaughterDirection = G4ThreeVector( SecondDaughterLV.X() ,
+ SecondDaughterLV.Y() ,
+ SecondDaughterLV.Z() );
+
+ Particle SecondDaughterParticle( m_DaughterNuclei[1],
+ SecondDaughterLV.E()-m_DaughterNuclei[1]->GetPDGMass(),
+ DaughterDirection.unit(),
+ decayingParticle.GetParticlePosition(),
+ m_shoot[1]);
+
+ ParticleStack::getInstance()->AddParticleToStack(FirstDaughterParticle);
+ ParticleStack::getInstance()->AddParticleToStack(SecondDaughterParticle);
}
+ }
+
+ // Case of a TGenPhaseSpace
+ else{
+ TLorentzVector NucleiLV( NucleiMomentum*Momentum.x()/GeV,
+ NucleiMomentum*Momentum.y()/GeV,
+ NucleiMomentum*Momentum.z()/GeV,
+ NucleiEnergy/GeV);
+
+ if( !m_TPhaseSpace.SetDecay(NucleiLV, m_DaughterNuclei.size(), m_Masses) )
+ cout << "Warning: Phase Space Decay forbiden by kinematic, or more than 18 particles, no particle emitted "<<endl;
- // Case of a TGenPhaseSpace
else{
- TLorentzVector NucleiLV( NucleiMomentum*Momentum.x()/GeV,
- NucleiMomentum*Momentum.y()/GeV,
- NucleiMomentum*Momentum.z()/GeV,
- NucleiEnergy/GeV);
+ m_TPhaseSpace.Generate();
+
+ TLorentzVector* daughterLV ;
+ double KineticEnergy;
+
+ for (unsigned int i = 0 ; i < m_DaughterNuclei.size(); i++) {
- if( !m_TPhaseSpace.SetDecay(NucleiLV, m_DaughterNuclei.size(), m_Masses) )
- cout << "Warning: Phase Space Decay forbiden by kinematic, or more than 18 particles, no particle emitted "<<endl;
+ daughterLV = m_TPhaseSpace.GetDecay(i);
+ G4ThreeVector daughterDirection = G4ThreeVector( daughterLV->X() ,
+ daughterLV->Y() ,
+ daughterLV->Z() );
- else{
- m_TPhaseSpace.Generate();
-
- TLorentzVector* daughterLV ;
- double KineticEnergy;
-
- for (unsigned int i = 0 ; i < m_DaughterNuclei.size(); i++) {
-
- daughterLV = m_TPhaseSpace.GetDecay(i);
- G4ThreeVector daughterDirection = G4ThreeVector( daughterLV->X() ,
- daughterLV->Y() ,
- daughterLV->Z() );
-
- KineticEnergy = daughterLV->E()-m_Masses[i] ;
-
- Particle daughterParticle( m_DaughterNuclei[i],
- KineticEnergy*GeV,
- daughterDirection.unit(),
- decayingParticle.GetParticlePosition(),
- m_shoot[i]);
- ParticleStack::getInstance()->AddParticleToStack(daughterParticle);
- }
- }
+ KineticEnergy = daughterLV->E()-m_Masses[i] ;
+
+ Particle daughterParticle( m_DaughterNuclei[i],
+ KineticEnergy*GeV,
+ daughterDirection.unit(),
+ decayingParticle.GetParticlePosition(),
+ m_shoot[i]);
+ ParticleStack::getInstance()->AddParticleToStack(daughterParticle);
+ }
}
+ }
}
void EventGeneratorParticleDecay::SetTarget(Target* Target){
- m_Target = Target;
+ m_Target = Target;
}
void EventGeneratorParticleDecay::SetDecay(vector<string> DaughterName, vector<bool> shoot, vector<double> ExcitationEnergy, string CSPath){
+
+ // Check the validity of the given data:
+ if (DaughterName.size() != shoot.size() || (DaughterName.size() != ExcitationEnergy.size() && ExcitationEnergy.size()!=0) ) {
+ G4cout << "ERROR : Missmatching information: Getting out " << G4endl ;
+ exit(1);
+ }
+
+ if ( DaughterName.size() != 2 && CSPath!="TGenPhaseSpace" ) {
+ G4cout << "ERROR 2: Missmatching information: Getting out " << G4endl ;
+ exit(1);
+ }
+
+ m_shoot = shoot ;
+ m_ExcitationEnergy= ExcitationEnergy ;
+
+ // If the Excitation Energy Token was omitted, then it is set to zero
+ if(m_ExcitationEnergy.size()==0)
+ for (unsigned int i = 0 ; i < DaughterName.size(); i++) {
+ m_ExcitationEnergy.push_back(0);
+ }
+
+ // used check for mass and charge conservation
+ Nucleus* myNucleus = new Nucleus(m_MotherNucleiName);
+ int InitialCharge = myNucleus->GetZ() ; int FinalCharge = 0 ;
+ int InitialMass = myNucleus->GetA() ; int FinalMass = 0 ;
+ delete myNucleus;
+ for (unsigned int i = 0 ; i< DaughterName.size(); i++) {
+ if(DaughterName[i] == "p"){
+ m_DaughterNuclei.push_back(G4ParticleTable::GetParticleTable()->FindParticle("proton"));
+ FinalMass++;FinalCharge++;
+ }
- // Check the validity of the given data:
- if (DaughterName.size() != shoot.size() || (DaughterName.size() != ExcitationEnergy.size() && ExcitationEnergy.size()!=0) ) {
- G4cout << "ERROR : Missmatching information: Getting out " << G4endl ;
- exit(1);
+ else if (DaughterName[i] == "n"){
+ m_DaughterNuclei.push_back(G4ParticleTable::GetParticleTable()->FindParticle("neutron"));
+ FinalMass++;
}
- if ( DaughterName.size() != 2 && CSPath!="TGenPhaseSpace" ) {
- G4cout << "ERROR 2: Missmatching information: Getting out " << G4endl ;
- exit(1);
+ else{
+ Nucleus* myNucleus = new Nucleus(DaughterName[i]);
+ m_DaughterNuclei.push_back(G4ParticleTable::GetParticleTable()->GetIon(myNucleus->GetZ(),
+ myNucleus->GetA(),
+ m_ExcitationEnergy[i]*MeV));
+ FinalMass+=myNucleus->GetA();
+ FinalCharge+=myNucleus->GetZ();
+ delete myNucleus;
}
+ }
+
+ // Check mass and charge conservation
+ if (InitialMass!=FinalMass || InitialCharge!=FinalCharge) {
+ G4cout << "ERROR: Mass and charge are not conserved." << G4endl;
+ exit(1);
+ }
+
+ m_DaughterName = DaughterName;
+
+ m_DifferentialCrossSection = CSPath;
+ if(CSPath!="TGenPhaseSpace") {
+ unsigned int CrossSectionSize = 0;
+ double* m_CrossSection = new double[CrossSectionSize] ;
+ double m_CrossSectionThetaMin = 0;
+ double m_CrossSectionThetaMax = 0;
- m_shoot = shoot ;
- m_ExcitationEnergy= ExcitationEnergy ;
+ string GlobalPath = getenv("NPTOOL");
+ string StandardPath = GlobalPath + "/Inputs/CrossSection/" + m_DifferentialCrossSection;
+ ifstream CSFile;
+ CSFile.open( StandardPath.c_str() );
- // If the Excitation Energy Token was omitted, then it is set to zero
- if(m_ExcitationEnergy.size()==0)
- for (unsigned int i = 0 ; i < DaughterName.size(); i++) {
- m_ExcitationEnergy.push_back(0);
- }
+ if(CSFile.is_open()) cout << "Reading Cross Section File " << m_DifferentialCrossSection << endl;
- // used check for mass and charge conservation
- Nucleus* myNucleus = new Nucleus(m_MotherNucleiName);
- int InitialCharge = myNucleus->GetZ() ; int FinalCharge = 0 ;
- int InitialMass = myNucleus->GetA() ; int FinalMass = 0 ;
- delete myNucleus;
- for (unsigned int i = 0 ; i< DaughterName.size(); i++) {
- if(DaughterName[i] == "p"){
- m_DaughterNuclei.push_back(G4ParticleTable::GetParticleTable()->FindParticle("proton"));
- FinalMass++;FinalCharge++;
- }
-
- else if (DaughterName[i] == "n"){
- m_DaughterNuclei.push_back(G4ParticleTable::GetParticleTable()->FindParticle("neutron"));
- FinalMass++;
- }
-
- else{
- Nucleus* myNucleus = new Nucleus(DaughterName[i]);
- m_DaughterNuclei.push_back(G4ParticleTable::GetParticleTable()->GetIon(myNucleus->GetZ(),
- myNucleus->GetA(),
- m_ExcitationEnergy[i]*MeV));
- FinalMass+=myNucleus->GetA();
- FinalCharge+=myNucleus->GetZ();
- delete myNucleus;
- }
+ // 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( m_DifferentialCrossSection.c_str() );
+ if(CSFile.is_open()) {
+ cout << "Reading Cross Section File " << m_DifferentialCrossSection << endl;
+ }
+
+ else {
+ cout << "ERROR : Cross Section File " << m_DifferentialCrossSection << " not found" << endl;
+ exit(1);
+ }
}
- // Check mass and charge conservation
- if (InitialMass!=FinalMass || InitialCharge!=FinalCharge) {
- G4cout << "ERROR: Mass and charge are not conserved." << G4endl;
- exit(1);
+ double CSBuffer,AngleBuffer;
+ vector<double> CrossSectionBuffer;
+ m_CrossSectionThetaMin = 200;
+ m_CrossSectionThetaMax = -10;
+ while(!CSFile.eof()) {
+ CSFile >> AngleBuffer;
+ CSFile >> CSBuffer;
+ double CSFinal = CSBuffer*sin(AngleBuffer*deg);
+ CrossSectionBuffer.push_back(CSFinal);
+ if (AngleBuffer < m_CrossSectionThetaMin) m_CrossSectionThetaMin = AngleBuffer;
+ if (AngleBuffer > m_CrossSectionThetaMax) m_CrossSectionThetaMax = AngleBuffer;
}
- m_DaughterName = DaughterName;
+ CSFile.close();
- m_DifferentialCrossSection = CSPath;
- if(CSPath!="TGenPhaseSpace") {
- unsigned int CrossSectionSize = 0;
- double* m_CrossSection = new double[CrossSectionSize] ;
- double m_CrossSectionThetaMin = 0;
- double m_CrossSectionThetaMax = 0;
-
- string GlobalPath = getenv("NPTOOL");
- string StandardPath = GlobalPath + "/Inputs/CrossSection/" + m_DifferentialCrossSection;
- ifstream CSFile;
- CSFile.open( StandardPath.c_str() );
-
- if(CSFile.is_open()) cout << "Reading Cross Section File " << m_DifferentialCrossSection << 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( m_DifferentialCrossSection.c_str() );
- if(CSFile.is_open()) {
- cout << "Reading Cross Section File " << m_DifferentialCrossSection << endl;
- }
-
- else {
- cout << "ERROR : Cross Section File " << m_DifferentialCrossSection << " not found" << endl;
- exit(1);
- }
- }
-
- double CSBuffer,AngleBuffer;
- vector<double> CrossSectionBuffer;
- m_CrossSectionThetaMin = 200;
- m_CrossSectionThetaMax = -10;
- while(!CSFile.eof()) {
- CSFile >> AngleBuffer;
- CSFile >> CSBuffer;
- double CSFinal = CSBuffer*sin(AngleBuffer*deg);
- CrossSectionBuffer.push_back(CSFinal);
- if (AngleBuffer < m_CrossSectionThetaMin) m_CrossSectionThetaMin = AngleBuffer;
- if (AngleBuffer > m_CrossSectionThetaMax) m_CrossSectionThetaMax = AngleBuffer;
- }
-
- CSFile.close();
-
- m_CrossSectionSize = CrossSectionBuffer.size();
- m_CrossSectionArray = new double[CrossSectionBuffer.size()] ;
-
- for(unsigned int i = 0 ; i < m_CrossSectionSize ; i++ ) {
- m_CrossSectionArray[i] = CrossSectionBuffer[i];
- }
-
- m_CrossSectionShoot = new RandGeneral(m_CrossSectionArray, m_CrossSectionSize);
-
+ m_CrossSectionSize = CrossSectionBuffer.size();
+ m_CrossSectionArray = new double[CrossSectionBuffer.size()] ;
+
+ for(unsigned int i = 0 ; i < m_CrossSectionSize ; i++ ) {
+ m_CrossSectionArray[i] = CrossSectionBuffer[i];
}
+ m_CrossSectionShoot = new RandGeneral(m_CrossSectionArray, m_CrossSectionSize);
- else{
- // Set up the array of masses
- m_Masses = new double[m_DaughterNuclei.size()];
-
- // Mass of the daugther nuclei are set once
- for (unsigned int i = 0 ; i < m_DaughterNuclei.size(); i++) {
- m_Masses[i] = m_DaughterNuclei[i]->GetPDGMass()/GeV;
- }
-
- }
+ }
+
+
+ else{
+ // Set up the array of masses
+ m_Masses = new double[m_DaughterNuclei.size()];
- // Change the name of the decaying nucleus to G4 standard
- m_MotherNucleiName = m_ParticleStack->ChangeNameToG4Standard(m_MotherNucleiName);
+ // Mass of the daugther nuclei are set once
+ for (unsigned int i = 0 ; i < m_DaughterNuclei.size(); i++) {
+ m_Masses[i] = m_DaughterNuclei[i]->GetPDGMass()/GeV;
+ }
+ }
+
+ // Change the name of the decaying nucleus to G4 standard
+ m_MotherNucleiName = m_ParticleStack->ChangeNameToG4Standard(m_MotherNucleiName);
+
}
\ No newline at end of file
diff --git a/NPSimulation/src/EventGeneratorPhaseSpace.cc b/NPSimulation/src/EventGeneratorPhaseSpace.cc
deleted file mode 100644
index e858d6f52c65cc3bccb284e9742a5a48dffbc35f..0000000000000000000000000000000000000000
--- a/NPSimulation/src/EventGeneratorPhaseSpace.cc
+++ /dev/null
@@ -1,437 +0,0 @@
-/*****************************************************************************
- * Copyright (C) 2009-2010 this file is part of the NPTool Project *
- * *
- * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
- * For the list of contributors see $NPTOOL/Licence/Contributors *
- *****************************************************************************/
-
-/*****************************************************************************
- * Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
- * *
- * Creation Date : Feb 2010 *
- * Last update : *
- *---------------------------------------------------------------------------*
- * Decription: *
- * This event Generator is used to simulated pure phase space event and *
- * especially evaluate the background contribution of different phase space *
- * channel during an experiment. Manage up to 18 body phase space *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
-
-// C++ headers
-#include <iostream>
-#include <fstream>
-#include <limits>
-
-// G4 header defining G4 types
-#include "globals.hh"
-
-// G4 headers
-#include "G4ParticleTable.hh"
-#include "G4ParticleGun.hh"
-#include "G4RotationMatrix.hh"
-
-// G4 headers including CLHEP headers
-// for generating random numbers
-#include "Randomize.hh"
-
-// NPTool headers
-#include "EventGeneratorPhaseSpace.hh"
-#include "RootOutput.h"
-
-//Root Headers
-#include "TGenPhaseSpace.h"
-
-using namespace std;
-using namespace CLHEP;
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-EventGeneratorPhaseSpace::EventGeneratorPhaseSpace()
-{
- //------------- Default Constructor -------------
- m_InitConditions = new TInitialConditions() ;
- m_Target = new Target() ;
- m_SigmaX = 0 ;
- m_SigmaY = 0 ;
- m_SigmaThetaX = 0 ;
- m_SigmaPhiY = 0 ;
- m_EventWeight = 0 ;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-EventGeneratorPhaseSpace::~EventGeneratorPhaseSpace()
-{
- //------------- Default Destructor ------------
- delete m_InitConditions ;
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorPhaseSpace::SetTarget(Target* Target)
- {
- if(Target!=0)
- m_Target = Target;
- }
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-EventGeneratorPhaseSpace::EventGeneratorPhaseSpace( G4ParticleDefinition* BeamParticle,
- G4ParticleDefinition* TargetParticle,
- vector<G4ParticleDefinition*> ReactionProducts,
- double BeamEnergy,
- double BeamEnergySpread,
- double SigmaX,
- double SigmaY,
- double SigmaThetaX,
- double SigmaPhiY)
-{
- //------------- Constructor with nuclei names and beam energy ------------
-
- SetEverything( BeamParticle,
- TargetParticle,
- ReactionProducts,
- BeamEnergy,
- BeamEnergySpread,
- SigmaX,
- SigmaY,
- SigmaThetaX,
- SigmaPhiY);
-
- m_EventWeight = 0;
-
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorPhaseSpace::InitializeRootOutput()
-{
- RootOutput *pAnalysis = RootOutput::getInstance();
- TTree *pTree = pAnalysis->GetTree();
- pTree->Branch("InitialConditions", "TInitialConditions", &m_InitConditions);
- pTree->Branch("EventWeight",&m_EventWeight,"EventWeigt/D");
-}
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// Inherit from VEventGenerator
-
-void EventGeneratorPhaseSpace::ReadConfiguration(string Path)
-{
-////////General Reading needs////////
- string LineBuffer;
- string DataBuffer;
-
-////////Reaction Setting needs///////
- string Beam, Target, Heavy, Light, CrossSectionPath ;
- G4double BeamEnergy = 0 , BeamEnergySpread = 0 , SigmaX = 0 , SigmaY = 0 , SigmaThetaX = 0 , SigmaPhiY=0;
- vector<G4ParticleDefinition*> ReactionProducts;
- G4ParticleDefinition* BeamParticle = NULL ;
- G4ParticleDefinition* TargetParticle = NULL ;
-
- bool ReadingStatus = false ;
- bool check_Beam = false ;
- bool check_Target = false ;
- bool check_ReactionProducts = false ;
- bool check_BeamEnergy = false ;
- bool check_BeamEnergySpread = false ;
- bool check_FWHMX = false ;
- bool check_FWHMY = false ;
- bool check_EmmitanceTheta = false ;
- bool check_EmmitancePhi = false ;
-
-//////////////////////////////////////////////////////////////////////////////////////////
- ifstream ReactionFile;
- ReactionFile.open(Path.c_str());
-
- if (ReactionFile.is_open()) {}
- else {
- return;
- }
-
- while (!ReactionFile.eof()) {
- //Pick-up next line
- getline(ReactionFile, LineBuffer);
-
- if (LineBuffer.compare(0, 10, "PhaseSpace") == 0) { ReadingStatus = true ;}
-
-
- 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.compare(0, 5, "Beam=") == 0) {
- check_Beam = true ;
- ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
- int A = atoi(DataBuffer.c_str());
- ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
- int Z = atoi(DataBuffer.c_str());
- ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
- int EXX = atoi(DataBuffer.c_str())*MeV;
-
- G4ParticleDefinition* ParticleBuffer = G4ParticleTable::GetParticleTable()->GetIon(Z, A, EXX) ;
- BeamParticle = ParticleBuffer ;
-
- G4cout << "Beam : A=" << A << " Z=" << Z << " Excitation Energy = " << EXX/MeV << "MeV" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 7, "Target=") == 0) {
- check_Target = true ;
- ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
- int A = atoi(DataBuffer.c_str());
- ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
- int Z = atoi(DataBuffer.c_str());
- ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
- int EXX = atoi(DataBuffer.c_str())*MeV;
-
- G4ParticleDefinition* ParticleBuffer = G4ParticleTable::GetParticleTable()->GetIon(Z, A, EXX) ;
- TargetParticle = ParticleBuffer ;
-
- G4cout << "Target : A=" << A << " Z=" << Z << " Excitation Energy = " << EXX/MeV << "MeV" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 11, "BeamEnergy=") == 0) {
- check_BeamEnergy = true ;
- ReactionFile >> DataBuffer;
- BeamEnergy = atof(DataBuffer.c_str()) * MeV;
- G4cout << "Beam Energy " << BeamEnergy / MeV << " MeV" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 17, "BeamEnergySpread=") == 0) {
- check_BeamEnergySpread = true ;
- ReactionFile >> DataBuffer;
- BeamEnergySpread = atof(DataBuffer.c_str()) * MeV;
- G4cout << "Beam Energy Spread " << BeamEnergySpread / MeV << " MeV" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 7, "SigmaX=") == 0) {
- check_FWHMX = true ;
- ReactionFile >> DataBuffer;
- SigmaX = atof(DataBuffer.c_str()) * mm;
- G4cout << "Beam FWHM X " << SigmaX << " mm" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 7, "SigmaY=") == 0) {
- check_FWHMY = true ;
- ReactionFile >> DataBuffer;
- SigmaY = atof(DataBuffer.c_str()) * mm;
- G4cout << "Beam FWHM Y " << SigmaX << " mm" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 19, "SigmaThetaX=") == 0) {
- check_EmmitanceTheta = true ;
- ReactionFile >> DataBuffer;
- SigmaThetaX = atof(DataBuffer.c_str()) * deg;
- G4cout << "Beam Emmitance Theta " << SigmaThetaX / deg << " deg" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 17, "SigmaPhiY=") == 0) {
- check_EmmitancePhi = true ;
- ReactionFile >> DataBuffer;
- SigmaPhiY = atof(DataBuffer.c_str()) * deg;
- G4cout << "Beam Emmitance Phi " << SigmaPhiY / deg << " deg" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 13, "DecayProduct=") == 0) {
- ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
- int A = atoi(DataBuffer.c_str());
- ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
- int Z = atoi(DataBuffer.c_str());
- ReactionFile >> DataBuffer; ReactionFile >> DataBuffer;
- int EXX = atoi(DataBuffer.c_str())*MeV;
-
- G4ParticleDefinition* ParticleBuffer ;
-
- if ( A==1 && Z==0 )
- ParticleBuffer = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
- else if ( A==1 && Z==1)
- ParticleBuffer = G4ParticleTable::GetParticleTable()->FindParticle("proton");
-// else if ( A==2 && Z==0)
-// ParticleBuffer = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
-// ParticleBuffer->SetPDGMass( ParticleBuffer->GetPDGMass*2);
- else
- ParticleBuffer = G4ParticleTable::GetParticleTable()->GetIon(Z, A, EXX) ;
-
- ReactionProducts.push_back(ParticleBuffer);
-
- G4cout << "Decay Product: A=" << A << " Z=" << Z << " Excitation Energy = " << EXX/MeV << "MeV" << G4endl;
- }
-
- else if (DataBuffer.compare(0, 21, "EndOfDecayProductList") == 0) {
- check_ReactionProducts = true ;
- }
-
-
- ///////////////////////////////////////////////////
- // 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_BeamEnergy && check_BeamEnergySpread && check_FWHMX && check_FWHMY && check_EmmitanceTheta
- && check_EmmitancePhi && check_ReactionProducts)
- ReadingStatus = false ;
-
- }
-
-
- }
-
- SetEverything( BeamParticle,
- TargetParticle,
- ReactionProducts,
- BeamEnergy,
- BeamEnergySpread,
- SigmaX,
- SigmaY,
- SigmaThetaX,
- SigmaPhiY);
-
- ReactionFile.close();
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorPhaseSpace::GenerateEvent(G4Event* anEvent , G4ParticleGun* particleGun)
-{
-
- // If first time, write the DeDx table
- if(anEvent->GetEventID()==0 && m_Target!=0)
- {
- m_Target->WriteDEDXTable(m_BeamParticle ,0, m_BeamEnergy+4*m_BeamEnergySpread);
- for(unsigned int k = 0 ; k < m_ReactionProducts.size() ; k++)
- {
- m_Target->WriteDEDXTable(m_ReactionProducts[k] ,0, m_BeamEnergy+4*m_BeamEnergySpread);
- }
- }
-
- // Clear contents of Precedent event (now stored in ROOTOutput)
- m_InitConditions->Clear() ;
- m_EventWeight = 0 ;
-
- ///////////////////////////////////////////////////////////////////////
- ///// Calculate the incident beam direction as well as the vertex /////
- ///// of interaction in target and Energy Loss of the beam within /////
- ///// the target. /////
- ///////////////////////////////////////////////////////////////////////
- G4ThreeVector InterCoord;
-
- G4double Beam_thetaX = 0, Beam_phiY = 0;
- G4double Beam_theta = 0, Beam_phi = 0;
- G4double FinalBeamEnergy = 0 ;
- G4double InitialBeamEnergy = RandGauss::shoot(m_BeamEnergy, m_BeamEnergySpread);
-
- m_Target->CalculateBeamInteraction( 0, m_SigmaX, 0, m_SigmaThetaX,
- 0, m_SigmaY, 0, m_SigmaPhiY,
- InitialBeamEnergy,
- m_BeamParticle,
- InterCoord, Beam_thetaX, Beam_phiY,
- Beam_theta, Beam_phi,
- FinalBeamEnergy);
-
- m_InitConditions->SetICIncidentEnergy(FinalBeamEnergy / MeV);
-
- // write vertex position to ROOT file
- G4double x0 = InterCoord.x();
- G4double y0 = InterCoord.y();
- G4double z0 = InterCoord.z();
- m_InitConditions->SetICPositionX(x0);
- m_InitConditions->SetICPositionY(y0);
- m_InitConditions->SetICPositionZ(z0);
-
- // write emittance angles to ROOT file
- m_InitConditions->SetICIncidentEmittanceTheta(Beam_thetaX / deg);
- m_InitConditions->SetICIncidentEmittancePhi(Beam_phiY / deg);
-
- // write angles to ROOT file
- m_InitConditions->SetICIncidentAngleTheta(Beam_theta / deg);
- m_InitConditions->SetICIncidentAnglePhi(Beam_phi / deg);
-
- //////////////////////////////////////////////////////////
- //////////////// Phase Space Calculation /////////////////
- //////////////////////////////////////////////////////////
-
- // Masses array of ecay products in GeV (unit used by ROOT)
- int NumberOfReactionProducts = m_ReactionProducts.size() ;
- double* masses = new double[NumberOfReactionProducts];
-
- for(int k = 0 ; k < NumberOfReactionProducts ; k++ )
- {
- masses[k] = m_ReactionProducts[k]-> GetPDGMass()/GeV ;
- }
-
- // Kinematics of reaction
- G4double M = m_BeamParticle -> GetPDGMass();
- G4double InitialE = FinalBeamEnergy + M ;
- G4double InitialMomentumX = sqrt( InitialE*InitialE - M*M) * sin(Beam_theta) * cos(Beam_phi);
- G4double InitialMomentumY = sqrt( InitialE*InitialE - M*M) * sin(Beam_theta) * sin(Beam_phi);
- G4double InitialMomentumZ = sqrt( InitialE*InitialE - M*M) * cos(Beam_theta);
-
- TLorentzVector Initial = TLorentzVector(InitialMomentumX/GeV, InitialMomentumY/GeV, InitialMomentumZ/GeV,InitialE/GeV) + TLorentzVector(0,0,0,m_TargetParticle -> GetPDGMass() / GeV);
-
-
- // Instentiate a Phase Space Generator, with flat distrution
- TGenPhaseSpace TPhaseSpace ;
-
- if( !TPhaseSpace.SetDecay(Initial, NumberOfReactionProducts , masses) ) cout << "Warning: Phase Space Decay forbiden by kinematic, or more than 18 particles "<<endl;
-
- // Generate event and store the associate weight. Think to use this weigt to get correcte spectrum
- m_EventWeight = TPhaseSpace.Generate();
-
-
- TLorentzVector* ProductLV ;
- for ( int u = 0; u < NumberOfReactionProducts ; u++)
- {
- ProductLV = TPhaseSpace.GetDecay(u);
-
- G4ThreeVector Momentum ( ProductLV->X()*GeV ,
- ProductLV->Y()*GeV ,
- ProductLV->Z()*GeV );
- Momentum.unit() ;
-
- G4double Energy = ProductLV->E()*GeV-masses[u]*GeV ;
-
- //Set the gun to shoot
- particleGun->SetParticleDefinition(m_ReactionProducts[u]);
- particleGun->SetParticleMomentumDirection(Momentum);
- particleGun->SetParticleEnergy(Energy);
- particleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0));
- // Shoot the Daugter
- particleGun->GeneratePrimaryVertex(anEvent) ;
-
- // get theta and phi in the world frame
- G4double theta_world = Momentum.theta();
- G4double phi_world = Momentum.phi();
- if (phi_world < 1e-6) phi_world += 2*pi;
- // write angles in ROOT file
- m_InitConditions->SetICEmittedAngleThetaLabWorldFrame(theta_world / deg);
- m_InitConditions->SetICEmittedAnglePhiWorldFrame(phi_world / deg);
- m_InitConditions->SetICEmittedEnergy(Energy);
-
- }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void EventGeneratorPhaseSpace::SetEverything( G4ParticleDefinition* BeamParticle,
- G4ParticleDefinition* TargetParticle,
- vector<G4ParticleDefinition*> ReactionProducts,
- double BeamEnergy,
- double BeamEnergySpread,
- double SigmaX,
- double SigmaY,
- double SigmaThetaX,
- double SigmaPhiY)
-{
- //------------- Constructor with nuclei names and beam energy ------------
- m_BeamParticle = BeamParticle;
- m_TargetParticle = TargetParticle;
- m_ReactionProducts = ReactionProducts;
- m_BeamEnergy = BeamEnergy;
- m_BeamEnergySpread = BeamEnergySpread;
- m_SigmaX = SigmaX;
- m_SigmaY = SigmaY;
- m_SigmaThetaX = SigmaThetaX;
- m_SigmaPhiY = SigmaPhiY;
-}
-
diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc
index 90c42c4b9322440750930cdcf97ca8299e4fcc98..585b98bec92233df892b3b748ed3939252c16e9b 100644
--- a/NPSimulation/src/EventGeneratorTransfert.cc
+++ b/NPSimulation/src/EventGeneratorTransfert.cc
@@ -449,6 +449,7 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent)
// Angles
RandGeneral CrossSectionShoot(m_Reaction->GetCrossSection(), m_Reaction->GetCrossSectionSize());
G4double ThetaCM = (m_Reaction->GetCrossSectionAngleMin() + CrossSectionShoot.shoot() * (m_Reaction->GetCrossSectionAngleMax() - m_Reaction->GetCrossSectionAngleMin())) * deg;
+
G4double phi = RandFlat::shoot() * 2*pi;
//////////////////////////////////////////////////
diff --git a/NPSimulation/src/ParticleStack.cc b/NPSimulation/src/ParticleStack.cc
index dde3fa0edcea5d6d6f2f51790186e97de224a70e..ba6a037852899eac52f74759f7c5ca8cbfeb5e6e 100644
--- a/NPSimulation/src/ParticleStack.cc
+++ b/NPSimulation/src/ParticleStack.cc
@@ -20,12 +20,14 @@
* *
*****************************************************************************/
-
+// NPS
#include"ParticleStack.hh"
// G4 headers
#include "G4ParticleTable.hh"
+// NPL
+#include "RootOutput.h"
ParticleStack* ParticleStack::instance = 0 ;
@@ -45,7 +47,10 @@ ParticleStack::ParticleStack(){
m_particleGun->SetParticlePosition(G4ThreeVector(0., 0., 0.));
m_particleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.));
-}
+
+ // Instantiate the TInitialConditions object and link it to the RootOutput tree
+ m_InitialConditions = new TInitialConditions();
+ RootOutput::getInstance()->GetTree()->Branch("InitialConditions","TInitialConditions",&m_InitialConditions);}
ParticleStack::~ParticleStack(){
@@ -61,6 +66,10 @@ void ParticleStack::SetParticleStack(vector<Particle> particle_stack){
void ParticleStack::AddParticleToStack(Particle particle){
m_ParticleStack.push_back(particle);
+ m_InitialConditions->SetICPositionX(particle.GetParticlePosition().x());
+ m_InitialConditions->SetICPositionY(particle.GetParticlePosition().y());
+ m_InitialConditions->SetICPositionZ(particle.GetParticlePosition().z());
+ m_InitialConditions->SetICEmittedAngleThetaLabWorldFrame(particle.GetParticleMomentumDirection().theta());
}
Particle ParticleStack::SearchAndRemoveParticle(string name){
diff --git a/NPSimulation/src/PrimaryGeneratorAction.cc b/NPSimulation/src/PrimaryGeneratorAction.cc
index 20d7b190234ab869dff1939285afc24f00240f77..ac0a641a312cd56bc06212d7501688807921de80 100644
--- a/NPSimulation/src/PrimaryGeneratorAction.cc
+++ b/NPSimulation/src/PrimaryGeneratorAction.cc
@@ -42,7 +42,6 @@
#include "EventGeneratorTransfert.hh"
#include "EventGeneratorIsotropic.hh"
#include "EventGeneratorBeam.hh"
-#include "EventGeneratorPhaseSpace.hh"
#include "EventGeneratorGammaDecay.hh"
#include "EventGeneratorParticleDecay.hh"
@@ -77,8 +76,6 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path)
{
- bool check_TransfertToResonance = false;
- bool check_PhaseSpace = false;
bool check_Isotropic = false;
bool check_Transfert = false;
bool check_Beam = false;
@@ -100,7 +97,6 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path)
cout << "Error, Event Generator file " << Path << " found" << endl;
}
-
while (!EventGeneratorFile.eof()) {
//Pick-up next line
getline(EventGeneratorFile, LineBuffer);
@@ -135,7 +131,7 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path)
}
//Search for Transfert
- else if (LineBuffer.compare(0, 9, "Transfert") == 0 && !check_Transfert && LineBuffer.compare(0, 11, "TransfertTo") != 0) {
+ else if (LineBuffer.compare(0, 9, "Transfert") == 0 && !check_Transfert) {
check_Transfert = true;
VEventGenerator* myEventGenerator = new EventGeneratorTransfert();
EventGeneratorFile.close();
@@ -178,20 +174,8 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path)
seenToken_ParticleDecay=0;
}
}
-
-
- //Search for Transfert To Resonance
- else if (LineBuffer.compare(0, 10, "PhaseSpace") == 0 && !check_PhaseSpace) {
- check_PhaseSpace = true;
- VEventGenerator* myEventGenerator = new EventGeneratorPhaseSpace();
- EventGeneratorFile.close();
- myEventGenerator->ReadConfiguration(Path);
- EventGeneratorFile.open(Path.c_str());
- myEventGenerator->InitializeRootOutput();
- myEventGenerator->SetTarget(m_detector->GetTarget());
- m_EventGenerator.push_back(myEventGenerator);
- }
}
+
EventGeneratorFile.close();
}