diff --git a/DPGAnalysis/MuonTools/BuildFile.xml b/DPGAnalysis/MuonTools/BuildFile.xml
new file mode 100644
index 0000000000000..6977478d3e4af
--- /dev/null
+++ b/DPGAnalysis/MuonTools/BuildFile.xml
@@ -0,0 +1,35 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/DPGAnalysis/MuonTools/interface/MuBaseFlatTableProducer.h b/DPGAnalysis/MuonTools/interface/MuBaseFlatTableProducer.h
new file mode 100644
index 0000000000000..4f9bf88862b7c
--- /dev/null
+++ b/DPGAnalysis/MuonTools/interface/MuBaseFlatTableProducer.h
@@ -0,0 +1,73 @@
+#ifndef Mu_MuBaseFlatTableProducer_h
+#define Mu_MuBaseFlatTableProducer_h
+
+/** \class MuBaseFlatTableProducer MuBaseFlatTableProducer.h DPGAnalysis/MuonTools/src/MuBaseFlatTableProducer.h
+ *
+ * Helper class defining the generic interface of a FlatTableProducer
+ *
+ * \author C. Battilana (INFN BO)
+ *
+ *
+ */
+
+#include "DPGAnalysis/MuonTools/interface/MuNtupleUtils.h"
+
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+
+#include "DataFormats/NanoAOD/interface/FlatTable.h"
+
+#include
+
+class MuBaseFlatTableProducer : public edm::stream::EDProducer<> {
+public:
+ /// Constructor
+ explicit MuBaseFlatTableProducer(const edm::ParameterSet &);
+
+ /// Configure event setup for each run
+ void beginRun(const edm::Run &run, const edm::EventSetup &config) final;
+
+ /// Fill ntuples event by event
+ void produce(edm::Event &, const edm::EventSetup &) final;
+
+ /// Empty, needed by interface
+ void endRun(const edm::Run &, const edm::EventSetup &) final {}
+
+protected:
+ /// The label name of the FlatTableProducer
+ std::string m_name;
+
+ /// Get info from the ES by run
+ virtual void getFromES(const edm::Run &run, const edm::EventSetup &environment) {}
+
+ /// Get info from the ES for a given event
+ virtual void getFromES(const edm::EventSetup &environment) {}
+
+ /// Fill ntuple
+ virtual void fillTable(edm::Event &ev) = 0;
+
+ /// Definition of default values for int variables
+ static constexpr int DEFAULT_INT_VAL{-999};
+
+ /// Definition of default values for int8 variables
+ static constexpr int8_t DEFAULT_INT8_VAL{-99};
+
+ /// Definition of default values for positive int variables
+ static constexpr int DEFAULT_INT_VAL_POS{-1};
+
+ /// Definition of default values for float variables
+ static constexpr double DEFAULT_DOUBLE_VAL{-999.0};
+
+ /// Definition of default values for positive float variables
+ static constexpr double DEFAULT_DOUBLE_VAL_POS{-1.0};
+
+ template
+ void addColumn(std::unique_ptr &table,
+ const std::string name,
+ const std::vector &vec,
+ const std::string descr) {
+ table->template addColumn>(name, vec, descr);
+ }
+};
+
+#endif
diff --git a/DPGAnalysis/MuonTools/interface/MuDigiBaseProducer.h b/DPGAnalysis/MuonTools/interface/MuDigiBaseProducer.h
new file mode 100644
index 0000000000000..27d161fa233c5
--- /dev/null
+++ b/DPGAnalysis/MuonTools/interface/MuDigiBaseProducer.h
@@ -0,0 +1,118 @@
+#ifndef MuonTools_MuDigiBaseProducer_h
+#define MuonTools_MuDigiBaseProducer_h
+
+/** \class MuDigiBaseProducer MuDigiBaseProducer.h DPGAnalysis/MuonTools/src/MuDigiBaseProducer.h
+ *
+ * Helper class defining the generic interface of a muon digi Producer
+ *
+ * \author C. Battilana (INFN BO)
+ *
+ *
+ */
+
+#include "DataFormats/MuonData/interface/MuonDigiCollection.h"
+#include "PhysicsTools/NanoAOD/interface/SimpleFlatTableProducer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include
+#include
+#include
+
+template
+class MuDigiBaseProducer : public SimpleFlatTableProducerBase> {
+ using COLLECTION = MuonDigiCollection;
+
+ using IntDetVar = FuncVariable, int>;
+ using UIntDetVar = FuncVariable, unsigned int>;
+ using Int8DetVar = FuncVariable, int8_t>;
+ using UInt8DetVar = FuncVariable, uint8_t>;
+
+ std::vector>> detIdVars_;
+
+public:
+ MuDigiBaseProducer(edm::ParameterSet const ¶ms) : SimpleFlatTableProducerBase(params) {
+ const auto &varCfgs = params.getParameter("detIdVariables");
+ const auto &varNames = varCfgs.getParameterNamesForType();
+
+ std::transform(varNames.begin(), varNames.end(), std::back_inserter(detIdVars_), [&](const auto &name) {
+ const edm::ParameterSet &varCfg = varCfgs.getParameter(name);
+ const std::string &type = varCfg.getParameter("type");
+
+ std::unique_ptr> detVarPtr;
+
+ if (type == "int") {
+ detVarPtr = std::move(std::make_unique(name, varCfg));
+ } else if (type == "uint") {
+ detVarPtr = std::move(std::make_unique(name, varCfg));
+ } else if (type == "int8") {
+ detVarPtr = std::move(std::make_unique(name, varCfg));
+ } else if (type == "uint8") {
+ detVarPtr = std::move(std::make_unique(name, varCfg));
+ } else {
+ throw cms::Exception("Configuration", "unsupported type " + type + " for variable " + name);
+ }
+
+ return detVarPtr;
+ });
+ }
+
+ ~MuDigiBaseProducer() override {}
+
+ static void fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
+ edm::ParameterSetDescription desc = SimpleFlatTableProducerBase::baseDescriptions();
+
+ edm::ParameterSetDescription variable;
+ edm::Comment comType{"the c++ type of the branch in the flat table"};
+ edm::Comment comPrecision{"the precision with which to store the value in the flat table"};
+
+ variable.add("expr")->setComment("a function to define the content of the branch in the flat table");
+ variable.add("doc")->setComment("few words description of the branch content");
+
+ variable.ifValue(edm::ParameterDescription("type", "int", true, comType),
+ edm::allowedValues("int", "uint", "int8", "uint8"));
+
+ edm::ParameterSetDescription variables;
+
+ variables.setComment("a parameters set to define all variable taken form detId to fill the flat table");
+
+ edm::ParameterWildcard variableWildCard{"*", edm::RequireZeroOrMore, true, variable};
+ variables.addNode(variableWildCard);
+
+ desc.add("detIdVariables", variables);
+
+ descriptions.addWithDefaultLabel(desc);
+ }
+
+ std::unique_ptr fillTable(const edm::Event &iEvent,
+ const edm::Handle &prod) const override {
+ std::vector digis;
+ std::vector detIds;
+ std::list detIdObjs; // CB needed to store DetIds (they are transient)
+
+ if (prod.isValid()) {
+ auto detIdIt = prod->begin();
+ auto detIdEnd = prod->end();
+
+ for (; detIdIt != detIdEnd; ++detIdIt) {
+ const auto &[detId, range] = (*detIdIt);
+ detIdObjs.push_back(detId);
+ std::fill_n(std::back_inserter(detIds), range.second - range.first, &detIdObjs.back());
+ std::transform(range.first, range.second, std::back_inserter(digis), [](const auto &digi) { return &digi; });
+ }
+ }
+
+ auto table = std::make_unique(digis.size(), this->name_, false, this->extension_);
+
+ for (const auto &var : this->vars_) {
+ var->fill(digis, *table);
+ }
+
+ for (const auto &var : detIdVars_) {
+ var->fill(detIds, *table);
+ }
+
+ return table;
+ }
+};
+
+#endif
diff --git a/DPGAnalysis/MuonTools/interface/MuLocalRecoBaseProducer.h b/DPGAnalysis/MuonTools/interface/MuLocalRecoBaseProducer.h
new file mode 100644
index 0000000000000..318ba3c60477e
--- /dev/null
+++ b/DPGAnalysis/MuonTools/interface/MuLocalRecoBaseProducer.h
@@ -0,0 +1,212 @@
+#ifndef MuonTools_MuRecObjBaseProducer_h
+#define MuonTools_MuRecObjBaseProducer_h
+
+/** \class MuRecObjBaseProducer MuRecObjBaseProducer.h DPGAnalysis/MuonTools/src/MuRecObjBaseProducer.h
+ *
+ * Helper class defining the generic interface of a muon digi Producer
+ *
+ * \author C. Battilana (INFN BO)
+ *
+ *
+ */
+
+#include "DataFormats/Common/interface/RangeMap.h"
+#include "DataFormats/Common/interface/ClonePolicy.h"
+#include "DataFormats/Common/interface/OwnVector.h"
+
+#include "PhysicsTools/NanoAOD/interface/SimpleFlatTableProducer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "DataFormats/TrackingRecHit/interface/RecSegment.h"
+
+#include "FWCore/Framework/interface/ESHandle.h"
+
+#include "FWCore/Utilities/interface/ESGetToken.h"
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
+#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
+
+#include
+#include
+#include
+#include
+
+template
+class MuRecObjBaseProducer
+ : public SimpleFlatTableProducerBase>> {
+ using COLLECTION = edm::RangeMap>;
+
+ edm::ESGetToken m_token;
+ edm::ESHandle m_geometry;
+
+ using IntDetVar = FuncVariable, int>;
+ using UIntDetVar = FuncVariable, unsigned int>;
+ using Int8DetVar = FuncVariable, int8_t>;
+ using UInt8DetVar = FuncVariable, uint8_t>;
+
+ std::vector>> detIdVars_;
+
+ using GlobalPosVar = FuncVariable, float>;
+ using GlobalDirVar = FuncVariable, float>;
+
+ std::vector>> globalPosVars_;
+ std::vector>> globalDirVars_;
+
+public:
+ MuRecObjBaseProducer(edm::ParameterSet const ¶ms)
+ : SimpleFlatTableProducerBase(params), m_token{this->template esConsumes()} {
+ auto varCfgs = params.getParameter("detIdVariables");
+ auto varNames = varCfgs.getParameterNamesForType();
+
+ std::transform(varNames.begin(), varNames.end(), std::back_inserter(detIdVars_), [&](const auto &name) {
+ const edm::ParameterSet &varCfg = varCfgs.getParameter(name);
+ const std::string &type = varCfg.getParameter("type");
+
+ std::unique_ptr> detVarPtr;
+
+ if (type == "int") {
+ detVarPtr = std::move(std::make_unique(name, varCfg)); // CB can improve?
+ } else if (type == "uint") {
+ detVarPtr = std::move(std::make_unique(name, varCfg));
+ } else if (type == "int8") {
+ detVarPtr = std::move(std::make_unique(name, varCfg));
+ } else if (type == "uint8") {
+ detVarPtr = std::move(std::make_unique(name, varCfg));
+ } else {
+ throw cms::Exception("Configuration", "unsupported type " + type + " for variable " + name);
+ }
+
+ return detVarPtr;
+ });
+
+ varCfgs = params.getParameter("globalPosVariables");
+ varNames = varCfgs.getParameterNamesForType();
+
+ std::transform(varNames.begin(), varNames.end(), std::back_inserter(globalPosVars_), [&](const auto &name) {
+ return std::make_unique(name, varCfgs.getParameter(name));
+ });
+
+ if constexpr (std::is_base_of_v) {
+ varCfgs = params.getParameter("globalDirVariables");
+ varNames = varCfgs.getParameterNamesForType();
+
+ std::transform(varNames.begin(), varNames.end(), std::back_inserter(globalDirVars_), [&](const auto &name) {
+ return std::make_unique(name, varCfgs.getParameter(name));
+ });
+ }
+ }
+
+ ~MuRecObjBaseProducer() override {}
+
+ static void fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
+ edm::ParameterSetDescription desc = SimpleFlatTableProducerBase::baseDescriptions();
+
+ auto baseDescription = []() {
+ edm::ParameterSetDescription varBase;
+
+ varBase.add("expr")->setComment("a function to define the content of the branch in the flat table");
+ varBase.add("doc")->setComment("few words description of the branch content");
+
+ return varBase;
+ };
+
+ auto fullDescription = [](auto const &var, std::string const label) {
+ edm::ParameterSetDescription fullDesc;
+
+ edm::ParameterWildcard detIdVarWildCard{"*", edm::RequireZeroOrMore, true, var};
+ fullDesc.setComment("a parameters set to define all " + label + " variables to the flat table");
+ fullDesc.addNode(detIdVarWildCard);
+
+ return fullDesc;
+ };
+
+ auto detIdVar{baseDescription()};
+ auto globalGeomVar{baseDescription()};
+
+ edm::Comment comType{"the c++ type of the branch in the flat table"};
+ detIdVar.ifValue(edm::ParameterDescription{"type", "int", true, comType},
+ edm::allowedValues("int", "uint", "int8", "uint8"));
+
+ edm::Comment comPrecision{"the precision with which to store the value in the flat table"};
+ globalGeomVar.addOptionalNode(edm::ParameterDescription{"precision", true, comPrecision}, false);
+
+ desc.add("detIdVariables", fullDescription(detIdVar, "DetId"));
+ desc.add("globalPosVariables", fullDescription(globalGeomVar, "Global Position"));
+
+ if constexpr (std::is_base_of_v) {
+ desc.add("globalDirVariables", fullDescription(globalGeomVar, "Global Direction"));
+ }
+
+ descriptions.addWithDefaultLabel(desc);
+ }
+
+ std::unique_ptr fillTable(const edm::Event &iEvent,
+ const edm::Handle &product) const override {
+ std::vector objs;
+ std::vector detIds;
+ std::vector globalPositions;
+ std::vector globalDirections;
+
+ // CB needed to store DetIds, global points and vectors (they are transient)
+ std::list detIdObjs;
+ std::list globalPointObjs;
+ std::list globalVectorObjs;
+
+ if (product.isValid()) {
+ auto detIdIt = product->id_begin();
+ const auto detIdEnd = product->id_end();
+
+ for (; detIdIt != detIdEnd; ++detIdIt) {
+ const auto &range = product->get(*detIdIt);
+ const GeomDet *geomDet = m_geometry->idToDet(*detIdIt);
+
+ detIdObjs.push_back(*detIdIt);
+ std::fill_n(std::back_inserter(detIds), range.second - range.first, &detIdObjs.back());
+
+ for (auto objIt{range.first}; objIt != range.second; ++objIt) {
+ objs.push_back(&(*objIt));
+ globalPointObjs.push_back(geomDet->toGlobal(objIt->localPosition()));
+ globalPositions.push_back(&globalPointObjs.back());
+ if constexpr (std::is_base_of_v) {
+ globalVectorObjs.push_back(geomDet->toGlobal(objIt->localDirection()));
+ globalDirections.push_back(&globalVectorObjs.back());
+ }
+ }
+ }
+ }
+
+ auto table = std::make_unique(objs.size(), this->name_, false, this->extension_);
+
+ for (const auto &var : this->vars_) {
+ var->fill(objs, *table);
+ }
+
+ for (const auto &var : detIdVars_) {
+ var->fill(detIds, *table);
+ }
+
+ for (const auto &var : globalPosVars_) {
+ var->fill(globalPositions, *table);
+ }
+
+ if constexpr (std::is_base_of_v) {
+ for (const auto &var : globalDirVars_) {
+ var->fill(globalDirections, *table);
+ }
+ }
+
+ return table;
+ }
+
+ void produce(edm::Event &event, const edm::EventSetup &environment) override {
+ edm::Handle src;
+ event.getByToken(this->src_, src);
+
+ m_geometry = environment.getHandle(m_token);
+ std::unique_ptr out = fillTable(event, src);
+ out->setDoc(this->doc_);
+
+ event.put(std::move(out));
+ }
+};
+
+#endif
diff --git a/DPGAnalysis/MuonTools/interface/MuNtupleUtils.h b/DPGAnalysis/MuonTools/interface/MuNtupleUtils.h
new file mode 100644
index 0000000000000..71116cf55956d
--- /dev/null
+++ b/DPGAnalysis/MuonTools/interface/MuNtupleUtils.h
@@ -0,0 +1,135 @@
+#ifndef MuNtuple_MuNtupleUtils_h
+#define MuNtuple_MuNtupleUtils_h
+
+/** \class MuNtupleUtils MuNtupleUtils.h MuDPGAnalysis/MuNtuples/src/MuNtupleUtils.h
+ *
+ * A set of helper classes class to handle :
+ * - Handing of InputTags and tokens
+ * - DB information from edm::EventSetup
+ * - Conversion between L1T trigger primitive coordinates and CMS global ones
+ *
+ * \author C. Battilana - INFN (BO)
+ * \author L. Giuducci - INFN (BO)
+ *
+ *
+ */
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+#include "Geometry/DTGeometry/interface/DTGeometry.h"
+#include "Geometry/Records/interface/MuonGeometryRecord.h"
+
+#include
+#include
+
+class L1MuDTChambPhDigi;
+class L1Phase2MuDTPhDigi;
+
+namespace nano_mu {
+
+ template
+ class EDTokenHandle {
+ public:
+ /// Constructor
+ EDTokenHandle(const edm::ParameterSet& config, edm::ConsumesCollector&& collector, std::string name)
+ : m_name{name}, m_inputTag{config.getParameter(name)} {
+ if (m_inputTag.label() != "none") {
+ m_token = collector.template consumes(m_inputTag);
+ }
+ }
+
+ /// Conditional getter
+ /// checks whether a token is valid and if
+ /// retireving the data collection succeded
+ auto conditionalGet(const edm::Event& ev) const {
+ edm::Handle collection;
+
+ if (!m_token.isUninitialized() && !ev.getByToken(m_token, collection))
+ edm::LogError("") << "[EDTokenHandle]::conditionalGet: " << m_inputTag.label()
+ << " collection does not exist !!!";
+
+ return collection;
+ }
+
+ private:
+ std::string m_name;
+ edm::InputTag m_inputTag;
+ edm::EDGetTokenT m_token;
+ };
+
+ template
+ class ESTokenHandle {
+ public:
+ /// Constructor
+ ESTokenHandle(edm::ConsumesCollector&& collector, const std::string& label = "")
+ : m_token{collector.template esConsumes(edm::ESInputTag{"", label})} {}
+
+ ///Get Handle from ES
+ void getFromES(const edm::EventSetup& environment) { m_handle = environment.getHandle(m_token); }
+
+ /// Check validity
+ bool isValid() { return m_handle.isValid(); }
+
+ /// Return handle
+ T const* operator->() { return m_handle.product(); }
+
+ private:
+ edm::ESGetToken m_token;
+ edm::ESHandle m_handle;
+ };
+
+ class DTTrigGeomUtils {
+ public:
+ struct chambCoord {
+ double pos{};
+ double dir{};
+ };
+
+ /// Constructor
+ DTTrigGeomUtils(edm::ConsumesCollector&& collector, bool dirInDeg = true);
+
+ /// Return local position and direction in chamber RF - legacy
+ chambCoord trigToReco(const L1MuDTChambPhDigi* trig);
+
+ /// Return local position and direction in chamber RF - analytical method
+ chambCoord trigToReco(const L1Phase2MuDTPhDigi* trig);
+
+ /// Checks id the chamber has positive RF;
+ bool hasPosRF(int wh, int sec) { return wh > 0 || (wh == 0 && sec % 4 > 1); };
+
+ /// Update EventSetup information
+ void getFromES(const edm::Run& run, const edm::EventSetup& environment) {
+ m_dtGeom.getFromES(environment);
+ for (int i_st = 0; i_st != 4; ++i_st) {
+ const DTChamberId chId(-2, i_st + 1, 4);
+ const DTChamber* chamb = m_dtGeom->chamber(chId);
+ const DTSuperLayer* sl1 = chamb->superLayer(DTSuperLayerId(chId, 1));
+ const DTSuperLayer* sl3 = chamb->superLayer(DTSuperLayerId(chId, 3));
+ m_zsl1[i_st] = chamb->surface().toLocal(sl1->position()).z();
+ m_zsl3[i_st] = chamb->surface().toLocal(sl3->position()).z();
+ m_zcn[i_st] = 0.5 * (m_zsl1[i_st] + m_zsl3[i_st]);
+ }
+ };
+
+ private:
+ ESTokenHandle m_dtGeom;
+
+ std::array m_zcn;
+ std::array m_zsl1;
+ std::array m_zsl3;
+
+ static constexpr double PH1_PHI_R = 4096.;
+ static constexpr double PH1_PHIB_R = 512.;
+
+ static constexpr double PH2_PHI_R = 65536. / 0.8;
+ static constexpr double PH2_PHIB_R = 2048. / 1.4;
+ };
+
+} // namespace nano_mu
+
+#endif
diff --git a/DPGAnalysis/MuonTools/plugins/BuildFile.xml b/DPGAnalysis/MuonTools/plugins/BuildFile.xml
new file mode 100644
index 0000000000000..ccd1c6c40c741
--- /dev/null
+++ b/DPGAnalysis/MuonTools/plugins/BuildFile.xml
@@ -0,0 +1,18 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/DPGAnalysis/MuonTools/plugins/MuCSCTnPFlatTableProducer.cc b/DPGAnalysis/MuonTools/plugins/MuCSCTnPFlatTableProducer.cc
new file mode 100644
index 0000000000000..368e2474733c8
--- /dev/null
+++ b/DPGAnalysis/MuonTools/plugins/MuCSCTnPFlatTableProducer.cc
@@ -0,0 +1,1031 @@
+/** \class MuCSCTnPFlatTableProducer MuCSCTnPFlatTableProducer.cc DPGAnalysis/MuonTools/plugins/MuCSCTnPFlatTableProducer.cc
+ *
+ * Helper class : the CSC Tag and probe segment efficiency filler
+ *
+ * \author M. Herndon (UW Madison)
+ *
+ *
+ */
+
+#include "DataFormats/DetId/interface/DetId.h"
+#include "DataFormats/CSCRecHit/interface/CSCSegment.h"
+#include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
+#include "DataFormats/MuonDetId/interface/CSCDetId.h"
+
+#include "DataFormats/MuonReco/interface/MuonSelectors.h"
+#include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
+
+#include "DataFormats/Math/interface/deltaR.h"
+
+#include "TString.h"
+#include "TRegexp.h"
+
+#include
+
+#include
+#include
+
+#include "DPGAnalysis/MuonTools/interface/MuBaseFlatTableProducer.h"
+
+#include "FWCore/Framework/interface/ESHandle.h"
+
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+
+#include "Geometry/CSCGeometry/interface/CSCGeometry.h"
+#include "Geometry/Records/interface/MuonGeometryRecord.h"
+
+#include "DataFormats/MuonReco/interface/Muon.h"
+#include "DataFormats/MuonReco/interface/MuonFwd.h"
+#include "DataFormats/MuonReco/interface/MuonIsolation.h"
+#include "DataFormats/MuonReco/interface/MuonPFIsolation.h"
+
+#include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/TrackReco/interface/TrackFwd.h"
+
+#include "DataFormats/VertexReco/interface/Vertex.h"
+#include "DataFormats/VertexReco/interface/VertexFwd.h"
+
+#include "DataFormats/Common/interface/TriggerResults.h"
+#include "DataFormats/HLTReco/interface/TriggerEvent.h"
+#include "DataFormats/HLTReco/interface/TriggerObject.h"
+#include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
+
+#include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
+#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
+#include "TrackingTools/Records/interface/TransientTrackRecord.h"
+#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
+#include "TrackingTools/GeomPropagators/interface/Propagator.h"
+
+class MuonServiceProxy;
+
+class MuCSCTnPFlatTableProducer : public MuBaseFlatTableProducer {
+public:
+ /// Constructor
+ MuCSCTnPFlatTableProducer(const edm::ParameterSet&);
+
+ /// Fill descriptors
+ static void fillDescriptions(edm::ConfigurationDescriptions&);
+
+protected:
+ /// Fill tree branches for a given events
+ void fillTable(edm::Event&) final;
+
+ /// Get info from the ES by run
+ void getFromES(const edm::Run&, const edm::EventSetup&) final;
+
+ /// Get info from the ES for a given event
+ void getFromES(const edm::EventSetup&) final;
+
+private:
+ static constexpr Float_t MEZ[6] = {601.3, 696.11, 696.11, 827.56, 936.44, 1025.9};
+
+ /// Tokens
+ nano_mu::EDTokenHandle m_muToken;
+ nano_mu::EDTokenHandle m_trackToken;
+
+ nano_mu::EDTokenHandle m_cscSegmentToken;
+
+ nano_mu::EDTokenHandle> m_primaryVerticesToken;
+
+ nano_mu::EDTokenHandle m_trigResultsToken;
+ nano_mu::EDTokenHandle m_trigEventToken;
+
+ /// Name of the triggers used by muon filler for trigger matching
+ std::string m_trigName;
+ std::string m_isoTrigName;
+
+ /// Handles to geometry, detector and specialized objects
+ /// CSC Geometry
+ nano_mu::ESTokenHandle m_cscGeometry;
+
+ /// Muon service proxy
+ std::unique_ptr m_muonSP;
+
+ /// Transient Track Builder
+ nano_mu::ESTokenHandle m_transientTrackBuilder;
+
+ // Extrapolator to cylinder
+ edm::ESHandle propagatorAlong;
+ edm::ESHandle propagatorOpposite;
+ edm::ESHandle theBField;
+
+ /// HLT config provider
+ HLTConfigProvider m_hltConfig;
+
+ /// Indices of the triggers used by muon filler for trigger matching
+ std::vector m_trigIndices;
+ std::vector m_isoTrigIndices;
+
+ /// Selection functions
+ //bool muonTagSelection(const reco::Muon & muon,edm::Handle> tracks);
+ bool trackProbeSelection(const reco::Track& track, edm::Handle>);
+ bool muonTagSelection(const reco::Muon&);
+ //bool trackProbeSelection(const reco::Track & track);
+ bool zSelection(const reco::Muon&, const reco::Track&);
+
+ // Calculation functions
+ double zMass(const reco::Track&, const reco::Muon&);
+ double calcDeltaR(double, double, double, double);
+ double iso(const reco::Track&, edm::Handle>);
+
+ // Track extrapolation and segment match functions
+ TrajectoryStateOnSurface surfExtrapTrkSam(const reco::Track&, double);
+ FreeTrajectoryState freeTrajStateMuon(const reco::Track&);
+
+ UChar_t ringCandidate(Int_t iiStation, Int_t station, Float_t feta, Float_t phi);
+ UChar_t thisChamberCandidate(UChar_t station, UChar_t ring, Float_t phi);
+
+ TrajectoryStateOnSurface* matchTTwithCSCSeg(const reco::Track&,
+ edm::Handle,
+ CSCSegmentCollection::const_iterator&,
+ CSCDetId&);
+ Float_t TrajectoryDistToSeg(TrajectoryStateOnSurface, CSCSegmentCollection::const_iterator);
+ std::vector GetEdgeAndDistToGap(const reco::Track&, CSCDetId&);
+ Float_t YDistToHVDeadZone(Float_t, Int_t);
+
+ /// The variables holding
+ /// the T&P related information
+
+ unsigned int m_nZCands; // the # of digis (size of all following vectors)
+
+ double _trackIso;
+ double _muonIso;
+ double _zMass;
+
+ bool hasTrigger(std::vector&,
+ const trigger::TriggerObjectCollection&,
+ edm::Handle&,
+ const reco::Muon&);
+
+ float computeTrkIso(const reco::MuonIsolation&, float);
+ float computePFIso(const reco::MuonPFIsolation&, float);
+};
+
+MuCSCTnPFlatTableProducer::MuCSCTnPFlatTableProducer(const edm::ParameterSet& config)
+ : MuBaseFlatTableProducer(config),
+ m_muToken{config, consumesCollector(), "muonSrc"},
+ m_trackToken{config, consumesCollector(), "trackSrc"},
+ m_cscSegmentToken{config, consumesCollector(), "cscSegmentSrc"},
+ m_primaryVerticesToken{config, consumesCollector(), "primaryVerticesSrc"},
+ m_trigResultsToken{config, consumesCollector(), "trigResultsSrc"},
+ m_trigEventToken{config, consumesCollector(), "trigEventSrc"},
+ m_trigName{config.getParameter("trigName")},
+ m_isoTrigName{config.getParameter("isoTrigName")},
+ m_cscGeometry{consumesCollector()},
+ m_muonSP{std::make_unique(config.getParameter("ServiceParameters"),
+ consumesCollector())},
+ m_transientTrackBuilder{consumesCollector(), "TransientTrackBuilder"} {
+ produces();
+}
+
+void MuCSCTnPFlatTableProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+ edm::ParameterSetDescription desc;
+
+ desc.add("name", "cscTnP");
+ desc.add("muonSrc", edm::InputTag{"muons"});
+ desc.add("trackSrc", edm::InputTag{"generalTracks"});
+ desc.add("cscSegmentSrc", edm::InputTag{"cscSegments"});
+ desc.add("primaryVerticesSrc", edm::InputTag{"offlinePrimaryVertices"});
+
+ desc.add("trigEventSrc", edm::InputTag{"hltTriggerSummaryAOD::HLT"});
+ desc.add("trigResultsSrc", edm::InputTag{"TriggerResults::HLT"});
+
+ desc.add("trigName", "none");
+ desc.add("isoTrigName", "HLT_IsoMu2*");
+
+ desc.setAllowAnything();
+
+ descriptions.addWithDefaultLabel(desc);
+}
+
+void MuCSCTnPFlatTableProducer::getFromES(const edm::Run& run, const edm::EventSetup& environment) {
+ m_cscGeometry.getFromES(environment);
+
+ bool changed{true};
+ m_hltConfig.init(run, environment, "HLT", changed);
+
+ const bool enableWildcard{true};
+
+ TString tName = TString(m_trigName);
+ TRegexp tNamePattern = TRegexp(tName, enableWildcard);
+
+ for (unsigned iPath = 0; iPath < m_hltConfig.size(); ++iPath) {
+ TString pathName = TString(m_hltConfig.triggerName(iPath));
+ if (pathName.Contains(tNamePattern))
+ m_trigIndices.push_back(static_cast(iPath));
+ }
+
+ tName = TString(m_isoTrigName);
+ tNamePattern = TRegexp(tName, enableWildcard);
+
+ for (unsigned iPath = 0; iPath < m_hltConfig.size(); ++iPath) {
+ TString pathName = TString(m_hltConfig.triggerName(iPath));
+ if (pathName.Contains(tNamePattern))
+ m_isoTrigIndices.push_back(static_cast(iPath));
+ }
+}
+
+void MuCSCTnPFlatTableProducer::getFromES(const edm::EventSetup& environment) {
+ m_transientTrackBuilder.getFromES(environment);
+ m_muonSP->update(environment);
+}
+
+void MuCSCTnPFlatTableProducer::fillTable(edm::Event& ev) {
+ unsigned int m_nZCands = 0; // the # of digis (size of all following vectors)
+
+ // Muon track tag variables
+ std::vector m_muonPt; // muon pT [GeV/c]
+ std::vector m_muonPhi; // muon phi [rad]
+ std::vector m_muonEta; // muon eta
+ std::vector m_muonPtError; // muon pT [GeV/c] error
+ std::vector m_muonPhiError; // muon phi [rad] error
+ std::vector m_muonEtaError; // muon eta error
+ std::vector m_muonCharge; // muon charge
+ std::vector m_muonDXY; // muon dXY
+ std::vector m_muonDZ; // muon dZ
+ std::vector m_muonTrkHits; // muon track Hits
+ std::vector m_muonChi2; // muon Chi2
+ std::vector m_muonTrigger; // muon trigger
+ std::vector m_muonIso; // track Iso
+
+ // Track probe variabes
+ std::vector m_trackPt; // track pT [GeV/c]
+ std::vector m_trackP; // track P [GeV/c]
+ std::vector m_trackPhi; // track phi [rad]
+ std::vector m_trackEta; // track eta
+ std::vector m_trackPtError; // track pT [GeV/c] error
+ std::vector m_trackPhiError; // track phi [rad] error
+ std::vector m_trackEtaError; // track eta error
+ std::vector m_trackCharge; // track charge
+ std::vector m_trackDXY; // track dXY
+ std::vector m_trackDZ; // track dZ
+ std::vector m_trackTrkHits; // track Hits
+ std::vector m_trackChi2; // track Chi2
+ std::vector m_trackIso; // track Iso
+
+ // Z and global variables
+ std::vector m_zMass; // z mass
+ std::vector m_dRTrackMuon; // dR between the track and muon
+ std::vector m_numberOfPrimaryVertices; // Number of primary Vertices
+
+ // CSC chamber information, station encoded in vector
+ std::vector m_chamberEndcap; // chamber endcap
+ // station encoded in array index
+ std::array, 4> m_chamberRing; // chamber Ring
+ std::array, 4> m_chamberChamber; // chamber Chamber
+ std::array, 4> m_chamberLayer; // Segment layer information
+
+ // Track intersection variables
+ std::array, 4> m_ttIntLocalX; // track trajector intersection local X on stations 1-4
+ std::array, 4> m_ttIntLocalY; // track trajector intersection local Y on stations 1-4
+ std::array, 4> m_ttIntLocalErrorX; // track trajector intersection local X on stations 1-4
+ std::array, 4> m_ttIntLocalErrorY; // track trajector intersection local Y on stations 1-4
+ std::array, 4> m_ttIntLocalW; // track trajector intersection local Wire on stations 1-4
+ std::array, 4> m_ttIntLocalS; // track trajector intersection local Strip on stations 1-4
+ std::array, 4> m_ttIntEta; // track trajector intersection Eta stations 1-4
+
+ // Track intersection fiducial information
+
+ std::array, 4>
+ m_ttDistToEdge; // track trajector intersection distance to edge, neg is with chamber, on stations 1-4
+ std::array, 4> m_ttDistToHVGap; // track trajector intersection distance to HV GAP on stations 1-4
+
+ // Segment location variables
+ std::array, 4> m_segLocalX; // segment local X on stations 1-4
+ std::array, 4> m_segLocalY; // segment local Y on stations 1-4
+ std::array, 4> m_segLocalErrorX; // segment local X error on stations 1-4
+ std::array, 4> m_segLocalErrorY; // segment local Y error on stations 1-4
+
+ // track intersection segment residuals variables
+ std::array, 4>
+ m_ttIntSegResidualLocalX; // track trajector intersection Segment residual local X on stations 1-4
+ std::array, 4>
+ m_ttIntSegResidualLocalY; // track trajector intersection Segment residuallocal Y on stations 1-4
+
+ auto&& propagator_along = m_muonSP->propagator("SteppingHelixPropagatorAlong");
+ auto&& propagator_opposite = m_muonSP->propagator("SteppingHelixPropagatorOpposite");
+
+ propagatorAlong = propagator_along;
+ propagatorOpposite = propagator_opposite;
+
+ theBField = m_muonSP->magneticField();
+
+ auto muons = m_muToken.conditionalGet(ev);
+ auto tracks = m_trackToken.conditionalGet(ev);
+ auto segments = m_cscSegmentToken.conditionalGet(ev);
+ auto primaryVertices = m_primaryVerticesToken.conditionalGet(ev);
+
+ auto triggerResults = m_trigResultsToken.conditionalGet(ev);
+ auto triggerEvent = m_trigEventToken.conditionalGet(ev);
+
+ if (muons.isValid() && tracks.isValid() && segments.isValid() && primaryVertices.isValid() &&
+ m_transientTrackBuilder.isValid()) {
+ for (const auto& muon : (*muons)) {
+ if (!muonTagSelection(muon))
+ continue;
+
+ bool muonTrigger = false;
+ if (triggerResults.isValid() && triggerEvent.isValid()) {
+ const auto& triggerObjects = triggerEvent->getObjects();
+ muonTrigger = (hasTrigger(m_isoTrigIndices, triggerObjects, triggerEvent, muon) ||
+ hasTrigger(m_trigIndices, triggerObjects, triggerEvent, muon));
+ }
+
+ for (const auto& track : (*tracks)) {
+ if (!trackProbeSelection(track, tracks))
+ continue;
+ if (!zSelection(muon, track))
+ continue;
+ //std::cout << "Z candidate found: " << _zMass << " track eta: " << track.eta() << std::endl;
+ //std::cout.flush();
+ m_nZCands++;
+
+ m_trackPt.push_back(track.pt());
+ m_trackP.push_back(track.p());
+ m_trackEta.push_back(track.eta());
+ m_trackPhi.push_back(track.phi());
+ m_trackPtError.push_back(track.pt());
+ m_trackEtaError.push_back(track.eta());
+ m_trackPhiError.push_back(track.phi());
+ m_trackCharge.push_back(track.charge());
+ m_trackDXY.push_back(track.dxy());
+ m_trackDZ.push_back(track.dz());
+ m_trackTrkHits.push_back(track.hitPattern().numberOfValidTrackerHits());
+ m_trackChi2.push_back(track.normalizedChi2());
+ m_trackIso.push_back(_trackIso);
+
+ m_muonPt.push_back(muon.track()->pt());
+ m_muonPhi.push_back(muon.track()->phi());
+ m_muonEta.push_back(muon.track()->eta());
+ m_muonPtError.push_back(muon.track()->ptError());
+ m_muonPhiError.push_back(muon.track()->phiError());
+ m_muonEtaError.push_back(muon.track()->etaError());
+ m_muonCharge.push_back(muon.charge());
+ m_muonDXY.push_back(muon.track()->dxy());
+ m_muonDZ.push_back(muon.track()->dz());
+ m_muonTrkHits.push_back(muon.track()->hitPattern().numberOfValidTrackerHits());
+ m_muonChi2.push_back(muon.track()->normalizedChi2());
+ m_muonIso.push_back(computeTrkIso(muon.isolationR03(), muon.pt()));
+ m_muonTrigger.push_back(muonTrigger);
+
+ m_zMass.push_back(_zMass);
+ double_t dR = calcDeltaR(track.eta(), muon.eta(), track.phi(), muon.phi());
+ //double_t dR = 1.0;
+ m_dRTrackMuon.push_back(dR);
+ const reco::VertexCollection& vertices = *primaryVertices.product();
+ m_numberOfPrimaryVertices.push_back(vertices.size());
+
+ bool ec = (track.eta() > 0);
+ UChar_t endcapCSC = ec ? 0 : 1;
+ m_chamberEndcap.push_back(endcapCSC * 1);
+
+ Int_t iiStationFail = 0;
+ for (int iiStationZ = 0; iiStationZ < 6; iiStationZ++) {
+ UChar_t stationCSC = iiStationZ > 2 ? iiStationZ - 2 : 0;
+ UChar_t ringCSC = 0;
+ TrajectoryStateOnSurface tsos = surfExtrapTrkSam(track, ec ? MEZ[iiStationZ] : -MEZ[iiStationZ]);
+
+ if (tsos.isValid()) {
+ Float_t trkEta = tsos.globalPosition().eta(), trkPhi = tsos.globalPosition().phi();
+ ringCSC = ringCandidate(iiStationZ, stationCSC + 1, trkEta, trkPhi);
+
+ if (ringCSC) {
+ UChar_t chamberCSC = thisChamberCandidate(stationCSC + 1, ringCSC, track.phi()) - 1;
+ CSCDetId Layer3id = CSCDetId(endcapCSC + 1, stationCSC + 1, ringCSC, chamberCSC + 1, 3);
+ CSCDetId Layer0Id = CSCDetId(endcapCSC + 1,
+ stationCSC + 1,
+ ringCSC,
+ chamberCSC + 1,
+ 0); //layer 0 is the mid point of the chamber. It is not a real layer.
+ // !!!!! need to fix Layer0Id problem with ME1/1 here
+
+ const BoundPlane& Layer3Surface = m_cscGeometry->idToDet(Layer3id)->surface();
+
+ tsos = surfExtrapTrkSam(track, Layer3Surface.position().z());
+
+ if (tsos.isValid()) {
+ // Fill track intersection denominator information
+ LocalPoint localTTIntPoint = Layer3Surface.toLocal(tsos.freeState()->position());
+ const CSCLayerGeometry* layerGeoma = m_cscGeometry->chamber(Layer0Id)->layer(3)->geometry();
+ const CSCLayerGeometry* layerGeomb = m_cscGeometry->chamber(Layer0Id)->layer(4)->geometry();
+
+ m_chamberRing[stationCSC].push_back(ringCSC);
+ m_chamberChamber[stationCSC].push_back(chamberCSC);
+ m_ttIntLocalX[stationCSC].push_back(localTTIntPoint.x());
+ m_ttIntLocalY[stationCSC].push_back(localTTIntPoint.y());
+ m_ttIntLocalW[stationCSC].push_back(
+ (layerGeoma->nearestWire(localTTIntPoint) + layerGeomb->nearestWire(localTTIntPoint)) / 2.0);
+ m_ttIntLocalS[stationCSC].push_back(
+ (layerGeoma->strip(localTTIntPoint) + layerGeomb->strip(localTTIntPoint)) / 2.0);
+ m_ttIntEta[stationCSC].push_back(trkEta);
+
+ // Errors are those of the track intersection, chosing the plane and exact geomentry is performed in the function
+ Float_t CSCProjEdgeDist = -9999.0;
+ Float_t ttIntLocalErrorX = -9999.0;
+ Float_t CSCDyProjHVGap = 9999.0;
+ Float_t ttIntLocalErrorY = -9999.0;
+ for (Int_t ly = 1; ly < 7; ly++) {
+ CSCDetId Layerid = CSCDetId(endcapCSC + 1, stationCSC + 1, ringCSC, chamberCSC + 1, ly);
+ std::vector EdgeAndDistToGap(GetEdgeAndDistToGap(
+ track, Layerid)); //values: 1-edge;2-err of edge;3-disttogap;4-err of dist to gap
+ if (EdgeAndDistToGap[0] > CSCProjEdgeDist) {
+ CSCProjEdgeDist = EdgeAndDistToGap[0];
+ ttIntLocalErrorX = EdgeAndDistToGap[1];
+ }
+ if (EdgeAndDistToGap[2] < CSCDyProjHVGap) {
+ CSCDyProjHVGap = EdgeAndDistToGap[2];
+ ttIntLocalErrorY = EdgeAndDistToGap[3];
+ }
+ }
+ m_ttDistToEdge[stationCSC].push_back(CSCProjEdgeDist);
+ m_ttDistToHVGap[stationCSC].push_back(CSCDyProjHVGap);
+ m_ttIntLocalErrorX[stationCSC].push_back(ttIntLocalErrorX);
+ m_ttIntLocalErrorY[stationCSC].push_back(ttIntLocalErrorY);
+
+ // now we have a local point for comparison to segments
+ CSCSegmentCollection::const_iterator cscSegOut;
+ TrajectoryStateOnSurface* TrajToSeg = matchTTwithCSCSeg(track, segments, cscSegOut, Layer3id);
+
+ if (TrajToSeg == nullptr) {
+ // fill Null Num
+ m_segLocalX[stationCSC].push_back(-9999.0);
+ m_segLocalY[stationCSC].push_back(-9999.0);
+ m_segLocalErrorX[stationCSC].push_back(0.0);
+ m_segLocalErrorY[stationCSC].push_back(0.0);
+
+ m_ttIntSegResidualLocalX[stationCSC].push_back(-9990.0);
+ m_ttIntSegResidualLocalY[stationCSC].push_back(-9990.0);
+
+ m_chamberLayer[stationCSC].push_back(-9);
+
+ continue;
+ }
+
+ LocalPoint localSegmentPoint = (*cscSegOut).localPosition();
+ LocalError localSegErr = (*cscSegOut).localPositionError();
+
+ m_segLocalX[stationCSC].push_back(localSegmentPoint.x());
+ m_segLocalY[stationCSC].push_back(localSegmentPoint.y());
+ m_segLocalErrorX[stationCSC].push_back(sqrt(localSegErr.xx()));
+ m_segLocalErrorY[stationCSC].push_back(sqrt(localSegErr.yy()));
+
+ m_ttIntSegResidualLocalX[stationCSC].push_back(localTTIntPoint.x() - localSegmentPoint.x());
+ m_ttIntSegResidualLocalY[stationCSC].push_back(localTTIntPoint.y() - localSegmentPoint.y());
+ /* Extract layers for hits */
+ int layers = 0;
+ for (std::vector::const_iterator itRH = cscSegOut->specificRecHits().begin();
+ itRH != cscSegOut->specificRecHits().end();
+ ++itRH) {
+ const CSCRecHit2D* recHit = &(*itRH);
+ int layer = recHit->cscDetId().layer();
+ layers |= 1 << (layer - 1);
+ }
+ m_chamberLayer[stationCSC].push_back(layers);
+
+ } // end preliminary tsos is valid
+
+ } // end found ring CSC
+
+ } // end refined tsos is valid
+
+ if ((!tsos.isValid()) || (ringCSC == 0)) {
+ // only fill Null denominator once for station 1, iiStation Z = 0,1,2
+ if (iiStationZ <= 2)
+ iiStationFail++;
+ if (iiStationZ > 2 || iiStationFail == 3) {
+ // fill Null Den Num
+ m_chamberRing[stationCSC].push_back(-9);
+ m_chamberChamber[stationCSC].push_back(-9);
+ m_ttIntLocalX[stationCSC].push_back(-9999.0);
+ m_ttIntLocalY[stationCSC].push_back(-9999.0);
+ m_ttIntLocalErrorX[stationCSC].push_back(0.0);
+ m_ttIntLocalErrorY[stationCSC].push_back(0.0);
+ m_ttIntLocalW[stationCSC].push_back(-9999.0);
+ m_ttIntLocalS[stationCSC].push_back(-9999.0);
+ m_ttIntEta[stationCSC].push_back(-9999.0);
+
+ m_ttDistToEdge[stationCSC].push_back(-9999.0);
+ m_ttDistToHVGap[stationCSC].push_back(-9999.9);
+
+ m_segLocalX[stationCSC].push_back(-9999.0);
+ m_segLocalY[stationCSC].push_back(-9999.0);
+ m_segLocalErrorX[stationCSC].push_back(0.0);
+ m_segLocalErrorY[stationCSC].push_back(0.0);
+
+ m_ttIntSegResidualLocalX[stationCSC].push_back(-9990.0);
+ m_ttIntSegResidualLocalY[stationCSC].push_back(-9990.0);
+
+ m_chamberLayer[stationCSC].push_back(-9);
+ }
+ }
+
+ } // end loop over CSC Z planes
+ } // endl loop over tracks
+ } // end loop over muons
+
+ } // End if good physics objects
+
+ // if (m_nZCands>0) {
+ auto table = std::make_unique(m_nZCands, m_name, false, false);
+
+ table->setDoc("CSC Tag & Probe segment efficiency information");
+
+ addColumn(table, "muonPt", m_muonPt, "muon pt [GeV/c]");
+ addColumn(table, "muonPhi", m_muonPhi, "muon phi [rad]");
+ addColumn(table, "muonEta", m_muonEta, "muon eta");
+ addColumn(table, "muonPtError", m_muonPtError, "muon pt error [GeV/c]");
+ addColumn(table, "muonPhiError", m_muonPhiError, "muon phi error [rad]");
+ addColumn(table, "muonEtaError", m_muonEtaError, "muon eta error");
+ addColumn(table, "muonCharge", m_muonCharge, "muon charge");
+ addColumn(table, "muonDXY", m_muonDXY, "muon dXY [cm]");
+ addColumn(table, "muonDZ", m_muonDZ, "muon dZ [cm]");
+ addColumn(table, "muonTrkHits", m_muonTrkHits, "muon track hits");
+ addColumn(table, "muonChi2", m_muonChi2, "muon chi2");
+ addColumn(table, "muonIso", m_trackIso, "muon relative iso");
+ addColumn(table, "muonTrigger", m_muonTrigger, "muon has trigger bool");
+
+ addColumn(table, "trackPt", m_trackPt, "track pt [GeV/c]");
+ addColumn(table, "trackP", m_trackPt, "track p [GeV/c]");
+ addColumn(table, "trackPhi", m_trackPhi, "track phi [rad]");
+ addColumn(table, "trackEta", m_trackEta, "track eta");
+ addColumn(table, "trackPtError", m_trackPtError, "track pt error [GeV/c]");
+ addColumn(table, "trackPhiError", m_trackPhiError, "track phi error [rad]");
+ addColumn(table, "trackEtaError", m_trackEtaError, "track eta error");
+ addColumn(table, "trackCharge", m_trackCharge, "track charge");
+ addColumn(table, "trackDXY", m_trackDXY, "track dXY [cm]");
+ addColumn(table, "trackDZ", m_trackDZ, "track dZ [cm]");
+ addColumn(table, "trackTrkHits", m_trackTrkHits, "track track hits");
+ addColumn(table, "trackChi2", m_trackChi2, "track chi2");
+ addColumn(table, "trackIso", m_trackIso, "track relative iso");
+
+ addColumn(table, "zMass", m_zMass, "Z mass [GeV/c^2]");
+ addColumn(table, "dRTrackMuon", m_dRTrackMuon, "dR between track and muon");
+ addColumn(table, "numberOfPrimaryVertidies", m_numberOfPrimaryVertices, "Number of PVs");
+
+ addColumn(table, "chamberEndcap", m_chamberEndcap, "");
+ addColumn(table, "chamberRing1", m_chamberRing[0], "");
+ addColumn(table, "chamberRing2", m_chamberRing[1], "");
+ addColumn(table, "chamberRing3", m_chamberRing[2], "");
+ addColumn(table, "chamberRing4", m_chamberRing[3], "");
+ addColumn(table, "chamberChamber1", m_chamberChamber[0], "");
+ addColumn(table, "chamberChamber2", m_chamberChamber[1], "");
+ addColumn(table, "chamberChamber3", m_chamberChamber[2], "");
+ addColumn(table, "chamberChamber4", m_chamberChamber[3], "");
+ addColumn(table, "chamberLayer1", m_chamberLayer[0], "");
+ addColumn(table, "chamberLayer2", m_chamberLayer[1], "");
+ addColumn(table, "chamberLayer3", m_chamberLayer[2], "");
+ addColumn(table, "chamberLayer4", m_chamberLayer[3], "");
+
+ addColumn(table, "ttIntLocalX1", m_ttIntLocalX[0], "");
+ addColumn(table, "ttIntLocalX2", m_ttIntLocalX[1], "");
+ addColumn(table, "ttIntLocalX3", m_ttIntLocalX[2], "");
+ addColumn(table, "ttIntLocalX4", m_ttIntLocalX[3], "");
+ addColumn(table, "ttIntLocalY1", m_ttIntLocalY[0], "");
+ addColumn(table, "ttIntLocalY2", m_ttIntLocalY[1], "");
+ addColumn(table, "ttIntLocalY3", m_ttIntLocalY[2], "");
+ addColumn(table, "ttIntLocalY4", m_ttIntLocalY[3], "");
+ addColumn(table, "ttIntLocalErrorX1", m_ttIntLocalErrorX[0], "");
+ addColumn(table, "ttIntLocalErrorX2", m_ttIntLocalErrorX[1], "");
+ addColumn(table, "ttIntLocalErrorX3", m_ttIntLocalErrorX[2], "");
+ addColumn(table, "ttIntLocalErrorX4", m_ttIntLocalErrorX[3], "");
+ addColumn(table, "ttIntLocalErrorY1", m_ttIntLocalErrorY[0], "");
+ addColumn(table, "ttIntLocalErrorY2", m_ttIntLocalErrorY[1], "");
+ addColumn(table, "ttIntLocalErrorY3", m_ttIntLocalErrorY[2], "");
+ addColumn(table, "ttIntLocalErrorY4", m_ttIntLocalErrorY[3], "");
+ addColumn(table, "ttIntLocalW1", m_ttIntLocalW[0], "");
+ addColumn(table, "ttIntLocalW2", m_ttIntLocalW[1], "");
+ addColumn(table, "ttIntLocalW3", m_ttIntLocalW[2], "");
+ addColumn(table, "ttIntLocalW4", m_ttIntLocalW[3], "");
+ addColumn(table, "ttIntLocalS1", m_ttIntLocalS[0], "");
+ addColumn(table, "ttIntLocalS2", m_ttIntLocalS[1], "");
+ addColumn(table, "ttIntLocalS3", m_ttIntLocalS[2], "");
+ addColumn(table, "ttIntLocalS4", m_ttIntLocalS[3], "");
+ addColumn(table, "ttIntEta1", m_ttIntEta[0], "");
+ addColumn(table, "ttIntEta2", m_ttIntEta[1], "");
+ addColumn(table, "ttIntEta3", m_ttIntEta[2], "");
+ addColumn(table, "ttIntEta4", m_ttIntEta[3], "");
+
+ addColumn(table, "ttDistToEdge1", m_ttDistToEdge[0], "");
+ addColumn(table, "ttDistToEdge2", m_ttDistToEdge[1], "");
+ addColumn(table, "ttDistToEdge3", m_ttDistToEdge[2], "");
+ addColumn(table, "ttDistToEdge4", m_ttDistToEdge[3], "");
+ addColumn(table, "ttDistToHVGap1", m_ttDistToHVGap[0], "");
+ addColumn(table, "ttDistToHVGap2", m_ttDistToHVGap[1], "");
+ addColumn(table, "ttDistToHVGap3", m_ttDistToHVGap[2], "");
+ addColumn(table, "ttDistToHVGap4", m_ttDistToHVGap[3], "");
+
+ addColumn(table, "segLocalX1", m_segLocalX[0], "");
+ addColumn(table, "segLocalX2", m_segLocalX[1], "");
+ addColumn(table, "segLocalX3", m_segLocalX[2], "");
+ addColumn(table, "segLocalX4", m_segLocalX[3], "");
+ addColumn(table, "segLocalY1", m_segLocalY[0], "");
+ addColumn(table, "segLocalY2", m_segLocalY[1], "");
+ addColumn(table, "segLocalY3", m_segLocalY[2], "");
+ addColumn(table, "segLocalY4", m_segLocalY[3], "");
+ addColumn(table, "segLocalErrorX1", m_segLocalErrorX[0], "");
+ addColumn(table, "segLocalErrorX2", m_segLocalErrorX[1], "");
+ addColumn(table, "segLocalErrorX3", m_segLocalErrorX[2], "");
+ addColumn(table, "segLocalErrorX4", m_segLocalErrorX[3], "");
+ addColumn(table, "segLocalErrorY1", m_segLocalErrorY[0], "");
+ addColumn(table, "segLocalErrorY2", m_segLocalErrorY[1], "");
+ addColumn(table, "segLocalErrorY3", m_segLocalErrorY[2], "");
+ addColumn(table, "segLocalErrorY4", m_segLocalErrorY[3], "");
+
+ addColumn(table, "ttIntSegResidualLocalX1", m_ttIntSegResidualLocalX[0], "");
+ addColumn(table, "ttIntSegResidualLocalX2", m_ttIntSegResidualLocalX[1], "");
+ addColumn(table, "ttIntSegResidualLocalX3", m_ttIntSegResidualLocalX[2], "");
+ addColumn(table, "ttIntSegResidualLocalX4", m_ttIntSegResidualLocalX[3], "");
+ addColumn(table, "ttIntSegResidualLocalY1", m_ttIntSegResidualLocalY[0], "");
+ addColumn(table, "ttIntSegResidualLocalY2", m_ttIntSegResidualLocalY[1], "");
+ addColumn(table, "ttIntSegResidualLocalY3", m_ttIntSegResidualLocalY[2], "");
+ addColumn(table, "ttIntSegResidualLocalY4", m_ttIntSegResidualLocalY[3], "");
+
+ ev.put(std::move(table));
+}
+
+float MuCSCTnPFlatTableProducer::computeTrkIso(const reco::MuonIsolation& isolation, float muonPt) {
+ return isolation.sumPt / muonPt;
+}
+
+float MuCSCTnPFlatTableProducer::computePFIso(const reco::MuonPFIsolation& pfIsolation, float muonPt) {
+ return (pfIsolation.sumChargedHadronPt +
+ std::max(0., pfIsolation.sumNeutralHadronEt + pfIsolation.sumPhotonEt - 0.5 * pfIsolation.sumPUPt)) /
+ muonPt;
+}
+
+bool MuCSCTnPFlatTableProducer::hasTrigger(std::vector& trigIndices,
+ const trigger::TriggerObjectCollection& trigObjs,
+ edm::Handle& trigEvent,
+ const reco::Muon& muon) {
+ float dRMatch = 999.;
+ for (int trigIdx : trigIndices) {
+ const std::vector trigModuleLabels = m_hltConfig.moduleLabels(trigIdx);
+
+ const unsigned trigModuleIndex =
+ std::find(trigModuleLabels.begin(), trigModuleLabels.end(), "hltBoolEnd") - trigModuleLabels.begin() - 1;
+ const unsigned hltFilterIndex = trigEvent->filterIndex(edm::InputTag(trigModuleLabels[trigModuleIndex], "", "HLT"));
+ if (hltFilterIndex < trigEvent->sizeFilters()) {
+ const trigger::Keys keys = trigEvent->filterKeys(hltFilterIndex);
+ const trigger::Vids vids = trigEvent->filterIds(hltFilterIndex);
+ const unsigned nTriggers = vids.size();
+
+ for (unsigned iTrig = 0; iTrig < nTriggers; ++iTrig) {
+ trigger::TriggerObject trigObj = trigObjs[keys[iTrig]];
+ float dR = deltaR(muon, trigObj);
+ if (dR < dRMatch)
+ dRMatch = dR;
+ }
+ }
+ }
+
+ return dRMatch < 0.1; //CB should get it programmable
+}
+
+//bool MuCSCTnPFlatTableProducer::muonTagSelection(const reco::Muon & muon,edm::Handle> tracks)
+bool MuCSCTnPFlatTableProducer::muonTagSelection(const reco::Muon& muon) {
+ float ptCut = 10.0;
+ int trackerHitsCut = 8;
+ float dxyCut = 2.0;
+ float dzCut = 24.0;
+ float chi2Cut = 4.0;
+
+ bool selected = false;
+ //_muonIso = iso(*muon.track(),tracks);
+ _muonIso = computePFIso(muon.pfIsolationR04(), muon.pt());
+
+ if (!muon.isTrackerMuon())
+ return false;
+ if (!muon.track().isNonnull())
+ return false;
+ selected =
+ ((muon.track()->pt() > ptCut) && (muon.track()->hitPattern().numberOfValidTrackerHits() >= trackerHitsCut) &&
+ (muon.track()->dxy() < dxyCut) && (std::abs(muon.track()->dz()) < dzCut) &&
+ (muon.track()->normalizedChi2() < chi2Cut) && _muonIso < 0.1);
+
+ return selected;
+}
+
+bool MuCSCTnPFlatTableProducer::trackProbeSelection(const reco::Track& track,
+ edm::Handle> tracks) {
+ float ptCut = 10.0;
+ int trackerHitsCut = 8;
+ float dxyCut = 2.0;
+ float dzCut = 24.0;
+ float chi2Cut = 4.0;
+
+ bool selected = false;
+ _trackIso = iso(track, tracks);
+
+ selected =
+ ((track.pt() > ptCut) && (std::abs(track.eta()) > 0.75) && (std::abs(track.eta()) < 2.55) &&
+ (track.numberOfValidHits() >= trackerHitsCut) && (track.dxy() < dxyCut) && (std::abs(track.dz()) < dzCut) &&
+ (track.normalizedChi2() > 0.0) && (track.normalizedChi2() < chi2Cut) && _trackIso < 0.1);
+
+ return selected;
+}
+
+bool MuCSCTnPFlatTableProducer::zSelection(const reco::Muon& muon, const reco::Track& track) {
+ bool selected = false;
+
+ _zMass = zMass(track, muon);
+ selected = (track.charge() * muon.charge() == -1 && (_zMass > 75.0) && (_zMass < 120.0));
+
+ return selected;
+}
+
+// get track position at a particular (xy) plane given its z
+TrajectoryStateOnSurface MuCSCTnPFlatTableProducer::surfExtrapTrkSam(const reco::Track& track, double z) {
+ Plane::PositionType pos(0, 0, z);
+ Plane::RotationType rot;
+ Plane::PlanePointer myPlane = Plane::build(pos, rot);
+
+ FreeTrajectoryState recoStart = freeTrajStateMuon(track);
+ TrajectoryStateOnSurface recoProp = propagatorAlong->propagate(recoStart, *myPlane);
+
+ if (!recoProp.isValid())
+ recoProp = propagatorOpposite->propagate(recoStart, *myPlane);
+
+ return recoProp;
+}
+
+FreeTrajectoryState MuCSCTnPFlatTableProducer::freeTrajStateMuon(const reco::Track& track) {
+ //no track extras in nanoaod so directly use vx and p
+ GlobalPoint innerPoint(track.vx(), track.vy(), track.vz());
+ GlobalVector innerVec(track.px(), track.py(), track.pz());
+
+ GlobalTrajectoryParameters gtPars(innerPoint, innerVec, track.charge(), &*theBField);
+
+ AlgebraicSymMatrix66 cov;
+ cov *= 1e-20;
+
+ CartesianTrajectoryError tCov(cov);
+
+ return (cov.kRows == 6 ? FreeTrajectoryState(gtPars, tCov) : FreeTrajectoryState(gtPars));
+}
+
+UChar_t MuCSCTnPFlatTableProducer::ringCandidate(Int_t iiStation, Int_t station, Float_t feta, Float_t phi) {
+ UChar_t ring = 0;
+
+ switch (station) {
+ case 1:
+ if (std::abs(feta) >= 0.85 && std::abs(feta) < 1.18) { //ME13
+ if (iiStation == 2)
+ ring = 3;
+ return ring;
+ }
+ if (std::abs(feta) >= 1.18 &&
+ std::abs(feta) <= 1.5) { //ME12 if(std::abs(feta)>1.18 && std::abs(feta)<1.7){//ME12
+ if (iiStation == 1)
+ ring = 2;
+ return ring;
+ }
+ if (std::abs(feta) > 1.5 && std::abs(feta) < 2.1) { //ME11
+ if (iiStation == 0)
+ ring = 1;
+ return ring;
+ }
+ if (std::abs(feta) >= 2.1 && std::abs(feta) < 2.45) { //ME11
+ if (iiStation == 0)
+ ring = 4;
+ return ring;
+ }
+ break;
+ case 2:
+ if (std::abs(feta) > 0.95 && std::abs(feta) < 1.6) { //ME22
+ ring = 2;
+ return ring;
+ }
+ if (std::abs(feta) > 1.55 && std::abs(feta) < 2.45) { //ME21
+ ring = 1;
+ return ring;
+ }
+ break;
+ case 3:
+ if (std::abs(feta) > 1.08 && std::abs(feta) < 1.72) { //ME32
+ ring = 2;
+ return ring;
+ }
+ if (std::abs(feta) > 1.69 && std::abs(feta) < 2.45) { //ME31
+ ring = 1;
+ return ring;
+ }
+ break;
+ case 4:
+ if (std::abs(feta) > 1.78 && std::abs(feta) < 2.45) { //ME41
+ ring = 1;
+ return ring;
+ }
+ if (std::abs(feta) > 1.15 && std::abs(feta) <= 1.78) { //ME42
+ ring = 2;
+ return ring;
+ }
+ break;
+ default:
+ edm::LogError("") << "Invalid station: " << station << std::endl;
+ break;
+ }
+ return 0;
+}
+
+UChar_t MuCSCTnPFlatTableProducer::thisChamberCandidate(UChar_t station, UChar_t ring, Float_t phi) {
+ // cout <<"\t\t TPTrackMuonSys::thisChamberCandidate..."< 1 && ring == 1) ? 18 : 36;
+ const Float_t ChamberSpan = 2 * M_PI / nVal;
+ Float_t dphi = phi + M_PI / 36;
+ while (dphi >= 2 * M_PI)
+ dphi -= 2 * M_PI;
+ while (dphi < 0)
+ dphi += 2 * M_PI;
+ UChar_t ChCand = floor(dphi / ChamberSpan) + 1;
+ return ChCand > nVal ? nVal : ChCand;
+}
+
+Float_t MuCSCTnPFlatTableProducer::TrajectoryDistToSeg(TrajectoryStateOnSurface TrajSuf,
+ CSCSegmentCollection::const_iterator segIt) {
+ if (!TrajSuf.isValid())
+ return 9999.;
+ const GeomDet* gdet = m_cscGeometry->idToDet((CSCDetId)(*segIt).cscDetId());
+ LocalPoint localTTPos = gdet->surface().toLocal(TrajSuf.freeState()->position());
+ LocalPoint localSegPos = (*segIt).localPosition();
+ Float_t CSCdeltaX = localSegPos.x() - localTTPos.x();
+ Float_t CSCdeltaY = localSegPos.y() - localTTPos.y();
+ return sqrt(pow(CSCdeltaX, 2) + pow(CSCdeltaY, 2));
+}
+
+TrajectoryStateOnSurface* MuCSCTnPFlatTableProducer::matchTTwithCSCSeg(const reco::Track& track,
+ edm::Handle cscSegments,
+ CSCSegmentCollection::const_iterator& cscSegOut,
+ CSCDetId& idCSC) {
+ TrajectoryStateOnSurface* TrajSuf = nullptr;
+ Float_t deltaCSCR = 9999.;
+ for (CSCSegmentCollection::const_iterator segIt = cscSegments->begin(); segIt != cscSegments->end(); segIt++) {
+ CSCDetId id = (CSCDetId)(*segIt).cscDetId();
+
+ if (idCSC.endcap() != id.endcap())
+ continue;
+ if (idCSC.station() != id.station())
+ continue;
+ if (idCSC.chamber() != id.chamber())
+ continue;
+
+ Bool_t ed1 =
+ (idCSC.station() == 1) && ((idCSC.ring() == 1 || idCSC.ring() == 4) && (id.ring() == 1 || id.ring() == 4));
+ Bool_t ed2 =
+ (idCSC.station() == 1) && ((idCSC.ring() == 2 && id.ring() == 2) || (idCSC.ring() == 3 && id.ring() == 3));
+ Bool_t ed3 = (idCSC.station() != 1) && (idCSC.ring() == id.ring());
+ Bool_t TMCSCMatch = (ed1 || ed2 || ed3);
+
+ if (!TMCSCMatch)
+ continue;
+
+ const CSCChamber* cscchamber = m_cscGeometry->chamber(id);
+
+ if (!cscchamber)
+ continue;
+
+ TrajectoryStateOnSurface TrajSuf_ = surfExtrapTrkSam(track, cscchamber->toGlobal((*segIt).localPosition()).z());
+ Float_t dR_ = std::abs(TrajectoryDistToSeg(TrajSuf_, segIt));
+ if (dR_ < deltaCSCR) {
+ delete TrajSuf;
+ TrajSuf = new TrajectoryStateOnSurface(TrajSuf_);
+ deltaCSCR = dR_;
+ cscSegOut = segIt;
+ }
+ } //loop over segments
+
+ return TrajSuf;
+}
+
+std::vector MuCSCTnPFlatTableProducer::GetEdgeAndDistToGap(const reco::Track& track, CSCDetId& detid) {
+ std::vector result(4, 9999.);
+ result[3] = -9999;
+ const GeomDet* gdet = m_cscGeometry->idToDet(detid);
+ TrajectoryStateOnSurface tsos = surfExtrapTrkSam(track, gdet->surface().position().z());
+ if (!tsos.isValid())
+ return result;
+ LocalPoint localTTPos = gdet->surface().toLocal(tsos.freeState()->position());
+ const CSCWireTopology* wireTopology = m_cscGeometry->layer(detid)->geometry()->wireTopology();
+ Float_t wideWidth = wireTopology->wideWidthOfPlane();
+ Float_t narrowWidth = wireTopology->narrowWidthOfPlane();
+ Float_t length = wireTopology->lengthOfPlane();
+ // If slanted, there is no y offset between local origin and symmetry center of wire plane
+ Float_t yOfFirstWire = std::abs(wireTopology->wireAngle()) > 1.E-06 ? -0.5 * length : wireTopology->yOfWire(1);
+ // y offset between local origin and symmetry center of wire plane
+ Float_t yCOWPOffset = yOfFirstWire + 0.5 * length;
+ // tangent of the incline angle from inside the trapezoid
+ Float_t tangent = (wideWidth - narrowWidth) / (2. * length);
+ // y position wrt bottom of trapezoid
+ Float_t yPrime = localTTPos.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_t halfWidthAtYPrime = 0.5 * narrowWidth + yPrime * tangent;
+ // x offset between local origin and symmetry center of wire plane is zero
+ // x offset of ME11s is also zero. x center of wire groups is not at zero, because it is not parallel to x. The wire groups of ME11s have a complex geometry, see the code in m_debug.
+ Float_t edgex = std::abs(localTTPos.x()) - halfWidthAtYPrime;
+ Float_t edgey = std::abs(localTTPos.y() - yCOWPOffset) - 0.5 * length;
+ LocalError localTTErr = tsos.localError().positionError();
+ if (edgex > edgey) {
+ result[0] = edgex;
+ result[1] = sqrt(localTTErr.xx());
+ //result[1] = sqrt(tsos.cartesianError().position().cxx());
+ } else {
+ result[0] = edgey;
+ result[1] = sqrt(localTTErr.yy());
+ //result[1] = sqrt(tsos.cartesianError().position().cyy());
+ }
+ result[2] = YDistToHVDeadZone(localTTPos.y(), detid.station() * 10 + detid.ring());
+ result[3] = sqrt(localTTErr.yy());
+ return result; //return values: 1-edge;2-err of edge;3-disttogap;4-err of dist to gap
+}
+
+//deadzone center is according to http://cmssdt.cern.ch/SDT/lxr/source/RecoLocalMuon/CSCEfficiency/src/CSCEfficiency.cc#605
+//wire spacing is according to CSCTDR
+Float_t MuCSCTnPFlatTableProducer::YDistToHVDeadZone(Float_t yLocal, Int_t StationAndRing) {
+ //the ME11 wires are not parallel to x, but no gap
+ //chamber edges are not included.
+ const Float_t deadZoneCenterME1_2[2] = {-32.88305, 32.867423};
+ const Float_t deadZoneCenterME1_3[2] = {-22.7401, 27.86665};
+ const Float_t deadZoneCenterME2_1[2] = {-27.47, 33.67};
+ const Float_t deadZoneCenterME3_1[2] = {-36.21, 23.68};
+ const Float_t deadZoneCenterME4_1[2] = {-26.14, 23.85};
+ const Float_t deadZoneCenterME234_2[4] = {-81.8744, -21.18165, 39.51105, 100.2939};
+ const Float_t* deadZoneCenter;
+ Float_t deadZoneHeightHalf = 0.32 * 7 / 2; // wire spacing * (wires missing + 1)/2
+ Float_t minY = 999999.;
+ UChar_t nGaps = 2;
+ switch (std::abs(StationAndRing)) {
+ case 11:
+ case 14:
+ return 162; //the height of ME11
+ break;
+ case 12:
+ deadZoneCenter = deadZoneCenterME1_2;
+ break;
+ case 13:
+ deadZoneCenter = deadZoneCenterME1_3;
+ break;
+ case 21:
+ deadZoneCenter = deadZoneCenterME2_1;
+ break;
+ case 31:
+ deadZoneCenter = deadZoneCenterME3_1;
+ break;
+ case 41:
+ deadZoneCenter = deadZoneCenterME4_1;
+ break;
+ default:
+ deadZoneCenter = deadZoneCenterME234_2;
+ nGaps = 4;
+ }
+ for (UChar_t iGap = 0; iGap < nGaps; iGap++) {
+ Float_t newMinY = yLocal < deadZoneCenter[iGap] ? deadZoneCenter[iGap] - deadZoneHeightHalf - yLocal
+ : yLocal - (deadZoneCenter[iGap] + deadZoneHeightHalf);
+ if (newMinY < minY)
+ minY = newMinY;
+ }
+ return minY;
+}
+
+double MuCSCTnPFlatTableProducer::iso(const reco::Track& track, edm::Handle> tracks) {
+ double isoSum = 0.0;
+ for (const auto& track2 : (*tracks)) {
+ double dR = calcDeltaR(track.eta(), track2.eta(), track.phi(), track2.phi());
+ if (track2.pt() > 1.0 && dR > 0.001 && dR < 0.3)
+ isoSum += track2.pt();
+ }
+ return isoSum / track.pt();
+}
+
+double MuCSCTnPFlatTableProducer::calcDeltaR(double eta1, double eta2, double phi1, double phi2) {
+ double deta = eta1 - eta2;
+ if (phi1 < 0)
+ phi1 += 2.0 * M_PI;
+ if (phi2 < 0)
+ phi2 += 2.0 * M_PI;
+ double dphi = phi1 - phi2;
+ if (dphi > M_PI)
+ dphi -= 2. * M_PI;
+ else if (dphi < -M_PI)
+ dphi += 2. * M_PI;
+ return std::sqrt(deta * deta + dphi * dphi);
+}
+
+double MuCSCTnPFlatTableProducer::zMass(const reco::Track& track, const reco::Muon& muon) {
+ double zMass = -99.0;
+ double mMu = 0.1134289256;
+
+ zMass = std::pow((std::sqrt(std::pow(muon.p(), 2) + mMu * mMu) + std::sqrt(std::pow(track.p(), 2) + mMu * mMu)), 2) -
+ (std::pow((muon.px() + track.px()), 2) + std::pow((muon.py() + track.py()), 2) +
+ std::pow((muon.pz() + track.pz()), 2));
+
+ return std::sqrt(zMass);
+}
+
+#include "FWCore/PluginManager/interface/ModuleDef.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+DEFINE_FWK_MODULE(MuCSCTnPFlatTableProducer);
diff --git a/DPGAnalysis/MuonTools/plugins/MuDTMuonExtTableProducer.cc b/DPGAnalysis/MuonTools/plugins/MuDTMuonExtTableProducer.cc
new file mode 100644
index 0000000000000..a1653d692407c
--- /dev/null
+++ b/DPGAnalysis/MuonTools/plugins/MuDTMuonExtTableProducer.cc
@@ -0,0 +1,396 @@
+/** \class MuDTMuonExtTableProducer MuDTMuonExtTableProducer.cc DPGAnalysis/MuonTools/plugins/MuDTMuonExtTableProducer.cc
+ *
+ * Helper class : the muon filler
+ *
+ * \author L. Lunerti (INFN BO)
+ *
+ *
+ */
+
+#include "DataFormats/DetId/interface/DetId.h"
+#include "DataFormats/MuonDetId/interface/DTChamberId.h"
+#include "DataFormats/MuonReco/interface/MuonChamberMatch.h"
+#include "DataFormats/MuonReco/interface/MuonSegmentMatch.h"
+#include "DataFormats/MuonReco/interface/MuonSelectors.h"
+
+#include "DataFormats/DTRecHit/interface/DTRecSegment4D.h"
+#include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
+#include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
+#include "DataFormats/MuonDetId/interface/DTChamberId.h"
+
+#include "DataFormats/Math/interface/deltaR.h"
+
+#include "TString.h"
+#include "TRegexp.h"
+
+#include
+#include
+
+#include "DPGAnalysis/MuonTools/interface/MuBaseFlatTableProducer.h"
+
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+
+#include "DataFormats/Common/interface/View.h"
+#include "DataFormats/MuonReco/interface/Muon.h"
+#include "DataFormats/VertexReco/interface/Vertex.h"
+#include "DataFormats/MuonReco/interface/MuonFwd.h"
+#include "DataFormats/MuonReco/interface/MuonIsolation.h"
+#include "DataFormats/MuonReco/interface/MuonPFIsolation.h"
+#include "DataFormats/DTRecHit/interface/DTRecSegment4D.h"
+
+#include "DataFormats/Common/interface/TriggerResults.h"
+#include "DataFormats/HLTReco/interface/TriggerEvent.h"
+#include "DataFormats/HLTReco/interface/TriggerObject.h"
+#include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
+
+class MuDTMuonExtTableProducer : public MuBaseFlatTableProducer {
+public:
+ /// Constructor
+ MuDTMuonExtTableProducer(const edm::ParameterSet&);
+
+ /// Fill descriptors
+ static void fillDescriptions(edm::ConfigurationDescriptions&);
+
+protected:
+ /// Fill tree branches for a given events
+ void fillTable(edm::Event&) final;
+
+ /// Get info from the ES by run
+ void getFromES(const edm::Run&, const edm::EventSetup&) final;
+
+private:
+ /// Tokens
+ nano_mu::EDTokenHandle> m_muToken;
+ nano_mu::EDTokenHandle m_dtSegmentToken;
+
+ nano_mu::EDTokenHandle m_trigResultsToken;
+ nano_mu::EDTokenHandle m_trigEventToken;
+
+ /// Fill matches table
+ bool m_fillMatches;
+
+ /// Name of the triggers used by muon filler for trigger matching
+ std::string m_trigName;
+ std::string m_isoTrigName;
+
+ /// DT Geometry
+ nano_mu::ESTokenHandle m_dtGeometry;
+
+ /// HLT config provider
+ HLTConfigProvider m_hltConfig;
+
+ /// Indices of the triggers used by muon filler for trigger matching
+ std::vector m_trigIndices;
+ std::vector m_isoTrigIndices;
+
+ bool hasTrigger(std::vector&,
+ const trigger::TriggerObjectCollection&,
+ edm::Handle&,
+ const reco::Muon&);
+};
+
+MuDTMuonExtTableProducer::MuDTMuonExtTableProducer(const edm::ParameterSet& config)
+ : MuBaseFlatTableProducer(config),
+ m_muToken{config, consumesCollector(), "src"},
+ m_dtSegmentToken{config, consumesCollector(), "dtSegmentSrc"},
+ m_trigResultsToken{config, consumesCollector(), "trigResultsSrc"},
+ m_trigEventToken{config, consumesCollector(), "trigEventSrc"},
+ m_fillMatches{config.getParameter("fillMatches")},
+ m_trigName{config.getParameter("trigName")},
+ m_isoTrigName{config.getParameter("isoTrigName")},
+ m_dtGeometry{consumesCollector()} {
+ produces();
+ if (m_fillMatches) {
+ produces("matches");
+ produces("staMatches");
+ }
+}
+
+void MuDTMuonExtTableProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+ edm::ParameterSetDescription desc;
+
+ desc.add("name", "muon");
+ desc.add("src", edm::InputTag{"patMuons"});
+ desc.add("dtSegmentSrc", edm::InputTag{"dt4DSegments"});
+ desc.add("trigEventSrc", edm::InputTag{"hltTriggerSummaryAOD::HLT"});
+ desc.add("trigResultsSrc", edm::InputTag{"TriggerResults::HLT"});
+
+ desc.add("fillMatches", true);
+
+ desc.add("trigName", "HLT_Mu55*");
+ desc.add("isoTrigName", "HLT_IsoMu2*");
+
+ descriptions.addWithDefaultLabel(desc);
+}
+
+void MuDTMuonExtTableProducer::getFromES(const edm::Run& run, const edm::EventSetup& environment) {
+ m_dtGeometry.getFromES(environment);
+
+ bool changed{true};
+ m_hltConfig.init(run, environment, "HLT", changed);
+
+ const bool enableWildcard{true};
+
+ TString tName = TString(m_trigName);
+ TRegexp tNamePattern = TRegexp(tName, enableWildcard);
+
+ for (unsigned iPath = 0; iPath < m_hltConfig.size(); ++iPath) {
+ TString pathName = TString(m_hltConfig.triggerName(iPath));
+ if (pathName.Contains(tNamePattern))
+ m_trigIndices.push_back(static_cast(iPath));
+ }
+
+ tName = TString(m_isoTrigName);
+ tNamePattern = TRegexp(tName, enableWildcard);
+
+ for (unsigned iPath = 0; iPath < m_hltConfig.size(); ++iPath) {
+ TString pathName = TString(m_hltConfig.triggerName(iPath));
+ if (pathName.Contains(tNamePattern))
+ m_isoTrigIndices.push_back(static_cast(iPath));
+ }
+}
+
+void MuDTMuonExtTableProducer::fillTable(edm::Event& ev) {
+ unsigned int nMuons{0};
+
+ std::vector firesIsoTrig;
+ std::vector firesTrig;
+
+ std::vector nMatches;
+ std::vector staMu_nMatchSeg;
+
+ std::vector matches_begin;
+ std::vector matches_end;
+
+ std::vector staMatches_begin;
+ std::vector staMatches_end;
+
+ std::vector matches_wheel;
+ std::vector