diff --git a/Inputs/EventGenerator/3He.source b/Inputs/EventGenerator/3He.source
index 2c4bf8b610a8ea35d9f21ecca33a4d4f7e939488..d85a5bce8286662497a458a9fad9ce991b45c15f 100644
--- a/Inputs/EventGenerator/3He.source
+++ b/Inputs/EventGenerator/3He.source
@@ -1,13 +1,13 @@
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% An Isotropic Source to be used as EventGenerator %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
- EnergyLow= 120
- EnergyHigh= 200
+ EnergyLow= 0
+ EnergyHigh= 520
HalfOpenAngleMin= 0
- HalfOpenAngleMax= 12
+ HalfOpenAngleMax= 90
x0= 0
y0= 0
z0= 0
diff --git a/Inputs/EventGenerator/EnergyDist_proton.root b/Inputs/EventGenerator/EnergyDist_proton.root
new file mode 100644
index 0000000000000000000000000000000000000000..7c3d13b3d50b0dd2be68f12e171589a6ab4598e0
Binary files /dev/null and b/Inputs/EventGenerator/EnergyDist_proton.root differ
diff --git a/Inputs/EventGenerator/ThetaDist_proton.root b/Inputs/EventGenerator/ThetaDist_proton.root
new file mode 100644
index 0000000000000000000000000000000000000000..897d7d38dd876a4abe2a48ab1f0e8bc52d4eb587
Binary files /dev/null and b/Inputs/EventGenerator/ThetaDist_proton.root differ
diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source
index 02b1c4dfeae0658076e3c8cff6632c4995e62693..efbad0a97245a4b50c5111e894a5886a989f3aa8 100644
--- a/Inputs/EventGenerator/alpha.source
+++ b/Inputs/EventGenerator/alpha.source
@@ -4,10 +4,10 @@
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
- EnergyLow= 5
- EnergyHigh= 5
+ EnergyLow= 0
+ EnergyHigh= 585
HalfOpenAngleMin= 0
- HalfOpenAngleMax= 5
+ HalfOpenAngleMax= 90
x0= 0
y0= 0
z0= 0
diff --git a/Inputs/EventGenerator/deuton.source b/Inputs/EventGenerator/deuton.source
index fd3680952d6534a10546409e9382b10f7139cc39..2730c05f1604e224d664a63851a9296419a93228 100644
--- a/Inputs/EventGenerator/deuton.source
+++ b/Inputs/EventGenerator/deuton.source
@@ -5,9 +5,9 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 0
- EnergyHigh= 9
+ EnergyHigh= 250
HalfOpenAngleMin= 0
- HalfOpenAngleMax= 45
+ HalfOpenAngleMax= 90
x0= 0
y0= 0
z0= 0
diff --git a/Inputs/EventGenerator/proton.pbuu b/Inputs/EventGenerator/proton.pbuu
new file mode 100644
index 0000000000000000000000000000000000000000..d903318385148f479fa0892a5a41fa9e14965ba5
--- /dev/null
+++ b/Inputs/EventGenerator/proton.pbuu
@@ -0,0 +1,14 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%% An Isotropic Source to be used as EventGenerator %%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Energy are given in MeV , Position in mm %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+pBUU
+ AngleHistPath= ThetaDist_proton.root spectrum2
+ EnergyHistPath= EnergyDist_proton.root spectrum1
+ x0= 0
+ y0= 0
+ z0= 0
+ Particle= proton
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Supported particle type: proton, neutron, deuton, triton, He3 , alpha
diff --git a/Inputs/EventGenerator/proton.source b/Inputs/EventGenerator/proton.source
index fe74556d683666950e9bdeb3785c9f507634b5af..9bb49875926c71e424ea7975009273cdd77649b2 100644
--- a/Inputs/EventGenerator/proton.source
+++ b/Inputs/EventGenerator/proton.source
@@ -5,7 +5,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 0
- EnergyHigh= 150
+ EnergyHigh= 190
HalfOpenAngleMin= 0
HalfOpenAngleMax= 90
x0= 0
diff --git a/Inputs/EventGenerator/triton.source b/Inputs/EventGenerator/triton.source
index 5a9aa2a5ac98c86b273827760b695cd43ae8986d..a3430584ab09a241a519a1340fb6b5a5cc1fc32b 100644
--- a/Inputs/EventGenerator/triton.source
+++ b/Inputs/EventGenerator/triton.source
@@ -5,9 +5,9 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 0
- EnergyHigh= 12
+ EnergyHigh= 280
HalfOpenAngleMin= 0
- HalfOpenAngleMax= 45
+ HalfOpenAngleMax= 90
x0= 0
y0= 0
z0= 0
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index 0869dae2d6623c9b523d51bc6ff5a59f7c7fe753..61dd7a2fc3de22b9ea3a140f33daccfb5e5864ed 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -9,23 +9,17 @@
/*****************************************************************************
* *
- * Original Author : Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * Original Author : Pierre MORFOUACE contact: morfouac@nscl.msu.edu *
* *
- * Creation Date : March 2009 *
- * Last update : January 2011 *
+ * Creation Date : April 2016 *
+ * Last update : April 2016 *
*---------------------------------------------------------------------------*
* Decription: *
- * This class deal with Two Body transfert Reaction *
+ * This class simulate particule with a give energy *
+ * and angular distirubtion in the lab. *
* Physical parameter (Nuclei mass) are loaded from the nubtab03.asc file *
* (2003 nuclear table of isotopes mass). *
* *
- * KineRelativistic: Used in NPSimulation *
- * 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. *
* *
*---------------------------------------------------------------------------*
* Comment: *
diff --git a/NPSimulation/Core/CMakeLists.txt b/NPSimulation/Core/CMakeLists.txt
index 0549e0593fe1a81ad50d083dcd8eac5c72c648a1..ade632670b93c175f115c1b1675c8d51cff067d1 100644
--- a/NPSimulation/Core/CMakeLists.txt
+++ b/NPSimulation/Core/CMakeLists.txt
@@ -1,2 +1,2 @@
-add_library(NPSCore SHARED CalorimeterScorers.cc EventAction.cc EventGeneratorParticleDecay.cc ObsoleteGeneralScorers.cc PrimaryGeneratorAction.cc Target.cc Chamber.cc EventGeneratorBeam.cc EventGeneratorTwoBodyReaction.cc Particle.cc PrimaryGeneratorActionMessenger.cc NPSVDetector.cc DetectorConstruction.cc EventGeneratorGammaDecay.cc MaterialManager.cc ParticleStack.cc SiliconScorers.cc PhotoDiodeScorers.cc VEventGenerator.cc DetectorMessenger.cc EventGeneratorIsotropic.cc MyMagneticField.cc PhysicsList.cc SteppingVerbose.cc NPSDetectorFactory.cc RunAction.cc)
+add_library(NPSCore SHARED CalorimeterScorers.cc EventAction.cc EventGeneratorParticleDecay.cc ObsoleteGeneralScorers.cc PrimaryGeneratorAction.cc Target.cc Chamber.cc EventGeneratorBeam.cc EventGeneratorTwoBodyReaction.cc EventGeneratorpBUU.cc Particle.cc PrimaryGeneratorActionMessenger.cc NPSVDetector.cc DetectorConstruction.cc EventGeneratorGammaDecay.cc MaterialManager.cc ParticleStack.cc SiliconScorers.cc PhotoDiodeScorers.cc VEventGenerator.cc DetectorMessenger.cc EventGeneratorIsotropic.cc MyMagneticField.cc PhysicsList.cc SteppingVerbose.cc NPSDetectorFactory.cc RunAction.cc)
target_link_libraries(NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPInitialConditions -lNPInteractionCoordinates)
diff --git a/NPSimulation/Core/EventGeneratorIsotropic.cc b/NPSimulation/Core/EventGeneratorIsotropic.cc
index e1cdcd423337d9b1bf85ab7c245dda57f5934d7a..9a4f8714e801625d9d29ebc140ffd48208245409 100644
--- a/NPSimulation/Core/EventGeneratorIsotropic.cc
+++ b/NPSimulation/Core/EventGeneratorIsotropic.cc
@@ -40,14 +40,16 @@ using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorIsotropic::EventGeneratorIsotropic(){
- m_EnergyLow = 0 ;
- m_EnergyHigh = 0 ;
- m_x0 = 0 ;
- m_y0 = 0 ;
- m_z0 = 0 ;
- m_particle = NULL;
- m_ParticleStack = ParticleStack::getInstance();
- m_ExcitationEnergy = 0 ;
+ m_EnergyLow = 0 ;
+ m_EnergyHigh = 0 ;
+ m_x0 = 0 ;
+ m_y0 = 0 ;
+ m_z0 = 0 ;
+ m_SigmaX = 0;
+ m_SigmaY = 0;
+ m_particle = NULL;
+ m_ParticleStack = ParticleStack::getInstance();
+ m_ExcitationEnergy = 0 ;
}
diff --git a/NPSimulation/Core/EventGeneratorpBUU.cc b/NPSimulation/Core/EventGeneratorpBUU.cc
new file mode 100644
index 0000000000000000000000000000000000000000..eeda164098f21852232d7625cddead6974ef5a99
--- /dev/null
+++ b/NPSimulation/Core/EventGeneratorpBUU.cc
@@ -0,0 +1,256 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2013 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: Pierre MORFOUACE contact address: morfouac@nscl.msu.edu *
+ * *
+ * Creation Date : April 2016 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This event Generator is used to simulated pBUU ion Source *
+ * Very usefull to figure out Geometric Efficacity of experimental Set-Up *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+// C++
+#include<limits>
+
+// G4 headers
+#include "G4ParticleTable.hh"
+#include "G4IonTable.hh"
+// G4 headers including CLHEP headers
+// for generating random numbers
+#include "Randomize.hh"
+
+// NPS headers
+#include "EventGeneratorpBUU.hh"
+
+// NPl headers
+#include "RootOutput.h"
+#include "NPNucleus.h"
+
+using namespace CLHEP;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+EventGeneratorpBUU::EventGeneratorpBUU(){
+ m_x0 = 0 ;
+ m_y0 = 0 ;
+ m_z0 = 0 ;
+ m_SigmaX = 0;
+ m_SigmaY = 0;
+ m_particle = NULL;
+ m_ParticleStack = ParticleStack::getInstance();
+ m_ExcitationEnergy = 0 ;
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+EventGeneratorpBUU::~EventGeneratorpBUU(){
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorpBUU::ReadConfiguration(string Path,int){
+ ////////General Reading needs////////
+ string LineBuffer;
+ string DataBuffer;
+
+ bool ReadingStatus = false ;
+ bool check_AngleHistPath = false;
+ bool check_EnergyHistPath = false;
+ bool check_x0 = false ;
+ bool check_y0 = false ;
+ bool check_z0 = false ;
+ bool check_particle = false ;
+ bool check_ExcitationEnergy = false ;
+
+ ////////Reaction Setting needs///////
+ string particle ;
+ //////////////////////////////////////////////////////////////////////////////////////////
+ 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, 4, "pBUU") == 0) {
+ G4cout << "///////////////////////////////////////////////////" << G4endl ;
+ G4cout << "pBUU Source Found" << G4endl ;
+ 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== "AngleHistPath=") {
+ check_AngleHistPath = true ;
+ TString FileName,HistName;
+ ReactionFile >> FileName >> HistName;
+ G4cout << "Reading Angle Distribution file: " << FileName << G4endl;
+
+ string GlobalPath = getenv("NPTOOL");
+ TString StandardPath = GlobalPath + "/Inputs/EventGenerator/" + FileName;
+
+ TFile *f1 = new TFile(StandardPath);
+ fAngleHist = (TH1F*) f1->FindObjectAny(HistName);
+ if(!fAngleHist){
+ G4cout << "Error: Histogramm " << HistName << " not found in file " << FileName << G4endl;
+ exit(1);
+ }
+ }
+
+ else if (DataBuffer== "EnergyHistPath=") {
+ check_EnergyHistPath = true ;
+ TString FileName,HistName;
+ ReactionFile >> FileName >> HistName;
+ G4cout << "Reading Energy Distribution file: " << FileName << G4endl;
+
+ string GlobalPath = getenv("NPTOOL");
+ TString StandardPath = GlobalPath + "/Inputs/EventGenerator/" + FileName;
+
+ TFile *f2 = new TFile(StandardPath);
+ fEnergyHist = (TH1F*) f2->FindObjectAny(HistName);
+ if(!fEnergyHist){
+ G4cout << "Error: Histogramm " << HistName << " not found in file " << FileName << G4endl;
+ exit(1);
+ }
+ }
+
+
+ else if (DataBuffer == "x0=") {
+ check_x0 = true ;
+ ReactionFile >> DataBuffer;
+ m_x0 = atof(DataBuffer.c_str()) * mm;
+ G4cout << "x0 " << m_x0 << " mm" << G4endl;
+ }
+
+ else if (DataBuffer == "y0=") {
+ check_y0 = true ;
+ ReactionFile >> DataBuffer;
+ m_y0 = atof(DataBuffer.c_str()) * mm;
+ G4cout << "y0 " << m_y0 << " mm" << G4endl;
+ }
+
+ else if (DataBuffer == "z0=" ) {
+ check_z0 = true ;
+ ReactionFile >> DataBuffer;
+ m_z0 = atof(DataBuffer.c_str()) * mm;
+ G4cout << "z0 " << m_z0 << " mm" << G4endl;
+ }
+
+ else if (DataBuffer == "SigmaX=" ) {
+ ReactionFile >> DataBuffer;
+ m_SigmaX = atof(DataBuffer.c_str()) * mm;
+ G4cout << "SigmaX " << m_SigmaX << " mm" << G4endl;
+ }
+
+ else if (DataBuffer == "SigmaY=" ) {
+ ReactionFile >> DataBuffer;
+ m_SigmaY = atof(DataBuffer.c_str()) * mm;
+ G4cout << "SigmaY " << m_SigmaY << " mm" << G4endl;
+ }
+
+ else if (DataBuffer=="Particle=" || DataBuffer=="particle=") {
+ check_particle = true ;
+ ReactionFile >> m_particleName;
+ G4cout << "Particle : " << m_particleName << G4endl ;
+
+ // Case of light particle
+ if(m_particleName=="proton"){ m_particleName="1H" ; check_ExcitationEnergy = true ;}
+ else if(m_particleName=="deuton"){ m_particleName="2H" ; check_ExcitationEnergy = true ;}
+ else if(m_particleName=="triton"){ m_particleName="3H" ; check_ExcitationEnergy = true ;}
+ else if(m_particleName=="alpha") { m_particleName="4He" ; check_ExcitationEnergy = true ;}
+ else if(m_particleName=="gamma") { check_ExcitationEnergy = true ;}
+ else if(m_particleName=="neutron") { check_ExcitationEnergy = true ;}
+ else { check_ExcitationEnergy = true ;}
+
+ }
+
+ else if (DataBuffer=="ExcitationEnergy=") {
+ check_ExcitationEnergy = true ;
+ ReactionFile >> DataBuffer;
+ m_ExcitationEnergy = atof(DataBuffer.c_str()) * MeV;
+ G4cout << "ExcitationEnergy : " << m_ExcitationEnergy << G4endl ;
+ }
+
+ // If no pBUU 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_EnergyHistPath && check_AngleHistPath && check_x0 && check_y0 && check_z0 && check_particle && check_ExcitationEnergy)
+ ReadingStatus = false ;
+
+ }
+
+ }
+
+ if( !check_AngleHistPath || !check_EnergyHistPath || !check_x0 || !check_y0 || !check_z0 || !check_particle )
+ {G4cout << "ERROR : Token Sequence Incomplete, pBUU definition can not be Fonctionnal" << G4endl ; exit(1);}
+
+
+
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorpBUU::GenerateEvent(G4Event*){
+
+ if(m_particle==NULL){
+ if(m_particleName=="gamma" || m_particleName=="neutron" || m_particleName=="opticalphoton"){
+ m_particle = G4ParticleTable::GetParticleTable()->FindParticle(m_particleName.c_str());
+ }
+ else{
+ NPL::Nucleus* N = new NPL::Nucleus(m_particleName);
+ m_particle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(N->GetZ(), N->GetA(),m_ExcitationEnergy);
+ delete N;
+ }
+ }
+
+ G4double theta = fAngleHist->GetRandom()*deg;
+ G4double particle_energy = fEnergyHist->GetRandom() / MeV;
+ G4double phi = RandFlat::shoot() * 2 * pi;
+
+ // Direction of particle, energy and laboratory angle
+ G4double momentum_x = sin(theta) * cos(phi) ;
+ G4double momentum_y = sin(theta) * sin(phi) ;
+ G4double momentum_z = cos(theta) ;
+
+ G4double x0 = RandGauss::shoot(m_x0,m_SigmaX);
+ G4double y0 = RandGauss::shoot(m_y0,m_SigmaY);
+
+ Particle particle(m_particle, theta,particle_energy,G4ThreeVector(momentum_x, momentum_y, momentum_z),G4ThreeVector(x0, y0, m_z0));
+
+
+ m_ParticleStack->AddParticleToStack(particle);
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorpBUU::InitializeRootOutput(){
+
+}
diff --git a/NPSimulation/Core/EventGeneratorpBUU.hh b/NPSimulation/Core/EventGeneratorpBUU.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ff4b6096dd03bb97b789a498b49276d19efdca0f
--- /dev/null
+++ b/NPSimulation/Core/EventGeneratorpBUU.hh
@@ -0,0 +1,66 @@
+#ifndef EventGeneratorpBUU_h
+#define EventGeneratorpBUU_h
+/*****************************************************************************
+ * Copyright (C) 2009-2013 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: Pierre MORFOUACE contact address: morfouac@nscl.msu.edu *
+ * *
+ * Creation Date : April 2016 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This event Generator is used to simulated ion Source from pBUU simulation *
+ * Very usefull to figure out Geometric Efficacity of experimental Set-Up *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+// C++ header
+#include <string>
+using namespace std;
+
+using namespace CLHEP;
+
+// G4 headers
+#include "G4Event.hh"
+
+// NPS headers
+#include "VEventGenerator.hh"
+#include "ParticleStack.hh"
+
+// ROOT Headers
+#include "TH1F.h"
+#include "TFile.h"
+
+class EventGeneratorpBUU : public NPS::VEventGenerator{
+public: // Constructor and destructor
+ EventGeneratorpBUU() ;
+ virtual ~EventGeneratorpBUU();
+
+public: // Inherit from VEventGenerator Class
+ void ReadConfiguration(string,int) ;
+ void GenerateEvent(G4Event*) ;
+ void InitializeRootOutput() ;
+
+private: // Source parameter from input file
+ G4double m_x0 ; // Vertex Position X
+ G4double m_y0 ; // Vertex Position Y
+ G4double m_z0 ; // Vertex Position Z
+ G4double m_SigmaX ;
+ G4double m_SigmaY ;
+ G4ParticleDefinition* m_particle ; // Kind of particle to shoot
+ G4double m_ExcitationEnergy ; // Excitation energy of the emitted particle
+ string m_particleName ;
+ ParticleStack* m_ParticleStack ;
+
+ TH1F* fAngleHist;
+ TH1F* fEnergyHist;
+
+};
+#endif
diff --git a/NPSimulation/Core/ParticleStack.cc b/NPSimulation/Core/ParticleStack.cc
index b93028c477e03dec42e8b6d4ae2c9a33b12f7608..09c012db238298ccc6f24a705199e69049fcdda0 100644
--- a/NPSimulation/Core/ParticleStack.cc
+++ b/NPSimulation/Core/ParticleStack.cc
@@ -96,7 +96,8 @@ void ParticleStack::AddBeamParticleToStack(Particle& particle){
m_ParticleStack.push_back(particle);
// Incident beam parameter
m_InitialConditions-> SetIncidentParticleName (particle.GetParticleDefinition()->GetParticleName());
- m_InitialConditions-> SetIncidentInitialKineticEnergy (particle. GetParticleThetaCM());
+ //m_InitialConditions-> SetIncidentInitialKineticEnergy (particle. GetParticleThetaCM());
+ m_InitialConditions-> SetIncidentInitialKineticEnergy (particle. GetParticleKineticEnergy());
G4ThreeVector U(1,0,0);
G4ThreeVector V(0,1,0);
diff --git a/NPSimulation/Core/PrimaryGeneratorAction.cc b/NPSimulation/Core/PrimaryGeneratorAction.cc
index 321070127e8a089c5332d75394e1dd0e765b56b9..a857be5a17059771749c109e98f61760c8f8b776 100644
--- a/NPSimulation/Core/PrimaryGeneratorAction.cc
+++ b/NPSimulation/Core/PrimaryGeneratorAction.cc
@@ -42,6 +42,7 @@
// Event Generator Class
#include "EventGeneratorTwoBodyReaction.hh"
#include "EventGeneratorIsotropic.hh"
+#include "EventGeneratorpBUU.hh"
#include "EventGeneratorBeam.hh"
#include "EventGeneratorGammaDecay.hh"
#include "EventGeneratorParticleDecay.hh"
@@ -79,9 +80,10 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path){
- bool check_Isotropic = false;
- bool check_TwoBodyReaction = false;
- bool check_Beam = false;
+ bool check_Isotropic = false;
+ bool check_pBUU = false;
+ bool check_TwoBodyReaction = false;
+ bool check_Beam = false;
// You can have more than one of those
int alreadyiInstantiate_GammaDecay = 0;
@@ -122,6 +124,17 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path){
m_EventGenerator.push_back(myEventGenerator);
}
+ //Search for pBUU source
+ else if (LineBuffer.compare(0, 4, "pBUU") == 0 && !check_pBUU) {
+ check_pBUU = true;
+ NPS::VEventGenerator* myEventGenerator = new EventGeneratorpBUU();
+ EventGeneratorFile.close();
+ myEventGenerator->ReadConfiguration(Path);
+ EventGeneratorFile.open(Path.c_str());
+ myEventGenerator->InitializeRootOutput();
+ m_EventGenerator.push_back(myEventGenerator);
+ }
+
//Search for Beam
else if (LineBuffer.compare(0, 4, "Beam") == 0 && !check_Beam) {
check_Beam = true;
diff --git a/NPSimulation/PhysicsListOption.txt b/NPSimulation/PhysicsListOption.txt
index 6e6fc7c040fc7fe15a6a81a2c6400d52d2ec1711..42fbf0b5f36097997a801301737c23b68c197e38 100644
--- a/NPSimulation/PhysicsListOption.txt
+++ b/NPSimulation/PhysicsListOption.txt
@@ -1,10 +1,10 @@
EmPhysicsList Option4
DefaultCutOff 1
-IonBinaryCascadePhysics 0
+IonBinaryCascadePhysics 1
EmExtraPhysics 0
HadronElasticPhysics 0
StoppingPhysics 0
OpticalPhysics 0
-HadronPhysicsQGSP_BIC_HP 0
-Decay 0
+HadronPhysicsQGSP_BIC_HP 1
+Decay 1
diff --git a/Projects/Lassa/Analysis.cxx b/Projects/Lassa/Analysis.cxx
index d8bb65d16c2bc3b7bd113b1caab3c4da9c3532a7..0c1783034329a45b6bd52b0ba04bd65c354276e1 100644
--- a/Projects/Lassa/Analysis.cxx
+++ b/Projects/Lassa/Analysis.cxx
@@ -40,10 +40,18 @@ void Analysis::Init(){
InitialConditions = new TInitialConditions();
InitOutputBranch();
InitInputBranch();
- EDelta = 0.;
totalEvents = 0;
detectedEvents = 0;
peakEvents = 0;
+
+ f_proton = new TF1("f_proton","1 - TMath::Exp(-(x-0.0918309)/27.7746)",0,10);
+ f_deuton = new TF1("f_deuton","1 - TMath::Exp(-(x-0.0434552)/21.134)",0,10);
+ f_triton = new TF1("f_triton","1 - TMath::Exp(-(x-0.0353106)/17.9059)",0,10);
+
+ // Energy loss table: the G4Table are generated by the simulation
+ Proton_CsI = EnergyLoss("proton_CsI.G4table","G4Table",100 );
+ Deuton_CsI = EnergyLoss("deuteron_CsI.G4table","G4Table",100 );
+ Triton_CsI = EnergyLoss("triton_CsI.G4table","G4Table",100 );
}
////////////////////////////////////////////////////////////////////////////////
@@ -53,18 +61,17 @@ void Analysis::TreatEvent(){
totalEvents++;
- double BeamEnergy = 120.0;
- //double BeamEnergy = InitialConditions->GetIncidentInitialKineticEnergy();
+
+ double ParticleEnergy = InitialConditions->GetKineticEnergy(0);
- //cout << BeamEnergy << endl;
- EDelta = 2.0;
+ double EDelta = 2.0;
- //cout << Lassa->ThickSi_E.size() << endl;
+ if(Lassa->ThickSi_E.size()>0) InitialEnergy = ParticleEnergy;
///////////////////////////LOOP on Lassa Hit//////////////////////////////////
if(Lassa->ThickSi_E.size() == 1){
-
detectedEvents++;
+
//for(unsigned int countLassa = 0 ; countLassa < Lassa->ThickSi_E.size(); countLassa++){
TelescopeNumber = Lassa->TelescopeNumber[0];
@@ -80,33 +87,58 @@ void Analysis::TreatEvent(){
double Y_target = InitialConditions->GetIncidentPositionY();
double Z_target = InitialConditions->GetIncidentPositionZ();
- TVector3 PositionOnTarget = TVector3(X_target,Y_target,Z_target);
+ TVector3 PositionOnTarget = TVector3(0,0,0);
+ //TVector3 PositionOnTarget = TVector3(X_target,Y_target,Z_target);
TVector3 HitDirection = PositionOnLassa-PositionOnTarget;
TVector3 HitDirectionUnit = HitDirection.Unit();
+ TVector3 BeamDirection = TVector3(0,0,1);
//TVector3 BeamDirection = InitialConditions->GetBeamDirection();
//double XBeam = BeamDirection.X();
//double YBeam = BeamDirection.Y();
//double ZBeam = BeamDirection.Z();
- //ThetaLab = BeamDirection.Angle(HitDirection);
- ThetaLab = 0;//ThetaLab/deg;
-
- PhiLab = HitDirection.Phi()/deg;
+ ThetaLab = BeamDirection.Angle(HitDirection);
+ PhiLab = HitDirection.Phi();
E_ThickSi = Lassa->ThickSi_E[0];
-
- if(Lassa->CsI_E.size()>=1){
-
- E_CsI = Lassa->CsI_E[0];
+ ELab = E_ThickSi;
+ ELab_nucl = E_ThickSi;
+
+ //R_alpha = (4.59+1.192*pow(E_ThickSi,1.724));//*cos(ThetaLab);//Lise
+ R_alpha = (4.90+1.17883*pow(E_ThickSi,1.73497));//*cos(ThetaLab);//Geant4
- ELab = E_ThickSi + E_CsI;
+ if(Lassa->CsI_E.size()==1){
+ E_CsI = Lassa->CsI_E[0];
+ ELab += E_CsI;
+ //Try to simulate the nuclear reaction loss
+ //ThicknessCsI = Proton_CsI.EvaluateMaterialThickness(0*MeV, Lassa->CsI_E[0]*MeV, 200*millimeter, 0.1*millimeter);
+ //ThicknessCsI = Deuton_CsI.EvaluateMaterialThickness(0*MeV, Lassa->CsI_E[0]*MeV, 200*millimeter, 0.1*millimeter);
+ //ThicknessCsI = Triton_CsI.EvaluateMaterialThickness(0*MeV, Lassa->CsI_E[0]*MeV, 200*millimeter, 0.1*millimeter);
+
+ //double eval = f_proton->Eval(ThicknessCsI/10);
+ //double eval = f_deuton->Eval(ThicknessCsI/10);
+ /*double eval = f_triton->Eval(ThicknessCsI/10);
+ double Random_value = Rand.Uniform(0,1);
+
+ if(Random_value>eval) ELab_nucl += E_CsI;
+ else ELab_nucl += Rand.Uniform(0,E_CsI);*/
+
+ }
- if(ELab > (BeamEnergy-EDelta) && ELab < (BeamEnergy+EDelta) && ELab > 0.){
- peakEvents++;
- }
+ if(fabs(ParticleEnergy-ELab)<EDelta){
+ peakEvents++;
}
+ else{
+ ELab = -100;
+ //E_CsI = -100;
+ }
+
+ if(fabs(ParticleEnergy-ELab_nucl)>EDelta) ELab_nucl = -100;
+
+ ThetaLab = ThetaLab/deg;
+ PhiLab = PhiLab/deg;
}
}
////////////////////////////////////////////////////////////////////////////////
@@ -128,13 +160,15 @@ void Analysis::End(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() {
- //RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThicknessCsI",&ThicknessCsI,"ThicknessCsI/D");
RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
+ RootOutput::getInstance()->GetTree()->Branch("ELab_nucl",&ELab_nucl,"ELab_nucl/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab,"PhiLab/D");
- // RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
- // RootOutput::getInstance()->GetTree()->Branch("totalEvents",&totalEvents,"totalEvents/I");
- // RootOutput::getInstance()->GetTree()->Branch("detectedEvents",&detectedEvents,"detectedEvents/I");
+ RootOutput::getInstance()->GetTree()->Branch("InitialEnergy",&InitialEnergy,"InitialEnergy/D");
+ RootOutput::getInstance()->GetTree()->Branch("E_ThickSi",&E_ThickSi,"E_ThickSi/D");
+ RootOutput::getInstance()->GetTree()->Branch("E_CsI",&E_CsI,"E_CsI/D");
+ RootOutput::getInstance()->GetTree()->Branch("R_alpha",&R_alpha,"R_alpha/D");
// RootOutput::getInstance()->GetTree()->Branch("peakEvents",&peakEvents,"peakEvents/I");
}
@@ -147,15 +181,19 @@ void Analysis::InitInputBranch(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){
- E_ThickSi = 0.;
- E_CsI = 0.;
- ELab = 0.;
- ThetaLab = -1000;
- PhiLab = -1000;
- X = -1000;
- Y = -1000;
- Z = -1000;
+ E_ThickSi = -100;
+ E_CsI =-100;
+ ELab = -100;
+ ELab_nucl = -100;
+ ThetaLab = -100;
+ PhiLab = -100;
+ X = -100;
+ Y = -100;
+ Z = -100;
TelescopeNumber = -1;
+ InitialEnergy = -100;
+ ThicknessCsI = -1;
+ R_alpha = -100;
}
////////////////////////////////////////////////////////////////////////////////
diff --git a/Projects/Lassa/Analysis.h b/Projects/Lassa/Analysis.h
index a4de214c71bde1f5f8a55d542b94251cb5127317..c60f8cb4d333276e26aaf7e7a6330272bc4db923 100644
--- a/Projects/Lassa/Analysis.h
+++ b/Projects/Lassa/Analysis.h
@@ -27,6 +27,7 @@
#include "NPEnergyLoss.h"
#include "NPReaction.h"
#include "TRandom3.h"
+#include "TF1.h"
class Analysis: public NPL::VAnalysis{
public:
Analysis();
@@ -43,14 +44,17 @@ class Analysis: public NPL::VAnalysis{
private:
double ELab;
+ double ELab_nucl;
double E_ThickSi;
double E_CsI;
- double EDelta;
double PhiLab;
double ThetaLab;
double X,Y,Z;
double TelescopeNumber;
double thresholdEnergy;
+ double InitialEnergy;
+ double ThicknessCsI;
+ double R_alpha;
int totalEvents;
int detectedEvents;
@@ -62,6 +66,15 @@ class Analysis: public NPL::VAnalysis{
// intermediate variable
TRandom3 Rand;
+ TF1* f_proton;
+ TF1* f_deuton;
+ TF1* f_triton;
+
+ //Energy loss table
+ NPL::EnergyLoss Proton_CsI;
+ NPL::EnergyLoss Deuton_CsI;
+ NPL::EnergyLoss Triton_CsI;
+
TLassaPhysics* Lassa;
TInitialConditions* InitialConditions;
};
diff --git a/Projects/Lassa/RunToTreat.txt b/Projects/Lassa/RunToTreat.txt
index a50a6ea6708811cf5ff58f7e738e4e9b6aa4d9e9..599a23270ad6b7aeec9cd5f184667e806de15bdb 100755
--- a/Projects/Lassa/RunToTreat.txt
+++ b/Projects/Lassa/RunToTreat.txt
@@ -1,5 +1,7 @@
TTreeName
SimulatedTree
RootFileName
- ../../Outputs/Simulation/dumb.root
+ ../../Outputs/Simulation/lassa_pid_p.root
+ ../../Outputs/Simulation/lassa_pid_d.root
+ ../../Outputs/Simulation/lassa_pid_t.root
diff --git a/Projects/macros/GeometricalEfficiency.C b/Projects/macros/GeometricalEfficiency.C
index 843843e63dc4bc81ac2ce3d4080fc581be289048..b39ffb055f7385053fc2faab86f4b8a1e8777319 100644
--- a/Projects/macros/GeometricalEfficiency.C
+++ b/Projects/macros/GeometricalEfficiency.C
@@ -48,7 +48,7 @@
using namespace std;
-void GeometricalEfficiency(const char * fname = "myResult"){
+void GeometricalEfficiency(const char * fname = "lassa_proton"){
// Open output ROOT file from NPTool simulation run
TString path = gSystem->Getenv("NPTOOL");
path += "/Outputs/Simulation/";