Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-07-25 23:16:07 

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