Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:07:26 

0001 #include "DQM/GEM/plugins/GEMEfficiencyAnalyzer.h"
0002 
0003 #include "FWCore/Framework/interface/ConsumesCollector.h"
0004 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0005 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0006 #include "DataFormats/GeometryCommonDetAlgo/interface/ErrorFrameTransformer.h"
0007 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0008 #include "DataFormats/Math/interface/deltaPhi.h"
0009 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0010 #include "Validation/MuonHits/interface/MuonHitHelper.h"
0011 #include "Validation/MuonGEMHits/interface/GEMValidationUtils.h"
0012 
0013 GEMEfficiencyAnalyzer::GEMEfficiencyAnalyzer(const edm::ParameterSet& ps)
0014     : GEMDQMEfficiencySourceBase(ps),
0015       kGEMGeometryTokenBeginRun_(esConsumes<edm::Transition::BeginRun>()),
0016       kTransientTrackBuilderToken_(
0017           esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))),
0018       kGEMRecHitCollectionToken_(consumes<GEMRecHitCollection>(ps.getUntrackedParameter<edm::InputTag>("recHitTag"))),
0019       kMuonViewToken_(consumes<edm::View<reco::Muon> >(ps.getUntrackedParameter<edm::InputTag>("muonTag"))),
0020       kMuonTrackTypeName_(ps.getUntrackedParameter<std::string>("muonTrackType")),
0021       kMuonTrackType_(getMuonTrackType(kMuonTrackTypeName_)),
0022       kMuonName_(TString(ps.getUntrackedParameter<std::string>("muonName"))),
0023       kFolder_(ps.getUntrackedParameter<std::string>("folder")),
0024       kScenario_(getScenarioOption(ps.getUntrackedParameter<std::string>("scenario"))),
0025       kStartingStateType_(getStartingStateType(ps.getUntrackedParameter<std::string>("startingStateType"))),
0026       kMuonSubdetForGEM_({
0027           ps.getUntrackedParameter<std::vector<int> >("muonSubdetForGE0"),
0028           ps.getUntrackedParameter<std::vector<int> >("muonSubdetForGE11"),
0029           ps.getUntrackedParameter<std::vector<int> >("muonSubdetForGE21"),
0030       }),
0031       kCSCForGEM_({
0032           ps.getUntrackedParameter<std::vector<int> >("cscForGE0"),
0033           ps.getUntrackedParameter<std::vector<int> >("cscForGE11"),
0034           ps.getUntrackedParameter<std::vector<int> >("cscForGE21"),
0035       }),
0036       kMuonSegmentMatchDRCut_(static_cast<float>(ps.getUntrackedParameter<double>("muonSegmentMatchDRCut"))),
0037       kMuonPtMinCuts_({
0038           ps.getUntrackedParameter<double>("muonPtMinCutGE0"),
0039           ps.getUntrackedParameter<double>("muonPtMinCutGE11"),
0040           ps.getUntrackedParameter<double>("muonPtMinCutGE21"),
0041       }),
0042       kMuonEtaMinCuts_({
0043           ps.getUntrackedParameter<double>("muonEtaMinCutGE0"),
0044           ps.getUntrackedParameter<double>("muonEtaMinCutGE11"),
0045           ps.getUntrackedParameter<double>("muonEtaMinCutGE21"),
0046       }),
0047       kMuonEtaMaxCuts_({
0048           ps.getUntrackedParameter<double>("muonEtaMaxCutGE0"),
0049           ps.getUntrackedParameter<double>("muonEtaMaxCutGE11"),
0050           ps.getUntrackedParameter<double>("muonEtaMaxCutGE21"),
0051       }),
0052       kPropagationErrorRCut_(static_cast<float>(ps.getUntrackedParameter<double>("propagationErrorRCut"))),
0053       kPropagationErrorPhiCut_(static_cast<float>(ps.getUntrackedParameter<double>("propagationErrorPhiCut"))),
0054       kBoundsErrorScale_(static_cast<float>(ps.getUntrackedParameter<double>("boundsErrorScale"))),
0055       kMatchingMetric_(getMatchingMetric(ps.getUntrackedParameter<std::string>("matchingMetric"))),
0056       kMatchingCut_(static_cast<float>(ps.getUntrackedParameter<double>("matchingCut"))),
0057       kMuonPtBins_(ps.getUntrackedParameter<std::vector<double> >("muonPtBins")),
0058       kMuonEtaNbins_({
0059           ps.getUntrackedParameter<int>("muonEtaNbinsGE0"),
0060           ps.getUntrackedParameter<int>("muonEtaNbinsGE11"),
0061           ps.getUntrackedParameter<int>("muonEtaNbinsGE21"),
0062       }),
0063       kMuonEtaLow_({
0064           ps.getUntrackedParameter<double>("muonEtaLowGE0"),
0065           ps.getUntrackedParameter<double>("muonEtaLowGE11"),
0066           ps.getUntrackedParameter<double>("muonEtaLowGE21"),
0067       }),
0068       kMuonEtaUp_({
0069           ps.getUntrackedParameter<double>("muonEtaUpGE0"),
0070           ps.getUntrackedParameter<double>("muonEtaUpGE11"),
0071           ps.getUntrackedParameter<double>("muonEtaUpGE21"),
0072       }),
0073       kModeDev_(ps.getUntrackedParameter<bool>("modeDev")) {
0074   muon_service_ =
0075       std::make_unique<MuonServiceProxy>(ps.getParameter<edm::ParameterSet>("ServiceParameters"), consumesCollector());
0076 }
0077 
0078 GEMEfficiencyAnalyzer::~GEMEfficiencyAnalyzer() {}
0079 
0080 void GEMEfficiencyAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0081   edm::ParameterSetDescription desc;
0082   // GEMDQMEfficiencySourceBase
0083   desc.addUntracked<edm::InputTag>("ohStatusTag", edm::InputTag("muonGEMDigis", "OHStatus"));
0084   desc.addUntracked<edm::InputTag>("vfatStatusTag", edm::InputTag("muonGEMDigis", "VFATStatus"));
0085   desc.addUntracked<bool>("monitorGE11", true);
0086   desc.addUntracked<bool>("monitorGE21", false);
0087   desc.addUntracked<bool>("monitorGE0", false);
0088   desc.addUntracked<bool>("maskChamberWithError", false);
0089   desc.addUntracked<std::string>("logCategory", "GEMEfficiencyAnalyzer");
0090 
0091   // GEMEfficiencyAnalyzer
0092   desc.addUntracked<edm::InputTag>("recHitTag", edm::InputTag("gemRecHits"));
0093   desc.addUntracked<edm::InputTag>("muonTag", edm::InputTag("muons"));
0094   desc.addUntracked<bool>("modeDev", false);
0095   desc.addUntracked<std::string>("muonTrackType", "OuterTrack");
0096   desc.addUntracked<std::string>("muonName", "STA Muon");
0097   desc.addUntracked<std::string>("folder", "GEM/Efficiency/muonSTA");
0098   desc.addUntracked<std::string>("scenario", "pp");
0099   //
0100   desc.addUntracked<std::string>("startingStateType", "OutermostMeasurementState");
0101   desc.addUntracked<double>("muonSegmentMatchDRCut", 5.0f);  // for cosmics, in cm,  TODO tune
0102   // muon pt cut
0103   desc.addUntracked<double>("muonPtMinCutGE0", 20.0f);
0104   desc.addUntracked<double>("muonPtMinCutGE11", 20.0f);
0105   desc.addUntracked<double>("muonPtMinCutGE21", 20.0f);
0106   // muon abs eta cut for GE11
0107   desc.addUntracked<double>("muonEtaMinCutGE11", 1.5);
0108   desc.addUntracked<double>("muonEtaMaxCutGE11", 2.2);
0109   // muon abs eta cut for GE21
0110   desc.addUntracked<double>("muonEtaMinCutGE21", 1.5);
0111   desc.addUntracked<double>("muonEtaMaxCutGE21", 2.5);
0112   // muon abs eta cut for GE0
0113   desc.addUntracked<double>("muonEtaMinCutGE0", 2.0);
0114   desc.addUntracked<double>("muonEtaMaxCutGE0", 3.0);
0115   // propagation error cuts
0116   desc.addUntracked<double>("propagationErrorRCut", 0.5);    // cm
0117   desc.addUntracked<double>("propagationErrorPhiCut", 0.2);  // degree
0118   //
0119   desc.addUntracked<double>("boundsErrorScale", -2.0);  // TODO tune
0120   // matching
0121   desc.addUntracked<std::string>("matchingMetric", "DeltaPhi");
0122   desc.addUntracked<double>("matchingCut", 0.2);  // DeltaPhi for pp, in degree TODO tune
0123   // for MinotorElement
0124   const std::vector<double> default_pt_bins{
0125       0, 5, 10., 20., 30., 40., 50., 60., 70., 80., 90., 100., 110.};  // actually edges
0126   desc.addUntracked<std::vector<double> >("muonPtBins", default_pt_bins);
0127   // GE11
0128   desc.addUntracked<int>("muonEtaNbinsGE11", 9);  // bin width = 0.1
0129   desc.addUntracked<double>("muonEtaLowGE11", 1.4);
0130   desc.addUntracked<double>("muonEtaUpGE11", 2.3);
0131   // GE21
0132   desc.addUntracked<int>("muonEtaNbinsGE21", 12);  // bin width = 0.1
0133   desc.addUntracked<double>("muonEtaLowGE21", 1.4);
0134   desc.addUntracked<double>("muonEtaUpGE21", 2.6);
0135   // GE0
0136   desc.addUntracked<int>("muonEtaNbinsGE0", 12);  // bin width = 0.1
0137   desc.addUntracked<double>("muonEtaLowGE0", 1.9);
0138   desc.addUntracked<double>("muonEtaUpGE0", 3.1);
0139 
0140   // MuonSubdetId's are listed in DataFormats/MuonDetId/interface/MuonSubdetId.h
0141   desc.addUntracked<std::vector<int> >("muonSubdetForGE0", {});  // allow all muon subdetectors. TODO optimzie.
0142   desc.addUntracked<std::vector<int> >("muonSubdetForGE11", {});
0143   desc.addUntracked<std::vector<int> >("muonSubdetForGE21", {});
0144   // INFO when muonTrackType is "CombinedTrack" or "OuterTrack"
0145   // https://github.com/cms-sw/cmssw/blob/CMSSW_12_4_0_pre3/DataFormats/MuonDetId/interface/CSCDetId.h#L187-L193
0146   // assumed to be the same area.
0147   desc.addUntracked<std::vector<int> >("cscForGE11", {1, 2});  // ME1a, ME1b
0148   desc.addUntracked<std::vector<int> >("cscForGE21", {});      // all CSCSegments allowed
0149   desc.addUntracked<std::vector<int> >("cscForGE0", {});       // all CSCSegments allowed
0150 
0151   // ServiceParameters for MuonServiceProxy
0152   // This will be initialized in the cfi file
0153   edm::ParameterSetDescription service_parameters;
0154   service_parameters.setAllowAnything();
0155   desc.add<edm::ParameterSetDescription>("ServiceParameters", service_parameters);
0156 
0157   descriptions.add("gemEfficiencyAnalyzerDefault", desc);
0158 }
0159 
0160 // convert a string to enum
0161 GEMEfficiencyAnalyzer::MatchingMetric GEMEfficiencyAnalyzer::getMatchingMetric(const std::string name) {
0162   MatchingMetric method;
0163 
0164   if (name == "DeltaPhi") {
0165     method = MatchingMetric::kDeltaPhi;
0166 
0167   } else if (name == "RdPhi") {
0168     method = MatchingMetric::kRdPhi;
0169 
0170   } else {
0171     edm::LogError(kLogCategory_) << "received an unexpected MatchingMetric: " << name
0172                                  << " -> MatchingMetric::kDeltaPhi will be used instead.";
0173     method = MatchingMetric::kDeltaPhi;
0174   }
0175 
0176   return method;
0177 }
0178 
0179 // convert a string to enum
0180 GEMEfficiencyAnalyzer::StartingStateType GEMEfficiencyAnalyzer::getStartingStateType(const std::string name) {
0181   StartingStateType type;
0182 
0183   if (name == "InnermostMeasurementState") {
0184     type = StartingStateType::kInnermostMeasurementState;
0185 
0186   } else if (name == "OutermostMeasurementState") {
0187     type = StartingStateType::kOutermostMeasurementState;
0188 
0189   } else if (name == "StateOnSurfaceWithCSCSegment") {
0190     type = StartingStateType::kStateOnSurfaceWithCSCSegment;
0191 
0192   } else if (name == "AlignmentStyle") {
0193     type = StartingStateType::kAlignmentStyle;
0194 
0195   } else {
0196     edm::LogError(kLogCategory_) << "received an unexpected StartingStateType: " << name
0197                                  << " -> StartingStateType::kOutermostMeasurementState will be used instead.";
0198     type = StartingStateType::kOutermostMeasurementState;
0199   }
0200 
0201   return type;
0202 }
0203 
0204 // convert a string to enum
0205 reco::Muon::MuonTrackType GEMEfficiencyAnalyzer::getMuonTrackType(const std::string name) {
0206   reco::Muon::MuonTrackType muon_track_type;
0207 
0208   // DO NOT ALLOW TYPO
0209   if (name == "InnerTrack") {
0210     muon_track_type = reco::Muon::MuonTrackType::InnerTrack;
0211 
0212   } else if (name == "OuterTrack") {
0213     muon_track_type = reco::Muon::MuonTrackType::OuterTrack;
0214 
0215   } else if (name == "CombinedTrack") {
0216     muon_track_type = reco::Muon::MuonTrackType::CombinedTrack;
0217 
0218   } else {
0219     edm::LogError(kLogCategory_) << "received an unexpected reco::Muon::MuonTrackType: " << name
0220                                  << " --> OuterTrack will be used instead.";
0221 
0222     muon_track_type = reco::Muon::MuonTrackType::OuterTrack;
0223   }
0224 
0225   return muon_track_type;
0226 }
0227 
0228 GEMEfficiencyAnalyzer::ScenarioOption GEMEfficiencyAnalyzer::getScenarioOption(const std::string name) {
0229   ScenarioOption scenario;
0230   if (name == "pp") {
0231     scenario = ScenarioOption::kPP;
0232 
0233   } else if (name == "cosmics") {
0234     scenario = ScenarioOption::kCosmics;
0235 
0236   } else if (name == "HeavyIons") {
0237     scenario = ScenarioOption::kHeavyIons;
0238 
0239     edm::LogInfo(kLogCategory_) << "The scenario is set to \"HeavyIons\""
0240                                 << " but there is no strategy dedicated to"
0241                                 << "\"HeavyIons\" scenario. The strategy for "
0242                                 << "the \"pp\" scenario will be used insteqad.";
0243 
0244   } else {
0245     scenario = ScenarioOption::kPP;
0246 
0247     edm::LogError(kLogCategory_) << "received an unexpected ScenarioOption: " << name
0248                                  << ". Choose from (\"pp\", \"cosmics\", \"HeavyIons\")"
0249                                  << " --> pp will be used instead.";
0250   }
0251 
0252   return scenario;
0253 }
0254 
0255 void GEMEfficiencyAnalyzer::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const& setup) {
0256   ibooker.setCurrentFolder(kFolder_);
0257 
0258   const GEMGeometry* gem = nullptr;
0259   if (auto handle = setup.getHandle(kGEMGeometryTokenBeginRun_)) {
0260     gem = handle.product();
0261   } else {
0262     edm::LogError(kLogCategory_ + "|bookHistograms") << "failed to get GEMGeometry";
0263     return;
0264   }
0265 
0266   for (const GEMStation* station : gem->stations()) {
0267     const int region_id = station->region();
0268     const int station_id = station->station();
0269 
0270     if (skipGEMStation(station_id)) {
0271       continue;
0272     }
0273 
0274     ////////////////////////////////////////////////////////////////////////////
0275     // Region-Station
0276     ////////////////////////////////////////////////////////////////////////////
0277     {  // shadowing to reuse short variable names
0278       const GEMDetId key = getReStKey(region_id, station_id);
0279       const TString suffix = GEMUtils::getSuffixName(region_id, station_id);
0280       const TString title = kMuonName_ + GEMUtils::getSuffixTitle(region_id, station_id);
0281 
0282       // sources for eff. vs muon pt
0283       TH1F* h_muon_pt = new TH1F("muon_pt" + suffix, title, kMuonPtBins_.size() - 1, &kMuonPtBins_[0]);
0284       me_muon_pt_[key] = ibooker.book1D(h_muon_pt->GetName(), h_muon_pt);
0285       me_muon_pt_[key]->setAxisTitle("Muon p_{T} [GeV]", 1);
0286       me_muon_pt_matched_[key] = bookNumerator1D(ibooker, me_muon_pt_[key]);
0287 
0288       // sources for eff. vs muon eta
0289       me_muon_eta_[key] = ibooker.book1D("muon_eta" + suffix,
0290                                          title,
0291                                          kMuonEtaNbins_.at(station_id),
0292                                          kMuonEtaLow_.at(station_id),
0293                                          kMuonEtaUp_.at(station_id));
0294       me_muon_eta_[key]->setAxisTitle("Muon |#eta|", 1);
0295       me_muon_eta_matched_[key] = bookNumerator1D(ibooker, me_muon_eta_[key]);
0296 
0297       // sources for eff. vs muon phi
0298       me_muon_phi_[key] = ibooker.book1D("muon_phi" + suffix, title, 36, -180, 180);
0299       me_muon_phi_[key]->setAxisTitle("Muon #phi [deg]");
0300       me_muon_phi_matched_[key] = bookNumerator1D(ibooker, me_muon_phi_[key]);
0301 
0302       if (kModeDev_) {
0303         // without cuts except the fiducial cut
0304         TH1F* h_muon_pt_all = new TH1F("muon_pt_all" + suffix, title, kMuonPtBins_.size() - 1, &kMuonPtBins_[0]);
0305         me_muon_pt_all_[key] = ibooker.book1D(h_muon_pt_all->GetName(), h_muon_pt_all);
0306         me_muon_pt_all_[key]->setAxisTitle("Muon p_{T} [GeV]", 1);
0307         me_muon_pt_all_matched_[key] = bookNumerator1D(ibooker, me_muon_pt_all_[key]);
0308 
0309         me_muon_eta_all_[key] = ibooker.book1D("muon_eta_all" + suffix,
0310                                                title,
0311                                                kMuonEtaNbins_.at(station_id),
0312                                                kMuonEtaLow_.at(station_id),
0313                                                kMuonEtaUp_.at(station_id));
0314         me_muon_eta_all_[key]->setAxisTitle("Muon |#eta|", 1);
0315         me_muon_eta_all_matched_[key] = bookNumerator1D(ibooker, me_muon_eta_all_[key]);
0316 
0317         me_muon_charge_[key] = ibooker.book1D("muon_charge" + suffix, title, 3, -1.5, 1.5);
0318         me_muon_charge_[key]->setAxisTitle("Muon charge", 1);
0319         me_muon_charge_matched_[key] = bookNumerator1D(ibooker, me_muon_charge_[key]);
0320       }
0321     }  // shadowing
0322 
0323     ////////////////////////////////////////////////////////////////////////////
0324     // Region - Station - Layer
0325     ////////////////////////////////////////////////////////////////////////////
0326     const std::vector<const GEMSuperChamber*> superchamber_vec = station->superChambers();
0327     if (not checkRefs(superchamber_vec)) {
0328       edm::LogError(kLogCategory_) << "got an invalid ptr from GEMStation::superChambers";
0329       return;
0330     }
0331 
0332     const std::vector<const GEMChamber*> chamber_vec = superchamber_vec.front()->chambers();
0333     if (not checkRefs(chamber_vec)) {
0334       edm::LogError(kLogCategory_) << "got an invalid ptr from GEMSuperChamber::chambers";
0335       return;
0336     }
0337 
0338     // we actually loop over layers
0339     for (const GEMChamber* chamber : chamber_vec) {
0340       const int layer_id = chamber->id().layer();
0341 
0342       {  // shadowing
0343         const GEMDetId key = getReStLaKey(chamber->id());
0344         const TString suffix = GEMUtils::getSuffixName(region_id, station_id, layer_id);
0345         const TString title = kMuonName_ + GEMUtils::getSuffixTitle(region_id, station_id, layer_id);
0346 
0347         me_chamber_ieta_[key] = bookChamberEtaPartition(ibooker, "chamber_ieta" + suffix, title, station);
0348         me_chamber_ieta_matched_[key] = bookNumerator2D(ibooker, me_chamber_ieta_[key]);
0349 
0350         if (kModeDev_) {
0351           me_prop_path_length_[key] = ibooker.book1D("prop_path_length" + suffix, title, 50, 0.0, 5.0);
0352           me_prop_path_length_[key]->setAxisTitle("Propagation path length [cm]", 1);
0353           me_prop_path_length_matched_[key] = bookNumerator1D(ibooker, me_prop_path_length_[key]);
0354 
0355           // prop. r error in the global coordinates
0356           me_prop_err_r_[key] = ibooker.book1D("prop_err_r" + suffix, title, 60, 0.0, 3.0);
0357           me_prop_err_r_[key]->setAxisTitle("Propagation global #sigma_{R} [cm]", 1);
0358           me_prop_err_r_matched_[key] = bookNumerator1D(ibooker, me_prop_err_r_[key]);
0359 
0360           // prop. r error in the global coordinates
0361           me_prop_err_phi_[key] = ibooker.book1D("prop_err_phi" + suffix, title, 50, 0.0, 1.0);
0362           me_prop_err_phi_[key]->setAxisTitle("Propagation's global #sigma_{#phi} [deg]", 1);
0363           me_prop_err_phi_matched_[key] = bookNumerator1D(ibooker, me_prop_err_phi_[key]);
0364 
0365           // cutflow
0366           me_cutflow_[key] = ibooker.book1D("cutflow" + suffix, title, 5, 0.5, 5.5);
0367           me_cutflow_[key]->setBinLabel(1, "All");
0368           me_cutflow_[key]->setBinLabel(2, Form("#sigma_{R} < %.3f cm", kPropagationErrorRCut_));
0369           me_cutflow_[key]->setBinLabel(3, Form("#sigma_{phi} < %.3f deg", kPropagationErrorPhiCut_));
0370           me_cutflow_[key]->setBinLabel(4, Form("p_{T} > %.1f GeV", kMuonPtMinCuts_.at(station_id)));
0371           me_cutflow_[key]->setBinLabel(
0372               5, Form("%.2f < |#eta| < %.2f", kMuonEtaMinCuts_.at(station_id), kMuonEtaMaxCuts_.at(station_id)));
0373 
0374           me_cutflow_matched_[key] = bookNumerator1D(ibooker, me_cutflow_.at(key));
0375         }
0376       }  // shadowing
0377     }    // GEMChamber
0378 
0379     ////////////////////////////////////////////////////////////////////////////
0380     // Region - Station - iEta
0381     ////////////////////////////////////////////////////////////////////////////
0382     const std::vector<const GEMEtaPartition*> eta_partition_vec = chamber_vec.front()->etaPartitions();
0383     if (not checkRefs(eta_partition_vec)) {
0384       edm::LogError(kLogCategory_) << "got an invalid ptr from GEMChamber::etaPartitions";
0385       continue;
0386     }
0387 
0388     for (const GEMEtaPartition* eta_partition : eta_partition_vec) {
0389       const int ieta = eta_partition->id().ieta();
0390 
0391       {  // shadowing
0392         const GEMDetId key = getReStEtKey(eta_partition->id());
0393         const TString gem_label = TString::Format("GE%d1-%c-E%d", station_id, (region_id > 0 ? 'P' : 'M'), ieta);
0394         const TString suffix = "_" + gem_label;
0395         const TString title = kMuonName_ + " " + gem_label;
0396 
0397         // FIXME
0398         const float dphi_up = (kMatchingMetric_ == MatchingMetric::kDeltaPhi) ? kMatchingCut_
0399                               : (kScenario_ == ScenarioOption::kCosmics)      ? 1.0
0400                                                                               : 0.2;
0401         me_residual_phi_[key] = ibooker.book1D("residual_phi" + suffix, title, 41, -dphi_up, dphi_up);
0402         me_residual_phi_[key]->setAxisTitle("Residual in global #phi [deg]", 1);
0403 
0404         if (kModeDev_) {
0405           // matching metric
0406           std::string matching_metric_x_title;
0407           if (kMatchingMetric_ == MatchingMetric::kDeltaPhi) {
0408             matching_metric_x_title = "#Delta#phi [deg]";
0409 
0410           } else if (kMatchingMetric_ == MatchingMetric::kRdPhi) {
0411             matching_metric_x_title = "R#Delta#phi [cm]";
0412 
0413           } else {
0414             matching_metric_x_title = "UNKNOWN METRIC";
0415           }
0416 
0417           // matching metrics without any cuts
0418           me_matching_metric_all_[key] =
0419               ibooker.book1D("matching_metric_all" + suffix, title, 101, -3 * kMatchingCut_, 3 * kMatchingCut_);
0420           me_matching_metric_all_[key]->setAxisTitle(matching_metric_x_title, 1);
0421 
0422           // matching metrics after cuts
0423           me_matching_metric_[key] =
0424               ibooker.book1D("matching_metric" + suffix, title, 101, -kMatchingCut_, kMatchingCut_);
0425           me_matching_metric_[key]->setAxisTitle(matching_metric_x_title, 1);
0426 
0427           // residuals in the global phi for muons (q < 0)
0428           me_residual_phi_muon_[key] =
0429               ibooker.book1D("residual_phi_muon" + suffix, title + " (#mu, q < 0)", 50, -0.5, 0.5);
0430           me_residual_phi_muon_[key]->setAxisTitle("Residual in global #phi [deg]", 1);
0431           me_residual_phi_muon_[key]->setAxisTitle("Number of muons", 2);
0432 
0433           // residuals in the global phi for anti-muons (q > 0)
0434           me_residual_phi_antimuon_[key] =
0435               ibooker.book1D("residual_phi_antimuon" + suffix, title + " (#tilde{#mu}, q > 0)", 50, -0.5, 0.5);
0436           me_residual_phi_antimuon_[key]->setAxisTitle("Residual in global #phi [deg]", 1);
0437           me_residual_phi_antimuon_[key]->setAxisTitle("Number of anti-muons", 2);
0438 
0439           // residuals in the local x
0440           me_residual_x_[key] = ibooker.book1D("residual_x" + suffix, title, 60, -1.5, 1.5);
0441           me_residual_x_[key]->setAxisTitle("Residual in local X [cm]", 1);
0442 
0443           // residuals in the local y
0444           me_residual_y_[key] = ibooker.book1D("residual_y" + suffix, title, 48, -12.0, 12.0);
0445           me_residual_y_[key]->setAxisTitle("Residual in local Y [cm]", 1);
0446 
0447           // the strip difference
0448           me_residual_strip_[key] = ibooker.book1D("residual_strip" + suffix, title, 21, -10.0, 10.0);
0449           me_residual_strip_[key]->setAxisTitle("propagation strip - hit strip", 1);
0450         }
0451       }  // shadowing
0452     }    // GEMEtaPartition
0453   }      // GEMStataion
0454 }
0455 
0456 // In the `cosmics` scenario, TODO doc
0457 bool GEMEfficiencyAnalyzer::isInsideOut(const reco::Track& track) {
0458   return track.innerPosition().mag2() > track.outerPosition().mag2();
0459 }
0460 
0461 //
0462 void GEMEfficiencyAnalyzer::buildGEMLayers(const GEMGeometry* gem) {
0463   std::map<GEMDetId, std::vector<const GEMChamber*> > chambers_per_layer;
0464 
0465   for (const GEMStation* station : gem->stations()) {
0466     const int region_id = station->region();
0467     const int station_id = station->station();
0468     const bool is_ge11 = station_id == 1;
0469 
0470     if (skipGEMStation(station_id)) {
0471       continue;
0472     }
0473 
0474     for (const GEMSuperChamber* superchamber : station->superChambers()) {
0475       // GE11: chamber == 0 for even chambers, chamber == 1 for odd chambers
0476       // GE21 and GE0: chamber == 0 for all chambers
0477       const int chamber_id = is_ge11 ? superchamber->id().chamber() % 2 : 0;
0478 
0479       for (const GEMChamber* chamber : superchamber->chambers()) {
0480         const int layer_id = chamber->id().layer();
0481 
0482         const GEMDetId key{region_id, 1, station_id, layer_id, chamber_id, 0};
0483 
0484         if (chambers_per_layer.find(key) == chambers_per_layer.end()) {
0485           chambers_per_layer.insert({key, std::vector<const GEMChamber*>()});
0486         }
0487         chambers_per_layer.at(key).push_back(chamber);
0488       }  // GEMChamber => iterate over layer ids
0489     }    // GEMSuperChamber => iterate over chamber ids
0490   }      // GEMStation
0491 
0492   gem_layers_.reserve(chambers_per_layer.size());
0493   for (auto [gem_id, chambers] : chambers_per_layer) {
0494     // layer position and rotation
0495     const float z_origin = chambers.front()->position().z();
0496     Surface::PositionType position{0.f, 0.f, z_origin};
0497     Surface::RotationType rotation;
0498 
0499     // eta partitions should have same R and Z spans.
0500     // XXX is it true?
0501     auto [r_min, r_max] = chambers.front()->surface().rSpan();
0502     auto [z_min, z_max] = chambers.front()->surface().zSpan();
0503 
0504     z_min -= z_origin;
0505     z_max -= z_origin;
0506 
0507     // the bounds from min and max R and Z in the local coordinates.
0508     SimpleDiskBounds* bounds = new SimpleDiskBounds(r_min, r_max, z_min, z_max);
0509     const Disk::DiskPointer layer = Disk::build(position, rotation, bounds);
0510 
0511     gem_layers_.emplace_back(layer, chambers, gem_id);
0512 
0513     LogDebug(kLogCategory_) << gem_id
0514                             << Form(" ==> (z_origin, z_min, z_max) = (%.2f, %.2f, %.2f)", z_origin, z_min, z_max);
0515   }  // ring
0516 }
0517 
0518 // TODO doc
0519 // See https://twiki.cern.ch/twiki/pub/CMS/GEMPPDOfflineDQM/check-muon-direction.pdf
0520 bool GEMEfficiencyAnalyzer::checkPropagationDirection(const reco::Track* track, const GEMLayer& layer) {
0521   const bool is_same_region = track->eta() * layer.id.region() > 0;
0522 
0523   bool skip = false;
0524   if (kScenario_ == ScenarioOption::kCosmics) {
0525     float p2_in = track->innerMomentum().mag2();
0526     float p2_out = track->outerMomentum().mag2();
0527     if (isInsideOut(*track))
0528       std::swap(p2_in, p2_out);
0529     const bool is_outgoing = p2_in > p2_out;
0530 
0531     skip = (is_outgoing xor is_same_region);
0532 
0533   } else {
0534     // beam scenario
0535     skip = not is_same_region;
0536   }
0537 
0538   return skip;
0539 }
0540 
0541 GEMEfficiencyAnalyzer::StartingState GEMEfficiencyAnalyzer::buildStartingState(
0542     const reco::Muon& muon, const reco::TransientTrack& transient_track, const GEMLayer& gem_layer) {
0543   bool found = false;
0544   TrajectoryStateOnSurface state;
0545   DetId det_id;
0546 
0547   switch (kStartingStateType_) {
0548     case StartingStateType::kOutermostMeasurementState: {
0549       std::tie(found, state, det_id) = getOutermostMeasurementState(transient_track);
0550       break;
0551     }
0552     case StartingStateType::kInnermostMeasurementState: {
0553       std::tie(found, state, det_id) = getInnermostMeasurementState(transient_track);
0554       break;
0555     }
0556     case StartingStateType::kStateOnSurfaceWithCSCSegment: {
0557       std::tie(found, state, det_id) = buildStateOnSurfaceWithCSCSegment(muon, transient_track, gem_layer);
0558       break;
0559     }
0560     case StartingStateType::kAlignmentStyle: {
0561       std::tie(found, state, det_id) = buildStartingStateAlignmentStyle(muon, transient_track, gem_layer);
0562       break;
0563     }
0564     default: {
0565       edm::LogError(kLogCategory_) << "got an unexpected StartingStateType";
0566       break;
0567     }
0568   }
0569 
0570   found &= state.isValid();
0571 
0572   if (found and (det_id.det() == DetId::Detector::Muon)) {
0573     found &= isMuonSubdetAllowed(det_id, gem_layer.id.station());
0574   }
0575 
0576   if (found) {
0577     if (MuonHitHelper::isGEM(det_id)) {
0578       const GEMDetId start_id{det_id};
0579 
0580       const bool are_same_region = gem_layer.id.region() == start_id.region();
0581       const bool are_same_station = gem_layer.id.station() == start_id.station();
0582       const bool are_same_layer = gem_layer.id.layer() == start_id.layer();
0583       if (are_same_region and are_same_station and are_same_layer) {
0584         LogDebug(kLogCategory_)
0585             << "The starting detector of the muon propagation is same with the destination. Skip this propagation.";
0586         found = false;
0587       }
0588     }  // isGEM
0589   }    // found
0590 
0591   return std::make_tuple(found, state, det_id);
0592 }
0593 
0594 // Use the innermost measurement state as an initial state for the muon propagation.
0595 // NOTE If the analyzer uses global or standalone muons and GEM hits are used in the
0596 // muon reconstruction, the result should be biased.
0597 // In 12_4_0_pre3, GEM hits are used in the pp scenario, but not in the cosmics scenario.
0598 // https://github.com/cms-sw/cmssw/blob/CMSSW_12_4_0_pre3/RecoMuon/StandAloneMuonProducer/python/standAloneMuons_cfi.py#L111-L127
0599 // https://github.com/cms-sw/cmssw/blob/CMSSW_12_4_0_pre3/RecoMuon/CosmicMuonProducer/python/cosmicMuons_cfi.py
0600 GEMEfficiencyAnalyzer::StartingState GEMEfficiencyAnalyzer::getInnermostMeasurementState(
0601     const reco::TransientTrack& transient_track) {
0602   TrajectoryStateOnSurface state;
0603   DetId det_id;
0604 
0605   const reco::Track& track = transient_track.track();
0606   // get real innermost state
0607   if (isInsideOut(track)) {
0608     state = transient_track.outermostMeasurementState();
0609     det_id = track.outerDetId();
0610 
0611   } else {
0612     state = transient_track.innermostMeasurementState();
0613     det_id = track.innerDetId();
0614   }
0615 
0616   return std::make_tuple(true, state, det_id);
0617 }
0618 
0619 // Use the outermost measurement state as an initial state for the muon propagation.
0620 GEMEfficiencyAnalyzer::StartingState GEMEfficiencyAnalyzer::getOutermostMeasurementState(
0621     const reco::TransientTrack& transient_track) {
0622   const reco::Track& track = transient_track.track();
0623 
0624   TrajectoryStateOnSurface state;
0625   DetId det_id;
0626 
0627   // get real innermost state
0628   if (isInsideOut(track)) {
0629     state = transient_track.innermostMeasurementState();
0630     det_id = track.innerDetId();
0631 
0632   } else {
0633     state = transient_track.outermostMeasurementState();
0634     det_id = track.outerDetId();
0635   }
0636 
0637   return std::make_tuple(true, state, det_id);
0638 }
0639 
0640 // Find the nearest CSC segment to the given GEM layer and then use a trajectory
0641 // state on the surface with the segment as an initial state.
0642 // XXX This method results in the residual phi distribution with two peaks
0643 // because the muon and antimuon make different peaks.
0644 GEMEfficiencyAnalyzer::StartingState GEMEfficiencyAnalyzer::buildStateOnSurfaceWithCSCSegment(
0645     const reco::Muon& muon, const reco::TransientTrack& transient_track, const GEMLayer& gem_layer) {
0646   bool found = false;
0647   TrajectoryStateOnSurface state;
0648   DetId det_id;
0649 
0650   if (const CSCSegment* csc_segment = findCSCSegment(muon, transient_track, gem_layer)) {
0651     const GeomDet* det = muon_service_->trackingGeometry()->idToDet(csc_segment->cscDetId());
0652     const GlobalPoint global_position = det->toGlobal(csc_segment->localPosition());
0653 
0654     found = true;
0655     state = transient_track.stateOnSurface(global_position);
0656     det_id = csc_segment->geographicalId();
0657   }
0658 
0659   return std::make_tuple(found, state, det_id);
0660 }
0661 
0662 // Find an ME11 segment and the build an initial state using the location and
0663 // direction of the ME11 segment. If the muon has an inner track, the outerP of
0664 // the inner track is used as the momentum magnitude. If not, the momentum
0665 // magnitude is set to 1 GeV.
0666 // https://github.com/gem-sw/alignment/blob/713e8fa/GEMCSCBendingAnalyzer/MuonAnalyser/plugins/analyser.cc#L435-L446
0667 GEMEfficiencyAnalyzer::StartingState GEMEfficiencyAnalyzer::buildStartingStateAlignmentStyle(
0668     const reco::Muon& muon, const reco::TransientTrack& transient_track, const GEMLayer& gem_layer) {
0669   bool found = false;
0670   TrajectoryStateOnSurface state;
0671   DetId det_id;
0672 
0673   if (const CSCSegment* csc_segment = findCSCSegment(muon, transient_track, gem_layer)) {
0674     found = true;
0675     det_id = csc_segment->geographicalId();
0676 
0677     // position
0678     const LocalPoint position = csc_segment->localPosition();
0679     // momentum
0680     const reco::TrackRef inner_track = muon.innerTrack();
0681     const float momentum_magnitude = inner_track.isNonnull() ? inner_track.get()->outerP() : 1.0f;
0682     const LocalVector momentum = momentum_magnitude * csc_segment->localDirection();
0683 
0684     // trajectory parameter
0685     const LocalTrajectoryParameters trajectory_parameters{position, momentum, muon.charge()};
0686 
0687     // trajectory error
0688     const LocalTrajectoryError trajectory_error =
0689         asSMatrix<5>(csc_segment->parametersError().similarityT(csc_segment->projectionMatrix()));
0690 
0691     // surface
0692     const Plane& surface = muon_service_->trackingGeometry()->idToDet(det_id)->surface();
0693 
0694     state =
0695         TrajectoryStateOnSurface{trajectory_parameters, trajectory_error, surface, &*muon_service_->magneticField()};
0696   }
0697 
0698   return std::make_tuple(found, state, det_id);
0699 }
0700 
0701 //
0702 // for beam scenario
0703 const CSCSegment* GEMEfficiencyAnalyzer::findCSCSegmentBeam(const reco::TransientTrack& transient_track,
0704                                                             const GEMLayer& gem_layer) {
0705   const CSCSegment* best_csc_segment = nullptr;
0706   double min_z_distance = std::numeric_limits<double>::infinity();  // in cm
0707 
0708   for (trackingRecHit_iterator tracking_rechit_iter = transient_track.recHitsBegin();
0709        tracking_rechit_iter != transient_track.recHitsEnd();
0710        tracking_rechit_iter++) {
0711     const TrackingRecHit* tracking_rechit = *tracking_rechit_iter;
0712     if (not tracking_rechit->isValid()) {
0713       LogDebug(kLogCategory_) << "got an invalid trackingRecHit_iterator from transient_track. skip it.";
0714       continue;
0715     }
0716 
0717     const DetId det_id = tracking_rechit->geographicalId();
0718     if (not MuonHitHelper::isCSC(det_id)) {
0719       continue;
0720     }
0721 
0722     if (tracking_rechit->dimension() != kCSCSegmentDimension_) {
0723       continue;
0724     }
0725 
0726     const CSCDetId csc_id{det_id};
0727     if (not isCSCAllowed(csc_id, gem_layer.id.station())) {
0728       continue;
0729     }
0730 
0731     if (auto csc_segment = dynamic_cast<const CSCSegment*>(tracking_rechit)) {
0732       const GeomDet* det = muon_service_->trackingGeometry()->idToDet(csc_id);
0733       if (det == nullptr) {
0734         edm::LogError(kLogCategory_) << "GlobalTrackingGeometry::idToDet returns nullptr; CSCDetId=" << csc_id;
0735         continue;
0736       }
0737       const GlobalPoint global_position = det->toGlobal(csc_segment->localPosition());
0738       const float z_distance = std::abs(gem_layer.disk->localZclamped(global_position));
0739 
0740       if (z_distance < min_z_distance) {
0741         best_csc_segment = csc_segment;
0742         min_z_distance = z_distance;
0743       }
0744 
0745     } else {
0746       edm::LogError(kLogCategory_)
0747           << "failed to perform the conversion from `const TrackingRechit*` to `const CSCSegment*`";
0748     }
0749   }  // trackingRecHit_iterator
0750 
0751   return best_csc_segment;
0752 }
0753 
0754 const CSCSegment* GEMEfficiencyAnalyzer::findCSCSegmentCosmics(const reco::Muon& muon, const GEMLayer& gem_layer) {
0755   const CSCSegment* best_csc_segment = nullptr;
0756 
0757   for (const reco::MuonChamberMatch& chamber_match : muon.matches()) {
0758     if (not MuonHitHelper::isCSC(chamber_match.id)) {
0759       continue;
0760     }
0761 
0762     const CSCDetId csc_id{chamber_match.id};
0763     if (not isCSCAllowed(csc_id, gem_layer.id.station())) {
0764       continue;
0765     }
0766 
0767     const float x_track = chamber_match.x;
0768     const float y_track = chamber_match.y;
0769 
0770     for (const reco::MuonSegmentMatch& segment_match : chamber_match.segmentMatches) {
0771       if (not segment_match.isMask(reco::MuonSegmentMatch::BestInStationByDR)) {
0772         continue;
0773       }
0774 
0775       const float dr = std::hypot(x_track - segment_match.x, y_track - segment_match.y);
0776       std::cout << kLogCategory_ << ": dr=" << dr << std::endl;
0777 
0778       if (dr > kMuonSegmentMatchDRCut_) {
0779         LogDebug(kLogCategory_) << "too large dR(muon, segment)";
0780         break;
0781       }
0782 
0783       if (segment_match.cscSegmentRef.isNonnull()) {
0784         best_csc_segment = segment_match.cscSegmentRef.get();
0785       }
0786     }  // MuonSegmentMatch
0787   }    // MuonChamberMatch
0788 
0789   return best_csc_segment;
0790 }
0791 
0792 // just thin wrapper
0793 const CSCSegment* GEMEfficiencyAnalyzer::findCSCSegment(const reco::Muon& muon,
0794                                                         const reco::TransientTrack& transient_track,
0795                                                         const GEMLayer& gem_layer) {
0796   if (kScenario_ == ScenarioOption::kCosmics) {
0797     return findCSCSegmentCosmics(muon, gem_layer);
0798   } else {
0799     // pp or HI
0800     return findCSCSegmentBeam(transient_track, gem_layer);
0801   }
0802 }
0803 
0804 bool GEMEfficiencyAnalyzer::isMuonSubdetAllowed(const DetId& det_id, const int gem_station) {
0805   if ((gem_station < 0) or (gem_station > 2)) {
0806     edm::LogError(kLogCategory_) << "got unexpected gem station " << gem_station;
0807     return false;
0808   }
0809 
0810   if (det_id.det() != DetId::Detector::Muon) {
0811     edm::LogError(kLogCategory_) << Form(
0812         "(Detector, Subdetector) = (%d, %d)", static_cast<int>(det_id.det()), det_id.subdetId());
0813     return false;
0814   }
0815 
0816   const std::vector<int> allowed = kMuonSubdetForGEM_.at(gem_station);
0817   return allowed.empty() or (std::find(allowed.begin(), allowed.end(), det_id.subdetId()) != allowed.end());
0818 }
0819 
0820 // Returns a bool value indicating whether or not the CSC detector can be used
0821 // as a start detector for a given GEM station.
0822 // See https://github.com/cms-sw/cmssw/blob/CMSSW_12_4_0_pre3/DataFormats/MuonDetId/interface/CSCDetId.h#L187-L193
0823 // This method is used when using `buildStateOnSurfaceWithCSCSegment` or
0824 // `buildStartingStateAlignmentStyle`
0825 bool GEMEfficiencyAnalyzer::isCSCAllowed(const CSCDetId& csc_id, const int gem_station) {
0826   if ((gem_station < 0) or (gem_station > 2)) {
0827     edm::LogError(kLogCategory_) << "got unexpected gem station " << gem_station;
0828     return false;
0829   }
0830 
0831   // unsigned short to int
0832   const int csc_chamber_type = static_cast<int>(csc_id.iChamberType());
0833 
0834   const std::vector<int> allowed = kCSCForGEM_.at(gem_station);
0835   return allowed.empty() or (std::find(allowed.begin(), allowed.end(), csc_chamber_type) != allowed.end());
0836 }
0837 
0838 bool GEMEfficiencyAnalyzer::checkBounds(const Plane& plane, const GlobalPoint& global_point) {
0839   const LocalPoint local_point = plane.toLocal(global_point);
0840   const LocalPoint local_point_2d(local_point.x(), local_point.y(), 0.0f);
0841   return plane.bounds().inside(local_point_2d);
0842 }
0843 
0844 // TODO comment on the scale
0845 // https://github.com/cms-sw/cmssw/blob/CMSSW_12_0_0_pre3/DataFormats/GeometrySurface/src/SimpleDiskBounds.cc#L20-L35
0846 bool GEMEfficiencyAnalyzer::checkBounds(const Plane& plane,
0847                                         const GlobalPoint& global_point,
0848                                         const GlobalError& global_error,
0849                                         const float scale) {
0850   const LocalPoint local_point = plane.toLocal(global_point);
0851   const LocalError local_error = ErrorFrameTransformer::transform(global_error, plane);
0852 
0853   const LocalPoint local_point_2d{local_point.x(), local_point.y(), 0.0f};
0854   return plane.bounds().inside(local_point_2d, local_error, scale);
0855 }
0856 
0857 const GEMEtaPartition* GEMEfficiencyAnalyzer::findEtaPartition(const GlobalPoint& global_point,
0858                                                                const GlobalError& global_error,
0859                                                                const std::vector<const GEMChamber*>& chamber_vector) {
0860   const GEMEtaPartition* bound = nullptr;
0861   for (const GEMChamber* chamber : chamber_vector) {
0862     if (not checkBounds(chamber->surface(), global_point, global_error, kBoundsErrorScale_)) {
0863       continue;
0864     }
0865 
0866     for (const GEMEtaPartition* eta_partition : chamber->etaPartitions()) {
0867       if (checkBounds(eta_partition->surface(), global_point, global_error, kBoundsErrorScale_)) {
0868         bound = eta_partition;
0869         break;
0870       }
0871     }  // GEMEtaPartition
0872   }    // GEMChamber
0873 
0874   return bound;
0875 }
0876 
0877 // Borrowed from https://github.com/gem-sw/alignment/blob/713e8fa/GEMCSCBendingAnalyzer/MuonAnalyser/plugins/analyser.cc#L321-L327
0878 float GEMEfficiencyAnalyzer::computeRdPhi(const GlobalPoint& prop_global_pos,
0879                                           const LocalPoint& hit_local_pos,
0880                                           const GEMEtaPartition* eta_partition) {
0881   const StripTopology& topology = eta_partition->specificTopology();
0882   const LocalPoint prop_local_pos = eta_partition->toLocal(prop_global_pos);
0883 
0884   const float dx = prop_local_pos.x() - hit_local_pos.x();
0885   const float dy = prop_local_pos.y() - hit_local_pos.y();
0886   const float hit_strip = eta_partition->strip(hit_local_pos);
0887   const float hit_phi = topology.stripAngle(hit_strip);
0888   const float rdphi = std::cos(hit_phi) * dx - std::sin(hit_phi) * dy;
0889   return rdphi;
0890 }
0891 
0892 // Returns a global delta phi between a propagated muon and a reconstructed hit.
0893 float GEMEfficiencyAnalyzer::computeDeltaPhi(const GlobalPoint& prop_global_pos,
0894                                              const LocalPoint& hit_local_pos,
0895                                              const GEMEtaPartition* eta_partition) {
0896   const GlobalPoint hit_global_pos = eta_partition->toGlobal(hit_local_pos);
0897   const float dphi = Geom::convertRadToDeg(prop_global_pos.phi() - hit_global_pos.phi());
0898   return dphi;
0899 }
0900 
0901 // a thin wrapper to hide a messy conditional statement
0902 float GEMEfficiencyAnalyzer::computeMatchingMetric(const GlobalPoint& prop_global_pos,
0903                                                    const LocalPoint& hit_local_pos,
0904                                                    const GEMEtaPartition* eta_partition) {
0905   float metric;
0906   switch (kMatchingMetric_) {
0907     case MatchingMetric::kDeltaPhi: {
0908       metric = computeDeltaPhi(prop_global_pos, hit_local_pos, eta_partition);
0909       break;
0910     }
0911     case MatchingMetric::kRdPhi: {
0912       metric = computeRdPhi(prop_global_pos, hit_local_pos, eta_partition);
0913       break;
0914     }
0915     default: {
0916       edm::LogError(kLogCategory_) << "unknown MatchingMetric.";  // TODO
0917       metric = std::numeric_limits<float>::infinity();
0918     }
0919   }
0920 
0921   return metric;
0922 }
0923 
0924 // This method finds the closest hit to a propagated muon in the eta partition
0925 // with that propagated muon. Adjacent eta partitions are excluded from the area
0926 // of interst to avoid ambiguity in defining the detection efficiency of each
0927 // eta partition.
0928 std::pair<const GEMRecHit*, float> GEMEfficiencyAnalyzer::findClosestHit(const GlobalPoint& prop_global_pos,
0929                                                                          const GEMRecHitCollection::range& rechit_range,
0930                                                                          const GEMEtaPartition* eta_partition) {
0931   const GEMRecHit* closest_hit = nullptr;
0932   float min_metric = std::numeric_limits<float>::infinity();
0933 
0934   for (auto hit = rechit_range.first; hit != rechit_range.second; ++hit) {
0935     const LocalPoint hit_local_pos = hit->localPosition();
0936 
0937     const float metric = computeMatchingMetric(prop_global_pos, hit_local_pos, eta_partition);
0938 
0939     if (std::abs(metric) < std::abs(min_metric)) {
0940       min_metric = metric;
0941       closest_hit = &(*hit);
0942     }
0943   }
0944 
0945   return std::make_pair(closest_hit, min_metric);
0946 }
0947 
0948 void GEMEfficiencyAnalyzer::dqmBeginRun(edm::Run const&, edm::EventSetup const& setup) {
0949   const GEMGeometry* gem = nullptr;
0950   if (auto handle = setup.getHandle(kGEMGeometryTokenBeginRun_)) {
0951     gem = handle.product();
0952   } else {
0953     edm::LogError(kLogCategory_ + "|dqmBeginRun") << "failed to get GEMGeometry";
0954     return;
0955   }
0956 
0957   buildGEMLayers(gem);
0958 }
0959 
0960 void GEMEfficiencyAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0961   //////////////////////////////////////////////////////////////////////////////
0962   // get data from Event
0963   //////////////////////////////////////////////////////////////////////////////
0964   const GEMRecHitCollection* rechit_collection = nullptr;
0965   if (auto handle = event.getHandle(kGEMRecHitCollectionToken_)) {
0966     rechit_collection = handle.product();
0967   } else {
0968     edm::LogError(kLogCategory_) << "failed to get GEMRecHitCollection";
0969     return;
0970   }
0971 
0972   const edm::View<reco::Muon>* muon_view = nullptr;
0973   if (auto handle = event.getHandle(kMuonViewToken_)) {
0974     muon_view = handle.product();
0975   } else {
0976     edm::LogError(kLogCategory_) << "failed to get View<Muon>";
0977     return;
0978   }
0979 
0980   const GEMOHStatusCollection* oh_status_collection = nullptr;
0981   const GEMVFATStatusCollection* vfat_status_collection = nullptr;
0982   if (kMaskChamberWithError_) {
0983     if (auto handle = event.getHandle(kGEMOHStatusCollectionToken_)) {
0984       oh_status_collection = handle.product();
0985     } else {
0986       edm::LogError(kLogCategory_) << "failed to get OHVFATStatusCollection";
0987       return;
0988     }
0989 
0990     if (auto handle = event.getHandle(kGEMVFATStatusCollectionToken_)) {
0991       vfat_status_collection = handle.product();
0992     } else {
0993       edm::LogError(kLogCategory_) << "failed to get GEMVFATStatusCollection";
0994       return;
0995     }
0996   }
0997 
0998   //////////////////////////////////////////////////////////////////////////////
0999   // get data from EventSetup
1000   //////////////////////////////////////////////////////////////////////////////
1001   const TransientTrackBuilder* transient_track_builder = nullptr;
1002   if (auto handle = setup.getHandle(kTransientTrackBuilderToken_)) {
1003     transient_track_builder = handle.product();
1004   } else {
1005     edm::LogError(kLogCategory_) << "failed to get TransientTrackBuilder";
1006     return;
1007   }
1008 
1009   //////////////////////////////////////////////////////////////////////////////
1010   // get more data from EventSetup using MuonServiceProxy
1011   //////////////////////////////////////////////////////////////////////////////
1012   muon_service_->update(setup);
1013 
1014   // TODO StraightLinePropagator if B < epsilon else SteppingHelixPropagatorAny
1015   const Propagator* propagator = nullptr;
1016   if (auto handle = muon_service_->propagator("SteppingHelixPropagatorAny")) {
1017     propagator = handle.product();
1018   } else {
1019     edm::LogError(kLogCategory_) << "failed to get Propagator";
1020     return;
1021   }
1022 
1023   //////////////////////////////////////////////////////////////////////////////
1024   //  Main loop
1025   //////////////////////////////////////////////////////////////////////////////
1026 
1027   for (const reco::Muon& muon : *muon_view) {
1028     const reco::Track* track = muon.muonTrack(kMuonTrackType_).get();
1029     if (track == nullptr) {
1030       LogDebug(kLogCategory_) << "failed to get a " << kMuonTrackTypeName_;
1031       continue;
1032     }
1033 
1034     const reco::TransientTrack transient_track = transient_track_builder->build(track);
1035     if (not transient_track.isValid()) {
1036       edm::LogError(kLogCategory_) << "failed to build TransientTrack";
1037       continue;
1038     }
1039 
1040     for (const GEMLayer& layer : gem_layers_) {
1041       if (checkPropagationDirection(track, layer)) {
1042         LogDebug(kLogCategory_) << "bad flight path. skip this propagation.";
1043         continue;
1044       }
1045 
1046       const auto [found_start_state, start_state, start_id] = buildStartingState(muon, transient_track, layer);
1047       if (not found_start_state) {
1048         LogDebug(kLogCategory_) << "propagation starting state not found";
1049         continue;
1050       }
1051 
1052       // the trajectory state on the destination surface
1053       const auto [propagated_state, prop_path_length] = propagator->propagateWithPath(start_state, *(layer.disk));
1054       if (not propagated_state.isValid()) {
1055         LogDebug(kLogCategory_) << "failed to propagate a muon from "
1056                                 << Form("(Detector, Subdetector) = (%d, %d)",
1057                                         static_cast<int>(start_id.det()),
1058                                         start_id.subdetId())
1059                                 << " to " << layer.id << ". The path length is " << prop_path_length;
1060         continue;
1061       }
1062 
1063       const GlobalPoint prop_global_pos = propagated_state.globalPosition();
1064       const GlobalError& prop_global_err =
1065           ErrorFrameTransformer::transform(propagated_state.localError().positionError(), *layer.disk);
1066 
1067       if (not checkBounds(*layer.disk, prop_global_pos, prop_global_err, kBoundsErrorScale_)) {
1068         LogDebug(kLogCategory_) << "failed to pass checkBounds";
1069         continue;
1070       }
1071 
1072       const GEMEtaPartition* eta_partition = findEtaPartition(prop_global_pos, prop_global_err, layer.chambers);
1073       if (eta_partition == nullptr) {
1074         LogDebug(kLogCategory_) << "failed to find an eta partition";
1075         continue;
1076       }
1077 
1078       const GEMDetId gem_id = eta_partition->id();
1079 
1080       if (kMaskChamberWithError_) {
1081         const bool has_error = maskChamberWithError(gem_id.chamberId(), oh_status_collection, vfat_status_collection);
1082         if (has_error) {
1083           LogDebug(kLogCategory_) << gem_id.chamberId() << " has an erorr. Skip this propagation.";
1084           continue;
1085         }
1086       }
1087 
1088       //////////////////////////////////////////////////////////////////////////
1089       //
1090       //////////////////////////////////////////////////////////////////////////
1091       const GEMDetId rs_key = getReStKey(gem_id);     // region-station
1092       const GEMDetId rsl_key = getReStLaKey(gem_id);  // region-station-layer
1093       const GEMDetId rse_key = getReStEtKey(gem_id);  // region-station-ieta
1094 
1095       const int station_id = gem_id.station();
1096       const int chamber_id = gem_id.chamber();
1097       const int ieta = gem_id.ieta();
1098 
1099       const double muon_pt = muon.pt();
1100       const double muon_eta = std::fabs(muon.eta());
1101       const double muon_phi = Geom::convertRadToDeg(muon.phi());
1102 
1103       const double prop_global_err_r = std::sqrt(prop_global_err.rerr(prop_global_pos));
1104       const double prop_global_err_phi = Geom::convertRadToDeg(std::sqrt(prop_global_err.phierr(prop_global_pos)));
1105 
1106       // cuts
1107       const bool passed_prop_err_r_cut = (prop_global_err_r < kPropagationErrorRCut_);
1108       const bool passed_prop_err_phi_cut = (prop_global_err_phi < kPropagationErrorPhiCut_);
1109       const bool passed_pt_cut = muon_pt > kMuonPtMinCuts_.at(station_id);
1110       const bool passed_eta_cut =
1111           (muon_eta > kMuonEtaMinCuts_.at(station_id)) and (muon_eta < kMuonEtaMaxCuts_.at(station_id));
1112 
1113       const bool passed_prop_err_cuts = passed_prop_err_r_cut and passed_prop_err_phi_cut;
1114       const bool passed_all_cuts = passed_prop_err_cuts and passed_pt_cut and passed_eta_cut;
1115 
1116       const int cutflow_last = not kModeDev_                 ? 0
1117                                : not passed_prop_err_r_cut   ? 1
1118                                : not passed_prop_err_phi_cut ? 2
1119                                : not passed_pt_cut           ? 3
1120                                : not passed_eta_cut          ? 4
1121                                                              : 5;
1122 
1123       //////////////////////////////////////////////////////////////////////////
1124       // Fill denominators
1125       //////////////////////////////////////////////////////////////////////////
1126       if (passed_eta_cut and passed_prop_err_cuts) {
1127         fillMEWithinLimits(me_muon_pt_, rs_key, muon_pt);
1128       }
1129 
1130       if (passed_pt_cut and passed_prop_err_cuts) {
1131         fillMEWithinLimits(me_muon_eta_, rs_key, muon_eta);
1132       }
1133 
1134       if (passed_all_cuts) {
1135         fillME(me_chamber_ieta_, rsl_key, chamber_id, ieta);
1136         fillME(me_muon_phi_, rs_key, muon_phi);
1137       }
1138 
1139       if (kModeDev_) {
1140         fillMEWithinLimits(me_prop_path_length_, rsl_key, prop_path_length);
1141         fillMEWithinLimits(me_prop_err_r_, rsl_key, prop_global_err_r);
1142         fillMEWithinLimits(me_prop_err_phi_, rsl_key, prop_global_err_phi);
1143 
1144         fillMEWithinLimits(me_muon_pt_all_, rs_key, muon_pt);
1145         fillMEWithinLimits(me_muon_eta_all_, rs_key, muon_eta);
1146 
1147         fillME(me_muon_charge_, rs_key, muon.charge());
1148 
1149         for (int bin = 1; bin <= cutflow_last; bin++) {
1150           fillME(me_cutflow_, rsl_key, bin);
1151         }
1152 
1153       }  // dev mode
1154 
1155       //////////////////////////////////////////////////////////////////////////
1156       // Find a closet hit
1157       //////////////////////////////////////////////////////////////////////////
1158       const auto [hit, matching_metric] =
1159           findClosestHit(prop_global_pos, rechit_collection->get(gem_id), eta_partition);
1160 
1161       if (hit == nullptr) {
1162         LogDebug(kLogCategory_) << "hit not found";
1163         continue;
1164       }
1165 
1166       if (kModeDev_) {
1167         fillMEWithinLimits(me_matching_metric_all_, rse_key, matching_metric);
1168       }
1169 
1170       if (std::abs(matching_metric) > kMatchingCut_) {
1171         LogDebug(kLogCategory_) << "failed to pass the residual rphi cut";
1172         continue;
1173       }
1174 
1175       //////////////////////////////////////////////////////////////////////////
1176       // Fill numerators
1177       //////////////////////////////////////////////////////////////////////////
1178       if (passed_eta_cut and passed_prop_err_cuts) {
1179         fillMEWithinLimits(me_muon_pt_matched_, rs_key, muon_pt);
1180       }
1181 
1182       if (passed_pt_cut and passed_prop_err_cuts) {
1183         fillMEWithinLimits(me_muon_eta_matched_, rs_key, muon_eta);
1184       }
1185 
1186       if (passed_all_cuts) {
1187         fillME(me_chamber_ieta_matched_, rsl_key, chamber_id, ieta);
1188         fillME(me_muon_phi_matched_, rs_key, muon_phi);
1189       }
1190 
1191       if (kModeDev_) {
1192         fillMEWithinLimits(me_prop_path_length_matched_, rsl_key, prop_path_length);
1193 
1194         fillMEWithinLimits(me_prop_err_r_matched_, rsl_key, prop_global_err_r);
1195         fillMEWithinLimits(me_prop_err_phi_matched_, rsl_key, prop_global_err_phi);
1196 
1197         fillMEWithinLimits(me_muon_pt_all_matched_, rs_key, muon_pt);
1198         fillMEWithinLimits(me_muon_eta_all_matched_, rs_key, muon_eta);
1199 
1200         fillME(me_muon_charge_matched_, rs_key, muon.charge());
1201 
1202         if (passed_all_cuts) {
1203           for (int bin = 1; bin <= cutflow_last; bin++) {
1204             fillME(me_cutflow_matched_, rsl_key, bin);
1205           }
1206         }
1207       }
1208 
1209       //////////////////////////////////////////////////////////////////////////
1210       // Fill resolutions
1211       //////////////////////////////////////////////////////////////////////////
1212       if (passed_all_cuts) {
1213         const LocalPoint hit_local_pos = hit->localPosition();
1214         const GlobalPoint& hit_global_pos = eta_partition->toGlobal(hit_local_pos);
1215         const float residual_phi = Geom::convertRadToDeg(prop_global_pos.phi() - hit_global_pos.phi());
1216 
1217         fillMEWithinLimits(me_residual_phi_, rse_key, residual_phi);
1218 
1219         if (kModeDev_) {
1220           const LocalPoint prop_local_pos = eta_partition->toLocal(prop_global_pos);
1221           const StripTopology& topology = eta_partition->specificTopology();
1222 
1223           const float residual_x = prop_local_pos.x() - hit_local_pos.x();
1224           const float residual_y = prop_local_pos.y() - hit_local_pos.y();
1225           const float residual_strip = topology.strip(prop_local_pos) - topology.strip(hit_local_pos);
1226 
1227           fillMEWithinLimits(me_matching_metric_, rse_key, matching_metric);
1228           fillMEWithinLimits(me_residual_x_, rse_key, residual_x);
1229           fillMEWithinLimits(me_residual_y_, rse_key, residual_y);
1230           fillMEWithinLimits(me_residual_strip_, rse_key, residual_strip);
1231 
1232           if (muon.charge() < 0) {
1233             fillMEWithinLimits(me_residual_phi_muon_, rse_key, residual_phi);
1234           } else {
1235             fillMEWithinLimits(me_residual_phi_antimuon_, rse_key, residual_phi);
1236           }
1237         }  // kModeDev_
1238       }    // passed_all_cuts
1239     }      // destination
1240   }        // Muon
1241 }  // analyze
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!