diff --git a/MuonAnalyser/plugins/MuonSimAnalyser.cc b/MuonAnalyser/plugins/MuonSimAnalyser.cc index 5e16ebeb870..96639cb03c6 100644 --- a/MuonAnalyser/plugins/MuonSimAnalyser.cc +++ b/MuonAnalyser/plugins/MuonSimAnalyser.cc @@ -71,6 +71,7 @@ #include "TLorentzVector.h" #include "TTree.h" #include "TDatabasePDG.h" +#include using namespace std; using namespace edm; @@ -90,6 +91,12 @@ class MuonSimAnalyser : public edm::EDAnalyzer { void initValue(); + enum detType {DET_NONE = 0, DET_ME0, DET_GEM, DET_CSC, DET_DT, DET_RPC}; + + bool isSimTrackGood(const SimTrack & , int = -1, detType = DET_NONE); + std::tuple matchSimHitToSimTrack(const PSimHit & sh, const edm::SimTrackContainer & sim_track, int pdg=-1, detType = DET_NONE); + int getMotherPdgId(const SimTrack & trk, const edm::PSimHitContainer & sim_hit, const edm::SimVertexContainer & sim_vertex); + // ----------member data --------------------------- edm::ParameterSet cfg_; edm::EDGetToken ME0SimHitsToken_; @@ -99,18 +106,21 @@ class MuonSimAnalyser : public edm::EDAnalyzer { edm::Service fs; + float GEM_minEta_ = 1.5; float GEM_maxEta_ = 2.5; + float ME0_minEta_ = 1.9; float ME0_maxEta_ = 3.; + TTree *t_event; int b_nME0SimHits, b_nGEMSimHits; /* ME */ TTree *t_ME0_simhit; - int b_ME0_SimHit_region, b_ME0_SimHit_chamber, b_ME0_SimHit_layer, b_ME0_SimHit_etaPartition, b_ME0_SimHit_pdgId; - float b_ME0_SimHit_pt, b_ME0_SimHit_eta, b_ME0_SimHit_phi; + int b_ME0_SimHit_region, b_ME0_SimHit_chamber, b_ME0_SimHit_layer, b_ME0_SimHit_etaPartition, b_ME0_SimHit_pdgId, b_ME0_SimHit_momPdgId; + float b_ME0_SimHit_pt, b_ME0_SimHit_eta, b_ME0_SimHit_phi, b_ME0_SimHit_E; /* GEM */ TTree *t_GEM_simhit; - int b_GEM_SimHit_region, b_GEM_SimHit_station, b_GEM_SimHit_ring, b_GEM_SimHit_chamber, b_GEM_SimHit_layer, b_GEM_SimHit_etaPartition, b_GEM_SimHit_pdgId; - float b_GEM_SimHit_pt, b_GEM_SimHit_eta, b_GEM_SimHit_phi; + int b_GEM_SimHit_region, b_GEM_SimHit_station, b_GEM_SimHit_ring, b_GEM_SimHit_chamber, b_GEM_SimHit_layer, b_GEM_SimHit_etaPartition, b_GEM_SimHit_pdgId, b_GEM_SimHit_momPdgId; + float b_GEM_SimHit_pt, b_GEM_SimHit_eta, b_GEM_SimHit_phi, b_GEM_SimHit_E; }; @@ -132,9 +142,11 @@ MuonSimAnalyser::MuonSimAnalyser(const edm::ParameterSet& iConfig) /*ME0*/ t_ME0_simhit = fs->make("ME0_SimHit", "ME0_SimHit"); t_ME0_simhit->Branch("SimHit_pdgId", &b_ME0_SimHit_pdgId, "SimHit_pdgId/I"); + t_ME0_simhit->Branch("SimHit_momPdgId", &b_ME0_SimHit_momPdgId, "SimHit_momPdgId/I"); t_ME0_simhit->Branch("SimHit_pt", &b_ME0_SimHit_pt, "SimHit_pt/F"); t_ME0_simhit->Branch("SimHit_eta", &b_ME0_SimHit_eta, "SimHit_eta/F"); t_ME0_simhit->Branch("SimHit_phi", &b_ME0_SimHit_phi, "SimHit_phi/F"); + t_ME0_simhit->Branch("SimHit_E", &b_ME0_SimHit_E, "SimHit_E/F"); t_ME0_simhit->Branch("SimHit_region", &b_ME0_SimHit_region, "SimHit_region/I"); t_ME0_simhit->Branch("SimHit_chamber", &b_ME0_SimHit_chamber, "SimHit_chaber/I"); t_ME0_simhit->Branch("SimHit_layer", &b_ME0_SimHit_layer, "SimHit_layer/I"); @@ -143,9 +155,11 @@ MuonSimAnalyser::MuonSimAnalyser(const edm::ParameterSet& iConfig) /*GEM*/ t_GEM_simhit = fs->make("GEM_SimHit", "GEM_SimHit"); t_GEM_simhit->Branch("SimHit_pdgId", &b_GEM_SimHit_pdgId, "SimHit_pdgId/I"); + t_GEM_simhit->Branch("SimHit_momPdgId", &b_GEM_SimHit_momPdgId, "SimHit_momPdgId/I"); t_GEM_simhit->Branch("SimHit_pt", &b_GEM_SimHit_pt, "SimHit_pt/F"); t_GEM_simhit->Branch("SimHit_eta", &b_GEM_SimHit_eta, "SimHit_eta/F"); t_GEM_simhit->Branch("SimHit_phi", &b_GEM_SimHit_phi, "SimHit_phi/F"); + t_GEM_simhit->Branch("SimHit_E", &b_GEM_SimHit_E, "SimHit_E/F"); t_GEM_simhit->Branch("SimHit_region", &b_GEM_SimHit_region, "SimHit_region/I"); t_GEM_simhit->Branch("SimHit_station", &b_GEM_SimHit_station, "SimHit_station/I"); t_GEM_simhit->Branch("SimHit_ring", &b_GEM_SimHit_ring, "SimHit_ring/I"); @@ -162,10 +176,15 @@ MuonSimAnalyser::~MuonSimAnalyser() void MuonSimAnalyser::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) { + /* ME0 Geometry */ + edm::ESHandle hME0Geom; + iSetup.get().get(hME0Geom); + const ME0Geometry* ME0Geometry_ = &*hME0Geom; + /* GEM Geometry */ - //edm::ESHandle hGEMGeom; - //iSetup.get().get(hGEMGeom); - //const GEMGeometry* GEMGeometry_ = &*hGEMGeom; + edm::ESHandle hGEMGeom; + iSetup.get().get(hGEMGeom); + const GEMGeometry* GEMGeometry_ = &*hGEMGeom; edm::Handle ME0SimHits; edm::Handle GEMSimHits; @@ -179,24 +198,36 @@ MuonSimAnalyser::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup initValue(); + const edm::SimTrackContainer & sim_tracks = *simTracks.product(); + const edm::SimVertexContainer & sim_vertices = *simVertices.product(); + + /* ME0 */ if (ME0SimHits.isValid()) { const edm::PSimHitContainer & ME0_sim_hits = *ME0SimHits.product(); for (auto& sh : ME0_sim_hits){ - ME0DetId det_id = sh.detUnitId(); - auto vec = sh.momentumAtEntry(); + ME0DetId det_id(sh.detUnitId()); + auto vec = sh.momentumAtEntry(); + auto shLp = sh.localPosition(); + const BoundPlane & surface = ME0Geometry_->idToDet(det_id)->surface(); + auto shGp = surface.toGlobal(shLp); + b_ME0_SimHit_pdgId = sh.particleType(); + auto sim_track_tuple = matchSimHitToSimTrack(sh, sim_tracks); + b_ME0_SimHit_momPdgId = getMotherPdgId(std::get<0>(sim_track_tuple), ME0_sim_hits, sim_vertices); + std::cout << " >>>>>>> TUPLE : [ " << std::get<0>(sim_track_tuple) << " ] / [ " << std::get<1>(sim_track_tuple) << " ] / [ " << std::get<2>(sim_track_tuple) << " ] / [ " << std::get<3>(sim_track_tuple) << " ] " << std::endl; TDatabasePDG *pdg = TDatabasePDG::Instance(); if (auto particle_pdg = pdg->GetParticle(b_ME0_SimHit_pdgId)) { - std::cout << particle_pdg->GetName() << " " << b_ME0_SimHit_pdgId << std::endl; + std::cout << particle_pdg->GetName() << " " << b_ME0_SimHit_pdgId << " mom pdg : " << b_ME0_SimHit_momPdgId << std::endl; } else { std::cout << "Cannot understand!!! " << b_ME0_SimHit_pdgId << std::endl; } b_ME0_SimHit_pt = vec.perp(); - b_ME0_SimHit_eta = vec.eta(); - b_ME0_SimHit_phi = sh.phiAtEntry(); + b_ME0_SimHit_eta = shGp.eta(); + b_ME0_SimHit_phi = shGp.phi(); + b_ME0_SimHit_E = sh.pabs(); b_ME0_SimHit_region = det_id.region(); b_ME0_SimHit_chamber = det_id.chamber(); @@ -213,21 +244,29 @@ MuonSimAnalyser::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup if (GEMSimHits.isValid()) { const edm::PSimHitContainer & GEM_sim_hits = *GEMSimHits.product(); for (auto& sh : GEM_sim_hits){ - GEMDetId det_id = sh.detUnitId(); - auto vec = sh.momentumAtEntry(); + GEMDetId det_id(sh.detUnitId()); + auto vec = sh.momentumAtEntry(); + auto shLp = sh.localPosition(); + const BoundPlane & surface = GEMGeometry_->idToDet(det_id)->surface(); + auto shGp = surface.toGlobal(shLp); + b_GEM_SimHit_pdgId = sh.particleType(); - + auto sim_track_tuple = matchSimHitToSimTrack(sh, sim_tracks); + b_GEM_SimHit_momPdgId = getMotherPdgId(std::get<0>(sim_track_tuple), GEM_sim_hits, sim_vertices); + std::cout << " >>>>>>> TUPLE : [ " << std::get<0>(sim_track_tuple) << " ] / [ " << std::get<1>(sim_track_tuple) << " ] / [ " << std::get<2>(sim_track_tuple) << " ] / [ " << std::get<3>(sim_track_tuple) << " ] " << std::endl; + TDatabasePDG *pdg = TDatabasePDG::Instance(); if (auto particle_pdg = pdg->GetParticle(b_GEM_SimHit_pdgId)) { - std::cout << particle_pdg->GetName() << " " << b_GEM_SimHit_pdgId << std::endl; + std::cout << particle_pdg->GetName() << " " << b_GEM_SimHit_pdgId << " mom pdg : " << b_GEM_SimHit_momPdgId << std::endl; } else { std::cout << "Cannot understand!!! " << b_GEM_SimHit_pdgId << std::endl; } b_GEM_SimHit_pt = vec.perp(); - b_GEM_SimHit_eta = vec.eta(); - b_GEM_SimHit_phi = sh.phiAtEntry(); - + b_GEM_SimHit_eta = shGp.eta(); + b_GEM_SimHit_phi = shGp.phi(); + b_GEM_SimHit_E = sh.pabs(); + b_GEM_SimHit_region = det_id.region(); b_GEM_SimHit_station = det_id.station(); b_GEM_SimHit_ring = det_id.ring(); @@ -253,12 +292,51 @@ void MuonSimAnalyser::initValue() { b_nME0SimHits = -1; b_nGEMSimHits = -1; /*ME0 */ - b_ME0_SimHit_region = -9; b_ME0_SimHit_chamber = -9; b_ME0_SimHit_layer = -9; b_ME0_SimHit_etaPartition = -9; b_ME0_SimHit_pdgId = -99; - b_ME0_SimHit_pt = -9; b_ME0_SimHit_eta = -9; b_ME0_SimHit_phi = -9; + b_ME0_SimHit_region = -9; b_ME0_SimHit_chamber = -9; b_ME0_SimHit_layer = -9; b_ME0_SimHit_etaPartition = -9; b_ME0_SimHit_pdgId = -99; b_ME0_SimHit_momPdgId = -99; + b_ME0_SimHit_pt = -9; b_ME0_SimHit_eta = -9; b_ME0_SimHit_phi = -9; b_ME0_SimHit_E = -9; /*GEM */ - b_GEM_SimHit_region = -9; b_GEM_SimHit_station = -9; b_GEM_SimHit_ring = -9; b_GEM_SimHit_chamber = -9; b_GEM_SimHit_layer = -9; b_GEM_SimHit_etaPartition = -9; b_GEM_SimHit_pdgId = -99; - b_GEM_SimHit_pt = -9; b_GEM_SimHit_eta = -9; b_GEM_SimHit_phi = -9; + b_GEM_SimHit_region = -9; b_GEM_SimHit_station = -9; b_GEM_SimHit_ring = -9; b_GEM_SimHit_chamber = -9; b_GEM_SimHit_layer = -9; b_GEM_SimHit_etaPartition = -9; b_GEM_SimHit_pdgId = -99; b_GEM_SimHit_momPdgId = -99; + b_GEM_SimHit_pt = -9; b_GEM_SimHit_eta = -9; b_GEM_SimHit_phi = -9; b_GEM_SimHit_E = -9; +} + +bool MuonSimAnalyser::isSimTrackGood(const SimTrack &t, int pdg, detType detName) { + // SimTrack selection + if (t.noVertex()) return false; + //if (t.noGenpart()) return false; + if (pdg != -1 ) { + if (t.type() != pdg) return false; // only interested in pdg matched particle + } + //if (t.momentum().pt() < minPt_ ) return false; + const float eta(std::abs(t.momentum().eta())); + if (detName == DET_GEM) { + if (eta > GEM_maxEta_ || eta < GEM_minEta_ ) return false; // no GEMs could be in such eta + } else if (detName == DET_ME0) { + if (eta > ME0_maxEta_ || eta < ME0_minEta_ ) return false; // no ME0 could be in such eta + } + return true; +} + +std::tuple MuonSimAnalyser::matchSimHitToSimTrack(const PSimHit & sh, const edm::SimTrackContainer & sim_track, int pdg, detType detName) { + int pdgId = sh.particleType(); + for (auto & trk : sim_track) { + if (!isSimTrackGood(trk, pdg, detName)) continue; + if (sh.trackId() != trk.trackId()) continue; + return std::make_tuple(trk, trk.vertIndex(), trk.genpartIndex(), pdgId); + } + return std::make_tuple(SimTrack(), -1, -1, pdgId); +} + +int MuonSimAnalyser::getMotherPdgId(const SimTrack & trk, const edm::PSimHitContainer & sim_hit, const edm::SimVertexContainer & sim_vertex){ + unsigned int p_vtx_id = 0; + for (auto & vtx : sim_vertex) { + if (vtx.noParent()) continue; + if((int)vtx.vertexId() == trk.vertIndex()) p_vtx_id = vtx.parentIndex(); + } + for (auto & sh : sim_hit) { + if (sh.trackId() == p_vtx_id) return sh.particleType(); + } + return -999; } //define this as a plug-in diff --git a/MuonAnalyser/plugins/MuonTrackAnalyser.cc b/MuonAnalyser/plugins/MuonTrackAnalyser.cc index a3bc90aaefc..e4dd2973c89 100644 --- a/MuonAnalyser/plugins/MuonTrackAnalyser.cc +++ b/MuonAnalyser/plugins/MuonTrackAnalyser.cc @@ -112,6 +112,8 @@ class MuonTrackAnalyser : public edm::EDAnalyzer { void initMuonValue(); void initValue(); + bool isME0MuonSelNew(reco::Muon muon, double dEtaCut, double dPhiCut, double dPhiBendCut); + // ----------member data --------------------------- edm::EDGetTokenT me0Segments_; @@ -132,6 +134,8 @@ class MuonTrackAnalyser : public edm::EDAnalyzer { edm::Service fs; + const ME0Geometry* ME0Geometry_; + TTree *t_event; int b_run, b_lumi, b_latency; int b_nME0Segments, b_nCSCSegments, b_nGEMSegments, b_nDT4DSegments; @@ -159,11 +163,11 @@ class MuonTrackAnalyser : public edm::EDAnalyzer { /* Muon */ TTree *t_genMuon; float b_genMuon_pt, b_genMuon_eta, b_genMuon_phi; - int b_genMuon_isTight, b_genMuon_isMedium, b_genMuon_isLoose, b_genMuon_isTracker, b_genMuon_isGlobal, b_genMuon_isME0, b_genMuon_nMatchedStationLayer; + int b_genMuon_isTight, b_genMuon_isMedium, b_genMuon_isLoose, b_genMuon_isTracker, b_genMuon_isGlobal, b_genMuon_isME0, b_genMuon_isLooseME0, b_genMuon_nMatchedStationLayer; TTree *t_Muon; float b_muon_pt, b_muon_eta, b_muon_phi; - int b_muon_isTight, b_muon_isMedium, b_muon_isLoose, b_muon_isTracker, b_muon_isGlobal, b_muon_isME0, b_muon_nMatchedStationLayer; + int b_muon_isTight, b_muon_isMedium, b_muon_isLoose, b_muon_isTracker, b_muon_isGlobal, b_muon_isME0, b_muon_isLooseME0, b_muon_nMatchedStationLayer; }; MuonTrackAnalyser::MuonTrackAnalyser(const edm::ParameterSet& iConfig) @@ -240,6 +244,7 @@ MuonTrackAnalyser::MuonTrackAnalyser(const edm::ParameterSet& iConfig) t_genMuon->Branch("isTracker", &b_genMuon_isTracker, "isTracker/I"); t_genMuon->Branch("isGlobal", &b_genMuon_isGlobal, "isGlobal/I"); t_genMuon->Branch("isME0", &b_genMuon_isME0, "isME0/I"); + t_genMuon->Branch("isLooseME0", &b_genMuon_isLooseME0, "isLooseME0/I"); t_genMuon->Branch("isTight", &b_genMuon_isTight, "isTight/I"); t_genMuon->Branch("isMedium", &b_genMuon_isMedium, "isMedium/I"); t_genMuon->Branch("isLoose", &b_genMuon_isLoose, "isLoose/I"); @@ -252,6 +257,7 @@ MuonTrackAnalyser::MuonTrackAnalyser(const edm::ParameterSet& iConfig) t_Muon->Branch("isTracker", &b_muon_isTracker, "isTracker/I"); t_Muon->Branch("isGlobal", &b_muon_isGlobal, "isGlobal/I"); t_Muon->Branch("isME0", &b_muon_isME0, "isME0/I"); + t_Muon->Branch("isLooseME0", &b_muon_isLooseME0, "isLooseME0/I"); t_Muon->Branch("isTight", &b_muon_isTight, "isTight/I"); t_Muon->Branch("isMedium", &b_muon_isMedium, "isMedium/I"); t_Muon->Branch("isLoose", &b_muon_isLoose, "isLoose/I"); @@ -267,7 +273,8 @@ MuonTrackAnalyser::analyze(const edm::Event& iEvent, const edm::EventSetup& iSet { edm::ESHandle hME0Geom; iSetup.get().get(hME0Geom); - const ME0Geometry* ME0Geometry_ = &*hME0Geom; + + ME0Geometry_ = &*hME0Geom; /* ME0 Geometry */ edm::Handle me0Segments; @@ -418,6 +425,12 @@ MuonTrackAnalyser::analyze(const edm::Event& iEvent, const edm::EventSetup& iSet initMuonValue(); edm::RefToBase muRef = muonHandle->refAt(i); const reco::Muon* mu = muRef.get(); + + double mom = mu->p(); + double dPhiCut_ = std::min(std::max(1.2/mom,1.2/100),0.056); + double dPhiBendCut_ = std::min(std::max(0.2/mom,0.2/100),0.0096); + b_muon_isLooseME0 = isME0MuonSelNew(*mu, 0.077, dPhiCut_, dPhiBendCut_); + b_muon_isTight = muon::isTightMuon(*mu, pv0); b_muon_isMedium = muon::isMediumMuon(*mu); b_muon_isLoose = muon::isLooseMuon(*mu); @@ -445,7 +458,7 @@ void MuonTrackAnalyser::endRun(Run const&, EventSetup const&){} void MuonTrackAnalyser::initMuonValue() { b_genMuon_pt = -9; b_genMuon_eta = -9; b_genMuon_phi = -9; - b_muon_isTight = -1; b_muon_isMedium = -1; b_muon_isLoose = -1; b_muon_isTracker = -1; b_muon_isGlobal = -1; b_muon_isME0 = -1; + b_muon_isTight = -1; b_muon_isMedium = -1; b_muon_isLoose = -1; b_muon_isTracker = -1; b_muon_isGlobal = -1; b_muon_isME0 = -1; b_muon_isLooseME0 = -1; b_muon_pt = -9; b_muon_eta = -9; b_muon_phi = -9; b_muon_nMatchedStationLayer = -9; @@ -474,5 +487,54 @@ void MuonTrackAnalyser::initValue() { b_DT_4DSeg_hasZed = -1; } +bool MuonTrackAnalyser::isME0MuonSelNew(reco::Muon muon, double dEtaCut, double dPhiCut, double dPhiBendCut) +{ + bool result = false; + bool isME0 = muon.isME0Muon(); + + if(isME0){ + + double deltaEta = 999; + double deltaPhi = 999; + double deltaPhiBend = 999; + + const std::vector& chambers = muon.matches(); + for( std::vector::const_iterator chamber = chambers.begin(); chamber != chambers.end(); ++chamber ){ + + if (chamber->detector() == 5){ + + for ( std::vector::const_iterator segment = chamber->me0Matches.begin(); segment != chamber->me0Matches.end(); ++segment ){ + + LocalPoint trk_loc_coord(chamber->x, chamber->y, 0); + LocalPoint seg_loc_coord(segment->x, segment->y, 0); + LocalVector trk_loc_vec(chamber->dXdZ, chamber->dYdZ, 1); + LocalVector seg_loc_vec(segment->dXdZ, segment->dYdZ, 1); + + const ME0Chamber * me0chamber = ME0Geometry_->chamber(chamber->id); + + GlobalPoint trk_glb_coord = me0chamber->toGlobal(trk_loc_coord); + GlobalPoint seg_glb_coord = me0chamber->toGlobal(seg_loc_coord); + + //double segDPhi = segment->me0SegmentRef->deltaPhi(); + // need to check if this works + double segDPhi = me0chamber->computeDeltaPhi(seg_loc_coord, seg_loc_vec); + double trackDPhi = me0chamber->computeDeltaPhi(trk_loc_coord, trk_loc_vec); + + deltaEta = std::abs(trk_glb_coord.eta() - seg_glb_coord.eta() ); + deltaPhi = std::abs(trk_glb_coord.phi() - seg_glb_coord.phi() ); + deltaPhiBend = std::abs(segDPhi - trackDPhi); + + if (deltaEta < dEtaCut && deltaPhi < dPhiCut && deltaPhiBend < dPhiBendCut) result = true; + + } + } + } + + } + + return result; + +} + //define this as a plug-in DEFINE_FWK_MODULE(MuonTrackAnalyser);