Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-09-07 04:37:47

0001 #include "RecoMuon/MuonIdentification/interface/MuonMesh.h"
0002 #include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
0003 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0004 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0005 
0006 #include <utility>
0007 #include <algorithm>
0008 
0009 MuonMesh::MuonMesh(const edm::ParameterSet& parm)
0010     : doME1a(parm.getParameter<bool>("ME1a")),
0011       doOverlaps(parm.getParameter<bool>("Overlap")),
0012       doClustering(parm.getParameter<bool>("Clustering")),
0013       OverlapDPhi(parm.getParameter<double>("OverlapDPhi")),
0014       OverlapDTheta(parm.getParameter<double>("OverlapDTheta")),
0015       ClusterDPhi(parm.getParameter<double>("ClusterDPhi")),
0016       ClusterDTheta(parm.getParameter<double>("ClusterDTheta")) {}
0017 
0018 void MuonMesh::fillMesh(std::vector<reco::Muon>* inputMuons) {
0019   for (std::vector<reco::Muon>::iterator muonIter1 = inputMuons->begin(); muonIter1 != inputMuons->end(); ++muonIter1) {
0020     if (muonIter1->isTrackerMuon()) {
0021       mesh_[&*muonIter1];  // create this entry if it's a tracker muon
0022       for (std::vector<reco::Muon>::iterator muonIter2 = inputMuons->begin(); muonIter2 != inputMuons->end();
0023            ++muonIter2) {
0024         if (muonIter2->isTrackerMuon()) {
0025           if (muonIter2 != muonIter1) {
0026             // now fill all the edges for muon1 based on overlaps with muon2
0027             for (std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = muonIter1->matches().begin();
0028                  chamberIter1 != muonIter1->matches().end();
0029                  ++chamberIter1) {
0030               for (std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
0031                    segmentIter1 != chamberIter1->segmentMatches.end();
0032                    ++segmentIter1) {
0033                 for (std::vector<reco::MuonChamberMatch>::iterator chamberIter2 = muonIter2->matches().begin();
0034                      chamberIter2 != muonIter2->matches().end();
0035                      ++chamberIter2) {
0036                   for (std::vector<reco::MuonSegmentMatch>::iterator segmentIter2 =
0037                            chamberIter2->segmentMatches.begin();
0038                        segmentIter2 != chamberIter2->segmentMatches.end();
0039                        ++segmentIter2) {
0040                     //if(segmentIter1->mask & 0x1e0000 && segmentIter2->mask &0x1e0000) {
0041 
0042                     bool addsegment(false);
0043 
0044                     if (segmentIter1->cscSegmentRef.isNonnull() && segmentIter2->cscSegmentRef.isNonnull()) {
0045                       if (doME1a && isDuplicateOf(segmentIter1->cscSegmentRef, segmentIter2->cscSegmentRef) &&
0046                           CSCDetId(chamberIter1->id).ring() == 4 && CSCDetId(chamberIter2->id).ring() == 4 &&
0047                           chamberIter1->id == chamberIter2->id) {
0048                         addsegment = true;
0049                         //std::cout << "\tME1/a sharing detected." << std::endl;
0050                       }
0051 
0052                       if (doOverlaps &&
0053                           isDuplicateOf(std::make_pair(CSCDetId(chamberIter1->id), segmentIter1->cscSegmentRef),
0054                                         std::make_pair(CSCDetId(chamberIter2->id), segmentIter2->cscSegmentRef))) {
0055                         addsegment = true;
0056                         //std::cout << "\tChamber Overlap sharing detected." << std::endl;
0057                       }
0058 
0059                       if (doClustering &&
0060                           isClusteredWith(std::make_pair(CSCDetId(chamberIter1->id), segmentIter1->cscSegmentRef),
0061                                           std::make_pair(CSCDetId(chamberIter2->id), segmentIter2->cscSegmentRef))) {
0062                         addsegment = true;
0063                         //std::cout << "\tCluster sharing detected." << std::endl;
0064                       }
0065                       //std::cout << std::endl;
0066                     }  // has valid csc segment ref
0067 
0068                     if (addsegment) {  // add segment if clusters/overlaps/replicant and doesn't already exist
0069 
0070                       if (find(mesh_[&*muonIter1].begin(),
0071                                mesh_[&*muonIter1].end(),
0072                                std::make_pair(&*muonIter2, std::make_pair(&*chamberIter2, &*segmentIter2))) ==
0073                           mesh_[&*muonIter1].end()) {
0074                         mesh_[&*muonIter1].push_back(
0075                             std::make_pair(&*muonIter2, std::make_pair(&*chamberIter2, &*segmentIter2)));
0076                       }  // find
0077                     }  // add segment?
0078                     //} // both segments won arbitration
0079                   }  // segmentIter 2
0080                 }  // chamberIter2
0081               }  //segmentIter1
0082             }  // chamberIter1
0083 
0084           }  // if different muon
0085         }  // is tracker muon
0086       }  // muonIter2
0087 
0088     }  // is tracker muon
0089   }  // muonIter1
0090 
0091   // special cases
0092 
0093   // one muon: mark all segments belonging to a muon as cleaned as there are no other muons to fight with
0094   if (mesh_.size() == 1) {
0095     for (std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = mesh_.begin()->first->matches().begin();
0096          chamberIter1 != mesh_.begin()->first->matches().end();
0097          ++chamberIter1) {
0098       for (std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
0099            segmentIter1 != chamberIter1->segmentMatches.end();
0100            ++segmentIter1) {
0101         segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0102       }  // segmentIter1
0103     }  // chamberIter1
0104   }  // if only one tracker muon set winner bit boosted arbitration
0105 
0106   // segments that are not shared amongst muons and the have won all segment arbitration flags need to be promoted
0107   // also promote DT segments
0108   if (mesh_.size() > 1) {
0109     for (MeshType::iterator i = mesh_.begin(); i != mesh_.end(); ++i) {
0110       for (std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = i->first->matches().begin();
0111            chamberIter1 != i->first->matches().end();
0112            ++chamberIter1) {
0113         for (std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
0114              segmentIter1 != chamberIter1->segmentMatches.end();
0115              ++segmentIter1) {
0116           bool shared(false);
0117 
0118           for (AssociationType::iterator j = i->second.begin(); j != i->second.end(); ++j) {
0119             if (segmentIter1->cscSegmentRef.isNonnull() && j->second.second->cscSegmentRef.isNonnull()) {
0120               if (chamberIter1->id.subdetId() == MuonSubdetId::CSC &&
0121                   j->second.first->id.subdetId() == MuonSubdetId::CSC) {
0122                 CSCDetId segIterId(chamberIter1->id), shareId(j->second.first->id);
0123 
0124                 if (doOverlaps && isDuplicateOf(std::make_pair(segIterId, segmentIter1->cscSegmentRef),
0125                                                 std::make_pair(shareId, j->second.second->cscSegmentRef)))
0126                   shared = true;
0127 
0128                 if (doME1a && isDuplicateOf(segmentIter1->cscSegmentRef, j->second.second->cscSegmentRef) &&
0129                     segIterId.ring() == 4 && shareId.ring() == 4 && segIterId == segIterId)
0130                   shared = true;
0131 
0132                 if (doClustering &&
0133                     isClusteredWith(std::make_pair(CSCDetId(chamberIter1->id), segmentIter1->cscSegmentRef),
0134                                     std::make_pair(CSCDetId(j->second.first->id), j->second.second->cscSegmentRef)))
0135                   shared = true;
0136               }  // in CSCs?
0137             }  // cscSegmentRef non null?
0138           }  // j
0139 
0140           // Promote segments which have won all arbitration and are not shared or are DT segments
0141           if (((segmentIter1->mask & 0x1e0000) == 0x1e0000 && !shared) ||
0142               (chamberIter1->id.subdetId() == MuonSubdetId::DT && (segmentIter1->mask & 0x1e0000)))
0143             segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0144 
0145         }  // segmentIter1
0146       }  // chamberIter1
0147     }  // i
0148   }  // if non-trivial case
0149 }
0150 
0151 void MuonMesh::pruneMesh() {
0152   for (MeshType::iterator i = mesh_.begin(); i != mesh_.end(); ++i) {
0153     for (AssociationType::iterator j = i->second.begin(); j != i->second.end(); ++j) {
0154       for (std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = i->first->matches().begin();
0155            chamberIter1 != i->first->matches().end();
0156            ++chamberIter1) {
0157         for (std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
0158              segmentIter1 != chamberIter1->segmentMatches.end();
0159              ++segmentIter1) {
0160           if (j->second.second->cscSegmentRef.isNonnull() && segmentIter1->cscSegmentRef.isNonnull()) {
0161             //UNUSED:       bool me1a(false), overlap(false), cluster(false);
0162 
0163             // remove physical overlap duplicates first
0164             if (doOverlaps &&
0165                 isDuplicateOf(std::make_pair(CSCDetId(chamberIter1->id), segmentIter1->cscSegmentRef),
0166                               std::make_pair(CSCDetId(j->second.first->id), j->second.second->cscSegmentRef))) {
0167               if (i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) >
0168                   j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000)) {
0169                 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByOvlClean);
0170                 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0171 
0172                 //UNUSED:       overlap = true;
0173               } else if (i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ==
0174                          j->first->numberOfMatches(
0175                              (reco::Muon::ArbitrationType)0x1e0000)) {  // muon with more matched stations wins
0176 
0177                 if ((segmentIter1->mask & 0x1e0000) >
0178                     (j->second.second->mask & 0x1e0000)) {  // segment with better match wins
0179 
0180                   segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByOvlClean);
0181                   segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0182 
0183                   //UNUSED:       overlap = true;
0184                 } else {  // ??
0185                   // leave this available for later
0186                 }
0187               }  // overlap duplicate resolution
0188             }  // is overlap duplicate
0189 
0190             // do ME1/a arbitration second since the tie breaker depends on other stations
0191             // Unlike the other cleanings this one removes the bits from segments associated to tracks which
0192             // fail cleaning. (Instead of setting bits for the segments which win.)
0193             if (doME1a && isDuplicateOf(segmentIter1->cscSegmentRef, j->second.second->cscSegmentRef) &&
0194                 CSCDetId(chamberIter1->id).ring() == 4 && CSCDetId(j->second.first->id).ring() == 4 &&
0195                 chamberIter1->id == j->second.first->id) {
0196               if (j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) <
0197                   i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000)) {
0198                 for (AssociationType::iterator AsscIter1 = i->second.begin(); AsscIter1 != i->second.end();
0199                      ++AsscIter1) {
0200                   if (AsscIter1->second.second->cscSegmentRef.isNonnull())
0201                     if (j->first == AsscIter1->first && j->second.first == AsscIter1->second.first &&
0202                         isDuplicateOf(segmentIter1->cscSegmentRef, AsscIter1->second.second->cscSegmentRef)) {
0203                       AsscIter1->second.second->mask &= ~reco::MuonSegmentMatch::BelongsToTrackByME1aClean;
0204                     }
0205                 }
0206 
0207                 //UNUSED:       me1a = true;
0208               } else if (j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ==
0209                          i->first->numberOfMatches(
0210                              (reco::Muon::ArbitrationType)0x1e0000)) {  // muon with best arbitration wins
0211 
0212                 bool bestArb(true);
0213 
0214                 for (AssociationType::iterator AsscIter1 = i->second.begin(); AsscIter1 != i->second.end();
0215                      ++AsscIter1) {
0216                   if (AsscIter1->second.second->cscSegmentRef.isNonnull())
0217                     if (j->first == AsscIter1->first && j->second.first == AsscIter1->second.first &&
0218                         isDuplicateOf(segmentIter1->cscSegmentRef, AsscIter1->second.second->cscSegmentRef) &&
0219                         (segmentIter1->mask & 0x1e0000) < (AsscIter1->second.second->mask & 0x1e0000))
0220                       bestArb = false;
0221                 }
0222 
0223                 if (bestArb) {
0224                   for (AssociationType::iterator AsscIter1 = i->second.begin(); AsscIter1 != i->second.end();
0225                        ++AsscIter1) {
0226                     if (AsscIter1->second.second->cscSegmentRef.isNonnull())
0227                       if (j->first == AsscIter1->first && j->second.first == AsscIter1->second.first &&
0228                           isDuplicateOf(segmentIter1->cscSegmentRef, AsscIter1->second.second->cscSegmentRef)) {
0229                         AsscIter1->second.second->mask &= ~reco::MuonSegmentMatch::BelongsToTrackByME1aClean;
0230                       }
0231                   }
0232                 }
0233                 //UNUSED        me1a = true;
0234 
0235               }  // ME1/a duplicate resolution
0236             }  // is ME1/aduplicate?
0237 
0238             if (doClustering &&
0239                 isClusteredWith(std::make_pair(CSCDetId(chamberIter1->id), segmentIter1->cscSegmentRef),
0240                                 std::make_pair(CSCDetId(j->second.first->id), j->second.second->cscSegmentRef))) {
0241               if (i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) >
0242                   j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000)) {
0243                 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
0244                 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0245 
0246                 //UNUSED:        cluster = true;
0247               } else if (i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) <
0248                          j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000)) {
0249                 j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
0250                 j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0251 
0252                 //UNUSED:       cluster = true;
0253               } else {  // muon with more matched stations wins
0254 
0255                 if ((segmentIter1->mask & 0x1e0000) >
0256                     (j->second.second->mask & 0x1e0000)) {  // segment with better match wins
0257 
0258                   segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
0259                   segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0260 
0261                   //UNUSED:        cluster = true;
0262                 } else if ((segmentIter1->mask & 0x1e0000) < (j->second.second->mask & 0x1e0000)) {  //
0263 
0264                   j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
0265                   j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0266 
0267                   //UNUSED:        cluster = true;
0268                 } else {
0269                 }
0270               }  // cluster sharing resolution
0271 
0272             }  // is clustered with?
0273 
0274           }  // csc ref nonnull
0275         }  // segmentIter1
0276       }  // chamberIter1
0277     }  // j, associated segments iterator
0278   }  // i, map iterator
0279 
0280   // final step: make sure everything that's won a cleaning flag has the "BelongsToTrackByCleaning" flag
0281 
0282   for (MeshType::iterator i = mesh_.begin(); i != mesh_.end(); ++i) {
0283     for (std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = i->first->matches().begin();
0284          chamberIter1 != i->first->matches().end();
0285          ++chamberIter1) {
0286       for (std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
0287            segmentIter1 != chamberIter1->segmentMatches.end();
0288            ++segmentIter1) {
0289         // set cleaning bit if initial no cleaning bit but there are cleaning algorithm bits set.
0290         if (!segmentIter1->isMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning) && segmentIter1->isMask(0xe00000))
0291           segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
0292       }  // segmentIter1
0293     }  // chamberIter1
0294   }  // i
0295 }
0296 
0297 bool MuonMesh::isDuplicateOf(const CSCSegmentRef& lhs, const CSCSegmentRef& rhs)
0298     const  // this isDuplicateOf() deals with duplicate segments in ME1/a
0299 {
0300   bool result(false);
0301 
0302   if (!lhs->isME11a_duplicate())
0303     return result;
0304 
0305   std::vector<CSCSegment> lhs_duplicates = lhs->duplicateSegments();
0306 
0307   if (fabs(lhs->localPosition().x() - rhs->localPosition().x()) < 1E-3 &&
0308       fabs(lhs->localPosition().y() - rhs->localPosition().y()) < 1E-3 &&
0309       fabs(lhs->localDirection().x() / lhs->localDirection().z() -
0310            rhs->localDirection().x() / rhs->localDirection().z()) < 1E-3 &&
0311       fabs(lhs->localDirection().y() / lhs->localDirection().z() -
0312            rhs->localDirection().y() / rhs->localDirection().z()) < 1E-3 &&
0313       fabs(lhs->localPositionError().xx() - rhs->localPositionError().xx()) < 1E-3 &&
0314       fabs(lhs->localPositionError().yy() - rhs->localPositionError().yy()) < 1E-3 &&
0315       fabs(lhs->localDirectionError().xx() - rhs->localDirectionError().xx()) < 1E-3 &&
0316       fabs(lhs->localDirectionError().yy() - rhs->localDirectionError().yy()) < 1E-3)
0317     result = true;
0318 
0319   for (std::vector<CSCSegment>::const_iterator segIter1 = lhs_duplicates.begin(); segIter1 != lhs_duplicates.end();
0320        ++segIter1) {  // loop over lhs duplicates
0321 
0322     if (fabs(segIter1->localPosition().x() - rhs->localPosition().x()) < 1E-3 &&
0323         fabs(segIter1->localPosition().y() - rhs->localPosition().y()) < 1E-3 &&
0324         fabs(segIter1->localDirection().x() / segIter1->localDirection().z() -
0325              rhs->localDirection().x() / rhs->localDirection().z()) < 1E-3 &&
0326         fabs(segIter1->localDirection().y() / segIter1->localDirection().z() -
0327              rhs->localDirection().y() / rhs->localDirection().z()) < 1E-3 &&
0328         fabs(segIter1->localPositionError().xx() - rhs->localPositionError().xx()) < 1E-3 &&
0329         fabs(segIter1->localPositionError().yy() - rhs->localPositionError().yy()) < 1E-3 &&
0330         fabs(segIter1->localDirectionError().xx() - rhs->localDirectionError().xx()) < 1E-3 &&
0331         fabs(segIter1->localDirectionError().yy() - rhs->localDirectionError().yy()) < 1E-3)
0332       result = true;
0333     /*
0334     if(fabs(segIter1->localPosition().x()        - rhs->localPosition().x()      ) < 2*sqrt(segIter1->localPositionError().xx()) &&
0335        fabs(segIter1->localPosition().y()        - rhs->localPosition().y()      ) < 2*sqrt(segIter1->localPositionError().yy()) &&
0336        fabs(segIter1->localDirection().x()/segIter1->localDirection().z()    - rhs->localDirection().x()/rhs->localDirection().z()   ) 
0337        < 2*std::sqrt(std::max(segIter1->localDirectionError().yy(),rhs->localDirectionError().xx())) &&
0338        fabs(segIter1->localDirection().y()/segIter1->localDirection().z()    - rhs->localDirection().y()/rhs->localDirection().z()   ) 
0339        < 2*std::sqrt(std::max(segIter1->localDirectionError().yy(),rhs->localDirectionError().yy())))
0340       result = true;
0341     */
0342 
0343   }  // loop over duplicates
0344 
0345   return result;
0346 }
0347 
0348 bool MuonMesh::isDuplicateOf(const std::pair<CSCDetId, CSCSegmentRef>& rhs,
0349                              const std::pair<CSCDetId, CSCSegmentRef>& lhs)
0350     const  // this isDuplicateOf() deals with "overlapping chambers" duplicates
0351 {
0352   bool result(false);
0353 
0354   // try to implement the simple case first just back-to-back segments without treatment of ME1/a ganging
0355   // ME1a should be a simple extension of this
0356 
0357   if (rhs.first.endcap() == lhs.first.endcap() && rhs.first.station() == lhs.first.station() &&
0358       rhs.first.ring() == lhs.first.ring()) {  // if same endcap,station,ring (minimal requirement for ovl candidate)
0359     /*
0360     std::cout << "Chamber 1: " << rhs.first << std::endl 
0361           << "Chamber 2: " << lhs.first << std::endl;
0362 
0363     std::cout << "Same endcap,ring,station." << std::endl;    
0364     */
0365     //create neighboring chamber labels, treat ring as (Z mod 36 or 18) + 1 number line: left, right defined accordingly.
0366     unsigned modulus = ((rhs.first.ring() != 1 || rhs.first.station() == 1) ? 36 : 18);
0367     int left_neighbor = (((rhs.first.chamber() - 1 + modulus) % modulus == 0)
0368                              ? modulus
0369                              : (rhs.first.chamber() - 1 + modulus) % modulus);  // + modulus to ensure positivity
0370     int right_neighbor = (((rhs.first.chamber() + 1) % modulus == 0) ? modulus : (rhs.first.chamber() + 1) % modulus);
0371 
0372     if (lhs.first.chamber() == left_neighbor ||
0373         lhs.first.chamber() == right_neighbor) {  // if this is a neighboring chamber then it can be an overlap
0374 
0375       std::vector<CSCSegment> thesegments;
0376       thesegments.push_back(*(lhs.second));
0377       /*
0378       if(lhs.second->isME11a_duplicate())
0379     thesegments.insert(thesegments.begin(),
0380                lhs.second->duplicateSegments().begin(),
0381                lhs.second->duplicateSegments().end());
0382       */
0383 
0384       //std::cout << "lhs is in neighoring chamber of rhs." << std::endl;
0385 
0386       // rhs local direction info
0387       /*
0388       double rhs_dydz = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).y()/
0389                     geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).z();
0390       double rhs_dxdz = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).x()/
0391                     geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).z();
0392       double rhs_dydz_err = rhs.second->localDirectionError().yy();
0393       double rhs_dxdz_err = rhs.second->localDirectionError().xx();
0394       */
0395 
0396       //rhs global position info
0397       double rhs_phi = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).phi();
0398       double rhs_theta = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).theta();
0399       double rhs_z = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).z();
0400 
0401       for (std::vector<CSCSegment>::const_iterator ilhs = thesegments.begin(); ilhs != thesegments.end(); ++ilhs) {
0402         // lhs local direction info
0403         /*
0404     double lhs_dydz = geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).y()/
0405                       geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).z();
0406     double lhs_dxdz = geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).x()/
0407                       geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).z();
0408     double lhs_dydz_err = ilhs->localDirectionError().yy();
0409     double lhs_dxdz_err = ilhs->localDirectionError().xx();
0410     */
0411 
0412         //lhs global position info
0413         double lhs_phi = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).phi();
0414         double lhs_theta = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).theta();
0415         double lhs_z = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).z();
0416         /*
0417       std::cout << "RHS Segment Parameters:" << std::endl
0418       << "\t RHS Phi   : " << rhs_phi << std::endl
0419       << "\t RHS Theta : " << rhs_theta << std::endl
0420       << "\t RHS dx/dz : " << rhs_dxdz << " +- " << rhs_dxdz_err << std::endl
0421       << "\t RHS dy/dz : " << rhs_dydz << " +- " << rhs_dydz_err << std::endl; 
0422       
0423       std::cout << "LHS Segment Parameters:" << std::endl
0424       << "\t LHS Phi   : " << lhs_phi << std::endl
0425       << "\t LHS Theta : " << lhs_theta << std::endl
0426       << "\t LHS dx/dz : " << lhs_dxdz << " +- " << lhs_dxdz_err << std::endl
0427       << "\t LHS dy/dz : " << lhs_dydz << " +- " << lhs_dydz_err << std::endl; 
0428     */
0429 
0430         double phidiff =
0431             (fabs(rhs_phi - lhs_phi) > 2 * M_PI ? fabs(rhs_phi - lhs_phi) - 2 * M_PI : fabs(rhs_phi - lhs_phi));
0432 
0433         if (phidiff < OverlapDPhi && fabs(rhs_theta - lhs_theta) < OverlapDTheta && fabs(rhs_z) < fabs(lhs_z) &&
0434             rhs_z * lhs_z > 0)  // phi overlap region is 3.5 degrees and rhs is infront of lhs
0435           result = true;
0436       }  // loop over duplicate segments
0437     }  // neighboring chamber
0438   }  // same endcap,station,ring
0439 
0440   return result;
0441 }
0442 
0443 bool MuonMesh::isClusteredWith(const std::pair<CSCDetId, CSCSegmentRef>& lhs,
0444                                const std::pair<CSCDetId, CSCSegmentRef>& rhs) const {
0445   bool result(false);
0446 
0447   // try to implement the simple case first just back-to-back segments without treatment of ME1/a ganging
0448   // ME1a should be a simple extension of this
0449 
0450   //std::cout << lhs.first << ' ' << rhs.first << std::endl;
0451 
0452   if (rhs.first.endcap() == lhs.first.endcap() && lhs.first.station() < rhs.first.station()) {
0453     std::vector<CSCSegment> thesegments;
0454     thesegments.push_back(*(lhs.second));
0455     /*
0456       if(lhs.second->isME11a_duplicate())
0457     thesegments.insert(thesegments.begin(),
0458                lhs.second->duplicateSegments().begin(),
0459                lhs.second->duplicateSegments().end());
0460       */
0461     //std::cout << "lhs is in neighoring chamber of rhs." << std::endl;
0462 
0463     // rhs local direction info
0464     /*
0465       double rhs_dydz = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).y()/
0466                     geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).z();
0467       double rhs_dxdz = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).x()/
0468                     geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).z();
0469       double rhs_dydz_err = rhs.second->localDirectionError().yy();
0470       double rhs_dxdz_err = rhs.second->localDirectionError().xx();
0471       */
0472 
0473     //rhs global position info
0474     double rhs_phi = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).phi();
0475     double rhs_theta = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).theta();
0476 
0477     for (std::vector<CSCSegment>::const_iterator ilhs = thesegments.begin(); ilhs != thesegments.end(); ++ilhs) {
0478       // lhs local direction info
0479       /*
0480     double lhs_dydz = geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).y()/
0481                       geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).z();
0482     double lhs_dxdz = geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).x()/
0483                       geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).z();
0484     double lhs_dydz_err = ilhs->localDirectionError().yy();
0485     double lhs_dxdz_err = ilhs->localDirectionError().xx();
0486     */
0487 
0488       //lhs global position info
0489       double lhs_phi = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).phi();
0490       double lhs_theta = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).theta();
0491       /*
0492       std::cout << "RHS Segment Parameters:" << std::endl
0493       << "\t RHS Phi   : " << rhs_phi << std::endl
0494       << "\t RHS Theta : " << rhs_theta << std::endl
0495       << "\t RHS dx/dz : " << rhs_dxdz << " +- " << rhs_dxdz_err << std::endl
0496       << "\t RHS dy/dz : " << rhs_dydz << " +- " << rhs_dydz_err << std::endl; 
0497       
0498       std::cout << "LHS Segment Parameters:" << std::endl
0499       << "\t LHS Phi   : " << lhs_phi << std::endl
0500       << "\t LHS Theta : " << lhs_theta << std::endl
0501       << "\t LHS dx/dz : " << lhs_dxdz << " +- " << lhs_dxdz_err << std::endl
0502       << "\t LHS dy/dz : " << lhs_dydz << " +- " << lhs_dydz_err << std::endl; 
0503     */
0504 
0505       double phidiff =
0506           (fabs(rhs_phi - lhs_phi) > 2 * M_PI ? fabs(rhs_phi - lhs_phi) - 2 * M_PI : fabs(rhs_phi - lhs_phi));
0507 
0508       if (phidiff < ClusterDPhi && fabs(rhs_theta - lhs_theta) < ClusterDTheta)  // phi overlap region is 37 degrees
0509         result = true;
0510     }  // loop over duplicate segments
0511   }  // same endcap,station,ring
0512 
0513   return result;
0514 }