diff --git a/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/AthenaHepMCInterface.h b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/AthenaHepMCInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..44ff779143baa649e57611468e1b02222a509b1e
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/AthenaHepMCInterface.h
@@ -0,0 +1,27 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef AthenaHepMCInterface_H
+#define AthenaHepMCInterface_H
+
+#include "FadsKinematics/FadsGenerator.h"
+#include <string>
+
+class StoreGateSvc;
+class MsgStream;
+
+class AthenaHepMCInterface: public FADS::FadsGenerator {
+public:
+ AthenaHepMCInterface(std::string);
+ ~AthenaHepMCInterface();
+ void Initialize();
+ void Terminate();
+ HepMC::GenEvent* GenerateAnEvent();
+private:
+ StoreGateSvc *p_sgSvc;
+ MsgStream * m_log;
+ MsgStream log();
+};
+
+#endif
diff --git a/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/AthenaStackingAction.h b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/AthenaStackingAction.h
new file mode 100644
index 0000000000000000000000000000000000000000..17a376942129455c0389759b27d2b51945ad2142
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/AthenaStackingAction.h
@@ -0,0 +1,31 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef AthenaStackingAction_H
+#define AthenaStackingAction_H
+
+#include "FadsActions/ApplicationStackingAction.h"
+
+class AthenaStackingAction:public FADS::ApplicationStackingAction {
+public:
+ AthenaStackingAction();
+ ~AthenaStackingAction(){}
+
+ G4ClassificationOfNewTrack ClassifyNewTrack(const G4Track* aTrack);
+ void NewStage();
+ void PrepareNewEvent();
+
+ inline void KillAllNeutrinos(const bool a){ p_killAllNeutrinos=a; }
+ inline bool KillAllNeutrinos() const { return p_killAllNeutrinos; }
+
+ inline void KillLowEPhotons(const double a){ p_stackEnergyCut=a; }
+ inline double KillLowEPhotons() const { return p_stackEnergyCut; }
+
+private:
+
+ bool p_killAllNeutrinos;
+ double p_stackEnergyCut;
+};
+
+#endif
diff --git a/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/AthenaTrackingAction.h b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/AthenaTrackingAction.h
new file mode 100644
index 0000000000000000000000000000000000000000..689b8690593a6d7659694faa5c6434c655bf6b87
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/AthenaTrackingAction.h
@@ -0,0 +1,21 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef AthenaTrackingAction_H
+#define AthenaTrackingAction_H
+
+#include "FadsActions/ApplicationTrackingAction.h"
+
+class AthenaTrackingAction: public FADS::ApplicationTrackingAction {
+public:
+ AthenaTrackingAction(): FADS::ApplicationTrackingAction()
+ {
+ }
+ void PreUserTrackingAction(const G4Track*);
+ void PostUserTrackingAction(const G4Track*);
+private:
+
+};
+
+#endif
diff --git a/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/G4AtlasAlg.h b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/G4AtlasAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..2d320eec5434b19564e48f8974ed017f3e712da2
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/G4AtlasAlg.h
@@ -0,0 +1,63 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef G4AtlasAlg_H
+#define G4AtlasAlg_H
+
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "AthenaKernel/IAtRndmGenSvc.h"
+
+#include <string>
+#include <map>
+
+class AthenaStackingAction;
+class AthenaTrackingAction;
+class G4AtlasRunManager;
+
+template <class ConcreteAlgorithm> class AlgFactory;
+
+
+class G4AtlasAlg : public AthAlgorithm {
+
+ friend class AlgFactory<G4AtlasAlg>;
+
+public:
+
+ G4AtlasAlg(const std::string& name, ISvcLocator* pSvcLocator);
+ virtual ~G4AtlasAlg() { };
+
+ StatusCode initialize();
+ StatusCode finalize();
+ StatusCode execute();
+
+ void initializeG4();
+
+private:
+
+ G4AtlasRunManager* p_runMgr;
+
+ // Properties for the jobOptions
+ std::string libList;
+ std::string physList;
+ std::string generator;
+ std::string fieldMap;
+ std::string rndmGen;
+ bool m_KillAllNeutrinos;
+ double m_KillLowEPhotons;
+ bool m_releaseGeoModel;
+ bool m_IncludeParentsInG4Event;
+ bool m_killAbortedEvents;
+ bool m_flagAbortedEvents;
+
+ std::map<std::string,std::string> m_verbosities;
+
+ AthenaStackingAction* stackingAction;
+ AthenaTrackingAction* trackingAction;
+
+ // Random number service
+ ServiceHandle<IAtRndmGenSvc> m_rndmGenSvc;
+};
+
+#endif
diff --git a/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/G4AtlasAlgDict.h b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/G4AtlasAlgDict.h
new file mode 100644
index 0000000000000000000000000000000000000000..2e6c94340a5c19590a00595b3811c06acd3d2be6
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/G4AtlasAlgDict.h
@@ -0,0 +1,10 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef G4ATLASALG_G4ATLASALGDICT_H
+#define G4ATLASALG_G4ATLASALGDICT_H 1
+
+#include "G4AtlasAlg/G4AtlasAlg.h"
+
+#endif
diff --git a/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/G4AtlasRunManager.h b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/G4AtlasRunManager.h
new file mode 100644
index 0000000000000000000000000000000000000000..2b4bd50bbec227b8de8877ded75f4cd0edba8a98
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/G4AtlasRunManager.h
@@ -0,0 +1,60 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+
+#ifndef G4AtlasRunManager_h
+#define G4AtlasRunManager_h 1
+
+#include "G4RunManager.hh"
+#include "G4VUserPhysicsList.hh"
+#include "GaudiKernel/MsgStream.h"
+
+#include <vector>
+
+class StoreGateSvc;
+class MsgStream;
+
+class G4AtlasRunManager: public G4RunManager {
+
+ friend class G4AtlasAlg;
+ friend class SimControl;
+
+public:
+
+ virtual ~G4AtlasRunManager() { if (m_log) delete m_log; }
+
+ static G4AtlasRunManager* GetG4AtlasRunManager();
+
+ void SetPhysicsList(G4VUserPhysicsList* p) { m_pl = p; }
+ G4Event* GenerateEvent(G4int i_event);
+ bool SimulateFADSEvent();
+ void RunTermination();
+ void SetCurrentG4Event(int);
+
+protected:
+
+ void InitializeGeometry();
+ void InitializePhysics();
+
+private:
+
+ G4AtlasRunManager();
+ void EndEvent();
+
+ void SetStoreGatePtr( StoreGateSvc *sgs) { m_sgSvc = sgs; }
+ void SetReleaseGeo(bool b) { m_releaseGeo = b; }
+ void SetLogLevel(int i) { log().setLevel(i); }
+
+private:
+
+ G4VUserPhysicsList* m_pl;
+ StoreGateSvc* m_sgSvc;
+ bool m_releaseGeo;
+ MsgStream * m_log;
+ MsgStream log();
+};
+
+
+#endif
+
diff --git a/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/PreEventAction.h b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/PreEventAction.h
new file mode 100644
index 0000000000000000000000000000000000000000..bccfc9e6cb291d3a043478fd7255ed3ec93613ae
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/PreEventAction.h
@@ -0,0 +1,20 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef PreEventAction_H
+#define PreEventAction_H
+
+#include <string>
+
+class PreEventAction {
+private:
+ std::string name;
+public:
+ PreEventAction(std::string);
+ virtual ~PreEventAction() {}
+ std::string GetName();
+ virtual void Execute() = 0;
+};
+
+#endif
diff --git a/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/PreEventActionManager.h b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/PreEventActionManager.h
new file mode 100644
index 0000000000000000000000000000000000000000..e389a2b5358725304cc289895084c3b4a5516ee2
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/PreEventActionManager.h
@@ -0,0 +1,27 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef PreEventActionManager_H
+#define PreEventActionManager_H
+
+#include <string>
+#include <map>
+
+class PreEventAction;
+typedef std::map<std::string, PreEventAction*, std::less<std::string> > preEvActionMap;
+
+
+class PreEventActionManager {
+public:
+ static PreEventActionManager* GetPreEventActionManager();
+ void Execute();
+ void RegisterAction(PreEventAction *);
+private:
+ static PreEventActionManager* thePointer;
+ PreEventActionManager();
+ PreEventActionManager(const PreEventActionManager&) {}
+ preEvActionMap theMap;
+};
+
+#endif
diff --git a/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/TruthHepMCEventConverter.h b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/TruthHepMCEventConverter.h
new file mode 100644
index 0000000000000000000000000000000000000000..bf2a982e84dabe65169d2f4d8ddd346cf5411294
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/TruthHepMCEventConverter.h
@@ -0,0 +1,38 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef TruthHepMCEventConverter_H
+#define TruthHepMCEventConverter_H
+
+#include "FadsKinematics/HepMCEventConverter.h"
+#include <vector>
+
+class G4PrimaryParticle;
+class G4PrimaryVertex;
+
+// Helper class to convert HepMC into G4. this is used by the Generator
+// Center, instead of the standard conversion helper
+
+class TruthHepMCEventConverter: public FADS::HepMCEventConverter {
+public:
+ TruthHepMCEventConverter();
+ ~TruthHepMCEventConverter();
+ void HepMC2G4(const HepMC::GenEvent* ,G4Event *, const bool) const;
+private:
+ TruthHepMCEventConverter (const TruthHepMCEventConverter&); //declared, but not defined
+ TruthHepMCEventConverter& operator= (const TruthHepMCEventConverter&); //declared, but not defined
+ // helper functions here
+ bool SeparateVertices(HepMC::GenVertex* , HepMC::GenVertex*) const ;
+ G4PrimaryVertex* ConstructG4PrimaryVertex(HepMC::GenVertex*) const ;
+ G4PrimaryParticle* ConstructG4PrimaryParticle(HepMC::GenParticle*) const;
+ HepMC::GenVertex* ConstructHepMCVertex(G4PrimaryVertex*) const ;
+ void TransformHepMCParticle(HepMC::GenParticle*,G4PrimaryParticle*) const;
+ void TransformHepMCVertex(HepMC::GenVertex*,G4PrimaryVertex*) const;
+ void ModifyVertex(G4PrimaryParticle *pp, HepMC::GenVertex *v) const;
+ bool ValidHepMCVertex(HepMC::GenVertex*) const;
+ int printLev;
+ std::vector<int> *m_modded;
+};
+
+#endif
diff --git a/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/selection.xml b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/selection.xml
new file mode 100644
index 0000000000000000000000000000000000000000..3c3359d5948cdb79a3f459055eb8371551277a5a
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/G4AtlasAlg/selection.xml
@@ -0,0 +1,3 @@
+<lcgdict>
+ <class name="G4AtlasAlg"/>
+</lcgdict>
diff --git a/Simulation/G4Atlas/G4AtlasAlg/cmt/requirements b/Simulation/G4Atlas/G4AtlasAlg/cmt/requirements
new file mode 100644
index 0000000000000000000000000000000000000000..b2d530fba9a0710417e884827ec1de3f2b70de0d
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/cmt/requirements
@@ -0,0 +1,35 @@
+package G4AtlasAlg
+
+author Andrea DellAcqua <dellacqu@mail.cern.ch>
+
+branches src cmt doc
+
+use AtlasPolicy AtlasPolicy-*
+use AthenaBaseComps AthenaBaseComps-* Control
+use AthenaKernel AthenaKernel-* Control
+use GaudiInterface GaudiInterface-* External
+use Geant4 Geant4-* External
+use FadsActions FadsActions-* Simulation/G4Sim/FADS
+use FadsKinematics FadsKinematics-* Simulation/G4Sim/FADS
+
+private
+use AtlasHepMC AtlasHepMC-* External
+use MCTruth MCTruth-* Simulation/G4Sim
+use FadsGeometry FadsGeometry-* Simulation/G4Sim/FADS
+use FadsPhysics FadsPhysics-* Simulation/G4Sim/FADS
+use GeneratorObjects GeneratorObjects-* Generators
+use SimHelpers SimHelpers-* Simulation/G4Sim
+use StoreGate StoreGate-* Control
+use SGTools SGTools-* Control
+use GeoModelInterfaces GeoModelInterfaces-* DetectorDescription/GeoModel
+use EventInfo EventInfo-* Event
+end_private
+
+# build a dual_use library
+apply_pattern dual_use_library files=*.cxx
+
+#private
+#apply_pattern lcgdict dict=G4AtlasAlg \
+# headerfiles="../G4AtlasAlg/G4AtlasAlg.h" \
+# selectionfile=selection.xml
+#end_private
diff --git a/Simulation/G4Atlas/G4AtlasAlg/src/AthenaHepMCInterface.cxx b/Simulation/G4Atlas/G4AtlasAlg/src/AthenaHepMCInterface.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..3e01a8109190e8c570d83f6ed0fbc03d611d6473
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/src/AthenaHepMCInterface.cxx
@@ -0,0 +1,90 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "G4AtlasAlg/AthenaHepMCInterface.h"
+#include "G4AtlasAlg/TruthHepMCEventConverter.h"
+
+#include "HepMC/GenEvent.h"
+#include "FadsKinematics/FadsGeneratorT.h"
+#include "FadsKinematics/GeneratorCenter.h"
+
+#include "StoreGate/StoreGate.h"
+#include "StoreGate/DataHandle.h"
+
+// to permit access to StoreGate
+#include "GaudiKernel/ISvcLocator.h"
+#include "GaudiKernel/Bootstrap.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/IMessageSvc.h"
+#include "StoreGate/StoreGateSvc.h"
+
+#include "GeneratorObjects/McEventCollection.h"
+
+FADS::FadsGeneratorT<AthenaHepMCInterface> gen("AthenaHepMCInterface");
+
+AthenaHepMCInterface::AthenaHepMCInterface(std::string s):FADS::FadsGenerator(s),p_sgSvc(0),m_log(0)
+{
+ FADS::GeneratorCenter::GetGeneratorCenter()->DeleteHepMCEvent(false);
+ FADS::GeneratorCenter::GetGeneratorCenter()->SetHepMCEventConverter(
+ new TruthHepMCEventConverter);
+}
+AthenaHepMCInterface::~AthenaHepMCInterface()
+{
+ log() << MSG::DEBUG << "\b the AthenaHepMCInterface is being destroyed "<<"\n"
+ << " should you wish to use it again type: "<<"\n"
+ << " /Generator/Select AthenaHepMCInterface "<<"\n"<<endreq;
+ if (m_log) delete m_log;
+}
+void AthenaHepMCInterface::Initialize()
+{}
+void AthenaHepMCInterface::Terminate()
+{}
+HepMC::GenEvent* AthenaHepMCInterface::GenerateAnEvent()
+{
+ // Lazy initialization
+ if (!p_sgSvc){
+ log() << MSG::DEBUG <<
+ " standard interface to StoreGate for retrieving HepMC events"<<endreq;
+ ISvcLocator* svcLocator = Gaudi::svcLocator(); // from Bootstrap
+ StatusCode sc=svcLocator->service("StoreGateSvc",p_sgSvc);
+ if (sc.isFailure())
+ {
+ log() << MSG::WARNING <<"AthenaHepMCInterface could not access StoreGate!!"<<endreq;
+ }
+ }
+
+ const DataHandle<McEventCollection> mcCollptr;
+ StatusCode status=p_sgSvc->retrieve(mcCollptr,"GEN_EVENT");
+ if (status.isSuccess() )
+ {
+ if (0==mcCollptr->size()){
+ log() << MSG::WARNING << "No event in MC Event Collection!" << endreq;
+ return 0;
+ }
+
+ McEventCollection::const_iterator it = mcCollptr->begin();
+
+ // getting only the first event here.
+ HepMC::GenEvent *p_evt = (*it);
+ return p_evt;
+
+ }
+ else
+ {
+ log() << MSG::WARNING << "no McEventCollection found." << endreq;
+ }
+ return 0;
+}
+
+MsgStream AthenaHepMCInterface::log()
+{
+ if (m_log) return *m_log;
+ ISvcLocator* svcLocator = Gaudi::svcLocator();
+ IMessageSvc* p_msgSvc = 0;
+ if (svcLocator->service("MessageSvc", p_msgSvc).isFailure() || !p_msgSvc)
+ std::cout << "AthenaHepMCInterface: Trouble getting the message service. Should never happen. Will crash now." << std::endl;
+ m_log = new MsgStream(p_msgSvc,"AthenaHepMCInterface");
+ return *m_log;
+}
+
diff --git a/Simulation/G4Atlas/G4AtlasAlg/src/AthenaStackingAction.cxx b/Simulation/G4Atlas/G4AtlasAlg/src/AthenaStackingAction.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..56114c2853b9687b9fea6ec45e4d32b0a843fa35
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/src/AthenaStackingAction.cxx
@@ -0,0 +1,105 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "GaudiKernel/ISvcLocator.h"
+#include "GaudiKernel/Bootstrap.h"
+
+#include "G4AtlasAlg/AthenaStackingAction.h"
+#include "MCTruth/PrimaryParticleInformation.h"
+#include "MCTruth/TrackInformation.h"
+#include "MCTruth/TrackBarcodeInfo.h"
+#include "MCTruth/EventInformation.h"
+#include <iostream>
+#include <string>
+
+#include "G4VProcess.hh"
+#include "G4Event.hh"
+#include "G4EventManager.hh"
+
+#include "G4NeutrinoE.hh"
+#include "G4NeutrinoMu.hh"
+#include "G4NeutrinoTau.hh"
+#include "G4AntiNeutrinoE.hh"
+#include "G4AntiNeutrinoMu.hh"
+#include "G4AntiNeutrinoTau.hh"
+
+#include "G4Gamma.hh"
+
+// static AthenaStackingAction stacking;
+
+AthenaStackingAction::AthenaStackingAction(): ApplicationStackingAction(), p_killAllNeutrinos(false), p_stackEnergyCut(-1){}
+
+G4ClassificationOfNewTrack AthenaStackingAction::ClassifyNewTrack(const G4Track* aTrack)
+{
+ // Neutrinos turned off at job options level
+ G4ParticleDefinition &particleType = *(aTrack->GetDefinition() );
+ if ( p_killAllNeutrinos &&
+ ( &particleType == G4NeutrinoE::NeutrinoEDefinition() ||
+ &particleType == G4AntiNeutrinoE::AntiNeutrinoEDefinition() ||
+ &particleType == G4NeutrinoMu::NeutrinoMuDefinition() ||
+ &particleType == G4AntiNeutrinoMu::AntiNeutrinoMuDefinition() ||
+ &particleType == G4NeutrinoTau::NeutrinoTauDefinition() ||
+ &particleType == G4AntiNeutrinoTau::AntiNeutrinoTauDefinition() ) ){
+ return fKill;
+ }
+
+ if ( &particleType == G4Gamma::Gamma() &&
+ aTrack->GetTotalEnergy() < 0.00005 ){ // Safe cut for ultra-low-energy photons
+ return fKill;
+ }
+
+
+ G4Track *inT=const_cast<G4Track*>(aTrack);
+
+ G4Event *ev=G4EventManager::GetEventManager()->GetNonconstCurrentEvent();
+ EventInformation* eventInfo __attribute__ ((unused)) =static_cast<EventInformation*>
+ (ev->GetUserInformation());
+
+ int pID=aTrack->GetParentID();
+ if (!pID) // this is a primary particle.
+ {
+ const G4PrimaryParticle *pp=0;
+ PrimaryParticleInformation *ppi=0;
+
+ const G4DynamicParticle *dp=aTrack->GetDynamicParticle();
+ if (dp) pp=(const G4PrimaryParticle *)dp->GetPrimaryParticle();
+ if (pp)
+ ppi=dynamic_cast
+ <PrimaryParticleInformation *> (pp->GetUserInformation());
+ if (ppi) {
+ const ISF::ISFParticle* isp = ppi->GetISFParticle();
+ const HepMC::GenParticle *part=ppi->GetHepMCParticle();
+ if (part) {// OK, we got back to HepMC
+ TrackInformation *ti=new TrackInformation(part);
+ ti->SetRegenerationNr(0);
+ ti->SetISFParticle(isp);
+ // regNr=0 and classify=Primary are default values anyway
+ inT->SetUserInformation(ti);
+ ti->SetClassification(Primary);
+
+ }
+ else if (ppi->GetParticleBarcode()>=0) {
+ // PrimaryParticleInformation should at least provide a barcode
+ TrackBarcodeInfo *bi=new TrackBarcodeInfo(ppi->GetParticleBarcode());
+ bi->SetISFParticle(isp);
+ inT->SetUserInformation(bi);
+ }
+ }
+ }
+ else // secondary track: see if it must be saved
+ {
+
+ // Time cut for particles stacking after a certain time...
+ if ( p_stackEnergyCut > 0 &&
+ &particleType == G4Gamma::Gamma() &&
+ aTrack->GetTotalEnergy() < p_stackEnergyCut ){
+ return fKill;
+ }
+
+ }
+ return fUrgent;
+}
+void AthenaStackingAction::NewStage() {;}
+void AthenaStackingAction::PrepareNewEvent() {;}
+
diff --git a/Simulation/G4Atlas/G4AtlasAlg/src/AthenaTrackingAction.cxx b/Simulation/G4Atlas/G4AtlasAlg/src/AthenaTrackingAction.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..27a7a28323261312eb43062dabed0868a673b3b9
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/src/AthenaTrackingAction.cxx
@@ -0,0 +1,52 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "G4AtlasAlg/AthenaTrackingAction.h"
+#include "MCTruth/PrimaryParticleInformation.h"
+#include "MCTruth/TrackInformation.h"
+#include "MCTruth/EventInformation.h"
+#include "MCTruth/TruthStrategyManager.h"
+#include <iostream>
+
+#include "G4DynamicParticle.hh"
+#include "G4PrimaryParticle.hh"
+#include "G4EventManager.hh"
+#include "MCTruth/TrackHelper.h"
+#include "FadsActions/FadsTrackingAction.h"
+#include "MCTruth/AtlasTrajectory.h"
+
+// static AthenaTrackingAction ta;
+
+void AthenaTrackingAction::PreUserTrackingAction(const G4Track* inTrack)
+{
+ static int ilevel=TruthStrategyManager::GetStrategyManager()->
+ GetSecondarySavingLevel();
+// std::cout<<" this is AthenaTrackingAction::PreUserTrackingAction"<<std::endl;
+ G4Track* inT = const_cast<G4Track*> (inTrack);
+ TrackHelper trackHelper(inT);
+ if (trackHelper.IsPrimary() || trackHelper.IsRegisteredSecondary())
+ {
+ HepMC::GenParticle* part=const_cast<HepMC::GenParticle*> (
+ trackHelper.GetTrackInformation()->
+ GetHepMCParticle());
+ EventInformation *eventInfo=TruthStrategyManager::GetStrategyManager()->
+ GetEventInformation();
+ if (trackHelper.IsPrimary()) eventInfo->SetCurrentPrimary(part);
+ eventInfo->SetCurrentlyTraced(part);
+ }
+ if (trackHelper.IsPrimary() ||
+ (trackHelper.IsRegisteredSecondary()&&ilevel>1) ||
+ (trackHelper.IsSecondary()&&ilevel>2))
+ {
+ G4VTrajectory *temp=new AtlasTrajectory(inTrack);
+ FADS::FadsTrackingAction::GetTrackingAction()->SetTraj(temp);
+ }
+}
+
+
+void AthenaTrackingAction::PostUserTrackingAction(const G4Track* /*outTrack*/)
+{
+// std::cout <<"AthenaTrackingAction::PostUserTrackingAction"<<std::endl;
+ FADS::FadsTrackingAction::GetTrackingAction()->ResetTraj();
+}
diff --git a/Simulation/G4Atlas/G4AtlasAlg/src/G4AtlasAlg.cxx b/Simulation/G4Atlas/G4AtlasAlg/src/G4AtlasAlg.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..240b847d41dc6031c03cd76fa7d1621b25ce431d
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/src/G4AtlasAlg.cxx
@@ -0,0 +1,237 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "G4AtlasAlg/G4AtlasAlg.h"
+
+#include "G4AtlasAlg/PreEventActionManager.h"
+#include "G4AtlasAlg/AthenaTrackingAction.h"
+#include "G4AtlasAlg/AthenaStackingAction.h"
+#include "G4AtlasAlg/G4AtlasRunManager.h"
+
+#include "GaudiKernel/AlgFactory.h"
+
+#include "FadsActions/FadsEventAction.h"
+#include "FadsActions/FadsSteppingAction.h"
+#include "FadsActions/FadsRunAction.h"
+#include "FadsActions/FadsTrackingAction.h"
+#include "FadsActions/FadsStackingAction.h"
+#include "FadsKinematics/GeneratorCenter.h"
+#include "FadsPhysics/PhysicsListCatalog.h"
+#include "FadsPhysics/PhysicsListSteering.h"
+#include "FadsGeometry/FadsDetectorConstruction.h"
+
+#include "G4TransportationManager.hh"
+#include "G4RunManagerKernel.hh"
+#include "G4EventManager.hh"
+#include "G4Navigator.hh"
+#include "G4PropagatorInField.hh"
+#include "G4TrackingManager.hh"
+#include "G4StackManager.hh"
+#include "G4UImanager.hh"
+
+#include "EventInfo/EventInfo.h"
+
+/////////////////////////////////////////////////////////////////////////////
+
+
+G4AtlasAlg::G4AtlasAlg(const std::string& name, ISvcLocator* pSvcLocator)
+ : AthAlgorithm(name, pSvcLocator),
+ m_rndmGenSvc("AtDSFMTGenSvc", name)
+{
+ ATH_MSG_DEBUG(std::endl << std::endl << std::endl);
+ ATH_MSG_INFO("++++++++++++ G4AtlasAlg created ++++++++++++" << std::endl << std::endl);
+ declareProperty( "Dll", libList="");
+ declareProperty( "Physics", physList="");
+ declareProperty( "Generator", generator="");
+ declareProperty( "FieldMap", fieldMap="");
+ declareProperty( "RandomGenerator", rndmGen="athena");
+ declareProperty( "KillAllNeutrinos",m_KillAllNeutrinos=true);
+ declareProperty( "KillLowEPhotons",m_KillLowEPhotons=-1.);
+ declareProperty( "ReleaseGeoModel",m_releaseGeoModel=true);
+ declareProperty( "IncludeParentsInG4Event",m_IncludeParentsInG4Event=false);
+ declareProperty( "KillAbortedEvents", m_killAbortedEvents=true);
+ declareProperty( "FlagAbortedEvents", m_flagAbortedEvents=false);
+
+ // Service instantiation
+ declareProperty("AtRndmGenSvc", m_rndmGenSvc);
+
+ p_runMgr = G4AtlasRunManager::GetG4AtlasRunManager();
+ // if (m_noGeomInit) p_runMgr->NoGeomInit();
+
+ G4VUserPhysicsList *thePL=FADS::PhysicsListCatalog::GetInstance()->GetMainPhysicsList();
+ p_runMgr->SetPhysicsList(thePL);
+ p_runMgr->SetUserInitialization(thePL);
+ p_runMgr->SetUserInitialization(new FADS::FadsDetectorConstruction);
+ trackingAction =new AthenaTrackingAction;
+ stackingAction =new AthenaStackingAction;
+ p_runMgr->SetUserAction(FADS::FadsRunAction::GetRunAction());
+ p_runMgr->SetUserAction(FADS::FadsEventAction::GetEventAction());
+ p_runMgr->SetUserAction(FADS::FadsSteppingAction::GetSteppingAction());
+ p_runMgr->SetUserAction(FADS::FadsTrackingAction::GetTrackingAction());
+ p_runMgr->SetUserAction(FADS::FadsStackingAction::GetStackingAction());
+
+ // Verbosities
+ m_verbosities.clear();
+ declareProperty( "Verbosities" , m_verbosities );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+StatusCode G4AtlasAlg::initialize()
+{
+ G4UImanager *ui = G4UImanager::GetUIpointer();
+
+ if (!libList.empty())
+ {
+ ATH_MSG_INFO("G4AtlasAlg specific libraries requested ") ;
+ std::string temp="/load "+libList;
+ ui->ApplyCommand(temp);
+ }
+ if (!physList.empty())
+ {
+ ATH_MSG_INFO("requesting a specific physics list "<< physList) ;
+ std::string temp="/Physics/GetPhysicsList "+physList;
+ ui->ApplyCommand(temp);
+ }
+ FADS::GeneratorCenter * gc = FADS::GeneratorCenter::GetGeneratorCenter();
+ gc->SetIncludeParentsInG4Event( m_IncludeParentsInG4Event );
+ if (!generator.empty())
+ {
+ ATH_MSG_INFO("requesting a specific generator "<< generator) ;
+ gc->SelectGenerator(generator);
+ } else {
+ // make sure that there is a default generator (i.e. HepMC interface)
+ gc->SelectGenerator("AthenaHepMCInterface");
+ }
+ if (!fieldMap.empty())
+ {
+ ATH_MSG_INFO("requesting a specific field map "<< fieldMap) ;
+ ATH_MSG_INFO("the field is initialized straight away") ;
+ std::string temp="/MagneticField/Select "+fieldMap;
+ ui->ApplyCommand(temp);
+ ui->ApplyCommand("/MagneticField/Initialize");
+ }
+
+ stackingAction->KillAllNeutrinos(m_KillAllNeutrinos);
+ stackingAction->KillLowEPhotons(m_KillLowEPhotons);
+ p_runMgr->SetReleaseGeo( m_releaseGeoModel );
+ p_runMgr->SetLogLevel( int(msg().level()) ); // Synch log levels
+
+ ATH_MSG_DEBUG(std::endl << std::endl << std::endl);
+ ATH_MSG_INFO("++++++++++++ G4AtlasAlg initialized ++++++++++++" << std::endl << std::endl);
+
+ ATH_MSG_INFO("Setting checkmode to true");
+ ui->ApplyCommand("/geometry/navigator/check_mode true");
+
+ if (rndmGen=="athena" || rndmGen=="ranecu") {
+ // Set the random number generator to AtRndmGen
+ if (m_rndmGenSvc.retrieve().isFailure()) {
+ ATH_MSG_ERROR("Could not initialize ATLAS Random Generator Service");
+ return StatusCode::FAILURE;
+ }
+ CLHEP::HepRandomEngine* engine = m_rndmGenSvc->GetEngine("AtlasG4");
+ CLHEP::HepRandom::setTheEngine(engine);
+ ATH_MSG_INFO("Random nr. generator is set to Athena");
+ }
+ else if (rndmGen=="geant4" || rndmGen.empty()) {
+ ATH_MSG_INFO("Random nr. generator is set to Geant4");
+ }
+
+ // G4 init moved to PyG4AtlasAlg / G4AtlasEngine
+ /// @todo Reinstate or delete?! This can't actually be called from the Py algs
+ //ATH_MSG_INFO("Firing initialization of G4!!!");
+ //initializeG4();
+
+ return StatusCode::SUCCESS;
+
+}
+
+void G4AtlasAlg::initializeG4()
+{
+
+ if (m_verbosities.size()>0){
+ G4TransportationManager *tm = G4TransportationManager::GetTransportationManager();
+ G4RunManagerKernel *rmk = G4RunManagerKernel::GetRunManagerKernel();
+ G4EventManager *em = G4EventManager::GetEventManager();
+
+ if (m_verbosities.find("Navigator") != m_verbosities.end()){
+ tm->GetNavigatorForTracking()->SetVerboseLevel( atof(m_verbosities.find("Navigator")->second.data()) );
+ }
+ if (m_verbosities.find("Propagator") != m_verbosities.end()){
+ tm->GetPropagatorInField()->SetVerboseLevel( atof(m_verbosities.find("Propagator")->second.data()) );
+ }
+ if (m_verbosities.find("Tracking") != m_verbosities.end()){
+ rmk->GetTrackingManager()->SetVerboseLevel( atof(m_verbosities.find("Tracking")->second.data()) );
+ }
+ if (m_verbosities.find("Stepping") != m_verbosities.end()){
+ rmk->GetTrackingManager()->GetSteppingManager()->SetVerboseLevel( atof(m_verbosities.find("Stepping")->second.data()) );
+ }
+ if (m_verbosities.find("Stacking") != m_verbosities.end()){
+ rmk->GetStackManager()->SetVerboseLevel( atof(m_verbosities.find("Stacking")->second.data()) );
+ }
+ if (m_verbosities.find("Event") != m_verbosities.end()){
+ em->SetVerboseLevel( atof(m_verbosities.find("Event")->second.data()) );
+ }
+ } // End of the setting of verbosities
+
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+StatusCode G4AtlasAlg::finalize() {
+ ATH_MSG_DEBUG(std::endl<<std::endl<<std::endl);
+ ATH_MSG_INFO("++++++++++++ G4AtlasAlg finalized ++++++++++++" <<std::endl<<std::endl);
+
+ ATH_MSG_DEBUG("\t terminating the current G4 run");
+
+ p_runMgr->RunTermination();
+
+ return StatusCode::SUCCESS;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+StatusCode G4AtlasAlg::execute()
+{
+ static int n_Event=0;
+ static PreEventActionManager *preEvent=PreEventActionManager::
+ GetPreEventActionManager();
+ ATH_MSG_DEBUG(std::endl<<std::endl<<std::endl);
+ ATH_MSG_INFO("++++++++++++ G4AtlasAlg execute ++++++++++++" <<std::endl<<std::endl);
+
+ n_Event += 1;
+
+ if (n_Event<=10 || (n_Event%100) == 0) {
+ ATH_MSG_ALWAYS("G4AtlasAlg: Event num. " << n_Event << " start processing");
+ }
+
+ preEvent->Execute();
+
+ ATH_MSG_DEBUG("Calling SimulateG4Event");
+
+ bool abort = p_runMgr->SimulateFADSEvent();
+
+ if (abort) {
+ ATH_MSG_WARNING("Event was aborted !! ");
+ ATH_MSG_WARNING("Simulation will now go on to the next event ");
+ if (m_killAbortedEvents){
+ ATH_MSG_WARNING("setFilterPassed is now False");
+ setFilterPassed(false);
+ }
+ if (m_flagAbortedEvents){
+ const DataHandle<EventInfo> eic = 0;
+ if ( sgSvc()->retrieve( eic ).isFailure() || !eic ){
+ ATH_MSG_WARNING( "Failed to retrieve EventInfo" );
+ } else {
+ // Gotta cast away the const... sadface
+ EventInfo *ei = const_cast< EventInfo * > (&(*eic));
+ ei->setErrorState(EventInfo::Core,EventInfo::Error);
+ ATH_MSG_WARNING( "Set error state in event info!" );
+ }
+ }
+ }
+
+ return StatusCode::SUCCESS;
+}
diff --git a/Simulation/G4Atlas/G4AtlasAlg/src/G4AtlasRunManager.cxx b/Simulation/G4Atlas/G4AtlasAlg/src/G4AtlasRunManager.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..5a4534ad0463a9f5cfa44711a7d21ab7067c21f5
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/src/G4AtlasRunManager.cxx
@@ -0,0 +1,227 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+
+#include "G4AtlasAlg/G4AtlasRunManager.h"
+
+#include "G4StateManager.hh"
+#include "G4UImanager.hh"
+#include "G4Timer.hh"
+#include "G4UserRunAction.hh"
+#include "G4Run.hh"
+
+#include "G4ios.hh"
+
+#include "G4Event.hh"
+
+#include "G4VUserDetectorConstruction.hh"
+#include "G4VUserPhysicsList.hh"
+
+#include "StoreGate/StoreGateSvc.h"
+#include "GaudiKernel/System.h"
+#include "GaudiKernel/IMessageSvc.h"
+
+#include "FadsKinematics/GeneratorCenter.h"
+#include "FadsPhysics/PhysicsListCatalog.h"
+#include "FadsPhysics/DummyPhysicsList.h"
+#include "G4VUserPrimaryGeneratorAction.hh"
+
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include "GeoModelInterfaces/IGeoModelSvc.h"
+#include "GaudiKernel/ISvcLocator.h"
+#include "GaudiKernel/Bootstrap.h"
+#include "GaudiKernel/IMessageSvc.h"
+
+#include "G4LogicalVolumeStore.hh"
+#include "G4RegionStore.hh"
+#include "G4ProductionCutsTable.hh"
+#include "G4GeometryManager.hh"
+
+
+G4AtlasRunManager::G4AtlasRunManager()
+ : G4RunManager(),
+ m_pl(NULL),
+ m_sgSvc(NULL),
+ m_releaseGeo(false),
+ m_log(0)
+{ }
+
+
+G4AtlasRunManager* G4AtlasRunManager::GetG4AtlasRunManager() {
+ static G4AtlasRunManager* thisManager=0;
+ if (!thisManager) thisManager=new G4AtlasRunManager;
+ return thisManager;
+}
+
+
+void G4AtlasRunManager::InitializeGeometry() {
+ G4LogicalVolumeStore *lvs = G4LogicalVolumeStore::GetInstance();
+ for (unsigned int i=0;i<lvs->size();++i){
+ if ( (*lvs)[i]->GetName() == "Muon::MuonSys" ){
+ (*lvs)[i]->SetSmartless( 0.1 );
+ log() << MSG::INFO << "Set smartlessness for Muon::MuonSys to 0.1" << endreq;
+ } else if ( (*lvs)[i]->GetName() == "LArMgr::LAr::EMB::STAC") {
+ (*lvs)[i]->SetSmartless( 0.5 );
+ log() << MSG::INFO << "Set smartlessness for LArMgr::LAr::EMB::STAC to 0.5" << endreq;
+ }
+ }
+
+ if (userDetector) {
+ G4RunManager::InitializeGeometry();
+ } else {
+ log() << MSG::WARNING <<" User Detector not set!!! Geometry NOT initialized!!!"<<endreq;
+ }
+}
+
+
+void G4AtlasRunManager::EndEvent() {
+ log() << MSG::DEBUG << "G4AtlasRunManager::EndEvent" << endreq;
+}
+
+
+void G4AtlasRunManager::InitializePhysics() {
+ kernel->InitializePhysics();
+ physicsInitialized = true;
+}
+
+
+G4Event* G4AtlasRunManager::GenerateEvent(G4int i_event)
+{
+ static FADS::GeneratorCenter* generatorCenter=FADS::GeneratorCenter::GetGeneratorCenter();
+
+ SetCurrentG4Event(i_event);
+
+ generatorCenter->GenerateAnEvent(currentEvent);
+
+ return currentEvent;
+}
+
+
+bool G4AtlasRunManager::SimulateFADSEvent()
+{
+
+// std::cout<<" SimulateFADSEvent : start simulating one event "<<std::endl;
+
+ G4StateManager* stateManager = G4StateManager::GetStateManager();
+ stateManager->SetNewState(G4State_GeomClosed);
+// stateManager->SetNewState(G4State_EventProc);
+
+ // Release GeoModel Geometry if necessary
+ if (m_releaseGeo){
+ ISvcLocator* svcLocator = Gaudi::svcLocator(); // from Bootstrap
+ StoreGateSvc* m_detStore;
+ StatusCode sc=svcLocator->service("DetectorStore",m_detStore);
+ if (sc.isFailure()){
+ log()<<MSG::ERROR<<"G4AtlasRunManager could not access the detector store - PANIC!!!!"<<endreq;
+ abort();
+ }
+
+ IGeoModelSvc* geoModel = 0;
+ StatusCode status=svcLocator->service("GeoModelSvc",geoModel);
+ if(status.isFailure()) {
+ log()<<MSG::WARNING<< " ----> Unable to retrieve GeoModelSvc" << endreq;
+ } else {
+ status = geoModel->clear();
+ if(status.isFailure()) {
+ log()<<MSG::WARNING<< " ----> GeoModelSvc::clear() failed" << endreq;
+ } else {
+ log()<<MSG::INFO<< " ----> GeoModelSvc::clear() succeeded " << endreq;
+ }
+ }
+ m_releaseGeo=false; // Don't do that again...
+ }
+
+
+ currentEvent = GenerateEvent(1);
+ if (currentEvent->IsAborted())
+ {
+ log()<<MSG::WARNING<<"G4AtlasRunManager::SimulateFADSEvent: "<<
+ "Event Aborted at Generator level"<<endreq;
+ currentEvent=0;
+ return true;
+ }
+
+ eventManager->ProcessOneEvent(currentEvent);
+ if (currentEvent->IsAborted())
+ {
+ log()<<MSG::WARNING<<"G4AtlasRunManager::SimulateFADSEvent: "<<
+ "Event Aborted at Detector Simulation level"<<endreq;
+ currentEvent=0;
+ return true;
+ }
+
+ AnalyzeEvent(currentEvent);
+ if (currentEvent->IsAborted())
+ {
+ log()<<MSG::WARNING<<"G4AtlasRunManager::SimulateFADSEvent: "<<
+ "Event Aborted at Analysis level"<<endreq;
+ currentEvent=0;
+ return true;
+ }
+
+// stateManager->SetNewState(G4State_GeomClosed);
+ StackPreviousEvent(currentEvent);
+ bool abort=currentEvent->IsAborted();
+ currentEvent = 0;
+ // std::cout<<" SimulateFADSEvent : done simulating one event "<<std::endl;
+ return abort;
+}
+
+
+void G4AtlasRunManager::RunTermination()
+{
+ // std::cout<<" this is G4AtlasRunManager::RunTermination() "<<std::endl;
+ for(size_t itr=0;itr<previousEvents->size();itr++)
+ { delete (*previousEvents)[itr]; }
+ previousEvents->clear();
+
+ if(userRunAction) userRunAction->EndOfRunAction(currentRun);
+
+ delete currentRun;
+ currentRun = 0;
+ runIDCounter++;
+
+ std::cout << "Changing the state..." << std::endl;
+ G4StateManager* stateManager = G4StateManager::GetStateManager();
+ stateManager->SetNewState(G4State_Idle);
+
+ std::cout << "Opening the geometry back up" << std::endl;
+ G4GeometryManager::GetInstance()->OpenGeometry();
+
+ std::cout << "Terminating the run... State is " << stateManager->GetStateString( stateManager->GetCurrentState() ) << std::endl;
+ kernel->RunTermination();
+ std::cout << "All done..." << std::endl;
+
+ // std::cout<<" setting all pointers in G4AtlasRunManager to 0"<<std::endl;
+ userRunAction=0;
+ userEventAction=0;
+ userSteppingAction=0;
+ userStackingAction=0;
+ userTrackingAction=0;
+ // physicsList=0;
+ userDetector=0;
+ userPrimaryGeneratorAction=0;
+
+ // std::cout<<" this is G4AtlasRunManager::RunTermination(): done "<<std::endl;
+}
+
+void G4AtlasRunManager::SetCurrentG4Event(int iEvent)
+{
+ currentEvent=new G4Event(iEvent);
+}
+
+MsgStream G4AtlasRunManager::log()
+{
+ if (m_log) return *m_log;
+ ISvcLocator* svcLocator = Gaudi::svcLocator();
+ IMessageSvc* p_msgSvc = 0;
+ if (svcLocator->service("MessageSvc", p_msgSvc).isFailure() || !p_msgSvc)
+ std::cout << "G4AtlasRunManager: Trouble getting the message service. Should never happen. Will crash now." << std::endl;
+ m_log = new MsgStream(p_msgSvc,"G4AtlasRunManager");
+ return *m_log;
+}
+
diff --git a/Simulation/G4Atlas/G4AtlasAlg/src/PreEventAction.cxx b/Simulation/G4Atlas/G4AtlasAlg/src/PreEventAction.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..2348c0c77bbde435d12565d8279e9754d7a533ef
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/src/PreEventAction.cxx
@@ -0,0 +1,16 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "G4AtlasAlg/PreEventAction.h"
+#include "G4AtlasAlg/PreEventActionManager.h"
+
+PreEventAction::PreEventAction(std::string n):name(n)
+{
+ PreEventActionManager::GetPreEventActionManager()->
+ RegisterAction(this);
+}
+std::string PreEventAction::GetName()
+{
+ return name;
+}
diff --git a/Simulation/G4Atlas/G4AtlasAlg/src/PreEventActionManager.cxx b/Simulation/G4Atlas/G4AtlasAlg/src/PreEventActionManager.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..589d36f30de9c857fa83ff17a1727498ced700d1
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/src/PreEventActionManager.cxx
@@ -0,0 +1,36 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "G4AtlasAlg/PreEventActionManager.h"
+#include "G4AtlasAlg/PreEventAction.h"
+#include <iostream>
+
+PreEventActionManager* PreEventActionManager::thePointer=0;
+
+PreEventActionManager::PreEventActionManager()
+{
+}
+
+PreEventActionManager* PreEventActionManager::GetPreEventActionManager()
+{
+ if (!thePointer) thePointer=new PreEventActionManager();
+ return thePointer;
+}
+
+void PreEventActionManager::Execute()
+{
+ preEvActionMap::const_iterator it;
+ for (it=theMap.begin();it!=theMap.end();it++)
+ (*it).second->Execute();
+}
+
+void PreEventActionManager::RegisterAction(PreEventAction *act)
+{
+ std::string n=act->GetName();
+ if (theMap.find(n)!=theMap.end())
+ std::cout<<"Attempt of re-defining an existing pre-event action"
+ <<std::endl;
+ else
+ theMap[n]=act;
+}
diff --git a/Simulation/G4Atlas/G4AtlasAlg/src/TruthHepMCEventConverter.cxx b/Simulation/G4Atlas/G4AtlasAlg/src/TruthHepMCEventConverter.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..32f2f5afc050f4fe42f0d28777de727267d52c00
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/src/TruthHepMCEventConverter.cxx
@@ -0,0 +1,468 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "G4AtlasAlg/TruthHepMCEventConverter.h"
+#include "G4Event.hh"
+
+#include "G4PrimaryParticle.hh"
+#include "G4PrimaryVertex.hh"
+#include "G4LorentzVector.hh"
+#include "G4Geantino.hh"
+#include "G4ChargedGeantino.hh"
+#include "G4EventManager.hh"
+#include "G4ParticleTable.hh"
+
+#include "FadsKinematics/VertexManipulator.h"
+#include "FadsKinematics/ParticleManipulator.h"
+#include "FadsKinematics/G4EventAnalyzer.h"
+#include "FadsKinematics/GeneratorCenter.h"
+
+#include "MCTruth/EventInformation.h"
+#include "MCTruth/TruthEvent.h"
+#include "MCTruth/PrimaryParticleInformation.h"
+
+#include "SimHelpers/ServiceAccessor.h"
+
+#include "GeneratorObjects/McEventCollection.h"
+
+#include "GaudiKernel/ISvcLocator.h"
+#include "GaudiKernel/Bootstrap.h"
+#include "StoreGate/StoreGateSvc.h"
+#include "SGTools/BuiltinsClids.h"
+
+#include <map>
+
+typedef HepMC::GenEvent Hevent;
+typedef HepMC::GenVertex Hvertex;
+typedef HepMC::GenParticle Hparticle;
+typedef HepMC::GenEvent::vertex_const_iterator Hv_iterator;
+typedef HepMC::GenVertex::particle_iterator Hp_iterator;
+// A note on unit conventions 26/10/2012
+//
+// HepMC objects store spatial coordinates in millimetres and time
+// coordinates are stored as c_light*t, also recorded in mm.
+//
+// G4PrimaryVertex objects store spatial coordinates in millimetres
+// and time coordinates are stored in ns.
+//
+// Therefore it is necessary to divide the time component by
+// CLHEP::c_light when going from HepMC to Geant4 objects and multiply
+// by CLHEP::c_light when going from Geant4 to HepMC objects.
+//
+// TODO: In the future we should include full conversion of units from
+// HepMC to G4 (CLHEP) i.e. G4_pos = HepMC_pos / HepMC::mm * CLHEP::mm
+// and G4_time = HepMC_time / HepMC::mm * CLHEP::mm / CLHEP::c_light.
+// Need to determine how best to get the HepMC::mm factor. Similary
+// for energies and momenta.
+
+
+class IsFinalState {
+// borrowed from one of the HepMC examples
+// this predicate returns true if the input has no decay vertex
+public:
+ bool operator()( const HepMC::GenParticle* p ) {
+ if ( !p->end_vertex() && p->status()==1 ) return 1;
+ return 0;
+ }
+};
+
+TruthHepMCEventConverter::TruthHepMCEventConverter()
+{
+ printLev=0;
+ m_modded = new std::vector<int>(10);
+}
+
+TruthHepMCEventConverter::~TruthHepMCEventConverter()
+{
+ if (m_modded) delete m_modded;
+}
+
+bool TruthHepMCEventConverter::ValidHepMCVertex(HepMC::GenVertex *testVertex) const
+{
+ bool vertexOK(true);
+ double coordinate(testVertex->position().x());
+ vertexOK&=!(std::isnan(coordinate) || std::isinf(coordinate));
+ coordinate = testVertex->position().y();
+ vertexOK&=!(std::isnan(coordinate) || std::isinf(coordinate));
+ coordinate = testVertex->position().z();
+ vertexOK&=!(std::isnan(coordinate) || std::isinf(coordinate));
+ coordinate = testVertex->position().t();
+ vertexOK&=!(std::isnan(coordinate) || std::isinf(coordinate));
+ return vertexOK;
+}
+
+bool TruthHepMCEventConverter::SeparateVertices(HepMC::GenVertex *v1,HepMC::GenVertex *v2) const
+{
+ const double distMin=0.01*CLHEP::mm;
+ G4LorentzVector vv1(v1->position().x(),v1->position().y(),v1->position().z(),v1->position().t());
+ G4LorentzVector vv2(v2->position().x(),v2->position().y(),v2->position().z(),v2->position().t());
+ return ((sqrt(pow(vv1.x()-vv2.x(),2)+pow(vv1.y()-vv2.y(),2)+
+ pow(vv1.z()-vv2.z(),2))>distMin) || (vv1.t()!=vv2.t()));
+}
+
+G4PrimaryVertex* TruthHepMCEventConverter::ConstructG4PrimaryVertex(HepMC::GenVertex *v) const
+{
+ if(!ValidHepMCVertex(v))
+ {
+ std::cout << "TruthHepMCEventConverter FATAL ConstructG4PrimaryVertex: Received invalid GenVertex object:" << *v
+ << " - Job will die now..." << std::endl;
+ throw;
+ }
+ G4LorentzVector vv(v->position().x(),v->position().y(),v->position().z(),v->position().t());
+ G4PrimaryVertex *vert= new G4PrimaryVertex (vv.x()*CLHEP::mm,
+ vv.y()*CLHEP::mm,
+ vv.z()*CLHEP::mm,
+ vv.t()*CLHEP::mm/CLHEP::c_light);
+ return vert;
+}
+
+void TruthHepMCEventConverter::TransformHepMCParticle(HepMC::GenParticle *hp, G4PrimaryParticle *gp) const
+{
+ hp->set_pdg_id(gp->GetPDGcode()); //this should not change at all, right?
+ CLHEP::Hep3Vector gpv=gp->GetMomentum();
+ const double pmass=gp->GetMass();
+ const double pe = sqrt(gpv.mag2() + pmass*pmass); // this does only change for boosts, etc.
+ hp->set_momentum(CLHEP::HepLorentzVector(gpv.x(),gpv.y(),gpv.z(),pe));
+}
+
+void TruthHepMCEventConverter::TransformHepMCVertex(HepMC::GenVertex *hv,G4PrimaryVertex *gv) const
+{
+ CLHEP::HepLorentzVector hvlv(hv->position().x(),hv->position().y(),hv->position().z(),hv->position().t());
+ if (hvlv.x()!=gv->GetX0()) hvlv.setX(gv->GetX0());
+ if (hvlv.y()!=gv->GetY0()) hvlv.setY(gv->GetY0());
+ if (hvlv.z()!=gv->GetZ0()) hvlv.setZ(gv->GetZ0());
+ const double cT0(gv->GetT0()*CLHEP::c_light);
+ if (hvlv.t()!=cT0) hvlv.setT(cT0);
+ hv->set_position(hvlv);
+}
+
+G4PrimaryParticle* TruthHepMCEventConverter::ConstructG4PrimaryParticle(HepMC::GenParticle *hp) const
+{
+ // check if stable
+ if(hp->status()!=1) return 0;
+
+ int pdgcode=hp->pdg_id();
+ CLHEP::HepLorentzVector lp(hp->momentum().px(),hp->momentum().py(),hp->momentum().pz(),hp->momentum().e());
+ G4ThreeVector& p=lp;
+
+ if (printLev>1)
+ {
+ std::cout << "GenParticle being converted :"<<std::endl;
+ std::cout << "momentum components (in MeV): "<<p<<std::endl;
+ std::cout << "Energy (in MeV): "<<lp.e()<<std::endl;
+ std::cout << "Mass (in MeV): "<<lp.m()<<std::endl;
+ }
+
+ G4PrimaryParticle *part;
+ if (pdgcode==999){
+ G4ParticleDefinition *pdef=
+ G4Geantino::GeantinoDefinition();
+ part=new G4PrimaryParticle(pdef,p.x(),p.y(),p.z());
+
+ } else if (pdgcode==998){
+ G4ParticleDefinition *pdef=
+ G4ChargedGeantino::ChargedGeantinoDefinition();
+ part=new G4PrimaryParticle(pdef,p.x(),p.y(),p.z());
+ } else {
+ part= new G4PrimaryParticle(pdgcode,
+ p.x(),
+ p.y(),
+ p.z());
+ }
+ return part;
+
+}
+
+HepMC::GenVertex* TruthHepMCEventConverter::ConstructHepMCVertex(G4PrimaryVertex *v) const
+{
+ G4ThreeVector pos=v->GetPosition();
+ double t0=v->GetT0();
+ G4LorentzVector lv(pos.x(),pos.y(),pos.z(),t0*CLHEP::c_light) ;
+ HepMC::GenVertex *hv=new HepMC::GenVertex(lv);
+ return hv;
+}
+
+struct CompositeVertex {
+ HepMC::GenVertex* primary;
+ std::vector<HepMC::GenVertex*> secondaries;
+};
+
+typedef std::vector<CompositeVertex> vertexVector;
+
+void TruthHepMCEventConverter::HepMC2G4(const HepMC::GenEvent* evt,G4Event * anEvent, const bool IncludeParentsInG4Event=false) const
+{
+
+ m_modded->clear();
+ int n_vertices=0;
+ int n_particles=0;
+
+ std::vector<G4PrimaryVertex* > vertexVector;
+ std::vector<G4PrimaryParticle* > particleVector;
+
+ if (printLev>3) std::cout<<" This is the Truth HepMC event converter "<<std::endl;
+
+ FADS::GeneratorCenter *gc=FADS::GeneratorCenter::GetGeneratorCenter();
+
+ FADS::vtxModifiers::const_iterator vit;
+ for (vit=gc->BeginVertexManipulator();vit!=gc->EndVertexManipulator();vit++)
+ if ((*vit).second->IsOn()) (*vit).second->EventInitialization();
+
+ FADS::particleModifiers::const_iterator pit;
+ for (pit=gc->BeginParticleManipulator();pit!=gc->EndParticleManipulator();pit++)
+ if ((*pit).second->IsOn()) (*pit).second->EventInitialization();
+
+ FADS::eventAnalyzers::const_iterator eit;
+ for (eit=gc->BeginEventAnalyzer();eit!=gc->EndEventAnalyzer();eit++)
+ if ((*eit).second->IsOn()) (*eit).second->EventInitialization();
+
+ if (printLev>=3) std::cout<<" creating the Truth event "<<std::endl;
+ HepMC::GenEvent* newEvt=new HepMC::GenEvent(*evt);
+
+ StoreGateSvc* storeGate=0;
+ ISvcLocator* svcLocator = Gaudi::svcLocator();
+ StatusCode st=svcLocator->service("StoreGateSvc", storeGate);
+ if (st.isFailure())
+ std::cout<<" storeGateSvc(); could not access StoreGateSvc"<<std::endl;
+ // Grab an instance of the particle table
+ G4ParticleTable *partTable = G4ParticleTable::GetParticleTable();
+ // If we failed, yell
+ if (partTable==0)
+ std::cout << "Could not find a particle table!" << std::endl;
+
+ // store the time in SG
+ double vtxTime = 0;
+ std::vector<HepMC::GenParticle*> incoming; incoming.clear();
+ for (HepMC::GenEvent::particle_iterator p = newEvt->particles_begin(); p != newEvt->particles_end(); ++p){
+ if ( ((*p)->pdg_id() >= 2212) && (!(*p)->production_vertex())) incoming.push_back( *p );
+ if (incoming.size() > 1) break;
+ }
+ if (incoming.size() == 2 &&
+ newEvt->vertices_size() != 0){
+ vtxTime = (*(newEvt->vertices_begin()))->position().z()/CLHEP::c_light;
+
+ // sign issue
+ if (incoming[0]->momentum().z() + incoming[1]->momentum().z() < 0) vtxTime = -vtxTime;
+ }
+ double *p_vtxTime = new double(vtxTime);
+ st=storeGate->record( p_vtxTime, "PrimarySimEventTime" );
+ if (st.isFailure())
+ std::cout << "TruthHepMCEventConverter: could not record primary vertex time to storeGate!"<<std::endl;
+
+// loop on the vertices
+// only those which have stable particles will be translated into G4PrimaryVertex
+// vertices will be displaced if their final position differs from the original
+
+// std::cout<<" primary vertex being processed "<<std::endl<<*(*(newEvt->vertices_begin()))<<std::endl;
+// std::cout<<" primary vertex time set at: "<<vtxInfo->GetPrimaryVertexTime()<<"ns"<<std::endl;
+ for ( Hv_iterator v = newEvt->vertices_begin();v != newEvt->vertices_end(); ++v )
+ {
+ if (printLev>4) std::cout<<" new vertex being processed "<<std::endl<<
+ *(*v)<<std::endl;
+
+ // Check if we've already modified this vertex
+ bool IsUsed=false;
+ for (unsigned int i=0;i<m_modded->size() && IncludeParentsInG4Event;i++){
+ if ( (*m_modded)[i] == (*v)->barcode() ){
+ IsUsed=true;
+ break;
+ }
+ }
+ if (IsUsed){
+ if (printLev>3) std::cout<< "Cutting a vertex because it's been used." << std::endl;
+ continue;
+ }
+
+ int nStableParticles=0;
+ IsFinalState isStable;
+
+ // create a G4PrimaryVertex - destroyed afterwards if not to be kept
+ //
+ G4PrimaryVertex *vert=ConstructG4PrimaryVertex((*v));
+ G4ThreeVector vtxPos1=vert->GetPosition();
+ FADS::vtxModifiers::const_iterator vit;
+ bool vGood=true;
+ for (vit=gc->BeginVertexManipulator();vit!=gc->EndVertexManipulator();++vit)
+ {
+ if ((*vit).second->IsOn())
+ vGood=vGood&&(*vit).second->EditVertex(vert);
+ if (!vGood) break;
+ }
+
+ for (vit=gc->BeginVertexManipulator();vit!=gc->EndVertexManipulator();++vit)
+ {
+ if ((*vit).second->GetName()=="VertexRangeChecker" && (*vit).second->IsOn()){
+ vGood=vGood&&(*vit).second->EditVertex(vert);
+ break;
+ }
+ }
+
+ if (!vGood)
+ {
+ // the vertex must be discarded
+ delete vert;
+ }
+
+
+ if (vGood)
+ {
+ TransformHepMCVertex((*v),vert);
+
+ std::map< int , HepMC::GenParticle*, std::less<int> > pMap;
+ std::map< int , HepMC::GenParticle*, std::less<int> >::const_iterator iit;
+
+ for (Hp_iterator it=(*v)->particles_begin(HepMC::children);
+ it!=(*v)->particles_end(HepMC::children);++it)
+ {
+ int bc=(*it)->barcode();
+ pMap[bc]=(*it);
+ }
+
+ for (iit=pMap.begin();iit!=pMap.end();iit++)
+ {
+ // If it is stable or it is known
+ bool IsKnown = partTable &&
+ partTable->FindParticle( ((*iit).second)->pdg_id() );
+ if (isStable((*iit).second) || (IsKnown && IncludeParentsInG4Event))
+ {
+ G4PrimaryParticle *part =ConstructG4PrimaryParticle((*iit).second);
+ if (0==part) continue;
+
+ bool pGood=true;
+ FADS::particleModifiers::const_iterator jt;
+
+ for (jt=gc->BeginParticleManipulator();jt!=gc->EndParticleManipulator();
+ ++jt)
+ {
+ if ((*jt).second->IsOn())
+ pGood=pGood&&(*jt).second->EditParticle(part);
+ if (!pGood) break;
+ }
+ if (pGood)
+ {
+ if (printLev>3) std::cout<<" particle accepted "<<std::endl;
+ TransformHepMCParticle((*iit).second, part);
+ nStableParticles++;
+ PrimaryParticleInformation* ppi=new PrimaryParticleInformation;
+
+ ppi->SetParticle((*iit).second);
+ ppi->SetRegenerationNr(0);
+ part->SetUserInformation(ppi);
+ vert->SetPrimary(part);
+ particleVector.push_back(part);
+
+ n_particles++;
+ if (IncludeParentsInG4Event &&
+ IsKnown &&
+ ( (*iit).second )->end_vertex() ){
+ ModifyVertex( &(*part), ((*iit).second)->end_vertex() );
+ }
+ if ( !((*iit).second)->end_vertex() && ((*iit).second)->status() != 1){
+ if (printLev>2)
+ std::cout << "Used a known(" << IsKnown << ") part with EV(" << ((*iit).second)->end_vertex()
+ << ") and status code " << ((*iit).second)->status() << std::endl;
+ }
+ }
+ else {
+ delete part;
+ }
+ }
+ }
+ if (nStableParticles){
+ n_vertices++;
+ anEvent->AddPrimaryVertex(vert);
+ vertexVector.push_back(vert);
+ } else {
+ delete vert;
+ }
+ }
+ }
+
+ EventInformation *eventInfo=new EventInformation;
+ eventInfo->SetNrOfPrimaryParticles(n_particles);
+ eventInfo->SetNrOfPrimaryVertices(n_vertices);
+ eventInfo->SetHepMCEvent(newEvt);
+ anEvent->SetUserInformation(eventInfo);
+
+// allow for event-by-event particle-level analysis
+ for (eit=gc->BeginEventAnalyzer();eit!=gc->EndEventAnalyzer();eit++)
+ if ((*eit).second->IsOn()) (*eit).second->EventAnalysis(vertexVector,particleVector);
+
+ // publish the truth to SG
+ if (printLev>=3) std::cout<<" recording the truth event to SG "<<std::endl;
+
+ std::string key = "TruthEvent";
+ McEventCollection* mcEvtColl;
+
+ if (storeGate->contains<McEventCollection>(key))
+ {
+ StatusCode status = storeGate->retrieve(mcEvtColl,key);
+ if (status.isSuccess())
+ std::cout << "found an McEventCollection in store" << std::endl;
+ }
+ else
+ {
+ // McCollection doesn't exist. Create it and fill it with the truth
+ mcEvtColl = new McEventCollection;
+ mcEvtColl->push_back(newEvt);
+ StatusCode status=storeGate->record( mcEvtColl, key );
+ if (status.isFailure())
+ std::cout << "TruthHepMCEventConverter: could not record MC evt collection to storeGate!"<<std::endl;
+ }
+
+}
+
+// Propagating a change using a primary particle and its decay vertex
+void TruthHepMCEventConverter::ModifyVertex(G4PrimaryParticle *pp, HepMC::GenVertex *v) const
+{
+ FADS::GeneratorCenter *gc= FADS::GeneratorCenter::GetGeneratorCenter();
+
+ // One incoming particle
+ G4LorentzVector lv(v->position().x(),
+ v->position().y(),
+ v->position().z(),
+ v->position().t());
+ if ( v->particles_in_size() < 1 ){
+ // Bark
+ std::cout << "ERROR: No particles going into a decay vertex!" << std::endl;
+ return;
+ }
+
+ G4LorentzVector lv0((*(v->particles_begin(HepMC::parents)))->production_vertex()->position().x(),
+ (*(v->particles_begin(HepMC::parents)))->production_vertex()->position().y(),
+ (*(v->particles_begin(HepMC::parents)))->production_vertex()->position().z(),
+ (*(v->particles_begin(HepMC::parents)))->production_vertex()->position().t());
+ pp->SetProperTime( (lv-lv0).mag()/CLHEP::c_light );
+
+ // Loop over daughters
+ for (HepMC::GenVertex::particle_iterator it = v->particles_begin(HepMC::children);
+ it != v->particles_end(HepMC::children);
+ it++){
+
+ // Either way we will add this daughter to the pp's daughters...
+ G4PrimaryParticle *daughter = ConstructG4PrimaryParticle( (*it) );
+ bool pGood=true;
+ FADS::particleModifiers::const_iterator jt;
+
+ for (jt=gc->BeginParticleManipulator();jt!=gc->EndParticleManipulator();
+ ++jt)
+ {
+ if ((*jt).second->IsOn())
+ pGood=pGood&&(*jt).second->EditParticle(daughter);
+ if (!pGood) break;
+ }
+ if (pGood){
+
+ pp->SetDaughter( daughter );
+
+ if ( (*it)->end_vertex() ) {
+ // Send the daughter with its end vertex back to this method
+ ModifyVertex( daughter, (*it)->end_vertex() );
+ }
+ }
+ }
+
+ // add the vertex to our list
+ m_modded->push_back( v->barcode() );
+
+}
diff --git a/Simulation/G4Atlas/G4AtlasAlg/src/components/G4AtlasAlg_entries.cxx b/Simulation/G4Atlas/G4AtlasAlg/src/components/G4AtlasAlg_entries.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..fea85a8f9282d83d1540edea1b4514aa4504fe7d
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/src/components/G4AtlasAlg_entries.cxx
@@ -0,0 +1,10 @@
+#include "GaudiKernel/DeclareFactoryEntries.h"
+
+#include "G4AtlasAlg/G4AtlasAlg.h"
+
+DECLARE_FACTORY_ENTRIES(G4AtlasAlg) {
+ DECLARE_ALGORITHM( G4AtlasAlg )
+}
+
+DECLARE_ALGORITHM_FACTORY( G4AtlasAlg )
+
diff --git a/Simulation/G4Atlas/G4AtlasAlg/src/components/G4AtlasAlg_load.cxx b/Simulation/G4Atlas/G4AtlasAlg/src/components/G4AtlasAlg_load.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..1931cdab48ca67195cafc827321471ca69b4c404
--- /dev/null
+++ b/Simulation/G4Atlas/G4AtlasAlg/src/components/G4AtlasAlg_load.cxx
@@ -0,0 +1,4 @@
+#include "GaudiKernel/LoadFactoryEntries.h"
+
+LOAD_FACTORY_ENTRIES(G4AtlasAlg)
+