Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

PFCluster truth selector #18106

Merged
merged 9 commits into from Apr 11, 2017
4 changes: 3 additions & 1 deletion RecoEgamma/EgammaMCTools/BuildFile.xml
@@ -1,7 +1,9 @@
<use name="SimDataFormats/Track"/>
<use name="SimDataFormats/Vertex"/>
<use name="clhep"/>
<use name="SimDataFormats/CrossingFrame"/>
<use name="DataFormats/Common"/>
<use name="DataFormats/HepMCCandidate"/>
<use name="clhep"/>
<export>
<lib name="1"/>
</export>
11 changes: 11 additions & 0 deletions RecoEgamma/EgammaMCTools/plugins/BuildFile.xml
@@ -0,0 +1,11 @@
<export>
</export>
<library name="PFClusterMatchedToPhotonsSelector_plugins" file="PFClusterMatchedToPhotonsSelector.cc">
<use name="FWCore/Framework"/>
<use name="FWCore/ParameterSet"/>
<use name="DataFormats/Common"/>
<use name="DataFormats/Math"/>
<use name="DataFormats/ParticleFlowReco"/>
<use name="DataFormats/HepMCCandidate"/>
<use name="SimDataFormats/TrackingAnalysis"/>
</library>
172 changes: 172 additions & 0 deletions RecoEgamma/EgammaMCTools/plugins/PFClusterMatchedToPhotonsSelector.cc
@@ -0,0 +1,172 @@
// -*- C++ -*-
//
// Package: CommonTools/RecoAlgos
// Class: PFClusterMatchedToPhotonsSelector
//
/**\class PFClusterMatchedToPhotonsSelector PFClusterMatchedToPhotonsSelector.cc CommonTools/RecoAlgos/plugins/PFClusterMatchedToPhotonsSelector.cc

Description: Matches ECAL PF clusters to photons that do not convert

Implementation:

*/
//
// Original Author: RCLSA
// Created: Wed, 22 Mar 2017 18:01:40 GMT
//
//


// system include files
#include <memory>
#include <iostream>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/Common/interface/ValueMap.h"
#include "DataFormats/Common/interface/Ref.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
#include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"

typedef reco::PFCluster::EEtoPSAssociation::value_type EEPSPair;
bool sortByKey(const EEPSPair& a, const EEPSPair& b) {
return a.first < b.first;
}

class PFClusterMatchedToPhotonsSelector : public edm::stream::EDProducer<> {
public:
PFClusterMatchedToPhotonsSelector(const edm::ParameterSet&);
static void fillDescriptions(edm::ConfigurationDescriptions&);

private:
virtual void produce(edm::Event&, const edm::EventSetup&);

edm::EDGetTokenT<reco::GenParticleCollection> genParticleToken_; // genParticles
edm::EDGetTokenT<reco::PFClusterCollection> particleFlowClusterECALToken_;
edm::EDGetTokenT<reco::PFCluster::EEtoPSAssociation> associationToken_;
edm::EDGetTokenT<TrackingParticleCollection> trackingParticleToken_;
double matchMaxDR2_;
double matchMaxDEDR2_;

double volumeZ_EB_;
double volumeRadius_EB_;
double volumeZ_EE_;
};

PFClusterMatchedToPhotonsSelector::PFClusterMatchedToPhotonsSelector(const edm::ParameterSet& iConfig) {
//now do what ever initialization is needed
particleFlowClusterECALToken_ = consumes<reco::PFClusterCollection>(iConfig.getParameter<edm::InputTag>("pfClustersTag"));
associationToken_ = consumes<reco::PFCluster::EEtoPSAssociation>(iConfig.getParameter<edm::InputTag>("pfClustersTag"));
trackingParticleToken_ = consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("trackingParticleTag"));
genParticleToken_ = consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("genParticleTag"));

matchMaxDR2_ = iConfig.getParameter<double>("maxDR2");
matchMaxDEDR2_ = iConfig.getParameter<double>("maxDEDR2");
volumeZ_EB_ = iConfig.getParameter<double>("volumeZ_EB");
volumeRadius_EB_ = iConfig.getParameter<double>("volumeRadius_EB");
volumeZ_EE_ = iConfig.getParameter<double>("volumeZ_EE");

produces<reco::PFClusterCollection>();
produces<reco::PFCluster::EEtoPSAssociation>();
produces<edm::ValueMap<reco::GenParticleRef> >();
}


