Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
57 commits
Select commit Hold shift + click to select a range
59538f0
update for 2018
jshlee Jul 17, 2018
33713ee
merge
yeckang Jul 19, 2018
2610bf8
add lumiScalers
yeckang Jul 19, 2018
826693f
add instant luminosity info
yeckang Jul 19, 2018
bc92587
add digicollection
yeckang Jul 23, 2018
6beba59
Merge pull request #1 from CPLUOS/master
yeckang Jul 23, 2018
75ce542
add nBxing info to output
yeckang Jul 24, 2018
87d873b
update
yeckang Aug 3, 2018
eebaf6f
checkpoint
yeckang Aug 16, 2018
e0b9516
add scripts for background plots
yeckang Aug 24, 2018
3450907
add cscrechits
yeckang Aug 31, 2018
98850c6
update latency
yeckang Sep 6, 2018
f1a85fa
update missing part
yeckang Sep 6, 2018
49e4161
add some histogram for drawing background plots
yeckang Sep 9, 2018
0434c17
checkpoint
yeckang Sep 10, 2018
cd88a2f
add muon branch and csc segment
yeckang Sep 11, 2018
a06024a
fix csc segments
yeckang Sep 11, 2018
0622afe
remove personal script
yeckang Sep 11, 2018
9b589f7
rename bkg analysis to avoid confilict
yeckang Sep 28, 2018
1d9c57c
merged from CPLUOS 28.Sep.2018
yeckang Sep 28, 2018
2634dd0
change name for consistency
yeckang Sep 28, 2018
4a2a3f0
clean up
yeckang Sep 28, 2018
dd3fe37
add more position information
yeckang Oct 9, 2018
d2cf249
Add GEMDebug and STASliceTestAnalysis
watson-ij Oct 11, 2018
5354987
Add scripts to run Gem Simulation, add output to STASliceTest
watson-ij Oct 11, 2018
f795490
Fix GemSim scripts
watson-ij Oct 11, 2018
bf430e3
Update STASliceTest for Jasons branch, add runners
watson-ij Oct 11, 2018
aa14c2b
Put runners in correct location
watson-ij Oct 11, 2018
c5a56e5
Merge pull request #3 from watson-ij/master
yeckang Oct 12, 2018
28b17b2
fix name
yeckang Oct 12, 2018
120676f
update bkg hit rate alnalyse
yeckang Dec 9, 2018
6743c6a
add scripts for changing of HGCal geometry
yeckang Jan 22, 2019
8ec86ac
Merge pull request #1 from yeckang/me0_StainlessSteel
jang00777 Jan 22, 2019
e0d88b3
Save b_etaPartition
jang00777 Jan 22, 2019
566f94c
Get me0Segment (not complete)
jang00777 Jan 22, 2019
75da58d
fix bug in segment loop
yeckang Jan 23, 2019
1c1bafa
add script for tenMu_modified samples
yeckang Jan 23, 2019
6085e34
Add CSC and GEM
jang00777 Jan 31, 2019
3eaa7e3
Save all stations and rings and add DT and RPC (cleaning up the code …
jang00777 Feb 3, 2019
21bf3c4
Save all of rings and stations and add DT and RPC
jang00777 Feb 3, 2019
a9e28b0
Add eta for digi and wire for DT
jang00777 Feb 5, 2019
5ccc3b2
Fix bug about variables(eta, wire) in seg, digi, rechit loop
jang00777 Feb 7, 2019
1818408
Add reco muons
jang00777 Feb 8, 2019
a241d3c
Add RPC region (I forgot that ... )
jang00777 Feb 11, 2019
e2998b4
Save hit counts and detector areas for ME0/GEM/CSC/RPC
jang00777 Feb 12, 2019
b2916cb
Add additional variables (phi for dets and isGlobal, isTracker and no…
jang00777 Feb 19, 2019
c98e514
Except slice test things and split HGCalSimTest into MuonDetHitAnalys…
jang00777 Feb 25, 2019
d0d6c9a
back
jang00777 Feb 25, 2019
40d08fb
Merge branch 'master' of https://github.com/CPLUOS/MuonPerformance in…
jang00777 Feb 26, 2019
a65d8eb
Add MuonSimAnalyser
jang00777 Feb 26, 2019
4069472
Fix bug for running HGCalSimTest
jang00777 Feb 28, 2019
7f8ac7a
Fix bug for running HGCalSimTest (+1)
jang00777 Feb 28, 2019
f7a72f8
Add ME0 simhit in MuonSimAnalyser
jang00777 Mar 4, 2019
58e0d25
Add isLooseME0 in MuonTrackAnalyser
jang00777 Mar 13, 2019
9491f90
Add analysis code for several detectors (#8)
jang00777 Mar 13, 2019
6b0d32b
Merge branch 'master' into Woojin
jang00777 Mar 14, 2019
4f5e8b3
Add function to extract track and vertex and get mom pdgId from simhit
jang00777 Mar 26, 2019
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
124 changes: 101 additions & 23 deletions MuonAnalyser/plugins/MuonSimAnalyser.cc
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,7 @@
#include "TLorentzVector.h"
#include "TTree.h"
#include "TDatabasePDG.h"
#include <tuple>

using namespace std;
using namespace edm;
Expand All @@ -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<SimTrack, unsigned int, unsigned int, int> 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_;
Expand All @@ -99,18 +106,21 @@ class MuonSimAnalyser : public edm::EDAnalyzer {

edm::Service<TFileService> 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;

};

Expand All @@ -132,9 +142,11 @@ MuonSimAnalyser::MuonSimAnalyser(const edm::ParameterSet& iConfig)
/*ME0*/
t_ME0_simhit = fs->make<TTree>("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");
Expand All @@ -143,9 +155,11 @@ MuonSimAnalyser::MuonSimAnalyser(const edm::ParameterSet& iConfig)
/*GEM*/
t_GEM_simhit = fs->make<TTree>("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");
Expand All @@ -162,10 +176,15 @@ MuonSimAnalyser::~MuonSimAnalyser()
void
MuonSimAnalyser::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
/* ME0 Geometry */
edm::ESHandle<ME0Geometry> hME0Geom;
iSetup.get<MuonGeometryRecord>().get(hME0Geom);
const ME0Geometry* ME0Geometry_ = &*hME0Geom;

/* GEM Geometry */
//edm::ESHandle<GEMGeometry> hGEMGeom;
//iSetup.get<MuonGeometryRecord>().get(hGEMGeom);
//const GEMGeometry* GEMGeometry_ = &*hGEMGeom;
edm::ESHandle<GEMGeometry> hGEMGeom;
iSetup.get<MuonGeometryRecord>().get(hGEMGeom);
const GEMGeometry* GEMGeometry_ = &*hGEMGeom;

edm::Handle<edm::PSimHitContainer> ME0SimHits;
edm::Handle<edm::PSimHitContainer> GEMSimHits;
Expand All @@ -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();
Expand All @@ -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();
Expand All @@ -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<SimTrack, unsigned int, unsigned int, int> 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
Expand Down
70 changes: 66 additions & 4 deletions MuonAnalyser/plugins/MuonTrackAnalyser.cc
Original file line number Diff line number Diff line change
Expand Up @@ -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<ME0SegmentCollection> me0Segments_;
Expand All @@ -132,6 +134,8 @@ class MuonTrackAnalyser : public edm::EDAnalyzer {

edm::Service<TFileService> fs;

const ME0Geometry* ME0Geometry_;

TTree *t_event;
int b_run, b_lumi, b_latency;
int b_nME0Segments, b_nCSCSegments, b_nGEMSegments, b_nDT4DSegments;
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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");
Expand All @@ -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");
Expand All @@ -267,7 +273,8 @@ MuonTrackAnalyser::analyze(const edm::Event& iEvent, const edm::EventSetup& iSet
{
edm::ESHandle<ME0Geometry> hME0Geom;
iSetup.get<MuonGeometryRecord>().get(hME0Geom);
const ME0Geometry* ME0Geometry_ = &*hME0Geom;

ME0Geometry_ = &*hME0Geom;

/* ME0 Geometry */
edm::Handle<ME0SegmentCollection> me0Segments;
Expand Down Expand Up @@ -418,6 +425,12 @@ MuonTrackAnalyser::analyze(const edm::Event& iEvent, const edm::EventSetup& iSet
initMuonValue();
edm::RefToBase<reco::Muon> 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);
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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<reco::MuonChamberMatch>& chambers = muon.matches();
for( std::vector<reco::MuonChamberMatch>::const_iterator chamber = chambers.begin(); chamber != chambers.end(); ++chamber ){

if (chamber->detector() == 5){

for ( std::vector<reco::MuonSegmentMatch>::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);