Skip to content

Commit

Permalink
Migrate to reco::GenParticle for real
Browse files Browse the repository at this point in the history
  • Loading branch information
@gpetruc
gpetruc committed Aug 28, 2013
1 parent 631bb66 commit a5c8f17
Show file tree
Hide file tree
Showing 4 changed files with 96 additions and 87 deletions.
85 changes: 27 additions & 58 deletions MuonAnalysis/MuonAssociators/plugins/MuonMCClassifier.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,27 +116,13 @@ class MuonMCClassifier : public edm::EDProducer {
}
}

const HepMC::GenParticle * getGpMother(const HepMC::GenParticle *gp) {
if (gp != 0) {
const HepMC::GenVertex *vtx = gp->production_vertex();
if (vtx != 0 && vtx->particles_in_size() > 0) {
return *vtx->particles_in_const_begin();
}
}
return 0;
}

/// Find the index of a genParticle given it's barcode. -1 if not found
int fetch(const edm::Handle<std::vector<int> > & genBarcodes, int barcode) const;

/// Convert TrackingParticle into GenParticle, save into output collection,
/// if mother is primary set reference to it,
/// return index in output collection
int convertAndPush(const TrackingParticle &tp,
reco::GenParticleCollection &out,
const TrackingParticleRef &momRef,
const edm::Handle<reco::GenParticleCollection> & genParticles,
const edm::Handle<std::vector<int> > & genBarcodes) const ;
const edm::Handle<reco::GenParticleCollection> & genParticles) const ;
};

MuonMCClassifier::MuonMCClassifier(const edm::ParameterSet &iConfig) :
Expand Down Expand Up @@ -196,10 +182,8 @@ MuonMCClassifier::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
iEvent.getByLabel(trackingParticles_, trackingParticles);

edm::Handle<reco::GenParticleCollection> genParticles;
edm::Handle<std::vector<int> > genBarcodes;
if (linkToGenParticles_) {
iEvent.getByLabel(genParticles_, genParticles);
iEvent.getByLabel(genParticles_, genBarcodes);
}

