diff --git a/DQMOffline/MuonDPG/plugins/GEMTnPEfficiencyTask.cc b/DQMOffline/MuonDPG/plugins/GEMTnPEfficiencyTask.cc new file mode 100644 index 0000000000000..170ec3d988060 --- /dev/null +++ b/DQMOffline/MuonDPG/plugins/GEMTnPEfficiencyTask.cc @@ -0,0 +1,1030 @@ +/* + * \file GEMTnPEfficiencyTask.cc + * \author Qianying + * + * \interited from the TnP framework of + * \author L. Lunerti - INFN Bologna + * + */ + +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "DataFormats/MuonReco/interface/MuonSegmentMatch.h" +#include "DataFormats/MuonReco/interface/MuonGEMHitMatch.h" + +#include "DQMOffline/MuonDPG/interface/BaseTnPEfficiencyTask.h" + +class GEMTnPEfficiencyTask : public BaseTnPEfficiencyTask { +public: + /// Constructor + GEMTnPEfficiencyTask(const edm::ParameterSet& config); + + /// Destructor + ~GEMTnPEfficiencyTask() override; + +protected: + std::string topFolder() const override; + + void bookHistograms(DQMStore::IBooker& iBooker, edm::Run const& run, edm::EventSetup const& context) override; + + /// Analyze + void analyze(const edm::Event& event, const edm::EventSetup& context) override; +}; + +GEMTnPEfficiencyTask::GEMTnPEfficiencyTask(const edm::ParameterSet& config) : BaseTnPEfficiencyTask(config) { + LogTrace("DQMOffline|MuonDPG|GEMTnPEfficiencyTask") << "[GEMTnPEfficiencyTask]: Constructor" << std::endl; +} + +GEMTnPEfficiencyTask::~GEMTnPEfficiencyTask() { + LogTrace("DQMOffline|MuonDPG|GEMTnPEfficiencyTask") + << "[GEMTnPEfficiencyTask]: analyzed " << m_nEvents << " events" << std::endl; +} + +void GEMTnPEfficiencyTask::bookHistograms(DQMStore::IBooker& iBooker, + edm::Run const& run, + edm::EventSetup const& context) { + BaseTnPEfficiencyTask::bookHistograms(iBooker, run, context); + + LogTrace("DQMOffline|MuonDPG|GEMTnPEfficiencyTask") << "[GEMTnPEfficiencyTask]: bookHistograms" << std::endl; + + auto baseDir = topFolder() + "Task/"; + iBooker.setCurrentFolder(baseDir); + + MonitorElement* me_GE11_pass_Ch_region = + iBooker.book2D("GE11_nPassingProbe_Ch_region", "GE11_nPassingProbe_Ch_region", 2, -1.5, 1.5, 36, 1, 37); + MonitorElement* me_GE11_fail_Ch_region = + iBooker.book2D("GE11_nFailingProbe_Ch_region", "GE11_nFailingProbe_Ch_region", 2, -1.5, 1.5, 36, 1, 37); + MonitorElement* me_GE21_pass_Ch_region = + iBooker.book2D("GE21_nPassingProbe_Ch_region", "GE21_nPassingProbe_Ch_region", 2, -1.5, 1.5, 36, 1, 37); + MonitorElement* me_GE21_fail_Ch_region = + iBooker.book2D("GE21_nFailingProbe_Ch_region", "GE21_nFailingProbe_Ch_region", 2, -1.5, 1.5, 36, 1, 37); + MonitorElement* me_GEM_pass_Ch_region_GE1 = + iBooker.book2D("GEM_nPassingProbe_Ch_region_GE1", "GEM_nPassingProbe_Ch_region_GE1", 4, 0, 4, 36, 1, 37); + MonitorElement* me_GEM_fail_Ch_region_GE1 = + iBooker.book2D("GEM_nFailingProbe_Ch_region_GE1", "GEM_nFailingProbe_Ch_region_GE1", 4, 0, 4, 36, 1, 37); + MonitorElement* me_GEM_pass_Ch_region_GE1_NoL = + iBooker.book2D("GEM_nPassingProbe_Ch_region_GE1_NoL", "GEM_nPassingProbe_Ch_region_GE1_NoL", 2, 0, 2, 36, 1, 37); + MonitorElement* me_GEM_fail_Ch_region_GE1_NoL = + iBooker.book2D("GEM_nFailingProbe_Ch_region_GE1_NoL", "GEM_nFailingProbe_Ch_region_GE1_NoL", 2, 0, 2, 36, 1, 37); + MonitorElement* me_GE11_pass_Ch_ieta = + iBooker.book2D("GE11_nPassingProbe_Ch_ieta", "GE11_nPassingProbe_Ch_ieta", 8, 1, 9, 36, 1, 37); + MonitorElement* me_GE11_fail_Ch_ieta = + iBooker.book2D("GE11_nFailingProbe_Ch_ieta", "GE11_nFailingProbe_Ch_ieta", 8, 1, 9, 36, 1, 37); + MonitorElement* me_GE11_pass_Ch_phi = iBooker.book2D( + "GE11_nPassingProbe_Ch_phi", "GE11_nPassingProbe_Ch_phi", 20, -TMath::Pi(), TMath::Pi(), 36, 1, 37); + MonitorElement* me_GE11_fail_Ch_phi = iBooker.book2D( + "GE11_nFailingProbe_Ch_phi", "GE11_nFailingProbe_Ch_phi", 20, -TMath::Pi(), TMath::Pi(), 36, 1, 37); + MonitorElement* me_GE11_pass_allCh_1D = + iBooker.book1D("GE11_nPassingProbe_allCh_1D", "GE11_nPassingProbe_allCh_1D", 2, -1.5, 1.5); + MonitorElement* me_GE11_fail_allCh_1D = + iBooker.book1D("GE11_nFailingProbe_allCh_1D", "GE11_nFailingProbe_allCh_1D", 2, -1.5, 1.5); + MonitorElement* me_GE11_pass_chamber_1D = + iBooker.book1D("GE11_nPassingProbe_chamber_1D", "GE11_nPassingProbe_chamber_1D", 36, 1, 37); + MonitorElement* me_GE11_fail_chamber_1D = + iBooker.book1D("GE11_nFailingProbe_chamber_1D", "GE11_nFailingProbe_chamber_1D", 36, 1, 37); + MonitorElement* me_GE21_pass_Ch_ieta = + iBooker.book2D("GE21_nPassingProbe_Ch_ieta", "GE21_nPassingProbe_Ch_ieta", 16, 1, 17, 18, 1, 19); + MonitorElement* me_GE21_fail_Ch_ieta = + iBooker.book2D("GE21_nFailingProbe_Ch_ieta", "GE21_nFailingProbe_Ch_ieta", 16, 1, 17, 18, 1, 19); + MonitorElement* me_GE21_pass_Ch_phi = iBooker.book2D( + "GE21_nPassingProbe_Ch_phi", "GE21_nPassingProbe_Ch_phi", 20, -TMath::Pi(), TMath::Pi(), 18, 1, 19); + MonitorElement* me_GE21_fail_Ch_phi = iBooker.book2D( + "GE21_nFailingProbe_Ch_phi", "GE21_nFailingProbe_Ch_phi", 20, -TMath::Pi(), TMath::Pi(), 18, 1, 19); + MonitorElement* me_GE21_pass_allCh_1D = + iBooker.book1D("GE21_nPassingProbe_allCh_1D", "GE21_nPassingProbe_allCh_1D", 2, -1.5, 1.5); + MonitorElement* me_GE21_fail_allCh_1D = + iBooker.book1D("GE21_nFailingProbe_allCh_1D", "GE21_nFailingProbe_allCh_1D", 2, -1.5, 1.5); + MonitorElement* me_GE21_pass_chamber_1D = + iBooker.book1D("GE21_nPassingProbe_chamber_1D", "GE21_nPassingProbe_chamber_1D", 18, 1, 19); + MonitorElement* me_GE21_fail_chamber_1D = + iBooker.book1D("GE21_nFailingProbe_chamber_1D", "GE21_nFailingProbe_chamber_1D", 18, 1, 19); + MonitorElement* me_GEM_pass_chamber_p1_1D = + iBooker.book1D("GEM_nPassingProbe_chamber_p1_1D", "GEM_nPassingProbe_chamber_p1_1D", 36, 1, 37); + MonitorElement* me_GEM_fail_chamber_p1_1D = + iBooker.book1D("GEM_nFailingProbe_chamber_p1_1D", "GEM_nFailingProbe_chamber_p1_1D", 36, 1, 37); + MonitorElement* me_GEM_pass_chamber_p2_1D = + iBooker.book1D("GEM_nPassingProbe_chamber_p2_1D", "GEM_nPassingProbe_chamber_p2_1D", 36, 1, 37); + MonitorElement* me_GEM_fail_chamber_p2_1D = + iBooker.book1D("GEM_nFailingProbe_chamber_p2_1D", "GEM_nFailingProbe_chamber_p2_1D", 36, 1, 37); + MonitorElement* me_GEM_pass_chamber_n1_1D = + iBooker.book1D("GEM_nPassingProbe_chamber_n1_1D", "GEM_nPassingProbe_chamber_n1_1D", 36, 1, 37); + MonitorElement* me_GEM_fail_chamber_n1_1D = + iBooker.book1D("GEM_nFailingProbe_chamber_n1_1D", "GEM_nFailingProbe_chamber_n1_1D", 36, 1, 37); + MonitorElement* me_GEM_pass_chamber_n2_1D = + iBooker.book1D("GEM_nPassingProbe_chamber_n2_1D", "GEM_nPassingProbe_chamber_n2_1D", 36, 1, 37); + MonitorElement* me_GEM_fail_chamber_n2_1D = + iBooker.book1D("GEM_nFailingProbe_chamber_n2_1D", "GEM_nFailingProbe_chamber_n2_1D", 36, 1, 37); + // + MonitorElement* me_GEM_pass_pt_1D = iBooker.book1D("GEM_nPassingProbe_pt_1D", "GEM_nPassingProbe_pt_1D", 20, 0, 100); + MonitorElement* me_GEM_fail_pt_1D = iBooker.book1D("GEM_nFailingProbe_pt_1D", "GEM_nFailingProbe_pt_1D", 20, 0, 100); + MonitorElement* me_GEM_pass_eta_1D = + iBooker.book1D("GEM_nPassingProbe_eta_1D", "GEM_nPassingProbe_eta_1D", 24, 0, 2.4); + MonitorElement* me_GEM_fail_eta_1D = + iBooker.book1D("GEM_nFailingProbe_eta_1D", "GEM_nFailingProbe_eta_1D", 24, 0, 2.4); + MonitorElement* me_GEM_pass_phi_1D = + iBooker.book1D("GEM_nPassingProbe_phi_1D", "GEM_nPassingProbe_phi_1D", 20, -TMath::Pi(), TMath::Pi()); + MonitorElement* me_GEM_fail_phi_1D = + iBooker.book1D("GEM_nFailingProbe_phi_1D", "GEM_nFailingProbe_phi_1D", 20, -TMath::Pi(), TMath::Pi()); + /// + MonitorElement* me_GEM_pass_pt_p1_1D = + iBooker.book1D("GEM_nPassingProbe_pt_p1_1D", "GEM_nPassingProbe_pt_p1_1D", 20, 0, 100); + MonitorElement* me_GEM_fail_pt_p1_1D = + iBooker.book1D("GEM_nFailingProbe_pt_p1_1D", "GEM_nFailingProbe_pt_p1_1D", 20, 0, 100); + MonitorElement* me_GEM_pass_eta_p1_1D = + iBooker.book1D("GEM_nPassingProbe_eta_p1_1D", "GEM_nPassingProbe_eta_p1_1D", 24, 0, 2.4); + MonitorElement* me_GEM_fail_eta_p1_1D = + iBooker.book1D("GEM_nFailingProbe_eta_p1_1D", "GEM_nFailingProbe_eta_p1_1D", 24, 0, 2.4); + MonitorElement* me_GEM_pass_phi_p1_1D = + iBooker.book1D("GEM_nPassingProbe_phi_p1_1D", "GEM_nPassingProbe_phi_p1_1D", 20, -TMath::Pi(), TMath::Pi()); + MonitorElement* me_GEM_fail_phi_p1_1D = + iBooker.book1D("GEM_nFailingProbe_phi_p1_1D", "GEM_nFailingProbe_phi_p1_1D", 20, -TMath::Pi(), TMath::Pi()); + MonitorElement* me_GEM_pass_pt_p2_1D = + iBooker.book1D("GEM_nPassingProbe_pt_p2_1D", "GEM_nPassingProbe_pt_p2_1D", 20, 0, 100); + MonitorElement* me_GEM_fail_pt_p2_1D = + iBooker.book1D("GEM_nFailingProbe_pt_p2_1D", "GEM_nFailingProbe_pt_p2_1D", 20, 0, 100); + MonitorElement* me_GEM_pass_eta_p2_1D = + iBooker.book1D("GEM_nPassingProbe_eta_p2_1D", "GEM_nPassingProbe_eta_p2_1D", 24, 0, 2.4); + MonitorElement* me_GEM_fail_eta_p2_1D = + iBooker.book1D("GEM_nFailingProbe_eta_p2_1D", "GEM_nFailingProbe_eta_p2_1D", 24, 0, 2.4); + MonitorElement* me_GEM_pass_phi_p2_1D = + iBooker.book1D("GEM_nPassingProbe_phi_p2_1D", "GEM_nPassingProbe_phi_p2_1D", 20, -TMath::Pi(), TMath::Pi()); + MonitorElement* me_GEM_fail_phi_p2_1D = + iBooker.book1D("GEM_nFailingProbe_phi_p2_1D", "GEM_nFailingProbe_phi_p2_1D", 20, -TMath::Pi(), TMath::Pi()); + MonitorElement* me_GEM_pass_pt_n1_1D = + iBooker.book1D("GEM_nPassingProbe_pt_n1_1D", "GEM_nPassingProbe_pt_n1_1D", 20, 0, 100); + MonitorElement* me_GEM_fail_pt_n1_1D = + iBooker.book1D("GEM_nFailingProbe_pt_n1_1D", "GEM_nFailingProbe_pt_n1_1D", 20, 0, 100); + MonitorElement* me_GEM_pass_eta_n1_1D = + iBooker.book1D("GEM_nPassingProbe_eta_n1_1D", "GEM_nPassingProbe_eta_n1_1D", 24, 0, 2.4); + MonitorElement* me_GEM_fail_eta_n1_1D = + iBooker.book1D("GEM_nFailingProbe_eta_n1_1D", "GEM_nFailingProbe_eta_n1_1D", 24, 0, 2.4); + MonitorElement* me_GEM_pass_phi_n1_1D = + iBooker.book1D("GEM_nPassingProbe_phi_n1_1D", "GEM_nPassingProbe_phi_n1_1D", 20, -TMath::Pi(), TMath::Pi()); + MonitorElement* me_GEM_fail_phi_n1_1D = + iBooker.book1D("GEM_nFailingProbe_phi_n1_1D", "GEM_nFailingProbe_phi_n1_1D", 20, -TMath::Pi(), TMath::Pi()); + MonitorElement* me_GEM_pass_pt_n2_1D = + iBooker.book1D("GEM_nPassingProbe_pt_n2_1D", "GEM_nPassingProbe_pt_n2_1D", 20, 0, 100); + MonitorElement* me_GEM_fail_pt_n2_1D = + iBooker.book1D("GEM_nFailingProbe_pt_n2_1D", "GEM_nFailingProbe_pt_n2_1D", 20, 0, 100); + MonitorElement* me_GEM_pass_eta_n2_1D = + iBooker.book1D("GEM_nPassingProbe_eta_n2_1D", "GEM_nPassingProbe_eta_n2_1D", 24, 0, 2.4); + MonitorElement* me_GEM_fail_eta_n2_1D = + iBooker.book1D("GEM_nFailingProbe_eta_n2_1D", "GEM_nFailingProbe_eta_n2_1D", 24, 0, 2.4); + MonitorElement* me_GEM_pass_phi_n2_1D = + iBooker.book1D("GEM_nPassingProbe_phi_n2_1D", "GEM_nPassingProbe_phi_n2_1D", 20, -TMath::Pi(), TMath::Pi()); + MonitorElement* me_GEM_fail_phi_n2_1D = + iBooker.book1D("GEM_nFailingProbe_phi_n2_1D", "GEM_nFailingProbe_phi_n2_1D", 20, -TMath::Pi(), TMath::Pi()); + //// + MonitorElement* me_ME0_pass_chamber_1D = + iBooker.book1D("ME0_nPassingProbe_chamber_1D", "ME0_nPassingProbe_chamber_1D", 18, 1, 19); + MonitorElement* me_ME0_fail_chamber_1D = + iBooker.book1D("ME0_nFailingProbe_chamber_1D", "ME0_nFailingProbe_chamber_1D", 18, 1, 19); + MonitorElement* me_GEM_pass_Ch_region_layer_phase2 = iBooker.book2D( + "GEM_nPassingProbe_Ch_region_layer_phase2", "GEM_nPassingProbe_Ch_region_layer_phase2", 10, 0, 10, 36, 1, 37); + MonitorElement* me_GEM_fail_Ch_region_layer_phase2 = iBooker.book2D( + "GEM_nFailingProbe_Ch_region_layer_phase2", "GEM_nFailingProbe_Ch_region_layer_phase2", 10, 0, 10, 36, 1, 37); + + me_GE11_pass_allCh_1D->setBinLabel(1, "GE-11", 1); + me_GE11_pass_allCh_1D->setBinLabel(2, "GE+11", 1); + me_GE11_pass_allCh_1D->setAxisTitle("Number of passing probes", 2); + + me_GE11_fail_allCh_1D->setBinLabel(1, "GE-11", 1); + me_GE11_fail_allCh_1D->setBinLabel(2, "GE+11", 1); + me_GE11_fail_allCh_1D->setAxisTitle("Number of failing probes", 2); + + me_GE11_pass_chamber_1D->setAxisTitle("Chamber", 1); + me_GE11_pass_chamber_1D->setAxisTitle("Number of passing probes", 2); + me_GE11_fail_chamber_1D->setAxisTitle("Chamber", 1); + me_GE11_fail_chamber_1D->setAxisTitle("Number of failing probes", 2); + + me_GE21_pass_allCh_1D->setBinLabel(1, "GE-21", 1); + me_GE21_pass_allCh_1D->setBinLabel(2, "GE+21", 1); + me_GE21_pass_allCh_1D->setAxisTitle("Number of passing probes", 2); + + me_GE21_fail_allCh_1D->setBinLabel(1, "GE-21", 1); + me_GE21_fail_allCh_1D->setBinLabel(2, "GE+21", 1); + me_GE21_fail_allCh_1D->setAxisTitle("Number of failing probes", 2); + + me_GE21_pass_chamber_1D->setAxisTitle("Chamber", 1); + me_GE21_pass_chamber_1D->setAxisTitle("Number of passing probes", 2); + me_GE21_fail_chamber_1D->setAxisTitle("Chamber", 1); + me_GE21_fail_chamber_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_chamber_p1_1D->setAxisTitle("Chamber", 1); + me_GEM_pass_chamber_p1_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_chamber_p1_1D->setAxisTitle("Chamber", 1); + me_GEM_fail_chamber_p1_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_chamber_p2_1D->setAxisTitle("Chamber", 1); + me_GEM_pass_chamber_p2_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_chamber_p2_1D->setAxisTitle("Chamber", 1); + me_GEM_fail_chamber_p2_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_chamber_n1_1D->setAxisTitle("Chamber", 1); + me_GEM_pass_chamber_n1_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_chamber_n1_1D->setAxisTitle("Chamber", 1); + me_GEM_fail_chamber_n1_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_chamber_n2_1D->setAxisTitle("Chamber", 1); + me_GEM_pass_chamber_n2_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_chamber_n2_1D->setAxisTitle("Chamber", 1); + me_GEM_fail_chamber_n2_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_pt_1D->setAxisTitle("P_{T}", 1); + me_GEM_pass_pt_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_pt_1D->setAxisTitle("P_{T}", 1); + me_GEM_fail_pt_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_eta_1D->setAxisTitle("#eta", 1); + me_GEM_pass_eta_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_eta_1D->setAxisTitle("#eta", 1); + me_GEM_fail_eta_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_phi_1D->setAxisTitle("#phi", 1); + me_GEM_pass_phi_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_phi_1D->setAxisTitle("#phi", 1); + me_GEM_fail_phi_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_pt_p1_1D->setAxisTitle("P_{T}", 1); + me_GEM_pass_pt_p1_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_pt_p1_1D->setAxisTitle("P_{T}", 1); + me_GEM_fail_pt_p1_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_eta_p1_1D->setAxisTitle("#eta", 1); + me_GEM_pass_eta_p1_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_eta_p1_1D->setAxisTitle("#eta", 1); + me_GEM_fail_eta_p1_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_phi_p1_1D->setAxisTitle("#phi", 1); + me_GEM_pass_phi_p1_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_phi_p1_1D->setAxisTitle("#phi", 1); + me_GEM_fail_phi_p1_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_pt_p2_1D->setAxisTitle("P_{T}", 1); + me_GEM_pass_pt_p2_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_pt_p2_1D->setAxisTitle("P_{T}", 1); + me_GEM_fail_pt_p2_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_eta_p2_1D->setAxisTitle("#eta", 1); + me_GEM_pass_eta_p2_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_eta_p2_1D->setAxisTitle("#eta", 1); + me_GEM_fail_eta_p2_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_phi_p2_1D->setAxisTitle("#phi", 1); + me_GEM_pass_phi_p2_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_phi_p2_1D->setAxisTitle("#phi", 1); + me_GEM_fail_phi_p2_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_pt_n1_1D->setAxisTitle("P_{T}", 1); + me_GEM_pass_pt_n1_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_pt_n1_1D->setAxisTitle("P_{T}", 1); + me_GEM_fail_pt_n1_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_eta_n1_1D->setAxisTitle("#eta", 1); + me_GEM_pass_eta_n1_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_eta_n1_1D->setAxisTitle("#eta", 1); + me_GEM_fail_eta_n1_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_phi_n1_1D->setAxisTitle("#phi", 1); + me_GEM_pass_phi_n1_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_phi_n1_1D->setAxisTitle("#phi", 1); + me_GEM_fail_phi_n1_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_pt_n2_1D->setAxisTitle("P_{T}", 1); + me_GEM_pass_pt_n2_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_pt_n2_1D->setAxisTitle("P_{T}", 1); + me_GEM_fail_pt_n2_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_eta_n2_1D->setAxisTitle("#eta", 1); + me_GEM_pass_eta_n2_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_eta_n2_1D->setAxisTitle("#eta", 1); + me_GEM_fail_eta_n2_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_pass_phi_n2_1D->setAxisTitle("#phi", 1); + me_GEM_pass_phi_n2_1D->setAxisTitle("Number of passing probes", 2); + me_GEM_fail_phi_n2_1D->setAxisTitle("#phi", 1); + me_GEM_fail_phi_n2_1D->setAxisTitle("Number of failing probes", 2); + + me_GE11_fail_Ch_region->setBinLabel(1, "GE-11", 1); + me_GE11_fail_Ch_region->setBinLabel(2, "GE+11", 1); + for (int i = 1; i < 37; ++i) { + me_GE11_fail_Ch_region->setBinLabel(i, std::to_string(i), 2); + } + me_GE11_fail_Ch_region->setAxisTitle("Chamber", 2); + me_GE11_fail_Ch_region->setAxisTitle("Number of failing probes", 3); + + me_GE11_pass_Ch_region->setBinLabel(1, "GE-11", 1); + me_GE11_pass_Ch_region->setBinLabel(2, "GE+11", 1); + for (int i = 1; i < 37; ++i) { + me_GE11_pass_Ch_region->setBinLabel(i, std::to_string(i), 2); + } + me_GE11_pass_Ch_region->setAxisTitle("Chamber", 2); + me_GE11_pass_Ch_region->setAxisTitle("Number of passing probes", 3); + + me_GE21_fail_Ch_region->setBinLabel(1, "GE-21", 1); + me_GE21_fail_Ch_region->setBinLabel(2, "GE+21", 1); + for (int i = 1; i < 19; ++i) { + me_GE21_fail_Ch_region->setBinLabel(i, std::to_string(i), 2); + } + me_GE21_fail_Ch_region->setAxisTitle("Chamber", 2); + me_GE21_fail_Ch_region->setAxisTitle("Number of failing probes", 3); + + me_GE21_pass_Ch_region->setBinLabel(1, "GE-21", 1); + me_GE21_pass_Ch_region->setBinLabel(2, "GE+21", 1); + for (int i = 1; i < 19; ++i) { + me_GE21_pass_Ch_region->setBinLabel(i, std::to_string(i), 2); + } + me_GE21_pass_Ch_region->setAxisTitle("Chamber", 2); + me_GE21_pass_Ch_region->setAxisTitle("Number of passing probes", 3); + + me_GEM_fail_Ch_region_GE1->setBinLabel(1, "GE-1/1_L2", 1); + me_GEM_fail_Ch_region_GE1->setBinLabel(2, "GE-1/1_L1", 1); + me_GEM_fail_Ch_region_GE1->setBinLabel(3, "GE+1/1_L1", 1); + me_GEM_fail_Ch_region_GE1->setBinLabel(4, "GE+1/1_L2", 1); + for (int i = 1; i < 37; ++i) { + me_GEM_fail_Ch_region_GE1->setBinLabel(i, std::to_string(i), 2); + } + me_GEM_fail_Ch_region_GE1->setAxisTitle("Chamber", 2); + me_GEM_fail_Ch_region_GE1->setAxisTitle("Number of passing probes", 3); + + me_GEM_pass_Ch_region_GE1->setBinLabel(1, "GE-1/1_L2", 1); + me_GEM_pass_Ch_region_GE1->setBinLabel(2, "GE-1/1_L1", 1); + me_GEM_pass_Ch_region_GE1->setBinLabel(3, "GE+1/1_L1", 1); + me_GEM_pass_Ch_region_GE1->setBinLabel(4, "GE+1/1_L2", 1); + for (int i = 1; i < 37; ++i) { + me_GEM_pass_Ch_region_GE1->setBinLabel(i, std::to_string(i), 2); + } + me_GEM_pass_Ch_region_GE1->setAxisTitle("Chamber", 2); + me_GEM_pass_Ch_region_GE1->setAxisTitle("Number of passing probes", 3); + + me_GEM_fail_Ch_region_GE1_NoL->setBinLabel(1, "GE-1", 1); + me_GEM_fail_Ch_region_GE1_NoL->setBinLabel(2, "GE+1", 1); + for (int i = 1; i < 37; ++i) { + me_GEM_fail_Ch_region_GE1_NoL->setBinLabel(i, std::to_string(i), 2); + } + me_GEM_fail_Ch_region_GE1_NoL->setAxisTitle("Chamber", 2); + me_GEM_fail_Ch_region_GE1_NoL->setAxisTitle("Number of passing probes", 3); + + me_GEM_pass_Ch_region_GE1_NoL->setBinLabel(1, "GE-1", 1); + me_GEM_pass_Ch_region_GE1_NoL->setBinLabel(2, "GE+1", 1); + for (int i = 1; i < 37; ++i) { + me_GEM_pass_Ch_region_GE1_NoL->setBinLabel(i, std::to_string(i), 2); + } + me_GEM_pass_Ch_region_GE1_NoL->setAxisTitle("Chamber", 2); + me_GEM_pass_Ch_region_GE1_NoL->setAxisTitle("Number of passing probes", 3); + /////////////// + for (int i = 1; i < 37; ++i) { + me_GE11_fail_Ch_ieta->setBinLabel(i, std::to_string(i), 2); + } + for (int i = 1; i < 9; ++i) { + me_GE11_fail_Ch_ieta->setBinLabel(i, std::to_string(i), 1); + } + me_GE11_fail_Ch_ieta->setAxisTitle("#ieta", 1); + me_GE11_fail_Ch_ieta->setAxisTitle("Chamber", 2); + me_GE11_fail_Ch_ieta->setAxisTitle("Number of failing probes", 3); + + for (int i = 1; i < 37; ++i) { + me_GE11_pass_Ch_ieta->setBinLabel(i, std::to_string(i), 2); + } + for (int i = 1; i < 9; ++i) { + me_GE11_pass_Ch_ieta->setBinLabel(i, std::to_string(i), 1); + } + me_GE11_pass_Ch_ieta->setAxisTitle("#ieta", 1); + me_GE11_pass_Ch_ieta->setAxisTitle("Chamber", 2); + me_GE11_pass_Ch_ieta->setAxisTitle("Number of passing probes", 3); + + for (int i = 1; i < 37; ++i) { + me_GE11_fail_Ch_phi->setBinLabel(i, std::to_string(i), 2); + } + me_GE11_fail_Ch_phi->setAxisTitle("#phi", 1); + me_GE11_fail_Ch_phi->setAxisTitle("Chamber", 2); + me_GE11_fail_Ch_phi->setAxisTitle("Number of failing probes", 3); + + for (int i = 1; i < 37; ++i) { + me_GE11_pass_Ch_phi->setBinLabel(i, std::to_string(i), 2); + } + me_GE11_pass_Ch_phi->setAxisTitle("#phi", 1); + me_GE11_pass_Ch_phi->setAxisTitle("Chamber", 2); + me_GE11_pass_Ch_phi->setAxisTitle("Number of passing probes", 3); + + for (int i = 1; i < 19; ++i) { + me_GE21_fail_Ch_ieta->setBinLabel(i, std::to_string(i), 2); + } + for (int i = 1; i < 17; ++i) { + me_GE21_fail_Ch_ieta->setBinLabel(i, std::to_string(i), 1); + } + me_GE21_fail_Ch_ieta->setAxisTitle("#ieta", 1); + me_GE21_fail_Ch_ieta->setAxisTitle("Chamber", 2); + me_GE21_fail_Ch_ieta->setAxisTitle("Number of failing probes", 3); + + for (int i = 1; i < 19; ++i) { + me_GE21_pass_Ch_ieta->setBinLabel(i, std::to_string(i), 2); + } + for (int i = 1; i < 17; ++i) { + me_GE21_pass_Ch_ieta->setBinLabel(i, std::to_string(i), 1); + } + me_GE21_pass_Ch_ieta->setAxisTitle("#ieta", 1); + me_GE21_pass_Ch_ieta->setAxisTitle("Chamber", 2); + me_GE21_pass_Ch_ieta->setAxisTitle("Number of passing probes", 3); + ///////////////////// + for (int i = 1; i < 19; ++i) { + me_GE21_fail_Ch_phi->setBinLabel(i, std::to_string(i), 2); + } + me_GE21_fail_Ch_phi->setAxisTitle("#phi", 1); + me_GE21_fail_Ch_phi->setAxisTitle("Chamber", 2); + me_GE21_fail_Ch_phi->setAxisTitle("Number of failing probes", 3); + + for (int i = 1; i < 19; ++i) { + me_GE21_pass_Ch_phi->setBinLabel(i, std::to_string(i), 2); + } + me_GE21_pass_Ch_phi->setAxisTitle("#phi", 1); + me_GE21_pass_Ch_phi->setAxisTitle("Chamber", 2); + me_GE21_pass_Ch_phi->setAxisTitle("Number of passing probes", 3); + + for (int i = 1; i < 19; ++i) { + me_ME0_pass_chamber_1D->setBinLabel(i, std::to_string(i), 1); + } + me_ME0_pass_chamber_1D->setAxisTitle("Chamber", 1); + me_ME0_pass_chamber_1D->setAxisTitle("Number of passing probes", 2); + for (int i = 1; i < 19; ++i) { + me_ME0_fail_chamber_1D->setBinLabel(i, std::to_string(i), 1); + } + me_ME0_fail_chamber_1D->setAxisTitle("Chamber", 1); + me_ME0_fail_chamber_1D->setAxisTitle("Number of failing probes", 2); + + me_GEM_fail_Ch_region_layer_phase2->setBinLabel(1, "GE-2/1_L2", 1); + me_GEM_fail_Ch_region_layer_phase2->setBinLabel(2, "GE-2/1_L1", 1); + me_GEM_fail_Ch_region_layer_phase2->setBinLabel(3, "GE-1/1_L2", 1); + me_GEM_fail_Ch_region_layer_phase2->setBinLabel(4, "GE-1/1_L1", 1); + me_GEM_fail_Ch_region_layer_phase2->setBinLabel(5, "GE0-", 1); + me_GEM_fail_Ch_region_layer_phase2->setBinLabel(6, "GE0+", 1); + me_GEM_fail_Ch_region_layer_phase2->setBinLabel(7, "GE+1/1_L1", 1); + me_GEM_fail_Ch_region_layer_phase2->setBinLabel(8, "GE+1/1_L2", 1); + me_GEM_fail_Ch_region_layer_phase2->setBinLabel(9, "GE+2/1_L1", 1); + me_GEM_fail_Ch_region_layer_phase2->setBinLabel(10, "GE+2/1_L2", 1); + for (int i = 1; i < 37; ++i) { + me_GEM_fail_Ch_region_layer_phase2->setBinLabel(i, std::to_string(i), 2); + } + me_GEM_fail_Ch_region_layer_phase2->setAxisTitle("Chamber", 2); + me_GEM_fail_Ch_region_layer_phase2->setAxisTitle("Number of passing probes", 3); + + me_GEM_pass_Ch_region_layer_phase2->setBinLabel(1, "GE-2/1_L2", 1); + me_GEM_pass_Ch_region_layer_phase2->setBinLabel(2, "GE-2/1_L1", 1); + me_GEM_pass_Ch_region_layer_phase2->setBinLabel(3, "GE-1/1_L2", 1); + me_GEM_pass_Ch_region_layer_phase2->setBinLabel(4, "GE-1/1_L1", 1); + me_GEM_pass_Ch_region_layer_phase2->setBinLabel(5, "GE0-", 1); + me_GEM_pass_Ch_region_layer_phase2->setBinLabel(6, "GE0+", 1); + me_GEM_pass_Ch_region_layer_phase2->setBinLabel(7, "GE+1/1_L1", 1); + me_GEM_pass_Ch_region_layer_phase2->setBinLabel(8, "GE+1/1_L2", 1); + me_GEM_pass_Ch_region_layer_phase2->setBinLabel(9, "GE+2/1_L1", 1); + me_GEM_pass_Ch_region_layer_phase2->setBinLabel(10, "GE+2/1_L2", 1); + + for (int i = 1; i < 37; ++i) { + me_GEM_pass_Ch_region_layer_phase2->setBinLabel(i, std::to_string(i), 2); + } + me_GEM_pass_Ch_region_layer_phase2->setAxisTitle("Chamber", 2); + me_GEM_pass_Ch_region_layer_phase2->setAxisTitle("Number of passing probes", 3); + + m_histos["GE11_nPassingProbe_Ch_region"] = me_GE11_pass_Ch_region; + m_histos["GE11_nFailingProbe_Ch_region"] = me_GE11_fail_Ch_region; + m_histos["GE21_nPassingProbe_Ch_region"] = me_GE21_pass_Ch_region; + m_histos["GE21_nFailingProbe_Ch_region"] = me_GE21_fail_Ch_region; + m_histos["GEM_nPassingProbe_Ch_region_GE1"] = me_GEM_pass_Ch_region_GE1; + m_histos["GEM_nFailingProbe_Ch_region_GE1"] = me_GEM_fail_Ch_region_GE1; + m_histos["GEM_nPassingProbe_Ch_region_GE1_NoL"] = me_GEM_pass_Ch_region_GE1_NoL; + m_histos["GEM_nFailingProbe_Ch_region_GE1_NoL"] = me_GEM_fail_Ch_region_GE1_NoL; + m_histos["GE11_nPassingProbe_Ch_ieta"] = me_GE11_pass_Ch_ieta; + m_histos["GE11_nFailingProbe_Ch_ieta"] = me_GE11_fail_Ch_ieta; + m_histos["GE11_nPassingProbe_Ch_phi"] = me_GE11_pass_Ch_phi; + m_histos["GE11_nFailingProbe_Ch_phi"] = me_GE11_fail_Ch_phi; + m_histos["GE21_nPassingProbe_Ch_ieta"] = me_GE21_pass_Ch_ieta; + m_histos["GE21_nFailingProbe_Ch_ieta"] = me_GE21_fail_Ch_ieta; + m_histos["GE21_nPassingProbe_Ch_phi"] = me_GE21_pass_Ch_phi; + m_histos["GE21_nFailingProbe_Ch_phi"] = me_GE21_fail_Ch_phi; + m_histos["GE11_nPassingProbe_allCh_1D"] = me_GE11_pass_allCh_1D; + m_histos["GE11_nFailingProbe_allCh_1D"] = me_GE11_fail_allCh_1D; + m_histos["GE21_nPassingProbe_allCh_1D"] = me_GE21_pass_allCh_1D; + m_histos["GE21_nFailingProbe_allCh_1D"] = me_GE21_fail_allCh_1D; + m_histos["GE11_nPassingProbe_chamber_1D"] = me_GE11_pass_chamber_1D; + m_histos["GE11_nFailingProbe_chamber_1D"] = me_GE11_fail_chamber_1D; + m_histos["GE21_nPassingProbe_chamber_1D"] = me_GE21_pass_chamber_1D; + m_histos["GE21_nFailingProbe_chamber_1D"] = me_GE21_fail_chamber_1D; + m_histos["GEM_nPassingProbe_chamber_p1_1D"] = me_GEM_pass_chamber_p1_1D; + m_histos["GEM_nFailingProbe_chamber_p1_1D"] = me_GEM_fail_chamber_p1_1D; + m_histos["GEM_nPassingProbe_chamber_p2_1D"] = me_GEM_pass_chamber_p2_1D; + m_histos["GEM_nFailingProbe_chamber_p2_1D"] = me_GEM_fail_chamber_p2_1D; + m_histos["GEM_nPassingProbe_chamber_n1_1D"] = me_GEM_pass_chamber_n1_1D; + m_histos["GEM_nFailingProbe_chamber_n1_1D"] = me_GEM_fail_chamber_n1_1D; + m_histos["GEM_nPassingProbe_chamber_n2_1D"] = me_GEM_pass_chamber_n2_1D; + m_histos["GEM_nFailingProbe_chamber_n2_1D"] = me_GEM_fail_chamber_n2_1D; + m_histos["GEM_nPassingProbe_pt_1D"] = me_GEM_pass_pt_1D; + m_histos["GEM_nFailingProbe_pt_1D"] = me_GEM_fail_pt_1D; + m_histos["GEM_nPassingProbe_eta_1D"] = me_GEM_pass_eta_1D; + m_histos["GEM_nFailingProbe_eta_1D"] = me_GEM_fail_eta_1D; + m_histos["GEM_nPassingProbe_phi_1D"] = me_GEM_pass_phi_1D; + m_histos["GEM_nFailingProbe_phi_1D"] = me_GEM_fail_phi_1D; + m_histos["GEM_nPassingProbe_pt_p1_1D"] = me_GEM_pass_pt_p1_1D; + m_histos["GEM_nFailingProbe_pt_p1_1D"] = me_GEM_fail_pt_p1_1D; + m_histos["GEM_nPassingProbe_eta_p1_1D"] = me_GEM_pass_eta_p1_1D; + m_histos["GEM_nFailingProbe_eta_p1_1D"] = me_GEM_fail_eta_p1_1D; + m_histos["GEM_nPassingProbe_phi_p1_1D"] = me_GEM_pass_phi_p1_1D; + m_histos["GEM_nFailingProbe_phi_p1_1D"] = me_GEM_fail_phi_p1_1D; + m_histos["GEM_nPassingProbe_pt_p2_1D"] = me_GEM_pass_pt_p2_1D; + m_histos["GEM_nFailingProbe_pt_p2_1D"] = me_GEM_fail_pt_p2_1D; + m_histos["GEM_nPassingProbe_eta_p2_1D"] = me_GEM_pass_eta_p2_1D; + m_histos["GEM_nFailingProbe_eta_p2_1D"] = me_GEM_fail_eta_p2_1D; + m_histos["GEM_nPassingProbe_phi_p2_1D"] = me_GEM_pass_phi_p2_1D; + m_histos["GEM_nFailingProbe_phi_p2_1D"] = me_GEM_fail_phi_p2_1D; + m_histos["GEM_nPassingProbe_pt_n1_1D"] = me_GEM_pass_pt_n1_1D; + m_histos["GEM_nFailingProbe_pt_n1_1D"] = me_GEM_fail_pt_n1_1D; + m_histos["GEM_nPassingProbe_eta_n1_1D"] = me_GEM_pass_eta_n1_1D; + m_histos["GEM_nFailingProbe_eta_n1_1D"] = me_GEM_fail_eta_n1_1D; + m_histos["GEM_nPassingProbe_phi_n1_1D"] = me_GEM_pass_phi_n1_1D; + m_histos["GEM_nFailingProbe_phi_n1_1D"] = me_GEM_fail_phi_n1_1D; + m_histos["GEM_nPassingProbe_pt_n2_1D"] = me_GEM_pass_pt_n2_1D; + m_histos["GEM_nFailingProbe_pt_n2_1D"] = me_GEM_fail_pt_n2_1D; + m_histos["GEM_nPassingProbe_eta_n2_1D"] = me_GEM_pass_eta_n2_1D; + m_histos["GEM_nFailingProbe_eta_n2_1D"] = me_GEM_fail_eta_n2_1D; + m_histos["GEM_nPassingProbe_phi_n2_1D"] = me_GEM_pass_phi_n2_1D; + m_histos["GEM_nFailingProbe_phi_n2_1D"] = me_GEM_fail_phi_n2_1D; + m_histos["ME0_nPassingProbe_chamber_1D"] = me_ME0_pass_chamber_1D; + m_histos["ME0_nFailingProbe_chamber_1D"] = me_ME0_fail_chamber_1D; + m_histos["GEM_nPassingProbe_Ch_region_layer_phase2"] = me_GEM_pass_Ch_region_layer_phase2; + m_histos["GEM_nFailingProbe_Ch_region_layer_phase2"] = me_GEM_fail_Ch_region_layer_phase2; + + std::string baseDir_ = topFolder() + "/detailed/"; + iBooker.setCurrentFolder(baseDir_); + m_histos["GEMseg_dx_ME0"] = iBooker.book1D("GEMseg_dx_ME0", "GEMseg_dx;probe dx [cm];Events", 100, 0., 20.); + m_histos["GEMhit_dx_GE1"] = iBooker.book1D("GEMhit_dx_GE1", "GEMhit_dx;probe dx [cm];Events", 100, 0., 10.); + m_histos["GEMhit_dx_GE2"] = iBooker.book1D("GEMhit_dx_GE2", "GEMhit_dx;probe dx [cm];Events", 100, 0., 10.); + + m_histos["GEMseg_x_ME0"] = iBooker.book1D("GEMhit_x_ME0", "GEMhit_x;probe x [cm];Events", 100, -10., 10.); + m_histos["GEMhit_x_GE1"] = iBooker.book1D("GEMhit_x_GE1", "GEMhit_x;probe x [cm];Events", 100, -10., 10.); + m_histos["GEMhit_x_GE2"] = iBooker.book1D("GEMhit_x_GE2", "GEMhit_x;probe x [cm];Events", 100, -10., 10.); + m_histos["Cham_x_ME0"] = iBooker.book1D("Cham_x_ME0", "Cham_x;probe x [cm];Events", 100, -10., 10.); + m_histos["Cham_x_GE1"] = iBooker.book1D("Cham_x_GE1", "Cham_x;probe x [cm];Events", 100, -10., 10.); + m_histos["Cham_x_GE2"] = iBooker.book1D("Cham_x_GE2", "Cham_x;probe x [cm];Events", 100, -10., 10.); +} + +void GEMTnPEfficiencyTask::analyze(const edm::Event& event, const edm::EventSetup& context) { + BaseTnPEfficiencyTask::analyze(event, context); + + edm::Handle muons; + event.getByToken(m_muToken, muons); + + //GE11 variables + std::vector> probe_coll_GE11_region; + std::vector> probe_coll_GE11_lay; + std::vector> probe_coll_GE11_chamber; + std::vector> probe_coll_GE11_pt; + std::vector> probe_coll_GE11_eta; + std::vector> probe_coll_GE11_ieta; + std::vector> probe_coll_GE11_phi; + std::vector> probe_coll_GE11_sta; + std::vector> probe_coll_GE11_dx; + + //GE21 variables + std::vector> probe_coll_GE21_region; + std::vector> probe_coll_GE21_lay; + std::vector> probe_coll_GE21_chamber; + std::vector> probe_coll_GE21_pt; + std::vector> probe_coll_GE21_eta; + std::vector> probe_coll_GE21_ieta; + std::vector> probe_coll_GE21_phi; + std::vector> probe_coll_GE21_sta; + std::vector> probe_coll_GE21_dx; + + std::vector probe_coll_GEM_staMatch; // ME0 to 0b0001, GE11 to 0b0010, GE21 to 0b0100 + + //ME0 variables + std::vector> probe_coll_ME0_region; + std::vector> probe_coll_ME0_roll; + std::vector> probe_coll_ME0_lay; + std::vector> probe_coll_ME0_chamber; + std::vector> probe_coll_ME0_pt; + std::vector> probe_coll_ME0_eta; + std::vector> probe_coll_ME0_ieta; + std::vector> probe_coll_ME0_phi; + std::vector> probe_coll_ME0_sta; + std::vector> probe_coll_ME0_dx; + + std::vector probe_indices; + if (!m_probeIndices.empty()) + probe_indices = m_probeIndices.back(); + + //Fill probe dx + subdetector coordinates + for (const auto i : probe_indices) { + //GE11 variables + std::vector probe_GE11_region; + std::vector probe_GE11_sta; + std::vector probe_GE11_lay; + std::vector probe_GE11_chamber; + std::vector probe_GE11_pt; + std::vector probe_GE11_eta; + std::vector probe_GE11_ieta; + std::vector probe_GE11_phi; + std::vector probe_GE11_dx; + //GE21 variables + std::vector probe_GE21_region; + std::vector probe_GE21_sta; + std::vector probe_GE21_lay; + std::vector probe_GE21_chamber; + std::vector probe_GE21_pt; + std::vector probe_GE21_eta; + std::vector probe_GE21_ieta; + std::vector probe_GE21_phi; + std::vector probe_GE21_dx; + //std::vector probe_GEM_dx_seg; + uint8_t GEM_stationMatching = 0; + //ME0 variables + std::vector probe_ME0_region; + std::vector probe_ME0_roll; + std::vector probe_ME0_sta; + std::vector probe_ME0_lay; + std::vector probe_ME0_chamber; + std::vector probe_ME0_pt; + std::vector probe_ME0_eta; + std::vector probe_ME0_ieta; + std::vector probe_ME0_phi; + std::vector probe_ME0_dx; + + bool gem_matched = false; // fill detailed plots only for probes matching GEM + + for (const auto& chambMatch : (*muons).at(i).matches()) { + // look in GEMs + bool hit_matched = false; // true if chambermatch has at least one hit (GE11, GE21) or segment (ME0) + if (chambMatch.detector() == MuonSubdetId::GEM) { + if (chambMatch.edgeX < m_borderCut && chambMatch.edgeY < m_borderCut) { + gem_matched = true; //fill detailed plots if at least one GEM probe match + + GEMDetId chId(chambMatch.id.rawId()); + + const int roll = chId.roll(); + const int region = chId.region(); + const int station = chId.station(); + const int layer = chId.layer(); + const int chamber = chId.chamber(); + const int ieta = chId.ieta(); + const float pt = (*muons).at(i).pt(); + const float eta = (*muons).at(i).eta(); + const float phi = (*muons).at(i).phi(); + GEM_stationMatching = GEM_stationMatching | (1 << station); + + if (station == 1 || station == 2) { + reco::MuonGEMHitMatch closest_matchedHit; + double smallestDx = 99999.; + double matched_GEMHit_x = 99999.; + + for (auto& gemHit : chambMatch.gemHitMatches) { + float dx = std::abs(chambMatch.x - gemHit.x); + if (dx < smallestDx) { + smallestDx = dx; + closest_matchedHit = gemHit; + matched_GEMHit_x = gemHit.x; + hit_matched = true; + } + } + + if (station == 1) { + probe_GE11_region.push_back(region); + probe_GE11_sta.push_back(station); + probe_GE11_lay.push_back(layer); + probe_GE11_chamber.push_back(chamber); + probe_GE11_ieta.push_back(ieta); + probe_GE11_pt.push_back(pt); + probe_GE11_eta.push_back(eta); + probe_GE11_phi.push_back(phi); + probe_GE11_dx.push_back(smallestDx); + } + + if (station == 2) { + probe_GE21_region.push_back(region); + probe_GE21_sta.push_back(station); + probe_GE21_lay.push_back(layer); + probe_GE21_chamber.push_back(chamber); + probe_GE21_ieta.push_back(ieta); + probe_GE21_pt.push_back(pt); + probe_GE21_eta.push_back(eta); + probe_GE21_phi.push_back(phi); + probe_GE21_dx.push_back(smallestDx); + } + + if (m_detailedAnalysis && hit_matched) { + if (station == 1) { + m_histos.find("GEMhit_dx_GE1")->second->Fill(smallestDx); + m_histos.find("GEMhit_x_GE1")->second->Fill(matched_GEMHit_x); + m_histos.find("Cham_x_GE1")->second->Fill(chambMatch.x); + } + if (station == 2) { + m_histos.find("GEMhit_dx_GE2")->second->Fill(smallestDx); + m_histos.find("GEMhit_x_GE2")->second->Fill(matched_GEMHit_x); + m_histos.find("Cham_x_GE2")->second->Fill(chambMatch.x); + } + } + } + + if (station == 0) { + reco::MuonSegmentMatch closest_matchedSegment; + double smallestDx_seg = 99999.; + + for (auto& seg : chambMatch.gemMatches) { + float dx_seg = std::abs(chambMatch.x - seg.x); + if (dx_seg < smallestDx_seg) { + smallestDx_seg = dx_seg; + closest_matchedSegment = seg; + hit_matched = true; + } + } + + probe_ME0_region.push_back(region); + probe_ME0_roll.push_back(roll); + probe_ME0_sta.push_back(station); + probe_ME0_lay.push_back(layer); + probe_ME0_chamber.push_back(chamber); + probe_ME0_ieta.push_back(ieta); + probe_ME0_pt.push_back(pt); + probe_ME0_eta.push_back(eta); + probe_ME0_phi.push_back(phi); + probe_ME0_dx.push_back(smallestDx_seg); + + if (m_detailedAnalysis && hit_matched) { + m_histos.find("GEMseg_dx_ME0")->second->Fill(smallestDx_seg); + m_histos.find("GEMseg_x_ME0")->second->Fill(closest_matchedSegment.x); + m_histos.find("Cham_x_ME0")->second->Fill(chambMatch.x); + } + } + } + } else + continue; + } //loop over chamber matches + + //Fill detailed plots + if (m_detailedAnalysis && gem_matched) { + m_histos.find("probeEta")->second->Fill((*muons).at(i).eta()); + m_histos.find("probePhi")->second->Fill((*muons).at(i).phi()); + m_histos.find("probeNumberOfMatchedStations")->second->Fill((*muons).at(i).numberOfMatchedStations()); + m_histos.find("probePt")->second->Fill((*muons).at(i).pt()); + //for(int ii=0; isecond->Fill(probe_GEM_dx[ii]); + // m_histos.find("GEMseg_dx")->second->Fill(probe_GEM_dx_seg[ii]); + //} + } + + //Fill GEM variables + probe_coll_GE11_region.push_back(probe_GE11_region); + probe_coll_GE11_sta.push_back(probe_GE11_sta); + probe_coll_GE11_lay.push_back(probe_GE11_lay); + probe_coll_GE11_chamber.push_back(probe_GE11_chamber); + probe_coll_GE11_ieta.push_back(probe_GE11_ieta); + probe_coll_GE11_pt.push_back(probe_GE11_pt); + probe_coll_GE11_eta.push_back(probe_GE11_eta); + probe_coll_GE11_phi.push_back(probe_GE11_phi); + probe_coll_GE11_dx.push_back(probe_GE11_dx); + + probe_coll_GEM_staMatch.push_back(GEM_stationMatching); + + //Fill GE21 variables + probe_coll_GE21_region.push_back(probe_GE21_region); + probe_coll_GE21_sta.push_back(probe_GE21_sta); + probe_coll_GE21_lay.push_back(probe_GE21_lay); + probe_coll_GE21_chamber.push_back(probe_GE21_chamber); + probe_coll_GE21_ieta.push_back(probe_GE21_ieta); + probe_coll_GE21_pt.push_back(probe_GE21_pt); + probe_coll_GE21_eta.push_back(probe_GE21_eta); + probe_coll_GE21_phi.push_back(probe_GE21_phi); + probe_coll_GE21_dx.push_back(probe_GE21_dx); + + //Fill ME0 variables + probe_coll_ME0_region.push_back(probe_ME0_region); + probe_coll_ME0_roll.push_back(probe_ME0_roll); // same as ieta + probe_coll_ME0_sta.push_back(probe_ME0_sta); + probe_coll_ME0_lay.push_back(probe_ME0_lay); + probe_coll_ME0_chamber.push_back(probe_ME0_chamber); + probe_coll_ME0_ieta.push_back(probe_ME0_ieta); + probe_coll_ME0_pt.push_back(probe_ME0_pt); + probe_coll_ME0_eta.push_back(probe_ME0_eta); + probe_coll_ME0_phi.push_back(probe_ME0_phi); + probe_coll_ME0_dx.push_back(probe_ME0_dx); + + } //loop over probe collection + + //Loop over probes + for (unsigned i = 0; i < probe_indices.size(); ++i) { + //uint8_t GEM_matchPatt = probe_coll_GEM_staMatch.at(i); // ME0 to 0b0001, GE11 to 0b0010, GE21 to 0b0100 + + //Loop over ME0 matches + unsigned nME0_matches = probe_coll_ME0_region.at(i).size(); + for (unsigned j = 0; j < nME0_matches; ++j) { + //ME0 variables + int ME0_region = probe_coll_ME0_region.at(i).at(j); + //int ME0_roll = probe_coll_ME0_roll.at(i).at(j); + //int ME0_sta = probe_coll_ME0_sta.at(i).at(j); + //int ME0_lay = probe_coll_ME0_lay.at(i).at(j); + int ME0_chamber = probe_coll_ME0_chamber.at(i).at(j); + //float ME0_pt = probe_coll_ME0_pt.at(i).at(j); + float ME0_dx = probe_coll_ME0_dx.at(i).at(j); + //float ME0_eta = probe_coll_ME0_eta.at(i).at(j); + //float ME0_phi = probe_coll_ME0_phi.at(i).at(j); + + if (ME0_dx < m_dxCut) { + m_histos.find("ME0_nPassingProbe_chamber_1D")->second->Fill(ME0_chamber); + if (ME0_region < 0) + m_histos.find("GEM_nPassingProbe_Ch_region_layer_phase2")->second->Fill(4, ME0_chamber); + else if (ME0_region > 0) + m_histos.find("GEM_nPassingProbe_Ch_region_layer_phase2")->second->Fill(5, ME0_chamber); + } else { + m_histos.find("ME0_nFailingProbe_chamber_1D")->second->Fill(ME0_chamber); + if (ME0_region < 0) + m_histos.find("GEM_nFailingProbe_Ch_region_layer_phase2")->second->Fill(4, ME0_chamber); + else if (ME0_region > 0) + m_histos.find("GEM_nFailingProbe_Ch_region_layer_phase2")->second->Fill(5, ME0_chamber); + } + } + // + + //Loop over GE11 matches + unsigned nGE11_matches = probe_coll_GE11_region.at(i).size(); + for (unsigned j = 0; j < nGE11_matches; ++j) { + //GEM variables + int GEM_region = probe_coll_GE11_region.at(i).at(j); + int GEM_sta = probe_coll_GE11_sta.at(i).at(j); + int GEM_lay = probe_coll_GE11_lay.at(i).at(j); + int GEM_chamber = probe_coll_GE11_chamber.at(i).at(j); + int GEM_ieta = probe_coll_GE11_ieta.at(i).at(j); + float GEM_pt = probe_coll_GE11_pt.at(i).at(j); + float GEM_dx = probe_coll_GE11_dx.at(i).at(j); + float GEM_eta = probe_coll_GE11_eta.at(i).at(j); + float GEM_phi = probe_coll_GE11_phi.at(i).at(j); + //Fill GEM plots + if (GEM_dx < m_dxCut) { + m_histos.find("GE11_nPassingProbe_Ch_region")->second->Fill(GEM_region, GEM_chamber); + m_histos.find("GE11_nPassingProbe_Ch_ieta")->second->Fill(GEM_ieta, GEM_chamber); + m_histos.find("GE11_nPassingProbe_Ch_phi")->second->Fill(GEM_phi, GEM_chamber); + m_histos.find("GE11_nPassingProbe_allCh_1D")->second->Fill(GEM_region); + m_histos.find("GE11_nPassingProbe_chamber_1D")->second->Fill(GEM_chamber); + if (GEM_region < 0) { + if (GEM_lay == 2) + m_histos.find("GEM_nPassingProbe_Ch_region_layer_phase2")->second->Fill(2, GEM_chamber); + else if (GEM_lay == 1) + m_histos.find("GEM_nPassingProbe_Ch_region_layer_phase2")->second->Fill(3, GEM_chamber); + } + if (GEM_region > 0) { + if (GEM_lay == 1) + m_histos.find("GEM_nPassingProbe_Ch_region_layer_phase2")->second->Fill(6, GEM_chamber); + else if (GEM_lay == 2) + m_histos.find("GEM_nPassingProbe_Ch_region_layer_phase2")->second->Fill(7, GEM_chamber); + } + if (GEM_region == -1) { + m_histos.find("GEM_nPassingProbe_Ch_region_GE1_NoL")->second->Fill(0, GEM_chamber); + } else if (GEM_region == 1) { + m_histos.find("GEM_nPassingProbe_Ch_region_GE1_NoL")->second->Fill(1, GEM_chamber); + } + + if (GEM_region == 1 && GEM_lay == 1) { + m_histos.find("GEM_nPassingProbe_chamber_p1_1D")->second->Fill(GEM_chamber); + m_histos.find("GEM_nPassingProbe_Ch_region_GE1")->second->Fill(2, GEM_chamber); + m_histos.find("GEM_nPassingProbe_pt_p1_1D")->second->Fill(GEM_pt); + m_histos.find("GEM_nPassingProbe_eta_p1_1D")->second->Fill(abs(GEM_eta)); + m_histos.find("GEM_nPassingProbe_phi_p1_1D")->second->Fill(GEM_phi); + } else if (GEM_region == 1 && GEM_lay == 2) { + m_histos.find("GEM_nPassingProbe_chamber_p2_1D")->second->Fill(GEM_chamber); + m_histos.find("GEM_nPassingProbe_Ch_region_GE1")->second->Fill(3, GEM_chamber); + m_histos.find("GEM_nPassingProbe_pt_p2_1D")->second->Fill(GEM_pt); + m_histos.find("GEM_nPassingProbe_eta_p2_1D")->second->Fill(abs(GEM_eta)); + m_histos.find("GEM_nPassingProbe_phi_p2_1D")->second->Fill(GEM_phi); + } else if (GEM_region == -1 && GEM_lay == 1) { + m_histos.find("GEM_nPassingProbe_chamber_n1_1D")->second->Fill(GEM_chamber); + m_histos.find("GEM_nPassingProbe_Ch_region_GE1")->second->Fill(1, GEM_chamber); + m_histos.find("GEM_nPassingProbe_pt_n1_1D")->second->Fill(GEM_pt); + m_histos.find("GEM_nPassingProbe_eta_n1_1D")->second->Fill(abs(GEM_eta)); + m_histos.find("GEM_nPassingProbe_phi_n1_1D")->second->Fill(GEM_phi); + } else if (GEM_region == -1 && GEM_lay == 2) { + m_histos.find("GEM_nPassingProbe_chamber_n2_1D")->second->Fill(GEM_chamber); + m_histos.find("GEM_nPassingProbe_Ch_region_GE1")->second->Fill(0, GEM_chamber); + m_histos.find("GEM_nPassingProbe_pt_n2_1D")->second->Fill(GEM_pt); + m_histos.find("GEM_nPassingProbe_eta_n2_1D")->second->Fill(abs(GEM_eta)); + m_histos.find("GEM_nPassingProbe_phi_n2_1D")->second->Fill(GEM_phi); + } + m_histos.find("GEM_nPassingProbe_pt_1D")->second->Fill(GEM_pt); + m_histos.find("GEM_nPassingProbe_eta_1D")->second->Fill(abs(GEM_eta)); + m_histos.find("GEM_nPassingProbe_phi_1D")->second->Fill(GEM_phi); + } else { + m_histos.find("GE11_nFailingProbe_Ch_region")->second->Fill(GEM_region, GEM_chamber); + m_histos.find("GE11_nFailingProbe_Ch_ieta")->second->Fill(GEM_ieta, GEM_chamber); + m_histos.find("GE11_nFailingProbe_Ch_phi")->second->Fill(GEM_phi, GEM_chamber); + m_histos.find("GE11_nFailingProbe_allCh_1D")->second->Fill(GEM_region); + m_histos.find("GE11_nFailingProbe_chamber_1D")->second->Fill(GEM_chamber); + if (GEM_region < 0) { + if (GEM_sta == 2 and GEM_lay == 2) + m_histos.find("GEM_nFailingProbe_Ch_region_layer_phase2")->second->Fill(0, GEM_chamber); + else if (GEM_sta == 2 and GEM_lay == 1) + m_histos.find("GEM_nFailingProbe_Ch_region_layer_phase2")->second->Fill(1, GEM_chamber); + else if (GEM_sta == 1 and GEM_lay == 2) + m_histos.find("GEM_nFailingProbe_Ch_region_layer_phase2")->second->Fill(2, GEM_chamber); + else if (GEM_sta == 1 and GEM_lay == 1) + m_histos.find("GEM_nFailingProbe_Ch_region_layer_phase2")->second->Fill(3, GEM_chamber); + } + if (GEM_region > 0) { + if (GEM_sta == 1 and GEM_lay == 1) + m_histos.find("GEM_nFailingProbe_Ch_region_layer_phase2")->second->Fill(6, GEM_chamber); + else if (GEM_sta == 1 and GEM_lay == 2) + m_histos.find("GEM_nFailingProbe_Ch_region_layer_phase2")->second->Fill(7, GEM_chamber); + else if (GEM_sta == 2 and GEM_lay == 1) + m_histos.find("GEM_nFailingProbe_Ch_region_layer_phase2")->second->Fill(8, GEM_chamber); + else if (GEM_sta == 2 and GEM_lay == 2) + m_histos.find("GEM_nFailingProbe_Ch_region_layer_phase2")->second->Fill(9, GEM_chamber); + } + if (GEM_region == -1) { + m_histos.find("GEM_nFailingProbe_Ch_region_GE1_NoL")->second->Fill(0, GEM_chamber); + } else if (GEM_region == 1) { + m_histos.find("GEM_nFailingProbe_Ch_region_GE1_NoL")->second->Fill(1, GEM_chamber); + } + // + if (GEM_region == 1 && GEM_lay == 1) { + m_histos.find("GEM_nFailingProbe_chamber_p1_1D")->second->Fill(GEM_chamber); + m_histos.find("GEM_nFailingProbe_Ch_region_GE1")->second->Fill(2, GEM_chamber); + m_histos.find("GEM_nFailingProbe_pt_p1_1D")->second->Fill(GEM_pt); + m_histos.find("GEM_nFailingProbe_eta_p1_1D")->second->Fill(abs(GEM_eta)); + m_histos.find("GEM_nFailingProbe_phi_p1_1D")->second->Fill(GEM_phi); + } else if (GEM_region == 1 && GEM_lay == 2) { + m_histos.find("GEM_nFailingProbe_chamber_p2_1D")->second->Fill(GEM_chamber); + m_histos.find("GEM_nFailingProbe_Ch_region_GE1")->second->Fill(3, GEM_chamber); + m_histos.find("GEM_nFailingProbe_pt_p2_1D")->second->Fill(GEM_pt); + m_histos.find("GEM_nFailingProbe_eta_p2_1D")->second->Fill(abs(GEM_eta)); + m_histos.find("GEM_nFailingProbe_phi_p2_1D")->second->Fill(GEM_phi); + } else if (GEM_region == -1 && GEM_lay == 1) { + m_histos.find("GEM_nFailingProbe_chamber_n1_1D")->second->Fill(GEM_chamber); + m_histos.find("GEM_nFailingProbe_Ch_region_GE1")->second->Fill(1, GEM_chamber); + m_histos.find("GEM_nFailingProbe_pt_n1_1D")->second->Fill(GEM_pt); + m_histos.find("GEM_nFailingProbe_eta_n1_1D")->second->Fill(abs(GEM_eta)); + m_histos.find("GEM_nFailingProbe_phi_n1_1D")->second->Fill(GEM_phi); + } else if (GEM_region == -1 && GEM_lay == 2) { + m_histos.find("GEM_nFailingProbe_chamber_n2_1D")->second->Fill(GEM_chamber); + m_histos.find("GEM_nFailingProbe_Ch_region_GE1")->second->Fill(0, GEM_chamber); + m_histos.find("GEM_nFailingProbe_pt_n2_1D")->second->Fill(GEM_pt); + m_histos.find("GEM_nFailingProbe_eta_n2_1D")->second->Fill(abs(GEM_eta)); + m_histos.find("GEM_nFailingProbe_phi_n2_1D")->second->Fill(GEM_phi); + } + m_histos.find("GEM_nFailingProbe_pt_1D")->second->Fill(GEM_pt); + m_histos.find("GEM_nFailingProbe_eta_1D")->second->Fill(abs(GEM_eta)); + m_histos.find("GEM_nFailingProbe_phi_1D")->second->Fill(GEM_phi); + } + } + + //Loop over GE21 matches + unsigned nGE21_matches = probe_coll_GE21_region.at(i).size(); + for (unsigned j = 0; j < nGE21_matches; ++j) { + //GEM variables + int GEM_region = probe_coll_GE21_region.at(i).at(j); + int GEM_lay = probe_coll_GE21_lay.at(i).at(j); + int GEM_chamber = probe_coll_GE21_chamber.at(i).at(j); + float GEM_ieta = probe_coll_GE21_ieta.at(i).at(j); + float GEM_dx = probe_coll_GE21_dx.at(i).at(j); + float GEM_phi = probe_coll_GE21_phi.at(i).at(j); + + //Fill GEM plots + if (GEM_dx < m_dxCut) { + m_histos.find("GE21_nPassingProbe_Ch_region")->second->Fill(GEM_region, GEM_chamber); + m_histos.find("GE21_nPassingProbe_Ch_ieta")->second->Fill(GEM_ieta, GEM_chamber); + m_histos.find("GE21_nPassingProbe_Ch_phi")->second->Fill(GEM_phi, GEM_chamber); + m_histos.find("GE21_nPassingProbe_allCh_1D")->second->Fill(GEM_region); + m_histos.find("GE21_nPassingProbe_chamber_1D")->second->Fill(GEM_chamber); + if (GEM_region < 0) { + if (GEM_lay == 2) + m_histos.find("GEM_nPassingProbe_Ch_region_layer_phase2")->second->Fill(0, GEM_chamber); + else if (GEM_lay == 1) + m_histos.find("GEM_nPassingProbe_Ch_region_layer_phase2")->second->Fill(1, GEM_chamber); + } + if (GEM_region > 0) { + if (GEM_lay == 1) + m_histos.find("GEM_nPassingProbe_Ch_region_layer_phase2")->second->Fill(8, GEM_chamber); + else if (GEM_lay == 2) + m_histos.find("GEM_nPassingProbe_Ch_region_layer_phase2")->second->Fill(9, GEM_chamber); + } + } else { + m_histos.find("GE21_nFailingProbe_Ch_region")->second->Fill(GEM_region, GEM_chamber); + m_histos.find("GE21_nFailingProbe_Ch_ieta")->second->Fill(GEM_ieta, GEM_chamber); + m_histos.find("GE21_nFailingProbe_Ch_phi")->second->Fill(GEM_phi, GEM_chamber); + m_histos.find("GE21_nFailingProbe_allCh_1D")->second->Fill(GEM_region); + m_histos.find("GE21_nFailingProbe_chamber_1D")->second->Fill(GEM_chamber); + if (GEM_region < 0) { + if (GEM_lay == 2) + m_histos.find("GEM_nFailingProbe_Ch_region_layer_phase2")->second->Fill(0, GEM_chamber); + else if (GEM_lay == 1) + m_histos.find("GEM_nFailingProbe_Ch_region_layer_phase2")->second->Fill(1, GEM_chamber); + } + if (GEM_region > 0) { + if (GEM_lay == 1) + m_histos.find("GEM_nFailingProbe_Ch_region_layer_phase2")->second->Fill(8, GEM_chamber); + else if (GEM_lay == 2) + m_histos.find("GEM_nFailingProbe_Ch_region_layer_phase2")->second->Fill(9, GEM_chamber); + } + } + } + } +} + +std::string GEMTnPEfficiencyTask::topFolder() const { return "GEM/Segment_TnP/"; }; + +DEFINE_FWK_MODULE(GEMTnPEfficiencyTask); \ No newline at end of file diff --git a/DQMOffline/MuonDPG/python/gemTnPEfficiencyClient_cfi.py b/DQMOffline/MuonDPG/python/gemTnPEfficiencyClient_cfi.py new file mode 100644 index 0000000000000..830db19dc7763 --- /dev/null +++ b/DQMOffline/MuonDPG/python/gemTnPEfficiencyClient_cfi.py @@ -0,0 +1,40 @@ +import FWCore.ParameterSet.Config as cms +from DQMServices.Core.DQMEDHarvester import DQMEDHarvester + +gemTnPEfficiencyClient = DQMEDHarvester("TnPEfficiencyClient", + #Histogram names listed as "passProbeHistoName:failProbeHistoName" + subsystem = cms.untracked.string("GEM"), + histoNames = cms.untracked.vstring("GE11_nPassingProbe_Ch_region:GE11_nFailingProbe_Ch_region", + "GE21_nPassingProbe_Ch_region:GE21_nFailingProbe_Ch_region", + "GEM_nPassingProbe_Ch_region_GE1:GEM_nFailingProbe_Ch_region_GE1", + "GEM_nPassingProbe_Ch_region_GE1_NoL:GEM_nFailingProbe_Ch_region_GE1_NoL", + "GE11_nPassingProbe_Ch_ieta:GE11_nFailingProbe_Ch_ieta", + "GE11_nPassingProbe_Ch_phi:GE11_nFailingProbe_Ch_phi", + "GE11_nPassingProbe_allCh_1D:GE11_nFailingProbe_allCh_1D", + "GE11_nPassingProbe_chamber_1D:GE11_nFailingProbe_chamber_1D", + "GE21_nPassingProbe_Ch_ieta:GE21_nFailingProbe_Ch_ieta", + "GE21_nPassingProbe_Ch_phi:GE21_nFailingProbe_Ch_phi", + "GE21_nPassingProbe_allCh_1D:GE21_nFailingProbe_allCh_1D", + "GE21_nPassingProbe_chamber_1D:GE21_nFailingProbe_chamber_1D", + "GEM_nPassingProbe_chamber_p1_1D:GEM_nFailingProbe_chamber_p1_1D", + "GEM_nPassingProbe_chamber_p2_1D:GEM_nFailingProbe_chamber_p2_1D", + "GEM_nPassingProbe_chamber_n1_1D:GEM_nFailingProbe_chamber_n1_1D", + "GEM_nPassingProbe_chamber_n2_1D:GEM_nFailingProbe_chamber_n2_1D", + "GEM_nPassingProbe_pt_1D:GEM_nFailingProbe_pt_1D", + "GEM_nPassingProbe_eta_1D:GEM_nFailingProbe_eta_1D", + "GEM_nPassingProbe_phi_1D:GEM_nFailingProbe_phi_1D", + "GEM_nPassingProbe_pt_p1_1D:GEM_nFailingProbe_pt_p1_1D", + "GEM_nPassingProbe_eta_p1_1D:GEM_nFailingProbe_eta_p1_1D", + "GEM_nPassingProbe_phi_p1_1D:GEM_nFailingProbe_phi_p1_1D", + "GEM_nPassingProbe_pt_p2_1D:GEM_nFailingProbe_pt_p2_1D", + "GEM_nPassingProbe_eta_p2_1D:GEM_nFailingProbe_eta_p2_1D", + "GEM_nPassingProbe_phi_p2_1D:GEM_nFailingProbe_phi_p2_1D", + "GEM_nPassingProbe_pt_n1_1D:GEM_nFailingProbe_pt_n1_1D", + "GEM_nPassingProbe_eta_n1_1D:GEM_nFailingProbe_eta_n1_1D", + "GEM_nPassingProbe_phi_n1_1D:GEM_nFailingProbe_phi_n1_1D", + "GEM_nPassingProbe_pt_n2_1D:GEM_nFailingProbe_pt_n2_1D", + "GEM_nPassingProbe_eta_n2_1D:GEM_nFailingProbe_eta_n2_1D", + "GEM_nPassingProbe_phi_n2_1D:GEM_nFailingProbe_phi_n2_1D", + "ME0_nPassingProbe_chamber_1D:ME0_nFailingProbe_chamber_1D", + "GEM_nPassingProbe_Ch_region_layer_phase2:GEM_nFailingProbe_Ch_region_layer_phase2"), + diagnosticPrescale = cms.untracked.int32(1)) \ No newline at end of file diff --git a/TrackingTools/TrackAssociator/src/TrackDetectorAssociator.cc b/TrackingTools/TrackAssociator/src/TrackDetectorAssociator.cc new file mode 100644 index 0000000000000..b31118fbdba1a --- /dev/null +++ b/TrackingTools/TrackAssociator/src/TrackDetectorAssociator.cc @@ -0,0 +1,1143 @@ +// -*- C++ -*- +// +// Package: TrackAssociator +// Class: TrackDetectorAssociator +// +/* + + Description: + + Implementation: + +*/ +// +// Original Author: Dmytro Kovalskyi +// Created: Fri Apr 21 10:59:41 PDT 2006 +// +// + +#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h" +#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" +#include "TrackingTools/TrackAssociator/interface/DetIdInfo.h" +#include "TrackingTools/Records/interface/DetIdAssociatorRecord.h" + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Utilities/interface/isFinite.h" +#include "DataFormats/Common/interface/Handle.h" + +#include "DataFormats/DetId/interface/DetId.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackExtra.h" +#include "DataFormats/CaloTowers/interface/CaloTower.h" +#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h" +#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h" +#include "DataFormats/EgammaReco/interface/SuperCluster.h" +#include "DataFormats/DTRecHit/interface/DTRecHitCollection.h" +#include "DataFormats/EcalDetId/interface/EBDetId.h" +#include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h" +#include "DataFormats/DTRecHit/interface/DTRecSegment2D.h" +#include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h" +#include "DataFormats/GEMRecHit/interface/GEMSegmentCollection.h" +#include "DataFormats/GEMRecHit/interface/ME0SegmentCollection.h" +#include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h" +#include "DataFormats/GEMRecHit/interface/GEMRecHitCollection.h" +#include "DataFormats/GEMRecHit/interface/ME0RecHitCollection.h" + +// calorimeter and muon infos +#include "Geometry/CommonDetUnit/interface/GeomDet.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "Geometry/Records/interface/MuonGeometryRecord.h" +#include "Geometry/DTGeometry/interface/DTGeometry.h" +#include "Geometry/CSCGeometry/interface/CSCGeometry.h" +#include "Geometry/GEMGeometry/interface/GEMGeometry.h" +#include "Geometry/GEMGeometry/interface/ME0Geometry.h" +#include "Geometry/GEMGeometry/interface/GEMEtaPartitionSpecs.h" + +#include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h" +#include +#include + +#include "DataFormats/Math/interface/LorentzVector.h" +#include "Math/VectorUtil.h" +#include + +#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/CaloHit/interface/PCaloHit.h" +#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h" + +using namespace reco; + +TrackDetectorAssociator::TrackDetectorAssociator() { + ivProp_ = nullptr; + useDefaultPropagator_ = false; +} + +TrackDetectorAssociator::~TrackDetectorAssociator() = default; + +void TrackDetectorAssociator::setPropagator(const Propagator* ptr) { + ivProp_ = ptr; + cachedTrajectory_.setPropagator(ivProp_); +} + +void TrackDetectorAssociator::useDefaultPropagator() { useDefaultPropagator_ = true; } + +void TrackDetectorAssociator::init(const edm::EventSetup& iSetup, const AssociatorParameters& parameters) { + // access the calorimeter geometry + theCaloGeometry_ = &iSetup.getData(parameters.theCaloGeometryToken); + + // get the tracking Geometry + theTrackingGeometry_ = &iSetup.getData(parameters.theTrackingGeometryToken); + + if (useDefaultPropagator_ && (!defProp_ || theMagneticFieldWatcher_.check(iSetup))) { + // setup propagator + const MagneticField* bField = &iSetup.getData(parameters.bFieldToken); + + auto prop = std::make_unique(bField, anyDirection); + prop->setMaterialMode(false); + prop->applyRadX0Correction(true); + // prop->setDebug(true); // tmp + defProp_ = std::move(prop); + setPropagator(defProp_.get()); + } + + ecalDetIdAssociator_ = &iSetup.getData(parameters.ecalDetIdAssociatorToken); + hcalDetIdAssociator_ = &iSetup.getData(parameters.hcalDetIdAssociatorToken); + hoDetIdAssociator_ = &iSetup.getData(parameters.hoDetIdAssociatorToken); + caloDetIdAssociator_ = &iSetup.getData(parameters.caloDetIdAssociatorToken); + muonDetIdAssociator_ = &iSetup.getData(parameters.muonDetIdAssociatorToken); + preshowerDetIdAssociator_ = &iSetup.getData(parameters.preshowerDetIdAssociatorToken); +} + +TrackDetMatchInfo TrackDetectorAssociator::associate(const edm::Event& iEvent, + const edm::EventSetup& iSetup, + const FreeTrajectoryState& fts, + const AssociatorParameters& parameters) { + return associate(iEvent, iSetup, parameters, &fts); +} + +TrackDetMatchInfo TrackDetectorAssociator::associate(const edm::Event& iEvent, + const edm::EventSetup& iSetup, + const AssociatorParameters& parameters, + const FreeTrajectoryState* innerState, + const FreeTrajectoryState* outerState) { + TrackDetMatchInfo info; + if (!parameters.useEcal && !parameters.useCalo && !parameters.useHcal && !parameters.useHO && !parameters.useMuon && + !parameters.usePreshower) + throw cms::Exception("ConfigurationError") + << "Configuration error! No subdetector was selected for the track association."; + + SteppingHelixStateInfo trackOrigin(*innerState); + info.stateAtIP = *innerState; + cachedTrajectory_.setStateAtIP(trackOrigin); + + init(iSetup, parameters); + // get track trajectory + // ECAL points (EB+EE) + // If the phi angle between a track entrance and exit points is more + // than 2 crystals, it is possible that the track will cross 3 crystals + // and therefore one has to check at least 3 points along the track + // trajectory inside ECAL. In order to have a chance to cross 4 crystalls + // in the barrel, a track should have P_t as low as 3 GeV or smaller + // If it's necessary, number of points along trajectory can be increased + + info.setCaloGeometry(theCaloGeometry_); + + cachedTrajectory_.reset_trajectory(); + // estimate propagation outer boundaries based on + // requested sub-detector information. For now limit + // propagation region only if muon matching is not + // requested. + double HOmaxR = hoDetIdAssociator_->volume().maxR(); + double HOmaxZ = hoDetIdAssociator_->volume().maxZ(); + double minR = ecalDetIdAssociator_->volume().minR(); + double minZ = preshowerDetIdAssociator_->volume().minZ(); + cachedTrajectory_.setMaxHORadius(HOmaxR); + cachedTrajectory_.setMaxHOLength(HOmaxZ * 2.); + cachedTrajectory_.setMinDetectorRadius(minR); + cachedTrajectory_.setMinDetectorLength(minZ * 2.); + + if (parameters.useMuon) { + double maxR = muonDetIdAssociator_->volume().maxR(); + double maxZ = muonDetIdAssociator_->volume().maxZ(); + cachedTrajectory_.setMaxDetectorRadius(maxR); + cachedTrajectory_.setMaxDetectorLength(maxZ * 2.); + } else { + cachedTrajectory_.setMaxDetectorRadius(HOmaxR); + cachedTrajectory_.setMaxDetectorLength(HOmaxZ * 2.); + } + + // If track extras exist and outerState is before HO maximum, then use outerState + if (outerState) { + if (outerState->position().perp() < HOmaxR && std::abs(outerState->position().z()) < HOmaxZ) { + LogTrace("TrackAssociator") << "Using outerState as trackOrigin at Rho=" << outerState->position().perp() + << " Z=" << outerState->position().z() << "\n"; + trackOrigin = SteppingHelixStateInfo(*outerState); + } else if (innerState) { + LogTrace("TrackAssociator") << "Using innerState as trackOrigin at Rho=" << innerState->position().perp() + << " Z=" << innerState->position().z() << "\n"; + trackOrigin = SteppingHelixStateInfo(*innerState); + } + } + + if (trackOrigin.momentum().mag() == 0) + return info; + if (edm::isNotFinite(trackOrigin.momentum().x()) or edm::isNotFinite(trackOrigin.momentum().y()) or + edm::isNotFinite(trackOrigin.momentum().z())) + return info; + if (!cachedTrajectory_.propagateAll(trackOrigin)) + return info; + + // get trajectory in calorimeters + cachedTrajectory_.findEcalTrajectory(ecalDetIdAssociator_->volume()); + cachedTrajectory_.findHcalTrajectory(hcalDetIdAssociator_->volume()); + cachedTrajectory_.findHOTrajectory(hoDetIdAssociator_->volume()); + cachedTrajectory_.findPreshowerTrajectory(preshowerDetIdAssociator_->volume()); + + info.trkGlobPosAtEcal = getPoint(cachedTrajectory_.getStateAtEcal().position()); + info.trkGlobPosAtHcal = getPoint(cachedTrajectory_.getStateAtHcal().position()); + info.trkGlobPosAtHO = getPoint(cachedTrajectory_.getStateAtHO().position()); + + info.trkMomAtEcal = cachedTrajectory_.getStateAtEcal().momentum(); + info.trkMomAtHcal = cachedTrajectory_.getStateAtHcal().momentum(); + info.trkMomAtHO = cachedTrajectory_.getStateAtHO().momentum(); + + if (parameters.useEcal) + fillEcal(iEvent, info, parameters); + if (parameters.useCalo) + fillCaloTowers(iEvent, info, parameters); + if (parameters.useHcal) + fillHcal(iEvent, info, parameters); + if (parameters.useHO) + fillHO(iEvent, info, parameters); + if (parameters.usePreshower) + fillPreshower(iEvent, info, parameters); + if (parameters.useMuon) + fillMuon(iEvent, info, parameters); + if (parameters.truthMatch) + fillCaloTruth(iEvent, info, parameters); + + return info; +} + +void TrackDetectorAssociator::fillEcal(const edm::Event& iEvent, + TrackDetMatchInfo& info, + const AssociatorParameters& parameters) { + const std::vector& trajectoryStates = cachedTrajectory_.getEcalTrajectory(); + + for (std::vector::const_iterator itr = trajectoryStates.begin(); + itr != trajectoryStates.end(); + itr++) + LogTrace("TrackAssociator") << "ECAL trajectory point (rho, z, phi): " << itr->position().perp() << ", " + << itr->position().z() << ", " << itr->position().phi(); + + std::vector coreTrajectory; + for (std::vector::const_iterator itr = trajectoryStates.begin(); + itr != trajectoryStates.end(); + itr++) + coreTrajectory.push_back(itr->position()); + + if (coreTrajectory.empty()) { + LogTrace("TrackAssociator") << "ECAL track trajectory is empty; moving on\n"; + info.isGoodEcal = false; + return; + } + info.isGoodEcal = true; + + // Find ECAL crystals + edm::Handle EBRecHits; + iEvent.getByToken(parameters.EBRecHitsToken, EBRecHits); + if (!EBRecHits.isValid()) + throw cms::Exception("FatalError") << "Unable to find EBRecHitCollection in the event!\n"; + + edm::Handle EERecHits; + iEvent.getByToken(parameters.EERecHitsToken, EERecHits); + if (!EERecHits.isValid()) + throw cms::Exception("FatalError") << "Unable to find EERecHitCollection in event!\n"; + + std::set ecalIdsInRegion; + if (parameters.accountForTrajectoryChangeCalo) { + // get trajectory change with respect to initial state + DetIdAssociator::MapRange mapRange = + getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToEcal), parameters.dREcalPreselection); + ecalIdsInRegion = ecalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], mapRange); + } else + ecalIdsInRegion = ecalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dREcalPreselection); + LogTrace("TrackAssociator") << "ECAL hits in the region: " << ecalIdsInRegion.size(); + if (parameters.dREcalPreselection > parameters.dREcal) + ecalIdsInRegion = ecalDetIdAssociator_->getDetIdsInACone(ecalIdsInRegion, coreTrajectory, parameters.dREcal); + LogTrace("TrackAssociator") << "ECAL hits in the cone: " << ecalIdsInRegion.size(); + info.crossedEcalIds = ecalDetIdAssociator_->getCrossedDetIds(ecalIdsInRegion, coreTrajectory); + const std::vector& crossedEcalIds = info.crossedEcalIds; + LogTrace("TrackAssociator") << "ECAL crossed hits " << crossedEcalIds.size(); + + // add EcalRecHits + for (std::vector::const_iterator itr = crossedEcalIds.begin(); itr != crossedEcalIds.end(); itr++) { + std::vector::const_iterator ebHit = (*EBRecHits).find(*itr); + std::vector::const_iterator eeHit = (*EERecHits).find(*itr); + if (ebHit != (*EBRecHits).end()) + info.crossedEcalRecHits.push_back(&*ebHit); + else if (eeHit != (*EERecHits).end()) + info.crossedEcalRecHits.push_back(&*eeHit); + else + LogTrace("TrackAssociator") << "Crossed EcalRecHit is not found for DetId: " << itr->rawId(); + } + for (std::set::const_iterator itr = ecalIdsInRegion.begin(); itr != ecalIdsInRegion.end(); itr++) { + std::vector::const_iterator ebHit = (*EBRecHits).find(*itr); + std::vector::const_iterator eeHit = (*EERecHits).find(*itr); + if (ebHit != (*EBRecHits).end()) + info.ecalRecHits.push_back(&*ebHit); + else if (eeHit != (*EERecHits).end()) + info.ecalRecHits.push_back(&*eeHit); + else + LogTrace("TrackAssociator") << "EcalRecHit from the cone is not found for DetId: " << itr->rawId(); + } +} + +void TrackDetectorAssociator::fillCaloTowers(const edm::Event& iEvent, + TrackDetMatchInfo& info, + const AssociatorParameters& parameters) { + // use ECAL and HCAL trajectories to match a tower. (HO isn't used for matching). + std::vector trajectory; + const std::vector& ecalTrajectoryStates = cachedTrajectory_.getEcalTrajectory(); + const std::vector& hcalTrajectoryStates = cachedTrajectory_.getHcalTrajectory(); + for (std::vector::const_iterator itr = ecalTrajectoryStates.begin(); + itr != ecalTrajectoryStates.end(); + itr++) + trajectory.push_back(itr->position()); + for (std::vector::const_iterator itr = hcalTrajectoryStates.begin(); + itr != hcalTrajectoryStates.end(); + itr++) + trajectory.push_back(itr->position()); + + if (trajectory.empty()) { + LogTrace("TrackAssociator") << "HCAL trajectory is empty; moving on\n"; + info.isGoodCalo = false; + return; + } + info.isGoodCalo = true; + + // find crossed CaloTowers + edm::Handle caloTowers; + iEvent.getByToken(parameters.caloTowersToken, caloTowers); + if (!caloTowers.isValid()) + throw cms::Exception("FatalError") << "Unable to find CaloTowers in event!\n"; + + std::set caloTowerIdsInRegion; + if (parameters.accountForTrajectoryChangeCalo) { + // get trajectory change with respect to initial state + DetIdAssociator::MapRange mapRange = + getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHcal), parameters.dRHcalPreselection); + caloTowerIdsInRegion = caloDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0], mapRange); + } else + caloTowerIdsInRegion = caloDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0], parameters.dRHcalPreselection); + + LogTrace("TrackAssociator") << "Towers in the region: " << caloTowerIdsInRegion.size(); + + auto caloTowerIdsInAConeBegin = caloTowerIdsInRegion.begin(); + auto caloTowerIdsInAConeEnd = caloTowerIdsInRegion.end(); + std::set caloTowerIdsInAConeTmp; + if (!caloDetIdAssociator_->selectAllInACone(parameters.dRHcal)) { + caloTowerIdsInAConeTmp = + caloDetIdAssociator_->getDetIdsInACone(caloTowerIdsInRegion, trajectory, parameters.dRHcal); + caloTowerIdsInAConeBegin = caloTowerIdsInAConeTmp.begin(); + caloTowerIdsInAConeEnd = caloTowerIdsInAConeTmp.end(); + } + LogTrace("TrackAssociator") << "Towers in the cone: " + << std::distance(caloTowerIdsInAConeBegin, caloTowerIdsInAConeEnd); + + info.crossedTowerIds = caloDetIdAssociator_->getCrossedDetIds(caloTowerIdsInRegion, trajectory); + const std::vector& crossedCaloTowerIds = info.crossedTowerIds; + LogTrace("TrackAssociator") << "Towers crossed: " << crossedCaloTowerIds.size(); + + // add CaloTowers + for (std::vector::const_iterator itr = crossedCaloTowerIds.begin(); itr != crossedCaloTowerIds.end(); itr++) { + CaloTowerCollection::const_iterator tower = (*caloTowers).find(*itr); + if (tower != (*caloTowers).end()) + info.crossedTowers.push_back(&*tower); + else + LogTrace("TrackAssociator") << "Crossed CaloTower is not found for DetId: " << (*itr).rawId(); + } + + for (std::set::const_iterator itr = caloTowerIdsInAConeBegin; itr != caloTowerIdsInAConeEnd; itr++) { + CaloTowerCollection::const_iterator tower = (*caloTowers).find(*itr); + if (tower != (*caloTowers).end()) + info.towers.push_back(&*tower); + else + LogTrace("TrackAssociator") << "CaloTower from the cone is not found for DetId: " << (*itr).rawId(); + } +} + +void TrackDetectorAssociator::fillPreshower(const edm::Event& iEvent, + TrackDetMatchInfo& info, + const AssociatorParameters& parameters) { + std::vector trajectory; + const std::vector& trajectoryStates = cachedTrajectory_.getPreshowerTrajectory(); + for (std::vector::const_iterator itr = trajectoryStates.begin(); + itr != trajectoryStates.end(); + itr++) + trajectory.push_back(itr->position()); + + if (trajectory.empty()) { + LogTrace("TrackAssociator") << "Preshower trajectory is empty; moving on\n"; + return; + } + + std::set idsInRegion = + preshowerDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0], parameters.dRPreshowerPreselection); + + LogTrace("TrackAssociator") << "Number of Preshower Ids in the region: " << idsInRegion.size(); + info.crossedPreshowerIds = preshowerDetIdAssociator_->getCrossedDetIds(idsInRegion, trajectory); + LogTrace("TrackAssociator") << "Number of Preshower Ids in crossed: " << info.crossedPreshowerIds.size(); +} + +void TrackDetectorAssociator::fillHcal(const edm::Event& iEvent, + TrackDetMatchInfo& info, + const AssociatorParameters& parameters) { + const std::vector& trajectoryStates = cachedTrajectory_.getHcalTrajectory(); + + std::vector coreTrajectory; + for (std::vector::const_iterator itr = trajectoryStates.begin(); + itr != trajectoryStates.end(); + itr++) + coreTrajectory.push_back(itr->position()); + + if (coreTrajectory.empty()) { + LogTrace("TrackAssociator") << "HCAL trajectory is empty; moving on\n"; + info.isGoodHcal = false; + return; + } + info.isGoodHcal = true; + + // find crossed Hcals + edm::Handle collection; + iEvent.getByToken(parameters.HBHEcollToken, collection); + if (!collection.isValid()) + throw cms::Exception("FatalError") << "Unable to find HBHERecHits in event!\n"; + + std::set idsInRegion; + if (parameters.accountForTrajectoryChangeCalo) { + // get trajectory change with respect to initial state + DetIdAssociator::MapRange mapRange = + getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHcal), parameters.dRHcalPreselection); + idsInRegion = hcalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], mapRange); + } else + idsInRegion = hcalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dRHcalPreselection); + + LogTrace("TrackAssociator") << "HCAL hits in the region: " << idsInRegion.size() << "\n" + << DetIdInfo::info(idsInRegion, nullptr); + + auto idsInAConeBegin = idsInRegion.begin(); + auto idsInAConeEnd = idsInRegion.end(); + std::set idsInAConeTmp; + if (!hcalDetIdAssociator_->selectAllInACone(parameters.dRHcal)) { + idsInAConeTmp = hcalDetIdAssociator_->getDetIdsInACone(idsInRegion, coreTrajectory, parameters.dRHcal); + idsInAConeBegin = idsInAConeTmp.begin(); + idsInAConeEnd = idsInAConeTmp.end(); + } + LogTrace("TrackAssociator") << "HCAL hits in the cone: " << std::distance(idsInAConeBegin, idsInAConeEnd) << "\n" + << DetIdInfo::info(std::set(idsInAConeBegin, idsInAConeEnd), nullptr); + info.crossedHcalIds = hcalDetIdAssociator_->getCrossedDetIds(idsInRegion, coreTrajectory); + const std::vector& crossedIds = info.crossedHcalIds; + LogTrace("TrackAssociator") << "HCAL hits crossed: " << crossedIds.size() << "\n" + << DetIdInfo::info(crossedIds, nullptr); + + // add Hcal + for (std::vector::const_iterator itr = crossedIds.begin(); itr != crossedIds.end(); itr++) { + HBHERecHitCollection::const_iterator hit = (*collection).find(*itr); + if (hit != (*collection).end()) + info.crossedHcalRecHits.push_back(&*hit); + else + LogTrace("TrackAssociator") << "Crossed HBHERecHit is not found for DetId: " << itr->rawId(); + } + for (std::set::const_iterator itr = idsInAConeBegin; itr != idsInAConeEnd; itr++) { + HBHERecHitCollection::const_iterator hit = (*collection).find(*itr); + if (hit != (*collection).end()) + info.hcalRecHits.push_back(&*hit); + else + LogTrace("TrackAssociator") << "HBHERecHit from the cone is not found for DetId: " << itr->rawId(); + } +} + +void TrackDetectorAssociator::fillHO(const edm::Event& iEvent, + TrackDetMatchInfo& info, + const AssociatorParameters& parameters) { + const std::vector& trajectoryStates = cachedTrajectory_.getHOTrajectory(); + + std::vector coreTrajectory; + for (std::vector::const_iterator itr = trajectoryStates.begin(); + itr != trajectoryStates.end(); + itr++) + coreTrajectory.push_back(itr->position()); + + if (coreTrajectory.empty()) { + LogTrace("TrackAssociator") << "HO trajectory is empty; moving on\n"; + info.isGoodHO = false; + return; + } + info.isGoodHO = true; + + // find crossed HOs + edm::Handle collection; + iEvent.getByToken(parameters.HOcollToken, collection); + if (!collection.isValid()) + throw cms::Exception("FatalError") << "Unable to find HORecHits in event!\n"; + + std::set idsInRegion; + if (parameters.accountForTrajectoryChangeCalo) { + // get trajectory change with respect to initial state + DetIdAssociator::MapRange mapRange = + getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHO), parameters.dRHcalPreselection); + idsInRegion = hoDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], mapRange); + } else + idsInRegion = hoDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dRHcalPreselection); + + LogTrace("TrackAssociator") << "idsInRegion.size(): " << idsInRegion.size(); + + auto idsInAConeBegin = idsInRegion.begin(); + auto idsInAConeEnd = idsInRegion.end(); + std::set idsInAConeTmp; + if (!hoDetIdAssociator_->selectAllInACone(parameters.dRHcal)) { + idsInAConeTmp = hoDetIdAssociator_->getDetIdsInACone(idsInRegion, coreTrajectory, parameters.dRHcal); + idsInAConeBegin = idsInAConeTmp.begin(); + idsInAConeEnd = idsInAConeTmp.end(); + } + LogTrace("TrackAssociator") << "idsInACone.size(): " << std::distance(idsInAConeBegin, idsInAConeEnd); + info.crossedHOIds = hoDetIdAssociator_->getCrossedDetIds(idsInRegion, coreTrajectory); + const std::vector& crossedIds = info.crossedHOIds; + + // add HO + for (std::vector::const_iterator itr = crossedIds.begin(); itr != crossedIds.end(); itr++) { + HORecHitCollection::const_iterator hit = (*collection).find(*itr); + if (hit != (*collection).end()) + info.crossedHORecHits.push_back(&*hit); + else + LogTrace("TrackAssociator") << "Crossed HORecHit is not found for DetId: " << itr->rawId(); + } + + for (std::set::const_iterator itr = idsInAConeBegin; itr != idsInAConeEnd; itr++) { + HORecHitCollection::const_iterator hit = (*collection).find(*itr); + if (hit != (*collection).end()) + info.hoRecHits.push_back(&*hit); + else + LogTrace("TrackAssociator") << "HORecHit from the cone is not found for DetId: " << itr->rawId(); + } +} + +FreeTrajectoryState TrackDetectorAssociator::getFreeTrajectoryState(const MagneticField* bField, + const SimTrack& track, + const SimVertex& vertex) { + GlobalVector vector(track.momentum().x(), track.momentum().y(), track.momentum().z()); + GlobalPoint point(vertex.position().x(), vertex.position().y(), vertex.position().z()); + + int charge = track.type() > 0 ? -1 : 1; // lepton convention + if (abs(track.type()) == 211 || // pion + abs(track.type()) == 321 || // kaon + abs(track.type()) == 2212) + charge = track.type() < 0 ? -1 : 1; + return getFreeTrajectoryState(bField, vector, point, charge); +} + +FreeTrajectoryState TrackDetectorAssociator::getFreeTrajectoryState(const MagneticField* bField, + const GlobalVector& momentum, + const GlobalPoint& vertex, + const int charge) { + GlobalTrajectoryParameters tPars(vertex, momentum, charge, bField); + + ROOT::Math::SMatrixIdentity id; + AlgebraicSymMatrix66 covT(id); + covT *= 1e-6; // initialize to sigma=1e-3 + CartesianTrajectoryError tCov(covT); + + return FreeTrajectoryState(tPars, tCov); +} + +FreeTrajectoryState TrackDetectorAssociator::getFreeTrajectoryState(const MagneticField* bField, + const reco::Track& track) { + GlobalVector vector(track.momentum().x(), track.momentum().y(), track.momentum().z()); + + GlobalPoint point(track.vertex().x(), track.vertex().y(), track.vertex().z()); + + GlobalTrajectoryParameters tPars(point, vector, track.charge(), bField); + + // FIX THIS !!! + // need to convert from perigee to global or helix (curvilinear) frame + // for now just an arbitrary matrix. + ROOT::Math::SMatrixIdentity id; + AlgebraicSymMatrix66 covT(id); + covT *= 1e-6; // initialize to sigma=1e-3 + CartesianTrajectoryError tCov(covT); + + return FreeTrajectoryState(tPars, tCov); +} + +DetIdAssociator::MapRange TrackDetectorAssociator::getMapRange(const std::pair& delta, const float dR) { + DetIdAssociator::MapRange mapRange; + mapRange.dThetaPlus = dR; + mapRange.dThetaMinus = dR; + mapRange.dPhiPlus = dR; + mapRange.dPhiMinus = dR; + if (delta.first > 0) + mapRange.dThetaPlus += delta.first; + else + mapRange.dThetaMinus += std::abs(delta.first); + if (delta.second > 0) + mapRange.dPhiPlus += delta.second; + else + mapRange.dPhiMinus += std::abs(delta.second); + LogTrace("TrackAssociator") << "Selection range: (dThetaPlus, dThetaMinus, dPhiPlus, dPhiMinus, dRPreselection): " + << mapRange.dThetaPlus << ", " << mapRange.dThetaMinus << ", " << mapRange.dPhiPlus + << ", " << mapRange.dPhiMinus << ", " << dR; + return mapRange; +} + +void TrackDetectorAssociator::getTAMuonChamberMatches(std::vector& matches, + const AssociatorParameters& parameters, + std::set occupancy) { + // Strategy: + // Propagate through the whole detector, estimate change in eta and phi + // along the trajectory, add this to dRMuon and find DetIds around this + // direction using the map. Then propagate fast to each surface and apply + // final matching criteria. + + // get the direction first + SteppingHelixStateInfo trajectoryPoint = cachedTrajectory_.getStateAtHcal(); + // If trajectory point at HCAL is not valid, try to use the outer most state of the + // trajectory instead. + if (!trajectoryPoint.isValid()) + trajectoryPoint = cachedTrajectory_.getOuterState(); + if (!trajectoryPoint.isValid()) { + LogTrace("TrackAssociator") + << "trajectory position at HCAL is not valid. Assume the track cannot reach muon detectors and skip it"; + return; + } + + GlobalVector direction = trajectoryPoint.momentum().unit(); + LogTrace("TrackAssociator") << "muon direction: " << direction + << "\n\t and corresponding point: " << trajectoryPoint.position() << "\n"; + + DetIdAssociator::MapRange mapRange = + getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::FullTrajectory), parameters.dRMuonPreselection); + + // and find chamber DetIds + + std::set muonIdsInRegion = muonDetIdAssociator_->getDetIdsCloseToAPoint(trajectoryPoint.position(), mapRange); + LogTrace("TrackAssociator") << "Number of chambers to check: " << muonIdsInRegion.size(); + + if (parameters.preselectMuonTracks) { + std::set muonIdsInRegionOccupied; + std::set_intersection(muonIdsInRegion.begin(), + muonIdsInRegion.end(), + occupancy.begin(), + occupancy.end(), + std::inserter(muonIdsInRegionOccupied, muonIdsInRegionOccupied.begin())); + if (muonIdsInRegionOccupied.empty()) + return; + } + + for (std::set::const_iterator detId = muonIdsInRegion.begin(); detId != muonIdsInRegion.end(); detId++) { + const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet(*detId); + TrajectoryStateOnSurface stateOnSurface = cachedTrajectory_.propagate(&geomDet->surface()); + if (!stateOnSurface.isValid()) { + LogTrace("TrackAssociator") << "Failed to propagate the track; moving on\n\t" + << "Element is not crosssed: " << DetIdInfo::info(*detId, nullptr) << "\n"; + continue; + } + LocalPoint localPoint = geomDet->surface().toLocal(stateOnSurface.freeState()->position()); + LocalError localError = stateOnSurface.localError().positionError(); + float distanceX = 0.f; + float distanceY = 0.f; + if (const CSCChamber* cscChamber = dynamic_cast(geomDet)) { + const CSCChamberSpecs* chamberSpecs = cscChamber->specs(); + if (!chamberSpecs) { + LogTrace("TrackAssociator") << "Failed to get CSCChamberSpecs from CSCChamber; moving on\n"; + continue; + } + const CSCLayerGeometry* layerGeometry = chamberSpecs->oddLayerGeometry(1); + if (!layerGeometry) { + LogTrace("TrackAssociator") << "Failed to get CSCLayerGeometry from CSCChamberSpecs; moving on\n"; + continue; + } + const CSCWireTopology* wireTopology = layerGeometry->wireTopology(); + if (!wireTopology) { + LogTrace("TrackAssociator") << "Failed to get CSCWireTopology from CSCLayerGeometry; moving on\n"; + continue; + } + + float wideWidth = wireTopology->wideWidthOfPlane(); + float narrowWidth = wireTopology->narrowWidthOfPlane(); + float length = wireTopology->lengthOfPlane(); + // If slanted, there is no y offset between local origin and symmetry center of wire plane + float yOfFirstWire = std::abs(wireTopology->wireAngle()) > 1.E-06f ? -0.5 * length : wireTopology->yOfWire(1); + // y offset between local origin and symmetry center of wire plane + float yCOWPOffset = yOfFirstWire + 0.5f * length; + + // tangent of the incline angle from inside the trapezoid + float tangent = (wideWidth - narrowWidth) / (2.f * length); + // y position wrt bottom of trapezoid + float yPrime = localPoint.y() + std::abs(yOfFirstWire); + // half trapezoid width at y' is 0.5 * narrowWidth + x side of triangle with the above tangent and side y' + float halfWidthAtYPrime = 0.5f * narrowWidth + yPrime * tangent; + distanceX = std::abs(localPoint.x()) - halfWidthAtYPrime; + distanceY = std::abs(localPoint.y() - yCOWPOffset) - 0.5f * length; + } else if (dynamic_cast(geomDet) || dynamic_cast(geomDet)) { + const TrapezoidalPlaneBounds* bounds = dynamic_cast(&geomDet->surface().bounds()); + + float wideWidth = bounds->width(); + float narrowWidth = 2.f * bounds->widthAtHalfLength() - wideWidth; + float length = bounds->length(); + float tangent = (wideWidth - narrowWidth) / (2.f * length); + float halfWidthAtY = tangent * localPoint.y() + 0.25f * (narrowWidth + wideWidth); + + distanceX = std::abs(localPoint.x()) - halfWidthAtY; + distanceY = std::abs(localPoint.y()) - 0.5f * length; + } else { + distanceX = std::abs(localPoint.x()) - 0.5f * geomDet->surface().bounds().width(); + distanceY = std::abs(localPoint.y()) - 0.5f * geomDet->surface().bounds().length(); + } + if ((distanceX < parameters.muonMaxDistanceX && distanceY < parameters.muonMaxDistanceY) || + (distanceX * distanceX < + localError.xx() * parameters.muonMaxDistanceSigmaX * parameters.muonMaxDistanceSigmaX && + distanceY * distanceY < + localError.yy() * parameters.muonMaxDistanceSigmaY * parameters.muonMaxDistanceSigmaY)) { + LogTrace("TrackAssociator") << "found a match: " << DetIdInfo::info(*detId, nullptr) << "\n"; + TAMuonChamberMatch match; + match.tState = stateOnSurface; + match.localDistanceX = distanceX; + match.localDistanceY = distanceY; + match.id = *detId; + matches.push_back(match); + } else { + LogTrace("TrackAssociator") << "chamber is too far: " << DetIdInfo::info(*detId, nullptr) + << "\n\tdistanceX: " << distanceX << "\t distanceY: " << distanceY + << "\t sigmaX: " << distanceX / sqrt(localError.xx()) + << "\t sigmaY: " << distanceY / sqrt(localError.yy()) << "\n"; + } + } +} + +void TrackDetectorAssociator::fillMuon(const edm::Event& iEvent, + TrackDetMatchInfo& info, + const AssociatorParameters& parameters) { + // Get the segments from the event + edm::Handle dtSegments; + iEvent.getByToken(parameters.dtSegmentsToken, dtSegments); + if (!dtSegments.isValid()) + throw cms::Exception("FatalError") << "Unable to find DTRecSegment4DCollection in event!\n"; + + edm::Handle cscSegments; + iEvent.getByToken(parameters.cscSegmentsToken, cscSegments); + if (!cscSegments.isValid()) + throw cms::Exception("FatalError") << "Unable to find CSCSegmentCollection in event!\n"; + + edm::Handle gemSegments; + if (parameters.useGEM) + iEvent.getByToken(parameters.gemSegmentsToken, gemSegments); + edm::Handle me0Segments; + if (parameters.useME0) + iEvent.getByToken(parameters.me0SegmentsToken, me0Segments); + + // Get the hits from the event only if track preselection is activated + // Get the chambers segments/hits in the events + std::set occupancy_set; + if (parameters.preselectMuonTracks) { + edm::Handle rpcRecHits; + iEvent.getByToken(parameters.rpcHitsToken, rpcRecHits); + if (!rpcRecHits.isValid()) + throw cms::Exception("FatalError") << "Unable to find RPCRecHitCollection in event!\n"; + + edm::Handle gemRecHits; + if (parameters.useGEM) + iEvent.getByToken(parameters.gemHitsToken, gemRecHits); + + edm::Handle me0RecHits; + if (parameters.useME0) + iEvent.getByToken(parameters.me0HitsToken, me0RecHits); + + for (const auto& dtSegment : *dtSegments) { + occupancy_set.insert(dtSegment.geographicalId()); + } + for (const auto& cscSegment : *cscSegments) { + occupancy_set.insert(cscSegment.geographicalId()); + } + for (const auto& rpcRecHit : *rpcRecHits) { + occupancy_set.insert(rpcRecHit.geographicalId()); + } + if (parameters.useGEM) { + for (const auto& gemSegment : *gemSegments) { + occupancy_set.insert(gemSegment.geographicalId()); + } + for (const auto& gemRecHit : *gemRecHits) { + occupancy_set.insert(gemRecHit.geographicalId()); + } + } + if (parameters.useME0) { + for (const auto& me0Segment : *me0Segments) { + occupancy_set.insert(me0Segment.geographicalId()); + } + for (const auto& me0RecHit : *me0RecHits) { + occupancy_set.insert(me0RecHit.geographicalId()); + } + } + if (occupancy_set.empty()) { + LogTrace("TrackAssociator") << "No segments or hits were found in the event: aborting track extrapolation" + << std::endl; + return; + } + } + + ///// get a set of DetId's in a given direction + + // check the map of available segments + // if there is no segments in a given direction at all, + // then there is no point to fly there. + // + // MISSING + // Possible solution: quick search for presence of segments + // for the set of DetIds + + // get a set of matches corresponding to muon chambers + std::vector matchedChambers; + LogTrace("TrackAssociator") << "Trying to Get ChamberMatches" << std::endl; + getTAMuonChamberMatches(matchedChambers, parameters, occupancy_set); + LogTrace("TrackAssociator") << "Chambers matched: " << matchedChambers.size() << "\n"; + + // Iterate over all chamber matches and fill segment matching + // info if it's available + for (std::vector::iterator matchedChamber = matchedChambers.begin(); + matchedChamber != matchedChambers.end(); + matchedChamber++) { + const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet((*matchedChamber).id); + // DT chamber + if (const DTChamber* chamber = dynamic_cast(geomDet)) { + // Get the range for the corresponding segments + DTRecSegment4DCollection::range range = dtSegments->get(chamber->id()); + // Loop over the segments of this chamber + for (DTRecSegment4DCollection::const_iterator segment = range.first; segment != range.second; segment++) { + if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) { + matchedChamber->segments.back().dtSegmentRef = DTRecSegment4DRef(dtSegments, segment - dtSegments->begin()); + } + } + } + // CSC Chamber + else if (const CSCChamber* chamber = dynamic_cast(geomDet)) { + // Get the range for the corresponding segments + CSCSegmentCollection::range range = cscSegments->get(chamber->id()); + // Loop over the segments + for (CSCSegmentCollection::const_iterator segment = range.first; segment != range.second; segment++) { + if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) { + matchedChamber->segments.back().cscSegmentRef = CSCSegmentRef(cscSegments, segment - cscSegments->begin()); + } + } + } else { + // GEM Chamber + if (parameters.useGEM) { + if (const GEMSuperChamber* chamber = dynamic_cast(geomDet)) { + // Get the range for the corresponding segments + GEMSegmentCollection::range range = gemSegments->get(chamber->id()); + // Loop over the segments + for (GEMSegmentCollection::const_iterator segment = range.first; segment != range.second; segment++) { + if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) { + matchedChamber->segments.back().gemSegmentRef = + GEMSegmentRef(gemSegments, segment - gemSegments->begin()); + } + } + } + } + // ME0 Chamber + if (parameters.useME0) { + if (const ME0Chamber* chamber = dynamic_cast(geomDet)) { + // Get the range for the corresponding segments + ME0SegmentCollection::range range = me0Segments->get(chamber->id()); + // Loop over the segments + for (ME0SegmentCollection::const_iterator segment = range.first; segment != range.second; segment++) { + if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) { + matchedChamber->segments.back().me0SegmentRef = + ME0SegmentRef(me0Segments, segment - me0Segments->begin()); + } + } + } + } + } + info.chambers.push_back(*matchedChamber); + } +} + +bool TrackDetectorAssociator::addTAMuonSegmentMatch(TAMuonChamberMatch& matchedChamber, + const RecSegment* segment, + const AssociatorParameters& parameters) { + LogTrace("TrackAssociator") << "Segment local position: " << segment->localPosition() << "\n" + << std::hex << segment->geographicalId().rawId() << "\n"; + + const GeomDet* chamber = muonDetIdAssociator_->getGeomDet(matchedChamber.id); + TrajectoryStateOnSurface trajectoryStateOnSurface = matchedChamber.tState; + GlobalPoint segmentGlobalPosition = chamber->toGlobal(segment->localPosition()); + + LogTrace("TrackAssociator") << "Segment global position: " << segmentGlobalPosition + << " \t (R_xy,eta,phi): " << segmentGlobalPosition.perp() << "," + << segmentGlobalPosition.eta() << "," << segmentGlobalPosition.phi() << "\n"; + + LogTrace("TrackAssociator") << "\teta hit: " << segmentGlobalPosition.eta() + << " \tpropagator: " << trajectoryStateOnSurface.freeState()->position().eta() << "\n" + << "\tphi hit: " << segmentGlobalPosition.phi() + << " \tpropagator: " << trajectoryStateOnSurface.freeState()->position().phi() + << std::endl; + + bool isGood = false; + bool isDTWithoutY = false; + const DTRecSegment4D* dtseg = dynamic_cast(segment); + if (dtseg && (!dtseg->hasZed())) + isDTWithoutY = true; + + float deltaPhi(std::abs(segmentGlobalPosition.phi() - trajectoryStateOnSurface.freeState()->position().phi())); + if (deltaPhi > float(M_PI)) + deltaPhi = std::abs(deltaPhi - float(M_PI) * 2.f); + float deltaEta = std::abs(segmentGlobalPosition.eta() - trajectoryStateOnSurface.freeState()->position().eta()); + + if (isDTWithoutY) { + isGood = deltaPhi < parameters.dRMuon; + // Be in chamber + isGood &= deltaEta < .3f; + } else + isGood = deltaEta * deltaEta + deltaPhi * deltaPhi < parameters.dRMuon * parameters.dRMuon; + + if (isGood) { + TAMuonSegmentMatch muonSegment; + muonSegment.segmentGlobalPosition = getPoint(segmentGlobalPosition); + muonSegment.segmentLocalPosition = getPoint(segment->localPosition()); + muonSegment.segmentLocalDirection = getVector(segment->localDirection()); + muonSegment.segmentLocalErrorXX = segment->localPositionError().xx(); + muonSegment.segmentLocalErrorYY = segment->localPositionError().yy(); + muonSegment.segmentLocalErrorXY = segment->localPositionError().xy(); + muonSegment.segmentLocalErrorDxDz = segment->localDirectionError().xx(); + muonSegment.segmentLocalErrorDyDz = segment->localDirectionError().yy(); + + // DANGEROUS - compiler cannot guaranty parameters ordering + // AlgebraicSymMatrix segmentCovMatrix = segment->parametersError(); + // muonSegment.segmentLocalErrorXDxDz = segmentCovMatrix[2][0]; + // muonSegment.segmentLocalErrorYDyDz = segmentCovMatrix[3][1]; + muonSegment.segmentLocalErrorXDxDz = -999.f; + muonSegment.segmentLocalErrorYDyDz = -999.f; + muonSegment.hasZed = true; + muonSegment.hasPhi = true; + + // timing information + muonSegment.t0 = 0.f; + if (dtseg) { + if ((dtseg->hasPhi()) && (!isDTWithoutY)) { + int phiHits = dtseg->phiSegment()->specificRecHits().size(); + // int zHits = dtseg->zSegment()->specificRecHits().size(); + int hits = 0; + double t0 = 0.; + // TODO: cuts on hit numbers not optimized in any way yet... + if (phiHits > 5 && dtseg->phiSegment()->ist0Valid()) { + t0 += dtseg->phiSegment()->t0() * phiHits; + hits += phiHits; + LogTrace("TrackAssociator") << " Phi t0: " << dtseg->phiSegment()->t0() << " hits: " << phiHits; + } + // the z segments seem to contain little useful information... + // if (zHits>3) { + // t0+=s->zSegment()->t0()*zHits; + // hits+=zHits; + // LogTrace("TrackAssociator") << " Z t0: " << s->zSegment()->t0() << " hits: " << zHits << std::endl; + // } + if (hits) + muonSegment.t0 = t0 / hits; + // LogTrace("TrackAssociator") << " --- t0: " << muonSegment.t0 << std::endl; + } else { + // check and set dimensionality + if (isDTWithoutY) + muonSegment.hasZed = false; + if (!dtseg->hasPhi()) + muonSegment.hasPhi = false; + } + } + matchedChamber.segments.push_back(muonSegment); + } + + return isGood; +} + +//********************** NON-CORE CODE ******************************// + +void TrackDetectorAssociator::fillCaloTruth(const edm::Event& iEvent, + TrackDetMatchInfo& info, + const AssociatorParameters& parameters) { + // get list of simulated tracks and their vertices + using namespace edm; + Handle simTracks; + iEvent.getByToken(parameters.simTracksToken, simTracks); + if (!simTracks.isValid()) + throw cms::Exception("FatalError") << "No simulated tracks found\n"; + + Handle simVertices; + iEvent.getByToken(parameters.simVerticesToken, simVertices); + if (!simVertices.isValid()) + throw cms::Exception("FatalError") << "No simulated vertices found\n"; + + // get sim calo hits + Handle simEcalHitsEB; + iEvent.getByToken(parameters.simEcalHitsEBToken, simEcalHitsEB); + if (!simEcalHitsEB.isValid()) + throw cms::Exception("FatalError") << "No simulated ECAL EB hits found\n"; + + Handle simEcalHitsEE; + iEvent.getByToken(parameters.simEcalHitsEEToken, simEcalHitsEE); + if (!simEcalHitsEE.isValid()) + throw cms::Exception("FatalError") << "No simulated ECAL EE hits found\n"; + + Handle simHcalHits; + iEvent.getByToken(parameters.simHcalHitsToken, simHcalHits); + if (!simHcalHits.isValid()) + throw cms::Exception("FatalError") << "No simulated HCAL hits found\n"; + + // find truth partner + SimTrackContainer::const_iterator simTrack = simTracks->begin(); + for (; simTrack != simTracks->end(); ++simTrack) { + math::XYZVector simP3(simTrack->momentum().x(), simTrack->momentum().y(), simTrack->momentum().z()); + math::XYZVector recoP3(info.stateAtIP.momentum().x(), info.stateAtIP.momentum().y(), info.stateAtIP.momentum().z()); + if (ROOT::Math::VectorUtil::DeltaR(recoP3, simP3) < 0.1) + break; + } + if (simTrack != simTracks->end()) { + info.simTrack = &(*simTrack); + float ecalTrueEnergy(0); + float hcalTrueEnergy(0); + + // loop over calo hits + for (PCaloHitContainer::const_iterator hit = simEcalHitsEB->begin(); hit != simEcalHitsEB->end(); ++hit) + if (hit->geantTrackId() == info.simTrack->genpartIndex()) + ecalTrueEnergy += hit->energy(); + + for (PCaloHitContainer::const_iterator hit = simEcalHitsEE->begin(); hit != simEcalHitsEE->end(); ++hit) + if (hit->geantTrackId() == info.simTrack->genpartIndex()) + ecalTrueEnergy += hit->energy(); + + for (PCaloHitContainer::const_iterator hit = simHcalHits->begin(); hit != simHcalHits->end(); ++hit) + if (hit->geantTrackId() == info.simTrack->genpartIndex()) + hcalTrueEnergy += hit->energy(); + + info.ecalTrueEnergy = ecalTrueEnergy; + info.hcalTrueEnergy = hcalTrueEnergy; + info.hcalTrueEnergyCorrected = hcalTrueEnergy; + if (std::abs(info.trkGlobPosAtHcal.eta()) < 1.3f) + info.hcalTrueEnergyCorrected = hcalTrueEnergy * 113.2f; + else if (std::abs(info.trkGlobPosAtHcal.eta()) < 3.0f) + info.hcalTrueEnergyCorrected = hcalTrueEnergy * 167.2f; + } +} + +TrackDetMatchInfo TrackDetectorAssociator::associate(const edm::Event& iEvent, + const edm::EventSetup& iSetup, + const reco::Track& track, + const AssociatorParameters& parameters, + Direction direction /*= Any*/) { + double currentStepSize = cachedTrajectory_.getPropagationStep(); + + const MagneticField* bField = &iSetup.getData(parameters.bFieldToken); + + if (track.extra().isNull()) { + if (direction != InsideOut) + throw cms::Exception("FatalError") << "No TrackExtra information is available and association is done with " + "something else than InsideOut track.\n" + << "Either change the parameter or provide needed data!\n"; + LogTrace("TrackAssociator") << "Track Extras not found\n"; + FreeTrajectoryState initialState = trajectoryStateTransform::initialFreeState(track, bField); + return associate(iEvent, iSetup, parameters, &initialState); // 5th argument is null pointer + } + + LogTrace("TrackAssociator") << "Track Extras found\n"; + FreeTrajectoryState innerState = trajectoryStateTransform::innerFreeState(track, bField); + FreeTrajectoryState outerState = trajectoryStateTransform::outerFreeState(track, bField); + FreeTrajectoryState referenceState = trajectoryStateTransform::initialFreeState(track, bField); + + LogTrace("TrackAssociator") << "inner track state (rho, z, phi):" << track.innerPosition().Rho() << ", " + << track.innerPosition().z() << ", " << track.innerPosition().phi() << "\n"; + LogTrace("TrackAssociator") << "innerFreeState (rho, z, phi):" << innerState.position().perp() << ", " + << innerState.position().z() << ", " << innerState.position().phi() << "\n"; + + LogTrace("TrackAssociator") << "outer track state (rho, z, phi):" << track.outerPosition().Rho() << ", " + << track.outerPosition().z() << ", " << track.outerPosition().phi() << "\n"; + LogTrace("TrackAssociator") << "outerFreeState (rho, z, phi):" << outerState.position().perp() << ", " + << outerState.position().z() << ", " << outerState.position().phi() << "\n"; + + // InsideOut first + if (crossedIP(track)) { + switch (direction) { + case InsideOut: + case Any: + return associate(iEvent, iSetup, parameters, &referenceState, &outerState); + break; + case OutsideIn: { + cachedTrajectory_.setPropagationStep(-std::abs(currentStepSize)); + TrackDetMatchInfo result = associate(iEvent, iSetup, parameters, &innerState, &referenceState); + cachedTrajectory_.setPropagationStep(currentStepSize); + return result; + break; + } + } + } else { + switch (direction) { + case InsideOut: + return associate(iEvent, iSetup, parameters, &innerState, &outerState); + break; + case OutsideIn: { + cachedTrajectory_.setPropagationStep(-std::abs(currentStepSize)); + TrackDetMatchInfo result = associate(iEvent, iSetup, parameters, &outerState, &innerState); + cachedTrajectory_.setPropagationStep(currentStepSize); + return result; + break; + } + case Any: { + // check if we deal with clear outside-in case + if (track.innerPosition().Dot(track.innerMomentum()) < 0 && + track.outerPosition().Dot(track.outerMomentum()) < 0) { + cachedTrajectory_.setPropagationStep(-std::abs(currentStepSize)); + TrackDetMatchInfo result; + if (track.innerPosition().Mag2() < track.outerPosition().Mag2()) + result = associate(iEvent, iSetup, parameters, &innerState, &outerState); + else + result = associate(iEvent, iSetup, parameters, &outerState, &innerState); + cachedTrajectory_.setPropagationStep(currentStepSize); + return result; + } + } + } + } + + // all other cases + return associate(iEvent, iSetup, parameters, &innerState, &outerState); +} + +TrackDetMatchInfo TrackDetectorAssociator::associate(const edm::Event& iEvent, + const edm::EventSetup& iSetup, + const SimTrack& track, + const SimVertex& vertex, + const AssociatorParameters& parameters) { + auto const* bField = &iSetup.getData(parameters.bFieldToken); + return associate(iEvent, iSetup, getFreeTrajectoryState(bField, track, vertex), parameters); +} + +TrackDetMatchInfo TrackDetectorAssociator::associate(const edm::Event& iEvent, + const edm::EventSetup& iSetup, + const GlobalVector& momentum, + const GlobalPoint& vertex, + const int charge, + const AssociatorParameters& parameters) { + auto const* bField = &iSetup.getData(parameters.bFieldToken); + return associate(iEvent, iSetup, getFreeTrajectoryState(bField, momentum, vertex, charge), parameters); +} + +bool TrackDetectorAssociator::crossedIP(const reco::Track& track) { + bool crossed = true; + crossed &= (track.innerPosition().rho() > 3); // something close to active volume + crossed &= (track.outerPosition().rho() > 3); // something close to active volume + crossed &= + ((track.innerPosition().x() * track.innerMomentum().x() + track.innerPosition().y() * track.innerMomentum().y() < + 0) != + (track.outerPosition().x() * track.outerMomentum().x() + track.outerPosition().y() * track.outerMomentum().y() < + 0)); + return crossed; +}