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

Adding digi-based shower tagging to reco::Muon #26369

Merged
merged 9 commits into from Apr 18, 2019
10 changes: 8 additions & 2 deletions DataFormats/MuonReco/interface/Muon.h
Expand Up @@ -185,7 +185,7 @@ namespace reco {

enum ArbitrationType { NoArbitration, SegmentArbitration, SegmentAndTrackArbitration, SegmentAndTrackArbitrationCleaned,
RPCHitAndTrackArbitration, GEMSegmentAndTrackArbitration, ME0SegmentAndTrackArbitration };

///
/// ====================== STANDARD SELECTORS ===========================
///
Expand Down Expand Up @@ -260,7 +260,13 @@ namespace reco {
unsigned int stationGapMaskDistance( float distanceCut = 10. ) const;
/// same as above for given number of sigmas
unsigned int stationGapMaskPull( float sigmaCut = 3. ) const;

/// # of digis in a given station layer
int nDigisInStation( int station, int muonSubdetId) const;
/// tag a shower in a given station layer
bool hasShowerInStation( int station, int muonSubdetId, int nDtDigisCut = 20, int nCscDigisCut = 36 ) const;
/// count the number of showers along a muon track
int numberOfShowers( int nDtDigisCut = 20, int nCscDigisCut = 36 ) const;

/// muon type - type of the algorithm that reconstructed this muon
/// multiple algorithms can reconstruct the same muon
static const unsigned int GlobalMuon = 1<<1;
Expand Down
24 changes: 13 additions & 11 deletions DataFormats/MuonReco/interface/MuonChamberMatch.h
Expand Up @@ -14,17 +14,19 @@ namespace reco {
std::vector<reco::MuonSegmentMatch> me0Matches; // segments matching propagated track trajectory
std::vector<reco::MuonSegmentMatch> truthMatches; // SimHit projection matching propagated track trajectory
std::vector<reco::MuonRPCHitMatch> rpcMatches; // rpc hits matching propagated track trajectory
float edgeX; // distance to closest edge in X (negative - inside, positive - outside)
float edgeY; // distance to closest edge in Y (negative - inside, positive - outside)
float x; // X position of the track
float y; // Y position of the track
float xErr; // propagation uncertainty in X
float yErr; // propagation uncertainty in Y
float dXdZ; // dX/dZ of the track
float dYdZ; // dY/dZ of the track
float dXdZErr; // propagation uncertainty in dX/dZ
float dYdZErr; // propagation uncertainty in dY/dZ
DetId id; // chamber ID
float edgeX; // distance to closest edge in X (negative - inside, positive - outside)
float edgeY; // distance to closest edge in Y (negative - inside, positive - outside)
float x; // X position of the track
float y; // Y position of the track
float xErr; // propagation uncertainty in X
float yErr; // propagation uncertainty in Y
float dXdZ; // dX/dZ of the track
float dYdZ; // dY/dZ of the track
float dXdZErr; // propagation uncertainty in dX/dZ
float dYdZErr; // propagation uncertainty in dY/dZ
DetId id; // chamber ID

int nDigisInRange; // # of DT/CSC digis in the chamber close-by to the propagated track

int detector() const { return id.subdetId(); }
int station() const;
Expand Down
72 changes: 72 additions & 0 deletions DataFormats/MuonReco/src/Muon.cc
Expand Up @@ -334,6 +334,78 @@ unsigned int Muon::stationGapMaskPull( float sigmaCut ) const
return totMask;
}

int Muon::nDigisInStation( int station, int muonSubdetId ) const
{
int nDigis(0);
std::map<int, int> me11DigisPerCh;

if ( muonSubdetId != MuonSubdetId::CSC &&
muonSubdetId != MuonSubdetId::DT )
return 0;

for ( auto & match : muMatches_ )
{
if ( match.detector() != muonSubdetId ||
match.station() != station )
continue;

int nDigisInCh = match.nDigisInRange;

if( muonSubdetId == MuonSubdetId::CSC && station == 1 )
{
CSCDetId id(match.id.rawId());

int chamber = id.chamber();
int ring = id.ring();

if ( ring == 1 || ring == 4 ) // merge ME1/1a and ME1/1b digis
{
if( me11DigisPerCh.find(chamber) == me11DigisPerCh.end() )
me11DigisPerCh[chamber] = 0;

me11DigisPerCh[chamber] += nDigisInCh;

continue;
}
}

if( nDigisInCh > nDigis )
nDigis = nDigisInCh;
}

for ( const auto & me11DigisInCh : me11DigisPerCh )
{
int nMe11DigisInCh = me11DigisInCh.second;
if ( nMe11DigisInCh > nDigis )
nDigis = nMe11DigisInCh;
}

return nDigis;
}

bool Muon::hasShowerInStation( int station, int muonSubdetId, int nDtDigisCut, int nCscDigisCut ) const
{
bool hasShower = muonSubdetId == MuonSubdetId::DT ?
nDigisInStation(station,muonSubdetId) >= nDtDigisCut :
Copy link
Contributor

Choose a reason for hiding this comment

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

if (muonSubdetId != MuonSubdetId::DT || muonSubdetId != MuonSubdetId::CSC) return false;
auto cut = muonSubdetId == MuonSubdetId::DT ? nDtDigisCut : nCscDigisCut;
return nDigisInStation(station,muonSubdetId) >= cut;

is perhaps shorter and more explicit about DT vs CSC.

Copy link
Author

Choose a reason for hiding this comment

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

Right.

Copy link
Author

Choose a reason for hiding this comment

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

But the OR should be an AND.

nDigisInStation(station,muonSubdetId) >= nCscDigisCut;

return hasShower;
}

int Muon::numberOfShowers( int nDtDigisCut, int nCscDigisCut ) const
{
int nShowers = 0;
for ( int iCh = 1; iCh < 5; ++iCh )
Copy link
Contributor

Choose a reason for hiding this comment

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

to reduce possible confusion when reading this, perhaps iSt instead of iCh because it's supposed to be a shorthand for a station

Copy link
Author

Choose a reason for hiding this comment

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

Changed in station (more consistent with the member functions above).

{
if ( hasShowerInStation(iCh,MuonSubdetId::DT,nDtDigisCut,nCscDigisCut) )
nShowers++;
if ( hasShowerInStation(iCh,MuonSubdetId::CSC,nDtDigisCut,nCscDigisCut) )
nShowers++;
}

return nShowers;
}

int Muon::numberOfSegments( int station, int muonSubdetId, ArbitrationType type ) const
{
int segments(0);
Expand Down
3 changes: 2 additions & 1 deletion DataFormats/MuonReco/src/classes_def.xml
Expand Up @@ -57,7 +57,8 @@ initial version number of a class which has never been stored before.
<class name="reco::HcalMuonRecHit" ClassVersion="3">
<version ClassVersion="3" checksum="2031556424"/>
</class>
<class name="reco::MuonChamberMatch" ClassVersion="12">
<class name="reco::MuonChamberMatch" ClassVersion="13">
<version ClassVersion="13" checksum="2674954213"/>
<version ClassVersion="12" checksum="2575855272"/>
<version ClassVersion="11" checksum="541727491"/>
<version ClassVersion="10" checksum="4050071853"/>
Expand Down
75 changes: 75 additions & 0 deletions RecoMuon/MuonIdentification/interface/MuonShowerDigiFiller.h
@@ -0,0 +1,75 @@
#ifndef MuonIdentification_MuonShowerDigiFiller_h
#define MuonIdentification_MuonShowerDigiFiller_h

// -*- C++ -*-
//
// Package: MuonShowerDigiFiller
// Class: MuonShowerDigiFiller
//
/**\class MuonShowerDigiFiller MuonShowerDigiFiller.h RecoMuon/MuonIdentification/interface/MuonShowerDigiFiller.h

Description: Class filling shower information using DT and CSC digis

Implementation:
<Notes on implementation>
*/
//
// Original Author: Carlo Battilana, INFN BO
// Created: Sat Mar 23 14:36:22 CET 2019
//
//

// system include files

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"

#include "DataFormats/DTDigi/interface/DTDigiCollection.h"
#include "DataFormats/CSCDigi/interface/CSCStripDigiCollection.h"

#include "Geometry/DTGeometry/interface/DTGeometry.h"
#include "Geometry/CSCGeometry/interface/CSCGeometry.h"

#include "DataFormats/MuonReco/interface/MuonChamberMatch.h"
#include "TrackingTools/TrackAssociator/interface/TAMuonChamberMatch.h"

//
// class decleration
//

class MuonShowerDigiFiller
{

public:

MuonShowerDigiFiller(const edm::ParameterSet&, edm::ConsumesCollector&& iC);

void getES( const edm::EventSetup& iSetup );
void getDigis( edm::Event& iEvent );

void fill( reco::MuonChamberMatch & muChMatch ) const;
void fillDefault( reco::MuonChamberMatch & muChMatch ) const;

private:

double m_digiMaxDistanceX;

edm::EDGetTokenT<DTDigiCollection> m_dtDigisToken;
edm::EDGetTokenT<CSCStripDigiCollection> m_cscDigisToken;

edm::ESHandle<DTGeometry> m_dtGeometry;
edm::ESHandle<CSCGeometry> m_cscGeometry;

edm::Handle<DTDigiCollection> m_dtDigis;
edm::Handle<CSCStripDigiCollection> m_cscDigis;

};

#endif
25 changes: 24 additions & 1 deletion RecoMuon/MuonIdentification/plugins/MuonIdProducer.cc
Expand Up @@ -59,6 +59,7 @@ MuonIdProducer::MuonIdProducer(const edm::ParameterSet& iConfig)
storeCrossedHcalRecHits_ = iConfig.getParameter<bool>("storeCrossedHcalRecHits");
fillMatching_ = iConfig.getParameter<bool>("fillMatching");
fillIsolation_ = iConfig.getParameter<bool>("fillIsolation");
fillShowerDigis_ = iConfig.getParameter<bool>("fillShowerDigis");
writeIsoDeposits_ = iConfig.getParameter<bool>("writeIsoDeposits");
fillGlobalTrackQuality_ = iConfig.getParameter<bool>("fillGlobalTrackQuality");
fillGlobalTrackRefits_ = iConfig.getParameter<bool>("fillGlobalTrackRefits");
Expand All @@ -80,7 +81,13 @@ MuonIdProducer::MuonIdProducer(const edm::ParameterSet& iConfig)
// Load parameters for the TimingFiller
edm::ParameterSet timingParameters = iConfig.getParameter<edm::ParameterSet>("TimingFillerParameters");
theTimingFiller_ = std::make_unique<MuonTimingFiller>(timingParameters,consumesCollector());


// Load parameters for the ShowerDigiFiller
if (fillShowerDigis_ && fillMatching_)
{
edm::ParameterSet showerDigiParameters = iConfig.getParameter<edm::ParameterSet>("ShowerDigiFillerParameters");
theShowerDigiFiller_ = std::make_unique<MuonShowerDigiFiller>(showerDigiParameters,consumesCollector());
}

if (fillCaloCompatibility_){
// Load MuonCaloCompatibility parameters
Expand Down Expand Up @@ -426,6 +433,8 @@ void MuonIdProducer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetu

meshAlgo_->setCSCGeometry(geomHandle.product());

if (fillShowerDigis_ && fillMatching_)
theShowerDigiFiller_->getES(iSetup);
}

bool validateGlobalMuonPair( const reco::MuonTrackLinks& goodMuon,
Expand All @@ -450,6 +459,9 @@ void MuonIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)

init(iEvent, iSetup);

if (fillShowerDigis_ && fillMatching_)
theShowerDigiFiller_->getDigis(iEvent);

// loop over input collections

// muons first - no cleaning, take as is.
Expand Down Expand Up @@ -876,6 +888,14 @@ void MuonIdProducer::fillMuonId(edm::Event& iEvent, const edm::EventSetup& iSetu
matchedChamber.edgeY = chamber.localDistanceY;

matchedChamber.id = chamber.id;

if (fillShowerDigis_) {
theShowerDigiFiller_->fill(matchedChamber);
}
else {
theShowerDigiFiller_->fillDefault(matchedChamber);
}

if ( ! chamber.segments.empty() ) ++nubmerOfMatchesAccordingToTrackAssociator;

// fill segments
Expand Down Expand Up @@ -958,6 +978,8 @@ void MuonIdProducer::fillMuonId(edm::Event& iEvent, const edm::EventSetup& iSetu
matchedChamber.edgeX = chamber.localDistanceX;
matchedChamber.edgeY = chamber.localDistanceY;

theShowerDigiFiller_->fillDefault(matchedChamber);

matchedChamber.id = chamber.id;

for ( const auto& rpcRecHit : *rpcHitHandle_ )
Expand Down Expand Up @@ -1324,6 +1346,7 @@ void MuonIdProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptio

desc.add<bool>("arbitrateTrackerMuons",false);
desc.add<bool>("storeCrossedHcalRecHits",false);
desc.add<bool>("fillShowerDigis",false);

edm::ParameterSetDescription descTrkAsoPar;
descTrkAsoPar.add<edm::InputTag>("GEMSegmentCollectionLabel",edm::InputTag("gemSegments"));
Expand Down
6 changes: 5 additions & 1 deletion RecoMuon/MuonIdentification/plugins/MuonIdProducer.h
Expand Up @@ -47,8 +47,9 @@
#include "RecoMuon/MuonIdentification/interface/MuonTimingFiller.h"
#include "RecoMuon/MuonIdentification/interface/MuonCaloCompatibility.h"
#include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractor.h"
// RPC-Muon stuffs
#include "RecoMuon/MuonIdentification/interface/MuonShowerDigiFiller.h"

// RPC-Muon stuffs
#include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
#include "DataFormats/RPCRecHit/interface/RPCRecHit.h"
#include "DataFormats/MuonReco/interface/MuonRPCHitMatch.h"
Expand Down Expand Up @@ -177,6 +178,8 @@ class MuonIdProducer : public edm::stream::EDProducer<> {

std::unique_ptr<MuonTimingFiller> theTimingFiller_;

std::unique_ptr<MuonShowerDigiFiller> theShowerDigiFiller_;

// selections
double minPt_;
double minP_;
Expand All @@ -196,6 +199,7 @@ class MuonIdProducer : public edm::stream::EDProducer<> {
bool fillEnergy_;
bool storeCrossedHcalRecHits_;
bool fillMatching_;
bool fillShowerDigis_;
bool fillIsolation_;
bool writeIsoDeposits_;
double ptThresholdToFillCandidateP4WithGlobalFit_;
Expand Down
11 changes: 11 additions & 0 deletions RecoMuon/MuonIdentification/python/MuonShowerDigiFiller_cfi.py
@@ -0,0 +1,11 @@
import FWCore.ParameterSet.Config as cms

MuonShowerDigiFillerBlock = cms.PSet(
ShowerDigiFillerParameters = cms.PSet(
digiMaxDistanceX = cms.double(25.0),
dtDigiCollectionLabel = cms.InputTag("muonDTDigis"),
cscDigiCollectionLabel = cms.InputTag("muonCSCDigis","MuonCSCStripDigi")
Copy link
Contributor

Choose a reason for hiding this comment

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

@battibass @slava77 it looks like this line is creating issues in wf 1102.0 step4 RECOFROMRECO:

   [7] Prefetching for module TrackProducer/'mergedDuplicateTracks'
   [8] Prefetching for module DuplicateTrackMerger/'duplicateTrackCandidates'
   [9] Prefetching for module TrackCollectionMerger/'preDuplicateMergingGeneralTracks'
   [10] Prefetching for module TrackProducer/'muonSeededTracksInOut'
   [11] Prefetching for module CkfTrackCandidateMaker/'muonSeededTrackCandidatesInOut'
   [12] Prefetching for module MuonReSeeder/'muonSeededSeedsInOut'
   [13] Calling method for module MuonIdProducer/'earlyMuons'
Exception Message:
Principal::getByToken: Found zero products matching all criteria
Looking for type: MuonDigiCollection&lt;CSCDetId,CSCStripDigi&gt;
Looking for module label: muonCSCDigis
Looking for productInstanceName: MuonCSCStripDigi

   Additional Info:
      [a] If you wish to continue processing events after a ProductNotFound exception,
add "SkipEvent = cms.untracked.vstring('ProductNotFound')" to the "options" PSet in the configuration.

----- End Fatal Exception -------------------------------------------------

This clearly needs to be fixed, my question is: if this is not affecting any other standard workflow as it seems, and it is not preventing the basic validation of pre4, can we consider to postpone its fix after pre4 is built? Or do you see further possible issues?

Switching this module off in the configuration allows also the step4 of 1102 to run without problems.

Copy link
Author

Choose a reason for hiding this comment

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

@fabiocos , @slava77 can surely reply more accurately than I can.
In principle this PR implies that digis are available (hence implies that unpacking was done starting from RAW). For any workflow for which this is the case there should be no problem.
The question is, are there RelVal workflows where this is not the case?
I'm happy to work on a patch asap, but I would need inputs to understand better how I should be testing it, as my check with runTheMatrix.py -l limited -i all clearly failed to spot this issue, as well as the one with TauAnalysis/MCEmbeddingTools/ that was instead spot during integration ...

Copy link
Contributor

Choose a reason for hiding this comment

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

@battibass indeed this is a bit special case, and only wf 1102.0 implementing it shows this issue. I see that @slava77 has already mentioned it in #26500 .

)
)


4 changes: 4 additions & 0 deletions RecoMuon/MuonIdentification/python/muons1stStep_cfi.py
Expand Up @@ -4,6 +4,7 @@
from RecoMuon.MuonIdentification.isolation_cff import *
from RecoMuon.MuonIdentification.caloCompatibility_cff import *
from RecoMuon.MuonIdentification.MuonTimingFiller_cfi import *
from RecoMuon.MuonIdentification.MuonShowerDigiFiller_cfi import *
from RecoMuon.MuonIdentification.TrackerKinkFinder_cfi import *
from TrackingTools.TrackAssociator.default_cfi import *
muons1stStep = cms.EDProducer("MuonIdProducer",
Expand All @@ -15,6 +16,8 @@
MIdIsoExtractorPSetBlock,
# MuonTiming
TimingFillerBlock,
# MuonShowerDigi
MuonShowerDigiFillerBlock,
# Kink finder
TrackerKinkFinderParametersBlock,

Expand Down Expand Up @@ -55,6 +58,7 @@
writeIsoDeposits = cms.bool(True),
minNumberOfMatches = cms.int32(1),
fillMatching = cms.bool(True),
fillShowerDigis = cms.bool(True),

# global fit for candidate p4 requirements
ptThresholdToFillCandidateP4WithGlobalFit = cms.double(200.0),
Expand Down