diff --git a/NPLib/Physics/CMakeLists.txt b/NPLib/Physics/CMakeLists.txt
index 02cf69825ff129124df7a827be3c897be0ba6bf6..3a3487d33d49fba7e12f3ea1810ca7d43b665977 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 9de2b666211aa8facae9ba279fafad1734c3437d..5d59b7e8e9983a5df5243a8b92f35ce040a69298 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 2f402568cf90ae561f8784163dd1a5754e428436..b19851ca6e57f04a99f0d9bd7eae04f764df1235 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 0ad7581521d6d225089109a97ef0f46fcb5f0c2e..0afde4ccda55ae51395e2116f19b584d1c115b59 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 340696b472d4a82c4a3d4273f0c7cda2eb0985cf..3914ff94fd3e9645610b3803ad753255f908140c 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 1ce1052ccc2a33fd99e676c858bbb9320ce067e1..0867793c216680a6d41bd5afac53abab79ccd6c9 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 0000000000000000000000000000000000000000..5ca389023640b3b1e55be32b9824a22f90ae2161
--- /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 0000000000000000000000000000000000000000..889da398f9625672f83468fdbbb03ef08a40449f
--- /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 3d096b55f6d1868699a974acbd2fb5274ed972a7..9493bc823950d5773d20b5d0f2575f2e6cdfa72f 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 9d52f4024450a9376cee50aef637832be20a52ec..9c30eba23af31d0ba6b9cf464f1d1ae55b4d26b8 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 8d964bcf221710db6308aece5842a07561773c11..eb2e65e8b5a5da67e25ffd4aeeab17434919334c 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 9a2b9b55a3c4ad0dfa88ebcfbb81ba68a0d80c5a..73ab8a18441fdd71533ccdb0de74ef68ea742c31 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 ef429ff280ca90d92fa593f704680920801d0d4b..189d397dc3e91dbeb312c14bdf742b88424c6673 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 e127c74c29cb28b2e59bab60f84186b4b21b5ded..9e24725e06d6db01aa515e28f850c5f785aec758 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 00f040f1ecfc58a711fa0c45cdd064ca7ca31954..4f09662829119114550a1fcd7462afbb6f3e62e3 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 f5b096c7e3b4ad217092fe8ff91be0a6fb8594e7..f6a3789eb562329e89e014bcbda160b965df7f5f 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;
};