diff --git a/MuonAnalyser/plugins/BuildFile.xml b/MuonAnalyser/plugins/BuildFile.xml index b5754384aa2..555883c5465 100644 --- a/MuonAnalyser/plugins/BuildFile.xml +++ b/MuonAnalyser/plugins/BuildFile.xml @@ -9,16 +9,21 @@ + + + + + diff --git a/MuonAnalyser/plugins/GEMDebug.cc b/MuonAnalyser/plugins/GEMDebug.cc new file mode 100644 index 00000000000..6e869292072 --- /dev/null +++ b/MuonAnalyser/plugins/GEMDebug.cc @@ -0,0 +1,259 @@ +// cd ../../.. && source /cvmfs/cms.cern.ch/cmsset_default.sh && eval `scramv1 runtime -sh` && eval `scramv1 runtime -sh` && scram b -j 10 + +// system include files +#include +#include +#include +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" + +#include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h" +#include "MagneticField/Engine/interface/MagneticField.h" + +#include "DataFormats/Scalers/interface/LumiScalers.h" +#include "DataFormats/MuonData/interface/MuonDigiCollection.h" +#include "DataFormats/GEMDigi/interface/GEMAMCStatusDigi.h" +#include "DataFormats/GEMDigi/interface/GEMVfatStatusDigi.h" +#include "DataFormats/GEMDigi/interface/GEMGEBStatusDigi.h" + +#include "DataFormats/PatCandidates/interface/Muon.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/GEMRecHit/interface/GEMRecHitCollection.h" +#include "DataFormats/MuonDetId/interface/GEMDetId.h" +#include "Geometry/GEMGeometry/interface/GEMGeometry.h" +#include "Geometry/GEMGeometry/interface/GEMEtaPartition.h" +#include "Geometry/GEMGeometry/interface/GEMEtaPartitionSpecs.h" +#include "Geometry/Records/interface/MuonGeometryRecord.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" +#include "Geometry/CommonTopologies/interface/StripTopology.h" + +#include "Geometry/Records/interface/MuonGeometryRecord.h" + +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/Run.h" + +#include "TH1D.h" +#include "TH2D.h" +#include "TString.h" +#include "TGraphAsymmErrors.h" +#include "TLorentzVector.h" +#include "TTree.h" + +#include "TVector2.h" + + +using namespace std; +using namespace edm; + +class GEMDebug : public edm::EDAnalyzer { +public: + explicit GEMDebug(const edm::ParameterSet&); + ~GEMDebug(); + +private: + virtual void analyze(const edm::Event&, const edm::EventSetup&); + virtual void beginJob() override; + virtual void endJob() override; + + virtual void beginRun(Run const&, EventSetup const&) override; + virtual void endRun(Run const&, EventSetup const&) override; + + // ----------member data --------------------------- + edm::EDGetTokenT gemRecHits_; + edm::EDGetTokenT > muons_; + edm::EDGetTokenT vertexCollection_; + edm::Service fs; + + MuonServiceProxy* theService_; + edm::ESHandle propagator_; + edm::ESHandle ttrackBuilder_; + edm::ESHandle bField_; + edm::EDGetTokenT lumiScalers_; + // edm::EDGetTokenT> gemDigis_; + // // edm::EDGetTokenT> vfatStatus_; + // edm::EDGetTokenT> gebStatus_; + + ULong64_t b_event; + int b_run, b_lumi; + float b_instLumi; + int b_firstStrip, b_nStrips, b_chamber, b_layer, b_etaPartition, b_muonQuality, b_bx; + int b_latency; + vector b_strips; + float b_x, b_y, b_z; + float b_mu_eta, b_mu_phi, b_mu_pt; + float b_pull_x, b_pull_y, b_res_x, b_res_y; + + int nEvents, nMuonTotal, nGEMFiducialMuon, nGEMTrackWithMuon, nGEMTotal; + int b_nMuons, b_nMuonsWithGEMHit, b_nMuonsInMuonTree, b_nHitsInHitTree; + int b_valid; + + int b_nGEMHits; + + // muon branches + int m_nhits, m_nvalidhits; + int m_nbounds; + int m_quality, m_charge; + float m_pt, m_eta, m_phi; + // GEMHits included in Muon + vector m_roll, m_chamber, m_layer; // hit info + vector m_resx, m_resy, m_pullx, m_pully; + // Propagation only information + vector m_in_vfat; + vector m_in_roll, m_in_chamber, m_in_layer; // propagation bound info + vector m_in_globalPhi, m_in_globalEta, m_in_nearGemPhi, m_in_nearGemEta; // global info + vector m_in_x, m_in_y, m_in_local_x, m_in_local_y, m_in_gemx, m_in_gemy, m_in_local_gemx, m_in_local_gemy, m_in_pullx, m_in_pully, m_in_resx, m_in_resy, m_in_trkextdx, m_in_gem_dx; + vector m_in_local_x_closetsos, m_in_local_y_closetsos, m_in_trkextdx_closetsos, m_in_trkextdx_inner; + vector m_in_local_x_inner, m_in_local_y_inner; + vector m_in_gemNStrips, m_in_gemFirstStrip, m_in_strip; + vector m_in_matchingGem; + + vector m_rec_roll, m_rec_chamber, m_rec_layer; + + vector> m_in_resx_tests; + + TTree *t_hit; + TTree *t_run; + TTree *t_muon; + TTree *t_event; +}; + +GEMDebug::GEMDebug(const edm::ParameterSet& iConfig) : + nEvents(0), + nMuonTotal(0), + nGEMFiducialMuon(0), + nGEMTrackWithMuon(0), + nGEMTotal(0) +{ + gemRecHits_ = consumes(iConfig.getParameter("gemRecHits")); + muons_ = consumes >(iConfig.getParameter("muons")); + edm::ParameterSet serviceParameters = iConfig.getParameter("ServiceParameters"); + theService_ = new MuonServiceProxy(serviceParameters); + lumiScalers_ = consumes(iConfig.getParameter("lumiScalers")); +} + +GEMDebug::~GEMDebug() +{ + std::cout << "::: GEM Slice Test Results :::" << std::endl; + std::cout << ": From " << nEvents << " events" << std::endl; + std::cout << std::endl; + std::cout << " # GEMHits " << nGEMTotal << std::endl; + std::cout << " # Muons " << nMuonTotal << std::endl; + std::cout << " # FidMu " << nGEMFiducialMuon << std::endl; + std::cout << " # GEMMu " << nGEMTrackWithMuon << std::endl; + std::cout << std::endl; +} + +#include + +void +GEMDebug::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + edm::Handle lumiScalers; + iEvent.getByToken(lumiScalers_, lumiScalers); + + edm::ESHandle GEMGeometry_; + iSetup.get().get(GEMGeometry_); + + edm::Handle gemRecHits; + iEvent.getByToken(gemRecHits_, gemRecHits); + + Handle > muons; + iEvent.getByToken(muons_, muons); + + iSetup.get().get("TransientTrackBuilder", ttrackBuilder_); + theService_->update(iSetup); + auto propagator = theService_->propagator("SteppingHelixPropagatorAny"); + bool output = false; + std::stringstream ss; + + ss << "\n\n ---------- \n\n" << "\n"; + ss << "nEvent: " << nEvents++ << " : " << iEvent.id().event() << "\n\n"; + ss << "GEM REC HITS" << "\n"; + for (auto & gem : *gemRecHits) { + auto detId = gem.gemId(); + + auto roll = GEMGeometry_->etaPartition(detId); + auto globalPosition = roll->toGlobal(gem.localPosition()); + b_x = globalPosition.x(); + b_y = globalPosition.y(); + b_z = globalPosition.z(); + + ss << " --- eta: " << globalPosition.eta() << " phi: " << globalPosition.phi() << " xyz: " << b_x << " " << b_y << " " << b_z << " ch:" << detId.chamber() << " l:" << detId.layer() << " p:" << detId.roll() << " s:" << gem.firstClusterStrip() << " [" << gem.clusterSize() << "]" << "\n"; + } + + ss << "\n\nMUONS" << "\n"; + for (auto & mu : *muons) { + ss << " -- Mu e:" << mu.eta() << " p:" << mu.phi() << " pt:" << mu.pt() << "\n"; + // only consider muons going in the right direction (toward the gem slice test) + if (mu.eta() > 0) continue; + + if (mu.passed(reco::Muon::Selector::CutBasedIdTight)) m_quality = 2; + else if (mu.passed(reco::Muon::Selector::CutBasedIdLoose)) m_quality = 1; + else m_quality = 0; + + + const reco::Track* muonTrack = 0; + if (mu.globalTrack().isNonnull()) muonTrack = mu.globalTrack().get(); + else if (mu.outerTrack().isNonnull()) muonTrack = mu.outerTrack().get(); + if (muonTrack) { + reco::TransientTrack ttTrack = ttrackBuilder_->build(muonTrack); + for (auto ch : GEMGeometry_->etaPartitions()) { + + TrajectoryStateOnSurface tsos = propagator->propagate(ttTrack.outermostMeasurementState(), + ch->surface()); + if (!tsos.isValid()) continue; + auto gemid = ch->id(); + if (gemid.chamber() == 1) continue; + + GlobalPoint tsosGP = tsos.globalPosition(); + const LocalPoint pos = ch->toLocal(tsosGP); + const LocalPoint pos2D(pos.x(), pos.y(), 0); + const BoundPlane& bps(ch->surface()); + + if (bps.bounds().inside(pos2D)) { + output = true; + ss << " ONgem: " << gemid.chamber() << " " << gemid.layer() << " " << gemid.roll() << " " << ch->strip(pos2D) << " - xyz: " << tsosGP.x() << " " << tsosGP.y() << " " << tsosGP.z() << "\n"; + } + } + + for (auto hit = muonTrack->recHitsBegin(); hit != muonTrack->recHitsEnd(); hit++) { + if ((*hit)->geographicalId().det() == DetId::Muon && + (*hit)->geographicalId().subdetId() == MuonSubdetId::GEM) { + output = true; + GEMDetId gemid((*hit)->geographicalId()); + + ss << " found hit " << gemid.chamber() << " " << gemid.layer() << " " << gemid.roll() << "\n"; + } + } + } + } + + ss << " ---------- " << "\n"; + if (output) { cout << ss.str() << endl; } +} + +void GEMDebug::beginJob(){} +void GEMDebug::endJob(){} + +void GEMDebug::beginRun(Run const& run, EventSetup const&){ +} +void GEMDebug::endRun(Run const&, EventSetup const&){} + +//define this as a plug-in +DEFINE_FWK_MODULE(GEMDebug); diff --git a/MuonAnalyser/plugins/MuonDetHitAnalyser.cc b/MuonAnalyser/plugins/MuonDetHitAnalyser.cc new file mode 100644 index 00000000000..b6dfae3c1fb --- /dev/null +++ b/MuonAnalyser/plugins/MuonDetHitAnalyser.cc @@ -0,0 +1,769 @@ +// cd /cms/ldap_home/iawatson/scratch/GEM/CMSSW_10_1_5/src/ && eval `scramv1 runtime -sh` && eval `scramv1 runtime -sh` && scram b -j 10 +// cd ../../.. && source /cvmfs/cms.cern.ch/cmsset_default.sh && eval `scramv1 runtime -sh` && eval `scramv1 runtime -sh` && scram b -j 10 +// system include files +#include +#include +#include +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" +// ME0 +#include "DataFormats/GEMDigi/interface/ME0DigiCollection.h" +#include "DataFormats/GEMRecHit/interface/ME0RecHitCollection.h" +#include "DataFormats/GEMRecHit/interface/ME0SegmentCollection.h" +#include "DataFormats/MuonDetId/interface/ME0DetId.h" +#include "Geometry/GEMGeometry/interface/ME0Geometry.h" +#include "Geometry/GEMGeometry/interface/ME0EtaPartition.h" +#include "Geometry/GEMGeometry/interface/ME0EtaPartitionSpecs.h" +// GEM +#include "DataFormats/GEMDigi/interface/GEMDigiCollection.h" +#include "DataFormats/GEMRecHit/interface/GEMRecHitCollection.h" +#include "DataFormats/GEMRecHit/interface/GEMSegmentCollection.h" +#include "DataFormats/MuonDetId/interface/GEMDetId.h" +#include "Geometry/GEMGeometry/interface/GEMGeometry.h" +#include "Geometry/GEMGeometry/interface/GEMEtaPartition.h" +#include "Geometry/GEMGeometry/interface/GEMEtaPartitionSpecs.h" +// CSC +#include "DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h" +#include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h" +#include "DataFormats/MuonDetId/interface/CSCDetId.h" +#include "Geometry/CSCGeometry/interface/CSCGeometry.h" +// DT +#include "DataFormats/DTDigi/interface/DTDigiCollection.h" +#include "DataFormats/DTRecHit/interface/DTRecHitCollection.h" +//#include "DataFormats/DTRecHit/interface/DTRecSegment2DCollection.h" +#include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h" +#include "DataFormats/MuonDetId/interface/DTBtiId.h" +#include "DataFormats/MuonDetId/interface/DTChamberId.h" +#include "DataFormats/MuonDetId/interface/DTLayerId.h" +#include "DataFormats/MuonDetId/interface/DTSectCollId.h" +#include "DataFormats/MuonDetId/interface/DTSuperLayerId.h" +#include "DataFormats/MuonDetId/interface/DTTracoId.h" +#include "DataFormats/MuonDetId/interface/DTWireId.h" +#include "Geometry/DTGeometry/interface/DTGeometry.h" +// RPC +#include "DataFormats/RPCDigi/interface/RPCDigiCollection.h" +#include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h" +#include "DataFormats/MuonDetId/interface/RPCDetId.h" +#include "Geometry/RPCGeometry/interface/RPCGeometry.h" +#include "Geometry/RPCGeometry/interface/RPCGeomServ.h" +// Muon +#include "DataFormats/PatCandidates/interface/Muon.h" +// L1 +#include "DataFormats/L1TMuon/interface/L1MuBMTrack.h" + +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" + +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" +#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" +#include "SimDataFormats/Associations/interface/MuonToTrackingParticleAssociator.h" +#include "SimDataFormats/Vertex/interface/SimVertex.h" +#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h" +#include "SimMuon/MCTruth/interface/MuonToSimAssociatorByHits.h" +#include "SimTracker/Common/interface/TrackingParticleSelector.h" +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h" + +#include "Geometry/Records/interface/MuonGeometryRecord.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" + +#include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h" +#include "Geometry/CommonTopologies/interface/RadialStripTopology.h" +#include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h" + +#include "Geometry/Records/interface/MuonGeometryRecord.h" + +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/Run.h" + +#include "TH1D.h" +#include "TH2D.h" +#include "TProfile.h" +#include "TString.h" +#include "TGraphAsymmErrors.h" +#include "TLorentzVector.h" +#include "TTree.h" + +using namespace std; +using namespace edm; + +class MuonDetHitAnalyser : public edm::EDAnalyzer { +public: + explicit MuonDetHitAnalyser(const edm::ParameterSet&); + ~MuonDetHitAnalyser(); + +private: + virtual void analyze(const edm::Event&, const edm::EventSetup&); + virtual void beginJob() override; + virtual void endJob() override; + + virtual void beginRun(Run const&, EventSetup const&) override; + virtual void endRun(Run const&, EventSetup const&) override; + + void initValue(); + + TProfile* hME0Area_; + TH1D* hME0Counts_; + + TProfile* hGEMArea_; + TH1D* hGEMCounts_; + + TProfile* hCSCArea_; + TH1D* hCSCCounts_; + + TProfile* hRPCArea_; + TH1D* hRPCCounts_; + + TH1D* hEvents_; + + // ----------member data --------------------------- + edm::EDGetTokenT me0Digis_; + edm::EDGetTokenT me0RecHits_; + + edm::EDGetTokenT csc2DRecHits_; + + edm::EDGetTokenT gemDigis_; + edm::EDGetTokenT gemRecHits_; + + edm::EDGetTokenT dtDigis_; + edm::EDGetTokenT dtRecHits_; + + edm::EDGetTokenT rpcDigis_; + edm::EDGetTokenT rpcRecHits_; + + edm::Service fs; + + TTree *t_event; + int b_run, b_lumi, b_latency; + int b_nME0Digis, b_nME0Segments, b_nME0RecHits, b_nCSCSegments, b_nCSC2DRecHits, b_nGEMDigis, b_nGEMSegments, b_nGEMRecHits, b_nDTDigis, b_nDT4DSegments, b_nDTRecHits, b_nRPCDigis, b_nRPCRecHits; + + /* ME0 */ + TTree *t_ME0_digi; + int b_ME0_Digi_region, b_ME0_Digi_chamber, b_ME0_Digi_layer, b_ME0_Digi_etaPartition, b_ME0_Digi_bx; + float b_ME0_Digi_eta, b_ME0_Digi_phi; + TTree *t_ME0_rec; + int b_ME0_RecHit_region, b_ME0_RecHit_chamber, b_ME0_RecHit_layer, b_ME0_RecHit_etaPartition; + float b_ME0_RecHit_eta, b_ME0_RecHit_phi; + + /* CSC */ + TTree *t_CSC_rec; + int b_CSC_2DRecHit_region, b_CSC_2DRecHit_station, b_CSC_2DRecHit_ring, b_CSC_2DRecHit_chamber, b_CSC_2DRecHit_layer; + float b_CSC_2DRecHit_eta, b_CSC_2DRecHit_phi; + + /* GEM */ + TTree *t_GEM_digi; + int b_GEM_Digi_region, b_GEM_Digi_station, b_GEM_Digi_ring, b_GEM_Digi_chamber, b_GEM_Digi_layer, b_GEM_Digi_etaPartition, b_GEM_Digi_bx; + float b_GEM_Digi_eta, b_GEM_Digi_phi; + TTree *t_GEM_rec; + int b_GEM_RecHit_region, b_GEM_RecHit_station, b_GEM_RecHit_ring, b_GEM_RecHit_chamber, b_GEM_RecHit_layer, b_GEM_RecHit_etaPartition, b_GEM_RecHit_bx; + float b_GEM_RecHit_eta, b_GEM_RecHit_phi; + + /* DT */ + TTree *t_DT_digi; + int b_DT_Digi_sector, b_DT_Digi_station, b_DT_Digi_wheel, b_DT_Digi_chamber, b_DT_Digi_superLayer, b_DT_Digi_layer, b_DT_Digi_wire; + float b_DT_Digi_eta, b_DT_Digi_phi; + TTree *t_DT_rec; + int b_DT_RecHit_sector, b_DT_RecHit_station, b_DT_RecHit_wheel, b_DT_RecHit_chamber, b_DT_RecHit_superLayer, b_DT_RecHit_layer, b_DT_RecHit_wire; + float b_DT_RecHit_eta, b_DT_RecHit_phi; + + /* RPC */ + TTree *t_RPC_digi; + int b_RPC_Digi_region, b_RPC_Digi_sector, b_RPC_Digi_subSector, b_RPC_Digi_station, b_RPC_Digi_ring, b_RPC_Digi_chamber, b_RPC_Digi_layer, b_RPC_Digi_roll, b_RPC_Digi_isIRPC; + float b_RPC_Digi_eta, b_RPC_Digi_phi; + TTree *t_RPC_rec; + int b_RPC_RecHit_region, b_RPC_RecHit_sector, b_RPC_RecHit_subSector, b_RPC_RecHit_station, b_RPC_RecHit_ring, b_RPC_RecHit_chamber, b_RPC_RecHit_layer, b_RPC_RecHit_roll, b_RPC_RecHit_isIRPC; + float b_RPC_RecHit_eta, b_RPC_RecHit_phi; +}; + +MuonDetHitAnalyser::MuonDetHitAnalyser(const edm::ParameterSet& iConfig) +{ + me0Digis_ = consumes(iConfig.getParameter("me0Digis")); + me0RecHits_ = consumes(iConfig.getParameter("me0RecHits")); + + csc2DRecHits_ = consumes(iConfig.getParameter("csc2DRecHits")); + + gemDigis_ = consumes(iConfig.getParameter("gemDigis")); + gemRecHits_ = consumes(iConfig.getParameter("gemRecHits")); + + dtDigis_ = consumes(iConfig.getParameter("dtDigis")); + dtRecHits_ = consumes(iConfig.getParameter("dtRecHits")); + + rpcDigis_ = consumes(iConfig.getParameter("rpcDigis")); + rpcRecHits_ = consumes(iConfig.getParameter("rpcRecHits")); + + t_event = fs->make("Event", "Event"); + t_event->Branch("nME0Digis", &b_nME0Digis, "nME0Digis/I"); + t_event->Branch("nME0RecHits", &b_nME0RecHits, "nME0RecHits/I"); + t_event->Branch("nCSC2DRecHits", &b_nCSC2DRecHits, "nCSC2DRecHits/I"); + t_event->Branch("nGEMDigis", &b_nGEMDigis, "nGEMDigis/I"); + t_event->Branch("nGEMRecHits", &b_nGEMRecHits, "nGEMRecHits/I"); + t_event->Branch("nDTDigis", &b_nDTDigis, "nDTDigis/I"); + t_event->Branch("nDTRecHits", &b_nDTRecHits, "nDTRecHits/I"); + t_event->Branch("nRPCDigis", &b_nRPCDigis, "nRPCDigis/I"); + t_event->Branch("nRPCRecHits", &b_nRPCRecHits, "nRPCRecHits/I"); + + /*ME0*/ + t_ME0_digi = fs->make("ME0_Hit", "ME0_Hit"); + t_ME0_digi->Branch("Digi_eta", &b_ME0_Digi_eta, "Digi_eta/F"); + t_ME0_digi->Branch("Digi_phi", &b_ME0_Digi_phi, "Digi_phi/F"); + t_ME0_digi->Branch("Digi_bx", &b_ME0_Digi_bx, "Digi_bx/I"); + t_ME0_digi->Branch("Digi_region", &b_ME0_Digi_region, "Digi_region/I"); + t_ME0_digi->Branch("Digi_chamber", &b_ME0_Digi_chamber, "Digi_chamber/I"); + t_ME0_digi->Branch("Digi_layer", &b_ME0_Digi_layer, "Digi_layer/I"); + t_ME0_digi->Branch("Digi_etaPartition", &b_ME0_Digi_etaPartition, "Digi_etaPartition/I"); + + t_ME0_rec = fs->make("ME0_RecHit", "ME0_RecHit"); + t_ME0_rec->Branch("RecHit_eta", &b_ME0_RecHit_eta, "RecHit_eta/F"); + t_ME0_rec->Branch("RecHit_phi", &b_ME0_RecHit_phi, "RecHit_phi/F"); + t_ME0_rec->Branch("RecHit_region", &b_ME0_RecHit_region, "RecHit_region/I"); + t_ME0_rec->Branch("RecHit_chamber", &b_ME0_RecHit_chamber, "RecHit_chaber/I"); + t_ME0_rec->Branch("RecHit_layer", &b_ME0_RecHit_layer, "RecHit_layer/I"); + t_ME0_rec->Branch("RecHit_etaPartition", &b_ME0_RecHit_etaPartition, "RecHit_etaParition/I"); + + /*CSC*/ + t_CSC_rec = fs->make("CSC_2DRecHit", "CSC_2DRecHit"); + t_CSC_rec->Branch("RecHit_eta", &b_CSC_2DRecHit_eta, "RecHit_eta/F"); + t_CSC_rec->Branch("RecHit_phi", &b_CSC_2DRecHit_phi, "RecHit_phi/F"); + t_CSC_rec->Branch("RecHit_region", &b_CSC_2DRecHit_region, "RecHit_region/I"); + t_CSC_rec->Branch("RecHit_station", &b_CSC_2DRecHit_station, "RecHit_station/I"); + t_CSC_rec->Branch("RecHit_ring", &b_CSC_2DRecHit_ring, "RecHit_ring/I"); + t_CSC_rec->Branch("RecHit_chamber", &b_CSC_2DRecHit_chamber, "RecHit_chaber/I"); + t_CSC_rec->Branch("RecHit_layer", &b_CSC_2DRecHit_layer, "RecHit_layer/I"); + + /*GEM*/ + t_GEM_digi = fs->make("GEM_Hit", "GEM_Hit"); + t_GEM_digi->Branch("Digi_eta", &b_GEM_Digi_eta, "Digi_eta/F"); + t_GEM_digi->Branch("Digi_phi", &b_GEM_Digi_phi, "Digi_phi/F"); + t_GEM_digi->Branch("Digi_bx", &b_GEM_Digi_bx, "Digi_bx/I"); + t_GEM_digi->Branch("Digi_region", &b_GEM_Digi_region, "Digi_region/I"); + t_GEM_digi->Branch("Digi_station", &b_GEM_Digi_station, "Digi_station/I"); + t_GEM_digi->Branch("Digi_ring", &b_GEM_Digi_ring, "Digi_ring/I"); + t_GEM_digi->Branch("Digi_chamber", &b_GEM_Digi_chamber, "Digi_chamber/I"); + t_GEM_digi->Branch("Digi_layer", &b_GEM_Digi_layer, "Digi_layer/I"); + t_GEM_digi->Branch("Digi_etaPartition", &b_GEM_Digi_etaPartition, "Digi_etaPartition/I"); + + t_GEM_rec = fs->make("GEM_RecHit", "GEM_RecHit"); + t_GEM_rec->Branch("RecHit_eta", &b_GEM_RecHit_eta, "RecHit_eta/F"); + t_GEM_rec->Branch("RecHit_phi", &b_GEM_RecHit_phi, "RecHit_phi/F"); + t_GEM_rec->Branch("RecHit_bx", &b_GEM_RecHit_bx, "RecHit_bx/I"); + t_GEM_rec->Branch("RecHit_region", &b_GEM_RecHit_region, "RecHit_region/I"); + t_GEM_rec->Branch("RecHit_station", &b_GEM_RecHit_station, "RecHit_station/I"); + t_GEM_rec->Branch("RecHit_ring", &b_GEM_RecHit_ring, "RecHit_ring/I"); + t_GEM_rec->Branch("RecHit_chamber", &b_GEM_RecHit_chamber, "RecHit_chaber/I"); + t_GEM_rec->Branch("RecHit_layer", &b_GEM_RecHit_layer, "RecHit_layer/I"); + t_GEM_rec->Branch("RecHit_etaPartition", &b_GEM_RecHit_etaPartition, "RecHit_etaParition/I"); + + /*DT*/ + t_DT_digi = fs->make("DT_Hit", "DT_Hit"); + t_DT_digi->Branch("Digi_eta", &b_DT_Digi_eta, "Digi_eta/F"); + t_DT_digi->Branch("Digi_phi", &b_DT_Digi_phi, "Digi_phi/F"); + t_DT_digi->Branch("Digi_sector", &b_DT_Digi_sector, "Digi_sector/I"); + t_DT_digi->Branch("Digi_station", &b_DT_Digi_station, "Digi_station/I"); + t_DT_digi->Branch("Digi_wheel", &b_DT_Digi_wheel, "Digi_wheel/I"); + t_DT_digi->Branch("Digi_chamber", &b_DT_Digi_chamber, "Digi_chamber/I"); + t_DT_digi->Branch("Digi_superLayer", &b_DT_Digi_superLayer, "Digi_superLayer/I"); + t_DT_digi->Branch("Digi_layer", &b_DT_Digi_layer, "Digi_layer/I"); + t_DT_digi->Branch("Digi_wire", &b_DT_Digi_wire, "Digi_wire/I"); + + t_DT_rec = fs->make("DT_RecHit", "DT_RecHit"); + t_DT_rec->Branch("RecHit_eta", &b_DT_RecHit_eta, "RecHit_eta/F"); + t_DT_rec->Branch("RecHit_phi", &b_DT_RecHit_phi, "RecHit_phi/F"); + t_DT_rec->Branch("RecHit_sector", &b_DT_RecHit_sector, "RecHit_sector/I"); + t_DT_rec->Branch("RecHit_station", &b_DT_RecHit_station, "RecHit_station/I"); + t_DT_rec->Branch("RecHit_wheel", &b_DT_RecHit_wheel, "RecHit_wheel/I"); + t_DT_rec->Branch("RecHit_chamber", &b_DT_RecHit_chamber, "RecHit_chaber/I"); + t_DT_rec->Branch("RecHit_superLayer", &b_DT_RecHit_superLayer, "RecHit_superLayer/I"); + t_DT_rec->Branch("RecHit_layer", &b_DT_RecHit_layer, "RecHit_layer/I"); + t_DT_rec->Branch("RecHit_wire", &b_DT_RecHit_wire, "RecHit_wire/I"); + + /*RPC*/ + t_RPC_digi = fs->make("RPC_Hit", "RPC_Hit"); + t_RPC_digi->Branch("Digi_eta", &b_RPC_Digi_eta, "Digi_eta/F"); + t_RPC_digi->Branch("Digi_phi", &b_RPC_Digi_phi, "Digi_phi/F"); + t_RPC_digi->Branch("Digi_isIRPC", &b_RPC_Digi_isIRPC, "Digi_isIRPC/I"); + t_RPC_digi->Branch("Digi_region", &b_RPC_Digi_region, "Digi_region/I"); + t_RPC_digi->Branch("Digi_sector", &b_RPC_Digi_sector, "Digi_sector/I"); + t_RPC_digi->Branch("Digi_subSector", &b_RPC_Digi_subSector, "Digi_subSector/I"); + t_RPC_digi->Branch("Digi_station", &b_RPC_Digi_station, "Digi_station/I"); + t_RPC_digi->Branch("Digi_ring", &b_RPC_Digi_ring, "Digi_ring/I"); + t_RPC_digi->Branch("Digi_chamber", &b_RPC_Digi_chamber, "Digi_chamber/I"); + t_RPC_digi->Branch("Digi_layer", &b_RPC_Digi_layer, "Digi_layer/I"); + t_RPC_digi->Branch("Digi_roll", &b_RPC_Digi_roll, "Digi_roll/I"); + t_RPC_digi->Branch("Digi_isIRPC", &b_RPC_Digi_isIRPC, "Digi_isIRPC/I"); + t_RPC_digi->Branch("Digi_eta", &b_RPC_Digi_eta, "Digi_eta/F"); + t_RPC_digi->Branch("Digi_phi", &b_RPC_Digi_phi, "Digi_phi/F"); + + t_RPC_rec = fs->make("RPC_RecHit", "RPC_RecHit"); + t_RPC_rec->Branch("RecHit_eta", &b_RPC_RecHit_eta, "RecHit_eta/F"); + t_RPC_rec->Branch("RecHit_phi", &b_RPC_RecHit_phi, "RecHit_phi/F"); + t_RPC_rec->Branch("RecHit_isIRPC", &b_RPC_RecHit_isIRPC, "RecHit_isIRPC/I"); + t_RPC_rec->Branch("RecHit_region", &b_RPC_RecHit_region, "RecHit_region/I"); + t_RPC_rec->Branch("RecHit_sector", &b_RPC_RecHit_sector, "RecHit_sector/I"); + t_RPC_rec->Branch("RecHit_subSector", &b_RPC_RecHit_subSector, "RecHit_subSector/I"); + t_RPC_rec->Branch("RecHit_station", &b_RPC_RecHit_station, "RecHit_station/I"); + t_RPC_rec->Branch("RecHit_ring", &b_RPC_RecHit_ring, "RecHit_ring/I"); + t_RPC_rec->Branch("RecHit_chamber", &b_RPC_RecHit_chamber, "RecHit_chaber/I"); + t_RPC_rec->Branch("RecHit_layer", &b_RPC_RecHit_layer, "RecHit_layer/I"); + t_RPC_rec->Branch("RecHit_roll", &b_RPC_RecHit_roll, "RecHit_roll/I"); +} + +MuonDetHitAnalyser::~MuonDetHitAnalyser() +{ +} + +void +MuonDetHitAnalyser::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + edm::ESHandle hME0Geom; + iSetup.get().get(hME0Geom); + const ME0Geometry* ME0Geometry_ = &*hME0Geom; + + /* ME0 Geometry */ + edm::Handle me0Digis; + iEvent.getByToken(me0Digis_, me0Digis); + + edm::Handle me0RecHits; + iEvent.getByToken(me0RecHits_, me0RecHits); + + /* CSC Geometry */ + edm::ESHandle hCSCGeom; + iSetup.get().get(hCSCGeom); + const CSCGeometry* CSCGeometry_ = &*hCSCGeom; + + edm::Handle csc2DRecHits; + iEvent.getByToken(csc2DRecHits_, csc2DRecHits); + + /* GEM Geometry */ + edm::ESHandle hGEMGeom; + iSetup.get().get(hGEMGeom); + const GEMGeometry* GEMGeometry_ = &*hGEMGeom; + + edm::Handle gemDigis; + iEvent.getByToken(gemDigis_, gemDigis); + + edm::Handle gemRecHits; + iEvent.getByToken(gemRecHits_, gemRecHits); + + /* DT Geometry */ + edm::ESHandle hDTGeom; + iSetup.get().get(hDTGeom); + const DTGeometry* DTGeometry_ = &*hDTGeom; + + edm::Handle dtDigis; + iEvent.getByToken(dtDigis_, dtDigis); + + edm::Handle dtRecHits; + iEvent.getByToken(dtRecHits_, dtRecHits); + + /* RPC Geometry */ + edm::ESHandle hRPCGeom; + iSetup.get().get(hRPCGeom); + const RPCGeometry* RPCGeometry_ = &*hRPCGeom; + + edm::Handle rpcDigis; + iEvent.getByToken(rpcDigis_, rpcDigis); + + edm::Handle rpcRecHits; + iEvent.getByToken(rpcRecHits_, rpcRecHits); + + initValue(); + + hEvents_->Fill(1); + + /* ME0 */ + for (auto roll : ME0Geometry_->etaPartitions()) { + ME0DetId rId = roll->id(); + // ME0 digi + auto digisRange = me0Digis->get(rId); + auto me0Digi = digisRange.first; + for (auto hit = me0Digi; hit != digisRange.second; ++hit) { + auto strip = hit->strip(); + auto digiLp = roll->centreOfStrip(strip); + auto gp = roll->toGlobal(digiLp); + b_ME0_Digi_eta = gp.eta(); + b_ME0_Digi_phi = gp.phi(); + b_ME0_Digi_bx = hit->bx(); + b_ME0_Digi_region = rId.region(); // Region id: 0 for Barrel Not in use, +/-1 For +/- Endcap + b_ME0_Digi_chamber = rId.chamber(); + b_ME0_Digi_layer = rId.layer(); + b_ME0_Digi_etaPartition = rId.roll(); + t_ME0_digi->Fill(); + b_nME0Digis++; + } + // ME0 RecHit + auto recRange = me0RecHits->get(rId); + auto me0Rec = recRange.first; + for (auto rec = me0Rec; rec != recRange.second; ++rec) { + auto recLd = rec->localPosition(); + auto gp = roll->toGlobal(recLd); + b_ME0_RecHit_eta = gp.eta(); + b_ME0_RecHit_phi = gp.phi(); + b_ME0_RecHit_region = rId.region(); + b_ME0_RecHit_chamber = rId.chamber(); + b_ME0_RecHit_layer = rId.layer(); + b_ME0_RecHit_etaPartition = rId.roll(); + b_nME0RecHits++; + t_ME0_rec->Fill(); + } + } + for ( auto me0Hit : *me0RecHits ) { + const auto detId = me0Hit.me0Id(); + auto region = std::to_string(detId.region()); + auto ch = std::to_string(detId.chamber()); + auto iEta = std::to_string(detId.roll()); + auto station = std::to_string(detId.station()); + auto ring = std::to_string(0);//std::to_string(detId.ring()); + auto layer = std::to_string(detId.layer()); + const string rollName = "Region_"+region+"_station_"+station+"_ring_"+ring+"_CH_"+ch+"_layer_"+layer+"_iEta_"+iEta; + const int idx = hME0Area_->GetXaxis()->FindBin(rollName.c_str()); + hME0Counts_->Fill(idx); + } + + /* CSC */ + for (auto ly : CSCGeometry_->layers()) { + CSCDetId lId = ly->id(); + // CSC rechit + auto recRange = csc2DRecHits->get(lId); + auto cscRec = recRange.first; + for (auto rec = cscRec; rec != recRange.second; ++rec) { + auto recLd = rec->localPosition(); + auto gp = ly->toGlobal(recLd); + b_CSC_2DRecHit_eta = gp.eta(); + b_CSC_2DRecHit_phi = gp.phi(); + b_CSC_2DRecHit_region = lId.endcap(); + b_CSC_2DRecHit_station = lId.station(); + b_CSC_2DRecHit_ring = lId.ring(); + b_CSC_2DRecHit_chamber = lId.chamber(); + b_CSC_2DRecHit_layer = lId.layer(); + b_nCSC2DRecHits++; + t_CSC_rec->Fill(); + } + } + for ( auto cscHit : *csc2DRecHits ) { + const auto detId = cscHit.cscDetId(); + auto endcap = std::to_string(detId.endcap()); // +1 : forward , -1 : backward + auto ch = std::to_string(detId.chamber()); + auto ring = std::to_string(detId.ring()); + auto station = std::to_string(detId.station()); + auto layer = std::to_string(detId.layer()); + const string name = "Endcap_"+endcap+"_station_"+station+"_ring_"+ring+"_CH_"+ch+"_layer_"+layer; + const int idx = hCSCArea_->GetXaxis()->FindBin(name.c_str()); + hCSCCounts_->Fill(idx); + } + + + /* GEM */ + for (auto roll : GEMGeometry_->etaPartitions()) { + GEMDetId rId = roll->id(); + // GEM digi + auto digisRange = gemDigis->get(rId); + auto gemDigi = digisRange.first; + for (auto hit = gemDigi; hit != digisRange.second; ++hit) { + auto strip = hit->strip(); + auto digiLp = roll->centreOfStrip(strip); + auto gp = roll->toGlobal(digiLp); + b_GEM_Digi_eta = gp.eta(); + b_GEM_Digi_phi = gp.phi(); + b_GEM_Digi_bx = hit->bx(); + b_GEM_Digi_region = rId.region(); + b_GEM_Digi_station = rId.station(); + b_GEM_Digi_ring = rId.ring(); + b_GEM_Digi_chamber = rId.chamber(); + b_GEM_Digi_layer = rId.layer(); + b_GEM_Digi_etaPartition = rId.roll(); + t_GEM_digi->Fill(); + b_nGEMDigis++; + } + // GEM rechit + auto recRange = gemRecHits->get(rId); + auto gemRec = recRange.first; + for (auto rec = gemRec; rec != recRange.second; ++rec) { + auto recLd = rec->localPosition(); + auto gp = roll->toGlobal(recLd); + b_GEM_RecHit_eta = gp.eta(); + b_GEM_RecHit_phi = gp.phi(); + b_GEM_RecHit_bx = rec->BunchX(); + b_GEM_RecHit_region = rId.region(); + b_GEM_RecHit_station = rId.station(); + b_GEM_RecHit_ring = rId.ring(); + b_GEM_RecHit_chamber = rId.chamber(); + b_GEM_RecHit_layer = rId.layer(); + b_GEM_RecHit_etaPartition = rId.roll(); + b_nGEMRecHits++; + t_GEM_rec->Fill(); + } + } + for ( auto gemHit : *gemRecHits ) { + const auto detId = gemHit.gemId(); + auto region = std::to_string(detId.region()); + auto ch = std::to_string(detId.chamber()); + auto iEta = std::to_string(detId.roll()); + auto station = std::to_string(detId.station()); + auto ring = std::to_string(detId.ring()); + auto layer = std::to_string(detId.layer()); + const string rollName = "Region_"+region+"_station_"+station+"_ring_"+ring+"_CH_"+ch+"_layer_"+layer+"_iEta_"+iEta; + + const int idx = hGEMArea_->GetXaxis()->FindBin(rollName.c_str()); + hGEMCounts_->Fill(idx); + } + + /* DT */ + for (auto ly : DTGeometry_->layers()) { + DTLayerId lId = ly->id(); + //auto ch = ly->chamber(); + //DTChamberId cId = ch->id(); + // DT digi + auto digisRange = dtDigis->get(lId); + auto dtDigi = digisRange.first; + for (auto hit = dtDigi; hit != digisRange.second; ++hit) { + b_DT_Digi_sector = lId.sector(); + b_DT_Digi_station = lId.station(); + b_DT_Digi_wheel = lId.wheel(); + //b_DT_Digi_chamber = cId.chamber(); + b_DT_Digi_superLayer = lId.superLayer(); /// Superlayers are numbered 1 (phi), 2 (Z), 3 (phi) + b_DT_Digi_layer = lId.layer(); + b_DT_Digi_wire = hit->wire(); + t_DT_digi->Fill(); + b_nDTDigis++; + } + // DT rechit + auto recRange = dtRecHits->get(lId); + auto dtRec = recRange.first; + for (auto rec = dtRec; rec != recRange.second; ++rec) { + auto recLd = rec->localPosition(); + auto gp = ly->toGlobal(recLd); + auto wireId = rec->wireId(); + b_DT_RecHit_eta = gp.eta(); + b_DT_RecHit_phi = gp.phi(); + b_DT_RecHit_sector = lId.sector(); + b_DT_RecHit_station = lId.station(); + b_DT_RecHit_wheel = lId.wheel(); + //b_DT_RecHit_chamber = cId.chamber() + b_DT_RecHit_superLayer = lId.superLayer(); + b_DT_RecHit_layer = lId.layer(); + b_DT_RecHit_wire = wireId.wire(); + b_nDTRecHits++; + t_DT_rec->Fill(); + } + } + + /* RPC */ + for (auto roll : RPCGeometry_->rolls()) { + RPCDetId rId = roll->id(); + //auto ch = roll->chamber(); + //auto cId = ch->id(); + // RPC digi + auto digisRange = rpcDigis->get(rId); + auto rpcDigi = digisRange.first; + for (auto hit = rpcDigi; hit != digisRange.second; ++hit) { + auto strip = hit->strip(); + auto digiLp = roll->centreOfStrip(strip); + auto gp = roll->toGlobal(digiLp); + b_RPC_Digi_eta = gp.eta(); + b_RPC_Digi_phi = gp.phi(); + b_RPC_Digi_isIRPC = (int)roll->isIRPC(); + b_RPC_Digi_region = rId.region(); // region : 0 for barrel +/-1 for +/- endcap + b_RPC_Digi_sector = rId.sector(); // Sector id: the group of chambers at same phi (and increasing r) + b_RPC_Digi_subSector = rId.subsector(); // subSector id : some sectors are divided along the phi direction in subsectors (from 1 to 4 in Barrel, from 1 to 6 in Endcap) + b_RPC_Digi_station = rId.station(); // Station id : For Barrel: the four groups of chambers at same r (distance from beam axis) and increasing phi / For Endcap: the three groups of chambers at same z (distance from interaction point), i.e. the disk + b_RPC_Digi_ring = rId.ring(); // Ring id: Wheel number in Barrel (from -2 to +2) Ring Number in Endcap (from 1 to 3) Ring has a different meaning in Barrel and Endcap! In Barrel it is wheel, in Endcap Ring has a different meaning in + //b_RPC_Digi_chamber = cId.chamber(); + b_RPC_Digi_layer = rId.layer(); // Layer id: each station can have two layers of chambers: layer 1 is the inner chamber and layer 2 is the outer chamber (when present) // Only in Barrel: RB1 and RB + b_RPC_Digi_roll = rId.roll(); // Roll id (also known as eta partition): each chamber is divided along the strip direction in - two or three parts (rolls) for Barrel and two, three or four parts for endcap + t_RPC_digi->Fill(); + b_nRPCDigis++; + } + // RPC rechit + auto recRange = rpcRecHits->get(rId); + auto rpcRec = recRange.first; + for (auto rec = rpcRec; rec != recRange.second; ++rec) { + auto recLd = rec->localPosition(); + auto gp = roll->toGlobal(recLd); + b_RPC_RecHit_eta = gp.eta(); + b_RPC_RecHit_phi = gp.phi(); + b_RPC_RecHit_isIRPC = (int)roll->isIRPC(); + b_RPC_RecHit_region = rId.region(); + b_RPC_RecHit_sector = rId.sector(); + b_RPC_RecHit_subSector = rId.subsector(); + b_RPC_RecHit_station = rId.station(); + b_RPC_RecHit_ring = rId.ring(); + //b_RPC_RecHit_chamber = cId.chamber(); + b_RPC_RecHit_layer = rId.layer(); + b_RPC_RecHit_roll = rId.roll(); + b_nRPCRecHits++; + t_RPC_rec->Fill(); + } + } + for ( auto rpcHit : *rpcRecHits ) { + const auto detId = rpcHit.rawId(); + const string rollName = RPCGeomServ(detId).name(); + const int idx = hRPCArea_->GetXaxis()->FindBin(rollName.c_str()); + hRPCCounts_->Fill(idx); + } + t_event->Fill(); +} + +void MuonDetHitAnalyser::beginJob(){} +void MuonDetHitAnalyser::endJob(){} + +void MuonDetHitAnalyser::beginRun(const edm::Run& run, const edm::EventSetup& iSetup){//(Run const& run, EventSetup const&){ + + hEvents_ = fs->make("hEvents", "hEvents;Types;Number of Events", 1, 1, 2); + + hME0Area_ = fs->make("hME0Area", "Roll area;Roll Name;Area [cm^{2}]", 5000, 1, 5001, 0, 100000); + hME0Counts_ = fs->make("hME0Counts", "Counts;Roll Index;Number of RecHits", 5000, 1, 5001); + + hGEMArea_ = fs->make("hGEMArea", "Roll area;Roll Name;Area [cm^{2}]", 5000, 1, 5001, 0, 100000); + hGEMCounts_ = fs->make("hGEMCounts", "Counts;Roll Index;Number of RecHits", 5000, 1, 5001); + + hCSCArea_ = fs->make("hCSCArea", "Roll area;Roll Name;Area [cm^{2}]", 5000, 1, 5001, 0, 100000); + hCSCCounts_ = fs->make("hCSCCounts", "Counts;Roll Index;Number of RecHits", 5000, 1, 5001); + + hRPCArea_ = fs->make("hRPCArea", "Roll area;Roll Name;Area [cm^{2}]", 5000, 1, 5001, 0, 100000); + hRPCCounts_ = fs->make("hRPCCounts", "Counts;Roll Index;Number of RecHits", 5000, 1, 5001); + + + edm::ESHandle me0Geom; + iSetup.get().get(me0Geom); + + edm::ESHandle gemGeom; + iSetup.get().get(gemGeom); + + edm::ESHandle cscGeom; + iSetup.get().get(cscGeom); + + edm::ESHandle rpcGeom; + iSetup.get().get(rpcGeom); + + /* ME0 */ + int iME0 = 0; + for ( auto roll : me0Geom->etaPartitions() ) { + ++iME0; + const auto detId = roll->id(); + auto region = std::to_string(detId.region()); + auto ch = std::to_string(detId.chamber()); + auto iEta = std::to_string(detId.roll()); + auto station = std::to_string(detId.station()); // for now, always 1 + auto ring = std::to_string(0);//std::to_string(detId.ring()); + auto layer = std::to_string(detId.layer()); + const string rollName = "Region_"+region+"_station_"+station+"_ring_"+ring+"_CH_"+ch+"_layer_"+layer+"_iEta_"+iEta; + const TrapezoidalStripTopology* top_(dynamic_cast(&(roll->topology()))); + const float striplength(top_->stripLength()); + const float pitch(roll->pitch()); + auto nstrip = top_->nstrips(); +// std::cout << " ME0 : " << rollName.c_str() << " / area : " << striplength*pitch*nstrip << std::endl; + hME0Area_->GetXaxis()->SetBinLabel(iME0, rollName.c_str()); + hME0Area_->Fill(iME0, striplength*pitch*nstrip); + } + + /* GEM */ + int iGEM = 0; + for ( auto roll : gemGeom->etaPartitions() ) { + ++iGEM; + const auto detId = roll->id(); + auto region = std::to_string(detId.region()); + auto ch = std::to_string(detId.chamber()); + auto iEta = std::to_string(detId.roll()); + auto station = std::to_string(detId.station()); + auto ring = std::to_string(detId.ring()); + auto layer = std::to_string(detId.layer()); + const string rollName = "Region_"+region+"_station_"+station+"_ring_"+ring+"_CH_"+ch+"_layer_"+layer+"_iEta_"+iEta; + const TrapezoidalStripTopology* top_(dynamic_cast(&(roll->topology()))); + const float striplength(top_->stripLength()); + const float pitch(roll->pitch()); + auto nstrip = top_->nstrips(); +// std::cout << " GEM : " << rollName.c_str() << " / area : " << striplength*pitch*nstrip << std::endl; + hGEMArea_->GetXaxis()->SetBinLabel(iGEM, rollName.c_str()); + hGEMArea_->Fill(iGEM, striplength*pitch*nstrip); + } + + /* CSC */ + int iCSC = 0; + for (auto ly : cscGeom->layers()) { + ++iCSC; + auto detId = ly->id(); + auto endcap = std::to_string(detId.endcap()); // +1 : forward , -1 : backward + auto ch = std::to_string(detId.chamber()); + auto ring = std::to_string(detId.ring()); + auto station = std::to_string(detId.station()); + auto layer = std::to_string(detId.layer()); + const string name = "Endcap_"+endcap+"_station_"+station+"_ring_"+ring+"_CH_"+ch+"_layer_"+layer; + const RadialStripTopology* top_(dynamic_cast(&(ly->topology()))); + const float striplength(top_->stripLength()); + const float pitch(ly->geometry()->stripPitch()); + const float nstrip = top_->nstrips(); + //std::cout << " CSC : " << name << " / area : " << striplength*pitch*nstrip << std::endl; + hCSCArea_->GetXaxis()->SetBinLabel(iCSC, name.c_str()); + hCSCArea_->Fill(iCSC, striplength*pitch*nstrip); + } + + /* RPC */ + int iRPC = 0; + for ( const RPCRoll* roll : rpcGeom->rolls() ) { + ++iRPC; + const auto detId = roll->id(); + const string rollName = RPCGeomServ(detId).name(); + const double width = roll->surface().bounds().width(); + const double height = roll->surface().bounds().length(); + //std::cout << " RPC : " << rollName.c_str() << " / area : " << width*height << std::endl; + hRPCArea_->GetXaxis()->SetBinLabel(iRPC, rollName.c_str()); + hRPCArea_->Fill(iRPC, width*height); + } + +} +void MuonDetHitAnalyser::endRun(Run const&, EventSetup const&){} + +void MuonDetHitAnalyser::initValue() { + b_nME0Digis = 0; b_nME0RecHits = 0; + b_nCSC2DRecHits = 0; + b_nGEMDigis = 0; b_nGEMRecHits = 0; + b_nDTDigis = 0; b_nDTRecHits = 0; + b_nRPCDigis = 0; b_nRPCRecHits = 0; + + /*ME0 digi*/ + b_ME0_Digi_region = -9; b_ME0_Digi_chamber = -9; b_ME0_Digi_layer = -9; b_ME0_Digi_etaPartition = -9; b_ME0_Digi_bx = -999; + b_ME0_Digi_eta = -9; b_ME0_Digi_phi = -9; + /*ME0 rechit*/ + b_ME0_RecHit_region = -9; b_ME0_RecHit_chamber = -9; b_ME0_RecHit_layer = -9; b_ME0_RecHit_etaPartition = -9; + b_ME0_RecHit_eta = -9; b_ME0_RecHit_phi = -9; + + /*CSC rechit*/ + b_CSC_2DRecHit_region = -9; b_CSC_2DRecHit_station = -9; b_CSC_2DRecHit_ring = -9; b_CSC_2DRecHit_chamber = -9; b_CSC_2DRecHit_layer = -9; + b_CSC_2DRecHit_eta = -9; b_CSC_2DRecHit_phi = -9; + + /*GEM digi*/ + b_GEM_Digi_region = -9; b_GEM_Digi_station = -9; b_GEM_Digi_ring = -9; b_GEM_Digi_chamber = -9; b_GEM_Digi_layer = -9; b_GEM_Digi_etaPartition = -9; b_GEM_Digi_bx = -999; + b_GEM_Digi_eta = -9; b_GEM_Digi_phi = -9; + /*GEM rechit*/ + b_GEM_RecHit_region = -9; b_GEM_RecHit_station = -9; b_GEM_RecHit_ring = -9; b_GEM_RecHit_chamber = -9; b_GEM_RecHit_layer = -9; b_GEM_RecHit_etaPartition = -9; b_GEM_RecHit_bx = -999; + b_GEM_RecHit_eta = -9; b_GEM_RecHit_phi = -9;; + + /* DT digi*/ + b_DT_Digi_sector = -9; b_DT_Digi_station = -9; b_DT_Digi_wheel = -9; b_DT_Digi_chamber = -9; b_DT_Digi_superLayer = -9; b_DT_Digi_layer = -9; b_DT_Digi_wire = -9; + b_DT_Digi_eta = -9; b_DT_Digi_phi = -9; + /* DT rechit*/ + b_DT_RecHit_sector = -9; b_DT_RecHit_station = -9; b_DT_RecHit_wheel = -9; b_DT_RecHit_chamber = -9; b_DT_RecHit_superLayer = -9; b_DT_RecHit_layer = -9; b_DT_RecHit_wire = -9; + b_DT_RecHit_eta = -9; b_DT_RecHit_phi = -9; + + /*RPC digi*/ + b_RPC_Digi_region = -9; b_RPC_Digi_sector = -9; b_RPC_Digi_subSector = -9; b_RPC_Digi_station = -9; b_RPC_Digi_ring = -9; b_RPC_Digi_chamber = -9; b_RPC_Digi_layer = -9; b_RPC_Digi_roll = -9; + b_RPC_Digi_eta = -9; b_RPC_Digi_phi = -9; + b_RPC_Digi_isIRPC = -1; + /*RPC rechit*/ + b_RPC_RecHit_region = -9; b_RPC_RecHit_sector = -9; b_RPC_RecHit_subSector = -9; b_RPC_RecHit_station = -9; b_RPC_RecHit_ring = -9; b_RPC_RecHit_chamber = -9; b_RPC_RecHit_layer = -9; b_RPC_RecHit_roll = -9; + b_RPC_RecHit_eta = -9; b_RPC_RecHit_phi = -9; + b_RPC_RecHit_isIRPC = -1; +} + +//define this as a plug-in +DEFINE_FWK_MODULE(MuonDetHitAnalyser); diff --git a/MuonAnalyser/plugins/MuonSimAnalyser.cc b/MuonAnalyser/plugins/MuonSimAnalyser.cc new file mode 100644 index 00000000000..5e16ebeb870 --- /dev/null +++ b/MuonAnalyser/plugins/MuonSimAnalyser.cc @@ -0,0 +1,265 @@ +// cd /cms/ldap_home/iawatson/scratch/GEM/CMSSW_10_1_5/src/ && eval `scramv1 runtime -sh` && eval `scramv1 runtime -sh` && scram b -j 10 +// cd ../../.. && source /cvmfs/cms.cern.ch/cmsset_default.sh && eval `scramv1 runtime -sh` && eval `scramv1 runtime -sh` && scram b -j 10 +// system include files +#include +#include +#include +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +//#include "DQMServices/Core/interface/DQMStore.h" +//#include "DQMServices/Core/interface/DQMEDAnalyzer.h" +//#include "DQMServices/Core/interface/MonitorElement.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" +// ME0 +#include "DataFormats/MuonDetId/interface/ME0DetId.h" +#include "Geometry/GEMGeometry/interface/ME0Geometry.h" +#include "Geometry/GEMGeometry/interface/ME0EtaPartition.h" +#include "Geometry/GEMGeometry/interface/ME0EtaPartitionSpecs.h" +// GEM +#include "DataFormats/MuonDetId/interface/GEMDetId.h" +#include "Geometry/GEMGeometry/interface/GEMGeometry.h" +#include "Geometry/GEMGeometry/interface/GEMEtaPartition.h" +#include "Geometry/GEMGeometry/interface/GEMEtaPartitionSpecs.h" + +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" + +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" +#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" +#include "SimDataFormats/Track/interface/SimTrackContainer.h" +#include "SimDataFormats/Vertex/interface/SimVertexContainer.h" +#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" +#include "SimDataFormats/Associations/interface/MuonToTrackingParticleAssociator.h" +#include "SimDataFormats/Vertex/interface/SimVertex.h" +#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h" +#include "SimMuon/MCTruth/interface/MuonToSimAssociatorByHits.h" +#include "SimTracker/Common/interface/TrackingParticleSelector.h" +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h" + +#include "Geometry/Records/interface/MuonGeometryRecord.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" + +#include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h" +#include "Geometry/CommonTopologies/interface/RadialStripTopology.h" +#include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h" + +#include "Geometry/Records/interface/MuonGeometryRecord.h" + +#include "FWCore/Utilities/interface/EDGetToken.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/Run.h" + +#include "TH1D.h" +#include "TH2D.h" +#include "TProfile.h" +#include "TString.h" +#include "TGraphAsymmErrors.h" +#include "TLorentzVector.h" +#include "TTree.h" +#include "TDatabasePDG.h" + +using namespace std; +using namespace edm; + +class MuonSimAnalyser : public edm::EDAnalyzer { +public: + explicit MuonSimAnalyser(const edm::ParameterSet&); + ~MuonSimAnalyser(); + +private: + virtual void analyze(const edm::Event&, const edm::EventSetup&); + virtual void beginJob() override; + virtual void endJob() override; + + virtual void beginRun(Run const&, EventSetup const&) override; + virtual void endRun(Run const&, EventSetup const&) override; + + void initValue(); + + // ----------member data --------------------------- + edm::ParameterSet cfg_; + edm::EDGetToken ME0SimHitsToken_; + edm::EDGetToken GEMSimHitsToken_; + edm::EDGetToken simTracksToken_; + edm::EDGetToken simVerticesToken_; + + edm::Service fs; + + 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; + + /* 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; + +}; + +MuonSimAnalyser::MuonSimAnalyser(const edm::ParameterSet& iConfig) +{ + + ME0SimHitsToken_ = consumes(iConfig.getParameter("ME0SimInputLabel")); + GEMSimHitsToken_ = consumes(iConfig.getParameter("GEMSimInputLabel")); + simTracksToken_ = consumes< edm::SimTrackContainer >(iConfig.getParameter("simTrackCollection")); + simVerticesToken_ = consumes< edm::SimVertexContainer >(iConfig.getParameter("simVertexCollection")); + cfg_ = iConfig; + + + + t_event = fs->make("Event", "Event"); + t_event->Branch("nME0SimHits", &b_nME0SimHits, "nME0SimHits/I"); + t_event->Branch("nGEMSimHits", &b_nGEMSimHits, "nGEMSimHits/I"); + + /*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_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_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"); + t_ME0_simhit->Branch("SimHit_etaPartition", &b_ME0_SimHit_etaPartition, "SimHit_etaParition/I"); + + /*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_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_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"); + t_GEM_simhit->Branch("SimHit_chamber", &b_GEM_SimHit_chamber, "SimHit_chaber/I"); + t_GEM_simhit->Branch("SimHit_layer", &b_GEM_SimHit_layer, "SimHit_layer/I"); + t_GEM_simhit->Branch("SimHit_etaPartition", &b_GEM_SimHit_etaPartition, "SimHit_etaParition/I"); + +} + +MuonSimAnalyser::~MuonSimAnalyser() +{ +} + +void +MuonSimAnalyser::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + /* GEM Geometry */ + //edm::ESHandle hGEMGeom; + //iSetup.get().get(hGEMGeom); + //const GEMGeometry* GEMGeometry_ = &*hGEMGeom; + + edm::Handle ME0SimHits; + edm::Handle GEMSimHits; + edm::Handle simTracks; + edm::Handle simVertices; + iEvent.getByToken(ME0SimHitsToken_, ME0SimHits); + iEvent.getByToken(GEMSimHitsToken_, GEMSimHits); + iEvent.getByToken(simTracksToken_, simTracks); + iEvent.getByToken(simVerticesToken_, simVertices); + //if ( !simHits.isValid() || !simTracks.isValid() || !simVertices.isValid()) return; + + initValue(); + + /* 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(); + b_ME0_SimHit_pdgId = sh.particleType(); + + 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; + } 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_region = det_id.region(); + b_ME0_SimHit_chamber = det_id.chamber(); + b_ME0_SimHit_layer = det_id.layer(); + b_ME0_SimHit_etaPartition = det_id.roll(); + + ++b_nME0SimHits; + t_ME0_simhit->Fill(); + } + } else return; + + + /* GEM */ + 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(); + b_GEM_SimHit_pdgId = sh.particleType(); + + 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; + } 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_region = det_id.region(); + b_GEM_SimHit_station = det_id.station(); + b_GEM_SimHit_ring = det_id.ring(); + b_GEM_SimHit_chamber = det_id.chamber(); + b_GEM_SimHit_layer = det_id.layer(); + b_GEM_SimHit_etaPartition = det_id.roll(); + + ++b_nGEMSimHits; + t_GEM_simhit->Fill(); + } + } else return; + t_event->Fill(); +} + +void MuonSimAnalyser::beginJob(){} +void MuonSimAnalyser::endJob(){} + +void MuonSimAnalyser::beginRun(const edm::Run& run, const edm::EventSetup& iSetup){//(Run const& run, EventSetup const&){ +} +void MuonSimAnalyser::endRun(Run const&, EventSetup const&){} + +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; + + /*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; +} + +//define this as a plug-in +DEFINE_FWK_MODULE(MuonSimAnalyser); diff --git a/MuonAnalyser/plugins/MuonTrackAnalyser.cc b/MuonAnalyser/plugins/MuonTrackAnalyser.cc new file mode 100644 index 00000000000..a3bc90aaefc --- /dev/null +++ b/MuonAnalyser/plugins/MuonTrackAnalyser.cc @@ -0,0 +1,478 @@ +// cd /cms/ldap_home/iawatson/scratch/GEM/CMSSW_10_1_5/src/ && eval `scramv1 runtime -sh` && eval `scramv1 runtime -sh` && scram b -j 10 +// cd ../../.. && source /cvmfs/cms.cern.ch/cmsset_default.sh && eval `scramv1 runtime -sh` && eval `scramv1 runtime -sh` && scram b -j 10 +// system include files +#include +#include +#include +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" +// ME0 +#include "DataFormats/GEMDigi/interface/ME0DigiCollection.h" +#include "DataFormats/GEMRecHit/interface/ME0RecHitCollection.h" +#include "DataFormats/GEMRecHit/interface/ME0SegmentCollection.h" +#include "DataFormats/MuonDetId/interface/ME0DetId.h" +#include "Geometry/GEMGeometry/interface/ME0Geometry.h" +#include "Geometry/GEMGeometry/interface/ME0EtaPartition.h" +#include "Geometry/GEMGeometry/interface/ME0EtaPartitionSpecs.h" +// GEM +#include "DataFormats/GEMDigi/interface/GEMDigiCollection.h" +#include "DataFormats/GEMRecHit/interface/GEMRecHitCollection.h" +#include "DataFormats/GEMRecHit/interface/GEMSegmentCollection.h" +#include "DataFormats/MuonDetId/interface/GEMDetId.h" +#include "Geometry/GEMGeometry/interface/GEMGeometry.h" +#include "Geometry/GEMGeometry/interface/GEMEtaPartition.h" +#include "Geometry/GEMGeometry/interface/GEMEtaPartitionSpecs.h" +// CSC +#include "DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h" +#include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h" +#include "DataFormats/MuonDetId/interface/CSCDetId.h" +#include "Geometry/CSCGeometry/interface/CSCGeometry.h" +// DT +#include "DataFormats/DTDigi/interface/DTDigiCollection.h" +#include "DataFormats/DTRecHit/interface/DTRecHitCollection.h" +//#include "DataFormats/DTRecHit/interface/DTRecSegment2DCollection.h" +#include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h" +#include "DataFormats/MuonDetId/interface/DTBtiId.h" +#include "DataFormats/MuonDetId/interface/DTChamberId.h" +#include "DataFormats/MuonDetId/interface/DTLayerId.h" +#include "DataFormats/MuonDetId/interface/DTSectCollId.h" +#include "DataFormats/MuonDetId/interface/DTSuperLayerId.h" +#include "DataFormats/MuonDetId/interface/DTTracoId.h" +#include "DataFormats/MuonDetId/interface/DTWireId.h" +#include "Geometry/DTGeometry/interface/DTGeometry.h" +// RPC +#include "DataFormats/RPCDigi/interface/RPCDigiCollection.h" +#include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h" +#include "DataFormats/MuonDetId/interface/RPCDetId.h" +#include "Geometry/RPCGeometry/interface/RPCGeometry.h" +#include "Geometry/RPCGeometry/interface/RPCGeomServ.h" +// Muon +#include "DataFormats/PatCandidates/interface/Muon.h" +// L1 +#include "DataFormats/L1TMuon/interface/L1MuBMTrack.h" + +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" + +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" +#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" +#include "SimDataFormats/Associations/interface/MuonToTrackingParticleAssociator.h" +#include "SimDataFormats/Vertex/interface/SimVertex.h" +#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h" +#include "SimMuon/MCTruth/interface/MuonToSimAssociatorByHits.h" +#include "SimTracker/Common/interface/TrackingParticleSelector.h" +#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h" + +#include "Geometry/Records/interface/MuonGeometryRecord.h" +#include "Geometry/CommonDetUnit/interface/GeomDet.h" + +#include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h" +#include "Geometry/CommonTopologies/interface/RadialStripTopology.h" +#include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h" + +#include "Geometry/Records/interface/MuonGeometryRecord.h" + +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/Run.h" + +#include "TH1D.h" +#include "TH2D.h" +#include "TProfile.h" +#include "TString.h" +#include "TGraphAsymmErrors.h" +#include "TLorentzVector.h" +#include "TTree.h" + +using namespace std; +using namespace edm; + +class MuonTrackAnalyser : public edm::EDAnalyzer { +public: + explicit MuonTrackAnalyser(const edm::ParameterSet&); + ~MuonTrackAnalyser(); + +private: + virtual void analyze(const edm::Event&, const edm::EventSetup&); + virtual void beginJob() override; + virtual void endJob() override; + + virtual void beginRun(Run const&, EventSetup const&) override; + virtual void endRun(Run const&, EventSetup const&) override; + + void initMuonValue(); + void initValue(); + + // ----------member data --------------------------- + edm::EDGetTokenT me0Segments_; + + edm::EDGetTokenT cscSegments_; + + edm::EDGetTokenT gemSegments_; + + edm::EDGetTokenT dt4DSegments_; + + edm::EDGetTokenT simToken_; + edm::EDGetTokenT > muonToken_; + +// edm::EDGetTokenT > muonL1Token_; + + const reco::VertexCollection* vertexes_; + edm::EDGetTokenT > vtxToken_; + + edm::Service fs; + + TTree *t_event; + int b_run, b_lumi, b_latency; + int b_nME0Segments, b_nCSCSegments, b_nGEMSegments, b_nDT4DSegments; + + /* ME0 */ + TTree *t_ME0_seg; + int b_ME0_Seg_region, b_ME0_Seg_chamber, b_ME0_Seg_layer, b_ME0_Seg_nRecHits; + float b_ME0_Seg_eta, b_ME0_Seg_phi; + + /* CSC */ + TTree *t_CSC_seg; + int b_CSC_Seg_region, b_CSC_Seg_station, b_CSC_Seg_ring, b_CSC_Seg_chamber, b_CSC_Seg_layer, b_CSC_Seg_nRecHits; + float b_CSC_Seg_eta, b_CSC_Seg_phi; + + /* GEM */ + TTree *t_GEM_seg; + int b_GEM_Seg_region, b_GEM_Seg_station, b_GEM_Seg_ring, b_GEM_Seg_chamber, b_GEM_Seg_layer, b_GEM_Seg_nRecHits; + float b_GEM_Seg_eta, b_GEM_Seg_phi; + + /* DT */ + TTree *t_DT_seg; + int b_DT_4DSeg_sector, b_DT_4DSeg_station, b_DT_4DSeg_wheel, b_DT_4DSeg_chamber, b_DT_4DSeg_nRecHits, b_DT_4DSeg_hasZed; + float b_DT_4DSeg_eta, b_DT_4DSeg_phi; + + /* 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; + + 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; +}; + +MuonTrackAnalyser::MuonTrackAnalyser(const edm::ParameterSet& iConfig) +{ + me0Segments_ = consumes(iConfig.getParameter("me0Segments")); + + cscSegments_ = consumes(iConfig.getParameter("cscSegments")); + + gemSegments_ = consumes(iConfig.getParameter("gemSegments")); + + dt4DSegments_ = consumes(iConfig.getParameter("dt4DSegments")); + + simToken_ = consumes(iConfig.getParameter("simLabel")); + muonToken_ = consumes >(iConfig.getParameter("muonLabel")); + +// muonL1Token_ = consumes >(iConfig.getParameter("muonL1Label")); + + vtxToken_ = consumes >(iConfig.getParameter("primaryVertex")); + + t_event = fs->make("Event", "Event"); + t_event->Branch("nME0Segments", &b_nME0Segments, "nME0Segments/I"); + t_event->Branch("nCSCSegments", &b_nCSCSegments, "nCSCSegments/I"); + t_event->Branch("nGEMSegments", &b_nGEMSegments, "nGEMSegments/I"); + t_event->Branch("nDT4DSegments", &b_nDT4DSegments, "nDT4DSegments/I"); + + /*ME0*/ + t_ME0_seg = fs->make("ME0_Segment", "ME0_Segment"); + t_ME0_seg->Branch("Seg_eta", &b_ME0_Seg_eta, "Seg_eta/F"); + t_ME0_seg->Branch("Seg_phi", &b_ME0_Seg_phi, "Seg_phi/F"); + t_ME0_seg->Branch("Seg_nRecHits", &b_ME0_Seg_nRecHits, "Seg_nRecHits/I"); + t_ME0_seg->Branch("Seg_region", &b_ME0_Seg_region, "Seg_region/I"); + t_ME0_seg->Branch("Seg_chamber", &b_ME0_Seg_chamber, "Seg_chamber/I"); + t_ME0_seg->Branch("Seg_layer", &b_ME0_Seg_layer, "Seg_layer/I"); + + /*CSC*/ + t_CSC_seg = fs->make("CSC_Segment", "CSC_Segment"); + t_CSC_seg->Branch("Seg_eta", &b_CSC_Seg_eta, "Seg_eta/F"); + t_CSC_seg->Branch("Seg_phi", &b_CSC_Seg_phi, "Seg_phi/F"); + t_CSC_seg->Branch("Seg_nRecHits", &b_CSC_Seg_nRecHits, "Seg_nRecHits/I"); + t_CSC_seg->Branch("Seg_region", &b_CSC_Seg_region, "Seg_region/I"); + t_CSC_seg->Branch("Seg_station", &b_CSC_Seg_station, "Seg_station/I"); + t_CSC_seg->Branch("Seg_ring", &b_CSC_Seg_ring, "Seg_ring/I"); + t_CSC_seg->Branch("Seg_chamber", &b_CSC_Seg_chamber, "Seg_chamber/I"); + t_CSC_seg->Branch("Seg_layer", &b_CSC_Seg_layer, "Seg_layer/I"); + + /*GEM*/ + t_GEM_seg = fs->make("GEM_Segment", "GEM_Segment"); + t_GEM_seg->Branch("Seg_eta", &b_GEM_Seg_eta, "Seg_eta/F"); + t_GEM_seg->Branch("Seg_phi", &b_GEM_Seg_phi, "Seg_phi/F"); + t_GEM_seg->Branch("Seg_nRecHits", &b_GEM_Seg_nRecHits, "Seg_nRecHits/I"); + t_GEM_seg->Branch("Seg_region", &b_GEM_Seg_region, "Seg_region/I"); + t_GEM_seg->Branch("Seg_station", &b_GEM_Seg_station, "Seg_station/I"); + t_GEM_seg->Branch("Seg_ring", &b_GEM_Seg_ring, "Seg_ring/I"); + t_GEM_seg->Branch("Seg_chamber", &b_GEM_Seg_chamber, "Seg_chamber/I"); + t_GEM_seg->Branch("Seg_layer", &b_GEM_Seg_layer, "Seg_layer/I"); + + /*DT*/ + t_DT_seg = fs->make("DT_4DSegment", "DT_4DSegment"); + t_DT_seg->Branch("Seg_eta", &b_DT_4DSeg_eta, "Seg_eta/F"); + t_DT_seg->Branch("Seg_phi", &b_DT_4DSeg_phi, "Seg_phi/F"); + t_DT_seg->Branch("Seg_nRecHits", &b_DT_4DSeg_nRecHits, "Seg_nRecHits/I"); + t_DT_seg->Branch("Seg_sector", &b_DT_4DSeg_sector, "Seg_sector/I"); + t_DT_seg->Branch("Seg_station", &b_DT_4DSeg_station, "Seg_station/I"); + t_DT_seg->Branch("Seg_wheel", &b_DT_4DSeg_wheel, "Seg_wheel/I"); + t_DT_seg->Branch("Seg_chamber", &b_DT_4DSeg_chamber, "Seg_chamber/I"); + t_DT_seg->Branch("Seg_hasZed", &b_DT_4DSeg_hasZed, "Seg_hasZed/I"); + + /*Muon*/ + t_genMuon = fs->make("gen_Muon", "gen_Muon"); + t_genMuon->Branch("pt", &b_genMuon_pt, "pt/F"); + t_genMuon->Branch("eta", &b_genMuon_eta, "eta/F"); + t_genMuon->Branch("phi", &b_genMuon_phi, "phi/F"); + t_genMuon->Branch("nMatchedStationLayer", &b_genMuon_nMatchedStationLayer, "nMatchedStationLayer/I"); + 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("isTight", &b_genMuon_isTight, "isTight/I"); + t_genMuon->Branch("isMedium", &b_genMuon_isMedium, "isMedium/I"); + t_genMuon->Branch("isLoose", &b_genMuon_isLoose, "isLoose/I"); + + t_Muon = fs->make("Muon", "Muon"); + t_Muon->Branch("pt", &b_muon_pt, "pt/F"); + t_Muon->Branch("eta", &b_muon_eta, "eta/F"); + t_Muon->Branch("phi", &b_muon_phi, "phi/F"); + t_Muon->Branch("nMatchedStationLayer", &b_muon_nMatchedStationLayer, "nMatchedStationLayer/I"); + 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("isTight", &b_muon_isTight, "isTight/I"); + t_Muon->Branch("isMedium", &b_muon_isMedium, "isMedium/I"); + t_Muon->Branch("isLoose", &b_muon_isLoose, "isLoose/I"); + +} + +MuonTrackAnalyser::~MuonTrackAnalyser() +{ +} + +void +MuonTrackAnalyser::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + edm::ESHandle hME0Geom; + iSetup.get().get(hME0Geom); + const ME0Geometry* ME0Geometry_ = &*hME0Geom; + + /* ME0 Geometry */ + edm::Handle me0Segments; + iEvent.getByToken(me0Segments_, me0Segments); + + /* CSC Geometry */ + edm::ESHandle hCSCGeom; + iSetup.get().get(hCSCGeom); + const CSCGeometry* CSCGeometry_ = &*hCSCGeom; + + edm::Handle cscSegments; + iEvent.getByToken(cscSegments_, cscSegments); + + /* GEM Geometry */ + edm::ESHandle hGEMGeom; + iSetup.get().get(hGEMGeom); + const GEMGeometry* GEMGeometry_ = &*hGEMGeom; + + edm::Handle gemSegments; + iEvent.getByToken(gemSegments_, gemSegments); + + /* DT Geometry */ + edm::ESHandle hDTGeom; + iSetup.get().get(hDTGeom); + const DTGeometry* DTGeometry_ = &*hDTGeom; + + edm::Handle dt4DSegments; + iEvent.getByToken(dt4DSegments_, dt4DSegments); + + /* Muon */ + Handle simHandle; + iEvent.getByToken(simToken_, simHandle); + + edm::Handle > muonHandle; + iEvent.getByToken(muonToken_, muonHandle); + +// edm::Handle > muonL1Handle; +// iEvent.getByToken(muonL1Token_, muonL1Handle); + + edm::Handle vertices; + iEvent.getByToken(vtxToken_, vertices); + + vertexes_ = vertices.product(); + reco::Vertex pv0 = vertexes_->at(0); + + initValue(); + + for (auto ch : ME0Geometry_->chambers()) { + /* ME0 seg */ + ME0DetId cId = ch->id(); + auto segsRange = me0Segments->get(cId); + auto me0Seg = segsRange.first; + for (auto seg = me0Seg; seg != segsRange.second; ++seg) { + auto segLd = seg->localPosition(); + auto gp = ch->toGlobal(segLd); + b_ME0_Seg_eta = gp.eta(); + b_ME0_Seg_phi = gp.phi(); + b_ME0_Seg_nRecHits = seg->nRecHits(); + b_ME0_Seg_region = cId.region(); + b_ME0_Seg_chamber = cId.chamber(); + b_ME0_Seg_layer = cId.layer(); + //std::cout << b_ME0_Seg_chamber << " ==> seg x : " << segLd.x() << " seg y : " << segLd.y() << " seg eta : " << gp.eta() << " nRecHits : " << b_ME0_Seg_nRecHits << std::endl; + b_nME0Segments++; + t_ME0_seg->Fill(); + } + } + + /* CSC */ + for (auto ch : CSCGeometry_->chambers()) { + // CSC seg + CSCDetId cId = ch->id(); + auto segsRange = cscSegments->get(cId); + auto cscSeg = segsRange.first; + for (auto seg = cscSeg; seg != segsRange.second; ++seg) { + auto segLd = seg->localPosition(); + auto gp = ch->toGlobal(segLd); + b_CSC_Seg_eta = gp.eta(); + b_CSC_Seg_phi = gp.phi(); + b_CSC_Seg_nRecHits = seg->nRecHits(); + b_CSC_Seg_region = cId.endcap(); + b_CSC_Seg_station = cId.station(); + b_CSC_Seg_ring = cId.ring(); + b_CSC_Seg_chamber = cId.chamber(); + b_CSC_Seg_layer = cId.layer(); + b_nCSCSegments++; + t_CSC_seg->Fill(); + } + } + + /* GEM */ + for (auto ch : GEMGeometry_->chambers()) { + // GEM + GEMDetId cId = ch->id(); + auto segsRange = gemSegments->get(cId); + auto gemSeg = segsRange.first; + for (auto seg = gemSeg; seg != segsRange.second; ++seg) { + auto segLd = seg->localPosition(); + auto gp = ch->toGlobal(segLd); + b_GEM_Seg_eta = gp.eta(); + b_GEM_Seg_phi = gp.phi(); + b_GEM_Seg_nRecHits = seg->nRecHits(); + b_GEM_Seg_region = cId.region(); + b_GEM_Seg_station = cId.station(); + b_GEM_Seg_ring = cId.ring(); + b_GEM_Seg_chamber = cId.chamber(); + b_GEM_Seg_layer = cId.layer(); + b_nGEMSegments++; + t_GEM_seg->Fill(); + } + } + + /* DT */ + for (auto ch : DTGeometry_->chambers()) { + // DT seg + DTChamberId cId = ch->id(); + //b_DT_4DSeg_chamber = ch->id(); + auto segsRange = dt4DSegments->get(cId); + auto dt4DSeg = segsRange.first; + for (auto seg = dt4DSeg; seg != segsRange.second; ++seg) { + auto segLd = seg->localPosition(); + auto gp = ch->toGlobal(segLd); + b_DT_4DSeg_eta = gp.eta(); + b_DT_4DSeg_phi = gp.phi(); + b_DT_4DSeg_nRecHits = seg->recHits().size(); + b_DT_4DSeg_sector = cId.sector(); + b_DT_4DSeg_station = cId.station(); + b_DT_4DSeg_wheel = cId.wheel(); + //b_DT_4DSeg_chamber = cId.chamber(); + b_DT_4DSeg_hasZed = seg->hasZed(); + b_nDT4DSegments++; + t_DT_seg->Fill(); + } + } + + /* gen Muon */ + for (TrackingParticleCollection::size_type i=0; isize(); i++) { + TrackingParticleRef simRef(simHandle, i); + const TrackingParticle* simTP = simRef.get(); + b_genMuon_pt = simTP->pt(); + b_genMuon_eta = simTP->eta(); + b_genMuon_phi = simTP->phi(); + bool isSignalMuon = abs(simTP->pdgId())==13 && !simTP->genParticles().empty() && (simTP->eventId().event() == 0) && (simTP->eventId().bunchCrossing() == 0); + if (!isSignalMuon) continue; + t_genMuon->Fill(); + } + /* reco Muon */ + for(size_t i = 0; i < muonHandle->size(); ++i) { + initMuonValue(); + edm::RefToBase muRef = muonHandle->refAt(i); + const reco::Muon* mu = muRef.get(); + b_muon_isTight = muon::isTightMuon(*mu, pv0); + b_muon_isMedium = muon::isMediumMuon(*mu); + b_muon_isLoose = muon::isLooseMuon(*mu); + b_muon_pt = mu->pt(); + b_muon_eta = mu->eta(); + b_muon_phi = mu->phi(); + b_muon_isGlobal = mu->isGlobalMuon(); + b_muon_isTracker = mu->isTrackerMuon(); + b_muon_isME0 = mu->isME0Muon(); + b_muon_nMatchedStationLayer = mu->numberOfMatchedStations(); + t_Muon->Fill(); + } + /*L1 trigger*/ + + t_event->Fill(); +} + +void MuonTrackAnalyser::beginJob(){} +void MuonTrackAnalyser::endJob(){} + +void MuonTrackAnalyser::beginRun(const edm::Run& run, const edm::EventSetup& iSetup){//(Run const& run, EventSetup const&){ +} +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_pt = -9; b_muon_eta = -9; b_muon_phi = -9; + b_muon_nMatchedStationLayer = -9; + +} +void MuonTrackAnalyser::initValue() { + b_nME0Segments = 0; + b_nCSCSegments = 0; + b_nGEMSegments = 0; + b_nDT4DSegments = 0; + + /*ME0 seg*/ + b_ME0_Seg_region = -9; b_ME0_Seg_chamber = -9; b_ME0_Seg_layer = -9; b_ME0_Seg_nRecHits = -9; + b_ME0_Seg_eta = -9; b_ME0_Seg_phi = -9; + + /*CSC seg*/ + b_CSC_Seg_region = -9; b_CSC_Seg_station = -9; b_CSC_Seg_ring = -9; b_CSC_Seg_chamber = -9; b_CSC_Seg_layer = -9; b_CSC_Seg_nRecHits = -9; + b_CSC_Seg_eta = -9; b_CSC_Seg_phi = -9; + + /*GEM seg*/ + b_GEM_Seg_region = -9; b_GEM_Seg_station = -9; b_GEM_Seg_ring = -9; b_GEM_Seg_chamber = -9; b_GEM_Seg_layer = -9; b_GEM_Seg_nRecHits = -9; + b_GEM_Seg_eta = -9; b_GEM_Seg_phi = -9; + + /* DT seg*/ + b_DT_4DSeg_sector = -9; b_DT_4DSeg_station = -9; b_DT_4DSeg_wheel = -9; b_DT_4DSeg_chamber = -9; b_DT_4DSeg_nRecHits = -9; + b_DT_4DSeg_eta = -9; b_DT_4DSeg_phi = -9; + b_DT_4DSeg_hasZed = -1; +} + +//define this as a plug-in +DEFINE_FWK_MODULE(MuonTrackAnalyser); diff --git a/MuonAnalyser/plugins/SliceTestAnalysis.cc b/MuonAnalyser/plugins/SliceTestAnalysis.cc deleted file mode 100644 index fa3dc613349..00000000000 --- a/MuonAnalyser/plugins/SliceTestAnalysis.cc +++ /dev/null @@ -1,545 +0,0 @@ -// cd /cms/ldap_home/iawatson/scratch/GEM/CMSSW_10_1_5/src/ && eval `scramv1 runtime -sh` && eval `scramv1 runtime -sh` && scram b -j 10 -// cd ../../.. && source /cvmfs/cms.cern.ch/cmsset_default.sh && eval `scramv1 runtime -sh` && eval `scramv1 runtime -sh` && scram b -j 10 -// system include files -#include -#include -#include -#include -#include - -// user include files -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/EDAnalyzer.h" - -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/MakerMacros.h" - -#include "FWCore/ParameterSet/interface/ParameterSet.h" -#include "FWCore/ServiceRegistry/interface/Service.h" -#include "CommonTools/UtilAlgos/interface/TFileService.h" - -#include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h" -#include "TrackingTools/GeomPropagators/interface/Propagator.h" -#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" -#include "TrackingTools/Records/interface/TransientTrackRecord.h" -#include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h" -#include "MagneticField/Engine/interface/MagneticField.h" - -#include "DataFormats/PatCandidates/interface/Muon.h" -#include "DataFormats/VertexReco/interface/Vertex.h" -#include "DataFormats/VertexReco/interface/VertexFwd.h" -#include "DataFormats/GEMRecHit/interface/GEMRecHitCollection.h" -#include "DataFormats/MuonDetId/interface/GEMDetId.h" -#include "Geometry/GEMGeometry/interface/GEMGeometry.h" -#include "Geometry/GEMGeometry/interface/GEMEtaPartition.h" -#include "Geometry/GEMGeometry/interface/GEMEtaPartitionSpecs.h" -#include "Geometry/Records/interface/MuonGeometryRecord.h" -#include "Geometry/CommonDetUnit/interface/GeomDet.h" -#include "Geometry/CommonTopologies/interface/StripTopology.h" - -#include "Geometry/Records/interface/MuonGeometryRecord.h" - -#include "FWCore/Framework/interface/ESHandle.h" -#include "FWCore/Framework/interface/Run.h" - -#include "TH1D.h" -#include "TH2D.h" -#include "TString.h" -#include "TGraphAsymmErrors.h" -#include "TLorentzVector.h" -#include "TTree.h" - -using namespace std; -using namespace edm; - -class SliceTestAnalysis : public edm::EDAnalyzer { -public: - explicit SliceTestAnalysis(const edm::ParameterSet&); - ~SliceTestAnalysis(); - -private: - virtual void analyze(const edm::Event&, const edm::EventSetup&); - virtual void beginJob() override; - virtual void endJob() override; - - virtual void beginRun(Run const&, EventSetup const&) override; - virtual void endRun(Run const&, EventSetup const&) override; - - // ----------member data --------------------------- - edm::EDGetTokenT gemRecHits_; - edm::EDGetTokenT > muons_; - edm::EDGetTokenT vertexCollection_; - edm::Service fs; - - MuonServiceProxy* theService_; - edm::ESHandle propagator_; - edm::ESHandle ttrackBuilder_; - edm::ESHandle bField_; - - TH2D* h_firstStrip[36][3]; - TH2D* h_allStrips[36][3]; - TH2D* h_globalPosOnGem; - TH1D* h_clusterSize, *h_totalStrips, *h_bxtotal; - TH1D* h_inEta[36][3]; - TH1D* h_hitEta[36][3]; - TH1D* h_trkEta[36][3]; - - TH1D* h_res_x, *h_res_y, *h_pull_x, *h_pull_y; - - TTree *t_hit; - int b_run, b_lumi, b_event; - int b_firstStrip, b_nStrips, b_chamber, b_layer, b_etaPartition, b_muonQuality, b_bx; - float b_x, b_y, b_z; - float b_mu_eta, b_mu_phi, b_mu_pt; - float b_pull_x, b_pull_y, b_res_x, b_res_y; - - int nEvents, nMuonTotal, nGEMFiducialMuon, nGEMTrackWithMuon; - int b_nMuons, b_nMuonsWithGEMHit; - int b_valid; - - int b_nGEMHits; - - // muon branches - int m_nhits, m_nvalidhits; - int m_nbounds; - int m_quality; - float m_pt, m_eta, m_phi; - vector m_roll, m_chamber, m_layer, m_strip; // hit info - vector m_resx, m_resy, m_pullx, m_pully; - vector m_in_roll, m_in_chamber, m_in_layer, m_in_strip; // propagation bound info - vector m_in_globalPhi, m_in_globalEta, m_in_nearGemPhi, m_in_nearGemEta, m_in_nearGemRoll, m_in_nearGemFirstStrip, m_in_nearGemClusterSize; // global info - - TTree *t_run; - TTree *t_muon; - TTree *t_event; -}; - -struct GEMMuonAssociation { - int muonQuality; - float mu_phi; - float mu_eta; - float mu_pt; - - float pull_x; - float pull_y; - float res_x; - float res_y; - - int valid; -}; - -SliceTestAnalysis::SliceTestAnalysis(const edm::ParameterSet& iConfig) : - nEvents(0), - nMuonTotal(0), - nGEMFiducialMuon(0), - nGEMTrackWithMuon(0) -{ - gemRecHits_ = consumes(iConfig.getParameter("gemRecHits")); - muons_ = consumes >(iConfig.getParameter("muons")); - vertexCollection_ = consumes(iConfig.getParameter("vertexCollection")); - edm::ParameterSet serviceParameters = iConfig.getParameter("ServiceParameters"); - theService_ = new MuonServiceProxy(serviceParameters); - - h_clusterSize=fs->make(Form("clusterSize"),"clusterSize",100,0,100); - h_totalStrips=fs->make(Form("totalStrips"),"totalStrips",200,0,200); - h_bxtotal=fs->make(Form("bx"),"bx",31,-15,15); - - h_globalPosOnGem = fs->make(Form("onGEM"), "onGEM", 100, -100, 100, 100, -100, 100); - - t_run = fs->make("Run", "Run"); - t_run->Branch("run", &b_run, "run/I"); - - t_event = fs->make("Event", "Event"); - t_event->Branch("nMuons", &b_nMuons, "nMuons/I"); - t_event->Branch("nMuonsWithGEMHit", &b_nMuonsWithGEMHit, "nMuonsWithGEMHit/I"); - t_event->Branch("nGEMHits", &b_nGEMHits, "nGEMHits/I"); - t_event->Branch("run", &b_run, "run/I"); - t_event->Branch("lumi", &b_lumi, "lumi/I"); - - h_res_x=fs->make(Form("res_x"),"res_x",100,-50,50); - h_res_y=fs->make(Form("res_y"),"res_y",100,-50,50); - h_pull_x=fs->make(Form("pull_x"),"pull_x",100,-50,50); - h_pull_y=fs->make(Form("pull_y"),"pull_y",100,-50,50); - - t_muon = fs->make("Muon", "Muon"); - t_muon->Branch("run", &b_run, "run/I"); - t_muon->Branch("lumi", &b_lumi, "lumi/I"); - t_muon->Branch("nhits", &m_nhits, "nhits/I")->SetTitle("n GEM hits associated to muon"); - t_muon->Branch("nvalidhits", &m_nvalidhits, "nvalidhits/I")->SetTitle("n GEM hits associated to muon, and muon can propagate to eta partition of hit"); - t_muon->Branch("nbounds", &m_nbounds, "nbounds/I")->SetTitle("times muon is in GEM eta partition bounds"); - t_muon->Branch("pt", &m_pt, "pt/F"); - t_muon->Branch("eta", &m_eta, "eta/F"); - t_muon->Branch("phi", &m_phi, "phi/F"); - t_muon->Branch("quality", &m_quality, "quality/I")->SetTitle("muon quality :: 0:noid 1:looseID 2:tightID"); - t_muon->Branch("roll", &m_roll); - t_muon->Branch("chamber", &m_chamber); - t_muon->Branch("layer", &m_layer); - t_muon->Branch("strip", &m_strip); - t_muon->Branch("resx", &m_resx); - t_muon->Branch("resy", &m_resy); - t_muon->Branch("pullx", &m_pullx); - t_muon->Branch("pully", &m_pully); - t_muon->Branch("in_roll", &m_in_roll); - t_muon->Branch("in_chamber", &m_in_chamber); - t_muon->Branch("in_layer", &m_in_layer); - t_muon->Branch("in_strip", &m_in_strip); - t_muon->Branch("in_globalPhi", &m_in_globalPhi); - t_muon->Branch("in_globalEta", &m_in_globalEta); - t_muon->Branch("in_nearGemPhi", &m_in_nearGemPhi); - t_muon->Branch("in_nearGemEta", &m_in_nearGemEta); - t_muon->Branch("in_nearGemRoll", &m_in_nearGemRoll); - t_muon->Branch("in_nearGemFirstStrip", &m_in_nearGemFirstStrip); - t_muon->Branch("in_nearGemClusterSize", &m_in_nearGemClusterSize); - - t_hit = fs->make("Hit", "Hit"); - t_hit->Branch("run", &b_run, "run/I"); - t_hit->Branch("lumi", &b_lumi, "lumi/I"); - - t_hit->Branch("bx", &b_bx, "bx/I"); - t_hit->Branch("firstStrip", &b_firstStrip, "firstStrip/I"); - t_hit->Branch("nStrips", &b_nStrips, "nStrips/I"); - t_hit->Branch("chamber", &b_chamber, "chamber/I"); - t_hit->Branch("layer", &b_layer, "layer/I"); - t_hit->Branch("etaPartition", &b_etaPartition, "etaPartition/I"); - t_hit->Branch("muonQuality", &b_muonQuality, "muonQuality/I")->SetTitle("muonQuality -1:none 0:noid 1:looseID 2:tightID"); - t_hit->Branch("x", &b_x, "x/F"); - t_hit->Branch("y", &b_y, "y/F"); - t_hit->Branch("z", &b_z, "z/F"); - t_hit->Branch("pull_x", &b_pull_x, "pull_x/F"); - t_hit->Branch("pull_y", &b_pull_y, "pull_y/F"); - t_hit->Branch("res_x", &b_res_x, "res_x/F"); - t_hit->Branch("res_y", &b_res_y, "res_y/F"); - t_hit->Branch("valid", &b_valid, "valid/I"); - t_hit->Branch("mu_phi", &b_mu_phi, "mu_phi/F"); - t_hit->Branch("mu_eta", &b_mu_eta, "mu_eta/F"); - t_hit->Branch("mu_pt", &b_mu_pt, "mu_pt/F"); - - for (int ichamber=0; ichamber<36;++ichamber) { - // for (int ichamber=27; ichamber<=30;++ichamber) { - for (int ilayer=1; ilayer<3;++ilayer) { - h_firstStrip[ichamber][ilayer] = fs->make(Form("firstStrip ch %i lay %i",ichamber, ilayer),"firstStrip",384,1,385,8,0.5,8.5); - h_firstStrip[ichamber][ilayer]->GetXaxis()->SetTitle("strip"); - h_firstStrip[ichamber][ilayer]->GetYaxis()->SetTitle("iEta"); - - h_allStrips[ichamber][ilayer] = fs->make(Form("allStrips ch %i lay %i",ichamber, ilayer),"allStrips",384,1,385,8,0.5,8.5); - h_allStrips[ichamber][ilayer]->GetXaxis()->SetTitle("strip"); - h_allStrips[ichamber][ilayer]->GetYaxis()->SetTitle("iEta"); - - h_inEta[ichamber][ilayer] = fs->make(Form("inEta ch %i lay %i",ichamber, ilayer),"inEta",8,0.5,8.5); - h_hitEta[ichamber][ilayer] = fs->make(Form("hitEta ch %i lay %i",ichamber, ilayer),"hitEta",8,0.5,8.5); - h_trkEta[ichamber][ilayer] = fs->make(Form("trkEta ch %i lay %i",ichamber, ilayer),"trkEta",8,0.5,8.5); - } - } -} - -SliceTestAnalysis::~SliceTestAnalysis() -{ - std::cout << "::: GEM Slice Test Results :::" << std::endl; - std::cout << ": From " << nEvents << " events" << std::endl; - std::cout << std::endl; - std::cout << " # Muons " << nMuonTotal << std::endl; - std::cout << " # FidMu " << nGEMFiducialMuon << std::endl; - std::cout << " # GEMMu " << nGEMTrackWithMuon << std::endl; - std::cout << std::endl; -} - -void -SliceTestAnalysis::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) -{ - nEvents++; - - b_nMuons = 0; - b_nMuonsWithGEMHit = 0; - b_nGEMHits = 0; - - b_run = iEvent.run(); - b_lumi = iEvent.luminosityBlock(); - - edm::ESHandle hGeom; - iSetup.get().get(hGeom); - const GEMGeometry* GEMGeometry_ = &*hGeom; - - iSetup.get().get("TransientTrackBuilder",ttrackBuilder_); - // iSetup.get().get("SteppingHelixPropagatorAny",propagator_); - // iSetup.get().get(bField_); - theService_->update(iSetup); - auto propagator = theService_->propagator("SteppingHelixPropagatorAny"); - - edm::Handle gemRecHits; - iEvent.getByToken(gemRecHits_, gemRecHits); - - edm::Handle vertexCollection; - iEvent.getByToken( vertexCollection_, vertexCollection ); - if(vertexCollection.isValid()) { - vertexCollection->size(); - // std::cout << "vertex->size() " << vertexCollection->size() < > muons; - iEvent.getByToken(muons_, muons); - - std::map muHits; - - for (size_t i = 0; i < muons->size(); ++i) { - b_nMuons++; - nMuonTotal++; - - m_nhits = 0; - m_nvalidhits = 0; - m_nbounds = 0; - m_roll.clear(); m_chamber.clear(); m_layer.clear(); m_strip.clear(); - m_resx.clear(); m_resy.clear(); m_pullx.clear(); m_pully.clear(); - m_in_roll.clear(); m_in_chamber.clear(); m_in_layer.clear(); m_in_strip.clear(); - m_in_globalPhi.clear(); m_in_globalEta.clear(); m_in_nearGemPhi.clear(); m_in_nearGemEta.clear(); m_in_nearGemFirstStrip.clear(); - - edm::RefToBase muRef = muons->refAt(i); - const reco::Muon* mu = muRef.get(); - - if (mu->passed(reco::Muon::Selector::CutBasedIdTight)) - m_quality = 2; - else if (mu->passed(reco::Muon::Selector::CutBasedIdLoose)) - m_quality = 1; - else - m_quality = 0; - - const reco::Track* muonTrack = 0; - if ( mu->globalTrack().isNonnull() ) muonTrack = mu->globalTrack().get(); - else if ( mu->outerTrack().isNonnull() ) muonTrack = mu->outerTrack().get(); - if (muonTrack) { - - std::set detLists; - - reco::TransientTrack ttTrack = ttrackBuilder_->build(muonTrack); - bool onDet = false; - - // prepropagate to near the GEM region, to speedup per/etapart prop. w/o loss of generatlity - // TrajectoryStateOnSurface tsos_ = propagator->propagate(ttTrack.outermostMeasurementState(), - // GEMGeometry_->etaPartitions()[0]->surface()); - // if (!tsos_.isValid()) continue; - - // GlobalPoint tsosGP = tsos.globalPosition(); - - for (auto ch : GEMGeometry_->etaPartitions()) { - TrajectoryStateOnSurface tsos = propagator->propagate(ttTrack.outermostMeasurementState(), - ch->surface()); - if (!tsos.isValid()) continue; - - GlobalPoint tsosGP = tsos.globalPosition(); - const LocalPoint pos = ch->toLocal(tsosGP); - const LocalPoint pos2D(pos.x(), pos.y(), 0); - const BoundPlane& bps(ch->surface()); - h_globalPosOnGem->Fill(tsosGP.x(), tsosGP.y()); - - if (bps.bounds().inside(pos2D)) { - m_nbounds++; - onDet = true; - auto gemid = ch->id(); - h_inEta[gemid.chamber()][gemid.layer()]->Fill(gemid.roll()); - m_in_roll.push_back(gemid.roll()); - m_in_chamber.push_back(gemid.chamber()); - m_in_layer.push_back(gemid.layer()); - m_in_strip.push_back((int)(ch->strip(pos))); - - m_in_globalPhi.push_back(tsosGP.phi()); - m_in_globalEta.push_back(tsosGP.eta()); - - float gemEta = +99.0, gemPhi = +99.0; - int gemRoll = -9, gemFirstStrip = -9, gemClusterSize = -9; - - - for (auto ch : GEMGeometry_->chambers()) { - for(auto roll : ch->etaPartitions()) { - GEMDetId rId = roll->id(); - if (rId != gemid) continue; - - auto recHitsRange = gemRecHits->get(rId); - auto gemRecHit = recHitsRange.first; - for (auto hit = gemRecHit; hit != recHitsRange.second; ++hit) { - auto gemGlob = ch->toGlobal(hit->localPosition()); - - // pick closest hit - if (fabs(gemGlob.phi() - tsosGP.phi()) < gemPhi) { - gemPhi = gemGlob.phi(); - gemEta = gemGlob.eta(); - gemRoll = rId.roll(); - gemFirstStrip = hit->firstClusterStrip(); - gemClusterSize = hit->clusterSize(); - } - - // h_hitEta[rId.chamber()][rId.layer()]->Fill(rId.roll()); - // h_firstStrip[rId.chamber()][rId.layer()]->Fill(hit->firstClusterStrip(), rId.roll()); - // h_clusterSize->Fill(hit->clusterSize()); - // h_bxtotal->Fill(hit->BunchX()); - } - } - } - m_in_nearGemPhi.push_back(gemPhi); - m_in_nearGemEta.push_back(gemEta); - m_in_nearGemRoll.push_back(gemRoll); - m_in_nearGemFirstStrip.push_back(gemFirstStrip); - m_in_nearGemClusterSize.push_back(gemClusterSize); - if (gemFirstStrip != -9) { - cout << " strip: " << gemFirstStrip << " VFAT: " << int(gemFirstStrip/128) << " X: " << ch->centreOfStrip(gemFirstStrip).x() << " phi: " << ch->toGlobal(ch->centreOfStrip(gemFirstStrip)).phi()<< endl; - cout << gemPhi << " " << tsosGP.phi() << endl; - } - - for (auto hit = muonTrack->recHitsBegin(); hit != muonTrack->recHitsEnd(); hit++) { - if ((*hit)->geographicalId().det() == 2 && (*hit)->geographicalId().subdetId() == 4) { - if ((*hit)->rawId() == ch->id().rawId() ) { - GEMDetId gemid((*hit)->geographicalId()); - auto etaPart = GEMGeometry_->etaPartition(gemid); - // cout << "found it "<< gemid - // << " lp " << (*hit)->localPosition() - // << " gp " << etaPart->toGlobal((*hit)->localPosition()) - // << endl; - } - } - } - } - } - - for (auto hit = muonTrack->recHitsBegin(); hit != muonTrack->recHitsEnd(); hit++) { - if ( (*hit)->geographicalId().det() == 2 && (*hit)->geographicalId().subdetId() == 4) { - GEMDetId gemid((*hit)->geographicalId()); - h_trkEta[gemid.chamber()][gemid.layer()]->Fill(gemid.roll()); - } - } - - if (onDet) ++nGEMFiducialMuon; - - if (muonTrack->hitPattern().numberOfValidMuonGEMHits()) { - ++b_nMuonsWithGEMHit; - ++nGEMTrackWithMuon; - for (auto hit = muonTrack->recHitsBegin(); hit != muonTrack->recHitsEnd(); hit++) { - if ( (*hit)->geographicalId().det() == 2 && (*hit)->geographicalId().subdetId() == 4) { - GEMDetId gemid((*hit)->geographicalId()); - auto etaPart = GEMGeometry_->etaPartition(gemid); - - m_nhits++; - - TrajectoryStateOnSurface tsos = propagator->propagate(ttTrack.outermostMeasurementState(),etaPart->surface()); - if (!tsos.isValid()) continue; - // GlobalPoint tsosGP = tsos.globalPosition(); - - m_nvalidhits++; - - LocalPoint && tsos_localpos = tsos.localPosition(); - LocalError && tsos_localerr = tsos.localError().positionError(); - LocalPoint && dethit_localpos = (*hit)->localPosition(); - LocalError && dethit_localerr = (*hit)->localPositionError(); - auto res_x = (dethit_localpos.x() - tsos_localpos.x()); - auto res_y = (dethit_localpos.y() - tsos_localpos.y()); - auto pull_x = (dethit_localpos.x() - tsos_localpos.x()) / - std::sqrt(dethit_localerr.xx() + tsos_localerr.xx()); - auto pull_y = (dethit_localpos.y() - tsos_localpos.y()) / - std::sqrt(dethit_localerr.yy() + tsos_localerr.yy()); - - h_res_x->Fill(res_x); - h_res_y->Fill(res_y); - h_pull_x->Fill(pull_x); - h_pull_y->Fill(pull_y); - - m_roll.push_back(gemid.roll()); - m_chamber.push_back(gemid.chamber()); - m_layer.push_back(gemid.layer()); - auto gemhit = static_cast(*hit); - m_strip.push_back((int) (etaPart->strip(dethit_localpos))); - - m_resx.push_back(res_x); - m_resy.push_back(res_y); - m_pullx.push_back(pull_x); - m_pully.push_back(pull_y); - - int isvalid = (*hit)->isValid(); - auto mup = tsos.globalMomentum(); - muHits[static_cast(*hit)] = {muonQuality: m_quality, - mu_phi: mup.phi(), mu_eta: mup.eta(), mu_pt: mup.perp(), - pull_x: pull_x, pull_y: pull_y, - res_x: res_x, res_y: res_y, valid: isvalid}; - } - } - } - } - - if (m_nhits > 0 or m_nbounds > 0) { - m_pt = mu->pt(); - m_eta = mu->eta(); - m_phi = mu->phi(); - t_muon->Fill(); - } - } // Muon Loop - - int totalStrips = 0; - - for (auto ch : GEMGeometry_->chambers()) { - for(auto roll : ch->etaPartitions()) { - GEMDetId rId = roll->id(); - - auto recHitsRange = gemRecHits->get(rId); - auto gemRecHit = recHitsRange.first; - for (auto hit = gemRecHit; hit != recHitsRange.second; ++hit) { - - h_hitEta[rId.chamber()][rId.layer()]->Fill(rId.roll()); - h_firstStrip[rId.chamber()][rId.layer()]->Fill(hit->firstClusterStrip(), rId.roll()); - h_clusterSize->Fill(hit->clusterSize()); - h_bxtotal->Fill(hit->BunchX()); - for (int nstrip = hit->firstClusterStrip(); nstrip < hit->firstClusterStrip()+hit->clusterSize(); ++nstrip) { - totalStrips++; - h_allStrips[rId.chamber()][rId.layer()]->Fill(nstrip, rId.roll()); - } - - b_firstStrip = hit->firstClusterStrip(); - b_bx = hit->BunchX(); - b_nStrips = hit->clusterSize(); - b_chamber = rId.chamber(); - b_layer = rId.layer(); - b_etaPartition = rId.roll(); - b_muonQuality = -1; - b_pull_x = -99; - b_pull_y = -99; - b_res_x = -99; - b_res_y = -99; - b_valid = -1; - b_mu_phi = b_mu_eta = b_mu_pt = -99; - for (const auto & kv : muHits) { - auto muHit = kv.first; - if (*hit == *muHit) { - auto assoc = kv.second; - b_muonQuality = assoc.muonQuality; - b_pull_x = assoc.pull_x; - b_pull_y = assoc.pull_y; - b_res_x = assoc.res_x; - b_res_y = assoc.res_y; - b_valid = assoc.valid; - b_mu_eta = assoc.mu_eta; - b_mu_phi = assoc.mu_phi; - b_mu_pt = assoc.mu_pt; - } - } - - auto globalPosition = roll->toGlobal(hit->localPosition()); - b_x = globalPosition.x(); - b_y = globalPosition.y(); - b_z = globalPosition.z(); - - t_hit->Fill(); - b_nGEMHits++; - } - } - } - h_totalStrips->Fill(totalStrips); - - t_event->Fill(); -} - -void SliceTestAnalysis::beginJob(){} -void SliceTestAnalysis::endJob(){} - -void SliceTestAnalysis::beginRun(Run const& run, EventSetup const&){ - b_run = run.run(); - t_run->Fill(); -} -void SliceTestAnalysis::endRun(Run const&, EventSetup const&){} - -//define this as a plug-in -DEFINE_FWK_MODULE(SliceTestAnalysis); diff --git a/MuonAnalyser/test/runHGCalSimTest.py b/MuonAnalyser/test/runHGCalSimTest.py new file mode 100644 index 00000000000..6230a238945 --- /dev/null +++ b/MuonAnalyser/test/runHGCalSimTest.py @@ -0,0 +1,74 @@ +import FWCore.ParameterSet.Config as cms +from Configuration.StandardSequences.Eras import eras + +process = cms.Process('HGCalSimTest',eras.Run2_2017,eras.run3_GEM) + +process.load("FWCore.MessageService.MessageLogger_cfi") +process.load('Configuration.Geometry.GeometryExtended2023D23Reco_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') +process.load('RecoLocalMuon.GEMRecHit.gemLocalReco_cff') +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, '101X_dataRun2_Prompt_v10', '') +process.MessageLogger.cerr.FwkReport.reportEvery = 5000 + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(-1) +) +#process.maxEvents.input = cms.untracked.int32(10) +# Input source +process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring()) +process.source.skipEvents = cms.untracked.uint32(0) + +#process.source.fileNames.append('/store/user/yekang/CRAB_PrivateMC/HGCalStainless_default_me0_singleNu_RECO/190129_065738/0000/singleNu_GEN-SIM-DIGI_10.root') +#process.source.fileNames.append('/store/user/yekang/CRAB_PrivateMC/HGCalStainless_default_MinBias/190213_140302/0000/step1_10.root') +#process.source.fileNames.append('root://uosaf0007.sscc.uos.ac.kr:1094//xrootd/store/user/yekang/CRAB_PrivateMC/HGCalStainless_default_me0_diMu_RECO/190225_172244/0000/diMu_GEN-SIM-DIGI_951.root') +process.source.fileNames.append('root://uosaf0007.sscc.uos.ac.kr:1094//xrootd/store/user/yekang/CRAB_PrivateMC/HGCalStainless_modified_me0_diMu_RECO/190225_172318/0000/diMu_GEN-SIM-DIGI_489.root') +#from glob import glob +#process.source.fileNames.append('file:/cms/ldap_home/yckang/me0/CMSSW_10_4_0/src/MuonPerformance/MuonAnalyser/test/tenMu_GEN-SIM-DIGI.root') + +#fname = 'singleMuon.txt' +#f = open(fname) +#for line in f: +# process.source.fileNames.append(line) + +process.options = cms.untracked.PSet( +) + +process.TFileService = cms.Service("TFileService",fileName = cms.string("histo.root")) +process.MuonDetHitAnalyser = cms.EDAnalyzer('MuonDetHitAnalyser', + me0Digis = cms.InputTag("simMuonME0Digis"), + me0RecHits = cms.InputTag("me0RecHits"), + csc2DRecHits = cms.InputTag("csc2DRecHits"), + gemDigis = cms.InputTag("simMuonGEMDigis"), + gemRecHits = cms.InputTag("gemRecHits"), + dtDigis = cms.InputTag("simMuonDTDigis"), + dtRecHits = cms.InputTag("dt1DRecHits"), + rpcDigis = cms.InputTag("simMuonRPCDigis"), + rpcRecHits = cms.InputTag("rpcRecHits"), +) + +process.MuonTrackAnalyser = cms.EDAnalyzer('MuonTrackAnalyser', + me0Segments = cms.InputTag("me0Segments"), + cscSegments = cms.InputTag("cscSegments"), + gemSegments = cms.InputTag("gemSegments"), + dt4DSegments = cms.InputTag("dt4DSegments"), + simLabel = cms.InputTag("mix","MergedTrackTruth"), + muonLabel = cms.InputTag("muons"), + primaryVertex = cms.InputTag('offlinePrimaryVertices'), +) + +process.MuonSimAnalyser = cms.EDAnalyzer('MuonSimAnalyser', +# verboseSimHit = cms.untracked.int32(1), + ME0SimInputLabel = cms.InputTag('g4SimHits', "MuonME0Hits"), + GEMSimInputLabel = cms.InputTag('g4SimHits', "MuonGEMHits"), +# simMuOnlyGEM = cms.untracked.bool(True), + #verboseSimHit = cms.untracked.int32(1), + #simInputLabel = cms.InputTag('g4SimHits',"MuonGEMHits"), + simTrackCollection = cms.InputTag('g4SimHits'), + simVertexCollection = cms.InputTag('g4SimHits') +) + +process.p = cms.Path(process.MuonDetHitAnalyser + process.MuonTrackAnalyser + process.MuonSimAnalyser) +#process.p = cms.Path(process.MuonSimAnalyser) + +process.schedule = cms.Schedule(process.p) diff --git a/MuonAnalyser/test/runHGCalSimTest_modified.py b/MuonAnalyser/test/runHGCalSimTest_modified.py new file mode 100644 index 00000000000..914c7c9494e --- /dev/null +++ b/MuonAnalyser/test/runHGCalSimTest_modified.py @@ -0,0 +1,57 @@ +import FWCore.ParameterSet.Config as cms +from Configuration.StandardSequences.Eras import eras + +process = cms.Process('HGCalSimTest',eras.Run2_2017,eras.run3_GEM) + +process.load("FWCore.MessageService.MessageLogger_cfi") +process.load('Configuration.Geometry.GeometryExtended2023D23Reco_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') +process.load('RecoLocalMuon.GEMRecHit.me0LocalReco_cff') +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, '101X_dataRun2_Prompt_v10', '') +process.MessageLogger.cerr.FwkReport.reportEvery = 5000 + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(-1) +) +#process.maxEvents.input = cms.untracked.int32(10) +# Input source +process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring()) +process.source.skipEvents = cms.untracked.uint32(0) + +process.source.fileNames.append('/store/user/yekang/me0/tenMu_modified/tenMu_GEN-SIM-DIGI_060.root') +#from glob import glob +#process.source.fileNames.append('file:/cms/ldap_home/yckang/me0/CMSSW_10_4_0/src/MuonPerformance/MuonAnalyser/test/tenMu_GEN-SIM-DIGI.root') + +#fname = 'singleMuon.txt' +#f = open(fname) +#for line in f: +# process.source.fileNames.append(line) + +process.options = cms.untracked.PSet( +) + +process.TFileService = cms.Service("TFileService",fileName = cms.string("histo.root")) + +process.reconstruction_step = cms.Path(process.me0LocalReco) +process.HGCalSimTest = cms.EDAnalyzer('HGCalSimTest', + me0Digis = cms.InputTag("simMuonME0Digis"), + me0Segments = cms.InputTag("me0Segments"), + me0RecHits = cms.InputTag("me0RecHits"), + cscSegments = cms.InputTag("cscSegments"), + csc2DRecHits = cms.InputTag("csc2DRecHits"), + gemDigis = cms.InputTag("simMuonGEMDigis"), + gemSegments = cms.InputTag("gemSegments"), + gemRecHits = cms.InputTag("gemRecHits"), + dtDigis = cms.InputTag("simMuonDTDigis"), + dt4DSegments = cms.InputTag("dt4DSegments"), + dtRecHits = cms.InputTag("dt1DRecHits"), + rpcDigis = cms.InputTag("simMuonRPCDigis"), + rpcRecHits = cms.InputTag("rpcRecHits"), + muonLabel = cms.InputTag("muons"), + primaryVertex = cms.InputTag('offlinePrimaryVertices'), + +) +process.p = cms.Path(process.HGCalSimTest) + +process.schedule = cms.Schedule(process.reconstruction_step, process.p) diff --git a/MuonAnalyser/test/runSliceTestAnalysis.py b/MuonAnalyser/test/runSliceTestAnalysis.py deleted file mode 100644 index 7c6c61668ef..00000000000 --- a/MuonAnalyser/test/runSliceTestAnalysis.py +++ /dev/null @@ -1,44 +0,0 @@ -import FWCore.ParameterSet.Config as cms -from Configuration.StandardSequences.Eras import eras - -process = cms.Process('SliceTestAnalysis',eras.Run2_2017,eras.run3_GEM) - -process.load("FWCore.MessageService.MessageLogger_cfi") -process.load('Configuration.StandardSequences.GeometryRecoDB_cff') -process.load('Configuration.StandardSequences.MagneticField_AutoFromDBCurrent_cff') -process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') -process.load('RecoMuon.TrackingTools.MuonServiceProxy_cff') -process.load('TrackingTools.TransientTrack.TransientTrackBuilder_cfi') -from Configuration.AlCa.GlobalTag import GlobalTag -process.GlobalTag = GlobalTag(process.GlobalTag, '101X_dataRun2_Prompt_v10', '') -process.MessageLogger.cerr.FwkReport.reportEvery = 5000 - -process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(-1) -) -#process.maxEvents.input = cms.untracked.int32(10) -# Input source -process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring()) -process.source.skipEvents = cms.untracked.uint32(0) - -#process.source.fileNames.append('/store/data/Run2018B/Cosmics/AOD/PromptReco-v1/000/317/428/00000/E4BC1D7B-3F6A-E811-9E05-FA163E57A064.root') -from glob import glob -process.source.fileNames.extend(['file:'+f for f in glob('/xrootd/store/user/jlee/SingleMuon/Run2017G-v1/RECOv2/step3*.root')][:5]) -#process.source.fileNames.append('/store/user/jlee/SingleMuon/Run2017G-v1/FEVTEvent/step3_313.root') - -#fname = 'singleMuon.txt' -#f = open(fname) -#for line in f: -# process.source.fileNames.append(line) - -process.options = cms.untracked.PSet() - -process.TFileService = cms.Service("TFileService",fileName = cms.string("histo.root")) - -process.SliceTestAnalysis = cms.EDAnalyzer('SliceTestAnalysis', - process.MuonServiceProxy, - gemRecHits = cms.InputTag("gemRecHits"), - muons = cms.InputTag("muons"), - vertexCollection = cms.InputTag("offlinePrimaryVertices"), -) -process.p = cms.Path(process.SliceTestAnalysis) diff --git a/MuonAnalyser/test/runSliceTestAnalysis2018.py b/MuonAnalyser/test/runSliceTestAnalysis2018.py deleted file mode 100644 index 5bfbdd1dbdc..00000000000 --- a/MuonAnalyser/test/runSliceTestAnalysis2018.py +++ /dev/null @@ -1,61 +0,0 @@ -import FWCore.ParameterSet.Config as cms -from Configuration.StandardSequences.Eras import eras - -process = cms.Process('SliceTestAnalysis',eras.Run2_2018,eras.run3_GEM) - -process.load("FWCore.MessageService.MessageLogger_cfi") -process.load('Configuration.StandardSequences.GeometryRecoDB_cff') -process.load('Configuration.StandardSequences.MagneticField_AutoFromDBCurrent_cff') -process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') -process.load('RecoMuon.TrackingTools.MuonServiceProxy_cff') -process.load('TrackingTools.TransientTrack.TransientTrackBuilder_cfi') -from Configuration.AlCa.GlobalTag import GlobalTag -process.GlobalTag = GlobalTag(process.GlobalTag, '101X_dataRun2_Prompt_v11', '') -process.MessageLogger.cerr.FwkReport.reportEvery = 5000 - -process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(-1) -) -#process.maxEvents.input = cms.untracked.int32(10) -# Input source -process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring()) -process.source.skipEvents = cms.untracked.uint32(0) - -### Files available as at 17.7.2018 -# /xrootd/store/user/jlee/SingleMuon/Run2018C-v1/RECOv1/step3_000.root 319337 -# /xrootd/store/user/jlee/SingleMuon/Run2018C-v1/RECOv1/step3_349.root 319347 -# /xrootd/store/user/jlee/SingleMuon/Run2018C-v1/RECOv1/step3_975.root 319348 -# /xrootd/store/user/jlee/SingleMuon/Run2018C-v1/RECOv1/step3_993.root 319349 -# /xrootd/store/user/jlee/SingleMuon/Run2018C-v1/RECOv1/step3_1140.root 319449 -# /xrootd/store/user/jlee/SingleMuon/Run2018C-v1/RECOv1/step3_1510.root 319450 -# /xrootd/store/user/jlee/SingleMuon/Run2018C-v1/RECOv1/step3_1518.root 319456 -# Error in : file /xrootd/store/user/jlee/SingleMuon/Run2018C-v1/RECOv1/step3_1808.root does not exist -# Skipping /xrootd/store/user/jlee/SingleMuon/Run2018C-v1/RECOv1/step3_1808.root -# /xrootd/store/user/jlee/SingleMuon/Run2018C-v1/RECOv1/step3_1857.root 319459 - - -#process.source.fileNames.append('/store/data/Run2018B/Cosmics/AOD/PromptReco-v1/000/317/428/00000/E4BC1D7B-3F6A-E811-9E05-FA163E57A064.root') -from glob import glob -process.source.fileNames.extend( - ['/store/user/jlee/SingleMuon/Run2018C-v1/RECOv2/step3_019.root'] - # ['file:/xrootd/store/user/iawatson/SingleMuon/Run2018C-v1/wGEM/RECOv1/step3_{0:03d}.root'.format(i) for i in range(682, 682+1)] - # ['file:'+f for f in glob('/xrootd/store/user/jlee/SingleMuon/Run2018C-v1/RECOv1/step3*.root')][:] -) -#process.source.fileNames.append('/store/user/jlee/SingleMuon/Run2017G-v1/FEVTEvent/step3_313.root') - -#fname = 'singleMuon.txt' -#f = open(fname) -#for line in f: -# process.source.fileNames.append(line) - -process.options = cms.untracked.PSet() - -process.TFileService = cms.Service("TFileService",fileName = cms.string("histo2018.root")) - -process.SliceTestAnalysis = cms.EDAnalyzer('SliceTestAnalysis', - process.MuonServiceProxy, - gemRecHits = cms.InputTag("gemRecHits"), - muons = cms.InputTag("muons"), - vertexCollection = cms.InputTag("offlinePrimaryVertices"), -) -process.p = cms.Path(process.SliceTestAnalysis) diff --git a/MuonAnalyser/test/step3_RAW2DIGI_L1Reco_RECO.py b/MuonAnalyser/test/step3_RAW2DIGI_L1Reco_RECO.py index 933721f86b6..b0368a661bb 100644 --- a/MuonAnalyser/test/step3_RAW2DIGI_L1Reco_RECO.py +++ b/MuonAnalyser/test/step3_RAW2DIGI_L1Reco_RECO.py @@ -2,12 +2,12 @@ # using: # Revision: 1.19 # Source: /local/reps/CMSSW/CMSSW/Configuration/Applications/python/ConfigBuilder.py,v -# with command line options: step3 --conditions 101X_dataRun2_Prompt_v10 -s RAW2DIGI,L1Reco,RECO --runUnscheduled --process reRECO --data --era Run2_2017 --eventcontent RECO --hltProcess reHLT --scenario pp --datatier RECO --customise Configuration/DataProcessing/RecoTLR.customisePostEra_Run2_2017 -n -1 --filein /store/data/Run2017G/SingleMuon/RAW/v1/000/306/801/00001/7ADC8664-3DCE-E711-ADEF-02163E019BF4.root --fileout file:step3.root +# with command line options: step3 --conditions 101X_dataRun2_Prompt_v10 -s RAW2DIGI,L1Reco,RECO --runUnscheduled --process reRECO --data --era Run2_2018 --eventcontent RECO --hltProcess reHLT --scenario pp --datatier RECO --customise Configuration/DataProcessing/RecoTLR.customisePostEra_Run2_2018 -n -1 --filein /store/data/Run2017G/SingleMuon/RAW/v1/000/306/801/00001/7ADC8664-3DCE-E711-ADEF-02163E019BF4.root --fileout file:step3.root import FWCore.ParameterSet.Config as cms from Configuration.StandardSequences.Eras import eras -process = cms.Process('reRECO',eras.Run2_2017,eras.run3_GEM) +process = cms.Process('reRECO',eras.Run2_2018,eras.run3_GEM) # import of standard configurations process.load('Configuration.StandardSequences.Services_cff') @@ -26,7 +26,7 @@ # Input source process.source = cms.Source("PoolSource", - fileNames = cms.untracked.vstring('file:/xrootd/store/data/Run2017G/SingleMuon/RAW/v1/000/306/826/00000/FE2B6447-A4CE-E711-8DD7-02163E019CD7.root'), + fileNames = cms.untracked.vstring('file:/xrootd/store/data/Run2018C/SingleMuon/RAW/v1/000/319/337/00000/EEE24158-BE82-E811-8419-FA163E00EF28.root'), secondaryFileNames = cms.untracked.vstring() ) #import FWCore.PythonUtilities.LumiList as LumiList @@ -58,9 +58,9 @@ # Additional output definition -# Other statements +# Other statements 101X_dataRun2_Prompt_v11 from Configuration.AlCa.GlobalTag import GlobalTag -process.GlobalTag = GlobalTag(process.GlobalTag, '101X_dataRun2_Prompt_v10', '') +process.GlobalTag = GlobalTag(process.GlobalTag, '101X_dataRun2_Prompt_v11', '') # Path and EndPath definitions process.raw2digi_step = cms.Path(process.RawToDigi) @@ -77,10 +77,10 @@ # customisation of the process. # Automatic addition of the customisation function from Configuration.DataProcessing.RecoTLR -from Configuration.DataProcessing.RecoTLR import customisePostEra_Run2_2017 +from Configuration.DataProcessing.RecoTLR import customisePostEra_Run2_2018 -#call to customisation function customisePostEra_Run2_2017 imported from Configuration.DataProcessing.RecoTLR -process = customisePostEra_Run2_2017(process) +#call to customisation function customisePostEra_Run2_2018 imported from Configuration.DataProcessing.RecoTLR +process = customisePostEra_Run2_2018(process) # End of customisation functions #do not add changes to your config after this point (unless you know what you are doing)