diff --git a/Validation/MuonHits/BuildFile.xml b/Validation/MuonHits/BuildFile.xml index 8896987508d9c..9dd05125a275c 100644 --- a/Validation/MuonHits/BuildFile.xml +++ b/Validation/MuonHits/BuildFile.xml @@ -10,7 +10,12 @@ + + + - + + + diff --git a/Validation/MuonHits/interface/CSCSimHitMatcher.h b/Validation/MuonHits/interface/CSCSimHitMatcher.h new file mode 100644 index 0000000000000..f8d9c59458a61 --- /dev/null +++ b/Validation/MuonHits/interface/CSCSimHitMatcher.h @@ -0,0 +1,79 @@ +#ifndef Validation_MuonHits_CSCSimHitMatcher_h +#define Validation_MuonHits_CSCSimHitMatcher_h + +/**\class CSCSimHitMatcher + + Description: Matching of CSC SimHit to SimTrack + + Author: Sven Dildick (TAMU), Tao Huang (TAMU) +*/ + +#include "Geometry/CSCGeometry/interface/CSCGeometry.h" +#include "Validation/MuonHits/interface/MuonSimHitMatcher.h" + +class CSCSimHitMatcher : public MuonSimHitMatcher { + public: + // constructor + CSCSimHitMatcher(const edm::ParameterSet& iPS, edm::ConsumesCollector&& iC); + + // destructor + ~CSCSimHitMatcher() {} + + // initialize the event + void init(const edm::Event& e, const edm::EventSetup& eventSetup); + + // do the matching + void match(const SimTrack& t, const SimVertex& v); + + // partitions' detIds with SimHits + std::set detIds(int type = MuonHitHelper::CSC_ALL) const; + + // chamber detIds with SimHits + std::set chamberIds(int type = MuonHitHelper::CSC_ALL) const; + + // CSC station detIds with SimHits + std::set chamberIdsStation(int station) const; + + // was there a hit in a particular CSC station? + bool hitStation(int, int) const; + + // number of stations with hits in at least X layers + int nStations(int nl = 4) const; + + // #layers with hits + int nLayersWithHitsInChamber(unsigned int) const; + + // How many CSC chambers with minimum number of layer with simhits did this + // simtrack get? + int nCoincidenceChambers(int min_n_layers = 4) const; + + // calculated the fitted position in a given layer for CSC simhits in a + // chamber + GlobalPoint simHitPositionKeyLayer(unsigned int chamberid) const; + + // local bending in a CSC chamber + float LocalBendingInChamber(unsigned int detid) const; + + // calculate average strip number for a provided collection of simhits + float simHitsMeanStrip(const edm::PSimHitContainer& sim_hits) const; + + // calculate average wg number for a provided collection of simhits (for CSC) + float simHitsMeanWG(const edm::PSimHitContainer& sim_hits) const; + + void chamberIdsToString(const std::set& set) const; + + // calculate the average position at the second station + GlobalPoint simHitsMeanPositionStation(int n) const; + + std::set hitStripsInDetId(unsigned int, int margin_n_strips = 0) const; + std::set hitWiregroupsInDetId(unsigned int, int margin_n_wg = 0) const; + + void camberIdsToString(const std::set&) const; + + private: + void matchSimHitsToSimTrack(); + + edm::ESHandle csc_geom_; +}; + +#endif diff --git a/Validation/MuonHits/interface/DTSimHitMatcher.h b/Validation/MuonHits/interface/DTSimHitMatcher.h new file mode 100644 index 0000000000000..98895b6b0a515 --- /dev/null +++ b/Validation/MuonHits/interface/DTSimHitMatcher.h @@ -0,0 +1,82 @@ +#ifndef Validation_MuonHits_DTSimHitMatcher_h +#define Validation_MuonHits_DTSimHitMatcher_h + +/**\class DTSimHitMatcher + + Description: Matching of DT SimHit to SimTrack + + Author: Sven Dildick (TAMU), Tao Huang (TAMU) +*/ + +#include "Geometry/DTGeometry/interface/DTGeometry.h" +#include "Validation/MuonHits/interface/MuonSimHitMatcher.h" + +class DTSimHitMatcher : public MuonSimHitMatcher { + public: + // constructor + DTSimHitMatcher(const edm::ParameterSet& iPS, edm::ConsumesCollector&& iC); + + // destructor + ~DTSimHitMatcher() {} + + // initialize the event + void init(const edm::Event& e, const edm::EventSetup& eventSetup); + + // do the matching + void match(const SimTrack& t, const SimVertex& v); + + // partitions' detIds with SimHits + std::set detIds(int type = MuonHitHelper::DT_ALL) const; + + // chamber detIds with SimHits + std::set chamberIds(int type = MuonHitHelper::DT_ALL) const; + + // DT station detIds with SimHits + std::set chamberIdsStation(int station) const; + + // DT layer detIds with SimHits + std::set layerIds() const; + + // DT super layer detIds with SimHits + std::set superlayerIds() const; + + // was there a hit in a particular DT/CSC station? + bool hitStation(int, int, int) const; + + // number of stations with hits in at least X layers + int nStations(int nsl = 1, int nl = 3) const; + + // access to DT hits + int nCellsWithHitsInLayer(unsigned int) const; + int nLayersWithHitsInSuperLayer(unsigned int) const; + int nSuperLayersWithHitsInChamber(unsigned int) const; + int nLayersWithHitsInChamber(unsigned int) const; + const edm::PSimHitContainer& hitsInLayer(unsigned int) const; + const edm::PSimHitContainer& hitsInSuperLayer(unsigned int) const; + const edm::PSimHitContainer& hitsInChamber(unsigned int) const; + + // calculate average wg number for a provided collection of simhits + float simHitsMeanWire(const edm::PSimHitContainer& sim_hits) const; + + // calculate the average position at the second station + GlobalPoint simHitsMeanPositionStation(int n) const; + + std::set hitWiresInDTLayerId( + unsigned int, int margin_n_wires = 0) const; // DT + std::set hitWiresInDTSuperLayerId( + unsigned int, int margin_n_wires = 0) const; // DT + std::set hitWiresInDTChamberId( + unsigned int, int margin_n_wires = 0) const; // DT + + void dtChamberIdsToString(const std::set&) const; + + private: + void matchSimHitsToSimTrack(); + + std::map layer_to_hits_; + std::map superlayer_to_hits_; + + edm::ESHandle dt_geom_; +}; + +#endif diff --git a/Validation/MuonHits/interface/GEMSimHitMatcher.h b/Validation/MuonHits/interface/GEMSimHitMatcher.h new file mode 100644 index 0000000000000..eda6fe813bbbd --- /dev/null +++ b/Validation/MuonHits/interface/GEMSimHitMatcher.h @@ -0,0 +1,90 @@ +#ifndef Validation_MuonHits_GEMSimHitMatcher_h +#define Validation_MuonHits_GEMSimHitMatcher_h + +/**\class GEMSimHitMatcher + + Description: Matching of GEM SimHit to SimTrack + + Author: Sven Dildick (TAMU), Tao Huang (TAMU) +*/ + +#include "Geometry/GEMGeometry/interface/GEMGeometry.h" +#include "Validation/MuonHits/interface/MuonSimHitMatcher.h" + +class GEMSimHitMatcher : public MuonSimHitMatcher { + public: + // constructor + GEMSimHitMatcher(const edm::ParameterSet& iPS, edm::ConsumesCollector&& iC); + + // destructor + ~GEMSimHitMatcher() {} + + // initialize the event + void init(const edm::Event& e, const edm::EventSetup& eventSetup); + + // do the matching + void match(const SimTrack& t, const SimVertex& v); + + // partitions' detIds with SimHits + std::set detIds(int gem_type = MuonHitHelper::GEM_ALL) const; + + // chamber detIds with SimHits + std::set chamberIds( + int gem_type = MuonHitHelper::GEM_ALL) const; + + // GEM detid's with hits in 2 layers of coincidence pads + // those are layer==1 only detid's + std::set detIdsCoincidences() const; + + // GEM superchamber detIds with SimHits + std::set superChamberIds() const; + + // GEM superchamber detIds with SimHits 2 layers of coincidence pads + std::set superChamberIdsCoincidences() const; + + // simhits from a particular superchamber + const edm::PSimHitContainer& hitsInSuperChamber(unsigned int) const; + + // was there a hit in a particular station? + bool hitStation(int, int) const; + + // number of stations with hits in at least X layers + int nStations(int nl = 2) const; + + // #layers with hits + int nLayersWithHitsInSuperChamber(unsigned int) const; + + // How many pads with simhits in GEM did this simtrack get? + int nPadsWithHits() const; + + // How many coincidence pads with simhits in GEM did this simtrack get? + int nCoincidencePadsWithHits() const; + + // transverse position in GEM + float simHitsGEMCentralPosition(const edm::PSimHitContainer& sim_hits) const; + + // calculate average strip number for a provided collection of simhits + float simHitsMeanStrip(const edm::PSimHitContainer& sim_hits) const; + + std::set hitStripsInDetId(unsigned int, int margin_n_strips = 0) const; + std::set hitPadsInDetId(unsigned int) const; + std::set hitCoPadsInDetId(unsigned int) const; + + // what unique partitions numbers were hit by this simtrack? + std::set hitPartitions() const; + + private: + void matchSimHitsToSimTrack(); + + edm::ESHandle gem_geom_; + + std::map superchamber_to_hits_; + + // detids with hits in pads + std::map > detids_to_pads_; + + // detids with hits in 2-layer pad coincidences + std::map > detids_to_copads_; +}; + +#endif diff --git a/Validation/MuonHits/interface/ME0SimHitMatcher.h b/Validation/MuonHits/interface/ME0SimHitMatcher.h new file mode 100644 index 0000000000000..38b433a65b1e9 --- /dev/null +++ b/Validation/MuonHits/interface/ME0SimHitMatcher.h @@ -0,0 +1,77 @@ +#ifndef Validation_MuonHits_ME0SimHitMatcher_h +#define Validation_MuonHits_ME0SimHitMatcher_h + +/**\class ME0SimHitMatcher + + Description: Matching of ME0 SimHit to SimTrack + + Author: Sven Dildick (TAMU), Tao Huang (TAMU) +*/ + +#include "Geometry/GEMGeometry/interface/ME0Geometry.h" +#include "Validation/MuonHits/interface/MuonSimHitMatcher.h" + +class ME0SimHitMatcher : public MuonSimHitMatcher { + public: + // constructor + ME0SimHitMatcher(const edm::ParameterSet& iPS, edm::ConsumesCollector&& iC); + + // destructor + ~ME0SimHitMatcher() {} + + // initialize the event + void init(const edm::Event& e, const edm::EventSetup& eventSetup); + + // do the matching + void match(const SimTrack& t, const SimVertex& v); + + // partitions' detIds with SimHits + std::set detIds() const; + + // chamber detIds with SimHits + std::set chamberIds() const; + + // ME0 superchamber detIds with SimHits + std::set superChamberIds() const; + + // simhits from a particular partition, chamber + const edm::PSimHitContainer& hitsInSuperChamber(unsigned int) const; + + // #layers with hits + int nLayersWithHitsInSuperChamber(unsigned int) const; + + // ME0 superchamber detIds with SimHits >=4 layers of coincidence pads + std::set superChamberIdsCoincidences( + int min_n_layers = 4) const; + + // How many ME0 chambers with minimum number of layer with simhits did this + // simtrack get? + int nCoincidenceChambers(int min_n_layers = 4) const; + + // calculate average strip for a provided collection of simhits + float simHitsMeanStrip(const edm::PSimHitContainer& sim_hits) const; + + std::set hitStripsInDetId(unsigned int, int margin_n_strips = 0) const; + std::set hitPadsInDetId(unsigned int) const; + + // what unique partitions numbers were hit by this simtrack? + std::set hitPartitions() const; + + // How many pads with simhits in ME0 did this simtrack get? + int nPadsWithHits() const; + + private: + void matchSimHitsToSimTrack(); + + edm::ESHandle me0_geom_; + + // detids with hits in pads + std::map > detids_to_pads_; + + // detids with hits in 2-layer pad coincidences + std::map > detids_to_copads_; + + std::map superChamber_to_hits_; +}; + +#endif diff --git a/Validation/MuonHits/interface/MuonHitHelper.h b/Validation/MuonHits/interface/MuonHitHelper.h new file mode 100644 index 0000000000000..fd97dfd9645c8 --- /dev/null +++ b/Validation/MuonHits/interface/MuonHitHelper.h @@ -0,0 +1,109 @@ +#ifndef Validation_MuonHits_MuonHitHelper_h +#define Validation_MuonHits_MuonHitHelper_h + +#include "DataFormats/MuonDetId/interface/CSCDetId.h" +#include "DataFormats/MuonDetId/interface/CSCTriggerNumbering.h" +#include "DataFormats/MuonDetId/interface/DTWireId.h" +#include "DataFormats/MuonDetId/interface/GEMDetId.h" +#include "DataFormats/MuonDetId/interface/ME0DetId.h" +#include "DataFormats/MuonDetId/interface/RPCDetId.h" + +class MuonHitHelper { + public: + /// CSC chamber types, according to CSCDetId::iChamberType() + enum CSCType { + CSC_ALL = 0, + CSC_ME11, + CSC_ME1a, + CSC_ME1b, + CSC_ME12, + CSC_ME13, + CSC_ME21, + CSC_ME22, + CSC_ME31, + CSC_ME32, + CSC_ME41, + CSC_ME42 + }; + + /// GEM chamber types + enum GEMType { GEM_ALL = 0, GEM_ME11, GEM_ME21 }; + + /// RPC endcap chamber types + enum RPCType { + RPC_ALL = 0, + RPC_ME12, + RPC_ME13, + RPC_ME22, + RPC_ME23, + RPC_ME31, + RPC_ME32, + RPC_ME33, + RPC_ME41, + RPC_ME42, + RPC_ME43, + RPC_MB01, + RPC_MB02, + RPC_MB03, + RPC_MB04, + RPC_MB11p, + RPC_MB12p, + RPC_MB13p, + RPC_MB14p, + RPC_MB21p, + RPC_MB22p, + RPC_MB23p, + RPC_MB24p, + RPC_MB11n, + RPC_MB12n, + RPC_MB13n, + RPC_MB14n, + RPC_MB21n, + RPC_MB22n, + RPC_MB23n, + RPC_MB24n + }; + + /// DT chamber types + enum DTType { + DT_ALL = 0, + DT_MB01, + DT_MB02, + DT_MB03, + DT_MB04, + DT_MB11p, + DT_MB12p, + DT_MB13p, + DT_MB14p, + DT_MB21p, + DT_MB22p, + DT_MB23p, + DT_MB24p, + DT_MB11n, + DT_MB12n, + DT_MB13n, + DT_MB14n, + DT_MB21n, + DT_MB22n, + DT_MB23n, + DT_MB24n + }; + + /// check detid type + static bool isDT(unsigned int detId); + static bool isGEM(unsigned int detId); + static bool isCSC(unsigned int detId); + static bool isRPC(unsigned int detId); + static bool isME0(unsigned int detId); + + // return MuonType for a particular DetId + static int toGEMType(int st, int ri); + static int toRPCType(int re, int st, int ri); + static int toDTType(int wh, int st); + static int toCSCType(int st, int ri); + + // get chamber number + static int chamber(const DetId& id); +}; + +#endif diff --git a/Validation/MuonHits/interface/MuonSimHitMatcher.h b/Validation/MuonHits/interface/MuonSimHitMatcher.h new file mode 100644 index 0000000000000..eb99675850064 --- /dev/null +++ b/Validation/MuonHits/interface/MuonSimHitMatcher.h @@ -0,0 +1,109 @@ +#ifndef Validation_MuonHits_MuonSimHitMatcher_h +#define Validation_MuonHits_MuonSimHitMatcher_h + +/**\class MuonSimHitMatcher + + Description: Matching of muon SimHit to SimTrack + + Author: Sven Dildick (TAMU), Tao Huang (TAMU) +*/ + +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include "DataFormats/Math/interface/deltaPhi.h" +#include "FWCore/Framework/interface/ConsumesCollector.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "Geometry/CommonDetUnit/interface/TrackingGeometry.h" +#include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h" +#include "Geometry/Records/interface/MuonGeometryRecord.h" +#include "SimDataFormats/Track/interface/SimTrackContainer.h" +#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" +#include "SimDataFormats/Vertex/interface/SimVertexContainer.h" +#include "Validation/MuonHits/interface/MuonHitHelper.h" + +#include +#include +#include + +class MuonSimHitMatcher { + public: + // constructor + MuonSimHitMatcher(const edm::ParameterSet& iPS, edm::ConsumesCollector&& iC); + + // destructor + ~MuonSimHitMatcher() {} + + // initialize the event + void init(const edm::Event& e, const edm::EventSetup& eventSetup); + + // do the matching + void match(const SimTrack& t, const SimVertex& v); + + // access to all the Muon SimHits (use MuonSubdetId::SubSystem) + const edm::PSimHitContainer& simHits(int) const; + + // partitions' detIds with SimHits + std::set detIds(int type = 0) const; + + // chamber detIds with SimHits + std::set chamberIds(int type = 0) const; + + // simhits from a particular partition, chamber + const edm::PSimHitContainer& hitsInDetId(unsigned int) const; + const edm::PSimHitContainer& hitsInChamber(unsigned int) const; + + // calculate Global average position for a provided collection of simhits + GlobalPoint simHitsMeanPosition(const edm::PSimHitContainer& sim_hits) const; + + // calculate Global average momentum for a provided collection of simhits in + // CSC + GlobalVector simHitsMeanMomentum(const edm::PSimHitContainer& sim_hits) const; + + // calculate the average position at the second station + GlobalPoint simHitsMeanPositionStation(int n) const; + + protected: + std::vector getIdsOfSimTrackShower( + unsigned trk_id, const edm::SimTrackContainer& simTracks, + const edm::SimVertexContainer& simVertices); + + bool verboseSimTrack_; + bool simMuOnly_; + bool discardEleHits_; + bool verbose_; + bool hasGeometry_; + + const TrackingGeometry* geometry_; + + edm::EDGetTokenT simVertexInput_; + edm::EDGetTokenT simTrackInput_; + edm::EDGetTokenT simHitInput_; + + edm::Handle simTracksH_; + edm::Handle simVerticesH_; + edm::Handle simHitsH_; + + edm::SimTrackContainer simTracks_; + edm::SimVertexContainer simVertices_; + // input collection + edm::PSimHitContainer simHits_; + + std::vector track_ids_; + std::map trkid_to_index_; + + edm::PSimHitContainer no_hits_; + + // selected hits + edm::PSimHitContainer hits_; + std::map detid_to_hits_; + std::map chamber_to_hits_; + + edm::ParameterSet simHitPSet_; +}; + +#endif diff --git a/Validation/MuonHits/interface/RPCSimHitMatcher.h b/Validation/MuonHits/interface/RPCSimHitMatcher.h new file mode 100644 index 0000000000000..bc078327d3498 --- /dev/null +++ b/Validation/MuonHits/interface/RPCSimHitMatcher.h @@ -0,0 +1,50 @@ +#ifndef Validation_MuonHits_RPCSimHitMatcher_h +#define Validation_MuonHits_RPCSimHitMatcher_h + +/**\class RPCSimHitMatcher + + Description: Matching of RPC SimHit to SimTrack + + Author: Sven Dildick (TAMU), Tao Huang (TAMU) +*/ + +#include "Geometry/RPCGeometry/interface/RPCGeometry.h" +#include "Validation/MuonHits/interface/MuonSimHitMatcher.h" + +class RPCSimHitMatcher : public MuonSimHitMatcher { + public: + // constructor + RPCSimHitMatcher(const edm::ParameterSet& iPS, edm::ConsumesCollector&& iC); + + // destructor + ~RPCSimHitMatcher() {} + + // initialize the event + void init(const edm::Event& e, const edm::EventSetup& eventSetup); + + // do the matching + void match(const SimTrack& t, const SimVertex& v); + + // partitions' detIds with SimHits + std::set detIds(int type = MuonHitHelper::RPC_ALL) const; + + // chamber detIds with SimHits + std::set chamberIds(int type = MuonHitHelper::RPC_ALL) const; + + bool hitStation(int st) const; + + // number of stations with hits + int nStations() const; + + // calculate average strip number for a provided collection of simhits + float simHitsMeanStrip(const edm::PSimHitContainer& sim_hits) const; + + std::set hitStripsInDetId(unsigned int, int margin_n_strips = 0) const; + + private: + void matchSimHitsToSimTrack(); + + edm::ESHandle rpc_geom_; +}; + +#endif diff --git a/Validation/MuonHits/plugins/BuildFile.xml b/Validation/MuonHits/plugins/BuildFile.xml new file mode 100644 index 0000000000000..027a53565f395 --- /dev/null +++ b/Validation/MuonHits/plugins/BuildFile.xml @@ -0,0 +1,12 @@ + + + + + + + + + + + + diff --git a/Validation/MuonHits/src/MuonSimHitsValidAnalyzer.cc b/Validation/MuonHits/plugins/MuonSimHitsValidAnalyzer.cc similarity index 99% rename from Validation/MuonHits/src/MuonSimHitsValidAnalyzer.cc rename to Validation/MuonHits/plugins/MuonSimHitsValidAnalyzer.cc index 355040194251c..c7cc0800f9624 100644 --- a/Validation/MuonHits/src/MuonSimHitsValidAnalyzer.cc +++ b/Validation/MuonHits/plugins/MuonSimHitsValidAnalyzer.cc @@ -1,4 +1,4 @@ -#include "Validation/MuonHits/src/MuonSimHitsValidAnalyzer.h" +#include "Validation/MuonHits/plugins/MuonSimHitsValidAnalyzer.h" #include "TFile.h" #include "TTree.h" @@ -504,6 +504,3 @@ void MuonSimHitsValidAnalyzer::fillDT(const edm::Event& iEvent, edm::LogInfo("MuonSimHitsValidAnalyzer::fillDT") << eventout << "\n"; return; } - - - diff --git a/Validation/MuonHits/src/MuonSimHitsValidAnalyzer.h b/Validation/MuonHits/plugins/MuonSimHitsValidAnalyzer.h similarity index 87% rename from Validation/MuonHits/src/MuonSimHitsValidAnalyzer.h rename to Validation/MuonHits/plugins/MuonSimHitsValidAnalyzer.h index 8f6c3426d12f8..3bbf5fd382201 100644 --- a/Validation/MuonHits/src/MuonSimHitsValidAnalyzer.h +++ b/Validation/MuonHits/plugins/MuonSimHitsValidAnalyzer.h @@ -1,5 +1,5 @@ -#ifndef MuonSimHitsValidAnalyzer_h -#define MuonSimHitsValidAnalyzer_h +#ifndef Validation_MuonHits_MuonSimHitsValidAnalyzer_h +#define Validation_MuonHits_MuonSimHitsValidAnalyzer_h /// framework & common header files #include "FWCore/Framework/interface/Frameworkfwd.h" @@ -12,13 +12,11 @@ #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/ServiceRegistry/interface/Service.h" -#include "Geometry/CommonDetUnit/interface/GeomDet.h" #include "DataFormats/GeometrySurface/interface/BoundPlane.h" -#include "DataFormats/DetId/interface/DetId.h" -#include -#include -#include +#include "DQMServices/Core/interface/DQMStore.h" +#include "DQMServices/Core/interface/MonitorElement.h" +#include "DQMServices/Core/interface/DQMEDAnalyzer.h" /// muon CSC, DT and RPC geometry info #include "Geometry/Records/interface/MuonGeometryRecord.h" @@ -27,22 +25,13 @@ #include "Geometry/RPCGeometry/interface/RPCGeometry.h" #include "DataFormats/MuonDetId/interface/MuonSubdetId.h" -/// muon CSC detector id -#include "DataFormats/MuonDetId/interface/CSCDetId.h" - /// data in edm::event -//#include "SimDataFormats/ValidationFormats/interface/PValidationFormats.h" #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" #include "SimDataFormats/Vertex/interface/SimVertexContainer.h" #include "SimDataFormats/Track/interface/SimTrackContainer.h" #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" -/// helper files -#include -#include "CLHEP/Units/GlobalSystemOfUnits.h" - #include - #include #include #include diff --git a/Validation/MuonHits/src/SealModule.cc b/Validation/MuonHits/plugins/SealModule.cc similarity index 75% rename from Validation/MuonHits/src/SealModule.cc rename to Validation/MuonHits/plugins/SealModule.cc index d1fc4046e5cce..19adcc2b71c64 100644 --- a/Validation/MuonHits/src/SealModule.cc +++ b/Validation/MuonHits/plugins/SealModule.cc @@ -2,8 +2,5 @@ #include "FWCore/Framework/interface/MakerMacros.h" #include "DQMServices/Core/interface/DQMStore.h" -#include "Validation/MuonHits/src/MuonSimHitsValidAnalyzer.h" - - - +#include "Validation/MuonHits/plugins/MuonSimHitsValidAnalyzer.h" DEFINE_FWK_MODULE(MuonSimHitsValidAnalyzer); diff --git a/Validation/MuonHits/python/muonSimHitMatcherPSet.py b/Validation/MuonHits/python/muonSimHitMatcherPSet.py new file mode 100644 index 0000000000000..42797c65be87b --- /dev/null +++ b/Validation/MuonHits/python/muonSimHitMatcherPSet.py @@ -0,0 +1,45 @@ +import FWCore.ParameterSet.Config as cms + +muonSimHitMatcherPSet = cms.PSet( + simTrack = cms.PSet( + verbose = cms.int32(0), + inputTag = cms.InputTag("g4SimHits"), + ), + simVertex = cms.PSet( + verbose = cms.int32(0), + inputTag = cms.InputTag("g4SimHits"), + ), + gemSimHit = cms.PSet( + verbose = cms.int32(0), + inputTag = cms.InputTag("g4SimHits", "MuonGEMHits"), + simMuOnly = cms.bool(True), + discardEleHits = cms.bool(True), + ), + me0SimHit = cms.PSet( + verbose = cms.int32(0), + inputTag = cms.InputTag("g4SimHits", "MuonME0Hits"), + simMuOnly = cms.bool(True), + discardEleHits = cms.bool(True), + minNHitsChamber = cms.int32(4), + ), + rpcSimHit = cms.PSet( + verbose = cms.int32(0), + inputTag = cms.InputTag("g4SimHits", "MuonRPCHits"), + simMuOnly = cms.bool(True), + discardEleHits = cms.bool(True), + ), + cscSimHit = cms.PSet( + verbose = cms.int32(0), + inputTag = cms.InputTag("g4SimHits", "MuonCSCHits"), + simMuOnly = cms.bool(True), + discardEleHits = cms.bool(True), + minNHitsChamber = cms.int32(4), + ), + dtSimHit = cms.PSet( + verbose = cms.int32(0), + inputTag = cms.InputTag("g4SimHits", "MuonDTHits"), + simMuOnly = cms.bool(True), + discardEleHits = cms.bool(True), + minNHitsChamber = cms.int32(4), + ) +) diff --git a/Validation/MuonHits/src/CSCSimHitMatcher.cc b/Validation/MuonHits/src/CSCSimHitMatcher.cc new file mode 100644 index 0000000000000..4332b42dc8fe8 --- /dev/null +++ b/Validation/MuonHits/src/CSCSimHitMatcher.cc @@ -0,0 +1,371 @@ +#include "Validation/MuonHits/interface/CSCSimHitMatcher.h" + +using namespace std; + +CSCSimHitMatcher::CSCSimHitMatcher(const edm::ParameterSet& ps, + edm::ConsumesCollector&& iC) + : MuonSimHitMatcher(ps, std::move(iC)) { + simHitPSet_ = ps.getParameterSet("cscSimHit"); + verbose_ = simHitPSet_.getParameter("verbose"); + simMuOnly_ = simHitPSet_.getParameter("simMuOnly"); + discardEleHits_ = simHitPSet_.getParameter("discardEleHits"); + + simHitInput_ = iC.consumes( + simHitPSet_.getParameter("inputTag")); +} + +/// initialize the event +void CSCSimHitMatcher::init(const edm::Event& iEvent, + const edm::EventSetup& iSetup) { + iSetup.get().get(csc_geom_); + if (csc_geom_.isValid()) { + geometry_ = &*csc_geom_; + } else { + hasGeometry_ = false; + edm::LogWarning("CSCSimHitMatcher") + << "+++ Info: CSC geometry is unavailable. +++\n"; + } + MuonSimHitMatcher::init(iEvent, iSetup); +} + +/// do the matching +void CSCSimHitMatcher::match(const SimTrack& track, const SimVertex& vertex) { + // instantiates the track ids and simhits + MuonSimHitMatcher::match(track, vertex); + + if (hasGeometry_) { + matchSimHitsToSimTrack(); + + if (verbose_) { + edm::LogInfo("CSCSimHitMatcher") + << "nTrackIds " << track_ids_.size() << " nSelectedCSCSimHits " + << hits_.size() << endl; + edm::LogInfo("CSCSimHitMatcher") + << "detids CSC " << detIds(0).size() << endl; + + for (const auto& id : detIds(0)) { + const auto& simhits = hitsInDetId(id); + const auto& simhits_gp = simHitsMeanPosition(simhits); + const auto& strips = hitStripsInDetId(id); + CSCDetId cscid(id); + if (cscid.station() == 1 and (cscid.ring() == 1 or cscid.ring() == 4)) { + edm::LogInfo("CSCSimHitMatcher") + << "cscdetid " << CSCDetId(id) << ": " << simhits.size() << " " + << simhits_gp.phi() << " " << detid_to_hits_[id].size() << endl; + edm::LogInfo("CSCSimHitMatcher") + << "nStrip " << strips.size() << endl; + edm::LogInfo("CSCSimHitMatcher") << "strips : "; + for (const auto& p : strips) { + edm::LogInfo("CSCSimHitMatcher") << p; + } + } + } + } + } +} + +void CSCSimHitMatcher::matchSimHitsToSimTrack() { + for (const auto& track_id : track_ids_) { + for (const auto& h : simHits_) { + if (h.trackId() != track_id) continue; + int pdgid = h.particleType(); + if (simMuOnly_ && std::abs(pdgid) != 13) continue; + // discard electron hits in the CSC chambers + if (discardEleHits_ && pdgid == 11) continue; + + const CSCDetId& layer_id(h.detUnitId()); + hits_.push_back(h); + detid_to_hits_[h.detUnitId()].push_back(h); + chamber_to_hits_[layer_id.chamberId().rawId()].push_back(h); + } + } +} + +std::set CSCSimHitMatcher::detIds(int csc_type) const { + std::set result; + for (const auto& p : detid_to_hits_) { + const auto& id = p.first; + if (csc_type > 0) { + CSCDetId detId(id); + if (MuonHitHelper::toCSCType(detId.station(), detId.ring()) != csc_type) + continue; + } + result.insert(id); + } + return result; +} + +std::set CSCSimHitMatcher::chamberIds(int csc_type) const { + std::set result; + for (const auto& p : chamber_to_hits_) { + const auto& id = p.first; + if (csc_type > 0) { + CSCDetId detId(id); + if (MuonHitHelper::toCSCType(detId.station(), detId.ring()) != csc_type) + continue; + } + result.insert(id); + } + return result; +} + +int CSCSimHitMatcher::nLayersWithHitsInChamber(unsigned int detid) const { + set layers_with_hits; + for (const auto& h : MuonSimHitMatcher::hitsInChamber(detid)) { + const CSCDetId& idd(h.detUnitId()); + layers_with_hits.insert(idd.layer()); + } + return layers_with_hits.size(); +} + +bool CSCSimHitMatcher::hitStation(int st, int nlayers) const { + int nst = 0; + for (const auto& ddt : chamberIds()) { + const CSCDetId id(ddt); + int ri(id.ring()); + if (id.station() != st) continue; + + // ME1/a case - check the ME1/b chamber + if (st == 1 and ri == 4) { + CSCDetId idME1b(id.endcap(), id.station(), 1, id.chamber(), 0); + const int nl1a(nLayersWithHitsInChamber(id.rawId())); + const int nl1b(nLayersWithHitsInChamber(idME1b.rawId())); + if (nl1a + nl1b < nlayers) continue; + ++nst; + } + // ME1/b case - check the ME1/a chamber + else if (st == 1 and ri == 1) { + CSCDetId idME1a(id.endcap(), id.station(), 4, id.chamber(), 0); + const int nl1a(nLayersWithHitsInChamber(idME1a.rawId())); + const int nl1b(nLayersWithHitsInChamber(id.rawId())); + if (nl1a + nl1b < nlayers) continue; + ++nst; + } + // default case + else { + const int nl(nLayersWithHitsInChamber(id.rawId())); + if (nl < nlayers) continue; + ++nst; + } + } + return nst; +} + +int CSCSimHitMatcher::nStations(int nlayers) const { + return (hitStation(1, nlayers) + hitStation(2, nlayers) + + hitStation(3, nlayers) + hitStation(4, nlayers)); +} + +float CSCSimHitMatcher::LocalBendingInChamber(unsigned int detid) const { + const CSCDetId cscid(detid); + if (nLayersWithHitsInChamber(detid) < 6) return -100; + float phi_layer1 = -10; + float phi_layer6 = 10; + + if (cscid.station() == 1 and (cscid.ring() == 1 or cscid.ring() == 4)) { + // phi in layer 1 + const CSCDetId cscid1a(cscid.endcap(), cscid.station(), 4, cscid.chamber(), + 1); + const CSCDetId cscid1b(cscid.endcap(), cscid.station(), 1, cscid.chamber(), + 1); + const edm::PSimHitContainer& hits1a = + MuonSimHitMatcher::hitsInDetId(cscid1a.rawId()); + const edm::PSimHitContainer& hits1b = + MuonSimHitMatcher::hitsInDetId(cscid1b.rawId()); + const GlobalPoint& gp1a = + simHitsMeanPosition(MuonSimHitMatcher::hitsInDetId(cscid1a.rawId())); + const GlobalPoint& gp1b = + simHitsMeanPosition(MuonSimHitMatcher::hitsInDetId(cscid1b.rawId())); + if (!hits1a.empty() and !hits1b.empty()) + phi_layer1 = (gp1a.phi() + gp1b.phi()) / 2.0; + else if (!hits1a.empty()) + phi_layer1 = gp1a.phi(); + else if (!hits1b.empty()) + phi_layer1 = gp1b.phi(); + + // phi in layer 6 + const CSCDetId cscid6a(cscid.endcap(), cscid.station(), 4, cscid.chamber(), + 6); + const CSCDetId cscid6b(cscid.endcap(), cscid.station(), 1, cscid.chamber(), + 6); + const edm::PSimHitContainer& hits6a = + MuonSimHitMatcher::hitsInDetId(cscid6a.rawId()); + const edm::PSimHitContainer& hits6b = + MuonSimHitMatcher::hitsInDetId(cscid6b.rawId()); + const GlobalPoint& gp6a = + simHitsMeanPosition(MuonSimHitMatcher::hitsInDetId(cscid6a.rawId())); + const GlobalPoint& gp6b = + simHitsMeanPosition(MuonSimHitMatcher::hitsInDetId(cscid6b.rawId())); + if (!hits6a.empty() and !hits6b.empty()) + phi_layer6 = (gp6a.phi() + gp6b.phi()) / 2.0; + else if (!hits6a.empty()) + phi_layer6 = gp6a.phi(); + else if (!hits6b.empty()) + phi_layer6 = gp6b.phi(); + + } else { + // phi in layer 1 + const CSCDetId cscid1(cscid.endcap(), cscid.station(), cscid.ring(), + cscid.chamber(), 1); + const GlobalPoint& gp1 = + simHitsMeanPosition(MuonSimHitMatcher::hitsInDetId(cscid1.rawId())); + phi_layer1 = gp1.phi(); + + // phi in layer 6 + const CSCDetId cscid6(cscid.endcap(), cscid.station(), cscid.ring(), + cscid.chamber(), 6); + const GlobalPoint& gp6 = + simHitsMeanPosition(MuonSimHitMatcher::hitsInDetId(cscid6.rawId())); + phi_layer6 = gp6.phi(); + } + return deltaPhi(phi_layer6, phi_layer1); +} + +float CSCSimHitMatcher::simHitsMeanStrip( + const edm::PSimHitContainer& sim_hits) const { + if (sim_hits.empty()) return -1.f; + + float sums = 0.f; + size_t n = 0; + for (const auto& h : sim_hits) { + const LocalPoint& lp = h.entryPoint(); + float s; + const auto& d = h.detUnitId(); + s = dynamic_cast(geometry_) + ->layer(d) + ->geometry() + ->strip(lp); + // convert to half-strip: + s *= 2.; + sums += s; + ++n; + } + if (n == 0) return -1.f; + return sums / n; +} + +float CSCSimHitMatcher::simHitsMeanWG( + const edm::PSimHitContainer& sim_hits) const { + if (sim_hits.empty()) return -1.f; + + float sums = 0.f; + size_t n = 0; + for (const auto& h : sim_hits) { + const LocalPoint& lp = h.entryPoint(); + float s; + const auto& d = h.detUnitId(); + // find nearest wire + const auto& layerG( + dynamic_cast(geometry_)->layer(d)->geometry()); + int nearestWire(layerG->nearestWire(lp)); + // then find the corresponding wire group + s = layerG->wireGroup(nearestWire); + sums += s; + ++n; + } + if (n == 0) return -1.f; + return sums / n; +} + +std::set CSCSimHitMatcher::hitStripsInDetId(unsigned int detid, + int margin_n_strips) const { + set result; + const auto& simhits = MuonSimHitMatcher::hitsInDetId(detid); + CSCDetId id(detid); + int max_nstrips = dynamic_cast(geometry_) + ->layer(id) + ->geometry() + ->numberOfStrips(); + for (const auto& h : simhits) { + const LocalPoint& lp = h.entryPoint(); + int central_strip = dynamic_cast(geometry_) + ->layer(id) + ->geometry() + ->nearestStrip(lp); + int smin = central_strip - margin_n_strips; + smin = (smin > 0) ? smin : 1; + int smax = central_strip + margin_n_strips; + smax = (smax <= max_nstrips) ? smax : max_nstrips; + for (int ss = smin; ss <= smax; ++ss) result.insert(ss); + } + return result; +} + +std::set CSCSimHitMatcher::hitWiregroupsInDetId(unsigned int detid, + int margin_n_wg) const { + set result; + const auto& simhits = MuonSimHitMatcher::hitsInDetId(detid); + CSCDetId id(detid); + int max_n_wg = dynamic_cast(geometry_) + ->layer(id) + ->geometry() + ->numberOfWireGroups(); + for (const auto& h : simhits) { + const LocalPoint& lp = h.entryPoint(); + const auto& layer_geo = + dynamic_cast(geometry_)->layer(id)->geometry(); + int central_wg = layer_geo->wireGroup(layer_geo->nearestWire(lp)); + int wg_min = central_wg - margin_n_wg; + wg_min = (wg_min > 0) ? wg_min : 1; + int wg_max = central_wg + margin_n_wg; + wg_max = (wg_max <= max_n_wg) ? wg_max : max_n_wg; + for (int wg = wg_min; wg <= wg_max; ++wg) result.insert(wg); + } + return result; +} + +int CSCSimHitMatcher::nCoincidenceChambers(int min_n_layers) const { + int result = 0; + const auto& chamber_ids = chamberIds(0); + for (const auto& id : chamber_ids) { + if (nLayersWithHitsInChamber(id) >= min_n_layers) result += 1; + } + return result; +} + +void CSCSimHitMatcher::chamberIdsToString( + const std::set& set) const { + for (const auto& p : set) { + CSCDetId detId(p); + edm::LogInfo("CSCSimHitMatcher") << " " << detId << "\n"; + } +} + +std::set CSCSimHitMatcher::chamberIdsStation(int station) const { + set result; + switch (station) { + case 1: { + const auto& p1(chamberIds(MuonHitHelper::CSC_ME1a)); + const auto& p2(chamberIds(MuonHitHelper::CSC_ME1b)); + const auto& p3(chamberIds(MuonHitHelper::CSC_ME12)); + const auto& p4(chamberIds(MuonHitHelper::CSC_ME13)); + result.insert(p1.begin(), p1.end()); + result.insert(p2.begin(), p2.end()); + result.insert(p3.begin(), p3.end()); + result.insert(p4.begin(), p4.end()); + break; + } + case 2: { + const auto& p1(chamberIds(MuonHitHelper::CSC_ME21)); + const auto& p2(chamberIds(MuonHitHelper::CSC_ME22)); + result.insert(p1.begin(), p1.end()); + result.insert(p2.begin(), p2.end()); + break; + } + case 3: { + const auto& p1(chamberIds(MuonHitHelper::CSC_ME31)); + const auto& p2(chamberIds(MuonHitHelper::CSC_ME32)); + result.insert(p1.begin(), p1.end()); + result.insert(p2.begin(), p2.end()); + break; + } + case 4: { + const auto& p1(chamberIds(MuonHitHelper::CSC_ME41)); + const auto& p2(chamberIds(MuonHitHelper::CSC_ME42)); + result.insert(p1.begin(), p1.end()); + result.insert(p2.begin(), p2.end()); + break; + } + }; + return result; +} diff --git a/Validation/MuonHits/src/DTSimHitMatcher.cc b/Validation/MuonHits/src/DTSimHitMatcher.cc new file mode 100644 index 0000000000000..b98bff00f48d2 --- /dev/null +++ b/Validation/MuonHits/src/DTSimHitMatcher.cc @@ -0,0 +1,374 @@ +#include "Validation/MuonHits/interface/DTSimHitMatcher.h" + +using namespace std; + +DTSimHitMatcher::DTSimHitMatcher(const edm::ParameterSet& ps, + edm::ConsumesCollector&& iC) + : MuonSimHitMatcher(ps, std::move(iC)) { + simHitPSet_ = ps.getParameterSet("dtSimHit"); + verbose_ = simHitPSet_.getParameter("verbose"); + simMuOnly_ = simHitPSet_.getParameter("simMuOnly"); + discardEleHits_ = simHitPSet_.getParameter("discardEleHits"); + + simHitInput_ = iC.consumes( + simHitPSet_.getParameter("inputTag")); +} + +/// initialize the event +void DTSimHitMatcher::init(const edm::Event& iEvent, + const edm::EventSetup& iSetup) { + iSetup.get().get(dt_geom_); + if (dt_geom_.isValid()) { + geometry_ = &*dt_geom_; + } else { + hasGeometry_ = false; + edm::LogWarning("DTSimHitMatcher") + << "+++ Info: DT geometry is unavailable. +++\n"; + } + MuonSimHitMatcher::init(iEvent, iSetup); +} + +/// do the matching +void DTSimHitMatcher::match(const SimTrack& track, const SimVertex& vertex) { + // instantiates the track ids and simhits + MuonSimHitMatcher::match(track, vertex); + + if (hasGeometry_) { + matchSimHitsToSimTrack(); + + if (verbose_) { + edm::LogInfo("DTSimHitMatcher") + << "nTrackIds " << track_ids_.size() << " nSelectedDTSimHits " + << hits_.size() << endl; + edm::LogInfo("DTSimHitMatcher") + << "detids DT " << detIds(0).size() << endl; + + const auto& dt_det_ids = detIds(0); + for (const auto& id : dt_det_ids) { + const auto& dt_simhits = MuonSimHitMatcher::hitsInDetId(id); + const auto& dt_simhits_gp = simHitsMeanPosition(dt_simhits); + edm::LogInfo("DTSimHitMatcher") + << "DTWireId " << DTWireId(id) << ": nHits " << dt_simhits.size() + << " eta " << dt_simhits_gp.eta() << " phi " << dt_simhits_gp.phi() + << " nCh " << chamber_to_hits_[id].size() << endl; + } + } + } +} + +void DTSimHitMatcher::matchSimHitsToSimTrack() { + for (const auto& track_id : track_ids_) { + for (const auto& h : simHits_) { + if (h.trackId() != track_id) continue; + int pdgid = h.particleType(); + if (simMuOnly_ && std::abs(pdgid) != 13) continue; + // discard electron hits in the DT chambers + if (discardEleHits_ && pdgid == 11) continue; + + const DTWireId layer_id(h.detUnitId()); + detid_to_hits_[h.detUnitId()].push_back(h); + hits_.push_back(h); + layer_to_hits_[layer_id.layerId().rawId()].push_back(h); + superlayer_to_hits_[layer_id.superlayerId().rawId()].push_back(h); + chamber_to_hits_[layer_id.chamberId().rawId()].push_back(h); + } + } +} + +std::set DTSimHitMatcher::detIds(int dt_type) const { + std::set result; + for (const auto& p : detid_to_hits_) { + const auto& id = p.first; + if (dt_type > 0) { + DTWireId detId(id); + if (MuonHitHelper::toDTType(detId.wheel(), detId.station()) != dt_type) + continue; + } + result.insert(id); + } + return result; +} + +std::set DTSimHitMatcher::chamberIds(int dt_type) const { + std::set result; + for (const auto& p : chamber_to_hits_) { + const auto& id = p.first; + if (dt_type > 0) { + DTChamberId detId(id); + if (MuonHitHelper::toDTType(detId.wheel(), detId.station()) != dt_type) + continue; + } + result.insert(id); + } + return result; +} + +std::set DTSimHitMatcher::layerIds() const { + std::set result; + for (const auto& p : layer_to_hits_) result.insert(p.first); + return result; +} + +std::set DTSimHitMatcher::superlayerIds() const { + std::set result; + for (const auto& p : superlayer_to_hits_) result.insert(p.first); + return result; +} + +const edm::PSimHitContainer& DTSimHitMatcher::hitsInLayer( + unsigned int detid) const { + if (!MuonHitHelper::isDT(detid)) return no_hits_; + + const DTWireId id(detid); + if (layer_to_hits_.find(id.layerId().rawId()) == layer_to_hits_.end()) + return no_hits_; + return layer_to_hits_.at(id.layerId().rawId()); +} + +const edm::PSimHitContainer& DTSimHitMatcher::hitsInSuperLayer( + unsigned int detid) const { + if (!MuonHitHelper::isDT(detid)) return no_hits_; + + const DTWireId id(detid); + if (superlayer_to_hits_.find(id.superlayerId().rawId()) == + superlayer_to_hits_.end()) + return no_hits_; + return superlayer_to_hits_.at(id.superlayerId().rawId()); +} + +const edm::PSimHitContainer& DTSimHitMatcher::hitsInChamber( + unsigned int detid) const { + if (!MuonHitHelper::isDT(detid)) return no_hits_; + + const DTWireId id(detid); + if (chamber_to_hits_.find(id.chamberId().rawId()) == chamber_to_hits_.end()) + return no_hits_; + return chamber_to_hits_.at(id.chamberId().rawId()); +} + +bool DTSimHitMatcher::hitStation(int st, int nsuperlayers, int nlayers) const { + int nst = 0; + for (const auto& ddt : chamberIds()) { + const DTChamberId id(ddt); + if (id.station() != st) continue; + + // require at least 1 superlayer + const int nsl(nSuperLayersWithHitsInChamber(id.rawId())); + if (nsl < nsuperlayers) continue; + + // require at least 3 layers hit per chamber + const int nl(nLayersWithHitsInChamber(id.rawId())); + if (nl < nlayers) continue; + ++nst; + } + return nst; +} + +int DTSimHitMatcher::nStations(int nsuperlayers, int nlayers) const { + return (hitStation(1, nsuperlayers, nlayers) + + hitStation(2, nsuperlayers, nlayers) + + hitStation(3, nsuperlayers, nlayers) + + hitStation(4, nsuperlayers, nlayers)); +} + +int DTSimHitMatcher::nCellsWithHitsInLayer(unsigned int detid) const { + set layers_with_hits; + const auto& hits = hitsInLayer(detid); + for (const auto& h : hits) { + if (MuonHitHelper::isDT(detid)) { + const DTWireId idd(h.detUnitId()); + layers_with_hits.insert(idd.wire()); + } + } + return layers_with_hits.size(); +} + +int DTSimHitMatcher::nLayersWithHitsInSuperLayer(unsigned int detid) const { + set layers_with_hits; + const auto& hits = hitsInSuperLayer(detid); + for (const auto& h : hits) { + if (MuonHitHelper::isDT(detid)) { + const DTLayerId idd(h.detUnitId()); + layers_with_hits.insert(idd.layer()); + } + } + return layers_with_hits.size(); +} + +int DTSimHitMatcher::nSuperLayersWithHitsInChamber(unsigned int detid) const { + set sl_with_hits; + const auto& hits = MuonSimHitMatcher::hitsInChamber(detid); + for (const auto& h : hits) { + if (MuonHitHelper::isDT(detid)) { + const DTSuperLayerId idd(h.detUnitId()); + sl_with_hits.insert(idd.superLayer()); + } + } + return sl_with_hits.size(); +} + +int DTSimHitMatcher::nLayersWithHitsInChamber(unsigned int detid) const { + int nLayers = 0; + const auto& superLayers(dynamic_cast(geometry_) + ->chamber(DTChamberId(detid)) + ->superLayers()); + for (const auto& sl : superLayers) { + nLayers += nLayersWithHitsInSuperLayer(sl->id().rawId()); + } + return nLayers; +} +float DTSimHitMatcher::simHitsMeanWire( + const edm::PSimHitContainer& sim_hits) const { + if (sim_hits.empty()) return -1.f; + + float sums = 0.f; + size_t n = 0; + for (const auto& h : sim_hits) { + const LocalPoint& lp = h.entryPoint(); + float s; + const auto& d = h.detUnitId(); + if (MuonHitHelper::isDT(d)) { + // find nearest wire + s = dynamic_cast(geometry_) + ->layer(DTLayerId(d)) + ->specificTopology() + .channel(lp); + } else + continue; + sums += s; + ++n; + } + if (n == 0) return -1.f; + return sums / n; +} + +std::set DTSimHitMatcher::hitWiresInDTLayerId( + unsigned int detid, int margin_n_wires) const { + set result; + if (MuonHitHelper::isDT(detid)) { + DTLayerId id(detid); + int max_nwires = dynamic_cast(geometry_) + ->layer(id) + ->specificTopology() + .channels(); + for (int wn = 0; wn <= max_nwires; ++wn) { + DTWireId wid(id, wn); + for (const auto& h : MuonSimHitMatcher::hitsInDetId(wid.rawId())) { + if (verbose_) + edm::LogInfo("DTSimHitMatcher") + << "central DTWireId " << wid << " simhit " << h << endl; + int smin = wn - margin_n_wires; + smin = (smin > 0) ? smin : 1; + int smax = wn + margin_n_wires; + smax = (smax <= max_nwires) ? smax : max_nwires; + for (int ss = smin; ss <= smax; ++ss) { + DTWireId widd(id, ss); + if (verbose_) + edm::LogInfo("DTSimHitMatcher") + << "\tadding DTWireId to collection " << widd << endl; + result.insert(widd.rawId()); + } + } + } + } + return result; +} + +std::set DTSimHitMatcher::hitWiresInDTSuperLayerId( + unsigned int detid, int margin_n_wires) const { + set result; + const auto& layers(dynamic_cast(geometry_) + ->superLayer(DTSuperLayerId(detid)) + ->layers()); + for (const auto& l : layers) { + if (verbose_) + edm::LogInfo("DTSimHitMatcher") + << "hitWiresInDTSuperLayerId::l id " << l->id() << endl; + const auto& p(hitWiresInDTLayerId(l->id().rawId(), margin_n_wires)); + result.insert(p.begin(), p.end()); + } + return result; +} + +std::set DTSimHitMatcher::hitWiresInDTChamberId( + unsigned int detid, int margin_n_wires) const { + set result; + const auto& superLayers(dynamic_cast(geometry_) + ->chamber(DTChamberId(detid)) + ->superLayers()); + for (const auto& sl : superLayers) { + if (verbose_) + edm::LogInfo("DTSimHitMatcher") + << "hitWiresInDTChamberId::sl id " << sl->id() << endl; + const auto& p(hitWiresInDTSuperLayerId(sl->id().rawId(), margin_n_wires)); + result.insert(p.begin(), p.end()); + } + return result; +} + +void DTSimHitMatcher::dtChamberIdsToString( + const std::set& set) const { + for (const auto& p : set) { + DTChamberId detId(p); + edm::LogInfo("DTSimHitMatcher") << " " << detId << "\n"; + } +} + +std::set DTSimHitMatcher::chamberIdsStation(int station) const { + set result; + switch (station) { + case 1: { + const auto& p1(chamberIds(MuonHitHelper::DT_MB21p)); + const auto& p2(chamberIds(MuonHitHelper::DT_MB11p)); + const auto& p3(chamberIds(MuonHitHelper::DT_MB01)); + const auto& p4(chamberIds(MuonHitHelper::DT_MB11n)); + const auto& p5(chamberIds(MuonHitHelper::DT_MB21n)); + result.insert(p1.begin(), p1.end()); + result.insert(p2.begin(), p2.end()); + result.insert(p3.begin(), p3.end()); + result.insert(p4.begin(), p4.end()); + result.insert(p5.begin(), p5.end()); + break; + } + case 2: { + const auto& p1(chamberIds(MuonHitHelper::DT_MB22p)); + const auto& p2(chamberIds(MuonHitHelper::DT_MB12p)); + const auto& p3(chamberIds(MuonHitHelper::DT_MB02)); + const auto& p4(chamberIds(MuonHitHelper::DT_MB12n)); + const auto& p5(chamberIds(MuonHitHelper::DT_MB22n)); + result.insert(p1.begin(), p1.end()); + result.insert(p2.begin(), p2.end()); + result.insert(p3.begin(), p3.end()); + result.insert(p4.begin(), p4.end()); + result.insert(p5.begin(), p5.end()); + break; + } + case 3: { + const auto& p1(chamberIds(MuonHitHelper::DT_MB23p)); + const auto& p2(chamberIds(MuonHitHelper::DT_MB13p)); + const auto& p3(chamberIds(MuonHitHelper::DT_MB03)); + const auto& p4(chamberIds(MuonHitHelper::DT_MB13n)); + const auto& p5(chamberIds(MuonHitHelper::DT_MB23n)); + result.insert(p1.begin(), p1.end()); + result.insert(p2.begin(), p2.end()); + result.insert(p3.begin(), p3.end()); + result.insert(p4.begin(), p4.end()); + result.insert(p5.begin(), p5.end()); + break; + } + case 4: { + const auto& p1(chamberIds(MuonHitHelper::DT_MB24p)); + const auto& p2(chamberIds(MuonHitHelper::DT_MB14p)); + const auto& p3(chamberIds(MuonHitHelper::DT_MB04)); + const auto& p4(chamberIds(MuonHitHelper::DT_MB14n)); + const auto& p5(chamberIds(MuonHitHelper::DT_MB24n)); + result.insert(p1.begin(), p1.end()); + result.insert(p2.begin(), p2.end()); + result.insert(p3.begin(), p3.end()); + result.insert(p4.begin(), p4.end()); + result.insert(p5.begin(), p5.end()); + break; + } + }; + return result; +} diff --git a/Validation/MuonHits/src/GEMSimHitMatcher.cc b/Validation/MuonHits/src/GEMSimHitMatcher.cc new file mode 100644 index 0000000000000..2f16bb04eac6b --- /dev/null +++ b/Validation/MuonHits/src/GEMSimHitMatcher.cc @@ -0,0 +1,360 @@ +#include "Validation/MuonHits/interface/GEMSimHitMatcher.h" + +using namespace std; + +GEMSimHitMatcher::GEMSimHitMatcher(const edm::ParameterSet& ps, + edm::ConsumesCollector&& iC) + : MuonSimHitMatcher(ps, std::move(iC)) { + simHitPSet_ = ps.getParameterSet("gemSimHit"); + verbose_ = simHitPSet_.getParameter("verbose"); + simMuOnly_ = simHitPSet_.getParameter("simMuOnly"); + discardEleHits_ = simHitPSet_.getParameter("discardEleHits"); + + simHitInput_ = iC.consumes( + simHitPSet_.getParameter("inputTag")); +} + +/// initialize the event +void GEMSimHitMatcher::init(const edm::Event& iEvent, + const edm::EventSetup& iSetup) { + iSetup.get().get(gem_geom_); + if (gem_geom_.isValid()) { + geometry_ = dynamic_cast(&*gem_geom_); + } else { + hasGeometry_ = false; + edm::LogWarning("GEMSimHitMatcher") + << "+++ Info: GEM geometry is unavailable. +++\n"; + } + MuonSimHitMatcher::init(iEvent, iSetup); +} + +/// do the matching +void GEMSimHitMatcher::match(const SimTrack& track, const SimVertex& vertex) { + // instantiates the track ids and simhits + MuonSimHitMatcher::match(track, vertex); + + if (hasGeometry_) { + matchSimHitsToSimTrack(); + + if (verbose_) { + edm::LogInfo("GEMSimHitMatcher") + << "nTrackIds " << track_ids_.size() << " nSelectedGEMSimHits " + << hits_.size() << endl; + edm::LogInfo("GEMSimHitMatcher") + << "detids GEM " << detIds(0).size() << endl; + + const auto& gem_ch_ids = detIds(); + for (const auto& id : gem_ch_ids) { + const auto& gem_simhits = MuonSimHitMatcher::hitsInDetId(id); + const auto& gem_simhits_gp = simHitsMeanPosition(gem_simhits); + edm::LogInfo("GEMSimHitMatcher") + << "gemchid " << GEMDetId(id) << ": nHits " << gem_simhits.size() + << " phi " << gem_simhits_gp.phi() << " nCh " + << chamber_to_hits_[id].size() << endl; + const auto& strips = hitStripsInDetId(id); + edm::LogInfo("GEMSimHitMatcher") << "nStrip " << strips.size() << endl; + edm::LogInfo("GEMSimHitMatcher") << "strips : "; + for (const auto& p : strips) { + edm::LogInfo("GEMSimHitMatcher") << p; + } + } + const auto& gem_sch_ids = superChamberIds(); + for (const auto& id : gem_sch_ids) { + const auto& gem_simhits = hitsInSuperChamber(id); + const auto& gem_simhits_gp = simHitsMeanPosition(gem_simhits); + edm::LogInfo("GEMSimHitMatcher") + << "gemschid " << GEMDetId(id) << ": " << nCoincidencePadsWithHits() + << " | " << gem_simhits.size() << " " << gem_simhits_gp.phi() << " " + << superchamber_to_hits_[id].size() << endl; + } + } + } +} + +void GEMSimHitMatcher::matchSimHitsToSimTrack() { + for (const auto& track_id : track_ids_) { + for (const auto& h : simHits_) { + if (h.trackId() != track_id) continue; + int pdgid = h.particleType(); + if (simMuOnly_ && std::abs(pdgid) != 13) continue; + // discard electron hits in the GEM chambers + if (discardEleHits_ && pdgid == 11) continue; + + const GEMDetId& p_id(h.detUnitId()); + detid_to_hits_[h.detUnitId()].push_back(h); + hits_.push_back(h); + chamber_to_hits_[p_id.chamberId().rawId()].push_back(h); + superchamber_to_hits_[p_id.superChamberId().rawId()].push_back(h); + } + } + + // find pads with hits + const auto& detids = detIds(); + // find 2-layer coincidence pads with hits + for (const auto& d : detids) { + GEMDetId id(d); + const auto& hits = hitsInDetId(d); + const auto& roll = + dynamic_cast(geometry_)->etaPartition(id); + // int max_npads = roll->npads(); + set pads; + for (const auto& h : hits) { + const LocalPoint& lp = h.entryPoint(); + pads.insert(1 + static_cast(roll->padTopology().channel(lp))); + } + detids_to_pads_[d] = pads; + } + + // find 2-layer coincidence pads with hits + for (const auto& d : detids) { + GEMDetId id1(d); + if (id1.layer() != 1) continue; + + // find pads with hits in layer1 + const auto& hits1 = hitsInDetId(d); + const auto& roll1 = + dynamic_cast(geometry_)->etaPartition(id1); + set pads1; + set pads2; + set copads; + + for (const auto& h : hits1) { + const LocalPoint& lp = h.entryPoint(); + pads1.insert(1 + static_cast(roll1->padTopology().channel(lp))); + if (verbose_) + edm::LogInfo("GEMSimHitMatcher") + << "GEMHits detid1 " << id1 << " pad1 " + << 1 + static_cast(roll1->padTopology().channel(lp)) + << std::endl; + } + + // find pads with hits in layer2 + for (const auto& d2 : detids) { + // staggered geometry???? improve here !! + GEMDetId id2(d2); + // does layer 2 has simhits? + if (id2.layer() != 2 or id2.region() != id1.region() or + id2.ring() != id1.ring() or id2.station() != id1.station() or + abs(id2.roll() - id1.roll()) > 1) + continue; + const auto& hits2 = hitsInDetId(id2()); + const auto& roll2 = + dynamic_cast(geometry_)->etaPartition(id2); + for (const auto& h : hits2) { + const LocalPoint& lp = h.entryPoint(); + pads2.insert(1 + static_cast(roll2->padTopology().channel(lp))); + if (verbose_) + edm::LogInfo("GEMSimHitMatcher") + << "GEMHits detid2 " << id2 << " pad2 " + << 1 + static_cast(roll2->padTopology().channel(lp)) + << std::endl; + } + } + + for (const auto& pad1 : pads1) { + for (const auto& pad2 : pads2) { + if (abs(pad1 - pad2) <= 2) { + if (copads.find(pad1) == copads.end()) copads.insert(pad1); + if (copads.find(pad2) == copads.end()) copads.insert(pad2); + } + } + } + + if (copads.empty()) continue; + + // detids here is layer1 id + detids_to_copads_[d] = copads; + } +} + +std::set GEMSimHitMatcher::detIds(int gem_type) const { + std::set result; + for (const auto& p : detid_to_hits_) { + const auto& id = p.first; + if (gem_type > 0) { + GEMDetId detId(id); + if (MuonHitHelper::toGEMType(detId.station(), detId.ring()) != gem_type) + continue; + } + result.insert(id); + } + return result; +} + +std::set GEMSimHitMatcher::detIdsCoincidences() const { + std::set result; + for (const auto& p : detids_to_copads_) result.insert(p.first); + return result; +} + +std::set GEMSimHitMatcher::chamberIds(int gem_type) const { + std::set result; + for (const auto& p : chamber_to_hits_) { + const auto& id = p.first; + if (gem_type > 0) { + GEMDetId detId(id); + if (MuonHitHelper::toGEMType(detId.station(), detId.ring()) != gem_type) + continue; + } + result.insert(id); + } + return result; +} + +std::set GEMSimHitMatcher::superChamberIds() const { + std::set result; + for (const auto& p : superchamber_to_hits_) result.insert(p.first); + return result; +} + +std::set GEMSimHitMatcher::superChamberIdsCoincidences() const { + std::set result; + for (const auto& p : detids_to_copads_) { + const GEMDetId& p_id(p.first); + result.insert(p_id.superChamberId().rawId()); + } + return result; +} + +const edm::PSimHitContainer& GEMSimHitMatcher::hitsInSuperChamber( + unsigned int detid) const { + if (MuonHitHelper::isGEM(detid)) { + const GEMDetId id(detid); + if (superchamber_to_hits_.find(id.chamberId().rawId()) == + superchamber_to_hits_.end()) + return no_hits_; + return superchamber_to_hits_.at(id.chamberId().rawId()); + } + + return no_hits_; +} + +int GEMSimHitMatcher::nLayersWithHitsInSuperChamber(unsigned int detid) const { + set layers_with_hits; + const auto& hits = hitsInSuperChamber(detid); + for (const auto& h : hits) { + const GEMDetId& idd(h.detUnitId()); + layers_with_hits.insert(idd.layer()); + } + return layers_with_hits.size(); +} + +bool GEMSimHitMatcher::hitStation(int st, int nlayers) const { + int nst = 0; + for (const auto& ddt : chamberIds()) { + const GEMDetId id(ddt); + if (id.station() != st) continue; + + const int nl(nLayersWithHitsInSuperChamber(id.rawId())); + if (nl < nlayers) continue; + ++nst; + } + return nst; +} + +int GEMSimHitMatcher::nStations(int nlayers) const { + return (hitStation(1, nlayers) + hitStation(2, nlayers)); +} + +float GEMSimHitMatcher::simHitsGEMCentralPosition( + const edm::PSimHitContainer& sim_hits) const { + if (sim_hits.empty()) return -0.0; // point "zero" + + float central = -0.0; + size_t n = 0; + for (const auto& h : sim_hits) { + LocalPoint lp(0., 0., 0.); // local central + GlobalPoint gp; + if (MuonHitHelper::isGEM(h.detUnitId())) { + gp = geometry_->idToDet(h.detUnitId())->surface().toGlobal(lp); + } + central = gp.perp(); + if (n >= 1) + edm::LogWarning("GEMSimHitMatcher") + << "warning! find more than one simhits in GEM chamber " << std::endl; + ++n; + } + + return central; +} + +float GEMSimHitMatcher::simHitsMeanStrip( + const edm::PSimHitContainer& sim_hits) const { + if (sim_hits.empty()) return -1.f; + + float sums = 0.f; + size_t n = 0; + for (const auto& h : sim_hits) { + const LocalPoint& lp = h.entryPoint(); + const auto& d = h.detUnitId(); + sums += + dynamic_cast(geometry_)->etaPartition(d)->strip(lp); + ++n; + } + if (n == 0) return -1.f; + return sums / n; +} + +std::set GEMSimHitMatcher::hitStripsInDetId(unsigned int detid, + int margin_n_strips) const { + set result; + const auto& simhits = MuonSimHitMatcher::hitsInDetId(detid); + GEMDetId id(detid); + int max_nstrips = + dynamic_cast(geometry_)->etaPartition(id)->nstrips(); + for (const auto& h : simhits) { + const LocalPoint& lp = h.entryPoint(); + int central_strip = + static_cast(dynamic_cast(geometry_) + ->etaPartition(id) + ->topology() + .channel(lp)); + int smin = central_strip - margin_n_strips; + smin = (smin > 0) ? smin : 1; + int smax = central_strip + margin_n_strips; + smax = (smax <= max_nstrips) ? smax : max_nstrips; + for (int ss = smin; ss <= smax; ++ss) result.insert(ss); + } + return result; +} + +std::set GEMSimHitMatcher::hitPadsInDetId(unsigned int detid) const { + set none; + if (detids_to_pads_.find(detid) == detids_to_pads_.end()) return none; + return detids_to_pads_.at(detid); +} + +std::set GEMSimHitMatcher::hitCoPadsInDetId(unsigned int detid) const { + set none; + if (detids_to_copads_.find(detid) == detids_to_copads_.end()) return none; + return detids_to_copads_.at(detid); +} + +std::set GEMSimHitMatcher::hitPartitions() const { + std::set result; + + const auto& detids = detIds(); + for (const auto& id : detids) { + GEMDetId idd(id); + result.insert(idd.roll()); + } + return result; +} + +int GEMSimHitMatcher::nPadsWithHits() const { + int result = 0; + const auto& pad_ids = detIds(); + for (const auto& id : pad_ids) { + result += hitPadsInDetId(id).size(); + } + return result; +} + +int GEMSimHitMatcher::nCoincidencePadsWithHits() const { + int result = 0; + const auto& copad_ids = detIdsCoincidences(); + for (const auto& id : copad_ids) { + result += hitCoPadsInDetId(id).size(); + } + return result; +} diff --git a/Validation/MuonHits/src/ME0SimHitMatcher.cc b/Validation/MuonHits/src/ME0SimHitMatcher.cc new file mode 100644 index 0000000000000..645e1bac17895 --- /dev/null +++ b/Validation/MuonHits/src/ME0SimHitMatcher.cc @@ -0,0 +1,215 @@ +#include "Validation/MuonHits/interface/ME0SimHitMatcher.h" + +using namespace std; + +ME0SimHitMatcher::ME0SimHitMatcher(const edm::ParameterSet& ps, + edm::ConsumesCollector&& iC) + : MuonSimHitMatcher(ps, std::move(iC)) { + simHitPSet_ = ps.getParameterSet("me0SimHit"); + verbose_ = simHitPSet_.getParameter("verbose"); + simMuOnly_ = simHitPSet_.getParameter("simMuOnly"); + discardEleHits_ = simHitPSet_.getParameter("discardEleHits"); + + simHitInput_ = iC.consumes( + simHitPSet_.getParameter("inputTag")); +} + +/// initialize the event +void ME0SimHitMatcher::init(const edm::Event& iEvent, + const edm::EventSetup& iSetup) { + iSetup.get().get(me0_geom_); + if (me0_geom_.isValid()) { + geometry_ = &*me0_geom_; + } else { + hasGeometry_ = false; + edm::LogWarning("ME0SimHitMatcher") + << "+++ Info: ME0 geometry is unavailable. +++\n"; + } + MuonSimHitMatcher::init(iEvent, iSetup); +} + +/// do the matching +void ME0SimHitMatcher::match(const SimTrack& track, const SimVertex& vertex) { + // instantiates the track ids and simhits + MuonSimHitMatcher::match(track, vertex); + + if (hasGeometry_) { + matchSimHitsToSimTrack(); + + if (verbose_) { + edm::LogInfo("ME0SimHitMatcher") + << "nTrackIds " << track_ids_.size() << " nSelectedME0SimHits " + << hits_.size() << endl; + edm::LogInfo("ME0SimHitMatcher") + << "detids ME0 " << detIds().size() << endl; + + const auto& me0_ch_ids = detIds(); + for (const auto& id : me0_ch_ids) { + const auto& me0_simhits = MuonSimHitMatcher::hitsInChamber(id); + const auto& me0_simhits_gp = simHitsMeanPosition(me0_simhits); + edm::LogInfo("ME0SimHitMatcher") + << "me0chid " << ME0DetId(id) << ": nHits " << me0_simhits.size() + << " phi " << me0_simhits_gp.phi() << " nCh " + << chamber_to_hits_[id].size() << endl; + const auto& strips = hitStripsInDetId(id); + edm::LogInfo("ME0SimHitMatcher") << "nStrip " << strips.size() << endl; + edm::LogInfo("ME0SimHitMatcher") << "strips : "; + for (const auto& p : strips) { + edm::LogInfo("ME0SimHitMatcher") << p; + } + } + } + } +} + +void ME0SimHitMatcher::matchSimHitsToSimTrack() { + for (const auto& track_id : track_ids_) { + for (const auto& h : simHits_) { + if (h.trackId() != track_id) continue; + int pdgid = h.particleType(); + if (simMuOnly_ && std::abs(pdgid) != 13) continue; + // discard electron hits in the ME0 chambers + if (discardEleHits_ && std::abs(pdgid) == 11) continue; + + const ME0DetId& layer_id(h.detUnitId()); + detid_to_hits_[h.detUnitId()].push_back(h); + hits_.push_back(h); + chamber_to_hits_[layer_id.layerId().rawId()].push_back(h); + superChamber_to_hits_[layer_id.chamberId().rawId()].push_back(h); + } + } + + // find pads with hits + const auto& detids = detIds(); + for (const auto& d : detids) { + ME0DetId id(d); + const auto& hits = hitsInDetId(d); + const auto& roll = + dynamic_cast(geometry_)->etaPartition(id); + // int max_npads = roll->npads(); + set pads; + for (const auto& h : hits) { + const LocalPoint& lp = h.entryPoint(); + pads.insert(1 + static_cast(roll->padTopology().channel(lp))); + } + detids_to_pads_[d] = pads; + } +} + +std::set ME0SimHitMatcher::detIds() const { + std::set result; + for (const auto& p : detid_to_hits_) result.insert(p.first); + return result; +} + +std::set ME0SimHitMatcher::chamberIds() const { + std::set result; + for (const auto& p : chamber_to_hits_) result.insert(p.first); + return result; +} + +std::set ME0SimHitMatcher::superChamberIds() const { + std::set result; + for (const auto& p : superChamber_to_hits_) result.insert(p.first); + return result; +} + +const edm::PSimHitContainer& ME0SimHitMatcher::hitsInSuperChamber( + unsigned int detid) const { + if (superChamber_to_hits_.find(detid) == superChamber_to_hits_.end()) + return no_hits_; + return superChamber_to_hits_.at(detid); +} + +int ME0SimHitMatcher::nLayersWithHitsInSuperChamber(unsigned int detid) const { + set layers_with_hits; + const auto& hits = hitsInSuperChamber(detid); + for (const auto& h : hits) { + const ME0DetId& idd(h.detUnitId()); + layers_with_hits.insert(idd.layer()); + } + return layers_with_hits.size(); +} + +std::set ME0SimHitMatcher::superChamberIdsCoincidences( + int min_n_layers) const { + set result; + // loop over the super chamber Ids + for (const auto& p : superChamberIds()) { + if (nLayersWithHitsInSuperChamber(p) >= min_n_layers) { + result.insert(p); + } + } + return result; +} + +int ME0SimHitMatcher::nCoincidenceChambers(int min_n_layers) const { + return superChamberIdsCoincidences(min_n_layers).size(); +} + +float ME0SimHitMatcher::simHitsMeanStrip( + const edm::PSimHitContainer& sim_hits) const { + if (sim_hits.empty()) return -1.f; + + float sums = 0.f; + size_t n = 0; + for (const auto& h : sim_hits) { + const LocalPoint& lp = h.entryPoint(); + float s; + const auto& d = h.detUnitId(); + s = dynamic_cast(geometry_)->etaPartition(d)->strip(lp); + sums += s; + ++n; + } + if (n == 0) return -1.f; + return sums / n; +} + +std::set ME0SimHitMatcher::hitStripsInDetId(unsigned int detid, + int margin_n_strips) const { + set result; + const auto& simhits = hitsInDetId(detid); + ME0DetId id(detid); + int max_nstrips = + dynamic_cast(geometry_)->etaPartition(id)->nstrips(); + for (const auto& h : simhits) { + const LocalPoint& lp = h.entryPoint(); + int central_strip = + 1 + static_cast(dynamic_cast(geometry_) + ->etaPartition(id) + ->topology() + .channel(lp)); + int smin = central_strip - margin_n_strips; + smin = (smin > 0) ? smin : 1; + int smax = central_strip + margin_n_strips; + smax = (smax <= max_nstrips) ? smax : max_nstrips; + for (int ss = smin; ss <= smax; ++ss) result.insert(ss); + } + return result; +} + +std::set ME0SimHitMatcher::hitPadsInDetId(unsigned int detid) const { + set none; + if (detids_to_pads_.find(detid) == detids_to_pads_.end()) return none; + return detids_to_pads_.at(detid); +} + +std::set ME0SimHitMatcher::hitPartitions() const { + std::set result; + + const auto& detids = detIds(); + for (const auto& id : detids) { + ME0DetId idd(id); + result.insert(idd.roll()); + } + return result; +} + +int ME0SimHitMatcher::nPadsWithHits() const { + int result = 0; + const auto& pad_ids = detIds(); + for (const auto& id : pad_ids) { + result += hitPadsInDetId(id).size(); + } + return result; +} diff --git a/Validation/MuonHits/src/MuonHitHelper.cc b/Validation/MuonHits/src/MuonHitHelper.cc new file mode 100644 index 0000000000000..ec3281ac78927 --- /dev/null +++ b/Validation/MuonHits/src/MuonHitHelper.cc @@ -0,0 +1,165 @@ +#include "Validation/MuonHits/interface/MuonHitHelper.h" + +bool MuonHitHelper::isDT(unsigned int detId) { + return (DetId(detId)).det() == DetId::Muon && + (DetId(detId)).subdetId() == MuonSubdetId::DT; +} + +bool MuonHitHelper::isGEM(unsigned int detId) { + return (DetId(detId)).det() == DetId::Muon && + (DetId(detId)).subdetId() == MuonSubdetId::GEM; +} + +bool MuonHitHelper::isCSC(unsigned int detId) { + return (DetId(detId)).det() == DetId::Muon && + (DetId(detId)).subdetId() == MuonSubdetId::CSC; +} + +bool MuonHitHelper::isRPC(unsigned int detId) { + return (DetId(detId)).det() == DetId::Muon && + (DetId(detId)).subdetId() == MuonSubdetId::RPC; +} + +bool MuonHitHelper::isME0(unsigned int detId) { + return (DetId(detId)).det() == DetId::Muon && + (DetId(detId)).subdetId() == MuonSubdetId::ME0; +} + +int MuonHitHelper::chamber(const DetId& id) { + if (id.det() != DetId::Detector::Muon) return -99; + int chamberN = 0; + switch (id.subdetId()) { + case MuonSubdetId::GEM: + chamberN = GEMDetId(id).chamber(); + break; + case MuonSubdetId::RPC: + // works only for endcap!! + chamberN = RPCDetId(id).sector(); + break; + case MuonSubdetId::CSC: + chamberN = CSCDetId(id).chamber(); + break; + case MuonSubdetId::ME0: + chamberN = ME0DetId(id).chamber(); + break; + case MuonSubdetId::DT: + chamberN = DTChamberId(id).sector(); + break; + }; + return chamberN; +} + +// return MuonType for a particular DetId +int MuonHitHelper::toGEMType(int st, int ri) { + if (st == 1) { + if (ri == 1) return GEM_ME11; + } else if (st == 2) { + if (ri == 1) return GEM_ME21; + } + return GEM_ALL; +} + +int MuonHitHelper::toRPCType(int re, int st, int ri) { + // endcap + if (std::abs(re) == 1) { + if (st == 1) { + if (ri == 2) return RPC_ME12; + if (ri == 3) return RPC_ME13; + } else if (st == 2) { + if (ri == 2) return RPC_ME22; + if (ri == 3) return RPC_ME23; + } else if (st == 3) { + if (ri == 1) return RPC_ME31; + if (ri == 2) return RPC_ME32; + if (ri == 3) return RPC_ME33; + } else if (st == 4) { + if (ri == 1) return RPC_ME41; + if (ri == 2) return RPC_ME42; + if (ri == 3) return RPC_ME43; + } + } + // Barrel + else { + if (ri == -2) { + if (st == 1) return RPC_MB21n; + if (st == 2) return RPC_MB22n; + if (st == 3) return RPC_MB23n; + if (st == 4) return RPC_MB24n; + } else if (ri == -1) { + if (st == 1) return RPC_MB11n; + if (st == 2) return RPC_MB12n; + if (st == 3) return RPC_MB13n; + if (st == 4) return RPC_MB14n; + } else if (ri == 0) { + if (st == 1) return RPC_MB01; + if (st == 2) return RPC_MB02; + if (st == 3) return RPC_MB03; + if (st == 4) return RPC_MB04; + } else if (ri == 1) { + if (st == 1) return RPC_MB11p; + if (st == 2) return RPC_MB12p; + if (st == 3) return RPC_MB13p; + if (st == 4) return RPC_MB14p; + } else if (ri == 2) { + if (st == 1) return RPC_MB21p; + if (st == 2) return RPC_MB22p; + if (st == 3) return RPC_MB23p; + if (st == 4) return RPC_MB24p; + } + } + return RPC_ALL; +} + +int MuonHitHelper::toDTType(int wh, int st) { + if (wh == -2) { + if (st == 1) return DT_MB21n; + if (st == 2) return DT_MB22n; + if (st == 3) return DT_MB23n; + if (st == 4) return DT_MB24n; + } + if (wh == -1) { + if (st == 1) return DT_MB11n; + if (st == 2) return DT_MB12n; + if (st == 3) return DT_MB13n; + if (st == 4) return DT_MB14n; + } + if (wh == 0) { + if (st == 1) return DT_MB01; + if (st == 2) return DT_MB02; + if (st == 3) return DT_MB03; + if (st == 4) return DT_MB04; + } + if (wh == 1) { + if (st == 1) return DT_MB11p; + if (st == 2) return DT_MB12p; + if (st == 3) return DT_MB13p; + if (st == 4) return DT_MB14p; + } + if (wh == 2) { + if (st == 1) return DT_MB21p; + if (st == 2) return DT_MB22p; + if (st == 3) return DT_MB23p; + if (st == 4) return DT_MB24p; + } + return DT_ALL; +} + +int MuonHitHelper::toCSCType(int st, int ri) { + if (st == 1) { + if (ri == 0) return CSC_ME11; + if (ri == 1) return CSC_ME1b; + if (ri == 2) return CSC_ME12; + if (ri == 3) return CSC_ME13; + if (ri == 4) return CSC_ME1a; + } else if (st == 2) { + if (ri == 1) return CSC_ME21; + if (ri == 2) return CSC_ME22; + } else if (st == 3) { + if (ri == 1) return CSC_ME31; + if (ri == 2) return CSC_ME32; + } else if (st == 4) { + if (ri == 1) return CSC_ME41; + if (ri == 2) return CSC_ME42; + } + return CSC_ALL; +} diff --git a/Validation/MuonHits/src/MuonSimHitMatcher.cc b/Validation/MuonHits/src/MuonSimHitMatcher.cc new file mode 100644 index 0000000000000..cf26b358c51e1 --- /dev/null +++ b/Validation/MuonHits/src/MuonSimHitMatcher.cc @@ -0,0 +1,139 @@ +#include "Validation/MuonHits/interface/MuonSimHitMatcher.h" + +#include + +using namespace std; + +MuonSimHitMatcher::MuonSimHitMatcher(const edm::ParameterSet& ps, + edm::ConsumesCollector&& iC) { + const auto& simVertex = ps.getParameterSet("simVertex"); + const auto& simTrack = ps.getParameterSet("simTrack"); + verboseSimTrack_ = simTrack.getParameter("verbose"); + + simVertexInput_ = iC.consumes( + simVertex.getParameter("inputTag")); + simTrackInput_ = iC.consumes( + simTrack.getParameter("inputTag")); +} + +/// initialize the event +void MuonSimHitMatcher::init(const edm::Event& iEvent, + const edm::EventSetup& iSetup) { + hasGeometry_ = true; + + iEvent.getByToken(simTrackInput_, simTracksH_); + iEvent.getByToken(simVertexInput_, simVerticesH_); + iEvent.getByToken(simHitInput_, simHitsH_); +} + +/// do the matching +void MuonSimHitMatcher::match(const SimTrack& track, const SimVertex& vertex) { + simTracks_ = *simTracksH_.product(); + simVertices_ = *simVerticesH_.product(); + simHits_ = *simHitsH_.product(); + + // fill trkId2Index association: + int no = 0; + trkid_to_index_.clear(); + for (const auto& t : simTracks_) { + trkid_to_index_[t.trackId()] = no; + no++; + } + + track_ids_ = + getIdsOfSimTrackShower(track.trackId(), simTracks_, simVertices_); + if (verboseSimTrack_) { + edm::LogInfo("MuonSimHitMatcher") << "Printing track_ids" << std::endl; + for (const auto& id : track_ids_) edm::LogInfo("MuonSimHitMatcher") << "id: " << id << std::endl; + } +} + +std::vector MuonSimHitMatcher::getIdsOfSimTrackShower( + unsigned int initial_trk_id, const edm::SimTrackContainer& simTracks, + const edm::SimVertexContainer& simVertices) { + vector result; + result.push_back(initial_trk_id); + + if (!simMuOnly_) return result; + + for (const auto& t : simTracks_) { + SimTrack last_trk = t; + // if (std::abs(t.type()) != 13) continue; + bool is_child = false; + while (true) { + if (last_trk.noVertex()) break; + if (simVertices_[last_trk.vertIndex()].noParent()) break; + + unsigned parentId = simVertices_[last_trk.vertIndex()].parentIndex(); + if (parentId == initial_trk_id) { + is_child = true; + break; + } + + const auto& association = trkid_to_index_.find(parentId); + if (association == trkid_to_index_.end()) break; + + last_trk = simTracks_[association->second]; + } + if (is_child) { + result.push_back(t.trackId()); + } + } + return result; +} + +const edm::PSimHitContainer& MuonSimHitMatcher::simHits(int sub) const { + return hits_; +} + +const edm::PSimHitContainer& MuonSimHitMatcher::hitsInDetId( + unsigned int detid) const { + if (detid_to_hits_.find(detid) == detid_to_hits_.end()) return no_hits_; + return detid_to_hits_.at(detid); +} + +const edm::PSimHitContainer& MuonSimHitMatcher::hitsInChamber( + unsigned int detid) const { + if (chamber_to_hits_.find(detid) == chamber_to_hits_.end()) return no_hits_; + return chamber_to_hits_.at(detid); +} + +GlobalPoint MuonSimHitMatcher::simHitsMeanPosition( + const edm::PSimHitContainer& sim_hits) const { + if (sim_hits.empty()) return GlobalPoint(); // point "zero" + + float sumx, sumy, sumz; + sumx = sumy = sumz = 0.f; + size_t n = 0; + for (const auto& h : sim_hits) { + const LocalPoint& lp = h.entryPoint(); + const GlobalPoint& gp = + geometry_->idToDet(h.detUnitId())->surface().toGlobal(lp); + sumx += gp.x(); + sumy += gp.y(); + sumz += gp.z(); + ++n; + } + if (n == 0) return GlobalPoint(); + return GlobalPoint(sumx / n, sumy / n, sumz / n); +} + +GlobalVector MuonSimHitMatcher::simHitsMeanMomentum( + const edm::PSimHitContainer& sim_hits) const { + if (sim_hits.empty()) return GlobalVector(); // point "zero" + + float sumx, sumy, sumz; + sumx = sumy = sumz = 0.f; + size_t n = 0; + for (const auto& h : sim_hits) { + const LocalVector& lv = h.momentumAtEntry(); + const GlobalVector& gv = + geometry_->idToDet(h.detUnitId())->surface().toGlobal(lv); + sumx += gv.x(); + sumy += gv.y(); + sumz += gv.z(); + ++n; + } + if (n == 0) return GlobalVector(); + return GlobalVector(sumx / n, sumy / n, sumz / n); +} diff --git a/Validation/MuonHits/src/RPCSimHitMatcher.cc b/Validation/MuonHits/src/RPCSimHitMatcher.cc new file mode 100644 index 0000000000000..1da78f99c8c25 --- /dev/null +++ b/Validation/MuonHits/src/RPCSimHitMatcher.cc @@ -0,0 +1,161 @@ +#include "Validation/MuonHits/interface/RPCSimHitMatcher.h" +#include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h" + +using namespace std; + +RPCSimHitMatcher::RPCSimHitMatcher(const edm::ParameterSet& ps, + edm::ConsumesCollector&& iC) + : MuonSimHitMatcher(ps, std::move(iC)) { + simHitPSet_ = ps.getParameterSet("rpcSimHit"); + verbose_ = simHitPSet_.getParameter("verbose"); + simMuOnly_ = simHitPSet_.getParameter("simMuOnly"); + discardEleHits_ = simHitPSet_.getParameter("discardEleHits"); + + simHitInput_ = iC.consumes( + simHitPSet_.getParameter("inputTag")); +} + +/// initialize the event +void RPCSimHitMatcher::init(const edm::Event& iEvent, + const edm::EventSetup& iSetup) { + iSetup.get().get(rpc_geom_); + if (rpc_geom_.isValid()) { + geometry_ = &*rpc_geom_; + } else { + hasGeometry_ = false; + edm::LogWarning("RPCSimHitMatcher") + << "+++ Info: RPC geometry is unavailable. +++\n"; + } + MuonSimHitMatcher::init(iEvent, iSetup); +} + +/// do the matching +void RPCSimHitMatcher::match(const SimTrack& track, const SimVertex& vertex) { + // instantiates the track ids and simhits + MuonSimHitMatcher::match(track, vertex); + + if (hasGeometry_) { + matchSimHitsToSimTrack(); + + if (verbose_) { + edm::LogInfo("RPCSimHitMatcher") + << "nSimHits " << simHits_.size() << " nTrackIds " + << track_ids_.size() << endl; + edm::LogInfo("RPCSimHitMatcher") + << "detids RPC " << detIds().size() << endl; + + const auto& ch_ids = chamberIds(); + for (const auto& id : ch_ids) { + const auto& simhits = MuonSimHitMatcher::hitsInChamber(id); + const auto& simhits_gp = simHitsMeanPosition(simhits); + edm::LogInfo("RPCSimHitMatcher") + << "RPCDetId " << RPCDetId(id) << ": nHits " << simhits.size() + << " eta " << simhits_gp.eta() << " phi " << simhits_gp.phi() + << " nCh " << chamber_to_hits_[id].size() << endl; + const auto& strips = hitStripsInDetId(id); + edm::LogInfo("RPCSimHitMatcher") << "nStrips " << strips.size() << endl; + edm::LogInfo("RPCSimHitMatcher") << "strips : "; + for (const auto& p : strips) { + edm::LogInfo("RPCSimHitMatcher") << p; + } + } + } + } +} + +void RPCSimHitMatcher::matchSimHitsToSimTrack() { + for (const auto& track_id : track_ids_) { + for (const auto& h : simHits_) { + if (h.trackId() != track_id) continue; + int pdgid = h.particleType(); + if (simMuOnly_ && std::abs(pdgid) != 13) continue; + // discard electron hits in the RPC chambers + if (discardEleHits_ && pdgid == 11) continue; + + const RPCDetId& layer_id(h.detUnitId()); + detid_to_hits_[h.detUnitId()].push_back(h); + hits_.push_back(h); + chamber_to_hits_[layer_id.chamberId().rawId()].push_back(h); + } + } +} + +std::set RPCSimHitMatcher::detIds(int type) const { + std::set result; + for (const auto& p : detid_to_hits_) { + const auto& id = p.first; + if (type > 0) { + RPCDetId detId(id); + if (MuonHitHelper::toRPCType(detId.region(), detId.station(), + detId.ring()) != type) + continue; + } + result.insert(id); + } + return result; +} + +std::set RPCSimHitMatcher::chamberIds(int type) const { + std::set result; + for (const auto& p : chamber_to_hits_) { + const auto& id = p.first; + if (type > 0) { + RPCDetId detId(id); + if (MuonHitHelper::toRPCType(detId.region(), detId.station(), + detId.ring()) != type) + continue; + } + result.insert(id); + } + return result; +} + +bool RPCSimHitMatcher::hitStation(int st) const { + int nst = 0; + for (const auto& ddt : chamberIds(0)) { + const RPCDetId id(ddt); + if (id.station() != st) continue; + ++nst; + } + return nst; +} + +int RPCSimHitMatcher::nStations() const { + return (hitStation(1) + hitStation(2) + hitStation(3) + hitStation(4)); +} + +float RPCSimHitMatcher::simHitsMeanStrip( + const edm::PSimHitContainer& sim_hits) const { + if (sim_hits.empty()) return -1.f; + + float sums = 0.f; + size_t n = 0; + for (const auto& h : sim_hits) { + const LocalPoint& lp = h.entryPoint(); + const auto& d = h.detUnitId(); + sums += dynamic_cast(geometry_)->roll(d)->strip(lp); + ++n; + } + if (n == 0) return -1.f; + return sums / n; +} + +std::set RPCSimHitMatcher::hitStripsInDetId(unsigned int detid, + int margin_n_strips) const { + set result; + RPCDetId id(detid); + for (const auto& roll : + dynamic_cast(geometry_)->chamber(id)->rolls()) { + int max_nstrips = roll->nstrips(); + for (const auto& h : MuonSimHitMatcher::hitsInDetId(roll->id().rawId())) { + const LocalPoint& lp = h.entryPoint(); + int central_strip = static_cast(roll->topology().channel(lp)); + int smin = central_strip - margin_n_strips; + smin = (smin > 0) ? smin : 1; + int smax = central_strip + margin_n_strips; + smax = (smax <= max_nstrips) ? smax : max_nstrips; + for (int ss = smin; ss <= smax; ++ss) result.insert(ss); + } + } + return result; +}