diff --git a/FullSimLight/fullSimLight.cc b/FullSimLight/fullSimLight.cc
index 90815c5afd03fe42daa096e860335680a9ed287c..db97f72be43370b6115cdddecd2b8289b0f0f4ba 100644
--- a/FullSimLight/fullSimLight.cc
+++ b/FullSimLight/fullSimLight.cc
@@ -36,6 +36,7 @@ static std::string parMacroFileName = "";
static std::string parPhysListName = "FTFP_BERT";
static bool parRunOverlapCheck = false;
static bool parSwitchField = false;
+static bool parStepStats = false;
void GetInputArguments(int argc, char** argv);
void Help();
@@ -48,12 +49,13 @@ int main(int argc, char** argv) {
G4cout
<< " =============== Running FullSimLight ================ " << G4endl
- << " Physics list name = " << parPhysListName << G4endl
- << " Geant4 macro = " << parMacroFileName << G4endl
- << " Performance mode = " << parIsPerformance << G4endl
- << " Geometry file = " << geometryFileName << G4endl
- << " Run Overlap Check = " << parRunOverlapCheck << G4endl
- << " Switch Field = " << parSwitchField << G4endl
+ << " Physics list name = " << parPhysListName << G4endl
+ << " Geant4 macro = " << parMacroFileName << G4endl
+ << " Performance mode = " << parIsPerformance << G4endl
+ << " Geometry file = " << geometryFileName << G4endl
+ << " Run Overlap Check = " << parRunOverlapCheck << G4endl
+ << " Switch Field = " << parSwitchField << G4endl
+ << " Collect Step Statistics = " << parStepStats << G4endl
<< " ===================================================== " << G4endl;
G4Timer myTotalCPUTimer;
@@ -104,7 +106,9 @@ int main(int argc, char** argv) {
}
// 3. User action
- runManager->SetUserInitialization(new MyActionInitialization(parIsPerformance));
+ G4String stepStatsFilename = "stepStats_FieldOn.root";
+ if (parSwitchField) stepStatsFilename = "stepStats_FieldOff.root";
+ runManager->SetUserInitialization(new MyActionInitialization(parIsPerformance, false, std::string(), parStepStats, stepStatsFilename));
// 4. Run the simulation in batch mode
G4UImanager* UI = G4UImanager::GetUIpointer();
@@ -137,6 +141,7 @@ static struct option options[] = {
{"pythia primary generator " , required_argument, 0, 'P'},
{"overlap geometry check " , no_argument , 0, 'o'},
{"switch field " , no_argument , 0, 's'},
+ {"collect step stats " , no_argument , 0, 't'},
{"help" , no_argument , 0, 'h'},
{0, 0, 0, 0}
};
@@ -171,7 +176,7 @@ void GetInputArguments(int argc, char** argv) {
}
while (true) {
int c, optidx = 0;
- c = getopt_long(argc, argv, "P:pm:f:sg:oh", options, &optidx);
+ c = getopt_long(argc, argv, "P:pm:f:stg:oh", options, &optidx);
if (c == -1)
break;
//
@@ -207,6 +212,9 @@ void GetInputArguments(int argc, char** argv) {
case 's':
parSwitchField = true;
break;
+ case 't':
+ parStepStats = true;
+ break;
case 'h':
Help();
exit(0);
diff --git a/FullSimLight/include/MyActionInitialization.hh b/FullSimLight/include/MyActionInitialization.hh
index 467ec7cb2853a88b1d7fd37df12d857d34aba2c3..a93d075f8510f740a8e8c12440618ee7aa5ac15b 100644
--- a/FullSimLight/include/MyActionInitialization.hh
+++ b/FullSimLight/include/MyActionInitialization.hh
@@ -7,7 +7,7 @@
class MyActionInitialization : public G4VUserActionInitialization {
public:
- MyActionInitialization(bool isperformance=false, bool createGeantinoMaps=false, G4String geantinoMapsFilename="geantinoMaps.root");
+ MyActionInitialization(bool isperformance=false, bool createGeantinoMaps=false, G4String geantinoMapsFilename="geantinoMaps.root", bool createStepStats=false, G4String stepStatsFilename="stepStats.root");
virtual ~MyActionInitialization();
virtual void BuildForMaster() const;
@@ -15,6 +15,7 @@ public:
void SetPerformanceModeFlag(bool val) { fIsPerformance = val; }
void SetCreateGeantinoMaps (bool val) { fCreateGeantinoMaps = val; }
+ void SetCreateStepStats (bool val) { fCreateStepStats = val; }
void SetCreateEtaPhiMaps (bool val) { fCreateEtaPhiMaps = val; }
void SetCreateDetectorsMaps(bool val) { fCreateDetectorsMaps = val; }
void SetCreateMaterialsMaps(bool val) { fCreateMaterialsMaps = val; }
@@ -32,6 +33,8 @@ private:
bool fCreateMaterialsMaps;
bool fCreateElementsMaps;
G4String fGeantinoMapsFilename;
+ bool fCreateStepStats;
+ G4String fStepStatsFilename;
G4double fRlimit;
G4double fZlimit;
G4double fXlimit;
diff --git a/FullSimLight/include/MyRunAction.hh b/FullSimLight/include/MyRunAction.hh
index 0b4ba95ffc31f530c9b3cfe49971b872156684fb..39b2471b8e473ec25118caf6dc95b057932149c9 100644
--- a/FullSimLight/include/MyRunAction.hh
+++ b/FullSimLight/include/MyRunAction.hh
@@ -16,14 +16,15 @@ class MyRunAction: public G4UserRunAction {
public:
- MyRunAction(bool isGeantino, G4String geantinoMapsFileName="geantinoMaps.root", bool isStepStats=false, G4String stepStatsFilename="stepStats.root");
+ MyRunAction(bool isGeantino, bool isStepStats, G4String geantinoMapsFileName="geantinoMaps.root", G4String stepStatsFilename="stepStats.root");
virtual ~MyRunAction();
virtual G4Run* GenerateRun();
virtual void BeginOfRunAction(const G4Run* aRun);
virtual void EndOfRunAction(const G4Run* aRun);
- void SetPerformanceFlag(G4bool val) { fIsPerformance=val; }
+ void SetPerformanceFlag(G4bool val) { fIsPerformance=val; }
+ bool getIsStepStats() { return fIsStepStats; }
private:
G4bool fIsPerformance;
@@ -49,6 +50,9 @@ public:
int fEtaInt_id;
// ID of the ZRSteps 2D Histogram
int fStepHeatMap_id;
+ // ID of the ESteps 1D Histogram
+ int fStepKinetic_id;
+
};
#endif
diff --git a/FullSimLight/include/MySteppingAction.hh b/FullSimLight/include/MySteppingAction.hh
index 3a93e5f569acff759b9c046ceb184a78f49c7c26..b86b06ee6a3319b5344eb6569baea927780008d6 100644
--- a/FullSimLight/include/MySteppingAction.hh
+++ b/FullSimLight/include/MySteppingAction.hh
@@ -5,20 +5,25 @@
#include "globals.hh"
#include "G4UserSteppingAction.hh"
+//G4AnalysisManager
+#include "MyAnalysis.hh"
+
class MyEventAction;
+class MyRunAction;
class G4Step;
class MySteppingAction : public G4UserSteppingAction {
public:
- MySteppingAction(MyEventAction* evtact);
+ MySteppingAction(MyEventAction* evtact, MyRunAction* runact);
virtual ~MySteppingAction();
virtual void UserSteppingAction(const G4Step*);
private:
- MyEventAction* fEventAction;
+ MyEventAction* fEventAction;
+ MyRunAction* fRunAction;
};
#endif
diff --git a/FullSimLight/src/MyActionInitialization.cc b/FullSimLight/src/MyActionInitialization.cc
index fd052b06ed7be2d30636b498b78cb76f494feb5b..8304e02daeb5f2d78c873b7ef333eea7a2f78c12 100644
--- a/FullSimLight/src/MyActionInitialization.cc
+++ b/FullSimLight/src/MyActionInitialization.cc
@@ -18,15 +18,15 @@
//const G4AnalysisManager* MyActionInitialization::fMasterAnalysisManager = nullptr;
-MyActionInitialization::MyActionInitialization(bool isperformance, bool createGeantinoMaps, G4String geantinoMapsFilename)
-: G4VUserActionInitialization(), fIsPerformance(isperformance), fCreateGeantinoMaps(createGeantinoMaps),fGeantinoMapsFilename(geantinoMapsFilename){}
+MyActionInitialization::MyActionInitialization(bool isperformance, bool createGeantinoMaps, G4String geantinoMapsFilename, bool createStepStats, G4String stepStatsFilename)
+: G4VUserActionInitialization(), fIsPerformance(isperformance), fCreateGeantinoMaps(createGeantinoMaps), fGeantinoMapsFilename(geantinoMapsFilename), fCreateStepStats(createStepStats), fStepStatsFilename(stepStatsFilename) {}
MyActionInitialization::~MyActionInitialization() {}
// called in case of MT
void MyActionInitialization::BuildForMaster() const {
- MyRunAction* masterRunAct = new MyRunAction(fCreateGeantinoMaps,fGeantinoMapsFilename);
+ MyRunAction* masterRunAct = new MyRunAction(fCreateGeantinoMaps, fCreateStepStats, fGeantinoMapsFilename, fStepStatsFilename);
masterRunAct->SetPerformanceFlag(fIsPerformance);
SetUserAction(masterRunAct);
}
@@ -47,21 +47,21 @@ void MyActionInitialization::Build() const {
// in sequential mode the BuildForMaster method is not called:
// - create the only one run action with perfomance flag true i.e. only time is measured
if (fIsPerformance) {
- MyRunAction* masterRunAct = new MyRunAction(fCreateGeantinoMaps, fGeantinoMapsFilename);
+ MyRunAction* masterRunAct = new MyRunAction(fCreateGeantinoMaps, fCreateStepStats, fGeantinoMapsFilename, fStepStatsFilename);
masterRunAct->SetPerformanceFlag(fIsPerformance);
SetUserAction(masterRunAct);
}
#endif
// do not create Run,Event,Stepping and Tracking actions in case of perfomance mode
if (!fIsPerformance) {
- MyRunAction* runact = new MyRunAction(fCreateGeantinoMaps, fGeantinoMapsFilename);
+ MyRunAction* runact = new MyRunAction(fCreateGeantinoMaps, fCreateStepStats, fGeantinoMapsFilename, fStepStatsFilename);
SetUserAction(runact);
if(!fCreateGeantinoMaps){
MyEventAction* evtact = new MyEventAction();
SetUserAction(evtact);
SetUserAction(new MyTrackingAction(evtact));
- SetUserAction(new MySteppingAction(evtact));
+ SetUserAction(new MySteppingAction(evtact, runact));
}
else
diff --git a/FullSimLight/src/MyRunAction.cc b/FullSimLight/src/MyRunAction.cc
index 0345bf643e723867752c0d76cea345ce7d08b417..328b50f7a9685b24dd9c37a840b3e7059c59109f 100644
--- a/FullSimLight/src/MyRunAction.cc
+++ b/FullSimLight/src/MyRunAction.cc
@@ -18,7 +18,7 @@
G4AnalysisManager* MyRunAction::fMasterAnalysisManager = nullptr;
-MyRunAction::MyRunAction(bool isGeantino, G4String geantinoMapFilename, bool isStepStats, G4String stepStatsFilename) : G4UserRunAction(), fIsPerformance(false), fIsGeantino(isGeantino), fIsStepStats(isStepStats), fRun(nullptr), fTimer(nullptr), fGeantinoMapsFilename(geantinoMapFilename), fStepStatsFilename(stepStatsFilename){
+MyRunAction::MyRunAction(bool isGeantino, bool isStepStats, G4String geantinoMapFilename, G4String stepStatsFilename) : G4UserRunAction(), fIsPerformance(false), fIsGeantino(isGeantino), fIsStepStats(isStepStats), fRun(nullptr), fTimer(nullptr), fGeantinoMapsFilename(geantinoMapFilename), fStepStatsFilename(stepStatsFilename){
}
@@ -106,7 +106,6 @@ void MyRunAction::BeginOfRunAction(const G4Run* /*aRun*/){
if (isMaster) {
fMasterAnalysisManager = analysisManager;
//G4cout<<"MyRunAction::BeginOfRunAction, created MASTER istance of the G4AnalysisManager: "<<fMasterAnalysisManager<<G4endl;
-
}
//else
// G4cout<<"MyRunAction::BeginOfRunAction, created WORKER istance of the G4AnalysisManager: "<<analysisManager<<G4endl;
@@ -116,19 +115,26 @@ void MyRunAction::BeginOfRunAction(const G4Run* /*aRun*/){
// Open output root file
G4cout<<"\n\nBeginOfRunAction: \n...create a root file using the G4AnalysisManager" << G4endl;
- if (!analysisManager->OpenFile(fStepStatsFilename)){
+ if (!analysisManager->OpenFile(fStepStatsFilename)) {
G4cout<<"\nBeginOfRunAction ERROR: File cannot be opened!"<<G4endl;
exit(-1);
- } else
+ }
+ else
G4cout<<"\n...output File "<<fStepStatsFilename<<" opened!"<<G4endl;
const char* stepHeatMapName = "ZRSteps";
if(analysisManager->GetH2Id(stepHeatMapName, false) < 0) {
fStepHeatMap_id = analysisManager->CreateH2(stepHeatMapName,stepHeatMapName,1000,-25000.,25000.,2000,0.,15000.);
- // G4cout<<"MyRunAction::BeginOfRunAction: G4AnalysisManager Created ZRSteps 2DHistogram with name: "<<stepHeatMapName<< " and with id: "<<fStepHeatMap_id<<G4endl;
- analysisManager->SetH2XAxisTitle(fStepHeatMap_id,"Z[mm]");
- analysisManager->SetH2YAxisTitle(fStepHeatMap_id,"R[mm]");
- analysisManager->SetH2ZAxisTitle(fStepHeatMap_id,"NumberOfSteps");
+ analysisManager->SetH2XAxisTitle(fStepHeatMap_id,"Z [mm]");
+ analysisManager->SetH2YAxisTitle(fStepHeatMap_id,"R [mm]");
+ analysisManager->SetH2ZAxisTitle(fStepHeatMap_id,"Steps");
+ }
+
+ const char* stepKineticName = "ESteps";
+ if(analysisManager->GetH1Id(stepKineticName, false) < 0) {
+ fStepKinetic_id = analysisManager->CreateH1(stepKineticName,stepKineticName,1000,-9,7);
+ analysisManager->SetH1XAxisTitle(fStepKinetic_id, "log(pre-step kinetic energy [MeV])");
+ analysisManager->SetH1YAxisTitle(fStepKinetic_id, "Steps");
}
}
@@ -180,7 +186,7 @@ void MyRunAction::BeginOfRunAction(const G4Run* /*aRun*/){
void MyRunAction::EndOfRunAction(const G4Run*) {
- if(fIsGeantino){
+ if(fIsGeantino || fIsStepStats){
auto analysisManager= G4AnalysisManager::Instance();
//Finalize analysisManager and Write out file
diff --git a/FullSimLight/src/MySteppingAction.cc b/FullSimLight/src/MySteppingAction.cc
index 5831fc401d7f86502da74c9cce9dfa09928e4a92..402d35f8eafec28689a7d6dc583e969a0813126b 100644
--- a/FullSimLight/src/MySteppingAction.cc
+++ b/FullSimLight/src/MySteppingAction.cc
@@ -2,14 +2,16 @@
#include "MySteppingAction.hh"
#include "MyEventAction.hh"
+#include "MyRunAction.hh"
#include "MyTrackInformation.hh"
#include "G4Step.hh"
#include "G4ParticleDefinition.hh"
#include "G4FieldManager.hh"
#include "G4Field.hh"
+#include <cmath>
-MySteppingAction::MySteppingAction(MyEventAction* evtact) : G4UserSteppingAction(), fEventAction(evtact) {}
+MySteppingAction::MySteppingAction(MyEventAction* evtact, MyRunAction* runact) : G4UserSteppingAction(), fEventAction(evtact), fRunAction(runact) {}
MySteppingAction::~MySteppingAction() {}
@@ -47,5 +49,34 @@ void MySteppingAction::UserSteppingAction(const G4Step* theStep) {
}
else
std::cout<<"Field manager does NOT exist!"<<std::endl;
+
*/
+
+ // Collect step statistics
+ G4bool doStepStats = fRunAction->getIsStepStats();
+ if (doStepStats) {
+
+ // pre-step point
+ const G4StepPoint* PreStepPoint = theStep->GetPreStepPoint();
+ // post-step point
+ const G4StepPoint* PostStepPoint = theStep->GetPostStepPoint();
+ // pre-step position
+ const G4ThreeVector position = PreStepPoint->GetPosition();
+ const G4double rXY = position.rho();
+ const G4double Z = position.z();
+ // pre-step kinetic energy
+ const G4double stepKinetic = PreStepPoint->GetKineticEnergy();
+
+ // Fill histograms
+ auto analysisManager = G4AnalysisManager::Instance();
+ {
+ // Protect concurrent histo filling
+ static std::mutex mutex_instance;
+ std::lock_guard<std::mutex> lock(mutex_instance);
+
+ analysisManager->FillH2(fRunAction->fStepHeatMap_id, Z, rXY);
+ analysisManager->FillH1(fRunAction->fStepKinetic_id, log10(stepKinetic));
+ }
+ }
+
}
\ No newline at end of file