diff --git a/FullSimLight/include/MySteppingAction.hh b/FullSimLight/include/MySteppingAction.hh
index 3a93e5f569acff759b9c046ceb184a78f49c7c26..f82b42d6ba29645a76966788ebe994ba63b09dee 100644
--- a/FullSimLight/include/MySteppingAction.hh
+++ b/FullSimLight/include/MySteppingAction.hh
@@ -4,21 +4,126 @@
#include "globals.hh"
#include "G4UserSteppingAction.hh"
+#include "MyRunAction.hh"
class MyEventAction;
class G4Step;
+
class MySteppingAction : public G4UserSteppingAction {
public:
- MySteppingAction(MyEventAction* evtact);
+ MySteppingAction(MyEventAction* evtact, MyRunAction* run);
virtual ~MySteppingAction();
virtual void UserSteppingAction(const G4Step*);
+
+ // maps to hold info per volume/process/material per particle type
+ //typedef std::map<G4String, TH1*> HistoMap_t;
+ typedef std::map<G4String, int> HistoMap_t;
+ typedef std::map<G4String, HistoMap_t> HistoMapMap_t;
+
+
+ /// this holds all the data from individual threads that needs to be merged at EoR
+ struct Report
+ {
+ // distributions per volume per particle type
+ HistoMapMap_t histoMapMap_vol_stepSize;
+ HistoMapMap_t histoMapMap_vol_stepKineticEnergy;
+ HistoMapMap_t histoMapMap_vol_postStepKineticEnergy;
+ HistoMapMap_t histoMapMap_vol_stepPseudorapidity;
+ HistoMapMap_t histoMapMap_vol_stepEnergyDeposit;
+ HistoMapMap_t histoMapMap_vol_stepEnergyNonIonDeposit;
+ HistoMapMap_t histoMapMap_vol_stepSecondaryKinetic;
+ HistoMapMap_t histoMapMap_vol_numberOfSteps;
+ HistoMapMap_t histoMapMap_vol_numberOfStepsPerInitialE;
+ HistoMapMap_t histoMapMap_vol_trackLengthPerInitialE;
+ HistoMapMap_t histoMapMap_vol_InitialE;
+
+ // distributions per material per particle type
+ HistoMapMap_t histoMapMap_mat_stepSize;
+ HistoMapMap_t histoMapMap_mat_stepKineticEnergy;
+ HistoMapMap_t histoMapMap_mat_postStepKineticEnergy;
+ HistoMapMap_t histoMapMap_mat_stepPseudorapidity;
+ HistoMapMap_t histoMapMap_mat_stepEnergyDeposit;
+ HistoMapMap_t histoMapMap_mat_stepEnergyNonIonDeposit;
+ HistoMapMap_t histoMapMap_mat_stepSecondaryKinetic;
+ HistoMapMap_t histoMapMap_mat_numberOfSteps;
+ HistoMapMap_t histoMapMap_mat_numberOfStepsPerInitialE;
+ HistoMapMap_t histoMapMap_mat_trackLengthPerInitialE;
+ HistoMapMap_t histoMapMap_mat_InitialE;
+
+ // distributions per process per particle type
+ HistoMapMap_t histoMapMap_prc_stepSize;
+ HistoMapMap_t histoMapMap_prc_stepKineticEnergy;
+ HistoMapMap_t histoMapMap_prc_postStepKineticEnergy;
+ HistoMapMap_t histoMapMap_prc_stepPseudorapidity;
+ HistoMapMap_t histoMapMap_prc_stepEnergyDeposit;
+ HistoMapMap_t histoMapMap_prc_stepEnergyNonIonDeposit;
+ HistoMapMap_t histoMapMap_prc_stepSecondaryKinetic;
+ HistoMapMap_t histoMapMap_prc_numberOfSteps;
+ HistoMapMap_t histoMapMap_prc_numberOfStepsPerInitialE;
+ HistoMapMap_t histoMapMap_prc_trackLengthPerInitialE;
+ HistoMapMap_t histoMapMap_prc_InitialE;
+
+ // all atlas
+ HistoMapMap_t histoMapMap_numberOfSteps;
+ HistoMapMap_t histoMapMap_numberOfStepsPerInitialE;
+ HistoMapMap_t histoMapMap_trackLengthPerInitialE;
+ HistoMapMap_t histoMapMap_InitialE;
+ HistoMapMap_t histoMapMap_stepKinetic;
+ HistoMapMap_t histoMapMap_postStepKinetic;
+
+ // 2D maps
+ HistoMapMap_t histoMapMap2D_vol_RZ;
+ HistoMapMap_t histoMapMap2D_mat_RZ;
+ HistoMapMap_t histoMapMap2D_prc_RZ;
+ HistoMapMap_t histoMapMap2D_vol_RZ_E;
+ HistoMapMap_t histoMapMap2D_mat_RZ_E;
+
+ // rather complicated function that merges two maps
+ void mergeMaps(HistoMapMap_t &selfMap, const HistoMapMap_t& refMap);
+
+ // function needed by ActionToolBaseReport base class
+ void merge(const Report & rep);
+ };
+
+ /// ctor
+ //StepHistogram(const Config&);
+
+ const Report& getReport() const { return m_report; }
+
+
private:
+ //Pointer to the run action
+ MyRunAction* m_run;
MyEventAction* fEventAction;
+ // report
+ Report m_report;
+
+ // initialize and fill histogram in a map
+ void InitializeFillHistogram2D(HistoMapMap_t &hMapMap, const char* suffix,
+ const G4String& pdgId, const G4String& vol,
+ int nbinsx, double xmin, double xmax,
+ int nbinsy, double ymin, double ymax,
+ double valuex, double valuey, double weight);
+
+ void InitializeFillHistogram(HistoMapMap_t &hMapMap, const char* suffix,
+ const G4String& pdgId, const G4String& vol,
+ int nbins, double xmin, double xmax, double value, double weight);
+
+ void InitializeFillHistogram(HistoMapMap_t &hMapMap, const char* suffix,
+ const G4String& pdgId, const G4String& vol,
+ int nbins, double *edges, double value, double weight);
+
+ float m_initialKineticEnergyOfStep;
+ G4String m_initialVolume;
+ G4String m_initialMaterial;
+ G4String m_initialProcess;
+ int m_trackID;
+
};
#endif
diff --git a/FullSimLight/src/MyActionInitialization.cc b/FullSimLight/src/MyActionInitialization.cc
index fd052b06ed7be2d30636b498b78cb76f494feb5b..4f51a68ca24bb69821b2966f76295c1d039b3c6e 100644
--- a/FullSimLight/src/MyActionInitialization.cc
+++ b/FullSimLight/src/MyActionInitialization.cc
@@ -61,7 +61,7 @@ void MyActionInitialization::Build() const {
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 1011db89eff3d3e0c485212d07ab2ecc725bb867..bb3c91618c3cbc4528d011e57c97342421e01cf6 100644
--- a/FullSimLight/src/MyRunAction.cc
+++ b/FullSimLight/src/MyRunAction.cc
@@ -23,7 +23,7 @@ MyRunAction::MyRunAction(bool isGeantino, G4String geantinoMapFilename) : G4User
}
MyRunAction::~MyRunAction() {
- if(fIsGeantino)
+ //if(fIsGeantino)
delete G4AnalysisManager::Instance();
}
@@ -97,11 +97,40 @@ void MyRunAction::BeginOfRunAction(const G4Run* /*aRun*/){
analysisManager->SetP1YAxisTitle(fEtaInt_id,"Interaction Length #lambda");
}
+ }
+ //Create step histograms only if not geantino simulation
+ else
+ {
+ //here create the "DEFAULT ONES"
+ // Create analysis manager
+ G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
+ if (isMaster) {
+ fMasterAnalysisManager = analysisManager;
+ G4cout<<"MyRunAction::BeginOfRunAction, created MASTER istance of the G4AnalysisManager: "<<fMasterAnalysisManager<<G4endl;
+
+ }
+
+ G4cout<<"MyRunAction::BeginOfRunAction, created WORKER istance of the G4AnalysisManager: "<<analysisManager<<G4endl;
+
+ G4cout << "Using G4AnalysisManager type: " << analysisManager->GetType() << G4endl;
+ analysisManager->SetVerboseLevel(1);
+
+ G4String histoFileName = "Histo.root"; //TO DO : customizable
+ // Open output root file
+ G4cout<<"\n\nBeginOfRunAction: \n...create a root file using the G4AnalysisManager" << G4endl;
+ if (!analysisManager->OpenFile(histoFileName)){
+ G4cout<<"\nBeginOfRunAction ERROR: File cannot be opened!"<<G4endl;
+ exit(-1);
+ } else
+ G4cout<<"\n...output File "<<histoFileName<<" opened!"<<G4endl;
+
+ ///There are no global/default histograms to create
+
}
if (isMaster) {
- //G4cout<<"\nBeginOfRunAction isMaster, and fMasterAnalysisManager: "<<fMasterAnalysisManager<<G4endl;
+ G4cout<<"\nBeginOfRunAction isMaster, and fMasterAnalysisManager: "<<fMasterAnalysisManager<<G4endl;
std::vector<G4Region*>* regionVect = G4RegionStore::GetInstance();
int numRegions = regionVect->size();
@@ -141,7 +170,7 @@ void MyRunAction::BeginOfRunAction(const G4Run* /*aRun*/){
fTimer = new G4Timer();
fTimer->Start();
}
- //else G4cout<<"BeginOfRunAction isWorker, and fMasterAnalysisManager: "<<fMasterAnalysisManager<<G4endl;
+ else G4cout<<"BeginOfRunAction isWorker, and fMasterAnalysisManager: "<<fMasterAnalysisManager<<G4endl;
}
@@ -159,6 +188,20 @@ void MyRunAction::EndOfRunAction(const G4Run*) {
G4cout<<"Output File successfully saved and closed! " << G4endl;
}
+ }
+ //Histogram stepping action TO DO: remove duplication of code?
+ else {
+
+ auto analysisManager= G4AnalysisManager::Instance();
+ //Finalize analysisManager and Write out file
+ if (analysisManager->IsOpenFile()){
+ G4cout<<"\n\n EndOfRunAction, writing output file: "<<analysisManager->GetFileName() << G4endl;
+ analysisManager->Write();
+ G4cout<<"... closing histogramming step file: "<<analysisManager->GetFileName() << G4endl;
+ analysisManager->CloseFile();
+ G4cout<<"Output File successfully saved and closed! " << G4endl;
+ }
+
}
if (isMaster) {
fTimer->Stop();
diff --git a/FullSimLight/src/MySteppingAction.cc b/FullSimLight/src/MySteppingAction.cc
index 6a46615c84e93f149fd77e952eaed530662d9e63..dfe2b5d6a4fae5f73fe455ee8fff53c63269d7a4 100644
--- a/FullSimLight/src/MySteppingAction.cc
+++ b/FullSimLight/src/MySteppingAction.cc
@@ -7,17 +7,417 @@
#include "G4ParticleDefinition.hh"
-MySteppingAction::MySteppingAction(MyEventAction* evtact) : G4UserSteppingAction(), fEventAction(evtact) {}
+//#include "G4StepPoint.hh"
+#include "G4VProcess.hh"
+#include "G4DebuggingHelper.hh"
+
+
+MySteppingAction::MySteppingAction(MyEventAction* evtact, MyRunAction* run) : G4UserSteppingAction(),
+m_run(run),
+fEventAction(evtact),
+m_initialKineticEnergyOfStep(),
+m_initialVolume(),
+m_initialMaterial(),
+m_initialProcess()
+{}
MySteppingAction::~MySteppingAction() {}
void MySteppingAction::UserSteppingAction(const G4Step* theStep) {
+ std::cout<<"************"<<std::endl;
+ std::cout<<"User stepping action: START"<<std::endl;
G4double edep = theStep->GetTotalEnergyDeposit();
G4double stepl = theStep->GetStepLength();
G4Track* theTrack = theStep->GetTrack();
G4bool isCharged = (theTrack->GetDefinition()->GetPDGCharge()!=0.);
G4int primTrackID = static_cast<MyTrackInformation*>(theTrack->GetUserInformation())->GetPrimaryTrackID();
fEventAction->AddData(edep, stepl, isCharged, primTrackID-1);
+
+
+ // track
+ G4Track *tr = theStep->GetTrack();
+
+ // pre-step point
+ G4StepPoint *PreStepPoint = theStep->GetPreStepPoint();
+
+ // post-step point
+ G4StepPoint *PostStepPoint = theStep->GetPostStepPoint();
+
+ // pre-step kinetic energy
+ double stepKinetic = PreStepPoint->GetKineticEnergy();
+
+ // post-step kinetic energy
+ double postStepKinetic = PostStepPoint->GetKineticEnergy();
+
+ // pre-step position
+ const G4ThreeVector& myPos = PreStepPoint->GetPosition();
+
+ // particle name
+ G4String particleName = "nucleus";
+
+ if (!(tr->GetDefinition()->GetParticleType() == "nucleus"))
+ particleName = G4DebuggingHelpers::ClassifyParticle(tr->GetParticleDefinition());
+
+ // pre-step volume
+ G4String volumeName = PreStepPoint->GetPhysicalVolume()->GetName();
+ volumeName= G4DebuggingHelpers::ClassifyVolume(volumeName);
+
+
+ // pre-step material
+ G4String materialName = PreStepPoint->GetMaterial()->GetName();
+ materialName = G4DebuggingHelpers::ClassifyMaterial(materialName);
+
+ // process name (uses post-step point)
+ G4String processName = PostStepPoint->GetProcessDefinedStep() ?
+ PostStepPoint->GetProcessDefinedStep()->GetProcessName() : "Unknown";
+
+ // secondaries
+ const std::vector<const G4Track*>* secondaries = theStep->GetSecondaryInCurrentStep();
+
+
+ std::cout<<"particleName: "<<particleName<<std::endl;
+ std::cout<<"volumeName: "<<volumeName<<std::endl;
+ std::cout<<"materialName: "<<materialName<<std::endl;
+ std::cout<<"processName: "<<processName<<std::endl;
+
+ // 2D map
+ bool do2DHistograms = false; //TO DO - configurable member
+ if (do2DHistograms) {
+ float DepositedE = theStep->GetTotalEnergyDeposit();
+ InitializeFillHistogram2D(m_report.histoMapMap2D_vol_RZ, "vol_RZ", particleName, volumeName,
+ 2000, -20000, 20000, 1000, 0, 5000, myPos.getZ(), myPos.perp(), 1.);
+ InitializeFillHistogram2D(m_report.histoMapMap2D_mat_RZ, "mat_RZ", particleName, materialName,
+ 2000, -20000, 20000, 1000, 0, 5000, myPos.getZ(), myPos.perp(), 1.);
+ InitializeFillHistogram2D(m_report.histoMapMap2D_prc_RZ, "prc_RZ", particleName, processName,
+ 2000, -20000, 20000, 1000, 0, 5000, myPos.getZ(), myPos.perp(), 1.);
+ InitializeFillHistogram2D(m_report.histoMapMap2D_vol_RZ_E, "vol_RZ_E", particleName, volumeName,
+ 2000, -20000, 20000, 1000, 0, 5000, myPos.getZ(), myPos.perp(), DepositedE);
+ InitializeFillHistogram2D(m_report.histoMapMap2D_mat_RZ_E, "mat_RZ_E", particleName, materialName,
+ 2000, -20000, 20000, 1000, 0, 5000, myPos.getZ(), myPos.perp(), DepositedE);
+ }
+
+ // step length
+ InitializeFillHistogram(m_report.histoMapMap_vol_stepSize, "vol_stepLength", particleName, volumeName,
+ 1000, -12, 4, log10(theStep->GetStepLength()), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_mat_stepSize, "mat_stepLength", particleName, materialName,
+ 1000, -12, 4, log10(theStep->GetStepLength()), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_prc_stepSize, "prc_stepLength", particleName, processName,
+ 1000, -12, 4, log10(theStep->GetStepLength()), 1.);
+
+
+ // step pseudorapidity
+ InitializeFillHistogram(m_report.histoMapMap_vol_stepPseudorapidity, "vol_stepPseudorapidity", particleName, volumeName,
+ 200, -10, 10, myPos.eta(), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_mat_stepPseudorapidity, "mat_stepPseudorapidity", particleName, materialName,
+ 200, -10, 10, myPos.eta(), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_prc_stepPseudorapidity, "prc_stepPseudorapidity", particleName, processName,
+ 200, -10, 10, myPos.eta(), 1.);
+
+ // step kinetic energy
+ InitializeFillHistogram(m_report.histoMapMap_vol_stepKineticEnergy, "vol_stepKineticEnergy", particleName, volumeName,
+ 1000, -9, 7, log10(stepKinetic), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_mat_stepKineticEnergy, "mat_stepKineticEnergy", particleName, materialName,
+ 1000, -9, 7, log10(stepKinetic), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_prc_stepKineticEnergy, "prc_stepKineticEnergy", particleName, processName,
+ 1000, -9, 7, log10(stepKinetic), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_stepKinetic, "stepKineticEnergy", particleName, "AllATLAS",
+ 1000, -9, 7, log10(stepKinetic), 1.);
+
+ // post step kinetic energy
+ InitializeFillHistogram(m_report.histoMapMap_vol_postStepKineticEnergy, "vol_postStepKineticEnergy", particleName, volumeName,
+ 1000, -9, 7, log10(postStepKinetic), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_mat_postStepKineticEnergy, "mat_postStepKineticEnergy", particleName, materialName,
+ 1000, -9, 7, log10(postStepKinetic), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_prc_postStepKineticEnergy, "prc_postStepKineticEnergy", particleName, processName,
+ 1000, -9, 7, log10(postStepKinetic), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_postStepKinetic, "postStepKineticEnergy", particleName, "AllATLAS",
+ 1000, -9, 7, log10(postStepKinetic), 1.);
+
+ // step energy deposit
+ InitializeFillHistogram(m_report.histoMapMap_vol_stepEnergyDeposit, "vol_stepEnergyDeposit", particleName, volumeName,
+ 1000, -11, 3, log10(theStep->GetTotalEnergyDeposit()), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_mat_stepEnergyDeposit, "mat_stepEnergyDeposit", particleName, materialName,
+ 1000, -11, 3, log10(theStep->GetTotalEnergyDeposit()), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_prc_stepEnergyDeposit, "prc_stepEnergyDeposit", particleName, processName,
+ 1000, -11, 3, log10(theStep->GetTotalEnergyDeposit()), 1.);
+
+ // step non-ionizing energy deposit
+ InitializeFillHistogram(m_report.histoMapMap_vol_stepEnergyNonIonDeposit, "vol_stepEnergyNonIonDeposit", particleName, volumeName,
+ 1000, -11, 1, log10(theStep->GetNonIonizingEnergyDeposit()), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_mat_stepEnergyNonIonDeposit, "mat_stepEnergyNonIonDeposit", particleName, materialName,
+ 1000, -11, 1, log10(theStep->GetNonIonizingEnergyDeposit()), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_prc_stepEnergyNonIonDeposit, "prc_stepEnergyNonIonDeposit", particleName, processName,
+ 1000, -11, 1, log10(theStep->GetNonIonizingEnergyDeposit()), 1.);
+
+ // secondary kinetic energy
+ for (const auto &track : *secondaries) {
+ G4String secondary_particleName = G4DebuggingHelpers::ClassifyParticle(track->GetParticleDefinition());
+ InitializeFillHistogram(m_report.histoMapMap_vol_stepSecondaryKinetic, "vol_stepSecondaryKinetic", secondary_particleName, volumeName,
+ 1000, -7, 5, log10(track->GetKineticEnergy()), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_mat_stepSecondaryKinetic, "mat_stepSecondaryKinetic", secondary_particleName, materialName,
+ 1000, -7, 5, log10(track->GetKineticEnergy()), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_prc_stepSecondaryKinetic, "prc_stepSecondaryKinetic", secondary_particleName, processName,
+ 1000, -7, 5, log10(track->GetKineticEnergy()), 1.);
+ }
+
+ bool doGeneralHistograms = false; //TO DO: configurable
+ // stop here if 'general' histograms not activated
+ // _______________________________________________
+ if (!doGeneralHistograms)
+ return;
+
+ // first step (after initial step)
+ if (tr->GetCurrentStepNumber()==1) {
+ // initial kinetic energy
+ m_initialKineticEnergyOfStep = stepKinetic;
+
+ // initial volume/material/processes
+ m_initialVolume = volumeName;
+ m_initialMaterial = materialName;
+ m_initialProcess = processName;
+
+ // save track ID for checking if we later have the same track
+ m_trackID = tr->GetTrackID();
+
+ // initial energy
+ InitializeFillHistogram(m_report.histoMapMap_InitialE, "InitialE", particleName, "AllATLAS",
+ 1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), 1.0);
+ InitializeFillHistogram(m_report.histoMapMap_vol_InitialE, "vol_InitialE", particleName, m_initialVolume,
+ 1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), 1.0);
+ InitializeFillHistogram(m_report.histoMapMap_mat_InitialE, "mat_InitialE", particleName, m_initialMaterial,
+ 1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), 1.0);
+ InitializeFillHistogram(m_report.histoMapMap_prc_InitialE, "prc_InitialE", particleName, m_initialProcess,
+ 1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), 1.0);
+ }
+
+ // last step
+ if ( tr->GetTrackStatus() == 2 ) {
+ // assert to check if we have the correct track
+ if (not (tr->GetTrackID() == m_trackID)) {
+ std::cerr<<"Track ID changed between the assumed first step and the last."<<std::endl;
+ throw std::exception();
+ }
+ // number of steps
+ int nSteps = tr->GetCurrentStepNumber() + 1;
+ InitializeFillHistogram(m_report.histoMapMap_numberOfSteps, "numberOfSteps", particleName, "AllATLAS",
+ 10000, 0.5, 10000.5, nSteps, 1.);
+ InitializeFillHistogram(m_report.histoMapMap_vol_numberOfSteps, "vol_numberOfSteps", particleName, m_initialVolume,
+ 10000, 0.5, 10000.5, nSteps, 1.);
+ InitializeFillHistogram(m_report.histoMapMap_mat_numberOfSteps, "mat_numberOfSteps", particleName, m_initialMaterial,
+ 10000, 0.5, 10000.5, nSteps, 1.);
+ InitializeFillHistogram(m_report.histoMapMap_prc_numberOfSteps, "prc_numberOfSteps", particleName, m_initialProcess,
+ 10000, 0.5, 10000.5, nSteps, 1.);
+ // number of steps vs initial energy
+ InitializeFillHistogram(m_report.histoMapMap_numberOfStepsPerInitialE, "numberOfStepsPerInitialE", particleName, "AllATLAS",
+ 1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), nSteps);
+ InitializeFillHistogram(m_report.histoMapMap_vol_numberOfStepsPerInitialE, "vol_numberOfStepsPerInitialE", particleName, m_initialVolume,
+ 1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), nSteps);
+ InitializeFillHistogram(m_report.histoMapMap_mat_numberOfStepsPerInitialE, "mat_numberOfStepsPerInitialE", particleName, m_initialMaterial,
+ 1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), nSteps);
+ InitializeFillHistogram(m_report.histoMapMap_prc_numberOfStepsPerInitialE, "prc_numberOfStepsPerInitialE", particleName, m_initialProcess,
+ 1000, -9, 7, std::log10(m_initialKineticEnergyOfStep), nSteps);
+ // track length vs initial energy
+ InitializeFillHistogram(m_report.histoMapMap_trackLengthPerInitialE, "trackLengthPerInitialE", particleName, "AllATLAS",
+ 1000, -9, 7, log10(tr->GetTrackLength()), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_vol_trackLengthPerInitialE, "vol_trackLengthPerInitialE", particleName, m_initialVolume,
+ 1000, -9, 7, log10(tr->GetTrackLength()), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_mat_trackLengthPerInitialE, "mat_trackLengthPerInitialE", particleName, m_initialMaterial,
+ 1000, -9, 7, log10(tr->GetTrackLength()), 1.);
+ InitializeFillHistogram(m_report.histoMapMap_prc_trackLengthPerInitialE, "prc_trackLengthPerInitialE", particleName, m_initialProcess,
+ 1000, -9, 7, log10(tr->GetTrackLength()), 1.);
+ }
+ std::cout<<"User Stepping action DONE"<<std::endl;
+ std::cout<<"************"<<std::endl;
+}
+
+void MySteppingAction::InitializeFillHistogram2D(HistoMapMap_t &hMapMap, const char* suffix,
+ const G4String& particleName, const G4String& vol,
+ int nbinsx, double xmin, double xmax,
+ int nbinsy, double ymin, double ymax,
+ double valuex, double valuey, double weight)
+{
+ if ( hMapMap.find(vol) == hMapMap.end() ) {
+ // initialize HistoMap_t if not yet exist
+ hMapMap.emplace(vol,HistoMap_t());
+ }
+ HistoMap_t &hMap = hMapMap[vol];
+
+ auto analysisManager = G4AnalysisManager::Instance();
+
+ if ( hMap.find(particleName) == hMap.end() ) {
+ // initialize histogram if not yet exist
+ std::ostringstream stringStream;
+ stringStream << vol << "_" << particleName << "_" << suffix;
+ //HERE I create the histogram if it doesn't exist
+ //DONE in what follows
+ //hMap[particleName] = new TH2F(stringStream.str().c_str(), stringStream.str().c_str(), nbinsx, xmin, xmax, nbinsy, ymin, ymax);
+ int id = m_run->fMasterAnalysisManager->GetH2Id(stringStream.str().c_str(), false);
+ if(id<0) {
+ //G4cout<<"Geant4: Eta rad profile of "<<pathEtaRL<<" didn't exist, creating P1 now on MASTER: "<<m_run->fMasterAnalysisManager<<G4endl;
+ //const std::string name(detName+"_RL");
+ //hMap[particleName] =
+ hMap[particleName] = m_run->fMasterAnalysisManager->CreateH2(stringStream.str().c_str(), stringStream.str().c_str(), nbinsx, xmin, xmax, nbinsy, ymin, ymax);
+// m_run->fMasterAnalysisManager->SetH2XAxisTitle(hMap[particleName], "#eta");
+// m_run->fMasterAnalysisManager->SetP1YAxisTitle(hMap[particleName], "%X0");
+
+ }
+ id=analysisManager->GetH2Id(stringStream.str().c_str(), false);
+ if(id<0)
+ {
+ //G4cout<<"Geant4: Eta rad profile of "<<pathEtaRL<<" didn't exist, creating P1 now on WORKER: "<<analysisManager<<G4endl;
+ hMap[particleName] = analysisManager->CreateH2(stringStream.str().c_str(), stringStream.str().c_str(), nbinsx, xmin, xmax, nbinsy, ymin, ymax);
+// analysisManager->SetP1XAxisTitle(hMap[particleName], "#eta");
+// analysisManager->SetP1YAxisTitle(hMap[particleName], "%X0");
+
+ }
+
+ }
+ //TO DO ((TH2*)hMap[particleName])->Fill(valuex, valuey, weight);
+ analysisManager->FillH2(hMap[particleName], valuex, valuey, weight);
+
+}
+
+void MySteppingAction::InitializeFillHistogram(HistoMapMap_t &hMapMap, const char* suffix,
+ const G4String& particleName, const G4String& vol,
+ int nbins, double xmin, double xmax, double value, double weight)
+{
+ if ( hMapMap.find(vol) == hMapMap.end() ) {
+ // initialize HistoMap_t if not yet exist
+ hMapMap.emplace(vol,HistoMap_t());
+ }
+ HistoMap_t &hMap = hMapMap[vol];
+ auto analysisManager = G4AnalysisManager::Instance();
+ if ( hMap.find(particleName) == hMap.end() ) {
+ // initialize histogram if not yet exist
+ std::ostringstream stringStream;
+ stringStream << vol << "_" << particleName << "_" << suffix;
+
+ //hMap[particleName] = new TH1F(stringStream.str().c_str(), stringStream.str().c_str(), nbins, xmin, xmax);
+
+ int id = m_run->fMasterAnalysisManager->GetH1Id(stringStream.str().c_str(), false);
+ if(id<0) {
+ hMap[particleName] = m_run->fMasterAnalysisManager->CreateH1(stringStream.str().c_str(), stringStream.str().c_str(), nbins, xmin, xmax);
+
+ }
+ id=analysisManager->GetH1Id(stringStream.str().c_str(), false);
+ if(id<0)
+ {
+ hMap[particleName] = analysisManager->CreateH1(stringStream.str().c_str(), stringStream.str().c_str(), nbins, xmin, xmax);
+ }
+
+ }
+ //TO DO hMap[particleName]->Fill(value, weight);
+ analysisManager->FillH1(hMap[particleName], value, weight);
+}
+
+void MySteppingAction::InitializeFillHistogram(HistoMapMap_t &hMapMap, const char* suffix,
+ const G4String& particleName, const G4String& vol,
+ int nbins, double *edges, double value, double weight)
+{
+ if ( hMapMap.find(vol) == hMapMap.end() ) {
+ // initialize HistoMap_t if not yet exist
+ hMapMap.emplace(vol,HistoMap_t());
+ }
+ HistoMap_t &hMap = hMapMap[vol];
+ auto analysisManager = G4AnalysisManager::Instance();
+ if ( hMap.find(particleName) == hMap.end() ) {
+ // initialize histogram if not yet exist
+ std::ostringstream stringStream;
+ stringStream << vol << "_" << particleName << "_" << suffix;
+
+ //hMap[particleName] = new TH1F(stringStream.str().c_str(), stringStream.str().c_str(), nbins, edges);
+ int id = m_run->fMasterAnalysisManager->GetH1Id(stringStream.str().c_str(), false);
+ if(id<0) {
+ hMap[particleName] = m_run->fMasterAnalysisManager->CreateH1(stringStream.str().c_str(), stringStream.str().c_str(), nbins, edges[0], edges[nbins]); //TO DO : check if this is proper and actually used
+
+ }
+ id=analysisManager->GetH2Id(stringStream.str().c_str(), false);
+ if(id<0)
+ {
+ hMap[particleName] = analysisManager->CreateH1(stringStream.str().c_str(), stringStream.str().c_str(), nbins, edges[0], edges[nbins]);
+ }
+
+ }
+ //hMap[particleName]->Fill(value, weight);
+ analysisManager->FillH1(hMap[particleName], value, weight);
+
}
+
+void MySteppingAction::Report::mergeMaps(HistoMapMap_t &selfMap, const HistoMapMap_t& refMap) {
+ for (auto const& ref : refMap)
+ {
+ if ( selfMap.find(ref.first) == selfMap.end() ) {
+ // HistoMap_t does not yet exist
+ selfMap.emplace(ref.first, ref.second);
+ }
+ else {
+ HistoMap_t &target = selfMap[ref.first];
+ for (auto const& hm : ref.second)
+ {
+ if ( target.find(hm.first) == target.end() ) {
+ // histogram does not yet exist
+ target.emplace(hm.first, hm.second);
+ }
+ else {
+ // add histograms
+ //TO DO target[hm.first]->Add(hm.second);
+ }
+ }
+ }
+ }
+}
+
+void MySteppingAction::Report::merge(const Report & rep) {
+ mergeMaps(histoMapMap_vol_stepSize, rep.histoMapMap_vol_stepSize);
+ mergeMaps(histoMapMap_vol_stepKineticEnergy, rep.histoMapMap_vol_stepKineticEnergy);
+ mergeMaps(histoMapMap_vol_postStepKineticEnergy, rep.histoMapMap_vol_postStepKineticEnergy);
+ mergeMaps(histoMapMap_vol_stepPseudorapidity, rep.histoMapMap_vol_stepPseudorapidity);
+ mergeMaps(histoMapMap_vol_stepEnergyDeposit, rep.histoMapMap_vol_stepEnergyDeposit);
+ mergeMaps(histoMapMap_vol_stepEnergyNonIonDeposit, rep.histoMapMap_vol_stepEnergyNonIonDeposit);
+ mergeMaps(histoMapMap_vol_stepSecondaryKinetic, rep.histoMapMap_vol_stepSecondaryKinetic);
+
+ mergeMaps(histoMapMap_mat_stepSize, rep.histoMapMap_mat_stepSize);
+ mergeMaps(histoMapMap_mat_stepKineticEnergy, rep.histoMapMap_mat_stepKineticEnergy);
+ mergeMaps(histoMapMap_mat_postStepKineticEnergy, rep.histoMapMap_mat_postStepKineticEnergy);
+ mergeMaps(histoMapMap_mat_stepPseudorapidity, rep.histoMapMap_mat_stepPseudorapidity);
+ mergeMaps(histoMapMap_mat_stepEnergyDeposit, rep.histoMapMap_mat_stepEnergyDeposit);
+ mergeMaps(histoMapMap_mat_stepEnergyNonIonDeposit, rep.histoMapMap_mat_stepEnergyNonIonDeposit);
+ mergeMaps(histoMapMap_mat_stepSecondaryKinetic, rep.histoMapMap_mat_stepSecondaryKinetic);
+
+ mergeMaps(histoMapMap_prc_stepSize, rep.histoMapMap_prc_stepSize);
+ mergeMaps(histoMapMap_prc_stepKineticEnergy, rep.histoMapMap_prc_stepKineticEnergy);
+ mergeMaps(histoMapMap_prc_postStepKineticEnergy, rep.histoMapMap_prc_postStepKineticEnergy);
+ mergeMaps(histoMapMap_prc_stepPseudorapidity, rep.histoMapMap_prc_stepPseudorapidity);
+ mergeMaps(histoMapMap_prc_stepEnergyDeposit, rep.histoMapMap_prc_stepEnergyDeposit);
+ mergeMaps(histoMapMap_prc_stepEnergyNonIonDeposit, rep.histoMapMap_prc_stepEnergyNonIonDeposit);
+ mergeMaps(histoMapMap_prc_stepSecondaryKinetic, rep.histoMapMap_prc_stepSecondaryKinetic);
+
+ mergeMaps(histoMapMap_numberOfSteps, rep.histoMapMap_numberOfSteps);
+ mergeMaps(histoMapMap_vol_numberOfSteps, rep.histoMapMap_vol_numberOfSteps);
+ mergeMaps(histoMapMap_mat_numberOfSteps, rep.histoMapMap_mat_numberOfSteps);
+ mergeMaps(histoMapMap_prc_numberOfSteps, rep.histoMapMap_prc_numberOfSteps);
+ mergeMaps(histoMapMap_numberOfStepsPerInitialE, rep.histoMapMap_numberOfStepsPerInitialE);
+ mergeMaps(histoMapMap_vol_numberOfStepsPerInitialE, rep.histoMapMap_vol_numberOfStepsPerInitialE);
+ mergeMaps(histoMapMap_mat_numberOfStepsPerInitialE, rep.histoMapMap_mat_numberOfStepsPerInitialE);
+ mergeMaps(histoMapMap_prc_numberOfStepsPerInitialE, rep.histoMapMap_prc_numberOfStepsPerInitialE);
+ mergeMaps(histoMapMap_trackLengthPerInitialE, rep.histoMapMap_trackLengthPerInitialE);
+ mergeMaps(histoMapMap_vol_trackLengthPerInitialE, rep.histoMapMap_vol_trackLengthPerInitialE);
+ mergeMaps(histoMapMap_mat_trackLengthPerInitialE, rep.histoMapMap_mat_trackLengthPerInitialE);
+ mergeMaps(histoMapMap_prc_trackLengthPerInitialE, rep.histoMapMap_prc_trackLengthPerInitialE);
+ mergeMaps(histoMapMap_InitialE, rep.histoMapMap_InitialE);
+ mergeMaps(histoMapMap_vol_InitialE, rep.histoMapMap_vol_InitialE);
+ mergeMaps(histoMapMap_mat_InitialE, rep.histoMapMap_mat_InitialE);
+ mergeMaps(histoMapMap_prc_InitialE, rep.histoMapMap_prc_InitialE);
+ mergeMaps(histoMapMap_stepKinetic, rep.histoMapMap_stepKinetic);
+ mergeMaps(histoMapMap_postStepKinetic, rep.histoMapMap_postStepKinetic);
+
+ mergeMaps(histoMapMap2D_vol_RZ, rep.histoMapMap2D_vol_RZ);
+ mergeMaps(histoMapMap2D_mat_RZ, rep.histoMapMap2D_mat_RZ);
+ mergeMaps(histoMapMap2D_prc_RZ, rep.histoMapMap2D_prc_RZ);
+ mergeMaps(histoMapMap2D_vol_RZ_E, rep.histoMapMap2D_vol_RZ_E);
+ mergeMaps(histoMapMap2D_mat_RZ_E, rep.histoMapMap2D_mat_RZ_E);
+}
+
+