Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:36:00 

0001 /** \class MuCSCTnPFlatTableProducer MuCSCTnPFlatTableProducer.cc DPGAnalysis/MuonTools/plugins/MuCSCTnPFlatTableProducer.cc
0002  *  
0003  * Helper class : the CSC Tag and probe segment efficiency  filler
0004  *
0005  * \author M. Herndon (UW Madison)
0006  *
0007  *
0008  */
0009 
0010 #include "DataFormats/DetId/interface/DetId.h"
0011 #include "DataFormats/CSCRecHit/interface/CSCSegment.h"
0012 #include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
0013 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0014 
0015 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0016 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0017 
0018 #include "DataFormats/Math/interface/deltaR.h"
0019 
0020 #include "TString.h"
0021 #include "TRegexp.h"
0022 
0023 #include <iostream>
0024 
0025 #include <numeric>
0026 #include <vector>
0027 
0028 #include "DPGAnalysis/MuonTools/interface/MuBaseFlatTableProducer.h"
0029 
0030 #include "FWCore/Framework/interface/ESHandle.h"
0031 
0032 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0033 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0034 
0035 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0036 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0037 
0038 #include "DataFormats/MuonReco/interface/Muon.h"
0039 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0040 #include "DataFormats/MuonReco/interface/MuonIsolation.h"
0041 #include "DataFormats/MuonReco/interface/MuonPFIsolation.h"
0042 
0043 #include "DataFormats/TrackReco/interface/Track.h"
0044 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0045 
0046 #include "DataFormats/VertexReco/interface/Vertex.h"
0047 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0048 
0049 #include "DataFormats/Common/interface/TriggerResults.h"
0050 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0051 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0052 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0053 
0054 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0055 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0056 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0057 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0058 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0059 
0060 class MuonServiceProxy;
0061 
0062 class MuCSCTnPFlatTableProducer : public MuBaseFlatTableProducer {
0063 public:
0064   /// Constructor
0065   MuCSCTnPFlatTableProducer(const edm::ParameterSet&);
0066 
0067   /// Fill descriptors
0068   static void fillDescriptions(edm::ConfigurationDescriptions&);
0069 
0070 protected:
0071   /// Fill tree branches for a given events
0072   void fillTable(edm::Event&) final;
0073 
0074   /// Get info from the ES by run
0075   void getFromES(const edm::Run&, const edm::EventSetup&) final;
0076 
0077   /// Get info from the ES for a given event
0078   void getFromES(const edm::EventSetup&) final;
0079 
0080 private:
0081   static constexpr Float_t MEZ[6] = {601.3, 696.11, 696.11, 827.56, 936.44, 1025.9};
0082 
0083   /// Tokens
0084   nano_mu::EDTokenHandle<reco::MuonCollection> m_muToken;
0085   nano_mu::EDTokenHandle<reco::TrackCollection> m_trackToken;
0086 
0087   nano_mu::EDTokenHandle<CSCSegmentCollection> m_cscSegmentToken;
0088 
0089   nano_mu::EDTokenHandle<std::vector<reco::Vertex>> m_primaryVerticesToken;
0090 
0091   nano_mu::EDTokenHandle<edm::TriggerResults> m_trigResultsToken;
0092   nano_mu::EDTokenHandle<trigger::TriggerEvent> m_trigEventToken;
0093 
0094   /// Name of the triggers used by muon filler for trigger matching
0095   std::string m_trigName;
0096   std::string m_isoTrigName;
0097 
0098   /// Handles to geometry, detector and specialized objects
0099   /// CSC Geometry
0100   nano_mu::ESTokenHandle<CSCGeometry, MuonGeometryRecord, edm::Transition::BeginRun> m_cscGeometry;
0101 
0102   /// Muon service proxy
0103   std::unique_ptr<MuonServiceProxy> m_muonSP;
0104 
0105   /// Transient Track Builder
0106   nano_mu::ESTokenHandle<TransientTrackBuilder, TransientTrackRecord> m_transientTrackBuilder;
0107 
0108   // Extrapolator to cylinder
0109   edm::ESHandle<Propagator> propagatorAlong;
0110   edm::ESHandle<Propagator> propagatorOpposite;
0111   edm::ESHandle<MagneticField> theBField;
0112 
0113   /// HLT config provider
0114   HLTConfigProvider m_hltConfig;
0115 
0116   /// Indices of the triggers used by muon filler for trigger matching
0117   std::vector<int> m_trigIndices;
0118   std::vector<int> m_isoTrigIndices;
0119 
0120   /// Selection functions
0121   //bool muonTagSelection(const reco::Muon & muon,edm::Handle<std::vector<reco::Track>> tracks);
0122   bool trackProbeSelection(const reco::Track& track, edm::Handle<std::vector<reco::Track>>);
0123   bool muonTagSelection(const reco::Muon&);
0124   //bool trackProbeSelection(const reco::Track & track);
0125   bool zSelection(const reco::Muon&, const reco::Track&);
0126 
0127   // Calculation functions
0128   double zMass(const reco::Track&, const reco::Muon&);
0129   double calcDeltaR(double, double, double, double);
0130   double iso(const reco::Track&, edm::Handle<std::vector<reco::Track>>);
0131 
0132   // Track extrapolation and segment match functions
0133   TrajectoryStateOnSurface surfExtrapTrkSam(const reco::Track&, double);
0134   FreeTrajectoryState freeTrajStateMuon(const reco::Track&);
0135 
0136   UChar_t ringCandidate(Int_t iiStation, Int_t station, Float_t feta, Float_t phi);
0137   UChar_t thisChamberCandidate(UChar_t station, UChar_t ring, Float_t phi);
0138 
0139   TrajectoryStateOnSurface* matchTTwithCSCSeg(const reco::Track&,
0140                                               edm::Handle<CSCSegmentCollection>,
0141                                               CSCSegmentCollection::const_iterator&,
0142                                               CSCDetId&);
0143   Float_t TrajectoryDistToSeg(TrajectoryStateOnSurface, CSCSegmentCollection::const_iterator);
0144   std::vector<Float_t> GetEdgeAndDistToGap(const reco::Track&, CSCDetId&);
0145   Float_t YDistToHVDeadZone(Float_t, Int_t);
0146 
0147   /// The variables holding
0148   /// the T&P related information
0149 
0150   unsigned int m_nZCands;  // the # of digis (size of all following vectors)
0151 
0152   double _trackIso;
0153   double _muonIso;
0154   double _zMass;
0155 
0156   bool hasTrigger(std::vector<int>&,
0157                   const trigger::TriggerObjectCollection&,
0158                   edm::Handle<trigger::TriggerEvent>&,
0159                   const reco::Muon&);
0160 
0161   float computeTrkIso(const reco::MuonIsolation&, float);
0162   float computePFIso(const reco::MuonPFIsolation&, float);
0163 };
0164 
0165 MuCSCTnPFlatTableProducer::MuCSCTnPFlatTableProducer(const edm::ParameterSet& config)
0166     : MuBaseFlatTableProducer(config),
0167       m_muToken{config, consumesCollector(), "muonSrc"},
0168       m_trackToken{config, consumesCollector(), "trackSrc"},
0169       m_cscSegmentToken{config, consumesCollector(), "cscSegmentSrc"},
0170       m_primaryVerticesToken{config, consumesCollector(), "primaryVerticesSrc"},
0171       m_trigResultsToken{config, consumesCollector(), "trigResultsSrc"},
0172       m_trigEventToken{config, consumesCollector(), "trigEventSrc"},
0173       m_trigName{config.getParameter<std::string>("trigName")},
0174       m_isoTrigName{config.getParameter<std::string>("isoTrigName")},
0175       m_cscGeometry{consumesCollector()},
0176       m_muonSP{std::make_unique<MuonServiceProxy>(config.getParameter<edm::ParameterSet>("ServiceParameters"),
0177                                                   consumesCollector())},
0178       m_transientTrackBuilder{consumesCollector(), "TransientTrackBuilder"} {
0179   produces<nanoaod::FlatTable>();
0180 }
0181 
0182 void MuCSCTnPFlatTableProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0183   edm::ParameterSetDescription desc;
0184 
0185   desc.add<std::string>("name", "cscTnP");
0186   desc.add<edm::InputTag>("muonSrc", edm::InputTag{"muons"});
0187   desc.add<edm::InputTag>("trackSrc", edm::InputTag{"generalTracks"});
0188   desc.add<edm::InputTag>("cscSegmentSrc", edm::InputTag{"cscSegments"});
0189   desc.add<edm::InputTag>("primaryVerticesSrc", edm::InputTag{"offlinePrimaryVertices"});
0190 
0191   desc.add<edm::InputTag>("trigEventSrc", edm::InputTag{"hltTriggerSummaryAOD::HLT"});
0192   desc.add<edm::InputTag>("trigResultsSrc", edm::InputTag{"TriggerResults::HLT"});
0193 
0194   desc.add<std::string>("trigName", "none");
0195   desc.add<std::string>("isoTrigName", "HLT_IsoMu2*");
0196 
0197   desc.setAllowAnything();
0198 
0199   descriptions.addWithDefaultLabel(desc);
0200 }
0201 
0202 void MuCSCTnPFlatTableProducer::getFromES(const edm::Run& run, const edm::EventSetup& environment) {
0203   m_cscGeometry.getFromES(environment);
0204 
0205   bool changed{true};
0206   m_hltConfig.init(run, environment, "HLT", changed);
0207 
0208   const bool enableWildcard{true};
0209 
0210   TString tName = TString(m_trigName);
0211   TRegexp tNamePattern = TRegexp(tName, enableWildcard);
0212 
0213   for (unsigned iPath = 0; iPath < m_hltConfig.size(); ++iPath) {
0214     TString pathName = TString(m_hltConfig.triggerName(iPath));
0215     if (pathName.Contains(tNamePattern))
0216       m_trigIndices.push_back(static_cast<int>(iPath));
0217   }
0218 
0219   tName = TString(m_isoTrigName);
0220   tNamePattern = TRegexp(tName, enableWildcard);
0221 
0222   for (unsigned iPath = 0; iPath < m_hltConfig.size(); ++iPath) {
0223     TString pathName = TString(m_hltConfig.triggerName(iPath));
0224     if (pathName.Contains(tNamePattern))
0225       m_isoTrigIndices.push_back(static_cast<int>(iPath));
0226   }
0227 }
0228 
0229 void MuCSCTnPFlatTableProducer::getFromES(const edm::EventSetup& environment) {
0230   m_transientTrackBuilder.getFromES(environment);
0231   m_muonSP->update(environment);
0232 }
0233 
0234 void MuCSCTnPFlatTableProducer::fillTable(edm::Event& ev) {
0235   unsigned int m_nZCands = 0;  // the # of digis (size of all following vectors)
0236 
0237   // Muon track tag variables
0238   std::vector<float> m_muonPt;        // muon pT [GeV/c]
0239   std::vector<float> m_muonPhi;       // muon phi [rad]
0240   std::vector<float> m_muonEta;       // muon eta
0241   std::vector<float> m_muonPtError;   // muon pT [GeV/c] error
0242   std::vector<float> m_muonPhiError;  // muon phi [rad] error
0243   std::vector<float> m_muonEtaError;  // muon eta error
0244   std::vector<int> m_muonCharge;      // muon charge
0245   std::vector<float> m_muonDXY;       // muon dXY
0246   std::vector<float> m_muonDZ;        // muon dZ
0247   std::vector<int> m_muonTrkHits;     // muon track Hits
0248   std::vector<float> m_muonChi2;      // muon Chi2
0249   std::vector<bool> m_muonTrigger;    // muon trigger
0250   std::vector<float> m_muonIso;       // track Iso
0251 
0252   // Track probe variabes
0253   std::vector<float> m_trackPt;        // track pT [GeV/c]
0254   std::vector<float> m_trackP;         // track P [GeV/c]
0255   std::vector<float> m_trackPhi;       // track phi [rad]
0256   std::vector<float> m_trackEta;       // track eta
0257   std::vector<float> m_trackPtError;   // track pT [GeV/c] error
0258   std::vector<float> m_trackPhiError;  // track phi [rad] error
0259   std::vector<float> m_trackEtaError;  // track eta error
0260   std::vector<int> m_trackCharge;      // track charge
0261   std::vector<float> m_trackDXY;       // track dXY
0262   std::vector<float> m_trackDZ;        // track dZ
0263   std::vector<int> m_trackTrkHits;     // track Hits
0264   std::vector<float> m_trackChi2;      // track Chi2
0265   std::vector<float> m_trackIso;       // track Iso
0266 
0267   // Z and global variables
0268   std::vector<float> m_zMass;                    // z mass
0269   std::vector<float> m_dRTrackMuon;              // dR between the track and muon
0270   std::vector<float> m_numberOfPrimaryVertices;  // Number of primary Vertices
0271 
0272   // CSC chamber information, station encoded in vector
0273   std::vector<int> m_chamberEndcap;                  // chamber endcap
0274                                                      // station encoded in array index
0275   std::array<std::vector<int>, 4> m_chamberRing;     // chamber Ring
0276   std::array<std::vector<int>, 4> m_chamberChamber;  // chamber Chamber
0277   std::array<std::vector<int>, 4> m_chamberLayer;    // Segment layer information
0278 
0279   // Track intersection variables
0280   std::array<std::vector<float>, 4> m_ttIntLocalX;       // track trajector intersection local X on stations 1-4
0281   std::array<std::vector<float>, 4> m_ttIntLocalY;       // track trajector intersection local Y on stations 1-4
0282   std::array<std::vector<float>, 4> m_ttIntLocalErrorX;  // track trajector intersection local X on stations 1-4
0283   std::array<std::vector<float>, 4> m_ttIntLocalErrorY;  // track trajector intersection local Y on stations 1-4
0284   std::array<std::vector<float>, 4> m_ttIntLocalW;       // track trajector intersection local Wire on stations 1-4
0285   std::array<std::vector<float>, 4> m_ttIntLocalS;       // track trajector intersection local Strip on stations 1-4
0286   std::array<std::vector<float>, 4> m_ttIntEta;          // track trajector intersection Eta stations 1-4
0287 
0288   // Track intersection fiducial information
0289 
0290   std::array<std::vector<float>, 4>
0291       m_ttDistToEdge;  // track trajector intersection distance to edge, neg is with chamber, on stations 1-4
0292   std::array<std::vector<float>, 4> m_ttDistToHVGap;  // track trajector intersection distance to HV GAP on stations 1-4
0293 
0294   // Segment location variables
0295   std::array<std::vector<float>, 4> m_segLocalX;       // segment local X on stations 1-4
0296   std::array<std::vector<float>, 4> m_segLocalY;       // segment local Y on stations 1-4
0297   std::array<std::vector<float>, 4> m_segLocalErrorX;  // segment local X error on stations 1-4
0298   std::array<std::vector<float>, 4> m_segLocalErrorY;  // segment local Y error on stations 1-4
0299 
0300   // track intersection segment residuals variables
0301   std::array<std::vector<float>, 4>
0302       m_ttIntSegResidualLocalX;  // track trajector intersection  Segment residual local X on stations 1-4
0303   std::array<std::vector<float>, 4>
0304       m_ttIntSegResidualLocalY;  // track trajector intersection  Segment residuallocal Y on stations 1-4
0305 
0306   auto&& propagator_along = m_muonSP->propagator("SteppingHelixPropagatorAlong");
0307   auto&& propagator_opposite = m_muonSP->propagator("SteppingHelixPropagatorOpposite");
0308 
0309   propagatorAlong = propagator_along;
0310   propagatorOpposite = propagator_opposite;
0311 
0312   theBField = m_muonSP->magneticField();
0313 
0314   auto muons = m_muToken.conditionalGet(ev);
0315   auto tracks = m_trackToken.conditionalGet(ev);
0316   auto segments = m_cscSegmentToken.conditionalGet(ev);
0317   auto primaryVertices = m_primaryVerticesToken.conditionalGet(ev);
0318 
0319   auto triggerResults = m_trigResultsToken.conditionalGet(ev);
0320   auto triggerEvent = m_trigEventToken.conditionalGet(ev);
0321 
0322   if (muons.isValid() && tracks.isValid() && segments.isValid() && primaryVertices.isValid() &&
0323       m_transientTrackBuilder.isValid()) {
0324     for (const auto& muon : (*muons)) {
0325       if (!muonTagSelection(muon))
0326         continue;
0327 
0328       bool muonTrigger = false;
0329       if (triggerResults.isValid() && triggerEvent.isValid()) {
0330         const auto& triggerObjects = triggerEvent->getObjects();
0331         muonTrigger = (hasTrigger(m_isoTrigIndices, triggerObjects, triggerEvent, muon) ||
0332                        hasTrigger(m_trigIndices, triggerObjects, triggerEvent, muon));
0333       }
0334 
0335       for (const auto& track : (*tracks)) {
0336         if (!trackProbeSelection(track, tracks))
0337           continue;
0338         if (!zSelection(muon, track))
0339           continue;
0340         //std::cout << "Z candidate found: " << _zMass << " track eta: " << track.eta() << std::endl;
0341         //std::cout.flush();
0342         m_nZCands++;
0343 
0344         m_trackPt.push_back(track.pt());
0345         m_trackP.push_back(track.p());
0346         m_trackEta.push_back(track.eta());
0347         m_trackPhi.push_back(track.phi());
0348         m_trackPtError.push_back(track.pt());
0349         m_trackEtaError.push_back(track.eta());
0350         m_trackPhiError.push_back(track.phi());
0351         m_trackCharge.push_back(track.charge());
0352         m_trackDXY.push_back(track.dxy());
0353         m_trackDZ.push_back(track.dz());
0354         m_trackTrkHits.push_back(track.hitPattern().numberOfValidTrackerHits());
0355         m_trackChi2.push_back(track.normalizedChi2());
0356         m_trackIso.push_back(_trackIso);
0357 
0358         m_muonPt.push_back(muon.track()->pt());
0359         m_muonPhi.push_back(muon.track()->phi());
0360         m_muonEta.push_back(muon.track()->eta());
0361         m_muonPtError.push_back(muon.track()->ptError());
0362         m_muonPhiError.push_back(muon.track()->phiError());
0363         m_muonEtaError.push_back(muon.track()->etaError());
0364         m_muonCharge.push_back(muon.charge());
0365         m_muonDXY.push_back(muon.track()->dxy());
0366         m_muonDZ.push_back(muon.track()->dz());
0367         m_muonTrkHits.push_back(muon.track()->hitPattern().numberOfValidTrackerHits());
0368         m_muonChi2.push_back(muon.track()->normalizedChi2());
0369         m_muonIso.push_back(computeTrkIso(muon.isolationR03(), muon.pt()));
0370         m_muonTrigger.push_back(muonTrigger);
0371 
0372         m_zMass.push_back(_zMass);
0373         double_t dR = calcDeltaR(track.eta(), muon.eta(), track.phi(), muon.phi());
0374         //double_t dR = 1.0;
0375         m_dRTrackMuon.push_back(dR);
0376         const reco::VertexCollection& vertices = *primaryVertices.product();
0377         m_numberOfPrimaryVertices.push_back(vertices.size());
0378 
0379         bool ec = (track.eta() > 0);
0380         UChar_t endcapCSC = ec ? 0 : 1;
0381         m_chamberEndcap.push_back(endcapCSC * 1);
0382 
0383         Int_t iiStationFail = 0;
0384         Int_t iiStation0Pass = 0;
0385         for (int iiStationZ = 0; iiStationZ < 6; iiStationZ++) {
0386           UChar_t stationCSC = iiStationZ > 2 ? iiStationZ - 2 : 0;
0387           UChar_t ringCSC = 0;
0388           // Could monitor this problem here if necessary;
0389           if (stationCSC == 0 && iiStation0Pass > 0)
0390             continue;
0391           TrajectoryStateOnSurface tsos = surfExtrapTrkSam(track, ec ? MEZ[iiStationZ] : -MEZ[iiStationZ]);
0392 
0393           if (tsos.isValid()) {
0394             Float_t trkEta = tsos.globalPosition().eta(), trkPhi = tsos.globalPosition().phi();
0395             ringCSC = ringCandidate(iiStationZ, stationCSC + 1, trkEta, trkPhi);
0396 
0397             if (ringCSC) {
0398               UChar_t chamberCSC = thisChamberCandidate(stationCSC + 1, ringCSC, track.phi()) - 1;
0399               CSCDetId Layer3id = CSCDetId(endcapCSC + 1, stationCSC + 1, ringCSC, chamberCSC + 1, 3);
0400               CSCDetId Layer0Id = CSCDetId(endcapCSC + 1,
0401                                            stationCSC + 1,
0402                                            ringCSC,
0403                                            chamberCSC + 1,
0404                                            0);  //layer 0 is the mid point of the chamber. It is not a real layer.
0405               // !!!!! need to fix Layer0Id problem with ME1/1 here
0406 
0407               const BoundPlane& Layer3Surface = m_cscGeometry->idToDet(Layer3id)->surface();
0408 
0409               tsos = surfExtrapTrkSam(track, Layer3Surface.position().z());
0410 
0411               if (tsos.isValid()) {
0412                 if (stationCSC == 0)
0413                   iiStation0Pass++;
0414                 // Fill track intersection denominator information
0415                 LocalPoint localTTIntPoint = Layer3Surface.toLocal(tsos.freeState()->position());
0416                 const CSCLayerGeometry* layerGeoma = m_cscGeometry->chamber(Layer0Id)->layer(3)->geometry();
0417                 const CSCLayerGeometry* layerGeomb = m_cscGeometry->chamber(Layer0Id)->layer(4)->geometry();
0418 
0419                 m_chamberRing[stationCSC].push_back(ringCSC);
0420                 m_chamberChamber[stationCSC].push_back(chamberCSC);
0421                 m_ttIntLocalX[stationCSC].push_back(localTTIntPoint.x());
0422                 m_ttIntLocalY[stationCSC].push_back(localTTIntPoint.y());
0423                 m_ttIntLocalW[stationCSC].push_back(
0424                     (layerGeoma->nearestWire(localTTIntPoint) + layerGeomb->nearestWire(localTTIntPoint)) / 2.0);
0425                 m_ttIntLocalS[stationCSC].push_back(
0426                     (layerGeoma->strip(localTTIntPoint) + layerGeomb->strip(localTTIntPoint)) / 2.0);
0427                 m_ttIntEta[stationCSC].push_back(trkEta);
0428 
0429                 // Errors are those of the track intersection, chosing the plane and exact geomentry is performed in the function
0430                 Float_t CSCProjEdgeDist = -9999.0;
0431                 Float_t ttIntLocalErrorX = -9999.0;
0432                 Float_t CSCDyProjHVGap = 9999.0;
0433                 Float_t ttIntLocalErrorY = -9999.0;
0434                 for (Int_t ly = 1; ly < 7; ly++) {
0435                   CSCDetId Layerid = CSCDetId(endcapCSC + 1, stationCSC + 1, ringCSC, chamberCSC + 1, ly);
0436                   std::vector<Float_t> EdgeAndDistToGap(GetEdgeAndDistToGap(
0437                       track, Layerid));  //values: 1-edge;2-err of edge;3-disttogap;4-err of dist to gap
0438                   if (EdgeAndDistToGap[0] > CSCProjEdgeDist) {
0439                     CSCProjEdgeDist = EdgeAndDistToGap[0];
0440                     ttIntLocalErrorX = EdgeAndDistToGap[1];
0441                   }
0442                   if (EdgeAndDistToGap[2] < CSCDyProjHVGap) {
0443                     CSCDyProjHVGap = EdgeAndDistToGap[2];
0444                     ttIntLocalErrorY = EdgeAndDistToGap[3];
0445                   }
0446                 }
0447                 m_ttDistToEdge[stationCSC].push_back(CSCProjEdgeDist);
0448                 m_ttDistToHVGap[stationCSC].push_back(CSCDyProjHVGap);
0449                 m_ttIntLocalErrorX[stationCSC].push_back(ttIntLocalErrorX);
0450                 m_ttIntLocalErrorY[stationCSC].push_back(ttIntLocalErrorY);
0451 
0452                 // now we have a local point for comparison to segments
0453                 CSCSegmentCollection::const_iterator cscSegOut;
0454                 TrajectoryStateOnSurface* TrajToSeg = matchTTwithCSCSeg(track, segments, cscSegOut, Layer3id);
0455 
0456                 if (TrajToSeg == nullptr) {
0457                   // fill Null Num
0458                   m_segLocalX[stationCSC].push_back(-9999.0);
0459                   m_segLocalY[stationCSC].push_back(-9999.0);
0460                   m_segLocalErrorX[stationCSC].push_back(0.0);
0461                   m_segLocalErrorY[stationCSC].push_back(0.0);
0462 
0463                   m_ttIntSegResidualLocalX[stationCSC].push_back(-9990.0);
0464                   m_ttIntSegResidualLocalY[stationCSC].push_back(-9990.0);
0465 
0466                   m_chamberLayer[stationCSC].push_back(-9);
0467 
0468                   continue;
0469                 }
0470 
0471                 LocalPoint localSegmentPoint = (*cscSegOut).localPosition();
0472                 LocalError localSegErr = (*cscSegOut).localPositionError();
0473 
0474                 m_segLocalX[stationCSC].push_back(localSegmentPoint.x());
0475                 m_segLocalY[stationCSC].push_back(localSegmentPoint.y());
0476                 m_segLocalErrorX[stationCSC].push_back(sqrt(localSegErr.xx()));
0477                 m_segLocalErrorY[stationCSC].push_back(sqrt(localSegErr.yy()));
0478 
0479                 m_ttIntSegResidualLocalX[stationCSC].push_back(localTTIntPoint.x() - localSegmentPoint.x());
0480                 m_ttIntSegResidualLocalY[stationCSC].push_back(localTTIntPoint.y() - localSegmentPoint.y());
0481                 /* Extract layers for hits */
0482                 int layers = 0;
0483                 for (std::vector<CSCRecHit2D>::const_iterator itRH = cscSegOut->specificRecHits().begin();
0484                      itRH != cscSegOut->specificRecHits().end();
0485                      ++itRH) {
0486                   const CSCRecHit2D* recHit = &(*itRH);
0487                   int layer = recHit->cscDetId().layer();
0488                   layers |= 1 << (layer - 1);
0489                 }
0490                 m_chamberLayer[stationCSC].push_back(layers);
0491 
0492               }  // end preliminary tsos is valid
0493 
0494             }  // end found ring CSC
0495 
0496           }  // end refined tsos is valid
0497 
0498           if ((!tsos.isValid()) || (ringCSC == 0)) {
0499             // only fill Null denominator once for station 1, iiStation Z = 0,1,2
0500             if (iiStationZ <= 2)
0501               iiStationFail++;
0502             if (iiStationZ > 2 || iiStationFail == 3) {
0503               // fill Null Den Num
0504               m_chamberRing[stationCSC].push_back(-9);
0505               m_chamberChamber[stationCSC].push_back(-9);
0506               m_ttIntLocalX[stationCSC].push_back(-9999.0);
0507               m_ttIntLocalY[stationCSC].push_back(-9999.0);
0508               m_ttIntLocalErrorX[stationCSC].push_back(0.0);
0509               m_ttIntLocalErrorY[stationCSC].push_back(0.0);
0510               m_ttIntLocalW[stationCSC].push_back(-9999.0);
0511               m_ttIntLocalS[stationCSC].push_back(-9999.0);
0512               m_ttIntEta[stationCSC].push_back(-9999.0);
0513 
0514               m_ttDistToEdge[stationCSC].push_back(-9999.0);
0515               m_ttDistToHVGap[stationCSC].push_back(-9999.9);
0516 
0517               m_segLocalX[stationCSC].push_back(-9999.0);
0518               m_segLocalY[stationCSC].push_back(-9999.0);
0519               m_segLocalErrorX[stationCSC].push_back(0.0);
0520               m_segLocalErrorY[stationCSC].push_back(0.0);
0521 
0522               m_ttIntSegResidualLocalX[stationCSC].push_back(-9990.0);
0523               m_ttIntSegResidualLocalY[stationCSC].push_back(-9990.0);
0524 
0525               m_chamberLayer[stationCSC].push_back(-9);
0526             }
0527           }
0528 
0529         }  // end loop over CSC Z planes
0530       }  // endl loop over tracks
0531     }  // end loop over muons
0532 
0533   }  // End if good physics objects
0534 
0535   //  if (m_nZCands>0) {
0536   auto table = std::make_unique<nanoaod::FlatTable>(m_nZCands, m_name, false, false);
0537 
0538   table->setDoc("CSC Tag & Probe segment efficiency  information");
0539 
0540   addColumn(table, "muonPt", m_muonPt, "muon pt [GeV/c]");
0541   addColumn(table, "muonPhi", m_muonPhi, "muon phi [rad]");
0542   addColumn(table, "muonEta", m_muonEta, "muon eta");
0543   addColumn(table, "muonPtError", m_muonPtError, "muon pt error [GeV/c]");
0544   addColumn(table, "muonPhiError", m_muonPhiError, "muon phi error [rad]");
0545   addColumn(table, "muonEtaError", m_muonEtaError, "muon eta error");
0546   addColumn(table, "muonCharge", m_muonCharge, "muon charge");
0547   addColumn(table, "muonDXY", m_muonDXY, "muon dXY [cm]");
0548   addColumn(table, "muonDZ", m_muonDZ, "muon dZ [cm]");
0549   addColumn(table, "muonTrkHits", m_muonTrkHits, "muon track hits");
0550   addColumn(table, "muonChi2", m_muonChi2, "muon chi2");
0551   addColumn(table, "muonIso", m_trackIso, "muon relative iso");
0552   addColumn(table, "muonTrigger", m_muonTrigger, "muon has trigger bool");
0553 
0554   addColumn(table, "trackPt", m_trackPt, "track pt [GeV/c]");
0555   addColumn(table, "trackP", m_trackPt, "track p [GeV/c]");
0556   addColumn(table, "trackPhi", m_trackPhi, "track phi [rad]");
0557   addColumn(table, "trackEta", m_trackEta, "track eta");
0558   addColumn(table, "trackPtError", m_trackPtError, "track pt error [GeV/c]");
0559   addColumn(table, "trackPhiError", m_trackPhiError, "track phi error [rad]");
0560   addColumn(table, "trackEtaError", m_trackEtaError, "track eta error");
0561   addColumn(table, "trackCharge", m_trackCharge, "track charge");
0562   addColumn(table, "trackDXY", m_trackDXY, "track dXY [cm]");
0563   addColumn(table, "trackDZ", m_trackDZ, "track dZ [cm]");
0564   addColumn(table, "trackTrkHits", m_trackTrkHits, "track track hits");
0565   addColumn(table, "trackChi2", m_trackChi2, "track chi2");
0566   addColumn(table, "trackIso", m_trackIso, "track relative iso");
0567 
0568   addColumn(table, "zMass", m_zMass, "Z mass [GeV/c^2]");
0569   addColumn(table, "dRTrackMuon", m_dRTrackMuon, "dR between track and muon");
0570   addColumn(table, "numberOfPrimaryVertidies", m_numberOfPrimaryVertices, "Number of PVs");
0571 
0572   addColumn(table, "chamberEndcap", m_chamberEndcap, "");
0573   addColumn(table, "chamberRing1", m_chamberRing[0], "");
0574   addColumn(table, "chamberRing2", m_chamberRing[1], "");
0575   addColumn(table, "chamberRing3", m_chamberRing[2], "");
0576   addColumn(table, "chamberRing4", m_chamberRing[3], "");
0577   addColumn(table, "chamberChamber1", m_chamberChamber[0], "");
0578   addColumn(table, "chamberChamber2", m_chamberChamber[1], "");
0579   addColumn(table, "chamberChamber3", m_chamberChamber[2], "");
0580   addColumn(table, "chamberChamber4", m_chamberChamber[3], "");
0581   addColumn(table, "chamberLayer1", m_chamberLayer[0], "");
0582   addColumn(table, "chamberLayer2", m_chamberLayer[1], "");
0583   addColumn(table, "chamberLayer3", m_chamberLayer[2], "");
0584   addColumn(table, "chamberLayer4", m_chamberLayer[3], "");
0585 
0586   addColumn(table, "ttIntLocalX1", m_ttIntLocalX[0], "");
0587   addColumn(table, "ttIntLocalX2", m_ttIntLocalX[1], "");
0588   addColumn(table, "ttIntLocalX3", m_ttIntLocalX[2], "");
0589   addColumn(table, "ttIntLocalX4", m_ttIntLocalX[3], "");
0590   addColumn(table, "ttIntLocalY1", m_ttIntLocalY[0], "");
0591   addColumn(table, "ttIntLocalY2", m_ttIntLocalY[1], "");
0592   addColumn(table, "ttIntLocalY3", m_ttIntLocalY[2], "");
0593   addColumn(table, "ttIntLocalY4", m_ttIntLocalY[3], "");
0594   addColumn(table, "ttIntLocalErrorX1", m_ttIntLocalErrorX[0], "");
0595   addColumn(table, "ttIntLocalErrorX2", m_ttIntLocalErrorX[1], "");
0596   addColumn(table, "ttIntLocalErrorX3", m_ttIntLocalErrorX[2], "");
0597   addColumn(table, "ttIntLocalErrorX4", m_ttIntLocalErrorX[3], "");
0598   addColumn(table, "ttIntLocalErrorY1", m_ttIntLocalErrorY[0], "");
0599   addColumn(table, "ttIntLocalErrorY2", m_ttIntLocalErrorY[1], "");
0600   addColumn(table, "ttIntLocalErrorY3", m_ttIntLocalErrorY[2], "");
0601   addColumn(table, "ttIntLocalErrorY4", m_ttIntLocalErrorY[3], "");
0602   addColumn(table, "ttIntLocalW1", m_ttIntLocalW[0], "");
0603   addColumn(table, "ttIntLocalW2", m_ttIntLocalW[1], "");
0604   addColumn(table, "ttIntLocalW3", m_ttIntLocalW[2], "");
0605   addColumn(table, "ttIntLocalW4", m_ttIntLocalW[3], "");
0606   addColumn(table, "ttIntLocalS1", m_ttIntLocalS[0], "");
0607   addColumn(table, "ttIntLocalS2", m_ttIntLocalS[1], "");
0608   addColumn(table, "ttIntLocalS3", m_ttIntLocalS[2], "");
0609   addColumn(table, "ttIntLocalS4", m_ttIntLocalS[3], "");
0610   addColumn(table, "ttIntEta1", m_ttIntEta[0], "");
0611   addColumn(table, "ttIntEta2", m_ttIntEta[1], "");
0612   addColumn(table, "ttIntEta3", m_ttIntEta[2], "");
0613   addColumn(table, "ttIntEta4", m_ttIntEta[3], "");
0614 
0615   addColumn(table, "ttDistToEdge1", m_ttDistToEdge[0], "");
0616   addColumn(table, "ttDistToEdge2", m_ttDistToEdge[1], "");
0617   addColumn(table, "ttDistToEdge3", m_ttDistToEdge[2], "");
0618   addColumn(table, "ttDistToEdge4", m_ttDistToEdge[3], "");
0619   addColumn(table, "ttDistToHVGap1", m_ttDistToHVGap[0], "");
0620   addColumn(table, "ttDistToHVGap2", m_ttDistToHVGap[1], "");
0621   addColumn(table, "ttDistToHVGap3", m_ttDistToHVGap[2], "");
0622   addColumn(table, "ttDistToHVGap4", m_ttDistToHVGap[3], "");
0623 
0624   addColumn(table, "segLocalX1", m_segLocalX[0], "");
0625   addColumn(table, "segLocalX2", m_segLocalX[1], "");
0626   addColumn(table, "segLocalX3", m_segLocalX[2], "");
0627   addColumn(table, "segLocalX4", m_segLocalX[3], "");
0628   addColumn(table, "segLocalY1", m_segLocalY[0], "");
0629   addColumn(table, "segLocalY2", m_segLocalY[1], "");
0630   addColumn(table, "segLocalY3", m_segLocalY[2], "");
0631   addColumn(table, "segLocalY4", m_segLocalY[3], "");
0632   addColumn(table, "segLocalErrorX1", m_segLocalErrorX[0], "");
0633   addColumn(table, "segLocalErrorX2", m_segLocalErrorX[1], "");
0634   addColumn(table, "segLocalErrorX3", m_segLocalErrorX[2], "");
0635   addColumn(table, "segLocalErrorX4", m_segLocalErrorX[3], "");
0636   addColumn(table, "segLocalErrorY1", m_segLocalErrorY[0], "");
0637   addColumn(table, "segLocalErrorY2", m_segLocalErrorY[1], "");
0638   addColumn(table, "segLocalErrorY3", m_segLocalErrorY[2], "");
0639   addColumn(table, "segLocalErrorY4", m_segLocalErrorY[3], "");
0640 
0641   addColumn(table, "ttIntSegResidualLocalX1", m_ttIntSegResidualLocalX[0], "");
0642   addColumn(table, "ttIntSegResidualLocalX2", m_ttIntSegResidualLocalX[1], "");
0643   addColumn(table, "ttIntSegResidualLocalX3", m_ttIntSegResidualLocalX[2], "");
0644   addColumn(table, "ttIntSegResidualLocalX4", m_ttIntSegResidualLocalX[3], "");
0645   addColumn(table, "ttIntSegResidualLocalY1", m_ttIntSegResidualLocalY[0], "");
0646   addColumn(table, "ttIntSegResidualLocalY2", m_ttIntSegResidualLocalY[1], "");
0647   addColumn(table, "ttIntSegResidualLocalY3", m_ttIntSegResidualLocalY[2], "");
0648   addColumn(table, "ttIntSegResidualLocalY4", m_ttIntSegResidualLocalY[3], "");
0649 
0650   ev.put(std::move(table));
0651 }
0652 
0653 float MuCSCTnPFlatTableProducer::computeTrkIso(const reco::MuonIsolation& isolation, float muonPt) {
0654   return isolation.sumPt / muonPt;
0655 }
0656 
0657 float MuCSCTnPFlatTableProducer::computePFIso(const reco::MuonPFIsolation& pfIsolation, float muonPt) {
0658   return (pfIsolation.sumChargedHadronPt +
0659           std::max(0., pfIsolation.sumNeutralHadronEt + pfIsolation.sumPhotonEt - 0.5 * pfIsolation.sumPUPt)) /
0660          muonPt;
0661 }
0662 
0663 bool MuCSCTnPFlatTableProducer::hasTrigger(std::vector<int>& trigIndices,
0664                                            const trigger::TriggerObjectCollection& trigObjs,
0665                                            edm::Handle<trigger::TriggerEvent>& trigEvent,
0666                                            const reco::Muon& muon) {
0667   float dRMatch = 999.;
0668   for (int trigIdx : trigIndices) {
0669     const std::vector<std::string> trigModuleLabels = m_hltConfig.moduleLabels(trigIdx);
0670 
0671     const unsigned trigModuleIndex =
0672         std::find(trigModuleLabels.begin(), trigModuleLabels.end(), "hltBoolEnd") - trigModuleLabels.begin() - 1;
0673     const unsigned hltFilterIndex = trigEvent->filterIndex(edm::InputTag(trigModuleLabels[trigModuleIndex], "", "HLT"));
0674     if (hltFilterIndex < trigEvent->sizeFilters()) {
0675       const trigger::Keys keys = trigEvent->filterKeys(hltFilterIndex);
0676       const trigger::Vids vids = trigEvent->filterIds(hltFilterIndex);
0677       const unsigned nTriggers = vids.size();
0678 
0679       for (unsigned iTrig = 0; iTrig < nTriggers; ++iTrig) {
0680         trigger::TriggerObject trigObj = trigObjs[keys[iTrig]];
0681         float dR = deltaR(muon, trigObj);
0682         if (dR < dRMatch)
0683           dRMatch = dR;
0684       }
0685     }
0686   }
0687 
0688   return dRMatch < 0.1;  //CB should get it programmable
0689 }
0690 
0691 //bool MuCSCTnPFlatTableProducer::muonTagSelection(const reco::Muon & muon,edm::Handle<std::vector<reco::Track>> tracks)
0692 bool MuCSCTnPFlatTableProducer::muonTagSelection(const reco::Muon& muon) {
0693   float ptCut = 10.0;
0694   int trackerHitsCut = 8;
0695   float dxyCut = 2.0;
0696   float dzCut = 24.0;
0697   float chi2Cut = 4.0;
0698 
0699   bool selected = false;
0700   //_muonIso = iso(*muon.track(),tracks);
0701   _muonIso = computePFIso(muon.pfIsolationR04(), muon.pt());
0702 
0703   if (!muon.isTrackerMuon())
0704     return false;
0705   if (!muon.track().isNonnull())
0706     return false;
0707   selected =
0708       ((muon.track()->pt() > ptCut) && (muon.track()->hitPattern().numberOfValidTrackerHits() >= trackerHitsCut) &&
0709        (muon.track()->dxy() < dxyCut) && (std::abs(muon.track()->dz()) < dzCut) &&
0710        (muon.track()->normalizedChi2() < chi2Cut) && _muonIso < 0.1);
0711 
0712   return selected;
0713 }
0714 
0715 bool MuCSCTnPFlatTableProducer::trackProbeSelection(const reco::Track& track,
0716                                                     edm::Handle<std::vector<reco::Track>> tracks) {
0717   float ptCut = 10.0;
0718   int trackerHitsCut = 8;
0719   float dxyCut = 2.0;
0720   float dzCut = 24.0;
0721   float chi2Cut = 4.0;
0722 
0723   bool selected = false;
0724   _trackIso = iso(track, tracks);
0725 
0726   selected =
0727       ((track.pt() > ptCut) && (std::abs(track.eta()) > 0.75) && (std::abs(track.eta()) < 2.55) &&
0728        (track.numberOfValidHits() >= trackerHitsCut) && (track.dxy() < dxyCut) && (std::abs(track.dz()) < dzCut) &&
0729        (track.normalizedChi2() > 0.0) && (track.normalizedChi2() < chi2Cut) && _trackIso < 0.1);
0730 
0731   return selected;
0732 }
0733 
0734 bool MuCSCTnPFlatTableProducer::zSelection(const reco::Muon& muon, const reco::Track& track) {
0735   bool selected = false;
0736 
0737   _zMass = zMass(track, muon);
0738   selected = (track.charge() * muon.charge() == -1 && (_zMass > 75.0) && (_zMass < 120.0));
0739 
0740   return selected;
0741 }
0742 
0743 // get track position at a particular (xy) plane given its z
0744 TrajectoryStateOnSurface MuCSCTnPFlatTableProducer::surfExtrapTrkSam(const reco::Track& track, double z) {
0745   Plane::PositionType pos(0, 0, z);
0746   Plane::RotationType rot;
0747   Plane::PlanePointer myPlane = Plane::build(pos, rot);
0748 
0749   FreeTrajectoryState recoStart = freeTrajStateMuon(track);
0750   TrajectoryStateOnSurface recoProp = propagatorAlong->propagate(recoStart, *myPlane);
0751 
0752   if (!recoProp.isValid())
0753     recoProp = propagatorOpposite->propagate(recoStart, *myPlane);
0754 
0755   return recoProp;
0756 }
0757 
0758 FreeTrajectoryState MuCSCTnPFlatTableProducer::freeTrajStateMuon(const reco::Track& track) {
0759   //no track extras in nanoaod so directly use vx and p
0760   GlobalPoint innerPoint(track.vx(), track.vy(), track.vz());
0761   GlobalVector innerVec(track.px(), track.py(), track.pz());
0762 
0763   GlobalTrajectoryParameters gtPars(innerPoint, innerVec, track.charge(), &*theBField);
0764 
0765   AlgebraicSymMatrix66 cov;
0766   cov *= 1e-20;
0767 
0768   CartesianTrajectoryError tCov(cov);
0769 
0770   return (cov.kRows == 6 ? FreeTrajectoryState(gtPars, tCov) : FreeTrajectoryState(gtPars));
0771 }
0772 
0773 UChar_t MuCSCTnPFlatTableProducer::ringCandidate(Int_t iiStation, Int_t station, Float_t feta, Float_t phi) {
0774   UChar_t ring = 0;
0775 
0776   switch (station) {
0777     case 1:
0778       if (std::abs(feta) >= 0.85 && std::abs(feta) < 1.18) {  //ME13
0779         if (iiStation == 2)
0780           ring = 3;
0781         return ring;
0782       }
0783       if (std::abs(feta) >= 1.18 &&
0784           std::abs(feta) <= 1.5) {  //ME12 if(std::abs(feta)>1.18 && std::abs(feta)<1.7){//ME12
0785         if (iiStation == 1)
0786           ring = 2;
0787         return ring;
0788       }
0789       if (std::abs(feta) > 1.5 && std::abs(feta) < 2.1) {  //ME11
0790         if (iiStation == 0)
0791           ring = 1;
0792         return ring;
0793       }
0794       if (std::abs(feta) >= 2.1 && std::abs(feta) < 2.45) {  //ME11
0795         if (iiStation == 0)
0796           ring = 4;
0797         return ring;
0798       }
0799       break;
0800     case 2:
0801       if (std::abs(feta) > 0.95 && std::abs(feta) < 1.6) {  //ME22
0802         ring = 2;
0803         return ring;
0804       }
0805       if (std::abs(feta) > 1.55 && std::abs(feta) < 2.45) {  //ME21
0806         ring = 1;
0807         return ring;
0808       }
0809       break;
0810     case 3:
0811       if (std::abs(feta) > 1.08 && std::abs(feta) < 1.72) {  //ME32
0812         ring = 2;
0813         return ring;
0814       }
0815       if (std::abs(feta) > 1.69 && std::abs(feta) < 2.45) {  //ME31
0816         ring = 1;
0817         return ring;
0818       }
0819       break;
0820     case 4:
0821       if (std::abs(feta) > 1.78 && std::abs(feta) < 2.45) {  //ME41
0822         ring = 1;
0823         return ring;
0824       }
0825       if (std::abs(feta) > 1.15 && std::abs(feta) <= 1.78) {  //ME42
0826         ring = 2;
0827         return ring;
0828       }
0829       break;
0830     default:
0831       edm::LogError("") << "Invalid station: " << station << std::endl;
0832       break;
0833   }
0834   return 0;
0835 }
0836 
0837 UChar_t MuCSCTnPFlatTableProducer::thisChamberCandidate(UChar_t station, UChar_t ring, Float_t phi) {
0838   //    cout <<"\t\t TPTrackMuonSys::thisChamberCandidate..."<<endl;
0839 
0840   //search for chamber candidate based on CMS IN-2007/024
0841   //10 deg chambers are ME1,ME22,ME32,ME42 chambers; 20 deg chambers are ME21,31,41 chambers
0842   //Chambers one always starts from approx -5 deg.
0843   const UChar_t nVal = (station > 1 && ring == 1) ? 18 : 36;
0844   const Float_t ChamberSpan = 2 * M_PI / nVal;
0845   Float_t dphi = phi + M_PI / 36;
0846   while (dphi >= 2 * M_PI)
0847     dphi -= 2 * M_PI;
0848   while (dphi < 0)
0849     dphi += 2 * M_PI;
0850   UChar_t ChCand = floor(dphi / ChamberSpan) + 1;
0851   return ChCand > nVal ? nVal : ChCand;
0852 }
0853 
0854 Float_t MuCSCTnPFlatTableProducer::TrajectoryDistToSeg(TrajectoryStateOnSurface TrajSuf,
0855                                                        CSCSegmentCollection::const_iterator segIt) {
0856   if (!TrajSuf.isValid())
0857     return 9999.;
0858   const GeomDet* gdet = m_cscGeometry->idToDet((CSCDetId)(*segIt).cscDetId());
0859   LocalPoint localTTPos = gdet->surface().toLocal(TrajSuf.freeState()->position());
0860   LocalPoint localSegPos = (*segIt).localPosition();
0861   Float_t CSCdeltaX = localSegPos.x() - localTTPos.x();
0862   Float_t CSCdeltaY = localSegPos.y() - localTTPos.y();
0863   return sqrt(pow(CSCdeltaX, 2) + pow(CSCdeltaY, 2));
0864 }
0865 
0866 TrajectoryStateOnSurface* MuCSCTnPFlatTableProducer::matchTTwithCSCSeg(const reco::Track& track,
0867                                                                        edm::Handle<CSCSegmentCollection> cscSegments,
0868                                                                        CSCSegmentCollection::const_iterator& cscSegOut,
0869                                                                        CSCDetId& idCSC) {
0870   TrajectoryStateOnSurface* TrajSuf = nullptr;
0871   Float_t deltaCSCR = 9999.;
0872   for (CSCSegmentCollection::const_iterator segIt = cscSegments->begin(); segIt != cscSegments->end(); segIt++) {
0873     CSCDetId id = (CSCDetId)(*segIt).cscDetId();
0874 
0875     if (idCSC.endcap() != id.endcap())
0876       continue;
0877     if (idCSC.station() != id.station())
0878       continue;
0879     if (idCSC.chamber() != id.chamber())
0880       continue;
0881 
0882     Bool_t ed1 =
0883         (idCSC.station() == 1) && ((idCSC.ring() == 1 || idCSC.ring() == 4) && (id.ring() == 1 || id.ring() == 4));
0884     Bool_t ed2 =
0885         (idCSC.station() == 1) && ((idCSC.ring() == 2 && id.ring() == 2) || (idCSC.ring() == 3 && id.ring() == 3));
0886     Bool_t ed3 = (idCSC.station() != 1) && (idCSC.ring() == id.ring());
0887     Bool_t TMCSCMatch = (ed1 || ed2 || ed3);
0888 
0889     if (!TMCSCMatch)
0890       continue;
0891 
0892     const CSCChamber* cscchamber = m_cscGeometry->chamber(id);
0893 
0894     if (!cscchamber)
0895       continue;
0896 
0897     TrajectoryStateOnSurface TrajSuf_ = surfExtrapTrkSam(track, cscchamber->toGlobal((*segIt).localPosition()).z());
0898     Float_t dR_ = std::abs(TrajectoryDistToSeg(TrajSuf_, segIt));
0899     if (dR_ < deltaCSCR) {
0900       delete TrajSuf;
0901       TrajSuf = new TrajectoryStateOnSurface(TrajSuf_);
0902       deltaCSCR = dR_;
0903       cscSegOut = segIt;
0904     }
0905   }  //loop over segments
0906 
0907   return TrajSuf;
0908 }
0909 
0910 std::vector<Float_t> MuCSCTnPFlatTableProducer::GetEdgeAndDistToGap(const reco::Track& track, CSCDetId& detid) {
0911   std::vector<Float_t> result(4, 9999.);
0912   result[3] = -9999;
0913   const GeomDet* gdet = m_cscGeometry->idToDet(detid);
0914   TrajectoryStateOnSurface tsos = surfExtrapTrkSam(track, gdet->surface().position().z());
0915   if (!tsos.isValid())
0916     return result;
0917   LocalPoint localTTPos = gdet->surface().toLocal(tsos.freeState()->position());
0918   const CSCWireTopology* wireTopology = m_cscGeometry->layer(detid)->geometry()->wireTopology();
0919   Float_t wideWidth = wireTopology->wideWidthOfPlane();
0920   Float_t narrowWidth = wireTopology->narrowWidthOfPlane();
0921   Float_t length = wireTopology->lengthOfPlane();
0922   // If slanted, there is no y offset between local origin and symmetry center of wire plane
0923   Float_t yOfFirstWire = std::abs(wireTopology->wireAngle()) > 1.E-06 ? -0.5 * length : wireTopology->yOfWire(1);
0924   // y offset between local origin and symmetry center of wire plane
0925   Float_t yCOWPOffset = yOfFirstWire + 0.5 * length;
0926   // tangent of the incline angle from inside the trapezoid
0927   Float_t tangent = (wideWidth - narrowWidth) / (2. * length);
0928   // y position wrt bottom of trapezoid
0929   Float_t yPrime = localTTPos.y() + std::abs(yOfFirstWire);
0930   // half trapezoid width at y' is 0.5 * narrowWidth + x side of triangle with the above tangent and side y'
0931   Float_t halfWidthAtYPrime = 0.5 * narrowWidth + yPrime * tangent;
0932   // x offset between local origin and symmetry center of wire plane is zero
0933   // 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.
0934   Float_t edgex = std::abs(localTTPos.x()) - halfWidthAtYPrime;
0935   Float_t edgey = std::abs(localTTPos.y() - yCOWPOffset) - 0.5 * length;
0936   LocalError localTTErr = tsos.localError().positionError();
0937   if (edgex > edgey) {
0938     result[0] = edgex;
0939     result[1] = sqrt(localTTErr.xx());
0940     //result[1] = sqrt(tsos.cartesianError().position().cxx());
0941   } else {
0942     result[0] = edgey;
0943     result[1] = sqrt(localTTErr.yy());
0944     //result[1] = sqrt(tsos.cartesianError().position().cyy());
0945   }
0946   result[2] = YDistToHVDeadZone(localTTPos.y(), detid.station() * 10 + detid.ring());
0947   result[3] = sqrt(localTTErr.yy());
0948   return result;  //return values: 1-edge;2-err of edge;3-disttogap;4-err of dist to gap
0949 }
0950 
0951 //deadzone center is according to http://cmssdt.cern.ch/SDT/lxr/source/RecoLocalMuon/CSCEfficiency/src/CSCEfficiency.cc#605
0952 //wire spacing is according to CSCTDR
0953 Float_t MuCSCTnPFlatTableProducer::YDistToHVDeadZone(Float_t yLocal, Int_t StationAndRing) {
0954   //the ME11 wires are not parallel to x, but no gap
0955   //chamber edges are not included.
0956   const Float_t deadZoneCenterME1_2[2] = {-32.88305, 32.867423};
0957   const Float_t deadZoneCenterME1_3[2] = {-22.7401, 27.86665};
0958   const Float_t deadZoneCenterME2_1[2] = {-27.47, 33.67};
0959   const Float_t deadZoneCenterME3_1[2] = {-36.21, 23.68};
0960   const Float_t deadZoneCenterME4_1[2] = {-26.14, 23.85};
0961   const Float_t deadZoneCenterME234_2[4] = {-81.8744, -21.18165, 39.51105, 100.2939};
0962   const Float_t* deadZoneCenter;
0963   Float_t deadZoneHeightHalf = 0.32 * 7 / 2;  // wire spacing * (wires missing + 1)/2
0964   Float_t minY = 999999.;
0965   UChar_t nGaps = 2;
0966   switch (std::abs(StationAndRing)) {
0967     case 11:
0968     case 14:
0969       return 162;  //the height of ME11
0970       break;
0971     case 12:
0972       deadZoneCenter = deadZoneCenterME1_2;
0973       break;
0974     case 13:
0975       deadZoneCenter = deadZoneCenterME1_3;
0976       break;
0977     case 21:
0978       deadZoneCenter = deadZoneCenterME2_1;
0979       break;
0980     case 31:
0981       deadZoneCenter = deadZoneCenterME3_1;
0982       break;
0983     case 41:
0984       deadZoneCenter = deadZoneCenterME4_1;
0985       break;
0986     default:
0987       deadZoneCenter = deadZoneCenterME234_2;
0988       nGaps = 4;
0989   }
0990   for (UChar_t iGap = 0; iGap < nGaps; iGap++) {
0991     Float_t newMinY = yLocal < deadZoneCenter[iGap] ? deadZoneCenter[iGap] - deadZoneHeightHalf - yLocal
0992                                                     : yLocal - (deadZoneCenter[iGap] + deadZoneHeightHalf);
0993     if (newMinY < minY)
0994       minY = newMinY;
0995   }
0996   return minY;
0997 }
0998 
0999 double MuCSCTnPFlatTableProducer::iso(const reco::Track& track, edm::Handle<std::vector<reco::Track>> tracks) {
1000   double isoSum = 0.0;
1001   for (const auto& track2 : (*tracks)) {
1002     double dR = calcDeltaR(track.eta(), track2.eta(), track.phi(), track2.phi());
1003     if (track2.pt() > 1.0 && dR > 0.001 && dR < 0.3)
1004       isoSum += track2.pt();
1005   }
1006   return isoSum / track.pt();
1007 }
1008 
1009 double MuCSCTnPFlatTableProducer::calcDeltaR(double eta1, double eta2, double phi1, double phi2) {
1010   double deta = eta1 - eta2;
1011   if (phi1 < 0)
1012     phi1 += 2.0 * M_PI;
1013   if (phi2 < 0)
1014     phi2 += 2.0 * M_PI;
1015   double dphi = phi1 - phi2;
1016   if (dphi > M_PI)
1017     dphi -= 2. * M_PI;
1018   else if (dphi < -M_PI)
1019     dphi += 2. * M_PI;
1020   return std::sqrt(deta * deta + dphi * dphi);
1021 }
1022 
1023 double MuCSCTnPFlatTableProducer::zMass(const reco::Track& track, const reco::Muon& muon) {
1024   double zMass = -99.0;
1025   double mMu = 0.1134289256;
1026 
1027   zMass = std::pow((std::sqrt(std::pow(muon.p(), 2) + mMu * mMu) + std::sqrt(std::pow(track.p(), 2) + mMu * mMu)), 2) -
1028           (std::pow((muon.px() + track.px()), 2) + std::pow((muon.py() + track.py()), 2) +
1029            std::pow((muon.pz() + track.pz()), 2));
1030 
1031   return std::sqrt(zMass);
1032 }
1033 
1034 #include "FWCore/PluginManager/interface/ModuleDef.h"
1035 #include "FWCore/Framework/interface/MakerMacros.h"
1036 
1037 DEFINE_FWK_MODULE(MuCSCTnPFlatTableProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!