diff --git a/FullSimLight/fullSimLight.cc b/FullSimLight/fullSimLight.cc
index c662a855cef33a5592d0f81abf18c8ced70ddd4f..04dde09d648ef1075973edc2316e2fea967212e9 100644
--- a/FullSimLight/fullSimLight.cc
+++ b/FullSimLight/fullSimLight.cc
@@ -41,10 +41,10 @@ void Help();
int main(int argc, char** argv) {
-
+
// Get input arguments
GetInputArguments(argc, argv);
-
+
G4cout
<< " =============== Running FullSimLight ================ " << G4endl
<< " Physics list name = " << parPhysListName << G4endl
@@ -53,10 +53,10 @@ int main(int argc, char** argv) {
<< " Geometry file = " << geometryFileName << G4endl
<< " Run Overlap Check = " << parRunOverlapCheck << G4endl
<< " ===================================================== " << G4endl;
-
+
G4Timer myTotalCPUTimer;
myTotalCPUTimer.Start();
-
+
//choose the Random engine: set to MixMax explicitely (default form 10.4)
G4Random::setTheEngine(new CLHEP::MixMaxRng);
// set seed and print info
@@ -67,7 +67,7 @@ int main(int argc, char** argv) {
<< " Initial seed = " << G4Random::getTheSeed() << G4endl
<< " ===================================================== " << G4endl
<< G4endl;
-
+
// Construct the default run manager
#ifdef G4MULTITHREADED
G4MTRunManager* runManager = new G4MTRunManager;
@@ -77,16 +77,16 @@ int main(int argc, char** argv) {
#else
G4RunManager* runManager = new G4RunManager;
#endif
-
+
// set mandatory initialization classes
// 1. Detector construction
MyDetectorConstruction* detector = new MyDetectorConstruction;
-
+
if (parRunOverlapCheck) detector->SetRunOverlapCheck(true);
-
+
detector->SetGeometryFileName (geometryFileName);
runManager->SetUserInitialization(detector);
-
+
// 2. Physics list
G4PhysListFactory factory;
if (factory.IsReferencePhysList(parPhysListName)) {
@@ -99,15 +99,15 @@ int main(int argc, char** argv) {
G4cerr << "ERROR: Physics List " << parPhysListName << " UNKNOWN!" << G4endl;
return -1;
}
-
+
// 3. User action
runManager->SetUserInitialization(new MyActionInitialization(parIsPerformance));
-
+
// 4. Run the simulation in batch mode
G4UImanager* UI = G4UImanager::GetUIpointer();
G4String command = "/control/execute ";
UI->ApplyCommand(command+parMacroFileName);
-
+
// Print out the final random number
G4cout << G4endl
<< " ================================================================= " << G4endl
@@ -150,7 +150,7 @@ void Help() {
<<" -P : use Pythia primary generator [config. available: ttbar/higgs/minbias or use a Pythia command input file]\n"
<<" -p : flag ==> run the application in performance mode i.e. no user actions \n"
<<" : - ==> run the application in NON performance mode i.e. with user actions (default) \n"<< std::endl;
-
+
std::cout <<"\nUsage: ./fullSimLight [OPTIONS] -m <MACRO_FILE>\n" <<std::endl;
for (int i=0; options[i].name!=NULL; i++) {
//printf("\t-%c --%s\t%s\n", options[i].val, options[i].name, options[i].has_arg ? options[i].name : "");
@@ -178,8 +178,6 @@ void GetInputArguments(int argc, char** argv) {
case 'P':
#if USE_PYTHIA
set_pythia_config(optarg);
- // Need to enable performance mode, as user actions require particle gun setup
- parIsPerformance = true;
#else
std::cerr << "Support for Pythia is not available. \nPlease visit the website http://home.thep.lu.se/Pythia/ to install it in your system." << std::endl;
exit(1);
diff --git a/FullSimLight/include/MyActionInitialization.hh b/FullSimLight/include/MyActionInitialization.hh
index 467ec7cb2853a88b1d7fd37df12d857d34aba2c3..a4eaab0910c2db3787711a2b19474c7baf8874c7 100644
--- a/FullSimLight/include/MyActionInitialization.hh
+++ b/FullSimLight/include/MyActionInitialization.hh
@@ -6,12 +6,14 @@
#include "G4String.hh"
class MyActionInitialization : public G4VUserActionInitialization {
+
public:
+
MyActionInitialization(bool isperformance=false, bool createGeantinoMaps=false, G4String geantinoMapsFilename="geantinoMaps.root");
- virtual ~MyActionInitialization();
+ ~MyActionInitialization() override;
- virtual void BuildForMaster() const;
- virtual void Build() const;
+ void BuildForMaster() const override;
+ void Build() const override;
void SetPerformanceModeFlag(bool val) { fIsPerformance = val; }
void SetCreateGeantinoMaps (bool val) { fCreateGeantinoMaps = val; }
@@ -25,12 +27,13 @@ public:
void SetYlimit(G4double y) { fYlimit = y; }
private:
- bool fIsPerformance;
- bool fCreateGeantinoMaps;
- bool fCreateEtaPhiMaps;
- bool fCreateDetectorsMaps;
- bool fCreateMaterialsMaps;
- bool fCreateElementsMaps;
+
+ bool fIsPerformance;
+ bool fCreateGeantinoMaps;
+ bool fCreateEtaPhiMaps;
+ bool fCreateDetectorsMaps;
+ bool fCreateMaterialsMaps;
+ bool fCreateElementsMaps;
G4String fGeantinoMapsFilename;
G4double fRlimit;
G4double fZlimit;
diff --git a/FullSimLight/include/MyEventAction.hh b/FullSimLight/include/MyEventAction.hh
index 2ac872dfa15590feaa8670c28d7a78034b9b73fa..f426fa8efe3f8a4454add5774dff9696ad9a1217 100644
--- a/FullSimLight/include/MyEventAction.hh
+++ b/FullSimLight/include/MyEventAction.hh
@@ -4,30 +4,28 @@
#include "G4UserEventAction.hh"
#include "globals.hh"
-#include "MyEventDataPerPrimary.hh"
-
-#include <vector>
+#include "MyEventData.hh"
class G4Event;
class G4Track;
-
class MyEventAction: public G4UserEventAction {
public:
- MyEventAction();
- virtual ~MyEventAction();
+ MyEventAction();
+ ~MyEventAction() override;
- virtual void BeginOfEventAction(const G4Event* evt);
- virtual void EndOfEventAction(const G4Event* evt);
+ void BeginOfEventAction(const G4Event* evt) override;
+ void EndOfEventAction(const G4Event* evt) override;
- void AddData(G4double edep, G4double length, G4bool ischarged, G4int primid);
- void AddSecondaryTrack(const G4Track* track, G4int primid);
+ void AddData(G4double edep, G4double length, G4bool ischarged);
+ void AddSecondaryTrack(const G4Track* track);
private:
- G4int fNumberOfPrimariesPerEvent;
- std::vector<MyEventDataPerPrimary> fEventDataPerPrimary;
+
+ MyEventData fEventData;
+
};
#endif
diff --git a/FullSimLight/include/MyEventDataPerPrimary.hh b/FullSimLight/include/MyEventData.hh
similarity index 78%
rename from FullSimLight/include/MyEventDataPerPrimary.hh
rename to FullSimLight/include/MyEventData.hh
index 4e020ead0233feb03b54ad9e5e56ab49ae41c295..f71ae9b5d7ad0c1f6367c105e76d4669fd69a52d 100644
--- a/FullSimLight/include/MyEventDataPerPrimary.hh
+++ b/FullSimLight/include/MyEventData.hh
@@ -1,19 +1,20 @@
-#ifndef MyEventDataPerPrimary_h
-#define MyEventDataPerPrimary_h 1
+#ifndef MyEventData_h
+#define MyEventData_h 1
#include "globals.hh"
#include <iostream>
-class MyEventDataPerPrimary {
+class MyEventData {
public:
- MyEventDataPerPrimary();
- ~MyEventDataPerPrimary();
+
+ MyEventData();
+ ~MyEventData();
void Clear();
- friend std::ostream& operator<<(std::ostream&, const MyEventDataPerPrimary&);
+ friend std::ostream& operator<<(std::ostream&, const MyEventData&);
G4double fEdep; // sum of energy deposit
G4double fTrackLCh; // sum of charged step length
diff --git a/FullSimLight/include/MyRun.hh b/FullSimLight/include/MyRun.hh
index 148da0f0d0c8f37b4fe1017aee1997bc2840b09b..0b7c2dbc086b6f0869bf6e2ecea79d0f4304ba32 100644
--- a/FullSimLight/include/MyRun.hh
+++ b/FullSimLight/include/MyRun.hh
@@ -4,29 +4,29 @@
#include "G4Run.hh"
#include "globals.hh"
-#include "MyRunDataPerPrimary.hh"
-
-#include <vector>
+#include "MyRunData.hh"
class G4Track;
-class MyEventDataPerPrimary;
+class MyEventData;
class MyRun : public G4Run {
+
public:
- MyRun();
- virtual ~MyRun();
- virtual void Merge(const G4Run*);
- void FillPerEventNumPrimaries(G4int primtypeindx) { fRunDataPerPrimary[primtypeindx].fNumPrimaries += 1.;}
- void FillPerEvent(const MyEventDataPerPrimary&, G4int primtypeindx);
+ MyRun();
+ ~MyRun() override;
+
+ void Merge(const G4Run*) override;
+
+ void FillPerEvent(const MyEventData&);
void EndOfRun();
- const MyRunDataPerPrimary& GetRunDataPerPrimary(G4int primtypeindx) const { return fRunDataPerPrimary[primtypeindx]; }
+ const MyRunData& GetRunData() const { return fRunData; }
private:
- G4int fNumPrimaryTypes;
- std::vector<MyRunDataPerPrimary> fRunDataPerPrimary;
+
+ MyRunData fRunData;
};
diff --git a/FullSimLight/include/MyRunAction.hh b/FullSimLight/include/MyRunAction.hh
index 3c0d1b13facbc4cdb62f7e0bd17ea9e15f5634c2..dd35c06a3686138fcb850fbf57dadd2e7e6d506f 100644
--- a/FullSimLight/include/MyRunAction.hh
+++ b/FullSimLight/include/MyRunAction.hh
@@ -17,25 +17,27 @@ class MyRunAction: public G4UserRunAction {
public:
MyRunAction(bool isGeantino, G4String geantinoMapsFileName="geantinoMaps.root");
- virtual ~MyRunAction();
+ ~MyRunAction();
- virtual G4Run* GenerateRun();
- virtual void BeginOfRunAction(const G4Run* aRun);
- virtual void EndOfRunAction(const G4Run* aRun);
-
- void SetPerformanceFlag(G4bool val) { fIsPerformance=val; }
+ G4Run* GenerateRun() override;
+ void BeginOfRunAction(const G4Run*) override;
+ void EndOfRunAction(const G4Run*) override;
+ void SetPerformanceFlag(G4bool val) { fIsPerformance = val; }
+ void SetPythiaConfig(G4String config) { fPythiaConfig = config; }
private:
+
G4bool fIsPerformance;
G4bool fIsGeantino;
MyRun* fRun;
G4Timer* fTimer;
G4String fGeantinoMapsFilename;
-
+ G4String fPythiaConfig;
+
//TO DO: make private and add Get methods
public:
-
+
static G4AnalysisManager* fMasterAnalysisManager;
// ID of the RadLen 2D Profile
int fRadName_id;
diff --git a/FullSimLight/include/MyRunDataPerPrimary.hh b/FullSimLight/include/MyRunData.hh
similarity index 81%
rename from FullSimLight/include/MyRunDataPerPrimary.hh
rename to FullSimLight/include/MyRunData.hh
index cce1c4793a7e5578ccee2d860cf3d6a5b2f031e4..3ba8d31e82ad711cff8676980cadfeec48648fc9 100644
--- a/FullSimLight/include/MyRunDataPerPrimary.hh
+++ b/FullSimLight/include/MyRunData.hh
@@ -1,19 +1,20 @@
-#ifndef MyRunDataPerPrimary_h
-#define MyRunDataPerPrimary_h 1
+#ifndef MyRunData_h
+#define MyRunData_h 1
#include "globals.hh"
-class MyRunDataPerPrimary {
+class MyRunData {
+
public:
- MyRunDataPerPrimary();
- ~MyRunDataPerPrimary();
+
+ MyRunData();
+ ~MyRunData();
void Clear();
- MyRunDataPerPrimary& operator+=(const MyRunDataPerPrimary& other);
+ MyRunData& operator+=(const MyRunData& other);
- G4double fNumPrimaries; // total number of primary particles
G4double fEdep; // sum of energy deposit (per event)
G4double fEdep2; // sum of energy deposit square
G4double fTrackLCh; // sum of charged step length (per event)
diff --git a/FullSimLight/include/MySteppingAction.hh b/FullSimLight/include/MySteppingAction.hh
index 3a93e5f569acff759b9c046ceb184a78f49c7c26..0af275c24e19c78f92dd2a405133e45afc37ed27 100644
--- a/FullSimLight/include/MySteppingAction.hh
+++ b/FullSimLight/include/MySteppingAction.hh
@@ -12,13 +12,15 @@ class MySteppingAction : public G4UserSteppingAction {
public:
- MySteppingAction(MyEventAction* evtact);
- virtual ~MySteppingAction();
+ MySteppingAction(MyEventAction*);
+ ~MySteppingAction() override;
- virtual void UserSteppingAction(const G4Step*);
+ void UserSteppingAction(const G4Step*) override;
private:
+
MyEventAction* fEventAction;
+
};
#endif
diff --git a/FullSimLight/include/MyTrackInformation.hh b/FullSimLight/include/MyTrackInformation.hh
deleted file mode 100644
index 6e8536d2c4ef332adac1866fe670339e540c0658..0000000000000000000000000000000000000000
--- a/FullSimLight/include/MyTrackInformation.hh
+++ /dev/null
@@ -1,21 +0,0 @@
-
-#ifndef MyTrackInformation_h
-#define MyTrackInformation_h 1
-
-#include "G4VUserTrackInformation.hh"
-#include "globals.hh"
-
-class MyTrackInformation : public G4VUserTrackInformation {
-public:
- MyTrackInformation(G4int id);
- virtual ~MyTrackInformation();
-
-
- G4int GetPrimaryTrackID() const { return fPrimaryTrackID; }
- void SetPrimaryTrackID(G4int id) { fPrimaryTrackID = id; }
-
-private:
- G4int fPrimaryTrackID;
-};
-
-#endif
diff --git a/FullSimLight/include/MyTrackingAction.hh b/FullSimLight/include/MyTrackingAction.hh
index 804a331758213e7576c67d6cb60eae1ff6c46454..262d5558269989113ac60e444aeec8c3f906c92c 100644
--- a/FullSimLight/include/MyTrackingAction.hh
+++ b/FullSimLight/include/MyTrackingAction.hh
@@ -9,16 +9,16 @@ class G4Track;
class MyEventAction;
class MyTrackingAction : public G4UserTrackingAction {
+
public:
- MyTrackingAction(MyEventAction*);
- virtual ~MyTrackingAction();
+ MyTrackingAction(MyEventAction*);
+ ~MyTrackingAction() override;
- virtual void PreUserTrackingAction(const G4Track*);
- virtual void PostUserTrackingAction(const G4Track*);
+ void PreUserTrackingAction(const G4Track*) override;
private:
-
+
MyEventAction *fEventAction;
};
diff --git a/FullSimLight/include/PythiaPrimaryGeneratorAction.hh b/FullSimLight/include/PythiaPrimaryGeneratorAction.hh
index 3ab2b3be06d8e924c0c5841187f0a7193a843392..4156103a4f2372521d5b1332f268806bc5433157 100644
--- a/FullSimLight/include/PythiaPrimaryGeneratorAction.hh
+++ b/FullSimLight/include/PythiaPrimaryGeneratorAction.hh
@@ -10,16 +10,22 @@ class G4Event;
bool use_pythia();
void set_pythia_config(const char*);
+const char* get_pythia_config();
class PythiaPrimaryGeneratorAction final : public G4VUserPrimaryGeneratorAction {
public:
- PythiaPrimaryGeneratorAction();
+ // seeding < 0 : fixed 1234 for all threads (used for performance measurements)
+ // seeding = 0 : re-seeding by the event ID in order to guarantee reproducibilty
+ PythiaPrimaryGeneratorAction(int seeding = -1);
void GeneratePrimaries(G4Event*) override;
private:
- Pythia8::Pythia pythia;
+
+ int fSeeding;
+ Pythia8::Pythia fPythia;
+
};
#endif
diff --git a/FullSimLight/src/MyActionInitialization.cc b/FullSimLight/src/MyActionInitialization.cc
index 95d6701660a53ee18c434e1bbd6f02f67548c244..3bcd4ed7680e6ca494230db0a914fac65fb84d98 100644
--- a/FullSimLight/src/MyActionInitialization.cc
+++ b/FullSimLight/src/MyActionInitialization.cc
@@ -31,6 +31,12 @@ MyActionInitialization::~MyActionInitialization() {}
void MyActionInitialization::BuildForMaster() const {
MyRunAction* masterRunAct = new MyRunAction(fCreateGeantinoMaps,fGeantinoMapsFilename);
masterRunAct->SetPerformanceFlag(fIsPerformance);
+#if USE_PYTHIA
+ if (use_pythia()) {
+ G4String str(get_pythia_config());
+ masterRunAct->SetPythiaConfig(str);
+ }
+#endif
SetUserAction(masterRunAct);
}
@@ -40,10 +46,14 @@ void MyActionInitialization::Build() const {
#if !USE_PYTHIA
SetUserAction(new MyPrimaryGeneratorAction());
#else
- if (use_pythia())
- SetUserAction(new PythiaPrimaryGeneratorAction());
- else
+ if (use_pythia()) {
+ // seed each generator/thread by 1234 if perfomance mode run and use the event
+ // ID+1 as seed otherwise (guaranted reproducibility while having different events)
+ G4int pythiaSeed = fIsPerformance ? -1 : 0;
+ SetUserAction(new PythiaPrimaryGeneratorAction(pythiaSeed));
+ } else {
SetUserAction(new MyPrimaryGeneratorAction());
+ }
#endif
#ifndef G4MULTITHREADED
@@ -52,6 +62,12 @@ void MyActionInitialization::Build() const {
if (fIsPerformance) {
MyRunAction* masterRunAct = new MyRunAction(fCreateGeantinoMaps, fGeantinoMapsFilename);
masterRunAct->SetPerformanceFlag(fIsPerformance);
+#if USE_PYTHIA
+ if (use_pythia()) {
+ G4String str(get_pythia_config());
+ masterRunAct->SetPythiaConfig(str);
+ }
+#endif
SetUserAction(masterRunAct);
}
#endif
@@ -59,14 +75,14 @@ void MyActionInitialization::Build() const {
if (!fIsPerformance) {
MyRunAction* runact = new MyRunAction(fCreateGeantinoMaps, fGeantinoMapsFilename);
SetUserAction(runact);
-
+
if(!fCreateGeantinoMaps){
MyEventAction* evtact = new MyEventAction();
SetUserAction(evtact);
SetUserAction(new MyTrackingAction(evtact));
SetUserAction(new MySteppingAction(evtact));
-
}
+
#if G4VERSION_NUMBER>=1040
else
{
@@ -85,7 +101,7 @@ void MyActionInitialization::Build() const {
myLenghtIntEventAct->SetCreateEtaPhiMaps(fCreateEtaPhiMaps);
SetUserAction(myLenghtIntEventAct);
SetUserAction(myLenghtIntSteppingAct);
-
+
}
#endif
//MultiEventActions?? TO DO?
diff --git a/FullSimLight/src/MyEventAction.cc b/FullSimLight/src/MyEventAction.cc
index 343a127a30b3be30a633c0f52c20fcf639ae0c06..275095007963d8c508f3a11ba23cec130b6ccc46 100644
--- a/FullSimLight/src/MyEventAction.cc
+++ b/FullSimLight/src/MyEventAction.cc
@@ -23,77 +23,46 @@
MyEventAction::MyEventAction() : G4UserEventAction() {
- fNumberOfPrimariesPerEvent = 10; // will be changed automatically if needed
- fEventDataPerPrimary.resize(fNumberOfPrimariesPerEvent);
- for (G4int id=0; id<fNumberOfPrimariesPerEvent; ++id) {
- fEventDataPerPrimary[id].Clear();
- }
+ fEventData.Clear();
}
-MyEventAction::~MyEventAction() {
- fEventDataPerPrimary.clear();
-}
+MyEventAction::~MyEventAction() { }
-void MyEventAction::BeginOfEventAction(const G4Event* evt) {
- G4int curNumPrims = evt->GetNumberOfPrimaryVertex();
- if (curNumPrims>fNumberOfPrimariesPerEvent) {
- fNumberOfPrimariesPerEvent = curNumPrims;
- fEventDataPerPrimary.resize(fNumberOfPrimariesPerEvent);
- }
- for (G4int id=0; id<curNumPrims; ++id) {
- fEventDataPerPrimary[id].Clear();
- }
+void MyEventAction::BeginOfEventAction(const G4Event*) {
+ fEventData.Clear();
}
-void MyEventAction::EndOfEventAction(const G4Event* evt) {
+void MyEventAction::EndOfEventAction(const G4Event*) {
//get the Run and add the data collected during the simulation of the event that has been completed
MyRun* run = static_cast<MyRun*>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());
- G4int nPrimPart = evt->GetNumberOfPrimaryVertex();
- G4cout << " \n================================================================= \n"
- << " === EndOfEventAction --- event = " << evt->GetEventID() << " with "<< nPrimPart << " primary:"<<G4endl;
- for (G4int ip=0; ip<nPrimPart; ++ip) {
- G4PrimaryParticle* pPart = evt->GetPrimaryVertex(ip)->GetPrimary(0);
- G4String pPartName = pPart->GetParticleDefinition()->GetParticleName();
- G4double pPartEkin = pPart->GetKineticEnergy();
- const G4ThreeVector& pPartDir = pPart->GetMomentumDirection();
- G4int primaryTypeIndx = MyPrimaryGeneratorAction::GetPrimaryTypeIndex(pPartName);
- run->FillPerEventNumPrimaries(primaryTypeIndx);
- run->FillPerEvent(fEventDataPerPrimary[ip],primaryTypeIndx);
- G4cout<< " Primary Particle: " << ip << " (type index = " << primaryTypeIndx << ")\n"
- << " Name = " << pPartName << " \n"
- << " Energy = " << pPartEkin/GeV << " [GeV]\n"
- << " Direction = " << pPartDir << " \n"
- << fEventDataPerPrimary[ip]<< G4endl;
- }
+ run->FillPerEvent(fEventData);
}
-void MyEventAction::AddData(G4double edep, G4double length, G4bool ischarged, G4int primid) {
- fEventDataPerPrimary[primid].fEdep += edep;
+// called from the SteppingAction
+void MyEventAction::AddData(G4double edep, G4double length, G4bool ischarged) {
+ fEventData.fEdep += edep;
if (ischarged) {
- fEventDataPerPrimary[primid].fTrackLCh += length;
- fEventDataPerPrimary[primid].fChargedStep += 1.;
+ fEventData.fTrackLCh += length;
+ fEventData.fChargedStep += 1.;
} else {
- fEventDataPerPrimary[primid].fTrackLNe += length;
- fEventDataPerPrimary[primid].fNeutralStep += 1.;
+ fEventData.fTrackLNe += length;
+ fEventData.fNeutralStep += 1.;
}
}
-void MyEventAction::AddSecondaryTrack(const G4Track* track, G4int primid) {
+// called from the TrackingAction
+void MyEventAction::AddSecondaryTrack(const G4Track* track) {
const G4ParticleDefinition* pdf = track->GetDefinition();
if (pdf==G4Gamma::Gamma()) {
- fEventDataPerPrimary[primid].fNGamma += 1.;
-
+ fEventData.fNGamma += 1.;
} else if (pdf==G4Electron::Electron()) {
- fEventDataPerPrimary[primid].fNElec += 1.;
-
+ fEventData.fNElec += 1.;
} else if (pdf==G4Positron::Positron()) {
- fEventDataPerPrimary[primid].fNPosit += 1.;
-
+ fEventData.fNPosit += 1.;
}
-
}
diff --git a/FullSimLight/src/MyEventDataPerPrimary.cc b/FullSimLight/src/MyEventData.cc
similarity index 72%
rename from FullSimLight/src/MyEventDataPerPrimary.cc
rename to FullSimLight/src/MyEventData.cc
index 9a2fa12601e1b17b9333e3e47822a1b5a19ef30d..658148ecc6ab685423f2ebc03b190b129ac1695f 100644
--- a/FullSimLight/src/MyEventDataPerPrimary.cc
+++ b/FullSimLight/src/MyEventData.cc
@@ -1,24 +1,17 @@
-#include "MyEventDataPerPrimary.hh"
+#include "MyEventData.hh"
#include "G4SystemOfUnits.hh"
-MyEventDataPerPrimary::MyEventDataPerPrimary() {
- fEdep = 0.;
- fTrackLCh = 0.;
- fTrackLNe = 0.;
- fChargedStep = 0.;
- fNeutralStep = 0.;
- fNGamma = 0.;
- fNElec = 0.;
- fNPosit = 0.;
+MyEventData::MyEventData() {
+ Clear();
}
-MyEventDataPerPrimary::~MyEventDataPerPrimary() {}
+MyEventData::~MyEventData() {}
-void MyEventDataPerPrimary::Clear() {
+void MyEventData::Clear() {
fEdep = 0.;
fTrackLCh = 0.;
fTrackLNe = 0.;
@@ -30,7 +23,7 @@ void MyEventDataPerPrimary::Clear() {
}
-std::ostream& operator<<(std::ostream& flux, const MyEventDataPerPrimary& evtdata) {
+std::ostream& operator<<(std::ostream& flux, const MyEventData& evtdata) {
std::ios::fmtflags mode = flux.flags();
flux.setf(std::ios::fixed,std::ios::floatfield);
long prec = flux.precision(3);
diff --git a/FullSimLight/src/MyRun.cc b/FullSimLight/src/MyRun.cc
index 69825892da4ba8b8647711bff233d7c7d2d2dbe4..15a6d33e3765c5a38bb4dee08383a27230da90c9 100644
--- a/FullSimLight/src/MyRun.cc
+++ b/FullSimLight/src/MyRun.cc
@@ -1,109 +1,89 @@
#include "MyRun.hh"
-#include "MyEventDataPerPrimary.hh"
-#include "MyPrimaryGeneratorAction.hh"
+#include "MyEventData.hh"
#include "G4SystemOfUnits.hh"
#include <iomanip>
MyRun::MyRun() : G4Run() {
- fNumPrimaryTypes = MyPrimaryGeneratorAction::GetNumberOfPrimaryTypes();
- fRunDataPerPrimary.resize(fNumPrimaryTypes);
- for (G4int ipt=0; ipt<fNumPrimaryTypes; ++ipt) {
- fRunDataPerPrimary[ipt].Clear();
- }
+ fRunData.Clear();
}
-MyRun::~MyRun() {
- fRunDataPerPrimary.clear();
-}
-
+MyRun::~MyRun() {}
-void MyRun::FillPerEvent(const MyEventDataPerPrimary& data, G4int primtypeindx) {
- MyRunDataPerPrimary& runData = fRunDataPerPrimary[primtypeindx];
- runData.fEdep += data.fEdep; runData.fEdep2 += data.fEdep*data.fEdep;
- runData.fTrackLCh += data.fTrackLCh; runData.fTrackLCh2 += data.fTrackLCh*data.fTrackLCh;
- runData.fTrackLNe += data.fTrackLNe; runData.fTrackLNe2 += data.fTrackLNe*data.fTrackLNe;
- runData.fChargedStep += data.fChargedStep; runData.fChargedStep2 += data.fChargedStep*data.fChargedStep;
- runData.fNeutralStep += data.fNeutralStep; runData.fNeutralStep2 += data.fNeutralStep*data.fNeutralStep;
- runData.fNGamma += data.fNGamma; runData.fNGamma2 += data.fNGamma*data.fNGamma;
- runData.fNElec += data.fNElec; runData.fNElec2 += data.fNElec*data.fNElec;
- runData.fNPosit += data.fNPosit; runData.fNPosit2 += data.fNPosit*data.fNPosit;
+void MyRun::FillPerEvent(const MyEventData& data) {
+ fRunData.fEdep += data.fEdep; fRunData.fEdep2 += data.fEdep*data.fEdep;
+ fRunData.fTrackLCh += data.fTrackLCh; fRunData.fTrackLCh2 += data.fTrackLCh*data.fTrackLCh;
+ fRunData.fTrackLNe += data.fTrackLNe; fRunData.fTrackLNe2 += data.fTrackLNe*data.fTrackLNe;
+ fRunData.fChargedStep += data.fChargedStep; fRunData.fChargedStep2 += data.fChargedStep*data.fChargedStep;
+ fRunData.fNeutralStep += data.fNeutralStep; fRunData.fNeutralStep2 += data.fNeutralStep*data.fNeutralStep;
+ fRunData.fNGamma += data.fNGamma; fRunData.fNGamma2 += data.fNGamma*data.fNGamma;
+ fRunData.fNElec += data.fNElec; fRunData.fNElec2 += data.fNElec*data.fNElec;
+ fRunData.fNPosit += data.fNPosit; fRunData.fNPosit2 += data.fNPosit*data.fNPosit;
}
void MyRun::Merge(const G4Run* run) {
const MyRun* localRun = static_cast<const MyRun*>(run);
if (localRun) {
- for (G4int ipt=0; ipt<fNumPrimaryTypes; ++ipt) {
- fRunDataPerPrimary[ipt] += localRun->GetRunDataPerPrimary(ipt);
- }
+ fRunData += localRun->GetRunData();
}
G4Run::Merge(run);
}
void MyRun::EndOfRun() {
- G4int numEvents = GetNumberOfEvent();
- G4int numPrimaries = 0;
- for (G4int ipt=0; ipt<fNumPrimaryTypes; ++ipt) {
- numPrimaries += fRunDataPerPrimary[ipt].fNumPrimaries;
- }
+ const G4int numEvents = GetNumberOfEvent();
+ //
std::ios::fmtflags mode = G4cout.flags();
G4int prec = G4cout.precision(2);
G4cout<< " \n ================================== Run summary ===================================== \n" << G4endl;
G4cout<< std::setprecision(4);
- G4cout<< " Number of events = " << numEvents << G4endl;
- G4cout<< " Total number of primary = " << numPrimaries << G4endl;
+ G4cout<< " Number of events = " << numEvents << G4endl;
G4cout<< " \n ---------------------------------------------------------------------------------------- \n" << G4endl;
- // compute and print run statistics per primary type per primary
- for (G4int ipt=0; ipt<fNumPrimaryTypes; ++ipt) {
- MyRunDataPerPrimary& runData = fRunDataPerPrimary[ipt];
- G4int nPrimaries = runData.fNumPrimaries;
- G4double norm = static_cast<G4double>(nPrimaries);
- if (norm>0.) {
- norm = 1./norm;
- } else {
- continue;
- }
- //compute and print statistic
- //
- G4String primName = MyPrimaryGeneratorAction::GetPrimaryName(ipt);
- G4double meanEdep = runData.fEdep*norm;
- G4double rmsEdep = std::sqrt(std::abs(runData.fEdep2*norm-meanEdep*meanEdep));
- G4double meanLCh = runData.fTrackLCh*norm;
- G4double rmsLCh = std::sqrt(std::abs(runData.fTrackLCh2*norm-meanLCh*meanLCh));
- G4double meanLNe = runData.fTrackLNe*norm;
- G4double rmsLNe = std::sqrt(std::abs(runData.fTrackLNe2*norm-meanLNe*meanLNe));
- G4double meanStpCh = runData.fChargedStep*norm;
- G4double rmsStpCh = std::sqrt(std::abs(runData.fChargedStep2*norm-meanStpCh*meanStpCh));
- G4double meanStpNe = runData.fNeutralStep*norm;
- G4double rmsStpNe = std::sqrt(std::abs(runData.fNeutralStep2*norm-meanStpNe*meanStpNe));
- G4double meanNGam = runData.fNGamma*norm;
- G4double rmsNGam = std::sqrt(std::abs(runData.fNGamma2*norm-meanNGam*meanNGam));
- G4double meanNElec = runData.fNElec*norm;
- G4double rmsNElec = std::sqrt(std::abs(runData.fNElec2*norm-meanNElec*meanNElec));
- G4double meanNPos = runData.fNPosit*norm;
- G4double rmsNPos = std::sqrt(std::abs(runData.fNPosit2*norm-meanNPos*meanNPos));
-
- G4cout<< " Number of primaries = " << nPrimaries << " " << primName <<G4endl;
- G4cout<< " Total energy deposit per primary = " << meanEdep/GeV << " +- " << rmsEdep/GeV << " [GeV]" <<G4endl;
- G4cout<< G4endl;
- G4cout<< " Total track length (charged) per primary = " << meanLCh/cm << " +- " << rmsLCh/cm << " [cm]" <<G4endl;
- G4cout<< " Total track length (neutral) per primary = " << meanLNe/cm << " +- " << rmsLNe/cm << " [cm]" <<G4endl;
- G4cout<< G4endl;
- G4cout<< " Number of steps (charged) per primary = " << meanStpCh << " +- " << rmsStpCh << G4endl;
- G4cout<< " Number of steps (neutral) per primary = " << meanStpNe << " +- " << rmsStpNe << G4endl;
- G4cout<< G4endl;
- G4cout<< " Number of secondaries per primary : " << G4endl
- << " Gammas = " << meanNGam << " +- " << rmsNGam << G4endl
- << " Electrons = " << meanNElec << " +- " << rmsNElec << G4endl
- << " Positrons = " << meanNPos << " +- " << rmsNPos << G4endl;
- G4cout<< " ......................................................................................... \n" << G4endl;
+ // compute and print run statistics per event
+ G4double norm = numEvents;
+ if (norm>0.) {
+ norm = 1./norm;
+ } else {
+ return;
}
+ //compute and print statistic
+ //
+ const G4double meanEdep = fRunData.fEdep*norm;
+ const G4double rmsEdep = std::sqrt(std::abs(fRunData.fEdep2*norm-meanEdep*meanEdep));
+ const G4double meanLCh = fRunData.fTrackLCh*norm;
+ const G4double rmsLCh = std::sqrt(std::abs(fRunData.fTrackLCh2*norm-meanLCh*meanLCh));
+ const G4double meanLNe = fRunData.fTrackLNe*norm;
+ const G4double rmsLNe = std::sqrt(std::abs(fRunData.fTrackLNe2*norm-meanLNe*meanLNe));
+ const G4double meanStpCh = fRunData.fChargedStep*norm;
+ const G4double rmsStpCh = std::sqrt(std::abs(fRunData.fChargedStep2*norm-meanStpCh*meanStpCh));
+ const G4double meanStpNe = fRunData.fNeutralStep*norm;
+ const G4double rmsStpNe = std::sqrt(std::abs(fRunData.fNeutralStep2*norm-meanStpNe*meanStpNe));
+ const G4double meanNGam = fRunData.fNGamma*norm;
+ const G4double rmsNGam = std::sqrt(std::abs(fRunData.fNGamma2*norm-meanNGam*meanNGam));
+ const G4double meanNElec = fRunData.fNElec*norm;
+ const G4double rmsNElec = std::sqrt(std::abs(fRunData.fNElec2*norm-meanNElec*meanNElec));
+ const G4double meanNPos = fRunData.fNPosit*norm;
+ const G4double rmsNPos = std::sqrt(std::abs(fRunData.fNPosit2*norm-meanNPos*meanNPos));
+
+ G4cout<< " Mean energy deposit per event = " << meanEdep/GeV << " +- " << rmsEdep/GeV << " [GeV]" <<G4endl;
+ G4cout<< G4endl;
+ G4cout<< " Mean track length (charged) per event = " << meanLCh/cm << " +- " << rmsLCh/cm << " [cm]" <<G4endl;
+ G4cout<< " Mean track length (neutral) per event = " << meanLNe/cm << " +- " << rmsLNe/cm << " [cm]" <<G4endl;
+ G4cout<< G4endl;
+ G4cout<< " Number of steps (charged) per event = " << meanStpCh << " +- " << rmsStpCh << G4endl;
+ G4cout<< " Number of steps (neutral) per event = " << meanStpNe << " +- " << rmsStpNe << G4endl;
+ G4cout<< G4endl;
+ G4cout<< " Number of secondaries per event : " << G4endl
+ << " Gammas = " << meanNGam << " +- " << rmsNGam << G4endl
+ << " Electrons = " << meanNElec << " +- " << rmsNElec << G4endl
+ << " Positrons = " << meanNPos << " +- " << rmsNPos << G4endl;
+ G4cout<< " ......................................................................................... \n" << G4endl;
+
G4cout<< " \n ======================================================================================== \n" << G4endl;
G4cout.setf(mode,std::ios::floatfield);
diff --git a/FullSimLight/src/MyRunAction.cc b/FullSimLight/src/MyRunAction.cc
index 135108acf1690945a700191e88787b8e76ec192f..be9b01f4b4aa70db38434848b438ff6ced1dc475 100644
--- a/FullSimLight/src/MyRunAction.cc
+++ b/FullSimLight/src/MyRunAction.cc
@@ -11,7 +11,7 @@
#include "globals.hh"
#include "MyRun.hh"
-#include "MyPrimaryGeneratorAction.hh"
+
#include "G4ProductionCutsTable.hh"
#include "G4Region.hh"
@@ -19,9 +19,10 @@
G4AnalysisManager* MyRunAction::fMasterAnalysisManager = nullptr;
-MyRunAction::MyRunAction(bool isGeantino, G4String geantinoMapFilename) : G4UserRunAction(), fIsPerformance(false), fIsGeantino(isGeantino), fRun(nullptr), fTimer(nullptr), fGeantinoMapsFilename(geantinoMapFilename){
-
-}
+MyRunAction::MyRunAction(bool isGeantino, G4String geantinoMapFilename)
+: G4UserRunAction(), fIsPerformance(false), fIsGeantino(isGeantino),
+ fRun(nullptr), fTimer(nullptr), fGeantinoMapsFilename(geantinoMapFilename),
+ fPythiaConfig("") { }
MyRunAction::~MyRunAction() {
if(fIsGeantino)
@@ -42,7 +43,7 @@ G4Run* MyRunAction::GenerateRun() {
void MyRunAction::BeginOfRunAction(const G4Run* /*aRun*/){
#if G4VERSION_NUMBER>=1040
-
+
if(fIsGeantino)
{
// Create analysis manager
@@ -50,14 +51,14 @@ 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;
-
+
G4cout << "Using G4AnalysisManager type: " << analysisManager->GetType() << G4endl;
analysisManager->SetVerboseLevel(1);
-
+
// Open output root file
G4cout<<"\n\nBeginOfRunAction: \n...create a root file using the G4AnalysisManager" << G4endl;
if (!analysisManager->OpenFile(fGeantinoMapsFilename)){
@@ -65,7 +66,7 @@ void MyRunAction::BeginOfRunAction(const G4Run* /*aRun*/){
exit(-1);
} else
G4cout<<"\n...output File "<<fGeantinoMapsFilename<<" opened!"<<G4endl;
-
+
const char* radName = "RZRadLen";
if(analysisManager->GetP2Id(radName, false) < 0){
fRadName_id = analysisManager->CreateP2(radName,radName,1000,-25000.,25000.,2000,0.,15000.);
@@ -73,7 +74,7 @@ void MyRunAction::BeginOfRunAction(const G4Run* /*aRun*/){
analysisManager->SetP2XAxisTitle(fRadName_id,"Z[mm]");
analysisManager->SetP2YAxisTitle(fRadName_id,"R[mm]");
analysisManager->SetP2ZAxisTitle(fRadName_id,"thickstepRL");
-
+
}
const char* intName = "RZIntLen";
if(analysisManager->GetP2Id(intName, false)< 0)
@@ -90,7 +91,7 @@ void MyRunAction::BeginOfRunAction(const G4Run* /*aRun*/){
fEtaRad_id = analysisManager->CreateP1(etaRadName,etaRadName,500,-6,6);
analysisManager->SetP1XAxisTitle(fEtaRad_id,"#eta");
analysisManager->SetP1YAxisTitle(fEtaRad_id,"Radiation Length %X0");
-
+
}
const char* etaIntName = "etaIntLen";
if(analysisManager->GetP1Id(etaIntName, false)< 0)
@@ -99,13 +100,13 @@ void MyRunAction::BeginOfRunAction(const G4Run* /*aRun*/){
analysisManager->SetP1XAxisTitle(fEtaInt_id,"#eta");
analysisManager->SetP1YAxisTitle(fEtaInt_id,"Interaction Length #lambda");
}
-
+
}
#endif
if (isMaster) {
-
+
//G4cout<<"\nBeginOfRunAction isMaster, and fMasterAnalysisManager: "<<fMasterAnalysisManager<<G4endl;
-
+
std::vector<G4Region*>* regionVect = G4RegionStore::GetInstance();
int numRegions = regionVect->size();
int sumNumMC = 0;
@@ -127,17 +128,20 @@ void MyRunAction::BeginOfRunAction(const G4Run* /*aRun*/){
sumNumMC += mcRVect[ir];
}
G4cout<< " === Total number of MC = " << sumNumMC << " vs " << numMC << " #regions = " << numRegions << G4endl;
-
-
+
+
+ G4String msg0 = (fPythiaConfig == "") ? "a Particle-Gun" : "PYTHIA (" + fPythiaConfig +")";
#ifdef G4MULTITHREADED
G4int nThreads = G4MTRunManager::GetMasterRunManager()->GetNumberOfThreads();
+ G4int nEvent = G4MTRunManager::GetMasterRunManager()->GetNumberOfEventsToBeProcessed();
G4cout << "\n =======================================================================================\n"
- << " Run started in MT mode with " << nThreads << " threads \n"
+ << " Run started in MT mode with " << nThreads << " threads, simulating " << nEvent << " events using " << msg0 << " \n"
<< " ======================================================================================= \n"
<< G4endl;
#else
+ G4int nEvent = G4RunManager::GetRunManager()->GetNumberOfEventsToBeProcessed();
G4cout << "\n =======================================================================================\n"
- << " Run started in sequential mode (i.e. with 1 thread) \n"
+ << " Run started in sequential mode (i.e. with 1 thread), simulating " << nEvent << " events using " << msg0 << " \n"
<< " ======================================================================================= \n"
<< G4endl;
#endif
@@ -153,7 +157,7 @@ void MyRunAction::EndOfRunAction(const G4Run*) {
#if G4VERSION_NUMBER>=1040
if(fIsGeantino){
-
+
auto analysisManager= G4AnalysisManager::Instance();
//Finalize analysisManager and Write out file
if (analysisManager->IsOpenFile()){
@@ -163,10 +167,10 @@ void MyRunAction::EndOfRunAction(const G4Run*) {
analysisManager->CloseFile();
G4cout<<"Output File successfully saved and closed! " << G4endl;
}
-
+
}
#endif
-
+
if (isMaster) {
fTimer->Stop();
// get number of events: even in case of perfomance mode when MyRun-s are not generated in GenerateRun()
@@ -183,8 +187,7 @@ void MyRunAction::EndOfRunAction(const G4Run*) {
G4cout << " Time: " << *fTimer << G4endl;
G4cout << " ======================================================================================= "<< G4endl;
delete fTimer;
- // print primary gun properties (not available at the begining of the run)
- MyPrimaryGeneratorAction::Print();
+
if (!fIsPerformance) { // otherwise we do not even create any MyRun objects so fRun is nullptr
fRun->EndOfRun();
}
diff --git a/FullSimLight/src/MyRunDataPerPrimary.cc b/FullSimLight/src/MyRunData.cc
similarity index 75%
rename from FullSimLight/src/MyRunDataPerPrimary.cc
rename to FullSimLight/src/MyRunData.cc
index d50d283e326accdf0fdfe078cc66992f7729191f..1030bd9929b8c0f4d053790019e192d18a37baf2 100644
--- a/FullSimLight/src/MyRunDataPerPrimary.cc
+++ b/FullSimLight/src/MyRunData.cc
@@ -1,17 +1,16 @@
-#include "MyRunDataPerPrimary.hh"
+#include "MyRunData.hh"
-MyRunDataPerPrimary::MyRunDataPerPrimary() {
+MyRunData::MyRunData() {
Clear();
}
-MyRunDataPerPrimary::~MyRunDataPerPrimary() {}
+MyRunData::~MyRunData() {}
-void MyRunDataPerPrimary::Clear() {
- fNumPrimaries = 0.;
+void MyRunData::Clear() {
fEdep = 0.;
fEdep2 = 0.;
fTrackLCh = 0.;
@@ -33,8 +32,7 @@ void MyRunDataPerPrimary::Clear() {
}
-MyRunDataPerPrimary& MyRunDataPerPrimary::operator+=(const MyRunDataPerPrimary& right) {
- fNumPrimaries += right.fNumPrimaries;
+MyRunData& MyRunData::operator+=(const MyRunData& right) {
fEdep += right.fEdep;
fEdep2 += right.fEdep2;
fTrackLCh += right.fTrackLCh;
diff --git a/FullSimLight/src/MySteppingAction.cc b/FullSimLight/src/MySteppingAction.cc
index 6a46615c84e93f149fd77e952eaed530662d9e63..e33e45022bbd9d45932662dfe9c16ba7a9426f09 100644
--- a/FullSimLight/src/MySteppingAction.cc
+++ b/FullSimLight/src/MySteppingAction.cc
@@ -2,7 +2,6 @@
#include "MySteppingAction.hh"
#include "MyEventAction.hh"
-#include "MyTrackInformation.hh"
#include "G4Step.hh"
#include "G4ParticleDefinition.hh"
@@ -16,8 +15,6 @@ MySteppingAction::~MySteppingAction() {}
void MySteppingAction::UserSteppingAction(const G4Step* theStep) {
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);
+ G4bool isCharged = (theStep->GetTrack()->GetDefinition()->GetPDGCharge() != 0.);
+ fEventAction->AddData(edep, stepl, isCharged);
}
diff --git a/FullSimLight/src/MyTrackInformation.cc b/FullSimLight/src/MyTrackInformation.cc
deleted file mode 100644
index 4e2ec6361dd3fb544e80d8dc8954d372d8a328bb..0000000000000000000000000000000000000000
--- a/FullSimLight/src/MyTrackInformation.cc
+++ /dev/null
@@ -1,6 +0,0 @@
-
-#include "MyTrackInformation.hh"
-
-MyTrackInformation::MyTrackInformation(G4int id) : fPrimaryTrackID(id) {}
-
-MyTrackInformation::~MyTrackInformation() {}
diff --git a/FullSimLight/src/MyTrackingAction.cc b/FullSimLight/src/MyTrackingAction.cc
index 33196227d369b101eafe6d85fe415024731593cb..9e3c591ac9354281061086a61afff2d9dd1f6e62 100644
--- a/FullSimLight/src/MyTrackingAction.cc
+++ b/FullSimLight/src/MyTrackingAction.cc
@@ -2,10 +2,7 @@
#include "MyTrackingAction.hh"
#include "G4Track.hh"
-#include "G4TrackVector.hh"
-#include "G4TrackingManager.hh"
#include "MyEventAction.hh"
-#include "MyTrackInformation.hh"
MyTrackingAction::MyTrackingAction(MyEventAction* evtact) : G4UserTrackingAction(), fEventAction(evtact) {}
@@ -15,28 +12,8 @@ MyTrackingAction::~MyTrackingAction() {}
void MyTrackingAction::PreUserTrackingAction(const G4Track* theTrack) {
- if (theTrack->GetParentID()==0) { // primary track so set its auxiliary track member: primary trackID
- theTrack->SetUserInformation(new MyTrackInformation(theTrack->GetTrackID()));
- } else {
- G4int primTrackID = static_cast<MyTrackInformation*>(theTrack->GetUserInformation())->GetPrimaryTrackID();
- fEventAction->AddSecondaryTrack(theTrack,primTrackID-1);
- }
-}
-
-
-void MyTrackingAction::PostUserTrackingAction(const G4Track* theTrack) {
- G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries();
- if (secondaries) {
- MyTrackInformation* info = static_cast<MyTrackInformation*>(theTrack->GetUserInformation());
- G4int primaryTrackID = info->GetPrimaryTrackID(); // get the primary trackID
- G4int currentTrackID = theTrack->GetTrackID(); // track ID of theTrack
- size_t nSecondaries = secondaries->size();
- for (size_t isec=0; isec<nSecondaries; ++isec) {
- // if the secondary was created by theTrack set its auxiliary info to store the primaryTrackID
- G4Track* secTrack = (*secondaries)[isec];
- if (secTrack->GetParentID()==currentTrackID) {
- secTrack->SetUserInformation(new MyTrackInformation(primaryTrackID));
- }
- }
+ // count secondaries
+ if (theTrack->GetParentID()!=0) {
+ fEventAction->AddSecondaryTrack(theTrack);
}
}
diff --git a/FullSimLight/src/PythiaPrimaryGeneratorAction.cc b/FullSimLight/src/PythiaPrimaryGeneratorAction.cc
index 37aecd0ab83bf294d3e7562f712a355e6c265560..0fba8446df1a7bc80e0b81aa2903cb7c8a2f3ce3 100644
--- a/FullSimLight/src/PythiaPrimaryGeneratorAction.cc
+++ b/FullSimLight/src/PythiaPrimaryGeneratorAction.cc
@@ -13,105 +13,110 @@
static const char *config = nullptr;
static const char *pythia_defaults[] = {
- "Print:quiet = on",
- "Beams:idA = 2212",
- "Beams:idB = 2212",
- "Beams:eCM = 13000.0",
- "Init:showProcesses = off",
- "Init:showMultipartonInteractions = off",
- "Init:showChangedParticleData = off",
- "Stat:showProcessLevel = off",
- "Stat:showErrors = off",
+ "Print:quiet = on",
+ "Beams:idA = 2212",
+ "Beams:idB = 2212",
+ "Beams:eCM = 13000.0",
+ "Init:showProcesses = off",
+ "Init:showMultipartonInteractions = off",
+ "Init:showChangedParticleData = off",
+ "Stat:showProcessLevel = off",
+ "Stat:showErrors = off",
};
static const char *pythia_minbias[] = {
- "HardQCD:all = on",
- "Beams:eCM = 13000.0",
- "PhaseSpace:pTHatMin = 20.0",
+ "HardQCD:all = on",
+ "Beams:eCM = 13000.0",
+ "PhaseSpace:pTHatMin = 20.0",
};
static const char *pythia_higgs[] = {
- "HiggsSM:all = on",
- "Beams:eCM = 13000.0",
- "PhaseSpace:pTHatMin = 20.0",
+ "HiggsSM:all = on",
+ "Beams:eCM = 13000.0",
+ "PhaseSpace:pTHatMin = 20.0",
};
static const char *pythia_ttbar[] = {
- "Top:gg2ttbar = on",
- "Top:qqbar2ttbar = on",
- "Beams:eCM = 13000.0",
- "PhaseSpace:pTHatMin = 20.0",
+ "Top:gg2ttbar = on",
+ "Top:qqbar2ttbar = on",
+ "Beams:eCM = 13000.0",
+ "PhaseSpace:pTHatMin = 20.0",
};
-bool use_pythia()
-{
+bool use_pythia() {
return config != nullptr;
}
-void set_pythia_config(const char *cfg)
-{
+void set_pythia_config(const char *cfg) {
config = cfg;
}
-PythiaPrimaryGeneratorAction::PythiaPrimaryGeneratorAction() {
- if (access(config, R_OK) == 0) {
- pythia.readFile(config);
- } else {
- /*
- flag name="Pythia:setSeed" default="off"
- Indicates whether a user-set seed should be used every time the Pythia::init routine is called. If off, the random number generator is initialized with its default seed at the beginning of the run, and never again. If on, each new Pythia::init call (should several be made in the same run) results in the random number being re-initialized, thereby possibly starting over with the same sequence, if you do not watch out.
-
- mode name="Pythia:seed" default="-1" max="900000000"
- The seed to be used, if setSeed is on.
- A negative value gives the default seed,
- a value 0 gives a random seed based on the time, and
- a value between 1 and 900,000,000 a unique different random number sequence.
- */
- pythia.readString("Random:setSeed = on");
- // use a reproducible seed: always the same results for benchmarks.
- pythia.readString("Random:seed = 1234");
- for (const auto str : pythia_defaults)
- pythia.readString(str);
-
- if (strcmp("minbias", config) == 0) {
- for (const auto str : pythia_minbias)
- pythia.readString(str);
- } else if (strcmp("higgs", config) == 0) {
- for (const auto str : pythia_higgs)
- pythia.readString(str);
- } else if (strcmp("ttbar", config) == 0) {
- for (const auto str : pythia_ttbar)
- pythia.readString(str);
- } else {
- errx(1, "unknown Pythia configuration: %s", config);
- }
- }
-
- pythia.init();
+const char* get_pythia_config() { return config; }
+
+PythiaPrimaryGeneratorAction::PythiaPrimaryGeneratorAction(int seeding) : fSeeding(seeding), fPythia("../share/Pythia8/xmldoc", false) {
+ if (access(config, R_OK) == 0) {
+ fPythia.readFile(config);
+ } else {
+ /*
+ flag name="Pythia:setSeed" default="off"
+ Indicates whether a user-set seed should be used every time the Pythia::init
+ routine is called. If off, the random number generator is initialized with its
+ default seed at the beginning of the run, and never again. If on, each new
+ Pythia::init call (should several be made in the same run) results in the
+ random number being re-initialized, thereby possibly starting over with the
+ same sequence, if you do not watch out.
+
+ mode name="Pythia:seed" default="-1" max="900000000"
+ The seed to be used, if setSeed is on.
+ A negative value gives the default seed,
+ a value 0 gives a random seed based on the time, and
+ a value between 1 and 900,000,000 a unique different random number sequence.
+ */
+ fPythia.readString("Random:setSeed = on");
+ // use a reproducible seed: always the same results for benchmarks.
+ fPythia.readString("Random:seed = 1234");
+ //
+ for (const auto str : pythia_defaults) fPythia.readString(str);
+ if (strcmp("minbias", config) == 0) {
+ for (const auto str : pythia_minbias) fPythia.readString(str);
+ } else if (strcmp("higgs", config) == 0) {
+ for (const auto str : pythia_higgs) fPythia.readString(str);
+ } else if (strcmp("ttbar", config) == 0) {
+ for (const auto str : pythia_ttbar) fPythia.readString(str);
+ } else {
+ errx(1, "unknown Pythia configuration: %s", config);
+ }
+ }
+ fPythia.init();
}
-void PythiaPrimaryGeneratorAction::GeneratePrimaries(G4Event *event)
-{
- static const G4double time = 0.0;
- static const G4ThreeVector position(0.0, 0.0, 0.0);
-
- G4PrimaryVertex* vertex = new G4PrimaryVertex(position, time);
-
- pythia.next();
-
- for (auto i = 1, n = pythia.event.size(); i < n; ++i) {
- const auto& particle = pythia.event[i];
-
- if (!particle.isFinal())
- continue;
-
- G4PrimaryParticle* p = new G4PrimaryParticle(particle.id());
- p->SetMass(particle.m() * GeV);
- p->SetMomentum(particle.px() * GeV, particle.py() * GeV, particle.pz() * GeV);
- vertex->SetPrimary(p);
- }
-
- event->AddPrimaryVertex(vertex);
+void PythiaPrimaryGeneratorAction::GeneratePrimaries(G4Event *event) {
+ static const G4double time = 0.0;
+ static const G4ThreeVector position(0.0, 0.0, 0.0);
+
+ G4PrimaryVertex* vertex = new G4PrimaryVertex(position, time);
+ //
+ // re-seed pythia with the event ID + 1 in order to ensure that the
+ // events are the same for a given event ID independently from any other
+ // run conditions
+ if (fSeeding == 0) {
+ fPythia.rndm.init(event->GetEventID()+1);
+ }
+ fPythia.next();
+
+ for (auto i = 1, n = fPythia.event.size(); i < n; ++i) {
+ const auto& particle = fPythia.event[i];
+
+ if (!particle.isFinal())
+ continue;
+
+ G4PrimaryParticle* p = new G4PrimaryParticle(particle.id());
+ p->SetMass(particle.m() * GeV);
+ p->SetMomentum(particle.px() * GeV, particle.py() * GeV, particle.pz() * GeV);
+ vertex->SetPrimary(p);
+ }
+
+ event->AddPrimaryVertex(vertex);
}
#endif