From 948135344f8102045cfaa79756d447b061f697d1 Mon Sep 17 00:00:00 2001
From: matta <matta@npt>
Date: Wed, 30 Jan 2013 15:43:46 +0000
Subject: [PATCH] * Checking how everything work together, - fixing couple of
bug - remove some warning - enhance performance
---
Inputs/EventGenerator/10He.reaction | 10 +-
.../InitialConditions/TInitialConditions.cxx | 123 ++++++------
NPLib/InitialConditions/TInitialConditions.h | 186 ++++++------------
NPLib/Physics/NPBeam.cxx | 22 ++-
NPLib/Physics/NPBeam.h | 10 +-
NPLib/Physics/NPFunction.cxx | 35 +---
NPLib/Physics/NPReaction.cxx | 18 +-
NPLib/Physics/NPReaction.h | 3 +-
NPSimulation/Simulation.cc | 5 +-
.../include/EventGeneratorIsotropic.hh | 3 -
NPSimulation/include/Particle.hh | 55 +++---
NPSimulation/include/ParticleStack.hh | 7 +-
NPSimulation/src/EventGeneratorBeam.cc | 18 +-
NPSimulation/src/EventGeneratorGammaDecay.cc | 3 +-
NPSimulation/src/EventGeneratorIsotropic.cc | 29 +--
.../src/EventGeneratorParticleDecay.cc | 7 +-
.../src/EventGeneratorTwoBodyReaction.cc | 28 +--
NPSimulation/src/Particle.cc | 65 ++++--
NPSimulation/src/ParticleStack.cc | 63 ++++--
19 files changed, 323 insertions(+), 367 deletions(-)
diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction
index 280f56a4f..4e7f72f0d 100644
--- a/Inputs/EventGenerator/10He.reaction
+++ b/Inputs/EventGenerator/10He.reaction
@@ -7,8 +7,8 @@ Beam
SigmaEnergy= 0
SigmaThetaX= 0.01
SigmaPhiY= 0.01
- SigmaX= 0.01
- SigmaY= 0.01
+ SigmaX= 5
+ SigmaY= 5
MeanThetaX= 0
MeanPhiY= 0
MeanX= 0
@@ -32,7 +32,7 @@ TwoBodyReaction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-ParticleDecay 10He
+%ParticleDecay 10He
Daughter= 9He n
ExcitationEnergy= 0.5 0
DifferentialCrossSection= flat.txt
@@ -40,7 +40,7 @@ ParticleDecay 10He
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-GammaDecay 9He
+%GammaDecay 9He
Cascade
BranchingRatio= 30
Energies= 0.1
@@ -51,7 +51,7 @@ GammaDecay 9He
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-ParticleDecay 9He
+%ParticleDecay 9He
Daughter= 8He n
DifferentialCrossSection= flat.txt
shoot= 1 1
diff --git a/NPLib/InitialConditions/TInitialConditions.cxx b/NPLib/InitialConditions/TInitialConditions.cxx
index 6dbd3411f..f19eb1df2 100644
--- a/NPLib/InitialConditions/TInitialConditions.cxx
+++ b/NPLib/InitialConditions/TInitialConditions.cxx
@@ -30,69 +30,72 @@ using namespace std;
ClassImp(TInitialConditions)
-TInitialConditions::TInitialConditions()
-{
- // Default constructor
-
- Clear();
+TInitialConditions::TInitialConditions(){
+ // Default constructor
+ Clear();
}
-TInitialConditions::~TInitialConditions()
-{}
-
-void TInitialConditions::Clear()
-{
- // Incident particle properties (before interactions in the target)
- // Vertex of interaction
- fIC_Position_X.clear();
- fIC_Position_Y.clear();
- fIC_Position_Z.clear();
- // Theta and Phi angles for the emittance
- fIC_Incident_Emittance_Theta.clear();
- fIC_Incident_Emittance_Phi.clear();
- // Incident particle angles
- fIC_Incident_Angle_Theta.clear();
- fIC_Incident_Angle_Phi.clear();
- // Incident particle energy
- fIC_Incident_Energy.clear();
-
- // Emitted particle angles
- fIC_Emitted_Angle_ThetaCM.clear();
- fIC_Emitted_Angle_ThetaLab_IncidentFrame.clear();
- fIC_Emitted_Angle_Phi_IncidentFrame.clear();
- fIC_Emitted_Angle_ThetaLab_WorldFrame.clear();
- fIC_Emitted_Angle_Phi_WorldFrame.clear();
- // Emitted particle energy
- fIC_Emitted_Energy.clear();
+TInitialConditions::~TInitialConditions(){
}
+void TInitialConditions::Clear(){
+ // Incident beam parameter
+ fIC_Incident_Particle_Name="";
+ fIC_Incident_Emittance_Phi=-1;
+ fIC_Incident_Emittance_Theta=-1;
+ fIC_Incident_Kinetic_Energy=-1;
+
+ // Beam status at the initial interaction point
+ fIC_Interaction_Kinetic_Energy=-1;
+ fIC_Interaction_Position_X=-1;
+ fIC_Interaction_Position_Y=-1;
+ fIC_Interaction_Position_Z=-1;
+
+ // emmitted particle
+ fIC_Particle_Name.clear();
+ fIC_Process_Name.clear();
+ fIC_ThetaCM.clear();
+ fIC_Kinetic_Energy.clear();
+ fIC_Momentum_Direction_X.clear();
+ fIC_Momentum_Direction_Y.clear();
+ fIC_Momentum_Direction_Z.clear();
+}
-
-void TInitialConditions::Dump() const
-{
- cout << "XXXXXXXXXXXXX Initial conditions XXXXXXXXXXXXXXXX" << endl;
-
- cout << "Vertex position : " << endl;
- cout << "\tX : " << fIC_Position_X[0] << endl;
- cout << "\tY : " << fIC_Position_Y[0] << endl;
- cout << "\tZ : " << fIC_Position_Z[0] << endl;
- cout << "Theta and Phi angles for the emittance : " << endl;
- cout << "\tTheta : " << fIC_Incident_Emittance_Theta[0] << endl;
- cout << "\tPhi : " << fIC_Incident_Emittance_Phi[0] << endl;
- cout << "Incident particle angles : " << endl;
- cout << "\tTheta : " << fIC_Incident_Angle_Theta[0] << endl;
- cout << "\tPhi : " << fIC_Incident_Angle_Phi[0] << endl;
- cout << "Incident energy : " << endl;
- cout << "\tEnergy : " << fIC_Incident_Energy[0] << endl;
- cout << endl;
- cout << "Emitted particle angles : " << endl;
- cout << "\tTheta CM : " << fIC_Emitted_Angle_ThetaCM[0] << endl;
- cout << " Incident frame :" << endl;
- cout << "\tTheta Lab : " << fIC_Emitted_Angle_ThetaLab_IncidentFrame[0] << endl;
- cout << "\tPhi : " << fIC_Emitted_Angle_Phi_IncidentFrame[0] << endl;
- cout << " World frame :" << endl;
- cout << "\tTheta Lab : " << fIC_Emitted_Angle_ThetaLab_WorldFrame[0] << endl;
- cout << "\tPhi : " << fIC_Emitted_Angle_Phi_WorldFrame[0] << endl;
- cout << "Emitted energy : " << endl;
- cout << "\tEnergy : " << fIC_Emitted_Energy[0] << endl;
+void TInitialConditions::Dump() const{
+ cout << "--------- Initial Condition Dump ---------" << endl ;
+
+ // Incident beam parameter
+ cout << "\t ---- Incident Beam ---- " << endl;
+ cout << "\t Particle Name: " << fIC_Incident_Particle_Name << endl;
+ cout << "\t Energy: " << fIC_Incident_Kinetic_Energy << endl;
+ cout << "\t Theta_X: " << fIC_Incident_Emittance_Theta << endl;
+ cout << "\t Phi_Y: " << fIC_Incident_Emittance_Phi << endl;
+
+
+
+ // Beam status at the initial interaction point
+ cout << "\t ---- Interaction Point ---- " << endl;
+ cout << "\t Energy: " << fIC_Interaction_Kinetic_Energy << endl;
+ cout << "\t Position: ( "
+ << fIC_Interaction_Position_X << " ; "
+ << fIC_Interaction_Position_Y << " ; "
+ << fIC_Interaction_Position_Z << ")" << endl;
+
+
+ // emmitted particle
+ unsigned int size = fIC_Particle_Name.size();
+ for(unsigned int i = 0 ; i < size; i ++){
+ cout << "\t ---- Particle " << i << " ---- " << endl;
+ cout << "\t Particle Name" << fIC_Particle_Name[i] << endl;
+ // cout << "\t Process Name" << fIC_Process_Name[i] << endl;
+ cout << "\t Theta CM" << fIC_ThetaCM[i] << endl;
+ cout << "\t Energy" << fIC_Kinetic_Energy[i] << endl;
+ cout << "\t Momentum Direction: ( "
+ << fIC_Momentum_Direction_X[i] << " ; "
+ << fIC_Momentum_Direction_Y[i] << " ; "
+ << fIC_Momentum_Direction_Z[i] << ")" << endl;
+
+ }
+
+
}
diff --git a/NPLib/InitialConditions/TInitialConditions.h b/NPLib/InitialConditions/TInitialConditions.h
index 699a3a4e9..4d8674cda 100644
--- a/NPLib/InitialConditions/TInitialConditions.h
+++ b/NPLib/InitialConditions/TInitialConditions.h
@@ -26,22 +26,26 @@
#ifndef __INITIALCONDITIONS__
#define __INITIALCONDITIONS__
+// STL Header
#include <vector>
-#include "TObject.h"
+#include <string>
+
using namespace std ;
+// Root Header
+#include "TObject.h"
-class TInitialConditions : public TObject
-{
+class TInitialConditions : public TObject{
private:
- /*
- // Incident beam parament
+
+ // Incident beam parameter
+ string fIC_Incident_Particle_Name;
double fIC_Incident_Emittance_Phi;
double fIC_Incident_Emittance_Theta;
- double fIC_Incident_Initial_Energy;
+ double fIC_Incident_Kinetic_Energy;
// Beam status at the initial interaction point
- double fIC_Interaction_Energy;
+ double fIC_Interaction_Kinetic_Energy;
double fIC_Interaction_Position_X;
double fIC_Interaction_Position_Y;
double fIC_Interaction_Position_Z;
@@ -55,122 +59,60 @@ private:
vector<double> fIC_Momentum_Direction_Y;
vector<double> fIC_Momentum_Direction_Z;
- */
-
- /*
- // 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;
- vector<Double_t> fIC_Position_Y;
- vector<Double_t> fIC_Position_Z;
- // Theta and Phi angles for the emittance
- vector<Double_t> fIC_Incident_Emittance_Theta;
- vector<Double_t> fIC_Incident_Emittance_Phi;
- // Incident particle angles
- vector<Double_t> fIC_Incident_Angle_Theta;
- vector<Double_t> fIC_Incident_Angle_Phi;
- // Incident particle energy
- vector<Double_t> fIC_Incident_Energy;
-
- // 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_ThetaLab_IncidentFrame;
- vector<Double_t> fIC_Emitted_Angle_Phi_IncidentFrame;
- // Emitted particle angles in the world frame
- vector<Double_t> fIC_Emitted_Angle_ThetaLab_WorldFrame;
- vector<Double_t> fIC_Emitted_Angle_Phi_WorldFrame;
- // Emittedparticle energy
- vector<Double_t> fIC_Emitted_Energy;
-
-
public:
- TInitialConditions();
- virtual ~TInitialConditions();
-
- void Clear();
- void Clear(const Option_t*) {};
- void Dump() const;
-
- ///////////////////// SETTERS ////////////////////////
- // Incident particle properties (before interactions in the target)
- // Vertex of interaction
- void SetICPositionX(Double_t PositionX) {fIC_Position_X.push_back(PositionX);}
- void SetICPositionY(Double_t PositionY) {fIC_Position_Y.push_back(PositionY);}
- void SetICPositionZ(Double_t PositionZ) {fIC_Position_Z.push_back(PositionZ);}
- // Theta and Phi angles for the emittance
- void SetICIncidentEmittanceTheta(Double_t Theta) {fIC_Incident_Emittance_Theta.push_back(Theta);}
- void SetICIncidentEmittancePhi(Double_t Phi) {fIC_Incident_Emittance_Phi.push_back(Phi);}
- // Incident particle angles
- void SetICIncidentAngleTheta(Double_t AngleTheta) {fIC_Incident_Angle_Theta.push_back(AngleTheta);}
- void SetICIncidentAnglePhi(Double_t AnglePhi) {fIC_Incident_Angle_Phi.push_back(AnglePhi);}
- // Incident particle energy
- void SetICIncidentEnergy(Double_t Energy) {fIC_Incident_Energy.push_back(Energy);}
-
- // Emitted particle angles
- // Center of mass
- void SetICEmittedAngleThetaCM(Double_t AngleTheta) {fIC_Emitted_Angle_ThetaCM.push_back(AngleTheta);}
- // Angles in the incident frame
- void SetICEmittedAngleThetaLabIncidentFrame(Double_t AngleTheta) {fIC_Emitted_Angle_ThetaLab_IncidentFrame.push_back(AngleTheta);}
- void SetICEmittedAnglePhiIncidentFrame(Double_t AnglePhi) {fIC_Emitted_Angle_Phi_IncidentFrame.push_back(AnglePhi);}
- // Angles in the world frame
- void SetICEmittedAngleThetaLabWorldFrame(Double_t AngleTheta) {fIC_Emitted_Angle_ThetaLab_WorldFrame.push_back(AngleTheta);}
- void SetICEmittedAnglePhiWorldFrame(Double_t AnglePhi) {fIC_Emitted_Angle_Phi_WorldFrame.push_back(AnglePhi);}
- // Emitted particle energy
- void SetICEmittedEnergy(Double_t Energy) {fIC_Emitted_Energy.push_back(Energy);}
-
-
- ///////////////////// GETTERS ////////////////////////
- // Incident particle properties (before interactions in the target)
- // Vertex of interaction
- Double_t GetICPositionX(Int_t i) {return fIC_Position_X.at(i);}
- Double_t GetICPositionY(Int_t i) {return fIC_Position_Y.at(i);}
- Double_t GetICPositionZ(Int_t i) {return fIC_Position_Z.at(i);}
- // Theta and Phi angles for the emittance
- Double_t GetICIncidentEmittanceTheta(Int_t i) {return fIC_Incident_Emittance_Theta.at(i);}
- Double_t GetICIncidentEmittancePhi(Int_t i) {return fIC_Incident_Emittance_Phi.at(i);}
- // Incident particle angles
- Double_t GetICIncidentAngleTheta(Int_t i) {return fIC_Incident_Angle_Theta.at(i);}
- Double_t GetICIncidentAnglePhi(Int_t i) {return fIC_Incident_Angle_Phi.at(i);}
- // Incident particle energy
- Double_t GetICIncidentEnergy(Int_t i) {return fIC_Incident_Energy.at(i);}
-
- // Emitted particle angles
- // Center of Mass
- Double_t GetICEmittedAngleThetaCM(Int_t i) {return fIC_Emitted_Angle_ThetaCM.at(i);}
- // Angles in the incident frame
- Double_t GetICEmittedAngleThetaLabIncidentFrame(Int_t i) {return fIC_Emitted_Angle_ThetaLab_IncidentFrame.at(i);}
- Double_t GetICEmittedAnglePhiIncidentFrame(Int_t i) {return fIC_Emitted_Angle_Phi_IncidentFrame.at(i);}
- // Angles in the world frame
- Double_t GetICEmittedAngleThetaLabWorldFrame(Int_t i) {return fIC_Emitted_Angle_ThetaLab_WorldFrame.at(i);}
- Double_t GetICEmittedAnglePhiWorldFrame(Int_t i) {return fIC_Emitted_Angle_Phi_WorldFrame.at(i);}
- // Emitted particle energy
- Double_t GetICEmittedEnergy(Int_t i) {return fIC_Emitted_Energy.at(i);}
-
- ClassDef(TInitialConditions, 1) // InitialConditions structure
+ TInitialConditions();
+ virtual ~TInitialConditions();
+
+ void Clear();
+ void Clear(const Option_t*) {};
+ void Dump() const;
+
+ ///////////////////// GetTERS ////////////////////////
+ // Incident beam parameter
+ void SetIncidentParticleName (string Incident_Particle_Name) {fIC_Incident_Particle_Name = Incident_Particle_Name;}
+ void SetIncidentKineticEnergy (double Incident_Kinetic_Energy) {fIC_Incident_Kinetic_Energy = Incident_Kinetic_Energy;}
+ void SetIncidentEmittanceTheta (double Incident_Emittance_Theta) {fIC_Incident_Emittance_Theta = Incident_Emittance_Theta;}
+ void SetIncidentEmittancePhi (double Incident_Emittance_Phi) {fIC_Incident_Emittance_Phi = Incident_Emittance_Phi;}
+
+ // Beam status at the initial interaction point
+ void SetInteractionKineticEnergy (double Interaction_Kinetic_Energy) {fIC_Interaction_Kinetic_Energy = Interaction_Kinetic_Energy;}
+ void SetInteractionPositionX (double Interaction_Position_X) {fIC_Interaction_Position_X = Interaction_Position_X;}
+ void SetInteractionPositionY (double Interaction_Position_Y) {fIC_Interaction_Position_Y = Interaction_Position_Y;}
+ void SetInteractionPositionZ (double Interaction_Position_Z) {fIC_Interaction_Position_Z = Interaction_Position_Z;}
+
+ // emmitted particle
+ void SetParticleName (string Particle_Name) {fIC_Particle_Name.push_back(Particle_Name);}
+ void SetProcessName (string Process_Name) {fIC_Process_Name.push_back(Process_Name);}
+ void SetThetaCM (double ThetaCM) {fIC_ThetaCM.push_back(ThetaCM);}
+ void SetKineticEnergy (double Kinetic_Energy) {fIC_Kinetic_Energy.push_back(Kinetic_Energy);}
+ void SetMomentumDirectionX (double Momentum_Direction_X) {fIC_Momentum_Direction_X.push_back(Momentum_Direction_X);}
+ void SetMomentumDirectionY (double Momentum_Direction_Y) {fIC_Momentum_Direction_Y.push_back(Momentum_Direction_Y);}
+ void SetMomentumDirectionZ (double Momentum_Direction_Z) {fIC_Momentum_Direction_Z.push_back(Momentum_Direction_Z);}
+ ///////////////////// GETTERS ////////////////////////
+ // Incident beam parameter
+ string GetIncidentParticleName () const {return fIC_Incident_Particle_Name ;}
+ double GetIncidentKineticEnergy () const {return fIC_Incident_Kinetic_Energy ;}
+ double GetIncidentEmittanceTheta () const {return fIC_Incident_Emittance_Theta ;}
+ double GetIncidentEmittancePhi () const {return fIC_Incident_Emittance_Phi ;}
+
+ // Beam status at the initial interaction point
+ double GetInteractionKineticEnergy () const {return fIC_Interaction_Kinetic_Energy ;}
+ double GetInteractionPositionX () const {return fIC_Interaction_Position_X ;}
+ double GetInteractionPositionY () const {return fIC_Interaction_Position_Y ;}
+ double GetInteractionPositionZ () const {return fIC_Interaction_Position_Z ;}
+
+ // emmitted particle
+ string GetParticleName (int i) const {return fIC_Particle_Name[i];}
+ string GetProcessName (int i) const {return fIC_Process_Name[i];}
+ double GetThetaCM (int i) const {return fIC_ThetaCM[i];}
+ double GetKineticEnergy (int i) const {return fIC_Kinetic_Energy[i];}
+ double GetMomentumDirectionX (int i) const {return fIC_Momentum_Direction_X[i];}
+ double GetMomentumDirectionY (int i) const {return fIC_Momentum_Direction_Y[i];}
+ double GetMomentumDirectionZ (int i) const {return fIC_Momentum_Direction_Z[i];}
+
+ unsigned int GetMult() const {return fIC_Particle_Name.size();}
+
+ ClassDef(TInitialConditions, 1) // InitialConditions structure
};
#endif
diff --git a/NPLib/Physics/NPBeam.cxx b/NPLib/Physics/NPBeam.cxx
index 56b62ced1..b4d8380f3 100644
--- a/NPLib/Physics/NPBeam.cxx
+++ b/NPLib/Physics/NPBeam.cxx
@@ -35,6 +35,9 @@
#include "NPFunction.h"
#include "NPOptionManager.h"
+// ROOT Header
+#include "TDirectory.h"
+
// Use CLHEP System of unit and Physical Constant
#include "CLHEP/Units/GlobalSystemOfUnits.h"
#include "CLHEP/Units/PhysicalConstants.h"
@@ -61,11 +64,16 @@ Beam::Beam(){
fTargetAngle = 0 ;
fTargetZ = 0 ;
fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
+
// case of user given distribution
- Global_BeamHistOffset++;
- fEnergyHist = new TH1F(Form("EnergyHist%i",Global_BeamHistOffset),"EnergyHist",1,0,1);
- fXThetaXHist = new TH2F(Form("XThetaXHis%i",Global_BeamHistOffset),"XThetaXHis",1,0,1,1,0,1);
- fYPhiYHist = new TH2F(Form("YPhiYHist%i",Global_BeamHistOffset),"YPhiYHist",1,0,1,1,0,1);
+ // do that to avoid warning from multiple Hist with same name...
+ int offset = 0;
+ while(gDirectory->FindObjectAny(Form("EnergyHist_%i",offset))!=0)
+ ++offset;
+
+ fEnergyHist = new TH1F(Form("EnergyHist_%i",offset),"EnergyHist",1,0,1);
+ fXThetaXHist = new TH2F(Form("XThetaXHis_%i",offset),"XThetaXHis",1,0,1,1,0,1);
+ fYPhiYHist = new TH2F(Form("YPhiYHist_%i",offset),"YPhiYHist",1,0,1,1,0,1);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -209,7 +217,7 @@ void Beam::ReadConfigurationFile(string Path){
string FileName,HistName;
ReactionFile >> FileName >> HistName;
if(fVerboseLevel==1) cout << "Reading Energy profile file: " << FileName << endl;
- fEnergyHist = Read1DProfile(FileName, HistName );
+ SetEnergyHist( Read1DProfile(FileName, HistName ));
}
else if (DataBuffer == "XThetaXProfilePath=") {
@@ -217,7 +225,7 @@ void Beam::ReadConfigurationFile(string Path){
string FileName,HistName;
ReactionFile >> FileName >> HistName;
if(fVerboseLevel==1) cout << "Reading X-ThetaX profile file: " << FileName << endl;
- fXThetaXHist = Read2DProfile(FileName, HistName );
+ SetXThetaXHist(Read2DProfile(FileName, HistName ) );
}
else if (DataBuffer == "YPhiYProfilePath=") {
@@ -225,7 +233,7 @@ void Beam::ReadConfigurationFile(string Path){
string FileName,HistName;
ReactionFile >> FileName >> HistName;
if(fVerboseLevel==1) cout << "Reading Y-ThetaY profile file: " << FileName << endl;
- fYPhiYHist = Read2DProfile(FileName, HistName );
+ SetYPhiYHist( Read2DProfile(FileName, HistName ));
}
diff --git a/NPLib/Physics/NPBeam.h b/NPLib/Physics/NPBeam.h
index 0e09f2cc1..63f4c4f4f 100755
--- a/NPLib/Physics/NPBeam.h
+++ b/NPLib/Physics/NPBeam.h
@@ -40,8 +40,6 @@ using namespace std;
using namespace NPL;
namespace NPL{
- // Needed to avoid warnig from multiple Th1 with same name
- static int Global_BeamHistOffset=0;
class Beam{
@@ -87,10 +85,10 @@ namespace NPL{
void SetMeanPhiY (double MeanPhiY) {fMeanPhiY=MeanPhiY;}
void SetSigmaThetaX (double SigmaThetaX) {fSigmaThetaX=SigmaThetaX;}
void SetSigmaPhiY (double SigmaPhiY) {fSigmaPhiY=SigmaPhiY;}
- void SetEnergyHist (TH1F* EnergyHist) {delete fEnergyHist; fEnergyHist = new TH1F(*EnergyHist);}
- void SetXThetaXHist (TH2F* XThetaXHist) {delete fXThetaXHist; fXThetaXHist = new TH2F(*XThetaXHist);}
- void SetYPhiYHist (TH2F* YPhiYHist) {delete fYPhiYHist; fYPhiYHist = new TH2F(*YPhiYHist);}
- void SetVerboseLevel(int verbose) {fVerboseLevel = verbose;}
+ void SetEnergyHist (TH1F* EnergyHist) {delete fEnergyHist; fEnergyHist = EnergyHist;}
+ void SetXThetaXHist (TH2F* XThetaXHist) {delete fXThetaXHist; fXThetaXHist = XThetaXHist;}
+ void SetYPhiYHist (TH2F* YPhiYHist) {delete fYPhiYHist; fYPhiYHist = YPhiYHist;}
+ void SetVerboseLevel(int verbose) {fVerboseLevel = verbose;}
// Get
Nucleus* GetNucleus () const {return fBeamNucleus;}
diff --git a/NPLib/Physics/NPFunction.cxx b/NPLib/Physics/NPFunction.cxx
index 9e3fc528a..fa29dde0c 100644
--- a/NPLib/Physics/NPFunction.cxx
+++ b/NPLib/Physics/NPFunction.cxx
@@ -58,28 +58,21 @@ TH1F* Read1DProfile(string filename,string HistName){
}
// Look for the step size, min and max of the distribution
- double min = 0 ;
- double max = 0 ;
- double previous = 0 ;
- int step =0;
+ double min ;
+ double max ;
unsigned int size = x.size();
if(size > 0){
min = x[0] ;
max = x[0] ;
- previous = x[0] ;
}
for(unsigned int i = 0 ; i < size ; i++){
if(x[i] > max) max = x[i] ;
if(x[i] < min) min = x[i] ;
-
- step= fabs(previous - x[i]) ;
- previous = x[i] ;
}
- h = new TH1F(HistName.c_str(),HistName.c_str(),step,min,max);
-
+ h = new TH1F(HistName.c_str(),HistName.c_str(),size,min,max);
for(unsigned int i = 0 ; i < size ; i++){
h->Fill(x[i],w[i]);
}
@@ -133,47 +126,35 @@ TH2F* Read2DProfile(string filename,string HistName){
}
// Look for the step size, min and max of the distribution
- double xmin = 0 ;
- double xmax = 0 ;
- double xprevious = 0 ;
- int xstep =0;
+ double xmin ;
+ double xmax ;
unsigned int xsize = x.size();
- double ymin = 0 ;
- double ymax = 0 ;
- double yprevious = 0 ;
- int ystep =0;
+ double ymin ;
+ double ymax ;
unsigned int ysize = y.size();
if(xsize > 0){
xmin = x[0] ;
xmax = x[0] ;
- xprevious = x[0] ;
}
if(ysize > 0){
ymin = y[0] ;
ymax = y[0] ;
- yprevious = y[0] ;
}
for(unsigned int i = 0 ; i < xsize ; i++){
if(x[i] > xmax) xmax = x[i] ;
if(x[i] < xmin) xmin = x[i] ;
-
- xstep= fabs(xprevious - x[i]) ;
- xprevious = x[i] ;
}
for(unsigned int i = 0 ; i < ysize ; i++){
if(y[i] > ymax) ymax = y[i] ;
if(y[i] < ymin) ymin = y[i] ;
-
- ystep= fabs(yprevious - y[i]) ;
- yprevious = y[i] ;
}
- h = new TH2F(HistName.c_str(),HistName.c_str(),xstep,xmin,xmax,ystep,ymin,ymax);
+ h = new TH2F(HistName.c_str(),HistName.c_str(),xsize,xmin,xmax,ysize,ymin,ymax);
for(unsigned int i = 0 ; i < xsize ; i++){
h->Fill(x[i],y[i],w[i]);
diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index fc7ba5379..e7bcce344 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -69,9 +69,13 @@ Reaction::Reaction(){
fQValue = 0;
fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
initializePrecomputeVariable();
- // Needed to avoid warning at compilation, not very clean...
- Global_BeamHistOffset=Global_BeamHistOffset;
- fCrossSectionHist = new TH1F(Form("Reaction_CS%i",Global_ReactionHistOffset++),"Reaction_CS",180,0,180);
+
+ // do that to avoid warning from multiple Hist with same name... int offset = 0;
+ int offset = 0;
+ while(gDirectory->FindObjectAny(Form("EnergyHist_%i",offset))!=0)
+ ++offset;
+
+ fCrossSectionHist = new TH1F(Form("EnergyHist_%i",offset),"Reaction_CS",1,0,180);
fshoot3=true;
fshoot4=true;
}
@@ -84,6 +88,7 @@ Reaction::~Reaction(){
delete fNuclei2;
delete fNuclei3;
delete fNuclei4;
+ delete fCrossSectionHist;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -181,7 +186,7 @@ double Reaction::EnergyLabToThetaCM(double EnergyLab, double ThetaLab){
fEnergyImpulsionCM_3 = fEnergyImpulsionLab_3;
fEnergyImpulsionCM_3.Boost(0,0,-BetaCM);
- double ThetaCM = CLHEP::pi - fEnergyImpulsionCM_1.Angle(fEnergyImpulsionCM_3.Vect());
+ double ThetaCM = M_PI - fEnergyImpulsionCM_1.Angle(fEnergyImpulsionCM_3.Vect());
return(ThetaCM);
}
@@ -300,7 +305,8 @@ void Reaction::ReadConfigurationFile(string Path){
string FileName,HistName;
ReactionFile >> FileName >> HistName;
if(fVerboseLevel==1) cout << "Reading Cross Section file: " << FileName << endl;
- fCrossSectionHist = Read1DProfile(FileName, HistName );
+ SetCrossSectionHist( Read1DProfile(FileName, HistName ));
+ cout << "cccc " << fCrossSectionHist->GetNbinsX() << endl;
}
else if (DataBuffer.compare(0, 17, "HalfOpenAngleMin=") == 0) {
@@ -352,7 +358,7 @@ void Reaction::ReadConfigurationFile(string Path){
// Pick up the beam energy from the Beam event generator
NPL::Beam* localBeam= new NPL::Beam();
- // localBeam->SetVerboseLevel(0);
+ localBeam->SetVerboseLevel(0);
localBeam->ReadConfigurationFile(Path);
// Modifiy the CS to shoot only within ]HalfOpenAngleMin,HalfOpenAngleMax[
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index 343c7d5e1..ee58d18b7 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -50,8 +50,6 @@
using namespace std;
namespace NPL{
- // Needed to avoid warnig from multiple Th1 with same name
- static int Global_ReactionHistOffset=0;
class Reaction{
@@ -91,6 +89,7 @@ namespace NPL{
void SetExcitationLight(double exci) {fExcitation3 = exci; initializePrecomputeVariable();}
void SetExcitationHeavy(double exci) {fExcitation4 = exci; initializePrecomputeVariable();}
void SetVerboseLevel(int verbose) {fVerboseLevel = verbose;}
+ void SetCrossSectionHist (TH1F* CrossSectionHist) {delete fCrossSectionHist; fCrossSectionHist = CrossSectionHist;}
double GetBeamEnergy() const {return fBeamEnergy;}
double GetThetaCM() const {return fThetaCM;}
diff --git a/NPSimulation/Simulation.cc b/NPSimulation/Simulation.cc
index 6306d9fe9..830f2c6de 100644
--- a/NPSimulation/Simulation.cc
+++ b/NPSimulation/Simulation.cc
@@ -131,9 +131,10 @@ int main(int argc, char** argv)
///////////////////////////////////////////////////////////////
////////////////////// Job termination ////////////////////////
///////////////////////////////////////////////////////////////
- RootOutput::getInstance()->Destroy();
- NPOptionManager::getInstance()->Destroy();
+ // delete primary; delete detector;
delete runManager;
+ NPOptionManager::getInstance()->Destroy();
+ RootOutput::getInstance()->Destroy();
return 0;
}
diff --git a/NPSimulation/include/EventGeneratorIsotropic.hh b/NPSimulation/include/EventGeneratorIsotropic.hh
index f6653d17a..f370a3ec4 100644
--- a/NPSimulation/include/EventGeneratorIsotropic.hh
+++ b/NPSimulation/include/EventGeneratorIsotropic.hh
@@ -55,8 +55,5 @@ private: // Source parameter from input file
G4double m_y0 ; // Vertex Position Y
G4double m_z0 ; // Vertex Position Z
G4ParticleDefinition* m_particle ; // Kind of particle to shoot isotropically
-
- private: // Output ROOT class storing characteristics of source
- TInitialConditions* m_InitConditions;
};
#endif
diff --git a/NPSimulation/include/Particle.hh b/NPSimulation/include/Particle.hh
index 7be0c5bf9..d3eef1924 100644
--- a/NPSimulation/include/Particle.hh
+++ b/NPSimulation/include/Particle.hh
@@ -28,33 +28,36 @@
#include"G4ThreeVector.hh"
class Particle{
-
+
public: // Constructor and Destructor
- // Empty constructor (return geantino at zero degree)
- Particle();
- ~Particle();
- // Constructor to be used
- Particle(G4ParticleDefinition* particle,double T,G4ThreeVector Direction, G4ThreeVector Position,bool ShootStatus=true);
-
-private: // Private Member
- G4ParticleDefinition* m_ParticleDefinition;
- double m_T ;
- G4ThreeVector m_Direction;
- G4ThreeVector m_Position;
- bool m_ShootStatus;
-
+ // Empty constructor (return geantino at zero degree)
+ Particle();
+ ~Particle();
+ // Constructor to be used
+ Particle(G4ParticleDefinition* particle,double ThetaCM,double T,G4ThreeVector Direction, G4ThreeVector Position,bool ShootStatus=true);
+
+private: // Private Member
+ G4ParticleDefinition* m_ParticleDefinition;
+ double m_ThetaCM;
+ double m_T ;
+ G4ThreeVector m_Direction;
+ G4ThreeVector m_Position;
+ bool m_ShootStatus;
+
public: // Setter and Getter
- // Getter
- G4ParticleDefinition* GetParticleDefinition();
- double GetParticleKineticEnergy();
- G4ThreeVector GetParticleMomentumDirection();
- G4ThreeVector GetParticlePosition();
- bool GetShootStatus();
- // Setter
- void SetParticleDefinition(G4ParticleDefinition*);
- void SetParticleKineticEnergy(double);
- void SetParticlePosition(G4ThreeVector);
- void SetParticleMomentumDirection(G4ThreeVector);
- void SetShootStatus(bool);
+ // Getter
+ G4ParticleDefinition* GetParticleDefinition();
+ double GetParticleThetaCM();
+ double GetParticleKineticEnergy();
+ G4ThreeVector GetParticleMomentumDirection();
+ G4ThreeVector GetParticlePosition();
+ bool GetShootStatus();
+ // Setter
+ void SetParticleDefinition(G4ParticleDefinition*);
+ void SetParticleThetaCM(double);
+ void SetParticleKineticEnergy(double);
+ void SetParticlePosition(G4ThreeVector);
+ void SetParticleMomentumDirection(G4ThreeVector);
+ void SetShootStatus(bool);
};
#endif
\ No newline at end of file
diff --git a/NPSimulation/include/ParticleStack.hh b/NPSimulation/include/ParticleStack.hh
index 5dabb7492..b76cdba5b 100644
--- a/NPSimulation/include/ParticleStack.hh
+++ b/NPSimulation/include/ParticleStack.hh
@@ -79,8 +79,11 @@ public: // Getter and Setter
public: // Particle management and shooting method
// EventGenerator use this method to add particle in the stack
- void AddParticleToStack(Particle);
-
+ void AddParticleToStack(Particle&);
+
+ // Add the Particle to the stack and fill Initial Condition Beam field
+ void AddBeamParticleToStack(Particle&);
+
// Search for a specific particle in the stack
Particle SearchAndRemoveParticle(string);
diff --git a/NPSimulation/src/EventGeneratorBeam.cc b/NPSimulation/src/EventGeneratorBeam.cc
index ffc9c4991..4a8b3f36d 100644
--- a/NPSimulation/src/EventGeneratorBeam.cc
+++ b/NPSimulation/src/EventGeneratorBeam.cc
@@ -52,6 +52,7 @@ EventGeneratorBeam::EventGeneratorBeam(){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorBeam::~EventGeneratorBeam(){
delete m_InitConditions;
+ delete m_Beam;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -99,21 +100,22 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent){
G4double Xdir = sin(Beam_thetaX);
G4double Ydir = sin(Beam_phiY);
G4double Zdir = cos(Beam_thetaX) + cos(Beam_phiY);
- G4ThreeVector BeamDir = G4ThreeVector(Xdir,Ydir,Zdir) ;
-
+ G4ThreeVector BeamDir(Xdir,Ydir,Zdir);
+ G4ThreeVector BeamPos(x0,y0,z0);
Beam_theta = BeamDir.theta() ;
Beam_phi = BeamDir.phi() ;
FinalBeamEnergy = m_Target->SlowDownBeam(m_particle, InitialBeamEnergy,z0,Beam_theta);
///////////////////////////////////////////////////////
///// Add the Beam particle to the particle Stack /////
///////////////////////////////////////////////////////
-
Particle BeamParticle( m_particle,
- FinalBeamEnergy,
- BeamDir.unit(),
- G4ThreeVector(x0,y0,z0),
- 1);
- m_ParticleStack->AddParticleToStack(BeamParticle);
+ InitialBeamEnergy,
+ FinalBeamEnergy,
+ BeamDir.unit(),
+ BeamPos,
+ 1);
+
+ m_ParticleStack->AddBeamParticleToStack(BeamParticle);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPSimulation/src/EventGeneratorGammaDecay.cc b/NPSimulation/src/EventGeneratorGammaDecay.cc
index bb34fe49a..d36b73476 100644
--- a/NPSimulation/src/EventGeneratorGammaDecay.cc
+++ b/NPSimulation/src/EventGeneratorGammaDecay.cc
@@ -236,6 +236,7 @@ void EventGeneratorGammaDecay::GenerateEvent(G4Event*){
= G4ParticleTable::GetParticleTable()->GetIon(decayingParticle.GetParticleDefinition()->GetAtomicNumber(), decayingParticle.GetParticleDefinition()->GetAtomicMass(), FinalExcitationEnergy*MeV);
Particle FinalParticle = Particle( FinalParticleDefition,
+ decayingParticle.GetParticleThetaCM(),
decayingParticle.GetParticleKineticEnergy(),
decayingParticle.GetParticleMomentumDirection(),
decayingParticle.GetParticlePosition(),
@@ -294,7 +295,7 @@ void EventGeneratorGammaDecay::GenerateEvent(G4Event*){
GammaLV.Py(),
GammaLV.Pz());
// Add the gamma to the stack
- Particle gammaParticle(gammaDefinition,GammaLV.E(),gammaDirection, decayingParticle.GetParticlePosition());
+ Particle gammaParticle(gammaDefinition,theta,GammaLV.E(),gammaDirection, decayingParticle.GetParticlePosition());
m_ParticleStack->AddParticleToStack(gammaParticle);
}
diff --git a/NPSimulation/src/EventGeneratorIsotropic.cc b/NPSimulation/src/EventGeneratorIsotropic.cc
index 3771fb438..a1f1d714a 100644
--- a/NPSimulation/src/EventGeneratorIsotropic.cc
+++ b/NPSimulation/src/EventGeneratorIsotropic.cc
@@ -40,7 +40,6 @@ using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorIsotropic::EventGeneratorIsotropic()
{
- m_InitConditions = new TInitialConditions();
m_EnergyLow = 0 ;
m_EnergyHigh = 0 ;
@@ -55,7 +54,6 @@ EventGeneratorIsotropic::EventGeneratorIsotropic()
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
EventGeneratorIsotropic::~EventGeneratorIsotropic()
{
- delete m_InitConditions;
}
@@ -187,7 +185,6 @@ void EventGeneratorIsotropic::ReadConfiguration(string Path)
void EventGeneratorIsotropic::GenerateEvent(G4Event* anEvent, G4ParticleGun* particleGun)
{
// Clear TInitialConditions
- m_InitConditions->Clear();
particleGun->SetParticleDefinition(m_particle);
@@ -209,28 +206,6 @@ void EventGeneratorIsotropic::GenerateEvent(G4Event* anEvent, G4ParticleGun* par
particleGun->SetParticleEnergy(particle_energy) ;
particleGun->SetParticlePosition(G4ThreeVector(m_x0, m_y0, m_z0)) ;
- // Fill TInitialConditions class
- // Interaction vertex
- m_InitConditions->SetICPositionX(m_x0 / mm);
- m_InitConditions->SetICPositionY(m_y0 / mm);
- m_InitConditions->SetICPositionZ(m_z0 / mm);
- // Incident "particles"
- // Everything is zero for a source
- m_InitConditions->SetICIncidentEmittanceTheta(0);
- m_InitConditions->SetICIncidentEmittancePhi(0);
- m_InitConditions->SetICIncidentAngleTheta(0);
- m_InitConditions->SetICIncidentAnglePhi(0);
- m_InitConditions->SetICIncidentEnergy(0);
- // Emitted particle angles
- m_InitConditions->SetICEmittedAngleThetaCM(theta / deg);
- m_InitConditions->SetICEmittedAngleThetaLabIncidentFrame(theta / deg);
- m_InitConditions->SetICEmittedAnglePhiIncidentFrame(phi / deg);
- m_InitConditions->SetICEmittedAngleThetaLabWorldFrame(theta / deg);
- m_InitConditions->SetICEmittedAnglePhiWorldFrame(phi / deg);
- // Emitted particle energy
- m_InitConditions->SetICEmittedEnergy(particle_energy / MeV);
-
-
//Shoot particle
particleGun->GeneratePrimaryVertex(anEvent) ;
}
@@ -240,7 +215,5 @@ void EventGeneratorIsotropic::GenerateEvent(G4Event* anEvent, G4ParticleGun* par
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void EventGeneratorIsotropic::InitializeRootOutput()
{
- RootOutput *pAnalysis = RootOutput::getInstance();
- TTree *pTree = pAnalysis->GetTree();
- pTree->Branch("InitialConditions", "TInitialConditions", &m_InitConditions);
+
}
diff --git a/NPSimulation/src/EventGeneratorParticleDecay.cc b/NPSimulation/src/EventGeneratorParticleDecay.cc
index 336058d2e..0ae6f688d 100644
--- a/NPSimulation/src/EventGeneratorParticleDecay.cc
+++ b/NPSimulation/src/EventGeneratorParticleDecay.cc
@@ -247,6 +247,7 @@ void EventGeneratorParticleDecay::GenerateEvent(G4Event* anEvent){
FirstDaughterLV.Z() );
Particle FirstDaughterParticle( m_DaughterNuclei[0],
+ ThetaCM,
FirstDaughterLV.E()-m_DaughterNuclei[0]->GetPDGMass(),
DaughterDirection.unit(),
decayingParticle.GetParticlePosition(),
@@ -257,6 +258,7 @@ void EventGeneratorParticleDecay::GenerateEvent(G4Event* anEvent){
SecondDaughterLV.Z() );
Particle SecondDaughterParticle( m_DaughterNuclei[1],
+ ThetaCM+M_PI,
SecondDaughterLV.E()-m_DaughterNuclei[1]->GetPDGMass(),
DaughterDirection.unit(),
decayingParticle.GetParticlePosition(),
@@ -288,12 +290,13 @@ void EventGeneratorParticleDecay::GenerateEvent(G4Event* anEvent){
daughterLV = m_TPhaseSpace.GetDecay(i);
G4ThreeVector daughterDirection = G4ThreeVector( daughterLV->X() ,
- daughterLV->Y() ,
- daughterLV->Z() );
+ daughterLV->Y() ,
+ daughterLV->Z() );
KineticEnergy = daughterLV->E()-m_Masses[i] ;
Particle daughterParticle( m_DaughterNuclei[i],
+ -1,
KineticEnergy*GeV,
daughterDirection.unit(),
decayingParticle.GetParticlePosition(),
diff --git a/NPSimulation/src/EventGeneratorTwoBodyReaction.cc b/NPSimulation/src/EventGeneratorTwoBodyReaction.cc
index c71dd4907..2136cc88f 100644
--- a/NPSimulation/src/EventGeneratorTwoBodyReaction.cc
+++ b/NPSimulation/src/EventGeneratorTwoBodyReaction.cc
@@ -186,44 +186,20 @@ void EventGeneratorTwoBodyReaction::GenerateEvent(G4Event* anEvent){
///////// Set up everything for shooting /////////
//////////////////////////////////////////////////
// Case of light particle
- // Instantiate a new particle
- Particle LightParticle;
-
- // Particle type
- LightParticle.SetParticleDefinition(LightName);
- // Particle energy
- LightParticle.SetParticleKineticEnergy(EnergyLight);
- // Particle vertex position
- LightParticle.SetParticlePosition(BeamParticle.GetParticlePosition());
// Particle direction
// Kinematical angles in the beam frame are transformed
// to the "world" frame
G4ThreeVector momentum_kine_world = BeamToWorld * momentum_kineLight_beam;
- //Set the Momentum Direction
- LightParticle.SetParticleMomentumDirection(momentum_kine_world);
- // Set the shoot status
- LightParticle.SetShootStatus(m_ShootLight) ;
//Add the particle to the particle stack
+ Particle LightParticle(LightName, ThetaCM , EnergyLight,momentum_kine_world, BeamParticle.GetParticlePosition(), m_ShootLight);
m_ParticleStack->AddParticleToStack(LightParticle);
// Case of heavy particle
- // Instantiate a new particle
- Particle HeavyParticle;
-
- // Particle type
- HeavyParticle.SetParticleDefinition(HeavyName);
- // Particle energy
- HeavyParticle.SetParticleKineticEnergy(EnergyHeavy);
- // Particle vertex position
- HeavyParticle.SetParticlePosition(BeamParticle.GetParticlePosition());
// Particle direction
// Kinematical angles in the beam frame are transformed
// to the "world" frame
momentum_kine_world = BeamToWorld * momentum_kineHeavy_beam;
- //Set the Momentum Direction
- HeavyParticle.SetParticleMomentumDirection(momentum_kine_world);
- // Set the shoot status
- HeavyParticle.SetShootStatus(m_ShootHeavy) ;
+ Particle HeavyParticle(HeavyName, ThetaCM + M_PI, EnergyHeavy,momentum_kine_world, BeamParticle.GetParticlePosition(), m_ShootHeavy);
//Add the particle to the particle stack
m_ParticleStack->AddParticleToStack(HeavyParticle);
}
diff --git a/NPSimulation/src/Particle.cc b/NPSimulation/src/Particle.cc
index 4cdb27810..4feffe6ef 100644
--- a/NPSimulation/src/Particle.cc
+++ b/NPSimulation/src/Particle.cc
@@ -23,60 +23,83 @@
#include "Particle.hh"
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Particle::Particle(){
- m_ParticleDefinition = NULL;
- m_T = 1*GeV;
- m_Direction = G4ThreeVector(0,0,1);
+ m_ParticleDefinition = NULL;
+ m_T = 1*GeV;
+ m_Direction = G4ThreeVector(0,0,1);
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Particle::~Particle(){
}
-Particle::Particle(G4ParticleDefinition* particle,double T,G4ThreeVector Direction, G4ThreeVector Position,bool ShootStatus){
- m_ParticleDefinition = particle;
- m_T = T;
- m_Direction = Direction;
- m_Position = Position;
- m_ShootStatus = ShootStatus;
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+Particle::Particle(G4ParticleDefinition* particle,double ThetaCM,double T,G4ThreeVector Direction, G4ThreeVector Position,bool ShootStatus){
+ m_ParticleDefinition = particle;
+ m_ThetaCM = ThetaCM;
+ m_T = T;
+ m_Direction = Direction;
+ m_Position = Position;
+ m_ShootStatus = ShootStatus;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4ParticleDefinition* Particle::GetParticleDefinition(){
- return m_ParticleDefinition;
+ return m_ParticleDefinition;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+double Particle::GetParticleThetaCM(){
+ return m_ThetaCM;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Particle::GetParticleKineticEnergy(){
- return m_T;
+ return m_T;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4ThreeVector Particle::GetParticleMomentumDirection(){
- return m_Direction;
+ return m_Direction;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4ThreeVector Particle::GetParticlePosition(){
- return m_Position;
+ return m_Position;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
bool Particle::GetShootStatus(){
- return m_ShootStatus;
+ return m_ShootStatus;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Particle::SetParticleDefinition(G4ParticleDefinition* particle){
- m_ParticleDefinition = particle;
+ m_ParticleDefinition = particle;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Particle::SetParticleThetaCM(double ThetaCM){
+ m_ThetaCM = ThetaCM ;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Particle::SetParticleKineticEnergy(double T){
- m_T = T ;
+ m_T = T ;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Particle::SetParticleMomentumDirection(G4ThreeVector Direction){
- m_Direction = Direction;
+ m_Direction = Direction;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Particle::SetParticlePosition(G4ThreeVector Position){
- m_Position = Position;
+ m_Position = Position;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Particle::SetShootStatus(bool ShootStatus){
- m_ShootStatus = ShootStatus;
-}
-
+ m_ShootStatus = ShootStatus;
+}
\ No newline at end of file
diff --git a/NPSimulation/src/ParticleStack.cc b/NPSimulation/src/ParticleStack.cc
index 9c72a59d1..67d3363a5 100644
--- a/NPSimulation/src/ParticleStack.cc
+++ b/NPSimulation/src/ParticleStack.cc
@@ -29,14 +29,17 @@
// NPL
#include "RootOutput.h"
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
ParticleStack* ParticleStack::instance = 0 ;
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
ParticleStack* ParticleStack::getInstance()
{
if (instance == 0) instance = new ParticleStack();
return instance ;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
ParticleStack::ParticleStack(){
m_particleGun = new G4ParticleGun(1);
@@ -56,37 +59,61 @@ ParticleStack::~ParticleStack(){
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
vector<Particle> ParticleStack::GetParticleStack(){
return m_ParticleStack;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ParticleStack::SetParticleStack(vector<Particle> particle_stack){
m_ParticleStack = particle_stack;
}
-void ParticleStack::AddParticleToStack(Particle particle){
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+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());
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void ParticleStack::AddBeamParticleToStack(Particle& particle){
+ m_InitialConditions->Clear();
+ m_ParticleStack.push_back(particle);
+
+ // Incident beam parameter
+ m_InitialConditions-> SetIncidentParticleName (particle.GetParticleDefinition()->GetParticleName());
+ m_InitialConditions-> SetIncidentKineticEnergy (particle. GetParticleThetaCM());
+
+ G4ThreeVector U(1,0,0);
+ G4ThreeVector V(0,1,0);
+
+ m_InitialConditions-> SetIncidentEmittanceTheta (particle.GetParticleMomentumDirection().angle(U)/deg);
+ m_InitialConditions-> SetIncidentEmittancePhi (particle.GetParticleMomentumDirection().angle(V)/deg);
+
+ // Beam status at the initial interaction point
+ m_InitialConditions-> SetInteractionKineticEnergy (particle. GetParticleKineticEnergy());
+ m_InitialConditions-> SetInteractionPositionX (particle. GetParticlePosition().x());
+ m_InitialConditions-> SetInteractionPositionY (particle. GetParticlePosition().y());
+ m_InitialConditions-> SetInteractionPositionZ (particle. GetParticlePosition().x());
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Particle ParticleStack::SearchAndRemoveParticle(string name){
- for(unsigned int i = 0 ; i < m_ParticleStack.size() ; i++){
+ unsigned int size = m_ParticleStack.size();
+
+ for(unsigned int i = 0 ; i < size ; i++){
string ParticleName = m_ParticleStack[i].GetParticleDefinition()->GetParticleName();
- if(ParticleName.compare(0, name.length(), name) == 0)
- {
+ if(ParticleName.compare(0, name.length(), name) == 0){
Particle my_Particule = m_ParticleStack[i];
m_ParticleStack.erase(m_ParticleStack.begin()+i);
return my_Particule;
- }
+ }
}
return Particle();
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
string ParticleStack::ChangeNameToG4Standard(string OriginalName){
string NumberOfMass ;
string Nucleid;
@@ -94,7 +121,7 @@ string ParticleStack::ChangeNameToG4Standard(string OriginalName){
for (unsigned int i = 0; i < OriginalName.length(); i++) {
ostringstream character;
character << OriginalName[i];
- if (character.str()=="0") NumberOfMass+="0";
+ if (character.str()=="0") NumberOfMass+="0";
else if (character.str()=="1") NumberOfMass+="1";
else if (character.str()=="2") NumberOfMass+="2";
else if (character.str()=="3") NumberOfMass+="3";
@@ -162,7 +189,7 @@ string ParticleStack::ChangeNameToG4Standard(string OriginalName){
// Special case for light particles
string FinalName=Nucleid+NumberOfMass;
- if (FinalName=="H1") FinalName="proton";
+ if (FinalName=="H1") FinalName="proton";
else if (FinalName=="H2") FinalName="deuteron";
else if (FinalName=="H3") FinalName="triton";
else if (FinalName=="He4") FinalName="alpha";
@@ -179,8 +206,11 @@ string ParticleStack::ChangeNameToG4Standard(string OriginalName){
return(FinalName);
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ParticleStack::ShootAllParticle(G4Event* anEvent){
- for(unsigned int i = 0 ; i < m_ParticleStack.size() ; i++){
+ unsigned int size = m_ParticleStack.size();
+
+ for(unsigned int i = 0 ; i < size ; i++){
if(m_ParticleStack[i].GetShootStatus()){
m_particleGun->SetParticleDefinition(m_ParticleStack[i].GetParticleDefinition());
@@ -188,8 +218,15 @@ void ParticleStack::ShootAllParticle(G4Event* anEvent){
m_particleGun->SetParticleMomentumDirection(m_ParticleStack[i].GetParticleMomentumDirection());
m_particleGun->SetParticlePosition(m_ParticleStack[i].GetParticlePosition());
m_particleGun->GeneratePrimaryVertex(anEvent);
+
+ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+ m_InitialConditions-> SetParticleName ( m_ParticleStack[i].GetParticleDefinition()->GetParticleName()) ;
+ m_InitialConditions-> SetThetaCM ( m_ParticleStack[i].GetParticleThetaCM()/deg);
+ m_InitialConditions-> SetKineticEnergy ( m_ParticleStack[i].GetParticleKineticEnergy());
+ m_InitialConditions-> SetMomentumDirectionX ( m_ParticleStack[i].GetParticleMomentumDirection().x());
+ m_InitialConditions-> SetMomentumDirectionY ( m_ParticleStack[i].GetParticleMomentumDirection().y());
+ m_InitialConditions-> SetMomentumDirectionZ ( m_ParticleStack[i].GetParticleMomentumDirection().z());
}
}
m_ParticleStack.clear();
- m_InitialConditions->Clear();
}
\ No newline at end of file
--
GitLab