edm::ESHandle<TrackAssociatorBase> associatorBase;
Expand Down Expand Up @@ -315,32 +299,25 @@ MuonMCClassifier::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
// Try to extract mother and grand mother of this muon.
// Unfortunately, SIM and GEN histories require diffent code :-(
if (!tp->genParticles().empty()) { // Muon is in GEN
#warning "This file has been modified just to get it to compile without any regard as to whether it still functions as intended"
#ifdef REMOVED_JUST_TO_GET_IT_TO_COMPILE__THIS_CODE_NEEDS_TO_BE_CHECKED
const HepMC::GenParticle * genMom = getGpMother(tp->genParticle()[0].get());
#else
const HepMC::GenParticle * genMom = NULL;
#endif
if (genMom) {
momPdgId[i] = tp->pdgId();
momStatus[i] = tp->status();
if (genMom->production_vertex()) {
const HepMC::ThreeVector & momVtx = genMom->production_vertex()->point3d();
momRho[i] = momVtx.perp() * 0.1; momZ[i] = momVtx.z() * 0.1; // HepMC is in mm!
}
reco::GenParticleRef genp = tp->genParticles()[0];
reco::GenParticleRef genMom = genp->numberOfMothers() > 0 ? genp->motherRef() : reco::GenParticleRef();
if (genMom.isNonnull()) {
momPdgId[i] = genMom->pdgId();
momStatus[i] = genMom->status();
momRho[i] = genMom->vertex().Rho(); momZ[i] = genMom->vz();
edm::LogVerbatim("MuonMCClassifier") << "\t Particle pdgId = "<<hitsPdgId[i] << " produced at rho = " << prodRho[i] << ", z = " << prodZ[i] << ", has GEN mother pdgId = " << momPdgId[i];
const HepMC::GenParticle * genGMom = getGpMother(genMom);
if (genGMom) {
gmomPdgId[i] = genGMom->pdg_id();
reco::GenParticleRef genGMom = genMom->numberOfMothers() > 0 ? genMom->motherRef() : reco::GenParticleRef();
if (genGMom.isNonnull()) {
gmomPdgId[i] = genGMom->pdgId();
edm::LogVerbatim("MuonMCClassifier") << "\t\t mother prod. vertex rho = " << momRho[i] << ", z = " << momZ[i] << ", grand-mom pdgId = " << gmomPdgId[i];
}
// in this case, we might want to know the heaviest mom flavour
for (const HepMC::GenParticle *nMom = genMom;
nMom != 0 && abs(nMom->pdg_id()) >= 100; // stop when we're no longer looking at hadrons or mesons
nMom = getGpMother(nMom)) {
int flav = flavour(nMom->pdg_id());
for (reco::GenParticleRef nMom = genMom;
nMom.isNonnull() && abs(nMom->pdgId()) >= 100; // stop when we're no longer looking at hadrons or mesons
nMom = nMom->numberOfMothers() > 0 ? nMom->motherRef() : reco::GenParticleRef()) {
int flav = flavour(nMom->pdgId());
if (hmomFlav[i] < flav) hmomFlav[i] = flav;
edm::LogVerbatim("MuonMCClassifier") << "\t\t backtracking flavour: mom pdgId = "<<nMom->pdg_id()<< ", flavour = " << flav << ", heaviest so far = " << hmomFlav[i];
edm::LogVerbatim("MuonMCClassifier") << "\t\t backtracking flavour: mom pdgId = "<<nMom->pdgId()<< ", flavour = " << flav << ", heaviest so far = " << hmomFlav[i];
}
}
} else { // Muon is in SIM Only
Expand All @@ -353,8 +330,8 @@ MuonMCClassifier::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
", has SIM mother pdgId = " << momPdgId[i] << " produced at rho = " << simMom->vertex().Rho() << ", z = " << simMom->vertex().Z();
if (!simMom->genParticles().empty()) {
momStatus[i] = simMom->genParticles()[0]->status();
#warning "This file has been modified just to get it to compile without any regard as to whether it still functions as intended"
gmomPdgId[i] = simMom->pdgId();
reco::GenParticleRef genGMom = (simMom->genParticles()[0]->numberOfMothers() > 0 ? simMom->genParticles()[0]->motherRef() : reco::GenParticleRef());
if (genGMom.isNonnull()) gmomPdgId[i] = genGMom->pdgId();
edm::LogVerbatim("MuonMCClassifier") << "\t\t SIM mother is in GEN (status " << momStatus[i] << "), grand-mom id = " << gmomPdgId[i];
} else {
momStatus[i] = -1;
Expand Down Expand Up @@ -418,14 +395,14 @@ MuonMCClassifier::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
if (linkToGenParticles_ && abs(ext[i]) >= 2) {
// Link to the genParticle if possible, but not decays in flight (in ppMuX they're in GEN block, but they have wrong parameters)
if (!tp->genParticles().empty() && abs(ext[i]) >= 5) {
#warning "This file has been modified just to get it to compile without any regard as to whether it still functions as intended"
#ifdef REMOVED_JUST_TO_GET_IT_TO_COMPILE__THIS_CODE_NEEDS_TO_BE_CHECKED
muToPrimary[i] = fetch(genBarcodes, tp->genParticle()[0]->barcode());
#endif
if (genParticles.id() != tp->genParticles().id()) {
throw cms::Exception("Configuration") << "Product ID mismatch between the genParticle collection (" << genParticles_ << ", id " << genParticles.id() << ") and the references in the TrackingParticles (id " << tp->genParticles().id() << ")\n";
}
muToPrimary[i] = tp->genParticles()[0].key();
} else {
// Don't put the same trackingParticle twice!
int &indexPlus1 = tpToSecondaries[tp]; // will create a 0 if the tp is not in the list already
if (indexPlus1 == 0) indexPlus1 = convertAndPush(*tp, *secondaries, getTpMother(tp), genParticles, genBarcodes) + 1;
if (indexPlus1 == 0) indexPlus1 = convertAndPush(*tp, *secondaries, getTpMother(tp), genParticles) + 1;
muToSecondary[i] = indexPlus1 - 1;
}
}
Expand Down Expand Up @@ -496,26 +473,18 @@ MuonMCClassifier::flavour(int pdgId) const {
return 0;
}

int MuonMCClassifier::fetch(const edm::Handle<std::vector<int> > & genBarcodes, int barcode) const
{
std::vector<int>::const_iterator it = std::find(genBarcodes->begin(), genBarcodes->end(), barcode);
return (it == genBarcodes->end()) ? -1 : (it - genBarcodes->begin());
}

// push secondary in collection.
// if it has a primary mother link to it.
int MuonMCClassifier::convertAndPush(const TrackingParticle &tp,
reco::GenParticleCollection &out,
const TrackingParticleRef & simMom,
const edm::Handle<reco::GenParticleCollection> & genParticles,
const edm::Handle<std::vector<int> > & genBarcodes) const {
const edm::Handle<reco::GenParticleCollection> & genParticles) const {
out.push_back(reco::GenParticle(tp.charge(), tp.p4(), tp.vertex(), tp.pdgId(), tp.status(), true));
if (simMom.isNonnull() && !simMom->genParticles().empty()) {
#warning "This file has been modified just to get it to compile without any regard as to whether it still functions as intended"
#ifdef REMOVED_JUST_TO_GET_IT_TO_COMPILE__THIS_CODE_NEEDS_TO_BE_CHECKED
int momIdx = fetch(genBarcodes, simMom->genParticle()[0]->barcode());
if (momIdx != -1) out.back().addMother(reco::GenParticleRef(genParticles, momIdx));
#endif
if (genParticles.id() != simMom->genParticles().id()) {
throw cms::Exception("Configuration") << "Product ID mismatch between the genParticle collection (" << genParticles_ << ", id " << genParticles.id() << ") and the references in the TrackingParticles (id " << simMom->genParticles().id() << ")\n";
}
out.back().addMother(simMom->genParticles()[0]);
}
return out.size()-1;
}
Expand Down
21 changes: 18 additions & 3 deletions MuonAnalysis/MuonAssociators/python/muonClassificationByHits_cfi.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,15 +5,30 @@
# SimMuon/MCTruth V02-05-00-03 (35X) or V02-06-00+ (37X+)

from SimGeneral.MixingModule.mixNoPU_cfi import *
from SimGeneral.TrackingAnalysis.trackingParticlesNoSimHits_cfi import *
trackingParticlesNoSimHits = mix.clone(
digitizers = cms.PSet(
mergedtruth = mix.digitizers.mergedtruth.clone(
simHitCollections = cms.PSet(
pixel = cms.VInputTag(),
tracker = cms.VInputTag(),
muon = cms.VInputTag(),
)
),
),
mixObjects = cms.PSet(
mixHepMC = mix.mixObjects.mixHepMC.clone(),
mixVertices = mix.mixObjects.mixVertices.clone(),
mixTracks = mix.mixObjects.mixTracks.clone(),
),
)
from SimMuon.MCTruth.MuonAssociatorByHitsESProducer_NoSimHits_cfi import *

classByHitsTM = cms.EDProducer("MuonMCClassifier",
muons = cms.InputTag("muons"),
muonPreselection = cms.string("isTrackerMuon"), #
#muonPreselection = cms.string("muonID('TrackerMuonArbitrated')"), # You might want this
trackType = cms.string("segments"), # or 'inner','outer','global'
trackingParticles = cms.InputTag("mergedtruthNoSimHits"),
trackingParticles = cms.InputTag("trackingParticlesNoSimHits","MergedTrackTruth"),
associatorLabel = cms.string("muonAssociatorByHits_NoSimHits"),
decayRho = cms.double(200), # to classifiy differently decay muons included in ppMuX
decayAbsZ = cms.double(400), # and decay muons that could not be in ppMuX
Expand All @@ -34,7 +49,7 @@


muonClassificationByHits = cms.Sequence(
mix +
#mix +
trackingParticlesNoSimHits +
( classByHitsTM +
classByHitsTMLSAT +
Expand Down
46 changes: 36 additions & 10 deletions MuonAnalysis/MuonAssociators/test/MuonMCClassifier/debugMuonsClassByHits_cfg.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,18 +5,37 @@
process.load("FWCore.MessageService.MessageLogger_cfi")
process.MessageLogger.cerr.FwkReport.reportEvery = 100

process.load("Configuration.StandardSequences.Geometry_cff")
process.load('Configuration.StandardSequences.Services_cff')
process.load("Configuration.StandardSequences.GeometryDB_cff")
process.load("Configuration.StandardSequences.MagneticField_38T_cff")
process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff")
process.GlobalTag.globaltag = 'START38_V12::All'
process.GlobalTag.globaltag = 'PRE_ST62_V8::All'

process.source = cms.Source("PoolSource",
fileNames = cms.untracked.vstring(
'root://pcmssd12.cern.ch//data/gpetrucc/7TeV/mu11/store/mc/Fall10/QCD_Pt_600to800_TuneZ2_7TeV_pythia6/GEN-SIM-RECO/E7TeV_ProbDist_2010Data_BX156_START38_V12-v1/0096/22B09051-20EA-DF11-84C6-00151796C1E8.root',
#'root://pcmssd12.cern.ch//data/gpetrucc/7TeV/mu11/store/mc/Fall10/QCD_Pt-20_MuEnrichedPt-15_TuneZ2_7TeV-pythia6/GEN-SIM-RECO/E7TeV_ProbDist_2010Data_BX156_START38_V12-v1/0070/F044B493-E3E5-DF11-ACBD-001A92811728.root'
'/store/relval/CMSSW_6_2_0/RelValTTbar/GEN-SIM-RECO/PU_PRE_ST62_V8-v2/00000/E03F79C5-A7EC-E211-A92E-003048F00520.root',
),
secondaryFileNames = cms.untracked.vstring(
'/store/relval/CMSSW_6_2_0/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG/PU_PRE_ST62_V8-v2/00000/DE640769-94EC-E211-ACE9-003048F23D8C.root',
'/store/relval/CMSSW_6_2_0/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG/PU_PRE_ST62_V8-v2/00000/D21DB243-94EC-E211-AEED-003048D2BC36.root',
'/store/relval/CMSSW_6_2_0/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG/PU_PRE_ST62_V8-v2/00000/BC74017C-94EC-E211-9F4A-BCAEC532971D.root',
'/store/relval/CMSSW_6_2_0/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG/PU_PRE_ST62_V8-v2/00000/90485044-94EC-E211-8BC7-003048F0111A.root',
'/store/relval/CMSSW_6_2_0/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG/PU_PRE_ST62_V8-v2/00000/88F032AD-94EC-E211-BA22-0030486730E8.root',
'/store/relval/CMSSW_6_2_0/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG/PU_PRE_ST62_V8-v2/00000/78F4633D-94EC-E211-A1FE-003048F11D46.root',
'/store/relval/CMSSW_6_2_0/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG/PU_PRE_ST62_V8-v2/00000/68ED2290-94EC-E211-A08E-003048F16B8E.root',
'/store/relval/CMSSW_6_2_0/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG/PU_PRE_ST62_V8-v2/00000/2EDDA656-94EC-E211-A208-003048F16F46.root',
'/store/relval/CMSSW_6_2_0/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG/PU_PRE_ST62_V8-v2/00000/2C1CD01D-94EC-E211-AA58-001E673982E1.root',
'/store/relval/CMSSW_6_2_0/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG/PU_PRE_ST62_V8-v2/00000/18B87680-9EEC-E211-AD34-003048F1C410.root',
'/store/relval/CMSSW_6_2_0/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG/PU_PRE_ST62_V8-v2/00000/16040F35-94EC-E211-A460-003048F1DBB6.root',
'/store/relval/CMSSW_6_2_0/RelValTTbar/GEN-SIM-DIGI-RAW-HLTDEBUG/PU_PRE_ST62_V8-v2/00000/04E987D4-94EC-E211-9AA2-003048C9CC70.root',
),
)


## IF twoFileSolution == TRUE: test with TrackingParticles from input secondary files
## IF twoFileSolution == FALSE: test making them from GEN-SIM-RECO
twoFileSolution = False

process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) )
process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True) )

Expand All @@ -30,12 +49,19 @@

process.classByHits = process.classByHitsGlb.clone(muons = "selMuons", muonPreselection = "")

process.go = cms.Path(
process.selMuons +
process.mix +
process.trackingParticlesNoSimHits +
process.classByHits
)
if twoFileSolution:
process.classByHits.trackingParticles = cms.InputTag("mix","MergedTrackTruth")
process.go = cms.Path(
process.selMuons +
process.classByHits
)
else:
del process.source.secondaryFileNames
process.go = cms.Path(
process.selMuons +
process.trackingParticlesNoSimHits +
process.classByHits
)

process.MessageLogger.categories += [ 'MuonMCClassifier' ]
process.MessageLogger.cerr.MuonMCClassifier = cms.untracked.PSet(
Expand Down
31 changes: 15 additions & 16 deletions MuonAnalysis/MuonAssociators/test/MuonMCClassifier/patMuonsWithClassByHits_cfg.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,23 +5,18 @@
process.load("FWCore.MessageService.MessageLogger_cfi")
process.MessageLogger.cerr.FwkReport.reportEvery = 100

process.load("Configuration.StandardSequences.Geometry_cff")
process.load('Configuration.StandardSequences.Services_cff')
process.load("Configuration.StandardSequences.GeometryDB_cff")
process.load("Configuration.StandardSequences.MagneticField_38T_cff")
process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff")
process.load("Configuration.StandardSequences.Reconstruction_cff")

process.load("SimGeneral.MixingModule.mixNoPU_cfi")
process.load("SimGeneral.TrackingAnalysis.trackingParticlesNoSimHits_cfi") # On RECO
process.load("SimMuon.MCTruth.MuonAssociatorByHitsESProducer_NoSimHits_cfi") # On RECO

process.GlobalTag.globaltag = 'START3X_V26A::All'
process.GlobalTag.globaltag = 'PRE_ST62_V8::All'

from Configuration.EventContent.EventContent_cff import *
process.source = cms.Source("PoolSource",
fileNames = cms.untracked.vstring(
#'rfio:/castor/cern.ch/user/g/gpetrucc/900GeV/MC/Feb9Skims/MC_CollisionEvents_MuonSkim_1.root',
#'rfio:/castor/cern.ch/user/g/gpetrucc/900GeV/MC/Feb9Skims/MC_CollisionEvents_MuonSkim_2.root',
#'rfio:/castor/cern.ch/user/g/gpetrucc/900GeV/MC/Feb9Skims/MC_CollisionEvents_MuonSkim_3.root',
'root://pcmssd12.cern.ch//data/gpetrucc/Feb9Skims/MC_CollisionEvents_MuonSkim_1.root'
'/store/relval/CMSSW_6_2_0/RelValTTbar/GEN-SIM-RECO/PU_PRE_ST62_V8-v2/00000/E03F79C5-A7EC-E211-A92E-003048F00520.root',
),
inputCommands = RECOSIMEventContent.outputCommands, # keep only RECO out of RAW+RECO, for tests
dropDescendantsOfDroppedBranches = cms.untracked.bool(False), # keep only RECO out of RAW+RECO, for tests
Expand All @@ -30,10 +25,6 @@
process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) )
process.options = cms.untracked.PSet( wantSummary = cms.untracked.bool(True) )

process.load('L1TriggerConfig.L1GtConfigProducers.L1GtTriggerMaskTechTrigConfig_cff')
from HLTrigger.HLTfilters.hltLevel1GTSeed_cfi import hltLevel1GTSeed
hltLevel1GTSeed.L1TechTriggerSeeding = cms.bool(True)
process.bscFilter = hltLevel1GTSeed.clone(L1SeedsLogicalExpression = cms.string('(40 OR 41) AND NOT (36 OR 37 OR 38 OR 39)'))
process.oneGoodVertexFilter = cms.EDFilter("VertexSelector",
src = cms.InputTag("offlinePrimaryVertices"),
cut = cms.string("!isFake && ndof >= 4 && abs(z) <= 15 && position.Rho <= 2"),
Expand All @@ -46,11 +37,11 @@
thresh = cms.untracked.double(0.25)
)
process.countCollisionEvents = cms.EDProducer("EventCountProducer")
process.preFilter = cms.Sequence(process.bscFilter * process.oneGoodVertexFilter * process.noScraping * process.countCollisionEvents )
process.preFilter = cms.Sequence(process.oneGoodVertexFilter * process.noScraping * process.countCollisionEvents )

process.mergedTruth = cms.EDProducer("GenPlusSimParticleProducer",
src = cms.InputTag("g4SimHits"), # use "famosSimHits" for FAMOS
setStatus = cms.int32(5), # set status = 8 for GEANT GPs
setStatus = cms.int32(5), # set status = 5 for GEANT GPs
filter = cms.vstring("pt > 0.0"), # just for testing (optional)
genParticles = cms.InputTag("genParticles") # original genParticle list
)
Expand All @@ -76,6 +67,9 @@
embedTrack = True,
embedCombinedMuon = True,
embedStandAloneMuon = True,
embedPFCandidate = False,
embedCaloMETMuonCorrs = cms.bool(False),
embedTcMETMuonCorrs = cms.bool(False),
# then switch off some features we don't need
#addTeVRefits = False, ## <<--- this doesn't work. PAT bug ??
embedPickyMuon = False,
Expand All @@ -85,6 +79,11 @@
addGenMatch = False,
embedGenMatch = False,
)
# Reset all these; the default in muonProducer_cfi is not empty, but wrong
process.patMuons.userData.userInts.src = []
process.patMuons.userData.userFloats.src = []
process.patMuons.userData.userCands.src = []
process.patMuons.userData.userClasses.src = []

process.load("MuonAnalysis.MuonAssociators.muonClassificationByHits_cfi")

Expand Down

0 comments on commit a5c8f17

Please sign in to comment.