void PFClusterMatchedToPhotonsSelector::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("pfClustersTag", edm::InputTag("particleFlowClusterECAL"));
desc.add<edm::InputTag>("trackingParticleTag", edm::InputTag("mix", "MergedTrackTruth"));
desc.add<edm::InputTag>("genParticleTag", edm::InputTag("genParticles"));
desc.add<double>("maxDR2", 0.1*0.1);
desc.add<double>("maxDEDR2", 0.5*0.5);
desc.add<double>("volumeZ_EB", 304.5);
desc.add<double>("volumeRadius_EB", 123.8);
desc.add<double>("volumeZ_EE", 317.0);
descriptions.add("pfClusterMatchedToPhotonsSelector", desc);

}

void PFClusterMatchedToPhotonsSelector::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {

edm::Handle<reco::PFClusterCollection> particleFlowClusterECALHandle_;
edm::Handle<reco::PFCluster::EEtoPSAssociation> associationHandle_;
edm::Handle<TrackingParticleCollection> trackingParticleHandle_;
edm::Handle<reco::GenParticleCollection> genParticleHandle_;
iEvent.getByToken(particleFlowClusterECALToken_, particleFlowClusterECALHandle_);
iEvent.getByToken(trackingParticleToken_, trackingParticleHandle_);
iEvent.getByToken(genParticleToken_, genParticleHandle_);
iEvent.getByToken(associationToken_, associationHandle_);

std::unique_ptr<reco::PFClusterCollection> out = std::make_unique<reco::PFClusterCollection>();
std::unique_ptr<reco::PFCluster::EEtoPSAssociation> association_out = std::make_unique<reco::PFCluster::EEtoPSAssociation>();
std::unique_ptr<edm::ValueMap<reco::GenParticleRef> > genmatching_out = std::make_unique<edm::ValueMap<reco::GenParticleRef> >();

std::vector<reco::GenParticleRef> genmatching;

size_t iP(0);
size_t iN(0);
for (auto&& pfCluster : *particleFlowClusterECALHandle_) {

bool isMatched = false;
reco::GenParticleRef::key_type matchedKey;

for (auto&& trackingParticle : *trackingParticleHandle_) {
if (trackingParticle.pdgId() != 22) continue;
if (trackingParticle.status() != 1) continue;
matchedKey = trackingParticle.genParticles().at(0).key();
float DR2 = reco::deltaR2(trackingParticle, pfCluster.position());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

CMS code rules ask lower case initials for local variables: what about theDR2 and theDE (just a suggestion)

float DE = 1. - trackingParticle.genParticles().at(0)->energy()/pfCluster.energy();
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it granted that pfCluster.energy() is different from zero (just asking, I don't know)? If there are possible pathological cases in which it can be zero, better to put a thresold...
For performance, better to compute ( 1./pfCluster.energy() ) outside the internal loop once for all, and multiply here that value.

if (DR2 > matchMaxDR2_) continue;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Move this right after the definition of DR2, so that you can avoid computing DE in most cases

if ((DR2 + DE*DE) > matchMaxDEDR2_) continue;

bool isConversion = false;
for (auto&& vertRef : trackingParticle.decayVertices()) {
if (vertRef->position().rho() > volumeRadius_EB_ && std::abs(vertRef->position().z()) < volumeZ_EB_) continue;
if (std::abs(vertRef->position().z()) > volumeZ_EE_) continue;

for(auto&& tpRef: vertRef->daughterTracks()) {
if(std::abs(tpRef->pdgId()) == 11) isConversion = true;
break;
}
if (isConversion) break;
}
if (isConversion) continue;

isMatched = true;
break;
}

if (isMatched) {
out->push_back(pfCluster);
for (size_t i=0; i<associationHandle_.product()->size(); i++) {
if (associationHandle_.product()->at(i).first == iP) {
association_out->push_back(std::make_pair(iN, associationHandle_.product()->at(i).second));
}
}
genmatching.push_back(edm::Ref<reco::GenParticleCollection>(genParticleHandle_,matchedKey));
}
iP++;
}

std::sort(association_out->begin(),association_out->end(),sortByKey);
edm::OrphanHandle<reco::PFClusterCollection> pfClusterHandle= iEvent.put(std::move(out));
iEvent.put(std::move(association_out));

edm::ValueMap<reco::GenParticleRef>::Filler mapFiller(*genmatching_out);
mapFiller.insert(pfClusterHandle, genmatching.begin(), genmatching.end());
mapFiller.fill();
iEvent.put(std::move(genmatching_out));
}

//define this as a plug-in
DEFINE_FWK_MODULE(PFClusterMatchedToPhotonsSelector);
10 changes: 10 additions & 0 deletions RecoEgamma/EgammaMCTools/src/classes.h
@@ -0,0 +1,10 @@
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
#include "DataFormats/Common/interface/ValueMap.h"

namespace DataFormats_PFClusterMatch {
struct dictionary {
edm::ValueMap<edm::Ref<std::vector<reco::GenParticle>,reco::GenParticle,edm::refhelper::FindUsingAdvance<std::vector<reco::GenParticle>,reco::GenParticle> > > vmg2;
edm::Wrapper<edm::ValueMap<edm::Ref<std::vector<reco::GenParticle>,reco::GenParticle,edm::refhelper::FindUsingAdvance<std::vector<reco::GenParticle>,reco::GenParticle> > > > wvmg2;
};
}
4 changes: 4 additions & 0 deletions RecoEgamma/EgammaMCTools/src/classes_def.xml
@@ -0,0 +1,4 @@
<lcgdict>
<class name="edm::ValueMap<edm::Ref<std::vector<reco::GenParticle>,reco::GenParticle,edm::refhelper::FindUsingAdvance<std::vector<reco::GenParticle>,reco::GenParticle> > >" />
<class name="edm::Wrapper<edm::ValueMap<edm::Ref<std::vector<reco::GenParticle>,reco::GenParticle,edm::refhelper::FindUsingAdvance<std::vector<reco::GenParticle>,reco::GenParticle> > > >" />
</lcgdict>
113 changes: 113 additions & 0 deletions RecoEgamma/EgammaMCTools/test/pfClusterForCalibration.py
@@ -0,0 +1,113 @@
# EGM skimmer
# Author: Rafael Lopes de Sa

import FWCore.ParameterSet.Config as cms

# Run with the 2017 detector
from Configuration.StandardSequences.Eras import eras
process = cms.Process('SKIM',eras.Run2_2017)

# Import the standard packages for reconstruction and digitization
process.load('Configuration.StandardSequences.Services_cff')
process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi')
process.load('FWCore.MessageService.MessageLogger_cfi')
process.load('Configuration.EventContent.EventContent_cff')
process.load('SimGeneral.MixingModule.mixNoPU_cfi')
process.load('Configuration.StandardSequences.Digi_cff')
process.load('Configuration.StandardSequences.GeometryRecoDB_cff')
process.load('Configuration.StandardSequences.MagneticField_cff')
process.load('Configuration.StandardSequences.RawToDigi_cff')
process.load('Configuration.StandardSequences.L1Reco_cff')
process.load('Configuration.StandardSequences.Reconstruction_cff')
process.load('Configuration.StandardSequences.EndOfProcess_cff')
process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff')
process.load('RecoEgamma.EgammaMCTools.pfClusterMatchedToPhotonsSelector_cfi')

# Global Tag configuration ... just using the same as in the RelVal
from Configuration.AlCa.GlobalTag import GlobalTag
process.GlobalTag = GlobalTag(process.GlobalTag, '81X_upgrade2017_realistic_v26', '')

process.MessageLogger.cerr.threshold = 'ERROR'
process.MessageLogger.cerr.FwkReport.reportEvery = 1000

process.options = cms.untracked.PSet( allowUnscheduled = cms.untracked.bool(True) )

# This is where users have some control.
# Define which collections to save and which dataformat we are using
savedCollections = cms.untracked.vstring('drop *',
# The commented ones are large collections that can be kept for debug
# 'keep EcalRecHitsSorted_*_*_*',
# 'keep recoPFClusters_*_*_*',
# 'keep recoCaloClusters_*_*_*',
# 'keep recoSuperClusters_*_*_*',
# 'keep recoGsfElectron*_*_*_*',
# 'keep recoPhoton*_*_*_*',
# 'keep *_mix_MergedTrackTruth_*',
'keep *_reducedEcalRecHits*_*_*',
'keep double_fixedGridRho*_*_*',
'keep recoGenParticles_*_*_*',
'keep GenEventInfoProduct_*_*_*',
'keep PileupSummaryInfos_*_*_*',
'keep *_ecalDigis_*_*',
'keep *_offlinePrimaryVertices_*_*',
'keep *_particleFlowCluster*_*_*')

process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(15))

