Skip to content
Snippets Groups Projects
Commit a74747a7 authored by Scott Snyder's avatar Scott Snyder Committed by scott snyder
Browse files

GeneratorFilters: const fixes.

DataVector has had a long-standing bug in which it has been possible
to get a non-const pointer from a const_iterator.
Fix instances where this package was relying on that to compile.

The class LeadingDiBjetFilter is further problematic as it was taking
advantage of this to modify the weights in an existing GenEvent object
in StoreGate.  However, there's nothing currently using this class,
so just remove it.
parent b8ee1a66
No related branches found
No related tags found
No related merge requests found
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
// --------------------------------------------------
//
// File: GeneratorFilters/LeadingDiBjetFilter.h
// AuthorList:
// Tim Barklow
//
// Based on Filter DiBjetFilter written by:
// Stephen Bieniek
//
#ifndef GENERATORFILTERSLEADINGDIBJETFILTER_H
#define GENERATORFILTERSLEADINGDIBJETFILTER_H
#include "GeneratorModules/GenFilter.h"
//Random number generator required for accepting light jets
class TRandom3;
class LeadingDiBjetFilter:public GenFilter {
public:
LeadingDiBjetFilter(const std::string& name, ISvcLocator* pSvcLocator);
virtual ~LeadingDiBjetFilter();
virtual StatusCode filterInitialize();
virtual StatusCode filterFinalize();
virtual StatusCode filterEvent();
private:
// Setable Properties:-
/* Variables for heavy flavour filtering */
double m_bottomPtMin;
double m_bottomEtaMax;
double m_deltaRFromTruth;
double m_jetPtMin;
double m_jetEtaMax;
/* Variables for cutting sample into pt slices */
double m_leadJet_ptMin;
double m_leadJet_ptMax;
std::string m_TruthJetContainerName;
/* Variables for light jet suppression */
double m_LightJetSuppressionFactor;
bool m_AcceptSomeLightEvents;
// Local Member Data:-
int m_NPass;
int m_Nevt;
double m_SumOfWeigths_Pass;
double m_SumOfWeigths_Evt;
TRandom3* m_ranNumGen;
bool isBwithWeakDK(const int pID) const;
};
#endif
......@@ -125,9 +125,7 @@ StatusCode DiBjetFilter::filterEvent() {
int bJetCounter = 0;
double weight = 1;
McEventCollection::const_iterator itr;
for (itr = events()->begin(); itr!=events()->end(); ++itr) {
const HepMC::GenEvent* genEvt = (*itr);
for (const HepMC::GenEvent* genEvt : *events()) {
weight = genEvt->weights().front();
HepMC::GenEvent::particle_const_iterator pitr;
std::vector< HepMC::GenEvent::particle_const_iterator > bHadrons;
......@@ -157,8 +155,7 @@ StatusCode DiBjetFilter::filterEvent() {
&& m_ranNumGen->Uniform() < (1.0 / m_LightJetSuppressionFactor )
&& passLeadJetCut ){
/* Modify event weight to account for light jet prescale */
for (itr = events()->begin(); itr!=events()->end(); ++itr) {
HepMC::GenEvent* genEvt = (*itr);
for (HepMC::GenEvent* genEvt : *events()) {
genEvt->weights().front() *= m_LightJetSuppressionFactor;
}
pass = true;
......
......@@ -247,10 +247,7 @@ GapJetFilter::filterEvent()
Int_t Clustag=0;
// Loop over all events in McEventCollection
McEventCollection::const_iterator itr;
for (itr = events()->begin(); itr != events()->end(); ++itr) {
HepMC::GenEvent* genEvt = *itr;
for (HepMC::GenEvent* genEvt : *events()) {
// Loop over all particles in event
HepMC::GenEvent::particle_const_iterator pitr;
......
......@@ -153,9 +153,7 @@ StatusCode JetForwardFilter::filterEvent() {
if (flagNJets != 0 && flag1stJet != 0 && flag2ndJet != 0 && flagJJ != 0 && flagW != 0) {
ATH_MSG_DEBUG("Pass filter");
// If passed then add a weight!
McEventCollection::const_iterator itr;
for (itr = events()->begin(); itr!=events()->end(); ++itr) {
HepMC::GenEvent* genEvt = (*itr);
for (HepMC::GenEvent* genEvt : *events()) {
// Store the inverse of the weighting for conventions sake
genEvt->weights().push_back(1.0/weighting);
}
......
......@@ -253,9 +253,7 @@ StatusCode JetIntervalFilter::filterEvent() {
ATH_MSG_DEBUG("Pass filter");
// If passed then calculate (and add) a weight!
if (m_weightingEvents) {
McEventCollection::const_iterator itr;
for (itr = events()->begin(); itr!=events()->end(); ++itr) {
HepMC::GenEvent* genEvt = *itr;
for (HepMC::GenEvent* genEvt : *events()) {
// Simple test to see if it is possible to put extra weights into Monte Carlo!
// Store the inverse of the weighting
genEvt->weights().push_back(1.0/weighting);
......
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
// -------------------------------------------------------------
// File: GeneratorFilters/LeadingDiBjetFilter.cxx
// Description:
// This is a modification of the DiBjetFilter filter. It differs in that
// b-hadrons must be matched to the two leading pT jets in order to be counted.
// Otherwise it is the same.
// It cuts on:
// - The presence of two b-hadron matched to the two leading pT jets.
// - The pT of the leading jet
// - It can also pass light jets and reweight the event accordingly to make an
// statisticaly inclusive sample of jets but with baised towards
// heavy flavour events
//
// AuthorList:
// S. Bieniek
// Header for this module:-
#include "GeneratorFilters/LeadingDiBjetFilter.h"
// Other classes used by this class:-
#include <math.h>
#include "GaudiKernel/SystemOfUnits.h"
#include "xAODJet/JetContainer.h"
#include "McParticleEvent/TruthParticle.h"
#include "CxxUtils/BasicTypes.h"
#include "TRandom3.h"
#include "TLorentzVector.h"
using HepMC::GenVertex;
using HepMC::GenParticle;
//--------------------------------------------------------------------------
LeadingDiBjetFilter::LeadingDiBjetFilter(const std::string& name,
ISvcLocator* pSvcLocator): GenFilter(name,pSvcLocator) {
//--------------------------------------------------------------------------
// Local Member Data:-
declareProperty("LeadJetPtMin",m_leadJet_ptMin=0*CLHEP::GeV);
declareProperty("LeadJetPtMax",m_leadJet_ptMax=50000*CLHEP::GeV);
declareProperty("BottomPtMin",m_bottomPtMin=5.0*CLHEP::GeV);
declareProperty("BottomEtaMax",m_bottomEtaMax=3.0);
declareProperty("JetPtMin",m_jetPtMin=15.0*CLHEP::GeV);
declareProperty("JetEtaMax",m_jetEtaMax=2.7);
declareProperty("DeltaRFromTruth",m_deltaRFromTruth=0.3);
declareProperty("TruthContainerName",m_TruthJetContainerName="AntiKt4TruthJets");
declareProperty("LightJetSuppressionFactor",m_LightJetSuppressionFactor=10);
declareProperty("AcceptSomeLightEvents",m_AcceptSomeLightEvents=false);
m_NPass = 0;
m_Nevt = 0;
m_SumOfWeigths_Pass = 0;
m_SumOfWeigths_Evt = 0;
m_ranNumGen = 0;
}
//--------------------------------------------------------------------------
LeadingDiBjetFilter::~LeadingDiBjetFilter(){
//--------------------------------------------------------------------------
}
//---------------------------------------------------------------------------
StatusCode LeadingDiBjetFilter::filterInitialize() {
//---------------------------------------------------------------------------
m_Nevt = 0;
m_NPass = 0;
m_SumOfWeigths_Pass = 0;
m_SumOfWeigths_Evt = 0;
if(m_AcceptSomeLightEvents){
m_ranNumGen = new TRandom3();
/* Set seed with respect to computer clock time */
m_ranNumGen->SetSeed(0);
}
ATH_MSG_INFO( "Initialized" );
return StatusCode::SUCCESS;
}
//---------------------------------------------------------------------------
StatusCode LeadingDiBjetFilter::filterFinalize() {
//---------------------------------------------------------------------------
if(m_AcceptSomeLightEvents){
delete m_ranNumGen;
}
ATH_MSG_INFO( m_NPass << " Events out of " << m_Nevt << " passed the filter" );
ATH_MSG_INFO( m_SumOfWeigths_Pass << " out of " << m_SumOfWeigths_Evt << " SumOfWeights counter, passed/total" );
return StatusCode::SUCCESS;
}
//---------------------------------------------------------------------------
StatusCode LeadingDiBjetFilter::filterEvent() {
//---------------------------------------------------------------------------
bool pass = false;
m_Nevt++;
// if(m_Nevt > 420) log << MSG::INFO << " m_Nevt= " << m_Nevt << " filterEvent point 001" << endmsg;
const xAOD::JetContainer* truthjetTES = 0;
// StatusCode sc=m_sgSvc->retrieve( truthjetTES, m_TruthJetContainerName);
StatusCode sc=evtStore()->retrieve( truthjetTES, m_TruthJetContainerName);
// if(m_Nevt > 420) log << MSG::INFO << " m_Nevt= " << m_Nevt << " filterEvent point 002" << endmsg;
if( sc.isFailure() || !truthjetTES ) {
ATH_MSG_WARNING( "No xAOD::JetContainer found in TDS " << m_TruthJetContainerName \
<< sc.isFailure() << " "<< !truthjetTES
);
return StatusCode::SUCCESS;
}
// if(m_Nevt > 420) log << MSG::INFO << " m_Nevt= " << m_Nevt << " filterEvent point 003" << endmsg;
bool passLeadJetCut = false;
xAOD::JetContainer::const_iterator jitr;
double lead_jet_pt = 0.0;
double lead_2nd_jet_pt = 0.0;
std::vector<xAOD::JetContainer::const_iterator> jets;
// if(m_Nevt > 420) log << MSG::INFO << " m_Nevt= " << m_Nevt << " filterEvent point 004" << endmsg;
for (jitr = (*truthjetTES).begin(); jitr !=(*truthjetTES).end(); ++jitr) {
if( (*jitr)->pt() < m_jetPtMin ) continue;
if( fabs( (*jitr)->eta() ) > m_jetEtaMax ) continue;
if(jets.size() == 0) {
jets.push_back(jitr);
}
else if(jets.size() ==1) {
if((*jitr)->pt()>(*jets[0])->pt()) {
jets.push_back(jets[0]);
jets[0]=jitr;
}
else {
jets.push_back(jitr);
}
}
else {
if((*jitr)->pt()>(*jets[0])->pt()) {
jets[1]=jets[0];
jets[0]=jitr;
}
else if((*jitr)->pt()>(*jets[1])->pt()) {
jets[1]=jitr;
}
}
if( (*jitr)->pt() > lead_jet_pt ){
lead_jet_pt = (*jitr)->pt();
}
}
// if(m_Nevt > 420) log << MSG::INFO << " m_Nevt= " << m_Nevt << " filterEvent point 005" << endmsg;
if(jets.size() > 0) lead_jet_pt = (*jets[0])->pt();
if(jets.size() > 1) lead_2nd_jet_pt=(*jets[1])->pt();
if( lead_jet_pt > m_leadJet_ptMin && lead_jet_pt < m_leadJet_ptMax ) passLeadJetCut = true;
// if(m_Nevt > 420) log << MSG::INFO << " m_Nevt= " << m_Nevt << " filterEvent point 006" << endmsg;
int bJetCounter = 0;
double weight = 1;
int n_events_in_collection=0;
int n_bHadrons_total=0;
// if(m_Nevt > 420) log << MSG::INFO << " m_Nevt= " << m_Nevt << " filterEvent point 007" << endmsg;
McEventCollection::const_iterator itr;
// if(m_Nevt > 420) log << MSG::INFO << " m_Nevt= " << m_Nevt << " filterEvent point 008" << endmsg;
for (itr = events()->begin(); itr!=events()->end(); ++itr) {
// if(m_Nevt > 420) log << MSG::INFO << " m_Nevt= " << m_Nevt << " filterEvent point 009" << endmsg;
n_events_in_collection++;
const HepMC::GenEvent* genEvt = (*itr);
weight = genEvt->weights().front();
HepMC::GenEvent::particle_const_iterator pitr;
std::vector< HepMC::GenEvent::particle_const_iterator > bHadrons;
// if(m_Nevt > 420) log << MSG::INFO << " m_Nevt= " << m_Nevt << " filterEvent point 010" << endmsg;
for(pitr=genEvt->particles_begin(); pitr!=genEvt->particles_end(); ++pitr ) {
if( !isBwithWeakDK( (*pitr)->pdg_id()) ) continue;
if( (*pitr)->momentum().perp() < m_bottomPtMin ) continue;
if( fabs( (*pitr)->momentum().pseudoRapidity() ) > m_bottomEtaMax) continue;
bHadrons.push_back(pitr);
}
// if(m_Nevt > 420) log << MSG::INFO << " m_Nevt= " << m_Nevt << " filterEvent point 011" << endmsg;
n_bHadrons_total += bHadrons.size();
for(uint i = 0; i < jets.size(); i++){
for(uint j = 0; j < bHadrons.size(); j++){
HepMC::FourVector tmp = (*bHadrons[j])->momentum();
TLorentzVector genpart(tmp.x(), tmp.y(), tmp.z(), tmp.t());
double dR = (*jets[i])->p4().DeltaR(genpart);
if(dR<m_deltaRFromTruth){
bJetCounter++;
break;
}
}
}
}
// if(m_Nevt > 420) log << MSG::INFO << " m_Nevt= " << m_Nevt << " filterEvent point 012" << endmsg;
m_SumOfWeigths_Evt += weight;
pass = (bJetCounter > 1) && passLeadJetCut;
if( (bJetCounter < 2) && m_AcceptSomeLightEvents
&& m_ranNumGen->Uniform() < (1.0 / m_LightJetSuppressionFactor )
&& passLeadJetCut ){
/* Modify event weight to account for light jet prescale */
for (itr = events()->begin(); itr!=events()->end(); ++itr) {
HepMC::GenEvent* genEvt = (*itr);
genEvt->weights().front() *= m_LightJetSuppressionFactor;
}
pass = true;
}
// if(m_Nevt > 420) log << MSG::INFO << " m_Nevt= " << m_Nevt << " filterEvent point 013" << endmsg;
if(pass){
m_NPass++;
m_SumOfWeigths_Pass += weight;
}
if(m_Nevt < 20 || m_Nevt%100 == 0 || pass) {
ATH_MSG_INFO( " m_Nevt= " << m_Nevt << " n_events_in_collection= " << n_events_in_collection << " 1st,2nd lead_jet_pt= " << lead_jet_pt << " " << lead_2nd_jet_pt << " passLeadJetCut= " << passLeadJetCut
<< " n_bHadrons_total= " << n_bHadrons_total << " bJetCounter= " << bJetCounter << " pass= " << pass << " m_NPass= " << m_NPass );
}
setFilterPassed(pass);
return StatusCode::SUCCESS;
}
bool LeadingDiBjetFilter::isBwithWeakDK(const int pID) const
{
int id = abs(pID);
if ( id == 511 || // B+
id == 521 || // B0
id == 531 || // Bs
id == 541 || // Bc
id == 5122 || // Lambda_B
id == 5132 || // Xi_b-
id == 5232 || // Xi_b0
id == 5332 ) // Omega_b-
return true;
else
return false;
}
......@@ -382,7 +382,7 @@ StatusCode TruthJetWeightFilter::saveWeight(double weight) {
// Prepare write access to put weights
McEventCollection* mymccoll;
CHECK(evtStore()->retrieve(mymccoll, "GEN_EVENT"));
McEventCollection::const_iterator itr = mymccoll->begin();
McEventCollection::iterator itr = mymccoll->begin();
for (int ii = 0; itr!=mymccoll->end(); ++itr, ii++) {
HepMC::GenEvent* genEvt = (*itr);
HepMC::WeightContainer& wgtsC = genEvt->weights();
......
......@@ -68,7 +68,6 @@
#include "GeneratorFilters/JetFilterWithTruthPhoton.h"
#include "GeneratorFilters/Boosted2DijetFilter.h"
#include "GeneratorFilters/DiBjetFilter.h"
#include "GeneratorFilters/LeadingDiBjetFilter.h"
#include "GeneratorFilters/HiggsFilter.h"
#include "GeneratorFilters/TTbarBoostCatFilter.h"
#include "GeneratorFilters/MultiParticleFilter.h"
......@@ -153,7 +152,6 @@ DECLARE_COMPONENT( VBFMjjIntervalFilter )
DECLARE_COMPONENT( JetFilterWithTruthPhoton )
DECLARE_COMPONENT( Boosted2DijetFilter )
DECLARE_COMPONENT( DiBjetFilter )
DECLARE_COMPONENT( LeadingDiBjetFilter )
DECLARE_COMPONENT( HiggsFilter )
DECLARE_COMPONENT( TTbarBoostCatFilter )
DECLARE_COMPONENT( MultiParticleFilter )
......@@ -167,4 +165,3 @@ DECLARE_COMPONENT( DecaysFinalStateFilter )
DECLARE_COMPONENT( HTFilter )
DECLARE_COMPONENT( MissingEtFilter )
DECLARE_COMPONENT( TrimuMassRangeFilter )
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