GEMCSCSegment/test/GEMCSCCoincidenceRateAnalyzer.cc

0001 #ifndef RecoLocalMuon_GEMCSCSegment_GEMCSCCoincidenceRateAnalyzer_h
0002 #define RecoLocalMuon_GEMCSCSegment_GEMCSCCoincidenceRateAnalyzer_h
0003 /** \class GEMCSCCoincidenceRateAnalyzer
0004  *
0005  * \author Seungjin Yang <seungjin.yang@cern.ch>
0006  */
0007 #include <algorithm>
0008 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014 #include "FWCore/Utilities/interface/InputTag.h"
0015 #include "DataFormats/GEMRecHit/interface/GEMCSCSegmentCollection.h"
0016 #include "DataFormats/GEMDigi/interface/GEMVFATStatusCollection.h"
0017 #include "DataFormats/GEMDigi/interface/GEMOHStatusCollection.h"
0018 #include "DataFormats/MuonReco/interface/Muon.h"
0019 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0020 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0021 #include "DataFormats/MuonDetId/interface/GEMDetId.h"
0022 #include "DataFormats/MuonReco/interface/MuonChamberMatch.h"
0023 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0024
0025 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0026 #include "FWCore/ServiceRegistry/interface/Service.h"
0027
0028 #include "TTree.h"
0029
0030 class GEMCSCCoincidenceRateAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0031 public:
0032   explicit GEMCSCCoincidenceRateAnalyzer(const edm::ParameterSet&);
0033   ~GEMCSCCoincidenceRateAnalyzer() override = default;
0034   static void fillDescriptions(edm::ConfigurationDescriptions&);
0035
0036 private:
0037   void analyze(const edm::Event&, const edm::EventSetup&) override;
0038
0039   void resetBranches();
0040   bool analyzeGEMCSCSegment(const GEMCSCSegment&,
0041                             const GEMOHStatusCollection*,
0042                             const GEMVFATStatusCollection*,
0043                             const std::vector<const CSCSegment*>&);
0044
0045   bool checkGEMChamberStatus(const GEMDetId& chamber_id, const GEMOHStatusCollection*, const GEMVFATStatusCollection*);
0046   std::vector<const CSCSegment*> findMuonSegments(const reco::MuonCollection*);
0047   bool isMuonSegment(const CSCSegment&, const std::vector<const CSCSegment*>);
0048   bool checkCSCChamberType(const CSCDetId&);
0049
0050   const edm::EDGetTokenT<GEMCSCSegmentCollection> gem_csc_segment_collection_token_;
0051   const edm::EDGetTokenT<GEMOHStatusCollection> gem_oh_status_collection_token_;
0052   const edm::EDGetTokenT<GEMVFATStatusCollection> gem_vfat_status_collection_token_;
0053   const edm::EDGetTokenT<reco::MuonCollection> muon_collection_token_;
0054   const std::string log_category_;
0055   const bool use_gem_daq_status_;
0056   const std::vector<unsigned int> csc_whitelist_;
0057
0058   //
0059   TTree* tree_;
0060   //// GEMCSCSegment
0061   int b_region_;
0062   int b_station_;
0063   int b_ring_;
0064   int b_chamber_;
0065   bool b_gem_chamber_has_error_;
0066   float b_norm_chi2_;
0067   //// CSCSegment
0068   float b_csc_norm_chi2_;
0069   int b_csc_num_hit_;
0070   bool b_csc_is_muon_;
0071   //// GEMRecHit
0072   int b_gem_num_hit;
0073   std::vector<int> b_gem_layer_;
0074   std::vector<int> b_gem_ieta_;
0075   std::vector<int> b_gem_first_strip_;
0076   std::vector<int> b_gem_cls_;
0077 };
0078
0079 GEMCSCCoincidenceRateAnalyzer::GEMCSCCoincidenceRateAnalyzer(const edm::ParameterSet& ps)
0080     : gem_csc_segment_collection_token_(
0081           consumes<GEMCSCSegmentCollection>(ps.getUntrackedParameter<edm::InputTag>("gemcscSegmentTag"))),
0082       gem_oh_status_collection_token_(
0083           consumes<GEMOHStatusCollection>(ps.getUntrackedParameter<edm::InputTag>("ohStatusTag"))),
0084       gem_vfat_status_collection_token_(
0085           consumes<GEMVFATStatusCollection>(ps.getUntrackedParameter<edm::InputTag>("vfatStatusTag"))),
0086       muon_collection_token_(consumes<reco::MuonCollection>(ps.getUntrackedParameter<edm::InputTag>("muonTag"))),
0087       log_category_(ps.getUntrackedParameter<std::string>("logCategory")),
0088       use_gem_daq_status_(ps.getUntrackedParameter<bool>("useGEMDAQStatus")),
0089       csc_whitelist_(ps.getUntrackedParameter<std::vector<unsigned int> >("cscWhitelist")) {
0090   usesResource(TFileService::kSharedResource);
0091
0092   edm::Service<TFileService> file_service;
0093   tree_ = file_service->make<TTree>("gemcsc", "gemcsc");
0094
0095   tree_->Branch("region", &b_region_);
0096   tree_->Branch("station", &b_station_);
0097   tree_->Branch("ring", &b_ring_);
0098   tree_->Branch("chamber", &b_chamber_);
0099   tree_->Branch("gem_chamber_has_error", &b_gem_chamber_has_error_, "gem_chamber_has_error/O");
0100
0101   tree_->Branch("norm_chi2", &b_norm_chi2_);
0102   tree_->Branch("csc_norm_chi2", &b_csc_norm_chi2_);
0103   tree_->Branch("csc_num_hit", &b_csc_num_hit_);
0104   tree_->Branch("csc_is_muon", &b_csc_is_muon_);
0105
0106   tree_->Branch("gem_num_hit", &b_gem_num_hit);
0107   tree_->Branch("gem_layer", &b_gem_layer_);
0108   tree_->Branch("gem_ieta", &b_gem_ieta_);
0109   tree_->Branch("gem_first_strip", &b_gem_first_strip_);
0110   tree_->Branch("gem_cls", &b_gem_cls_);
0111 }
0112
0113 void GEMCSCCoincidenceRateAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0114   edm::ParameterSetDescription desc;
0115   desc.addUntracked<edm::InputTag>("ohStatusTag", edm::InputTag("muonGEMDigis", "OHStatus"));
0116   desc.addUntracked<edm::InputTag>("vfatStatusTag", edm::InputTag("muonGEMDigis", "VFATStatus"));
0117   desc.addUntracked<std::string>("logCategory", "GEMCSCCoincidenceRateAnalyzer");
0118   desc.addUntracked<edm::InputTag>("gemcscSegmentTag", edm::InputTag("gemcscSegments"));
0119   desc.addUntracked<edm::InputTag>("muonTag", edm::InputTag("muons"));
0120   desc.addUntracked<bool>("useGEMDAQStatus", false);
0121   desc.addUntracked<std::vector<unsigned int> >("cscWhitelist",
0122                                                 {CSCDetId::iChamberType(1, 1), CSCDetId::iChamberType(2, 1)});
0123   descriptions.addWithDefaultLabel(desc);
0124 }
0125
0126 void GEMCSCCoincidenceRateAnalyzer::analyze(const edm::Event& event, const edm::EventSetup&) {
0127   /////////////////////////////////////////////////////////////////////////////
0128   // get data from edm::Event
0129   /////////////////////////////////////////////////////////////////////////////
0130   const GEMCSCSegmentCollection* gem_csc_segment_collection = nullptr;
0131   if (const edm::Handle<GEMCSCSegmentCollection> handle = event.getHandle(gem_csc_segment_collection_token_)) {
0132     gem_csc_segment_collection = handle.product();
0133
0134   } else {
0135     edm::LogError(log_category_) << "failed to get GEMCSCSegmentCollection";
0136     return;
0137   }
0138
0139   const GEMOHStatusCollection* gem_oh_status_collection = nullptr;
0140   const GEMVFATStatusCollection* gem_vfat_status_collection = nullptr;
0141   if (use_gem_daq_status_) {
0142     if (auto handle = event.getHandle(gem_oh_status_collection_token_)) {
0143       gem_oh_status_collection = handle.product();
0144
0145     } else {
0146       edm::LogError(log_category_) << "failed to get OHVFATStatusCollection";
0147       return;
0148     }
0149
0150     if (auto handle = event.getHandle(gem_vfat_status_collection_token_)) {
0151       gem_vfat_status_collection = handle.product();
0152
0153     } else {
0154       edm::LogError(log_category_) << "failed to get GEMVFATStatusCollection";
0155       return;
0156     }
0157   }
0158
0159   const reco::MuonCollection* muon_collection = nullptr;
0160   if (const edm::Handle<reco::MuonCollection> handle = event.getHandle(muon_collection_token_)) {
0161     muon_collection = handle.product();
0162
0163   } else {
0164     edm::LogError(log_category_) << "failed to get reco::MuonCollection";
0165     return;
0166   }
0167
0168   /////////////////////////////////////////////////////////////////////////////
0169   // analyze
0170   /////////////////////////////////////////////////////////////////////////////
0171   const std::vector<const CSCSegment*> muon_segment_vec = findMuonSegments(muon_collection);
0172
0173   for (edm::OwnVector<GEMCSCSegment>::const_iterator iter = gem_csc_segment_collection->begin();
0174        iter != gem_csc_segment_collection->end();
0175        iter++) {
0176     const GEMCSCSegment& gem_csc_segment = *iter;
0177     if (not checkCSCChamberType(gem_csc_segment.cscDetId())) {
0178       continue;
0179     }
0180     resetBranches();
0181     if (analyzeGEMCSCSegment(gem_csc_segment, gem_oh_status_collection, gem_vfat_status_collection, muon_segment_vec)) {
0182       tree_->Fill();
0183     }
0184   }  // GEMCSCSegment
0185 }
0186
0187 void GEMCSCCoincidenceRateAnalyzer::resetBranches() {
0188   // detector
0189   b_region_ = 0;
0190   b_station_ = 0;
0191   b_ring_ = 0;
0192   b_chamber_ = 0;
0193   b_gem_chamber_has_error_ = false;
0194   // GEMCSCSegment
0195   b_norm_chi2_ = 0.0f;
0196   // CSCSegment
0197   b_csc_norm_chi2_ = 0.0f;
0198   b_csc_num_hit_ = 0;
0199   b_csc_is_muon_ = false;
0200   // GEMRecHit
0201   b_gem_num_hit = 0;
0202   b_gem_layer_.clear();
0203   b_gem_ieta_.clear();
0204   b_gem_first_strip_.clear();
0205   b_gem_cls_.clear();
0206 }
0207
0208 // Main function to analyzer GEMCSCSegment.
0209 // Returns a boolean indicating whether the GEMCSCSegment analysis is successful.
0210 bool GEMCSCCoincidenceRateAnalyzer::analyzeGEMCSCSegment(const GEMCSCSegment& gem_csc_segment,
0211                                                          const GEMOHStatusCollection* gem_oh_status_collection,
0212                                                          const GEMVFATStatusCollection* gem_vfat_status_collection,
0213                                                          const std::vector<const CSCSegment*>& muon_segment_vec) {
0214   const CSCDetId csc_id = gem_csc_segment.cscDetId();
0215   const CSCSegment csc_segment = gem_csc_segment.cscSegment();
0216   const std::vector<GEMRecHit>& gem_hit_vec = gem_csc_segment.gemRecHits();
0217   // GEMDetId(int region, int ring, int station, int layer, int chamber, int ieta)
0218   const GEMDetId gem_chamber_id{csc_id.zendcap(), 1, csc_id.station(), 0, csc_id.chamber(), 0};
0219
0220   // detector
0221   b_region_ = csc_id.zendcap();
0222   b_station_ = csc_id.station();
0223   b_ring_ = csc_id.ring();
0224   b_chamber_ = csc_id.chamber();
0225   b_gem_chamber_has_error_ =
0226       use_gem_daq_status_ ? checkGEMChamberStatus(gem_chamber_id, gem_oh_status_collection, gem_vfat_status_collection)
0227                           : false;
0228   // GEMCSCSegment
0229   b_norm_chi2_ = static_cast<float>(gem_csc_segment.chi2() / gem_csc_segment.nRecHits());
0230   // CSCSegment
0231   b_csc_norm_chi2_ = static_cast<float>(csc_segment.chi2() / csc_segment.nRecHits());
0232   b_csc_num_hit_ = csc_segment.nRecHits();
0233   b_csc_is_muon_ = isMuonSegment(csc_segment, muon_segment_vec);
0234   // GEMRecHit
0235   b_gem_num_hit = gem_hit_vec.size();
0236
0237   b_gem_layer_.reserve(b_gem_num_hit);
0238   b_gem_ieta_.reserve(b_gem_num_hit);
0239   b_gem_first_strip_.reserve(b_gem_num_hit);
0240   b_gem_cls_.reserve(b_gem_num_hit);
0241
0242   for (const GEMRecHit& gem_hit : gem_hit_vec) {
0243     const GEMDetId gem_id = gem_hit.gemId();
0244
0245     b_gem_layer_.push_back(gem_id.layer());
0246     b_gem_ieta_.push_back(gem_id.ieta());
0247     b_gem_first_strip_.push_back(gem_hit.firstClusterStrip());
0248     b_gem_cls_.push_back(gem_hit.clusterSize());
0249   }  // GEMRecHit
0250
0251   return true;
0252 }
0253
0254 // Returns a vector of CSCSegments, which are matched with muons.
0255 std::vector<const CSCSegment*> GEMCSCCoincidenceRateAnalyzer::findMuonSegments(
0256     const reco::MuonCollection* muon_collection) {
0257   std::vector<const CSCSegment*> muon_segment_vec;
0258
0259   for (const reco::Muon& muon : *muon_collection) {
0260     for (const reco::MuonChamberMatch& chamber_match : muon.matches()) {
0261       if (chamber_match.detector() != MuonSubdetId::CSC) {
0262         continue;
0263       }
0264
0265       if (checkCSCChamberType(chamber_match.id)) {
0266         for (const reco::MuonSegmentMatch& segment_match : chamber_match.segmentMatches) {
0267           if (segment_match.isMask(reco::MuonSegmentMatch::BestInStationByDR)) {
0268             muon_segment_vec.push_back(segment_match.cscSegmentRef.get());
0269             break;
0270           }
0271         }  // MuonSegmentMatch
0272       }  // checkCSCChamberType
0273     }  // MuonChamberMatch
0274   }  // MuonCollection
0275
0276   return muon_segment_vec;
0277 }
0278
0279 // Returns a boolean indicating whether the given CSCSegment is matched with a
0280 // muon.
0281 // TODO better segment comparison
0282 bool GEMCSCCoincidenceRateAnalyzer::isMuonSegment(const CSCSegment& csc_segment,
0283                                                   const std::vector<const CSCSegment*> muon_segment_vec) {
0284   bool found = false;
0285
0286   const CSCDetId csc_id = csc_segment.cscDetId();
0287   for (const CSCSegment* muon_segment : muon_segment_vec) {
0288     if (csc_id != muon_segment->cscDetId())
0289       continue;
0290     if (csc_segment.localPosition().x() != muon_segment->localPosition().x())
0291       continue;
0292     if (csc_segment.localPosition().y() != muon_segment->localPosition().y())
0293       continue;
0294     if (csc_segment.localPosition().z() != muon_segment->localPosition().z())
0295       continue;
0296     if (csc_segment.time() != muon_segment->time())
0297       continue;
0298
0299     found = true;
0300   }
0301
0302   return found;
0303 }
0304
0305 // Returns a boolean indicating whether the GEMChamber has any DAQ error.
0306 bool GEMCSCCoincidenceRateAnalyzer::checkGEMChamberStatus(const GEMDetId& chamber_id,
0307                                                           const GEMOHStatusCollection* oh_status_collection,
0308                                                           const GEMVFATStatusCollection* vfat_status_collection) {
0309   const bool has_error = true;
0310   if (not use_gem_daq_status_) {
0311     edm::LogError(log_category_)
0312         << "useGEMDAQStatus is false but checkGEMChamberStatus is called. gem_chamber_has_error will be set to false";
0313     return not has_error;
0314   }
0315
0316   for (auto iter = oh_status_collection->begin(); iter != oh_status_collection->end(); iter++) {
0317     const auto [oh_id, range] = (*iter);
0318     if (chamber_id != oh_id) {
0319       continue;
0320     }
0321
0322     for (auto oh_status = range.first; oh_status != range.second; oh_status++) {
0323       if (oh_status->isBad()) {
0324         // GEMOHStatus is bad. Mask this chamber.
0325         return has_error;
0326       }  // isBad
0327     }  // range
0328   }  // collection
0329
0330   for (auto iter = vfat_status_collection->begin(); iter != vfat_status_collection->end(); iter++) {
0331     const auto [vfat_id, range] = (*iter);
0332     if (chamber_id != vfat_id.chamberId()) {
0333       continue;
0334     }
0335     for (auto vfat_status = range.first; vfat_status != range.second; vfat_status++) {
0336       if (vfat_status->isBad()) {
0337         return has_error;
0338       }
0339     }  // range
0340   }  // collection
0341
0342   return not has_error;
0343 }
0344
0345 // Returns a boolean indicating whether to allow a CSCDetId.
0346 bool GEMCSCCoincidenceRateAnalyzer::checkCSCChamberType(const CSCDetId& csc_id) {
0347   const auto csc_chamber_type = static_cast<unsigned int>(csc_id.iChamberType());
0348   return std::find(csc_whitelist_.begin(), csc_whitelist_.end(), csc_chamber_type) != csc_whitelist_.end();
0349 }
0350
0351 DEFINE_FWK_MODULE(GEMCSCCoincidenceRateAnalyzer);
0352 #endif  // RecoLocalMuon_GEMCSCSegment_GEMCSCCoincidenceRateAnalyzer_h