process.source = cms.Source("PoolSource",
fileNames = cms.untracked.vstring(
'/store/mc/PhaseIFall16DR/GluGluHToGG_M-125_13TeV_powheg_pythia8/AODSIM/FlatPU28to62HcalNZSRAW_81X_upgrade2017_realistic_v26-v1/100000/005AB6CE-27ED-E611-98CA-E0071B7A8590.root'
),
secondaryFileNames = cms.untracked.vstring(
'/store/mc/PhaseIFall16DR/GluGluHToGG_M-125_13TeV_powheg_pythia8/GEN-SIM-RAW/FlatPU28to62HcalNZSRAW_81X_upgrade2017_realistic_v26-v1/100000/0416D6B7-04ED-E611-B342-E0071B7A8550.root',
'/store/mc/PhaseIFall16DR/GluGluHToGG_M-125_13TeV_powheg_pythia8/GEN-SIM-RAW/FlatPU28to62HcalNZSRAW_81X_upgrade2017_realistic_v26-v1/100000/14829DD8-04ED-E611-8049-A0000420FE80.root',
'/store/mc/PhaseIFall16DR/GluGluHToGG_M-125_13TeV_powheg_pythia8/GEN-SIM-RAW/FlatPU28to62HcalNZSRAW_81X_upgrade2017_realistic_v26-v1/100000/54AFE9C4-04ED-E611-952D-A0000420FE80.root',
'/store/mc/PhaseIFall16DR/GluGluHToGG_M-125_13TeV_powheg_pythia8/GEN-SIM-RAW/FlatPU28to62HcalNZSRAW_81X_upgrade2017_realistic_v26-v1/100000/5A32C6B9-04ED-E611-B1EB-E0071B7A8550.root',
'/store/mc/PhaseIFall16DR/GluGluHToGG_M-125_13TeV_powheg_pythia8/GEN-SIM-RAW/FlatPU28to62HcalNZSRAW_81X_upgrade2017_realistic_v26-v1/100000/60E162B8-04ED-E611-898D-E0071B7A58F0.root',
'/store/mc/PhaseIFall16DR/GluGluHToGG_M-125_13TeV_powheg_pythia8/GEN-SIM-RAW/FlatPU28to62HcalNZSRAW_81X_upgrade2017_realistic_v26-v1/100000/6A47DD1A-FEEC-E611-81EB-A0000420FE80.root',
'/store/mc/PhaseIFall16DR/GluGluHToGG_M-125_13TeV_powheg_pythia8/GEN-SIM-RAW/FlatPU28to62HcalNZSRAW_81X_upgrade2017_realistic_v26-v1/100000/92B923B6-04ED-E611-9DC9-24BE05C48821.root',
'/store/mc/PhaseIFall16DR/GluGluHToGG_M-125_13TeV_powheg_pythia8/GEN-SIM-RAW/FlatPU28to62HcalNZSRAW_81X_upgrade2017_realistic_v26-v1/100000/B40E77B4-04ED-E611-9E30-E0071B7A45D0.root',
'/store/mc/PhaseIFall16DR/GluGluHToGG_M-125_13TeV_powheg_pythia8/GEN-SIM-RAW/FlatPU28to62HcalNZSRAW_81X_upgrade2017_realistic_v26-v1/100000/C48157B5-04ED-E611-BEC1-E0071B7A45D0.root',
'/store/mc/PhaseIFall16DR/GluGluHToGG_M-125_13TeV_powheg_pythia8/GEN-SIM-RAW/FlatPU28to62HcalNZSRAW_81X_upgrade2017_realistic_v26-v1/100000/CAED3A16-FEEC-E611-8262-24BE05CEFB41.root'
)
)
process.PFCLUSTERoutput = cms.OutputModule("PoolOutputModule",
dataset = cms.untracked.PSet(dataTier = cms.untracked.string('RECO'),
filterName = cms.untracked.string('')
),
eventAutoFlushCompressedSize = cms.untracked.int32(5242880),
fileName = cms.untracked.string('skimEGMobjects_fromRAW.root'),
outputCommands = savedCollections,
splitLevel = cms.untracked.int32(0)
)

# Run the digitizer to make the trackingparticles
process.mix.digitizers = cms.PSet(process.theDigitizersValid)
process.trackingtruth_step = cms.Path(process.pdigi_valid)

# Remake the PFClusters
process.pfclusters_step = cms.Path(process.bunchSpacingProducer *
process.ecalDigis *
process.ecalPreshowerDigis *
process.ecalPreshowerRecHit *
process.ecalMultiFitUncalibRecHit *
process.ecalDetIdToBeRecovered *
process.ecalRecHit *
process.particleFlowRecHitPS *
process.particleFlowRecHitECAL *
process.particleFlowClusterECALUncorrected *
process.particleFlowClusterPS *
process.particleFlowClusterECAL)

# Select the PFClusters we want to calibrate
process.particleFlowClusterECALMatchedToPhotons = process.pfClusterMatchedToPhotonsSelector.clone()
process.selection_step = cms.Path(process.particleFlowClusterECALMatchedToPhotons)

# Ends job and writes our output
process.endjob_step = cms.EndPath(process.endOfProcess)
process.output_step = cms.EndPath(process.PFCLUSTERoutput)

# Schedule definition, rebuilding rechits
process.schedule = cms.Schedule(process.trackingtruth_step,process.pfclusters_step,process.selection_step,process.endjob_step,process.output_step)