Skip to content
Snippets Groups Projects
Commit a049c98e authored by Graeme Stewart's avatar Graeme Stewart
Browse files

Simulation/G4Extensions/PionDecayer deleted from 21.0

Former-commit-id: 8b03ff63
parent 26283f8c
No related branches found
No related tags found
No related merge requests found
################################################################################
# Package: PionDecayer
################################################################################
# Declare the package name:
atlas_subdir( PionDecayer )
# Declare the package's dependencies:
atlas_depends_on_subdirs( PUBLIC
Simulation/G4Sim/FADS/FadsKinematics
PRIVATE
Simulation/G4Sim/MCTruth )
# External dependencies:
find_package( CLHEP )
find_package( Geant4 )
find_package( ROOT COMPONENTS Core Tree MathCore Hist RIO pthread )
find_package( XercesC )
# this line failed automatic conversion in cmt2cmake :
# macro_prepend PionDecayerDict_shlibflags " -lPionDecayer "
# Component(s) in the package:
atlas_add_component( PionDecayer
src/*.cxx
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${XERCESC_INCLUDE_DIRS} ${GEANT4_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS}
LINK_LIBRARIES ${ROOT_LIBRARIES} ${XERCESC_LIBRARIES} ${GEANT4_LIBRARIES} ${CLHEP_LIBRARIES} FadsKinematics MCTruth )
atlas_add_dictionary( PionDecayerDict
PionDecayer/PionDecayerDict.h
PionDecayer/selection.xml
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${XERCESC_INCLUDE_DIRS} ${GEANT4_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS}
LINK_LIBRARIES ${ROOT_LIBRARIES} ${XERCESC_LIBRARIES} ${GEANT4_LIBRARIES} ${CLHEP_LIBRARIES} FadsKinematics MCTruth )
# Install files from the package:
atlas_install_headers( PionDecayer )
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#ifndef PionDecayer_H
#define PionDecayer_H
#include "FadsKinematics/G4EventAnalyzer.h"
#include <utility>
#include <vector>
class PionDecayer: public FADS::G4EventAnalyzer {
private:
G4int verboseLevel;
static double alenz;
static double cvrad;
static double ptmuon;
static bool m_useCustom;
int m_motherId;
std::vector<int> *m_childrenId;
/* select pions and kaons coming from a decay mode specified through jobOptions */
std::vector<G4PrimaryParticle*> customDecayMode( std::vector<G4PrimaryParticle*>& part );
public:
PionDecayer(std::string s);
~PionDecayer();
inline void switchOn(){TurnOn();};
inline void switchOff(){TurnOff();};
void SetVerboseLevel(G4int i) { verboseLevel=i; }
G4int GetVerboseLevel() const { return verboseLevel; }
bool GetStatus() {return IsOn();}
// for the C++ layer only
void EventInitialization();
void EventAnalysis(std::vector<G4PrimaryVertex*>&, std::vector<G4PrimaryParticle*>&);
static void SetIDETlength(double);
static void SetIDETradius(double);
static void SetPtmuon(double);
static void UseCustomDecayMode();
void SetMother(int);
void AddDecayProduct(int);
private:
PionDecayer(const PionDecayer&); //declaring, but not defining as an object of this class should never be copied.
PionDecayer& operator= (const PionDecayer&); //declaring, but not defining as an object of this class should never be copied.
};
#endif
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#include "PionDecayer/PionDecayer.h"
<lcgdict>
<class name="PionDecayer"/>
<exclusion>
<class name="PionDecayer">
<method name="EventInitialization"/>
<method name="EventAnalysis"/>
</class>
</exclusion>
</lcgdict>
package PionDecayer
version PionDecayer-00-00-00
branches src cmt
use AtlasPolicy AtlasPolicy-*
use FadsKinematics FadsKinematics-* Simulation/G4Sim/FADS
include_dirs "$(PionDecayer_root)" "$(PionDecayer_root)/dict"
library PionDecayer *.cxx
apply_pattern component_library
private
use AtlasReflex AtlasReflex-* External
apply_pattern lcgdict dict=PionDecayer \
headerfiles="../PionDecayer/PionDecayerDict.h" \
selectionfile=selection.xml
macro_prepend PionDecayerDict_shlibflags " -lPionDecayer "
use AtlasCLHEP AtlasCLHEP-* External
use Geant4 Geant4-* External
use MCTruth MCTruth-* Simulation/G4Sim
end_private
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
/**
@mainpage PionDecayer
@author Andrea Dell'Acqua (dellacqu@mail.cern.ch)
@section PionDecayer Introduction
This package was introduced to ensure a pion decays inside the inner detector.
@section PionDecayer Class Overview
There is only one class in the package:
- PionDecayer : This controls the decays based on several input parameters.
@ref used_PionDecayer
@ref requirements_PionDecayer
*/
/**
@page used_PionDecayer Used Packages
@htmlinclude used_packages.html
*/
/**
@page requirements_PionDecayer Requirements
@include requirements
*/
###############################################################
#
# jobOptions example for using PionDecayer in decay mode with
# decay specified to Bs0 -> K- mu+ nu
#
###############################################################
#--- Detector flags -------------------------------------------
from AthenaCommon.DetFlags import DetFlags
# - Select detectors
DetFlags.ID_setOn()
DetFlags.Calo_setOn()
DetFlags.Muon_setOn()
DetFlags.Truth_setOn()
#--- Simulation flags -----------------------------------------
from G4AtlasApps.SimFlags import SimFlags
SimFlags.load_atlas_flags()
SimFlags.SimLayout.set_On()
## Global conditions tag
from AthenaCommon.GlobalFlags import globalflags
globalflags.ConditionsTag = "OFLCOND-SDR-BS7T-04-03"
from AthenaCommon.AthenaCommonFlags import athenaCommonFlags
athenaCommonFlags.PoolHitsOutput = 'G4Hits.pool.root'
athenaCommonFlags.PoolEvgenInput = ['/afs/cern.ch/atlas/offline/ProdData/15.6.11.3/mu_E200_eta0-60-10000.evgen.pool.root']
athenaCommonFlags.EvtMax = 10
def load_pionDecayer():
from G4AtlasApps import AtlasG4Eng
AtlasG4Eng.G4Eng.load_Dict("PionDecayerDict")
pionDecayer = AtlasG4Eng.G4Eng.gbl.PionDecayer("PionDecayer")
pionDecayer.switchOn()
pionDecayer.UseCustomDecayMode()
pionDecayer.SetMother(531)
pionDecayer.AddDecayProduct(-321)
pionDecayer.AddDecayProduct(-13)
pionDecayer.AddDecayProduct(14)
SimFlags.initFunctions.add_function("preInitPhysics", load_pionDecayer)
## Populate alg sequence
from AthenaCommon.AlgSequence import AlgSequence
topSeq = AlgSequence()
from G4AtlasApps.PyG4Atlas import PyG4AtlasAlg
topSeq += PyG4AtlasAlg()
###############################################################
#
# jobOptions example for using PionDecayer in threshold mode
# with Pt threshold at 2GeV
#
###############################################################
#--- Detector flags -------------------------------------------
from AthenaCommon.DetFlags import DetFlags
# - Select detectors
DetFlags.ID_setOn()
DetFlags.Calo_setOn()
DetFlags.Muon_setOn()
DetFlags.Truth_setOn()
#--- Simulation flags -----------------------------------------
from G4AtlasApps.SimFlags import SimFlags
SimFlags.load_atlas_flags()
SimFlags.SimLayout.set_On()
## Global conditions tag
from AthenaCommon.GlobalFlags import globalflags
globalflags.ConditionsTag = "OFLCOND-SDR-BS7T-04-03"
from AthenaCommon.AthenaCommonFlags import athenaCommonFlags
athenaCommonFlags.PoolHitsOutput = 'G4Hits.pool.root'
athenaCommonFlags.PoolEvgenInput = ['/afs/cern.ch/atlas/offline/ProdData/15.6.11.3/mu_E200_eta0-60-10000.evgen.pool.root']
athenaCommonFlags.EvtMax = 10
def load_pionDecayer():
from G4AtlasApps import AtlasG4Eng
AtlasG4Eng.G4Eng.load_Dict("PionDecayerDict")
pionDecayer = AtlasG4Eng.G4Eng.gbl.PionDecayer("PionDecayer")
pionDecayer.switchOn()
pionDecayer.SetPtmuon(2*GeV)
# Uncomment to change decay volume
#pionDecayer.SetIDETlength(340*cm)
#pionDecayer.SetIDETradius(110*cm)
SimFlags.initFunctions.add_function("preInitPhysics", load_pionDecayer)
## Populate alg sequence
from AthenaCommon.AlgSequence import AlgSequence
topSeq = AlgSequence()
from G4AtlasApps.PyG4Atlas import PyG4AtlasAlg
topSeq += PyG4AtlasAlg()
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#include "PionDecayer/PionDecayer.h"
#include "G4LorentzVector.hh"
#include "G4ParticleDefinition.hh"
#include "G4RunManager.hh"
#include "G4EventManager.hh"
#include "G4Event.hh"
#include <CLHEP/Random/RandFlat.h>
#include "MCTruth/PrimaryParticleInformation.h"
static PionDecayer filter("PionDecayer");
double PionDecayer::alenz=350.*CLHEP::cm;
double PionDecayer::cvrad=115.*CLHEP::cm;
double PionDecayer::ptmuon=1.*CLHEP::GeV;
bool PionDecayer::m_useCustom=false;
PionDecayer::PionDecayer(std::string s): FADS::G4EventAnalyzer(s),
verboseLevel(1),
m_motherId(0)
{
// TurnOn();
m_childrenId = new std::vector<int>;
}
PionDecayer::~PionDecayer()
{
if (m_childrenId) delete m_childrenId;
}
void PionDecayer::SetIDETlength(double l)
{
alenz=l;
}
void PionDecayer::SetIDETradius(double r)
{
cvrad=r;
}
void PionDecayer::SetPtmuon(double s)
{
ptmuon=s;
}
void PionDecayer::UseCustomDecayMode()
{
m_useCustom=true;
}
void PionDecayer::SetMother(int id)
{
m_motherId = id;
}
void PionDecayer::AddDecayProduct(int id)
{
m_childrenId->push_back(id);
}
void PionDecayer::EventInitialization()
{
std::cout << " Pion Decayer :EventInitialization() called here "<<std::endl;
}
void PionDecayer::EventAnalysis(std::vector<G4PrimaryVertex*>& /*vert*/, std::vector<G4PrimaryParticle*>& part) {
std::cout <<" Pion Decayer :EventAnalysis called here "<<std::endl;
std::cout <<" ptmuon cut set to "<<ptmuon<< " MeV "<<std::endl;
std::cout <<" fiducial volume dimensions "<<alenz<<" "<<cvrad<<" mm "<<std::endl;
std::vector<G4PrimaryParticle*> particles;
std::vector<G4PrimaryParticle*> particles2;
bool caseSelect=false;
if (part.size()) {
std::cout << "++++ number of particles in the event ++++ "<<part.size()<<std::endl;
/// first loop to take care of those cases where the particle to decay was already
/// chosen at the generator level
if ( m_useCustom ) {
std::cout << "++++ Selected decay channel ++++" << std::endl;
std::cout <<" Mother id: " << m_motherId << std::endl;
for (unsigned int i=0 ; i<m_childrenId->size() ; i++ )
std::cout <<" Child id: " << (*m_childrenId)[i] << std::endl;
particles2 = customDecayMode( part );
if (particles2.size()) caseSelect=true;
}
/// second loop
/// if no particles were selected on the generator level, pick a particle to decay
if (!caseSelect) {
for (unsigned int i=0;i<part.size();i++) {
G4PrimaryParticle *p=part[i];
int pdgcode=p->GetPDGcode();
double p_X=p->GetPx();
double p_Y=p->GetPy();
double ptmod=sqrt(p_X*p_X+p_Y*p_Y);
if((ptmod>ptmuon) && ((abs(pdgcode) ==211) || (abs(pdgcode) ==321)))
particles.push_back(part[i]);
}
if (!particles.size()) {
std::cout << " no particles found that pass the threshold: aborting the event!"<<std::endl;
G4Event *ev=const_cast<G4Event *>(G4RunManager::GetRunManager()->GetCurrentEvent());
ev->SetEventAborted();
return;
}
std::cout << "++++ number of muons and kaons in the event ++++ "<<particles.size()<<std::endl;
double rn=CLHEP::RandFlat::shoot();
unsigned int i=int(rn*particles.size());
particles2.push_back(particles[i]);
}
/// Do the decaying
std::cout << "Number of particles to decay: " << particles2.size() << std::endl;
for (unsigned int i=0;i<particles2.size();i++) {
G4PrimaryParticle *pa=particles2[i];
/// flag particle by setting status=9999
HepMC::GenParticle *hPart=const_cast<HepMC::GenParticle *>(((PrimaryParticleInformation*)(pa->GetUserInformation()))->GetHepMCParticle());
hPart->set_status(9999);
double p_X=pa->GetPx();
double p_Y=pa->GetPy();
double p_Z=pa->GetPz();
double mass=pa->GetG4code()->GetPDGMass();
double pmod=sqrt(p_X*p_X+p_Y*p_Y+p_Z*p_Z);
double ptmod=sqrt(p_X*p_X+p_Y*p_Y);
double E=sqrt(p_X*p_X+p_Y*p_Y+p_Z*p_Z + mass*mass);
G4LorentzVector vL(G4ThreeVector(p_X,p_Y,p_Z),E);
double beta = vL.beta();
double gamma = vL.gamma();
double betaGamma=beta*gamma;
double thetac=atan(cvrad/alenz);
double thetap=atan(ptmod/fabs(p_Z));
double decayL=std::min(cvrad/sin(thetap),alenz/cos(thetap));
double randNum=CLHEP::RandFlat::shoot();
decayL = decayL*randNum;
std::cout <<" mass "<<mass<<std::endl;
std::cout <<" pmod "<<pmod<<std::endl;
std::cout <<" decayL "<<decayL<<std::endl;
std::cout <<" thetac "<<thetac<<" decayL "<<decayL<<" ptmod "<<ptmod<<std::endl;
std::cout <<" ratio "<<mass/pmod<<std::endl;
std::cout <<" beta,gamma,betaGamma "<<beta<<" "<<gamma<<" "<<betaGamma<<std::endl;
std::cout <<" set c*tau to (decayL/betaGamma) "<<decayL/betaGamma<<std::endl;
pa->SetProperTime(decayL/(betaGamma*CLHEP::c_light));
}
}
}
std::vector<G4PrimaryParticle*> PionDecayer::customDecayMode( std::vector<G4PrimaryParticle*>& part) {
std::vector<G4PrimaryParticle*> selectedParticles;
bool particleInList=false;
for (unsigned int i=0;i<part.size();i++) {
G4PrimaryParticle *p=part[i];
/// Check if the current particle is in the list of children pdgIds
int pdgId=((PrimaryParticleInformation*)(p->GetUserInformation()))->GetHepMCParticle()->pdg_id();
for ( unsigned int j=0 ; j<m_childrenId->size() ; j++ )
if ( pdgId==(*m_childrenId)[j] ) particleInList=true;
if ( particleInList ) {
const HepMC::GenVertex* vertex = ((PrimaryParticleInformation*)(p->GetUserInformation()))->GetHepMCParticle()->production_vertex();
HepMC::GenVertex::particles_in_const_iterator mother = vertex->particles_in_const_begin();
if ( *mother!=NULL ) {
/// Check if the particle come from decay of the correct particle.
/// If it does, check that the rest of the decay products are the correct ones.
int id = (*mother)->pdg_id();
if (id==m_motherId) {
unsigned int nbrMatch = 0;
HepMC::GenVertex::particles_out_const_iterator partItr = vertex->particles_out_const_begin();
HepMC::GenVertex::particles_out_const_iterator partItrE = vertex->particles_out_const_end();
for ( ; partItr!=partItrE ; partItr++ ) {
for ( unsigned int j=0 ; j<m_childrenId->size() ; j++ ) {
if ( (*partItr)->pdg_id()==(*m_childrenId)[j] ) {
nbrMatch++;
}
}
}
/// If all decay products are matched to specified children, this is the correct decay.
/// If this was a pion or a kaon, put the particle in the list of particles which should
/// be decayed by PionDecayer
if ( nbrMatch==m_childrenId->size() && ( std::abs(pdgId)==211 || std::abs(pdgId)==321 ) ) {
std::cout << "adding " << pdgId << " to decay list" << std::endl;
selectedParticles.push_back(part[i]);
}
}
}
}
}
return selectedParticles;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment