Skip to content
Snippets Groups Projects
Commit b701f87b authored by James Beacham's avatar James Beacham Committed by John Chapman
Browse files

Merge branch 'ActuallyAddRadLengthActionClass_21.0' into '21.0'

Add RadLengthAction UserAction

See merge request atlas/athena!14486
parent 3a209a8a
No related branches found
No related tags found
No related merge requests found
/*
Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
*/
#ifndef G4UserActions_RadLengthAction_H
#define G4UserActions_RadLengthAction_H
#include <map>
#include <vector>
#include <string>
#include "G4UserRunAction.hh"
#include "G4UserEventAction.hh"
#include "G4UserSteppingAction.hh"
#include "G4VPhysicalVolume.hh"
#include "G4VSensitiveDetector.hh"
class TTree;
namespace G4UA
{
class RadLengthAction final : public G4UserRunAction,
public G4UserEventAction,
public G4UserSteppingAction
{
public:
struct Config
{
// depth of volume, how many daughters should be taken
// into account, value recieved from jobOptions
int VolumeDepthLevel=1;
};
RadLengthAction(const Config& config);
virtual void BeginOfRunAction(const G4Run*) override;
virtual void EndOfRunAction(const G4Run*) override;
virtual void BeginOfEventAction(const G4Event*) override;
virtual void EndOfEventAction(const G4Event*) override;
virtual void UserSteppingAction(const G4Step*) override;
private:
Config m_config;
// decision if muon chamber or trigger was already hit
// set to false every BeginOfEvent and changed in Stepping
bool MuChamberPassed;
bool MuTriggerPassed;
// primary variables obtained at BeginOfEvent
double etaPrimary, phiPrimary, chargePrimary;
// map of volumes initialized at BeginOfRun and used
// for comparison in Stepping
std::map<std::string,G4VPhysicalVolume*> topvolmap;
// map of trees named like the volumes initalized during
// BeginOfRun and filled at EndOfEvent
std::map<std::string,TTree*> treeMap;
// map of vector which is initialized at BeginOfRun, components
// of the vector are referenced in the branches of the tree with
// specific volume name, in FillVariables the vectors get values
std::map<std::string,std::vector<double> > variables;
// sensitive detectors retrieved at BeginOfRun and need in Stepping
// see ./atlas/MuonSpectrometer/MuonG4/MuonG4SD/share/muonSD.mac
G4VSensitiveDetector* m_SDMDT;
G4VSensitiveDetector* m_SDTGC;
G4VSensitiveDetector* m_SDCSC;
G4VSensitiveDetector* m_SDRPC;
// methode to fill vector stored in variables map (index volume name)
void fillVariables(std::vector<double> varvec, std::string name);
}; //class RadLengthAction
} //namespace G4UA
#endif // G4UserActions_RadLengthAction_H
/*
Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
*/
#ifndef G4USERACTIONS_G4UA__RADLENGTHACTIONTOOL_H
#define G4USERACTIONS_G4UA__RADLENGTHACTIONTOOL_H
// Infrastructure includes
#include "G4AtlasTools/UserActionToolBase.h"
// Local includes
#include "G4UserActions/RadLengthAction.h"
namespace G4UA{
/// @class RadLengthActionTool
/// @brief A tool to manage RadLengthAction actions
///
/// creates one RadLengthAction instance per thread
class RadLengthActionTool: public UserActionToolBase<RadLengthAction>
{
public:
/// standard tool ctor
RadLengthActionTool(const std::string& type, const std::string& name,const IInterface* parent);
protected:
/// creates the action instances
virtual std::unique_ptr<RadLengthAction>
makeAndFillAction(G4AtlasUserActions&) override final;
private:
/// holds the python configuration
RadLengthAction::Config m_config;
}; // class RadLengthActionTool
} // namespace G4UA
#endif
......@@ -91,3 +91,20 @@ def getFluxRecorderTool(name="G4UA::FluxRecorderTool", **kwargs):
def getRadiationMapsMakerTool(name="G4UA::RadiationMapsMakerTool", **kwargs):
return CfgMgr.G4UA__RadiationMapsMakerTool(name, **kwargs)
def getStoppedParticleActionTool(name="G4UA::StoppedParticleActionTool", **kwargs):
# Just have to set the stopping condition
from G4AtlasApps.SimFlags import simFlags
# example custom configuration
if name in simFlags.UserActionConfig.get_Value().keys():
for prop,value in simFlags.UserActionConfig.get_Value()[name].iteritems():
kwargs.setdefault(prop,value)
return CfgMgr.G4UA__StoppedParticleActionTool(name, **kwargs)
def getRadLengthActionTool(name="G4UA::RadLengthActionTool", **kwargs):
from G4AtlasApps.SimFlags import simFlags
# example custom configuration
if name in simFlags.UserActionConfig.get_Value().keys():
for prop,value in simFlags.UserActionConfig.get_Value()[name].iteritems():
kwargs.setdefault(prop,value)
return CfgMgr.G4UA__RadLengthActionTool(name, **kwargs)
......@@ -8,7 +8,6 @@ addTool("G4UserActions.G4UserActionsConf.G4UA__CosmicPerigeeActionTool", "G4UA::
addTool("G4UserActions.G4UserActionsConf.G4UA__G4TrackCounterTool", "G4UA::G4TrackCounterTool")
addTool("G4UserActions.G4UserActionsConf.G4UA__LengthIntegratorTool", "G4UA::LengthIntegratorTool")
addTool("G4UserActions.G4UserActionsConf.G4UA__PhotonKillerTool", "G4UA::PhotonKillerTool")
addTool("G4UserActions.G4UserActionsConf.G4UA__StoppedParticleActionTool", "G4UA::StoppedParticleActionTool")
addTool("G4UserActions.G4UserActionsConfig.getFastIDKillerTool", "G4UA::FastIDKillerTool")
addTool("G4UserActions.G4UserActionsConfig.getFastMBKillerTool", "G4UA::FastMBKillerTool")
......@@ -24,5 +23,6 @@ addTool("G4UserActions.G4UserActionsConfig.getScoringVolumeTrackKillerTool", "G4
addTool("G4UserActions.G4UserActionsConfig.getFluxRecorderTool", "G4UA::FluxRecorderTool")
addTool("G4UserActions.G4UserActionsConfig.getScoringPlaneTool", "G4UA::ScoringPlaneTool")
addTool("G4UserActions.G4UserActionsConfig.getRadiationMapsMakerTool", "G4UA::RadiationMapsMakerTool")
addTool("G4UserActions.G4UserActionsConfig.getStoppedParticleActionTool", "G4UA::StoppedParticleActionTool")
addTool("G4UserActions.G4UserActionsConfig.getRadLengthActionTool", "G4UA::RadLengthActionTool")
/*
Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
*/
#include "G4UserActions/RadLengthAction.h"
#include "SimHelpers/ServiceAccessor.h"
#include "MCTruth/TrackHelper.h"
#include <iostream>
#include <string>
#include <map>
#include <vector>
#include <math.h>
#include <TTree.h>
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/ITHistSvc.h"
#include "GaudiKernel/MsgStream.h"
#include "G4PrimaryParticle.hh"
#include "G4PrimaryVertex.hh"
#include "G4Event.hh"
#include "G4TouchableHistory.hh"
#include "G4Step.hh"
#include "G4StepPoint.hh"
#include "G4Track.hh"
#include "G4SDManager.hh"
#include "G4VSensitiveDetector.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4LogicalVolume.hh"
#include "G4VPhysicalVolume.hh"
#include "G4Material.hh"
namespace G4UA
{
RadLengthAction::RadLengthAction(const Config& config)
: m_config(config)
, m_SDMDT(nullptr)
, m_SDTGC(nullptr)
, m_SDCSC(nullptr)
, m_SDRPC(nullptr)
{
}
// BeginOfRunAction used for initialization of the most crucial maps,
// registering trees and reading depth level from jobOptions
void RadLengthAction::BeginOfRunAction(const G4Run*)
{
// getting level of depth from jobOptions
// depth = atof( theProperties["VolumeDepthLevel"].c_str() );
// creating reference to sensitive detectors of muonsystem
// later on used for comparison in Stepping
G4SDManager *SDM=G4SDManager::GetSDMpointer();
m_SDMDT = (G4VSensitiveDetector*) SDM->FindSensitiveDetector("AMS");
m_SDTGC = (G4VSensitiveDetector*) SDM->FindSensitiveDetector("TMS");
m_SDCSC = (G4VSensitiveDetector*) SDM->FindSensitiveDetector("CMS");
m_SDRPC = (G4VSensitiveDetector*) SDM->FindSensitiveDetector("RMS");
// clearing the maps with comparison volumes, trees and variables
topvolmap.clear();
treeMap.clear();
variables.clear();
// creating volume store
G4LogicalVolumeStore *store = G4LogicalVolumeStore::GetInstance();
// getting most top volume "Atlas::Atlas" by hand
G4LogicalVolume* atlas = store->GetVolume("Atlas::Atlas");
// logical volume vector, initialized with atlas only
// later on filled with all logicals of depth level
std::vector<G4LogicalVolume*> logvec;
logvec.push_back(atlas);
// physical volume vector later on used for topvolmap and names
std::vector<G4VPhysicalVolume*> physvec;
// counter for loop up to depth level
int counter = 0;
// loop to fill physvec with all physical volumes of depth level
while(counter != m_config.VolumeDepthLevel){
// counter increase up to depth
counter++;
// create tmplogvec to be filled with (logical) daughters of counter depth level
std::vector<G4LogicalVolume*> tmplogvec;
// iterator for logical
std::vector<G4LogicalVolume*>::iterator logit;
// looping on logical vector of previous depth level
for(logit=logvec.begin() ; logit < logvec.end(); logit++ ){
for(int k=0; k<(*logit)->GetNoDaughters();k++){
// dumping all logical daughters for next level
tmplogvec.push_back((*logit)->GetDaughter(k)->GetLogicalVolume());
// dumping the corresponding physical if not made of gas or last depth level where all are dumped
if(((*logit)->GetDaughter(k)->GetLogicalVolume()->GetMaterial()->GetRadlen()<100*CLHEP::cm) || (counter == m_config.VolumeDepthLevel)){
physvec.push_back((*logit)->GetDaughter(k));
}
}
}
// logical vector is put to the tmplogvec (with daughters) for next stage
// to physvec all not gas daughters of this level where already added
logvec=tmplogvec;
}
// loop on physvec to put the volumes in the topvolmap using their name
// for comparison in Stepping and initializing all other maps
std::vector<G4VPhysicalVolume*>::iterator physit;
for( physit=physvec.begin() ; physit < physvec.end(); physit++ ){
std::string fulldaughtername = (*physit)->GetName();
std::string daughtername;
std::string::size_type npos;
npos=fulldaughtername.find("::");
// nicer naming for depth level =1
if(npos!=std::string::npos && m_config.VolumeDepthLevel==1) daughtername = fulldaughtername.substr(0,npos);
else daughtername = fulldaughtername;
topvolmap[daughtername]= (*physit);
}
// adding by hand two additional names to topvolmap (without volume)
// due to initialization purpose of topvolmap
topvolmap["ToMuChamber"] = 0;
topvolmap["ToMuTrigger"] = 0;
// using topvolmap to initialize treemap (trees with volumenames)
// and variables map of vectors with FIXED SIZE TO NUMBER OF VARIABLES TO BE DUMPED
std::map<std::string,G4VPhysicalVolume*>::iterator volit;
for(volit=topvolmap.begin(); volit!=topvolmap.end(); volit++){
treeMap[(*volit).first]=new TTree(TString((*volit).first), TString((*volit).first));
variables[(*volit).first].resize(12);
}
// procedure to get the THistoSvc
StatusCode status;
static ITHistSvc* hSvc=0;
if(!hSvc){
ISvcLocator* svcLocator = Gaudi::svcLocator();
status = svcLocator->service("THistSvc", hSvc, true);
if( status.isFailure() ){
return;
}
}
// using already initialized treeMap to register the trees with volname
// and branches which are REFERENCED to the components of the corresponding
// entry in variables map
std::map<std::string,TTree*>::iterator it;
for(it=treeMap.begin(); it!=treeMap.end(); it++){
std::string filename= "/RadLengthAction/";
std::string treepath= filename+(*it).first;
if (hSvc) hSvc->regTree(treepath.c_str(), treeMap[(*it).first]);
//if (!hSvc) log()<< MSG::ERROR << "Cannot register Tree!" << (*it).first << endreq;
treeMap[(*it).first]->Branch("EnergyLoss", &variables[(*it).first].at(0), "EnergyLoss/D");
treeMap[(*it).first]->Branch("RadLength", &variables[(*it).first].at(1), "RadLength/D");
treeMap[(*it).first]->Branch("Intlength", &variables[(*it).first].at(2), "Intlength/D");
treeMap[(*it).first]->Branch("InitialEta", &variables[(*it).first].at(3), "InitialEta/D");
treeMap[(*it).first]->Branch("InitialPhi", &variables[(*it).first].at(4), "InitialPhi/D");
treeMap[(*it).first]->Branch("PointEta", &variables[(*it).first].at(5), "PointEta/D");
treeMap[(*it).first]->Branch("PointPhi", &variables[(*it).first].at(6), "PointPhi/D");
treeMap[(*it).first]->Branch("PointX", &variables[(*it).first].at(7), "PointX/D");
treeMap[(*it).first]->Branch("PointY", &variables[(*it).first].at(8), "PointY/D");
treeMap[(*it).first]->Branch("PointZ", &variables[(*it).first].at(9), "PointZ/D");
treeMap[(*it).first]->Branch("PointR", &variables[(*it).first].at(10), "PointR/D");
treeMap[(*it).first]->Branch("Charge", &variables[(*it).first].at(11), "Charge/D");
}
}
// empty EndOfRunAction
void RadLengthAction::EndOfRunAction(const G4Run*)
{
}
// BeginOfEventAction used to get variables of primary particle and reinitialization
void RadLengthAction::BeginOfEventAction(const G4Event* anEvent)
{
// getting variables of the primary particle
G4PrimaryVertex *vert=anEvent->GetPrimaryVertex(0);
// vert->SetPosition(0.0,0.0,0.0);
G4PrimaryParticle *part=vert->GetPrimary();
G4ThreeVector mom=part->GetMomentum();
etaPrimary=mom.eta();
phiPrimary=mom.phi();
chargePrimary=part->GetCharge();
G4cout<<"Begin Of Event"<<" "<<(double)vert->GetX0()<<" "<<(double)vert->GetY0()<<" "<<(double)vert->GetZ0()<< G4endl;
G4cout<<"Begin Of Event"<<" "<<(double)mom[0]<<" "<<(double)mom[1]<<" "<<(double)mom[2]<< G4endl;
// putting passflags for muon chambers/trigger (back) to false
MuChamberPassed = false;
MuTriggerPassed = false;
// reinitialize the variables vector for this event
std::map<std::string,std::vector<double> >::iterator it;
for(it=variables.begin(); it!=variables.end(); it++){
for(unsigned int i=0; i< variables[(*it).first].size(); i++) variables[(*it).first].at(i)=0;
}
}
// EndOfEvent used for filling all trees (branches are already refering to variables)
void RadLengthAction::EndOfEventAction(const G4Event*)
{
//for the case there was never a chamber passed set eloss to 1 radlength/intlength to -1
if(!MuChamberPassed){
variables["ToMuChamber"].at(0) = 1.;
variables["ToMuChamber"].at(1) = -1.;
variables["ToMuChamber"].at(2) = -1.;
}
//for the case there was never a trigger passed set eloss to 1 radlength/intlength to -1
if(!MuChamberPassed){
variables["ToMuTrigger"].at(0) = 1.;
variables["ToMuTrigger"].at(1) = -1.;
variables["ToMuTrigger"].at(2) = -1.;
}
std::map<std::string,TTree*>::iterator it;
for(it=treeMap.begin(); it!=treeMap.end(); it++){
treeMap[(*it).first]->Fill();
}
}
// Stepping to get the variables and make volume association
void RadLengthAction::UserSteppingAction(const G4Step* aStep)
{
// killing secondaries to save computing time, perhaps
if(TrackHelper(aStep->GetTrack()).IsSecondary()) aStep->GetTrack()->SetTrackStatus(fStopAndKill);
// entering primary particle
if(TrackHelper(aStep->GetTrack()).IsPrimary()) {
// get touchable history since used often
G4TouchableHistory* touchHist =(G4TouchableHistory*)aStep->GetPreStepPoint()->GetTouchable();
// get point before Stepping was started
G4ThreeVector xyz = (G4ThreeVector) aStep->GetPreStepPoint()->GetPosition();
// get radlength and intlength of volumes passed with this Stepping
// double radl=touchHist->GetVolume()->GetLogicalVolume()->GetMaterial()->GetRadlen();
// double intl=touchHist->GetVolume()->GetLogicalVolume()->GetMaterial()->GetNuclearInterLength();
double radl=aStep->GetPreStepPoint()->GetMaterial()->GetRadlen();
double intl=aStep->GetPreStepPoint()->GetMaterial()->GetNuclearInterLength();
// get the steplength for calculation X0 and intlength
double stepl=aStep->GetStepLength();
G4cout<<aStep->GetPreStepPoint()->GetMaterial()->GetName()<<" "<<radl<<" "<<stepl<<" "<<stepl/radl<<" "<<(double) xyz[0]<<" "<<(double) xyz[1]<<" "<< (double) xyz[2]<<G4endl;
// define and fill vector for the FIXED NUMBER of variables which should be dumped later on
const unsigned int nVariablesToSave(12);
std::vector<double> varvec(nVariablesToSave);
varvec.at(0) = (double) aStep->GetDeltaEnergy();
varvec.at(1) = (double) stepl/radl;
varvec.at(2) = (double) stepl/intl;
varvec.at(3) = (double) etaPrimary;
varvec.at(4) = (double) phiPrimary;
varvec.at(5) = (double) xyz.pseudoRapidity();
varvec.at(6) = (double) xyz.phi();
varvec.at(7) = (double) xyz[0];
varvec.at(8) = (double) xyz[1];
varvec.at(9) = (double) xyz[2];
varvec.at(10) = (double) sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
varvec.at(11) = (double) chargePrimary;
// loop on topvolmap to search for mother of current volume on volume depth level
// dump the variables in corresponding component of variables map
std::map<std::string,G4VPhysicalVolume*>::iterator it;
for(it=topvolmap.begin(); it!=topvolmap.end(); it++){
if((*it).second == touchHist->GetVolume(touchHist->GetHistoryDepth()-m_config.VolumeDepthLevel) ){
this->fillVariables(varvec, (*it).first);
}
}
// get sensitive detector of current volume
G4VSensitiveDetector* SD = (G4VSensitiveDetector*) touchHist->GetVolume()->GetLogicalVolume()->GetSensitiveDetector();
// set muon chamber pass flag true up to end of event and dump
// variables collected up to this point in variables map
if(SD==m_SDCSC || SD==m_SDMDT) MuChamberPassed = true;
if(!MuChamberPassed) this->fillVariables(varvec, "ToMuChamber");
// set muon trigger pass flag true up to end of event and dump
// variables collected up to this point in variables map
if(SD==m_SDRPC || SD==m_SDTGC) MuTriggerPassed = true;
if(!MuTriggerPassed) this->fillVariables(varvec, "ToMuTrigger");
}//end primary particle
}
// how variables should be dumped in variables map
void RadLengthAction::fillVariables(std::vector<double> varvec, std::string name){
// first three components should be added (deltaenergy, radlength, intlength)
for(unsigned int i = 0; i<3; i++) variables[name].at(i) += varvec.at(i);
// other components should be overwritten
// (last value inside the volume will be kept by change and dumped in EndOfEvent)
for(unsigned int j = 3; j<12; j++) variables[name].at(j) = varvec.at(j);
}
} //namespace G4UA
/*
Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
*/
#include "G4UserActions/RadLengthActionTool.h"
namespace G4UA{
RadLengthActionTool::RadLengthActionTool(const std::string& type,
const std::string& name,
const IInterface* parent)
: UserActionToolBase<RadLengthAction>(type, name, parent)
, m_config()
{
declareProperty("VolumeDepthLevel",m_config.VolumeDepthLevel);
}
std::unique_ptr<RadLengthAction>
RadLengthActionTool::makeAndFillAction(G4AtlasUserActions& actionList)
{
ATH_MSG_DEBUG("Constructing a RadLengthAction action");
auto action = std::make_unique<RadLengthAction>(m_config);
actionList.runActions.push_back( action.get() );
actionList.eventActions.push_back( action.get() );
actionList.steppingActions.push_back( action.get() );
return action;
}
} // namespace G4UA
......@@ -14,6 +14,7 @@
#include "G4UserActions/FluxRecorderTool.h"
#include "G4UserActions/ScoringPlaneTool.h"
#include "G4UserActions/RadiationMapsMakerTool.h"
#include "G4UserActions/RadLengthActionTool.h"
#include "../TestActionTool.h"
DECLARE_COMPONENT( G4UA::G4SimTimerTool )
......@@ -32,4 +33,5 @@ DECLARE_COMPONENT( G4UA::StoppedParticleActionTool )
DECLARE_COMPONENT( G4UA::FluxRecorderTool )
DECLARE_COMPONENT( G4UA::ScoringPlaneTool )
DECLARE_COMPONENT( G4UA::RadiationMapsMakerTool )
DECLARE_COMPONENT( G4UA::RadLengthActionTool )
DECLARE_COMPONENT( G4UA::TestActionTool )
# Include the RadLengthAction UserAction
# Used for production of X0 maps in the muon system
from G4AtlasApps.SimFlags import simFlags
simFlags.OptionalUserActionList.addAction('G4UA::RadLengthActionTool',['Event','Step','Run'])
#simFlags.UserActionConfig.addConfig('G4UA::RadLengthActionTool','VolumeDepthLevel',1)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment