From 378759115f4a6729144130e8f3ee60d60740df94 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Val=C3=A9rian=20Alcindor?= <valerian.alcindor@hotmail.fr>
Date: Thu, 9 Sep 2021 13:46:05 +0200
Subject: [PATCH] Several updates. 1) enable possiblity to record more
information from the tracks using Stepping Action. 2) Add to the p2p reaction
the possibility for the heavy nucleus to decay using the ABLA deexcitation
code included in G4.
---
NPLib/Physics/CMakeLists.txt | 7 +-
NPLib/Physics/NPFunction.cxx | 139 +++++++++++---------
NPLib/Physics/NPQFS.cxx | 5 +-
NPLib/Physics/NPQFS.h | 2 +
NPLib/Physics/TInteractionCoordinates.cxx | 8 ++
NPLib/Physics/TInteractionCoordinates.h | 62 ++++++++-
NPLib/Physics/TTrackInfo.cxx | 83 ++++++++++++
NPLib/Physics/TTrackInfo.h | 144 +++++++++++++++++++++
NPSimulation/Core/CMakeLists.txt | 2 +-
NPSimulation/Core/ParticleStack.cc | 2 +-
NPSimulation/Core/SteppingAction.cc | 144 +++++++++++++++++----
NPSimulation/Core/SteppingAction.hh | 59 ++++++---
NPSimulation/Process/BeamReaction.cc | 59 +++++++--
NPSimulation/Process/BeamReaction.hh | 4 +
NPSimulation/Scorers/InteractionScorers.cc | 40 ++++--
NPSimulation/Scorers/InteractionScorers.hh | 85 +++++++++++-
16 files changed, 706 insertions(+), 139 deletions(-)
create mode 100644 NPLib/Physics/TTrackInfo.cxx
create mode 100644 NPLib/Physics/TTrackInfo.h
diff --git a/NPLib/Physics/CMakeLists.txt b/NPLib/Physics/CMakeLists.txt
index 02cf69825..3a3487d33 100644
--- a/NPLib/Physics/CMakeLists.txt
+++ b/NPLib/Physics/CMakeLists.txt
@@ -14,6 +14,8 @@ add_custom_command(OUTPUT NPEnergyLossDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scrip
add_custom_command(OUTPUT TInitialConditionsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TInitialConditions.h TInitialConditionsDict.cxx TInitialConditions.rootmap libNPInitialConditions.so DEPENDS TInitialConditions.h)
+add_custom_command(OUTPUT TTrackInfoDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TTrackInfo.h TTrackInfoDict.cxx TTrackInfo.rootmap libNPTrackInfo.so DEPENDS TTrackInfo.h)
+
add_custom_command(OUTPUT TInteractionCoordinatesDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TInteractionCoordinates.h TInteractionCoordinatesDict.cxx TInteractionCoordinates.rootmap libNPInteractionCoordinates.so DEPENDS TInteractionCoordinates.h)
add_custom_command(OUTPUT TReactionConditionsDict.cxx COMMAND ${CMAKE_BINARY_DIR}/scripts/build_dict.sh TReactionConditions.h TReactionConditionsDict.cxx TReactionConditions.rootmap libNPReactionConditions.so DEPENDS TReactionConditions.h)
@@ -26,6 +28,9 @@ target_link_libraries(NPPhysics ${ROOT_LIBRARIES} Physics NPCore)
add_library(NPInitialConditions SHARED TInitialConditions.cxx TInitialConditionsDict.cxx )
target_link_libraries(NPInitialConditions ${ROOT_LIBRARIES} )
+add_library(NPTrackInfo SHARED TTrackInfo.cxx TTrackInfoDict.cxx )
+target_link_libraries(NPTrackInfo ${ROOT_LIBRARIES} )
+
add_library(NPInteractionCoordinates SHARED TInteractionCoordinates.cxx TInteractionCoordinatesDict.cxx)
target_link_libraries(NPInteractionCoordinates ${ROOT_LIBRARIES} )
@@ -36,4 +41,4 @@ add_library(NPFissionConditions SHARED TFissionConditions.cxx TFissionConditions
target_link_libraries(NPFissionConditions ${ROOT_LIBRARIES} )
-install(FILES GEF.h NPFissionDecay.h NPDecay.h NPBeam.h NPEnergyLoss.h NPFunction.h NPParticle.h NPNucleus.h NPReaction.h NPQFS.h TInitialConditions.h TInteractionCoordinates.h TReactionConditions.h TFissionConditions.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+install(FILES GEF.h NPFissionDecay.h NPDecay.h NPBeam.h NPEnergyLoss.h NPFunction.h NPParticle.h NPNucleus.h NPReaction.h NPQFS.h TInitialConditions.h TTrackInfo.h TInteractionCoordinates.h TReactionConditions.h TFissionConditions.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/Physics/NPFunction.cxx b/NPLib/Physics/NPFunction.cxx
index 9de2b6662..5d59b7e8e 100644
--- a/NPLib/Physics/NPFunction.cxx
+++ b/NPLib/Physics/NPFunction.cxx
@@ -374,73 +374,88 @@ namespace NPL{
string NumberOfMass ;
string Nucleid;
+ bool g4_exc = false;
for (unsigned int i = 0; i < OriginalName.length(); i++) {
ostringstream character;
character << OriginalName[i];
- 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";
- else if (character.str()=="4") NumberOfMass+="4";
- else if (character.str()=="5") NumberOfMass+="5";
- else if (character.str()=="6") NumberOfMass+="6";
- else if (character.str()=="7") NumberOfMass+="7";
- else if (character.str()=="8") NumberOfMass+="8";
- else if (character.str()=="9") NumberOfMass+="9";
+ // std::cout << OriginalName << " "<< OriginalName[i]<< std::endl;
- else if (character.str()=="A") Nucleid+="A";
- else if (character.str()=="B") Nucleid+="B";
- else if (character.str()=="C") Nucleid+="C";
- else if (character.str()=="D") Nucleid+="D";
- else if (character.str()=="E") Nucleid+="E";
- else if (character.str()=="F") Nucleid+="F";
- else if (character.str()=="G") Nucleid+="G";
- else if (character.str()=="H") Nucleid+="H";
- else if (character.str()=="I") Nucleid+="I";
- else if (character.str()=="J") Nucleid+="J";
- else if (character.str()=="K") Nucleid+="K";
- else if (character.str()=="L") Nucleid+="L";
- else if (character.str()=="M") Nucleid+="M";
- else if (character.str()=="N") Nucleid+="N";
- else if (character.str()=="O") Nucleid+="O";
- else if (character.str()=="P") Nucleid+="P";
- else if (character.str()=="Q") Nucleid+="Q";
- else if (character.str()=="R") Nucleid+="R";
- else if (character.str()=="S") Nucleid+="S";
- else if (character.str()=="T") Nucleid+="T";
- else if (character.str()=="U") Nucleid+="U";
- else if (character.str()=="V") Nucleid+="V";
- else if (character.str()=="W") Nucleid+="W";
- else if (character.str()=="X") Nucleid+="X";
- else if (character.str()=="Y") Nucleid+="Y";
- else if (character.str()=="Z") Nucleid+="Z";
+ // To remove Nuclei excitation from the name
+ if (character.str()=="["){
+ g4_exc=true;
+ }
+ if (character.str()=="]"){
+ g4_exc=false;
+ }
- else if (character.str()=="a") Nucleid+="a";
- else if (character.str()=="b") Nucleid+="b";
- else if (character.str()=="c") Nucleid+="c";
- else if (character.str()=="d") Nucleid+="d";
- else if (character.str()=="e") Nucleid+="e";
- else if (character.str()=="f") Nucleid+="f";
- else if (character.str()=="g") Nucleid+="g";
- else if (character.str()=="h") Nucleid+="h";
- else if (character.str()=="i") Nucleid+="i";
- else if (character.str()=="j") Nucleid+="j";
- else if (character.str()=="k") Nucleid+="k";
- else if (character.str()=="l") Nucleid+="l";
- else if (character.str()=="m") Nucleid+="m";
- else if (character.str()=="n") Nucleid+="n";
- else if (character.str()=="o") Nucleid+="o";
- else if (character.str()=="p") Nucleid+="p";
- else if (character.str()=="q") Nucleid+="q";
- else if (character.str()=="r") Nucleid+="r";
- else if (character.str()=="s") Nucleid+="s";
- else if (character.str()=="t") Nucleid+="t";
- else if (character.str()=="u") Nucleid+="u";
- else if (character.str()=="v") Nucleid+="v";
- else if (character.str()=="w") Nucleid+="w";
- else if (character.str()=="x") Nucleid+="x";
- else if (character.str()=="y") Nucleid+="y";
- else if (character.str()=="z") Nucleid+="z";
+ if(g4_exc == true){
+ continue;
+ }
+
+ 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";
+ else if (character.str()=="4") NumberOfMass+="4";
+ else if (character.str()=="5") NumberOfMass+="5";
+ else if (character.str()=="6") NumberOfMass+="6";
+ else if (character.str()=="7") NumberOfMass+="7";
+ else if (character.str()=="8") NumberOfMass+="8";
+ else if (character.str()=="9") NumberOfMass+="9";
+
+ else if (character.str()=="A") Nucleid+="A";
+ else if (character.str()=="B") Nucleid+="B";
+ else if (character.str()=="C") Nucleid+="C";
+ else if (character.str()=="D") Nucleid+="D";
+ else if (character.str()=="E") Nucleid+="E";
+ else if (character.str()=="F") Nucleid+="F";
+ else if (character.str()=="G") Nucleid+="G";
+ else if (character.str()=="H") Nucleid+="H";
+ else if (character.str()=="I") Nucleid+="I";
+ else if (character.str()=="J") Nucleid+="J";
+ else if (character.str()=="K") Nucleid+="K";
+ else if (character.str()=="L") Nucleid+="L";
+ else if (character.str()=="M") Nucleid+="M";
+ else if (character.str()=="N") Nucleid+="N";
+ else if (character.str()=="O") Nucleid+="O";
+ else if (character.str()=="P") Nucleid+="P";
+ else if (character.str()=="Q") Nucleid+="Q";
+ else if (character.str()=="R") Nucleid+="R";
+ else if (character.str()=="S") Nucleid+="S";
+ else if (character.str()=="T") Nucleid+="T";
+ else if (character.str()=="U") Nucleid+="U";
+ else if (character.str()=="V") Nucleid+="V";
+ else if (character.str()=="W") Nucleid+="W";
+ else if (character.str()=="X") Nucleid+="X";
+ else if (character.str()=="Y") Nucleid+="Y";
+ else if (character.str()=="Z") Nucleid+="Z";
+
+ else if (character.str()=="a") Nucleid+="a";
+ else if (character.str()=="b") Nucleid+="b";
+ else if (character.str()=="c") Nucleid+="c";
+ else if (character.str()=="d") Nucleid+="d";
+ else if (character.str()=="e") Nucleid+="e";
+ else if (character.str()=="f") Nucleid+="f";
+ else if (character.str()=="g") Nucleid+="g";
+ else if (character.str()=="h") Nucleid+="h";
+ else if (character.str()=="i") Nucleid+="i";
+ else if (character.str()=="j") Nucleid+="j";
+ else if (character.str()=="k") Nucleid+="k";
+ else if (character.str()=="l") Nucleid+="l";
+ else if (character.str()=="m") Nucleid+="m";
+ else if (character.str()=="n") Nucleid+="n";
+ else if (character.str()=="o") Nucleid+="o";
+ else if (character.str()=="p") Nucleid+="p";
+ else if (character.str()=="q") Nucleid+="q";
+ else if (character.str()=="r") Nucleid+="r";
+ else if (character.str()=="s") Nucleid+="s";
+ else if (character.str()=="t") Nucleid+="t";
+ else if (character.str()=="u") Nucleid+="u";
+ else if (character.str()=="v") Nucleid+="v";
+ else if (character.str()=="w") Nucleid+="w";
+ else if (character.str()=="x") Nucleid+="x";
+ else if (character.str()=="y") Nucleid+="y";
+ else if (character.str()=="z") Nucleid+="z";
}
string FinalName=NumberOfMass+Nucleid;
diff --git a/NPLib/Physics/NPQFS.cxx b/NPLib/Physics/NPQFS.cxx
index 2f402568c..b19851ca6 100644
--- a/NPLib/Physics/NPQFS.cxx
+++ b/NPLib/Physics/NPQFS.cxx
@@ -87,7 +87,7 @@ QFS::QFS(){
fPerpMomentumHist = NULL;
fParMomentumHist = NULL;
-
+ fDeexcitation = false;
}
////////////////////////////////////////////////////////////////////////////////
@@ -174,6 +174,9 @@ void QFS::ReadConfigurationFile(NPL::InputParser parser){
TH1F* Partemp = Read1DProfile(file_par[0], file_par[1]);
SetParMomentumHist(Partemp);
}
+ if(blocks[i]->HasToken("Deexcitation")){
+ fDeexcitation = blocks[i]->GetInt("Deexcitation");
+ }
}
cout << "\033[0m" ;
diff --git a/NPLib/Physics/NPQFS.h b/NPLib/Physics/NPQFS.h
index 0ad758152..0afde4ccd 100644
--- a/NPLib/Physics/NPQFS.h
+++ b/NPLib/Physics/NPQFS.h
@@ -93,6 +93,7 @@ namespace NPL{
bool fshoot1; // shoot light ejectile &
bool fshoot2; // shoot light ejectile 2
bool fUseExInGeant4;
+ bool fDeexcitation;
public:
Particle GetParticle(string name, NPL::InputParser parser);
@@ -182,6 +183,7 @@ namespace NPL{
bool GetShoot1() const {return fshoot1;}
bool GetShoot2() const {return fshoot2;}
bool GetShootB() const {return fshootB;}
+ bool GetDeexcitation() const {return fDeexcitation;}
double GetThetaCM() const {return fThetaCM;}
double GetPhiCM() const {return fPhiCM;}
double GetMomentumSigma() const {return fMomentumSigma;}
diff --git a/NPLib/Physics/TInteractionCoordinates.cxx b/NPLib/Physics/TInteractionCoordinates.cxx
index 340696b47..3914ff94f 100644
--- a/NPLib/Physics/TInteractionCoordinates.cxx
+++ b/NPLib/Physics/TInteractionCoordinates.cxx
@@ -35,6 +35,7 @@ TInteractionCoordinates::~TInteractionCoordinates()
{}
void TInteractionCoordinates::Clear(){
+ fDetected_Index.clear();
// Incident particle properties (before interactions in the target)
// Energy and Time
fDetected_Energy.clear();
@@ -46,6 +47,13 @@ void TInteractionCoordinates::Clear(){
// Incident particle angles
fDetected_Angle_Theta.clear();
fDetected_Angle_Phi.clear();
+
+ fDetected_Particle_Name.clear();
+ fDetected_A.clear();
+ fDetected_Z.clear();
+ fDetected_Mass.clear();
+ fDetected_Charge.clear();
+ fDetected_Brho.clear();
}
diff --git a/NPLib/Physics/TInteractionCoordinates.h b/NPLib/Physics/TInteractionCoordinates.h
index 1ce1052cc..0867793c2 100644
--- a/NPLib/Physics/TInteractionCoordinates.h
+++ b/NPLib/Physics/TInteractionCoordinates.h
@@ -11,7 +11,8 @@
* Original Author: N. de Sereville contact address: deserevi@ipno.in2p3.fr *
* *
* Creation Date : 10/06/09 *
- * Last update : *
+ * Last update : 01/09/2021 Valerian Alcindor adding particle name to *
+ * interaction coordinate for easier g4 analysis simulation *
*---------------------------------------------------------------------------*
* Decription: This class mainly records the coordinates of interaction *
* between a particle and a detector. *
@@ -25,14 +26,18 @@
+#include "NPFunction.h"
#include <vector>
#include "TObject.h"
#include <iostream>
+
using namespace std ;
class TInteractionCoordinates : public TObject{
private:
+ // TrackID or index for correlations
+ vector<int> fDetected_Index;
// Detected particle properties (before interactions in the target)
// Energy and Time
vector<double> fDetected_Energy;
@@ -44,6 +49,13 @@ class TInteractionCoordinates : public TObject{
// Particle angles
vector<double> fDetected_Angle_Theta;
vector<double> fDetected_Angle_Phi;
+ // Particle characteristics
+ vector<std::string> fDetected_Particle_Name;
+ vector<int> fDetected_A;
+ vector<int> fDetected_Z;
+ vector<double> fDetected_Mass;
+ vector<int> fDetected_Charge;
+ vector<double> fDetected_Brho;
public:
TInteractionCoordinates();
@@ -64,6 +76,37 @@ class TInteractionCoordinates : public TObject{
fDetected_Angle_Phi.push_back(Phi);
}
+ void SetInteraction(const int& Index, const double& Energy, const double&Time, const double& PositionX, const double& PositionY, const double& PositionZ,const double& Theta, const double& Phi){
+ fDetected_Index.push_back(Index);
+ fDetected_Energy.push_back(Energy);
+ fDetected_Time.push_back(Time);
+ fDetected_Position_X.push_back(PositionX);
+ fDetected_Position_Y.push_back(PositionY);
+ fDetected_Position_Z.push_back(PositionZ);
+ fDetected_Angle_Theta.push_back(Theta);
+ fDetected_Angle_Phi.push_back(Phi);
+ }
+
+ void SetInteraction(const int& Index, const double& Energy, const double&Time, const double& PositionX, const double& PositionY, const double& PositionZ,const double& Theta, const double& Phi, const std::string &ParticleName, const int &A, const int &Z, const double &Mass, const int &Charge, const double &Brho){
+ fDetected_Index.push_back(Index);
+ fDetected_Energy.push_back(Energy);
+ fDetected_Time.push_back(Time);
+ fDetected_Position_X.push_back(PositionX);
+ fDetected_Position_Y.push_back(PositionY);
+ fDetected_Position_Z.push_back(PositionZ);
+ fDetected_Angle_Theta.push_back(Theta);
+ fDetected_Angle_Phi.push_back(Phi);
+ if(ParticleName != "e-" && ParticleName != "e+")
+ fDetected_Particle_Name.push_back(NPL::ChangeNameFromG4Standard(ParticleName));
+ else
+ fDetected_Particle_Name.push_back(ParticleName);
+ fDetected_A.push_back(A);
+ fDetected_Z.push_back(Z);
+ fDetected_Mass.push_back(Mass);
+ fDetected_Charge.push_back(Charge);
+ fDetected_Brho.push_back(Brho);
+ }
+
///////////////////// SETTERS ////////////////////////
// Incident particle properties (before interactions in the target)
// Vertex of interaction
@@ -74,11 +117,18 @@ class TInteractionCoordinates : public TObject{
void SetDetectedAngleTheta(const double& AngleTheta) {fDetected_Angle_Theta.push_back(AngleTheta);}//!
void SetDetectedAnglePhi(const double& AnglePhi) {fDetected_Angle_Phi.push_back(AnglePhi);}//!
+ void SetDetectedParticleName(const std::string& ParticleName) {fDetected_Particle_Name.push_back(ParticleName);}//!
+ void SetDetectedA(const int& A) {fDetected_A.push_back(A);}//!
+ void SetDetectedZ(const int& Z) {fDetected_Z.push_back(Z);}//!
+ void SetDetectedMass(const double& Mass) {fDetected_Mass.push_back(Mass);}//!
+ void SetDetectedCharge(const int& Charge) {fDetected_Charge.push_back(Charge);}//!
+ void SetDetectedBrho(const double& Brho) {fDetected_Brho.push_back(Brho);}//!
+
///////////////////// GETTERS ////////////////////////
// Number of interactions (multiplicity)
int GetDetectedMultiplicity() const {return fDetected_Position_X.size();}
// Incident particle properties (before interactions in the target)
- // Enery and Time
+ // Energy and Time
double GetEnergy(const int& i) const {return fDetected_Energy[i];}//!
double GetTime(const int& i) const {return fDetected_Time[i];}//!
// Vertex of interaction
@@ -89,6 +139,14 @@ class TInteractionCoordinates : public TObject{
double GetDetectedAngleTheta(const int& i) const {return fDetected_Angle_Theta[i];}//!
double GetDetectedAnglePhi(const int& i) const {return fDetected_Angle_Phi[i];}//!
+ std::string GetParticleName(const int& i) const {return fDetected_Particle_Name[i];}//!
+
+ int GetA(const int& i) const {return fDetected_A[i];}//!
+ int GetZ(const int& i) const {return fDetected_Z[i];}//!
+ double GetMass(const int& i) const {return fDetected_Mass[i];}//!
+ int GetCharge(const int& i) const {return fDetected_Charge[i];}//!
+ double GetBrho(const int& i) const {return fDetected_Brho[i];}//!
+
ClassDef(TInteractionCoordinates, 2) // InteractionCoordinates structure
};
diff --git a/NPLib/Physics/TTrackInfo.cxx b/NPLib/Physics/TTrackInfo.cxx
new file mode 100644
index 000000000..5ca389023
--- /dev/null
+++ b/NPLib/Physics/TTrackInfo.cxx
@@ -0,0 +1,83 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2016 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: N. de Sereville contact address: deserevi@ipno.in2p3.fr *
+ * *
+ * Creation Date : 10/06/09 *
+ * Last update : 04/09/09 *
+ *---------------------------------------------------------------------------*
+ * Decription: This class records all the information concerning the event *
+ * generators, e.g. vertex of interaction, angles of emitted *
+ * particles... *
+ * This class derives from TObject (ROOT) and its aim is to be *
+ * stored in the output TTree of the G4 simulation *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * + 04/09/09: Add private members for emittance (N. de Sereville) *
+ * *
+ * *
+ *****************************************************************************/
+
+#include <iostream>
+using namespace std;
+
+#include "TTrackInfo.h"
+
+ClassImp(TTrackInfo)
+////////////////////////////////////////////////////////////////////////////////
+TTrackInfo::TTrackInfo(){
+}
+////////////////////////////////////////////////////////////////////////////////
+TTrackInfo::~TTrackInfo(){
+}
+////////////////////////////////////////////////////////////////////////////////
+void TTrackInfo::Clear(){
+
+ // emmitted particle
+ fTI_Particle_Name.clear();
+ fTI_Volume_Name.clear();
+ fTI_Kinetic_Energy.clear();
+ fTI_Mass.clear();
+ fTI_Charge.clear();
+ fTI_Z.clear();
+ fTI_A.clear();
+ fTI_Brho.clear();
+ fTI_Momentum.clear();
+ fTI_Momentum_X.clear();
+ fTI_Momentum_Y.clear();
+ fTI_Momentum_Z.clear();
+ fTI_PositionX.clear();
+ fTI_PositionY.clear();
+ fTI_PositionZ.clear();
+ fTI_Index.clear();
+ fTI_Theta.clear();
+ fTI_Phi.clear();
+}
+////////////////////////////////////////////////////////////////////////////////
+void TTrackInfo::Dump() const{
+ cout << "--------- Initial Condition Dump ---------" << endl ;
+
+ // emmitted particle
+ unsigned int size = fTI_Particle_Name.size();
+ for(unsigned int i = 0 ; i < size; i ++){
+ cout << "\t ---- Particle " << i << " ---- " << endl;
+ cout << "\t Particle Name" << fTI_Particle_Name[i] << endl;
+ cout << "\t Energy" << fTI_Kinetic_Energy[i] << endl;
+ cout << "\t Momentum Direction: ( "
+ << fTI_Momentum_X[i] << " ; "
+ << fTI_Momentum_Y[i] << " ; "
+ << fTI_Momentum_Z[i] << ")" << endl;
+ }
+}
+////////////////////////////////////////////////////////////////////////////////
+TVector3 TTrackInfo::GetParticleDirection (const int &i) const {
+ return TVector3( fTI_Momentum_X[i],
+ fTI_Momentum_Y[i],
+ fTI_Momentum_Z[i]).Unit();
+}
+
diff --git a/NPLib/Physics/TTrackInfo.h b/NPLib/Physics/TTrackInfo.h
new file mode 100644
index 000000000..889da398f
--- /dev/null
+++ b/NPLib/Physics/TTrackInfo.h
@@ -0,0 +1,144 @@
+#ifndef __TRACKINFO__
+#define __TRACKINFO__
+/*****************************************************************************
+ * Copyright (C) 2009-2016 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: N. de Sereville contact address: deserevi@ipno.in2p3.fr *
+ * *
+ * Creation Date : 10/06/09 *
+ * Last update : 04/09/09 *
+ *---------------------------------------------------------------------------*
+ * Decription: This class records all the information concerning the event *
+ * generators, e.g. vertex of interaction, angles of emitted *
+ * particles... *
+ * This class derives from TObject (ROOT) and its aim is to be *
+ * stored in the output TTree of the G4 simulation *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * + 04/09/09: Add private members for emittance (N. de Sereville) *
+ * *
+ * *
+ *****************************************************************************/
+
+// STL Header
+#include <cmath>
+#include <string>
+#include <vector>
+using namespace std;
+
+// Root Header
+#include "TObject.h"
+#include "TVector3.h"
+
+// NPTOOL headers
+#include "NPGlobalSystemOfUnits.h"
+#include "NPPhysicalConstants.h"
+#ifdef NP_SYSTEM_OF_UNITS_H
+using namespace NPUNITS;
+#endif
+
+class TTrackInfo : public TObject {
+private:
+ // Track info
+
+public:
+ // Particles info
+ vector<string> fTI_Particle_Name;
+ vector<string> fTI_Volume_Name;
+ vector<double> fTI_Kinetic_Energy;
+ vector<double> fTI_Theta;
+ vector<double> fTI_Phi;
+
+ vector<double> fTI_Mass;
+ vector<double> fTI_Charge;
+ vector<double> fTI_Z;
+ vector<double> fTI_A;
+ vector<double> fTI_Brho;
+
+ // vector<double> fTI_Spin;
+
+ vector<TVector3> fTI_Momentum;
+ vector<double> fTI_Momentum_X;
+ vector<double> fTI_Momentum_Y;
+ vector<double> fTI_Momentum_Z;
+
+ vector<double> fTI_PositionX;
+ vector<double> fTI_PositionY;
+ vector<double> fTI_PositionZ;
+
+ vector<int> fTI_Index;
+
+ // vector<double> fTI_;
+
+public:
+ TTrackInfo();
+ virtual ~TTrackInfo();
+
+ void Clear();
+ void Clear(const Option_t*) { Clear(); };
+ void Dump() const;
+
+ // emmitted particle
+ void SetParticleName(const string& Particle_Name) {
+ fTI_Particle_Name.push_back(Particle_Name);
+ }
+ void SetVolumeName(const string& Volume_Name) {
+ fTI_Volume_Name.push_back(Volume_Name);
+ }
+ void SetKineticEnergy(const double& Kinetic_Energy) {
+ fTI_Kinetic_Energy.push_back(Kinetic_Energy);
+ }
+ void SetTheta(const double& Theta) {
+ fTI_Theta.push_back(Theta);
+ }
+ void SetPhi(const double& Phi) {
+ fTI_Phi.push_back(Phi);
+ }
+ void SetMass(const double& Mass) { fTI_Mass.push_back(Mass); }
+ void SetCharge(const double& Charge) { fTI_Charge.push_back(Charge); }
+ void SetA(const double& A) { fTI_A.push_back(A); }
+ void SetZ(const double& Z) { fTI_Z.push_back(Z); }
+ void SetBrho(const double& Brho) { fTI_Brho.push_back(Brho); }
+ void SetMomentum(const TVector3& Momentum) {
+ fTI_Momentum.push_back(Momentum);
+ }
+ void SetMomentumX(const double& Momentum_X) {
+ fTI_Momentum_X.push_back(Momentum_X);
+ }
+ void SetMomentumY(const double& Momentum_Y) {
+ fTI_Momentum_Y.push_back(Momentum_Y);
+ }
+ void SetMomentumZ(const double& Momentum_Z) {
+ fTI_Momentum_Z.push_back(Momentum_Z);
+ }
+ void SetPositionX(const double& X) { fTI_PositionX.push_back(X); }
+ void SetPositionY(const double& Y) { fTI_PositionY.push_back(Y); }
+ void SetPositionZ(const double& Z) { fTI_PositionZ.push_back(Z); }
+ void SetIndex(const int& Index) { fTI_Index.push_back(Index); }
+
+ // emmitted particle
+ int GetParticleMultiplicity() const { return fTI_Kinetic_Energy.size(); }
+ string GetParticleName(const int& i) const { return fTI_Particle_Name[i]; }
+ double GetKineticEnergy(const int& i) const { return fTI_Kinetic_Energy[i]; }
+ TVector3 GetMomentum(const int& i) const { return fTI_Momentum[i]; }
+ double GetMomentumX(const int& i) const { return fTI_Momentum_X[i]; }
+ double GetMomentumY(const int& i) const { return fTI_Momentum_Y[i]; }
+ double GetMomentumZ(const int& i) const { return fTI_Momentum_Z[i]; }
+
+ TVector3 GetParticleDirection(const int& i) const;
+
+ double GetThetaLab_WorldFrame(const int& i) const {
+ return (GetParticleDirection(i).Angle(TVector3(0, 0, 1))) / deg;
+ }
+
+ unsigned int GetEmittedMult() const { return fTI_Particle_Name.size(); }
+
+ ClassDef(TTrackInfo, 1) // TrackInfo structure
+};
+
+#endif
diff --git a/NPSimulation/Core/CMakeLists.txt b/NPSimulation/Core/CMakeLists.txt
index 3d096b55f..9493bc823 100644
--- a/NPSimulation/Core/CMakeLists.txt
+++ b/NPSimulation/Core/CMakeLists.txt
@@ -1,2 +1,2 @@
add_library(NPSCore SHARED $<TARGET_OBJECTS:NPSEventGenerator> $<TARGET_OBJECTS:NPSProcess> SteppingAction.cc EventAction.cc PrimaryGeneratorAction.cc Target.cc Chamber.cc PrimaryGeneratorActionMessenger.cc NPSVDetector.cc DetectorConstruction.cc MaterialManager.cc DetectorMessenger.cc MyMagneticField.cc SteppingVerbose.cc NPSDetectorFactory.cc RunAction.cc Particle.cc ParticleStack.cc NPSFunction.cc)
-target_link_libraries(NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} NPSScorers NPInitialConditions NPInteractionCoordinates NPReactionConditions NPFissionConditions )
+target_link_libraries(NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} NPSScorers NPInitialConditions NPInteractionCoordinates NPReactionConditions NPFissionConditions NPTrackInfo)
diff --git a/NPSimulation/Core/ParticleStack.cc b/NPSimulation/Core/ParticleStack.cc
index 9d52f4024..9c30eba23 100644
--- a/NPSimulation/Core/ParticleStack.cc
+++ b/NPSimulation/Core/ParticleStack.cc
@@ -90,7 +90,7 @@ void ParticleStack::AddParticleToStack(Particle& particle){
// If the particle is the first one to be added, then the IC are cleared
if(m_First)
- m_InitialConditions->Clear();
+ m_InitialConditions->Clear();
m_ParticleStack.push_back(particle);
diff --git a/NPSimulation/Core/SteppingAction.cc b/NPSimulation/Core/SteppingAction.cc
index 8d964bcf2..eb2e65e8b 100644
--- a/NPSimulation/Core/SteppingAction.cc
+++ b/NPSimulation/Core/SteppingAction.cc
@@ -6,12 +6,12 @@
*****************************************************************************/
/*****************************************************************************
- * Original Author: Adrien MATTA contact address: matta@lpccaen.in2p3.fr *
+ * Original Author: ValerianAlcindor contact address: valcindor@@ikp.tu-darmstadt.de
* *
- * Creation Date : January 2021 *
+ * Creation Date : September 2021 *
* Last update : *
*---------------------------------------------------------------------------*
- * Decription: *
+ * Decription: [] *
* A quite Standard Geant4 EventAction class. *
* Call the Fill method of the output tree. *
*---------------------------------------------------------------------------*
@@ -20,35 +20,125 @@
* *
*****************************************************************************/
#include "SteppingAction.hh"
+#include "G4UnitsTable.hh"
+#include "NPFunction.h"
#include "NPOptionManager.h"
#include "RootOutput.h"
- //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-SteppingAction::SteppingAction(){
- m_record_track=NPOptionManager::getInstance()->GetRecordTrack();
- m_tree = RootOutput::getInstance()->GetTree();
- // Attach track info class to m_tree
- }
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+SteppingAction::SteppingAction() {
+ m_record_track = false;
+ m_tree = RootOutput::getInstance()->GetTree();
+
+ m_First = true;
+ ParticleName = "";
+ KineticEnergy = -1000;
+ Theta = -1000;
+ Phi = -1000;
+ Mass = -1000;
+ Charge = -1000;
+ Z = -1000;
+ A = -1000;
+ Momentum.setX(-1000);
+ Momentum.setY(-1000);
+ Momentum.setZ(-1000);
+ Position.setX(-1000);
+ Position.setY(-1000);
+ Position.setZ(-1000);
+ LastStepNumber = -1000;
+ VolumeName = "";
+ nClear = 0;
+ TrackID = -1000;
- //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void SteppingAction::UserSteppingAction (const G4Step* /*step*/){
- if(m_record_track){
- //TrackRecording(step);
- }
+ m_TrackInfo = new TTrackInfo();
+ AttachTrackInfo();
+ if (!RootOutput::getInstance()->GetTree()->FindBranch("TrackInfo"))
+ RootOutput::getInstance()->GetTree()->Branch("TrackInfo", "TTrackInfo",
+ &m_TrackInfo);
+
+ // Attach track info class to m_tree
}
-
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void SteppingAction::TrackRecording(const G4Step* /*step*/){
-// Clear the tracks class
-// m_Tracks->Clear();
-// unsigned int size = traj->size();
- // for(unsigned int i = 0 ; i < size ; i++){
- // FILL
- // Particle name
- // Interaction points
- // Momentum
- // Kinetic energy
- // Volume Name
- // }
+void SteppingAction::AttachTrackInfo() {
+ // Reasssigned the branch address
+ if (RootOutput::getInstance()->GetTree()->FindBranch("TrackInfo"))
+ RootOutput::getInstance()->GetTree()->SetBranchAddress("TrackInfo",
+ &m_TrackInfo);
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void SteppingAction::UserSteppingAction(const G4Step* step) {
+ if (m_record_track) {
+ TrackRecording(step);
}
+}
+
+// FIXME Still underdeveloppement
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void SteppingAction::TrackRecording(const G4Step* step) {
+ if (m_First)
+ m_TrackInfo->Clear();
+ m_First = false;
+
+ G4Track* track = step->GetTrack();
+ int StepNumber = track->GetCurrentStepNumber();
+
+ if (eventID
+ < G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID()) {
+ m_TrackInfo->Clear();
+ nClear++;
+ }
+
+ if (StepNumber == 1) {
+ ParticleName = track->GetParticleDefinition()->GetParticleName();
+ KineticEnergy = track->GetDynamicParticle()->GetKineticEnergy();
+ Mass = track->GetDynamicParticle()->GetMass();
+ Charge = track->GetDynamicParticle()->GetCharge();
+ Z = track->GetParticleDefinition()->GetAtomicNumber();
+ A = track->GetParticleDefinition()->GetAtomicMass();
+ Momentum = track->GetDynamicParticle()->GetMomentum();
+ Theta = Momentum.theta() * 180. / M_PI;
+ Phi = Momentum.phi() * 180. / M_PI;
+ Position = track->GetPosition();
+ G4VPhysicalVolume* volume = track->GetVolume();
+ VolumeName = volume->GetName();
+ TrackID = track->GetTrackID();
+
+ double c_light = 299.792458; // To go from T.m to MeV/e
+ double Brho = sqrt(KineticEnergy * KineticEnergy + 2 * KineticEnergy * Mass)
+ / (c_light * Charge);
+
+ m_TrackInfo->SetKineticEnergy(KineticEnergy);
+ m_TrackInfo->SetTheta(Theta);
+ m_TrackInfo->SetPhi(Phi);
+ m_TrackInfo->SetMass(Mass);
+ m_TrackInfo->SetCharge(Charge);
+ m_TrackInfo->SetZ(Z);
+ m_TrackInfo->SetA(A);
+ m_TrackInfo->SetBrho(Brho);
+ TVector3 Mom;
+ Mom.SetX(Momentum.x());
+ Mom.SetY(Momentum.y());
+ Mom.SetZ(Momentum.z());
+
+ m_TrackInfo->SetMomentum(Mom);
+
+ m_TrackInfo->SetMomentumX(Momentum.x());
+ m_TrackInfo->SetMomentumY(Momentum.y());
+ m_TrackInfo->SetMomentumZ(Momentum.z());
+
+ m_TrackInfo->SetPositionX(Position.x());
+ m_TrackInfo->SetPositionY(Position.y());
+ m_TrackInfo->SetPositionZ(Position.z());
+
+ m_TrackInfo->SetVolumeName(VolumeName);
+ m_TrackInfo->SetIndex(eventID + TrackID);
+
+ if (ParticleName != "e-" && ParticleName != "e+")
+ m_TrackInfo->SetParticleName(NPL::ChangeNameFromG4Standard(ParticleName));
+ else
+ m_TrackInfo->SetParticleName(ParticleName);
+ }
+ eventID = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
+}
diff --git a/NPSimulation/Core/SteppingAction.hh b/NPSimulation/Core/SteppingAction.hh
index 9a2b9b55a..73ab8a184 100644
--- a/NPSimulation/Core/SteppingAction.hh
+++ b/NPSimulation/Core/SteppingAction.hh
@@ -24,27 +24,56 @@
// G4 header defining G4 types
#include "globals.hh"
+// NPL
+#include "TTrackInfo.h"
+
// G4 header
+#include "G4RunManager.hh"
+#include "G4Track.hh"
#include "G4UserSteppingAction.hh"
-//Root
+// Root
#include "TTree.h"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-class SteppingAction : public G4UserSteppingAction{
- public:
- SteppingAction();
- ~SteppingAction(){};
-
- public:
- void TrackRecording(const G4Step* step);
- void UserSteppingAction (const G4Step* step);
-
- private: // tree
- TTree* m_tree;
-
- private:
- bool m_record_track;
+class SteppingAction : public G4UserSteppingAction {
+public:
+ SteppingAction();
+ ~SteppingAction(){};
+
+public:
+ void TrackRecording(const G4Step* step);
+ void UserSteppingAction(const G4Step* step);
+
+ std::string ParticleName;
+ double KineticEnergy;
+ double Theta;
+ double Phi;
+ double Mass;
+ double Charge;
+ double A;
+ double Z;
+ G4ThreeVector Momentum;
+ G4ThreeVector Position;
+ std::string VolumeName;
+ int nClear;
+ int eventID;
+ int TrackID;
+
+private: // tree
+ TTree* m_tree;
+ bool m_First;
+ unsigned int LastStepNumber;
+
+ // Host the Initial conditions TObject
+ TTrackInfo* m_TrackInfo;
+
+private:
+ bool m_record_track;
+
+public:
+ // Attach the TrackInfo object to the tree
+ void AttachTrackInfo();
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index ef429ff28..189d397dc 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -47,6 +47,8 @@ NPS::BeamReaction::BeamReaction(G4String modelName, G4Region* envelope)
m_shoot=false;
m_rand=0;
m_Z=0;
+
+ ABLA = new G4AblaInterface();
}
////////////////////////////////////////////////////////////////////////////////
@@ -149,17 +151,17 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
bool is_first = (to_entrance==0);
if(is_first && m_shoot){
- /* Does occur rarely when event is tangent to the target surface and scatters out
- std::cout << "Something went wrong in beam reaction, m_shoot and is_first variables cannot be true simultaneously" << std::endl;
- std::cout << "m_shoot: " << m_shoot << std::endl;
- std::cout << "rand: " << m_rand << std::endl;
- std::cout << "to_entrance: " << to_entrance << std::endl;
- std::cout << "to_exit: " << to_exit << std::endl;
- std::cout << "length: " << m_length << std::endl;
- std::cout << "step: " << m_StepSize << std::endl;
- std::cout << "Z: " << m_Z << std::endl;
- std::cout << "S: " << m_S << std::endl;
- */
+ /* Does occur rarely when event is tangent to the target surface and scatters out
+ std::cout << "Something went wrong in beam reaction, m_shoot and is_first variables cannot be true simultaneously" << std::endl;
+ std::cout << "m_shoot: " << m_shoot << std::endl;
+ std::cout << "rand: " << m_rand << std::endl;
+ std::cout << "to_entrance: " << to_entrance << std::endl;
+ std::cout << "to_exit: " << to_exit << std::endl;
+ std::cout << "length: " << m_length << std::endl;
+ std::cout << "step: " << m_StepSize << std::endl;
+ std::cout << "Z: " << m_Z << std::endl;
+ std::cout << "S: " << m_S << std::endl;
+ */
m_shoot = false;
}
@@ -556,8 +558,39 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
fastStep.CreateSecondaryTrack(particle2, localPosition, time);
}
if (m_QFS.GetShootB()) {
- G4DynamicParticle particleB(HeavyName, momentum_kineB_world, TKEB);
- fastStep.CreateSecondaryTrack(particleB, localPosition, time);
+ if(m_QFS.GetDeexcitation()){
+ double P_B2 = P_B->Px()*P_B->Px() + P_B->Py()*P_B->Py() + P_B->Pz()*P_B->Pz();
+ double scaling = 1;
+ if(P_B2>0.0){
+ double P_Bnew2 = TKEB*TKEB + 2 * TKEB * m_QFS.GetParticleB()->Mass();
+ scaling = sqrt(P_Bnew2)/sqrt(P_B2);
+ }
+ G4LorentzVector G4HeavyMomentum;
+ G4HeavyMomentum.setE(P_B->Energy());
+ G4HeavyMomentum.setPx(P_B->X());
+ G4HeavyMomentum.setPy(P_B->Y());
+ G4HeavyMomentum.setPz(P_B->Z());
+ G4Fragment HeavyFragment(HeavyName->GetAtomicMass(), HeavyName->GetAtomicNumber(), G4HeavyMomentum);
+ G4ReactionProductVector* HeavyDeexcitation = ABLA->DeExcite(HeavyFragment);
+
+ unsigned int sizeFrag = HeavyDeexcitation->size();
+ G4ThreeVector Momsum(0,0,0);
+ for(unsigned int i = 0 ; i < sizeFrag ; i++){
+ const G4ParticleDefinition* HeavyFragName = HeavyDeexcitation->at(i)->GetDefinition();
+ G4ThreeVector HeavyFragMomentum = HeavyDeexcitation->at(i)->GetMomentum();
+ double TKEFrag = HeavyDeexcitation->at(i)->GetKineticEnergy();
+ double HeavyFragModule = sqrt(TKEFrag*TKEFrag + 2*TKEFrag*HeavyDeexcitation->at(i)->GetMass());
+ G4ThreeVector HeavyFragMomentumDirection(HeavyFragMomentum.x()/HeavyFragModule,
+ HeavyFragMomentum.y()/HeavyFragModule,
+ HeavyFragMomentum.z()/HeavyFragModule);
+ G4DynamicParticle particleFrag(HeavyFragName, HeavyFragMomentumDirection, TKEFrag);
+ fastStep.CreateSecondaryTrack(particleFrag, localPosition, time);
+ }
+ }
+ else{
+ G4DynamicParticle particleB(HeavyName, momentum_kineB_world, TKEB);
+ fastStep.CreateSecondaryTrack(particleB, localPosition, time);
+ }
}
///////////////////////////////////
diff --git a/NPSimulation/Process/BeamReaction.hh b/NPSimulation/Process/BeamReaction.hh
index e127c74c2..9e24725e0 100644
--- a/NPSimulation/Process/BeamReaction.hh
+++ b/NPSimulation/Process/BeamReaction.hh
@@ -25,6 +25,9 @@
#define BeamReaction_h
#include "G4VFastSimulationModel.hh"
+#include "G4Abla.hh"
+#include "G4AblaInterface.hh"
+#include "G4Fragment.hh"
#include "PhysicsList.hh"
#include "NPReaction.h"
#include "NPQFS.h"
@@ -54,6 +57,7 @@ namespace NPS{
NPL::QFS m_QFS;
string m_BeamName;
int m_ReactionType;
+ G4AblaInterface* ABLA;
bool m_active;// is the process active
bool m_shoot;
diff --git a/NPSimulation/Scorers/InteractionScorers.cc b/NPSimulation/Scorers/InteractionScorers.cc
index 00f040f1e..4f0966282 100644
--- a/NPSimulation/Scorers/InteractionScorers.cc
+++ b/NPSimulation/Scorers/InteractionScorers.cc
@@ -24,7 +24,7 @@ using namespace InteractionScorers ;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
vector<InteractionData>::iterator InteractionDataVector::find(const unsigned int& index){
- for(vector<InteractionData>::iterator it= m_Data.begin() ; it !=m_Data.end() ; it++){
+ for(vector<InteractionData>::iterator it= m_Data.begin() ; it !=m_Data.end() ; it++){
if((*it).GetIndex()==index)
return it;
}
@@ -35,8 +35,8 @@ vector<InteractionData>::iterator InteractionDataVector::find(const unsigned int
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PS_Interactions::PS_Interactions(G4String name,TInteractionCoordinates* Inter, int depth) :G4VPrimitiveScorer(name, depth){
- m_Level = depth;
- m_InterractionCoordinates=Inter;
+ m_Level = depth;
+ m_InterractionCoordinates=Inter;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -46,26 +46,34 @@ G4bool PS_Interactions::ProcessHits(G4Step* aStep, G4TouchableHistory*){
t_Energy = aStep->GetTotalEnergyDeposit();
t_Time = point->GetGlobalTime();
t_Position = point->GetPosition();
-
+ t_Index = aStep->GetTrack()->GetTrackID();
+ t_ParticleName = aStep->GetTrack()->GetParticleDefinition()->GetParticleName();
+ t_A = aStep->GetTrack()->GetParticleDefinition()->GetAtomicMass();
+ t_Z = aStep->GetTrack()->GetParticleDefinition()->GetAtomicNumber();
+ t_Mass = aStep->GetTrack()->GetDynamicParticle()->GetMass();
+ t_Charge = aStep->GetTrack()->GetDynamicParticle()->GetCharge();
+ double KineticEnergy = aStep->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
+ t_Brho = sqrt(KineticEnergy * KineticEnergy + 2 * KineticEnergy * t_Mass)
+ / (c_light * t_Charge);
+
// add it to check the theta of momentum
// MOMENT = aStep->GetPreStepPoint()->GetMomentumDirection();
-
- t_Index = aStep->GetTrack()->GetTrackID();
+
vector<InteractionData>::iterator it;
it = m_DataVector.find(t_Index);
if(it!=m_DataVector.end())
it->Add(t_Energy);
else
- { m_DataVector.Set(t_Index,t_Energy,t_Time,t_Position.x(),t_Position.y(),t_Position.z(), /*MOMENT*/ t_Position.theta(),t_Position.phi());
-
-}
+ { m_DataVector.Set(t_Index,t_Energy,t_Time,t_Position.x(),t_Position.y(),t_Position.z(), /*MOMENT*/ t_Position.theta(),t_Position.phi(), t_ParticleName, t_A, t_Z, t_Mass, t_Charge, t_Brho);
+
+ }
return TRUE;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PS_Interactions::Initialize(G4HCofThisEvent*){
- // Clear is called by EventAction
+ // Clear is called by EventAction
clear();
}
@@ -73,9 +81,15 @@ void PS_Interactions::Initialize(G4HCofThisEvent*){
void PS_Interactions::EndOfEvent(G4HCofThisEvent*){
unsigned int size = m_DataVector.size();
- for(unsigned int i = 0 ; i < size ; i++)
- m_InterractionCoordinates->SetInteraction(m_DataVector[i]->GetEnergy(),m_DataVector[i]->GetTime(),m_DataVector[i]->GetPositionX(),m_DataVector[i]->GetPositionY(),m_DataVector[i]->GetPositionZ(),m_DataVector[i]->GetTheta()/deg,m_DataVector[i]->GetPhi()/deg);
-
+ for(unsigned int i = 0 ; i < size ; i++){
+ int eventID = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
+ int Index = t_Index + eventID;
+ m_InterractionCoordinates->SetInteraction(Index,m_DataVector[i]->GetEnergy(),m_DataVector[i]->GetTime(),
+ m_DataVector[i]->GetPositionX(),m_DataVector[i]->GetPositionY(),m_DataVector[i]->GetPositionZ(),
+ m_DataVector[i]->GetTheta()/deg,m_DataVector[i]->GetPhi()/deg, m_DataVector[i]->GetParticleName(),
+ m_DataVector[i]->GetA(), m_DataVector[i]->GetZ(), m_DataVector[i]->GetMass(), m_DataVector[i]->GetCharge(),
+ m_DataVector[i]->GetBrho());
+ }
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPSimulation/Scorers/InteractionScorers.hh b/NPSimulation/Scorers/InteractionScorers.hh
index f5b096c7e..f6a3789eb 100644
--- a/NPSimulation/Scorers/InteractionScorers.hh
+++ b/NPSimulation/Scorers/InteractionScorers.hh
@@ -21,6 +21,7 @@
* *
*****************************************************************************/
#include "G4VPrimitiveScorer.hh"
+#include "G4RunManager.hh"
#include "TInteractionCoordinates.h"
#include "NPImage.h"
#include <map>
@@ -44,6 +45,27 @@ namespace InteractionScorers {
m_Phi = Phi;
}
+ InteractionData(const unsigned int& Index ,const double& Energy, const double& Time ,
+ const double& PositionX, const double& PositionY, const double& PositionZ,
+ const double& Theta, const double& Phi, const std::string &ParticleName,
+ const int &A, const int &Z, const double &Mass, const int &Charge, const double &Brho){
+ m_Index = Index;
+ m_Energy = Energy;
+ m_Time = Time;
+ m_PositionX = PositionX;;
+ m_PositionY = PositionY;
+ m_PositionZ = PositionZ;
+ m_Theta = Theta;
+ m_Phi = Phi;
+ m_ParticleName = ParticleName;
+ m_A = A;
+ m_Z = Z;
+ m_Mass = Mass;
+ m_Charge = Charge;
+ m_Brho = Brho;
+
+ }
+
~InteractionData(){};
private:
@@ -55,7 +77,12 @@ namespace InteractionScorers {
double m_PositionZ;
double m_Theta;
double m_Phi;
-
+ std::string m_ParticleName;
+ int m_A;
+ int m_Z;
+ double m_Mass;
+ int m_Charge;
+ double m_Brho;
public:
unsigned int GetIndex() const{return m_Index;};
@@ -66,10 +93,19 @@ namespace InteractionScorers {
double GetPositionZ() const{return m_PositionZ;};
double GetTheta() const{return m_Theta;};
double GetPhi() const{return m_Phi;};
+ std::string GetParticleName() const{return m_ParticleName;};
+ int GetA() const{return m_A;};
+ int GetZ() const{return m_Z;};
+ double GetMass() const{return m_Mass;};
+ int GetCharge() const{return m_Charge;};
+ double GetBrho() const{return m_Brho;};
public:
- void Set(const unsigned int& Index, const double& Energy, const double& Time , const double& PositionX, const double& PositionY, const double& PositionZ, const double& Theta, const double& Phi){
+ void Set(const unsigned int& Index, const double& Energy, const double& Time ,
+ const double& PositionX, const double& PositionY, const double& PositionZ,
+ const double& Theta, const double& Phi){
+
m_Index = Index;
m_Energy = Energy;
m_Time = Time;
@@ -79,6 +115,28 @@ namespace InteractionScorers {
m_Theta = Theta;
m_Phi = Phi;
}
+
+ void Set(const unsigned int& Index, const double& Energy, const double& Time ,
+ const double& PositionX, const double& PositionY, const double& PositionZ,
+ const double& Theta, const double& Phi, const std::string &ParticleName,
+ const int &A, const int &Z, const double &Mass, const int &Charge, const double &Brho){
+
+ m_Index = Index;
+ m_Energy = Energy;
+ m_Time = Time;
+ m_PositionX = PositionX;;
+ m_PositionY = PositionY;
+ m_PositionZ = PositionZ;
+ m_Theta = Theta;
+ m_Phi = Phi;
+ m_ParticleName = ParticleName;
+ m_A = A;
+ m_Z = Z;
+ m_Mass = Mass;
+ m_Charge = Charge;
+ m_Brho = Brho;
+
+ }
void Add(const double& Energy){m_Energy+=Energy;};
unsigned int GetIndex(){return m_Index;};
};
@@ -100,7 +158,22 @@ namespace InteractionScorers {
vector<InteractionData>::iterator begin() {return m_Data.begin();};
unsigned int size() {return m_Data.size();};
void Add(const unsigned int& index,const double& Energy) {find(index)->Add(Energy);};
- void Set(const unsigned int& index,const double& Energy, const double& Time , const double& PositionX, const double& PositionY, const double& PositionZ, const double& Theta, const double& Phi) {m_Data.push_back(InteractionData(index,Energy,Time,PositionX,PositionY,PositionZ,Theta,Phi));};
+
+ void Set(const unsigned int& index,const double& Energy, const double& Time ,
+ const double& PositionX, const double& PositionY, const double& PositionZ,
+ const double& Theta, const double& Phi) {
+
+ m_Data.push_back(InteractionData(index,Energy,Time,PositionX,PositionY,PositionZ,Theta,Phi));
+
+ };
+
+ void Set(const unsigned int& index,const double& Energy, const double& Time ,
+ const double& PositionX, const double& PositionY, const double& PositionZ,
+ const double& Theta, const double& Phi, const std::string &ParticleName,
+ const int &A, const int &Z, const double &Mass, const int &Charge, const double &Brho){
+ m_Data.push_back(InteractionData(index,Energy,Time,PositionX,PositionY,PositionZ,Theta,Phi,ParticleName, A, Z, Mass, Charge, Brho));
+ };
+
InteractionData* operator[](const unsigned int& i){return &m_Data[i];};
};
@@ -132,6 +205,12 @@ namespace InteractionScorers {
double t_Energy;
double t_Time;
G4ThreeVector t_Position;
+ std::string t_ParticleName;
+ int t_A;
+ int t_Z;
+ double t_Mass;
+ int t_Charge;
+ double t_Brho;
//G4ThreeVector MOMENT;
};
--
GitLab