Commit 647264f1 authored by James Catmore's avatar James Catmore Committed by Graeme Stewart
Browse files

Adding TruthEvent to HepMC SG keys for which xAOD truth building runs...

Adding TruthEvent to HepMC SG keys for which xAOD truth building runs (DerivationFrameworkMCTruth-00-01-28)
parent 1608dfc6
///////////////////////// -*- C++ -*- /////////////////////////////
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
// CompactHardTruth.h
// Header file for class CompactHardTruth
// Author: Frank Paige <paige@bnl.gov>
//
// Algorithm to construct small, consistent hard-scattering HepMC truth
// record from GEN_EVENT. Only final partons have physical momenta, so
// must recluster them. Main steps:
// (1) Identify hadronization vertices with partons in, hadrons out.
// (2) Detach their incoming partons and delete their descendants.
// (3) Iteratively recombine 1->2 branchings of final partons with
// mass M < Mcut. Split 2->1 nonperturbative color recombinations.
// (4) Remove partons not connected to hard process. May make apparent
// anomalies, e.g., dropping soft q from incoming g -> g q.
//
// New final partons have correct momenta and ancestors/descendants.
// CompactHardTruth changes particles/vertices, so it cannot be used directly
// with ThinningSvc, but it can be used in parallel.
//
///////////////////////////////////////////////////////////////////
#ifndef DERIVATIONFRAMEWORK_COMPACTHARDTRUTH_H
#define DERIVATIONFRAMEWORK_COMPACTHARDTRUTH_H 1
// STL includes
#include <string>
// FrameWork includes
#include "AthenaBaseComps/AthAlgorithm.h"
#include "GaudiKernel/ToolHandle.h"
#include "HepMC/GenEvent.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
// fwd declares
class IMcVtxFilterTool;
class TH1F;
namespace DerivationFramework {
class CompactHardTruth
: public ::AthAlgorithm
{
///////////////////////////////////////////////////////////////////
// Public methods:
///////////////////////////////////////////////////////////////////
public:
// Copy constructor:
/// Constructor with parameters:
CompactHardTruth( const std::string& name, ISvcLocator* pSvcLocator );
/// Destructor:
virtual ~CompactHardTruth();
// Assignment operator:
//CompactHardTruth &operator=(const CompactHardTruth &alg);
// Athena algorithm Hooks
virtual StatusCode initialize();
virtual StatusCode execute();
virtual StatusCode finalize();
///////////////////////////////////////////////////////////////////
// Const methods:
///////////////////////////////////////////////////////////////////
// Parton is quark or gluon
virtual bool isParton( const HepMC::GenParticle* );
// Final parton is quark or gluon ending in vertex giving !isParton
virtual bool isFinalParton( const HepMC::GenParticle* );
// Hadron excludes leptons and BSM particles
virtual bool isHadron( const HepMC::GenParticle* p );
// Total in/out FourVector for vertex
virtual HepMC::FourVector vtxInMom(HepMC::GenVertex*);
virtual HepMC::FourVector vtxOutMom(HepMC::GenVertex*);
///////////////////////////////////////////////////////////////////
// Non-const methods:
///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
// Private data:
///////////////////////////////////////////////////////////////////
private:
/// Default constructor:
CompactHardTruth();
/// Containers
std::string m_mcEventsName;
std::string m_thinnedMcEventsName;
// Variables
int m_evtCount;
int m_missCount;
float m_partonCut;
float m_hardCut;
int m_dangleFound;
int m_dangleRemoved;
float m_danglePtMax;
float m_danglePtCut;
int m_maxCount;
int m_thinParticles;
int m_thinVertices;
};
// I/O operators
//////////////////////
///////////////////////////////////////////////////////////////////
// Inline methods:
///////////////////////////////////////////////////////////////////
} //> end namespace DerivationFramework
#endif //> !DERIVATIONFRAMEWORK_COMPACTHARDTRUTH_H
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
// ======================================================================
// DecayGraphHelper.h
// James.Catmore@cern.ch
// Helper struct which implements two simple recursive functions;
// descends or ascends through a decay from a head particle provided by
// the user, and records all particles and vertices it encounters in the
// masks
//
// Also contains functions shared by TruthDressing tool and
// TruthIsolationTool.
// ======================================================================
#pragma once
#include "xAODTruth/TruthParticleContainer.h"
#include "xAODTruth/TruthVertexContainer.h"
#include "HepPID/ParticleIDMethods.hh"
#include<unordered_set>
namespace DerivationFramework {
struct DecayGraphHelper {
// descendants: starting from particle (simple listing version)
void descendants(const xAOD::TruthParticle* pHead,
std::vector<int> &particleList,
std::unordered_set<int> &encounteredBarcodes) {
// Check that this barcode hasn't been seen before (e.g. we are in a loop)
//if ( find(encounteredBarcodes.begin(),encounteredBarcodes.end(),pHead->barcode()) != encounteredBarcodes.end()) return;
std::unordered_set<int>::const_iterator found = encounteredBarcodes.find(pHead->barcode());
if (found!=encounteredBarcodes.end()) return;
encounteredBarcodes.insert(pHead->barcode());
// Get the decay vertex
const xAOD::TruthVertex* decayVtx(0);
if (pHead->hasDecayVtx()) {decayVtx = pHead->decayVtx();}
else {return;}
// Save the PDG ID number of the particle
particleList.push_back(pHead->pdgId());
// Get children particles and self-call
int nChildren = decayVtx->nOutgoingParticles();
for (int i=0; i<nChildren; ++i) descendants(decayVtx->outgoingParticle(i),particleList,encounteredBarcodes);
return;
}
// descendants: starting from particle
void descendants(const xAOD::TruthParticle* pHead,
std::vector<bool> &particleMask,
std::vector<bool> &vertexMask,
std::unordered_set<int> &encounteredBarcodes,
bool includeGeant) {
// Check that this barcode hasn't been seen before (e.g. we are in a loop)
//if ( find(encounteredBarcodes.begin(),encounteredBarcodes.end(),pHead->barcode()) != encounteredBarcodes.end()) return;
std::unordered_set<int>::const_iterator found = encounteredBarcodes.find(pHead->barcode());
if (found!=encounteredBarcodes.end()) return;
encounteredBarcodes.insert(pHead->barcode());
// Save the particle position in the mask
// Note: is this safe???
int headIndex = pHead->index();
particleMask[headIndex] = true;
// Get the decay vertex
const xAOD::TruthVertex* decayVtx(0);
if (pHead->hasDecayVtx()) {decayVtx = pHead->decayVtx();}
else {return;}
// If user doesn't want Geant, check and reject
if (!includeGeant && pHead->barcode()>200000) return;
// Save the decay vertex
int vtxIndex = decayVtx->index();
vertexMask[vtxIndex] = true;
// Get children particles and self-call
int nChildren = decayVtx->nOutgoingParticles();
for (int i=0; i<nChildren; ++i) descendants(decayVtx->outgoingParticle(i),particleMask,vertexMask,encounteredBarcodes,includeGeant);
return;
}
// ancestors: starting from particle
void ancestors(const xAOD::TruthParticle* pHead,
std::vector<bool> &particleMask,
std::vector<bool> &vertexMask,
std::unordered_set<int> &encounteredBarcodes ) {
// Check that this barcode hasn't been seen before (e.g. we are in a loop)
//if ( find(encounteredBarcodes.begin(),encounteredBarcodes.end(),pHead->barcode()) != encounteredBarcodes.end()) return;
std::unordered_set<int>::const_iterator found = encounteredBarcodes.find(pHead->barcode());
if (found!=encounteredBarcodes.end()) return;
encounteredBarcodes.insert(pHead->barcode());
// Save particle position in the mask
int headIndex = pHead->index();
particleMask[headIndex] = true;
// Get the production vertex
const xAOD::TruthVertex* prodVtx(0);
if (pHead->hasProdVtx()) {prodVtx = pHead->prodVtx();}
else {return;}
// Save the production vertex
int vtxIndex = prodVtx->index();
vertexMask[vtxIndex] = true;
// Get children particles and self-call
int nParents = prodVtx->nIncomingParticles();
for (int i=0; i<nParents; ++i) ancestors(prodVtx->incomingParticle(i),particleMask,vertexMask,encounteredBarcodes);
return;
}
void constructListOfFinalParticles(const xAOD::TruthParticleContainer* allParticles,
std::vector<const xAOD::TruthParticle*> &selectedlist,
const std::vector<int> &pdgId, bool usePhotonsFromHadrons = false,
bool chargedOnly = false) const
{
//fill the vector selectedlist with only particles with abs(ID) in the list
//pdgID that are 'good' for dressing: status==1, barcode < 2e5 (for the
//photons), not from hadron decay,
//skip pdgId check if pdgId is an empty vector,
//ignore photons coming from hadrons if usePhotonsFromHadrons=false,
//only use charged particles if chargedOnly=true
bool skipPdgCheck = (pdgId.size()==0);
//bypass the hadron veto?
for (xAOD::TruthParticleContainer::const_iterator pItr=allParticles->begin(); pItr!=allParticles->end(); ++pItr) {
const xAOD::TruthParticle *particle = *pItr;
if (particle->status() != 1) continue;
if (!skipPdgCheck && find(pdgId.begin(), pdgId.end(), abs(particle->pdgId())) == pdgId.end()) continue;
//ensure particles are not from GEANT
if (particle->barcode() >= 2e5) continue;
//check if we have a neutral particle (threeCharge returns int)
if (chargedOnly && HepPID::threeCharge(particle->pdgId()) == 0) continue;
//if we have a photon from hadron decay, and usePhotonsFromHadrons=false, skip this particle
if (!usePhotonsFromHadrons && fromHadronDecay(particle) && particle->pdgId()==22) continue;
//good particle, add to list
selectedlist.push_back(particle);
}
}
bool fromHadronDecay(const xAOD::TruthParticle* particle) const
{
//does this particle come from a hadron decay?
//copied and modified from
//MenuTruthDressingTool::isLeptonFromTau()
const xAOD::TruthVertex* prod = particle->prodVtx();
if(!prod) return false; // no parent.
// Simple loop catch
if (prod==particle->decayVtx()) return false;
// Loop over the parents of this particle.
unsigned int nIncoming = prod->nIncomingParticles();
for(unsigned int itr = 0; itr<nIncoming; ++itr){
int parentId = prod->incomingParticle(itr)->pdgId();
if(HepPID::isHadron(parentId)) {
return true;
}
//note from tancredi: should we search up the chain for leptons as
//well? B->l->photon should be rejected?
if(parentId == particle->pdgId()) {
//continue searching up photon->photon chains
if(fromHadronDecay(prod->incomingParticle(itr))) {
return true;
}
}
}
return false;
}
float calculateDeltaR(const xAOD::IParticle *p1, const xAOD::IParticle *p2) const
{
float phi1 = p1->phi();
float phi2 = p2->phi();
float eta1 = p1->eta();
float eta2 = p2->eta();
float deltaPhi = fabs(phi1-phi2);
if (deltaPhi>TMath::Pi()) deltaPhi = 2.0*TMath::Pi() - deltaPhi;
float deltaPhiSq = (phi1-phi2)*(phi1-phi2);
float deltaEtaSq = (eta1-eta2)*(eta1-eta2);
float deltaR = sqrt(deltaPhiSq+deltaEtaSq);
return deltaR;
}
};
}
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
///////////////////////////////////////////////////////////////////
// GenericTruthThinning.h, (c) ATLAS Detector software
///////////////////////////////////////////////////////////////////
#ifndef DERIVATIONFRAMEWORK_GENERICTRUTHTHINNING_H
#define DERIVATIONFRAMEWORK_GENERICTRUTHTHINNING_H
#include <string>
#include "AthenaBaseComps/AthAlgTool.h"
#include "DerivationFrameworkInterfaces/IThinningTool.h"
#include "DerivationFrameworkMCTruth/DecayGraphHelper.h"
#include "GaudiKernel/ToolHandle.h"
namespace ExpressionParsing {
class ExpressionParser;
}
class IThinningSvc;
namespace DerivationFramework {
class GenericTruthThinning : public AthAlgTool, public IThinningTool {
public:
GenericTruthThinning(const std::string& t, const std::string& n, const IInterface* p);
~GenericTruthThinning();
StatusCode initialize();
StatusCode finalize();
virtual StatusCode doThinning() const;
private:
ServiceHandle<IThinningSvc> m_thinningSvc;
ExpressionParsing::ExpressionParser *m_vertParser;
ExpressionParsing::ExpressionParser *m_partParser;
mutable unsigned int m_ntotvtx, m_ntotpart, m_npassvtx, m_npasspart;
std::string m_particlesKey, m_verticesKey, m_eventsKey;
std::string m_partString;
//std::string m_vtxString;
bool m_preserveDescendants;
bool m_preserveGeneratorDescendants;
bool m_preserveAncestors;
bool m_tauHandling;
};
}
#endif // DERIVATIONFRAMEWORK_GENERICTRUTHTHINNING_H
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
///////////////////////////////////////////////////////////////////
// HardTruthThinning, (c) ATLAS Detector software
// Author: Frank Paige
// ThinningSvc tool to use with CompactHardTruth and truth jets.
// Selects:
// (1) Stable TruthEvent particles matching HardTruth barcodes.
// (2) Constituents of truth jets.
// (3) Selected particles (e.g., B's) and their decays.
// (4) Stable particles within cone around hard leptons/photons.
// Treatment of Geant particles not yet implemented.
//
// Required parameter Type MCGN2 key
// EventInfo EventInfo "EventInfo"
// TruthParticles TruthParticleContainer "TruthParticle"
// TruthVertices TruthVertexContainer key "TruthVertex"
// HardParticles TruthParticleContainer "TruthHardParticle"
//
// Parameters that turn on features:
// JetName If not empty, save constituents with cuts
// KeepIds If pdgId list not empty, save particles and decay chains
// IsolRadius If positive, save stable particles in isolation cones
//
///////////////////////////////////////////////////////////////////
#ifndef DERIVATIONFRAMEWORK_HARDTRUTHTHINNINGTOOL_H
#define DERIVATIONFRAMEWORK_HARDTRUTHTHINNINGTOOL_H
#include <string>
#include "AthenaBaseComps/AthAlgTool.h"
#include "DerivationFrameworkInterfaces/IThinningTool.h"
#include "DerivationFrameworkMCTruth/DecayGraphHelper.h"
#include "xAODTruth/TruthParticleContainer.h"
#include "GaudiKernel/ToolHandle.h"
class IThinningSvc;
namespace DerivationFramework {
class HardTruthThinning : public AthAlgTool, public IThinningTool {
public:
HardTruthThinning(const std::string& t, const std::string& n, const IInterface* p);
~HardTruthThinning();
StatusCode initialize();
StatusCode finalize();
virtual StatusCode doThinning() const;
int getDescendants(const xAOD::TruthParticle* p,
std::vector<const xAOD::TruthParticle*>& d) const;
void printxAODTruth(long long evnum,
const xAOD::TruthParticleContainer* truths) const;
private:
// handle to the thinning service
ServiceHandle<IThinningSvc> m_thinningSvc;
// TruthParticle container names
std::string m_eventInfoName;
std::string m_truthParticleName;
std::string m_truthVertexName;
std::string m_hardParticleName;
// TruthJet name and parameters
std::string m_jetName;
float m_jetPtCut;
float m_jetEtaCut;
float m_jetConstPtCut;
float m_jetPhotonPtCut;
float m_isolR;
float m_isolPtCut;
// Counters
mutable int m_evtCount;
int m_maxCount;
mutable int m_errCount;
// Special particles to keep (with descendants)
std::vector<int> m_keepIds;
};
} //end namespace
#endif // DERIVATIONFRAMEWORK_HARDTRUTHTHINNING_H
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
///////////////////////////////////////////////////////////////////
// MenuTruthThinning, (c) ATLAS Detector software
///////////////////////////////////////////////////////////////////
#ifndef DERIVATIONFRAMEWORK_MENUTRUTHTHINNINGTOOL_H
#define DERIVATIONFRAMEWORK_MENUTRUTHTHINNINGTOOL_H
#include <string>
#include <unordered_set>
#include "AthenaBaseComps/AthAlgTool.h"
#include "DerivationFrameworkInterfaces/IThinningTool.h"
#include "DerivationFrameworkMCTruth/DecayGraphHelper.h"
#include "xAODTruth/TruthParticleContainer.h"
#include "GaudiKernel/ToolHandle.h"
class IThinningSvc;
namespace DerivationFramework {
class MenuTruthThinning : public AthAlgTool, public IThinningTool {
public:
MenuTruthThinning(const std::string& t, const std::string& n, const IInterface* p);
~MenuTruthThinning();
StatusCode initialize();
StatusCode finalize();
virtual StatusCode doThinning() const;
bool isAccepted(const xAOD::TruthParticle*) const;
bool matchHadronIncTau(const xAOD::TruthParticle* part) const;
bool matchQuarkIncTau(const xAOD::TruthParticle* part) const;
bool isOrphanIncTau(const xAOD::TruthParticle* part) const;
bool matchGenParticle(const xAOD::TruthParticle* part,
std::vector<int> &targetIDs, std::vector<int> &intermediateIDs,
bool targetsAreRange) const;
bool isLeptonFromTau(const xAOD::TruthParticle*) const;
bool isFromTau(const xAOD::TruthParticle*) const;
bool isBSM(const xAOD::TruthParticle*) const;
bool isBoson(const xAOD::TruthParticle*) const;
bool isFsrFromLepton(const xAOD::TruthParticle*) const;
private:
// THE MENU
/// Names of the truth container?
std::string m_particlesKey;
std::string m_verticesKey;
std::string m_eventsKey;
/// Parameter: Keep partons?
bool m_writePartons;
/// Parameter: Keep hadrons?
bool m_writeHadrons;
/// Parameter: Keep b-hadrons?
bool m_writeBHadrons;
/// Parameter: Keep geant particles?
bool m_writeGeant;
/// Parameter: Write Geant photons with Pt above this threshold.
/// Set to < 0 to not write any.
float m_geantPhotonPtThresh;
/// Parameter: Keep hadronic tau decays?
bool m_writeTauHad;
/// Parameter: Keep BSM particles?
bool m_writeBSM;
/// Parameter: Keep bosons?
bool m_writeBosons;
/// Parameter: Write partons with Pt above this threshold.