Commit c8ac7ad4 authored by Zach Marshall's avatar Zach Marshall Committed by Graeme Stewart
Browse files

Protecting against empty MCEventCollections (G4AtlasAlg-00-03-10)

parent 0decae88
/*
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
/*
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
/*
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
/*
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
/*
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
/*
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
/*
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
/*
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
/*
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
<lcgdict>
<class name="G4AtlasAlg"/>
</lcgdict>
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
/*
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;
}
/*
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() {;}
/*
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);
}
}