MuonMesh.cc

Go to the documentation of this file.

#include "RecoMuon/MuonIdentification/interface/MuonMesh.h"
#include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
#include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
#include "Geometry/CSCGeometry/interface/CSCGeometry.h"

#include <utility>
#include <algorithm>

MuonMesh::MuonMesh(const edm::ParameterSet& parm) : doME1a(parm.getParameter<bool>("ME1a")),
                            doOverlaps(parm.getParameter<bool>("Overlap")),
                            doClustering(parm.getParameter<bool>("Clustering")),
                            OverlapDPhi(parm.getParameter<double>("OverlapDPhi")), 
                            OverlapDTheta(parm.getParameter<double>("OverlapDTheta")), 
                            ClusterDPhi(parm.getParameter<double>("ClusterDPhi")), 
                            ClusterDTheta(parm.getParameter<double>("ClusterDTheta")) {
}

void MuonMesh::fillMesh(std::vector<reco::Muon>* inputMuons) {

  for(std::vector<reco::Muon>::iterator muonIter1 = inputMuons->begin();
      muonIter1 != inputMuons->end();
      ++muonIter1) {
    if(muonIter1->isTrackerMuon()){

      mesh_[&*muonIter1]; // create this entry if it's a tracker muon
      for(std::vector<reco::Muon>::iterator muonIter2 = inputMuons->begin();
      muonIter2 != inputMuons->end();
      ++muonIter2) {
    if(muonIter2->isTrackerMuon()) {
      if(muonIter2 != muonIter1) {
        
        // now fill all the edges for muon1 based on overlaps with muon2
        for(std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = muonIter1->matches().begin();
        chamberIter1 != muonIter1->matches().end();
        ++chamberIter1) {         
          for(std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
          segmentIter1 != chamberIter1->segmentMatches.end();
          ++segmentIter1) {
        for(std::vector<reco::MuonChamberMatch>::iterator chamberIter2 = muonIter2->matches().begin();
            chamberIter2 != muonIter2->matches().end();
            ++chamberIter2) {         
          for(std::vector<reco::MuonSegmentMatch>::iterator segmentIter2 = chamberIter2->segmentMatches.begin();
              segmentIter2 != chamberIter2->segmentMatches.end();
              ++segmentIter2) {
                          

            //if(segmentIter1->mask & 0x1e0000 && segmentIter2->mask &0x1e0000) {

              bool addsegment(false);
              
              if( segmentIter1->cscSegmentRef.isNonnull() && segmentIter2->cscSegmentRef.isNonnull() ) {
                        
            if( doME1a && 
                isDuplicateOf(segmentIter1->cscSegmentRef,segmentIter2->cscSegmentRef) &&
                CSCDetId(chamberIter1->id).ring() == 4 && CSCDetId(chamberIter2->id).ring() == 4 &&
                chamberIter1->id == chamberIter2->id ) {
              addsegment = true;
              //std::cout << "\tME1/a sharing detected." << std::endl;
            }
            
            if( doOverlaps &&
                isDuplicateOf(std::make_pair(CSCDetId(chamberIter1->id),segmentIter1->cscSegmentRef),
                     std::make_pair(CSCDetId(chamberIter2->id),segmentIter2->cscSegmentRef)) ) {
              addsegment = true;
              //std::cout << "\tChamber Overlap sharing detected." << std::endl;
            }
            
            if( doClustering &&
                isClusteredWith(std::make_pair(CSCDetId(chamberIter1->id),segmentIter1->cscSegmentRef),
                       std::make_pair(CSCDetId(chamberIter2->id),segmentIter2->cscSegmentRef)) ) {
              addsegment = true;
              //std::cout << "\tCluster sharing detected." << std::endl;
            }
            //std::cout << std::endl;
              } // has valid csc segment ref
              
              if(addsegment) { // add segment if clusters/overlaps/replicant and doesn't already exist
            
            if(find(mesh_[&*muonIter1].begin(),
                mesh_[&*muonIter1].end(),
                std::make_pair(&*muonIter2,
                           std::make_pair(&*chamberIter2,
                                  &*segmentIter2)
                           )
                ) == mesh_[&*muonIter1].end()) {
              mesh_[&*muonIter1].push_back(std::make_pair(&*muonIter2,
                                      std::make_pair(&*chamberIter2,
                                             &*segmentIter2)
                                      )
                               );
            } // find
              } // add segment?
              //} // both segments won arbitration
          } // segmentIter 2                          
        } // chamberIter2
          } //segmentIter1
        } // chamberIter1   
        
      } // if different muon    
    } // is tracker muon  
      } // muonIter2
      
    } // is tracker muon
  } // muonIter1

  // special cases

  // one muon: mark all segments belonging to a muon as cleaned as there are no other muons to fight with
  if(mesh_.size() == 1) {
    for(std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = mesh_.begin()->first->matches().begin();
    chamberIter1 != mesh_.begin()->first->matches().end();
    ++chamberIter1) {         
      for(std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
      segmentIter1 != chamberIter1->segmentMatches.end();
      ++segmentIter1) {
    segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);    
      } // segmentIter1
    } // chamberIter1
  } // if only one tracker muon set winner bit boosted arbitration

  // segments that are not shared amongst muons and the have won all segment arbitration flags need to be promoted
  // also promote DT segments
  if(mesh_.size() > 1) {
    for( MeshType::iterator i = mesh_.begin(); i != mesh_.end(); ++i ) {
      for( std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = i->first->matches().begin();
       chamberIter1 != i->first->matches().end();
       ++chamberIter1 ) {   
    for(std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
        segmentIter1 != chamberIter1->segmentMatches.end();
        ++segmentIter1) {

      bool shared(false);
      
      for( AssociationType::iterator j = i->second.begin();
           j != i->second.end();
           ++j ) {

        if( segmentIter1->cscSegmentRef.isNonnull() && 
        j->second.second->cscSegmentRef.isNonnull() ) {               
          if(chamberIter1->id.subdetId() == MuonSubdetId::CSC &&
         j->second.first->id.subdetId() == MuonSubdetId::CSC ) {
        CSCDetId segIterId(chamberIter1->id), shareId(j->second.first->id);

        if( doOverlaps &&
            isDuplicateOf(std::make_pair(segIterId,segmentIter1->cscSegmentRef),
                  std::make_pair(shareId,j->second.second->cscSegmentRef)) )
          shared = true;
        
        if( doME1a && 
            isDuplicateOf(segmentIter1->cscSegmentRef,j->second.second->cscSegmentRef) &&
            segIterId.ring() == 4 && shareId.ring() == 4 &&
            segIterId == segIterId)
          shared = true;

        if( doClustering && 
            isClusteredWith(std::make_pair(CSCDetId(chamberIter1->id),segmentIter1->cscSegmentRef),
                    std::make_pair(CSCDetId(j->second.first->id),j->second.second->cscSegmentRef)) )
          shared = true;
          } // in CSCs?
        } // cscSegmentRef non null?
      } // j
      
      // Promote segments which have won all arbitration and are not shared or are DT segments
      if( ((segmentIter1->mask&0x1e0000) == 0x1e0000 && !shared) ||
          (chamberIter1->id.subdetId() == MuonSubdetId::DT && (segmentIter1->mask&0x1e0000)) ) 
          segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);

    } // segmentIter1
      } // chamberIter1
    }// i
  } // if non-trivial case

}

void MuonMesh::pruneMesh() {
  
  for( MeshType::iterator i = mesh_.begin(); i != mesh_.end(); ++i ) {

        for( AssociationType::iterator j = i->second.begin();
     j != i->second.end();
     ++j ) {

      for( std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = i->first->matches().begin();
       chamberIter1 != i->first->matches().end();
       ++chamberIter1 ) {   
    for(std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
        segmentIter1 != chamberIter1->segmentMatches.end();
        ++segmentIter1) {

      if(j->second.second->cscSegmentRef.isNonnull() && segmentIter1->cscSegmentRef.isNonnull()) {
        
        //UNUSED:       bool me1a(false), overlap(false), cluster(false);

        // remove physical overlap duplicates first
        if( doOverlaps &&
        isDuplicateOf(std::make_pair(CSCDetId(chamberIter1->id),segmentIter1->cscSegmentRef),
                  std::make_pair(CSCDetId(j->second.first->id),j->second.second->cscSegmentRef)) ) {

          if ( i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) >
           j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ) {
        
        segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByOvlClean);
        segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
        
        //UNUSED:       overlap = true;
          } else if ( i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ==
              j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ) { // muon with more matched stations wins
        
        if((segmentIter1->mask & 0x1e0000) > (j->second.second->mask & 0x1e0000)) { // segment with better match wins 
          
          segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByOvlClean);
          segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
          
          //UNUSED:       overlap = true;                  
        } else { // ??
          // leave this available for later
        }        
          } // overlap duplicate resolution
        } // is overlap duplicate
        
        // do ME1/a arbitration second since the tie breaker depends on other stations
        // Unlike the other cleanings this one removes the bits from segments associated to tracks which
        // fail cleaning. (Instead of setting bits for the segments which win.)
        if( doME1a && 
        isDuplicateOf(segmentIter1->cscSegmentRef,j->second.second->cscSegmentRef) &&
        CSCDetId(chamberIter1->id).ring() == 4 && CSCDetId(j->second.first->id).ring() == 4 &&
        chamberIter1->id ==  j->second.first->id ) {          

          if( j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) < 
          i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ) {
                
        for(AssociationType::iterator AsscIter1 = i->second.begin();
            AsscIter1 != i->second.end();
            ++AsscIter1) {
          if(AsscIter1->second.second->cscSegmentRef.isNonnull())
            if(j->first == AsscIter1->first &&
               j->second.first == AsscIter1->second.first &&
               isDuplicateOf(segmentIter1->cscSegmentRef,AsscIter1->second.second->cscSegmentRef)) {
              AsscIter1->second.second->mask &= ~reco::MuonSegmentMatch::BelongsToTrackByME1aClean;
            }       
        }  

        //UNUSED:       me1a = true;
          } else if ( j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) == 
              i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ) { // muon with best arbitration wins
                        
        bool bestArb(true);
        
        for(AssociationType::iterator AsscIter1 = i->second.begin();
            AsscIter1 != i->second.end();
            ++AsscIter1) {
          if(AsscIter1->second.second->cscSegmentRef.isNonnull())
            if(j->first == AsscIter1->first &&
               j->second.first == AsscIter1->second.first &&
               isDuplicateOf(segmentIter1->cscSegmentRef,AsscIter1->second.second->cscSegmentRef) && 
               (segmentIter1->mask & 0x1e0000) < (AsscIter1->second.second->mask & 0x1e0000))
              bestArb = false;
        }
        
        if(bestArb) {
          
          for(AssociationType::iterator AsscIter1 = i->second.begin();
              AsscIter1 != i->second.end();
              ++AsscIter1) {
            if(AsscIter1->second.second->cscSegmentRef.isNonnull())
              if(j->first == AsscIter1->first &&
             j->second.first == AsscIter1->second.first &&
             isDuplicateOf(segmentIter1->cscSegmentRef,AsscIter1->second.second->cscSegmentRef)) {
            AsscIter1->second.second->mask &= ~reco::MuonSegmentMatch::BelongsToTrackByME1aClean;
              }
          }
          
        }    
        //UNUSED        me1a = true;                

          } // ME1/a duplicate resolution                 
        } // is ME1/aduplicate?
        
        if(doClustering && 
           isClusteredWith(std::make_pair(CSCDetId(chamberIter1->id),segmentIter1->cscSegmentRef),
                   std::make_pair(CSCDetId(j->second.first->id),j->second.second->cscSegmentRef))) {

          if (i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) >
          j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ) {

         segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
         segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);

         //UNUSED:       cluster = true;
          } else if (i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) <
             j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000)) {

        j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
        j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
        
        //UNUSED:       cluster = true;
          } else { // muon with more matched stations wins
         
         if((segmentIter1->mask & 0x1e0000) > (j->second.second->mask & 0x1e0000)) { // segment with better match wins 
           
           segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
           segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
           
           //UNUSED:           cluster = true;                 
         } else if ((segmentIter1->mask & 0x1e0000) < (j->second.second->mask & 0x1e0000)){ //
           
           j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
           j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
           
           //UNUSED:           cluster = true;                 
         } else {
         }
           } // cluster sharing resolution
          
        } // is clustered with?
        
      } // csc ref nonnull    
    } // segmentIter1
      } // chamberIter1       
    } // j, associated segments iterator
  } // i, map iterator

  // final step: make sure everything that's won a cleaning flag has the "BelongsToTrackByCleaning" flag

   for( MeshType::iterator i = mesh_.begin(); i != mesh_.end(); ++i ) {
     for( std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = i->first->matches().begin();
      chamberIter1 != i->first->matches().end();
      ++chamberIter1 ) {    
       for(std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
       segmentIter1 != chamberIter1->segmentMatches.end();
       ++segmentIter1) {
     // set cleaning bit if initial no cleaning bit but there are cleaning algorithm bits set.
     if( !segmentIter1->isMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning) &&
         segmentIter1->isMask(0xe00000) )
       segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);  
       }// segmentIter1
     } // chamberIter1
   } // i
  
}

bool MuonMesh::isDuplicateOf(const CSCSegmentRef& lhs, const CSCSegmentRef& rhs) const // this isDuplicateOf() deals with duplicate segments in ME1/a
{
  bool result(false);
 
  if(!lhs->isME11a_duplicate())
    return result;

  std::vector<CSCSegment> lhs_duplicates = lhs->duplicateSegments();
      
  if(fabs(lhs->localPosition().x()        - rhs->localPosition().x()      ) < 1E-3 &&
     fabs(lhs->localPosition().y()        - rhs->localPosition().y()      ) < 1E-3 &&
     fabs(lhs->localDirection().x()/lhs->localDirection().z()    - rhs->localDirection().x()/rhs->localDirection().z()   ) < 1E-3 &&
     fabs(lhs->localDirection().y()/lhs->localDirection().z()    - rhs->localDirection().y()/rhs->localDirection().z()   ) < 1E-3 &&
     fabs(lhs->localPositionError().xx()  - rhs->localPositionError().xx()   ) < 1E-3 &&
     fabs(lhs->localPositionError().yy()  - rhs->localPositionError().yy()   ) < 1E-3 &&
     fabs(lhs->localDirectionError().xx() - rhs->localDirectionError().xx()) < 1E-3 &&
     fabs(lhs->localDirectionError().yy() - rhs->localDirectionError().yy()) < 1E-3)
    result = true;

  for( std::vector<CSCSegment>::const_iterator segIter1 = lhs_duplicates.begin();
       segIter1 != lhs_duplicates.end();
       ++segIter1 ) { // loop over lhs duplicates

    if(fabs(segIter1->localPosition().x()        - rhs->localPosition().x()      ) < 1E-3 &&
       fabs(segIter1->localPosition().y()        - rhs->localPosition().y()      ) < 1E-3 &&
       fabs(segIter1->localDirection().x()/segIter1->localDirection().z()    - rhs->localDirection().x()/rhs->localDirection().z()   ) < 1E-3 &&
       fabs(segIter1->localDirection().y()/segIter1->localDirection().z()    - rhs->localDirection().y()/rhs->localDirection().z()   ) < 1E-3 &&
       fabs(segIter1->localPositionError().xx()  - rhs->localPositionError().xx()   ) < 1E-3 &&
       fabs(segIter1->localPositionError().yy()  - rhs->localPositionError().yy()   ) < 1E-3 &&
       fabs(segIter1->localDirectionError().xx() - rhs->localDirectionError().xx()) < 1E-3 &&
       fabs(segIter1->localDirectionError().yy() - rhs->localDirectionError().yy()) < 1E-3)
      result = true;
    /*
    if(fabs(segIter1->localPosition().x()        - rhs->localPosition().x()      ) < 2*sqrt(segIter1->localPositionError().xx()) &&
       fabs(segIter1->localPosition().y()        - rhs->localPosition().y()      ) < 2*sqrt(segIter1->localPositionError().yy()) &&
       fabs(segIter1->localDirection().x()/segIter1->localDirection().z()    - rhs->localDirection().x()/rhs->localDirection().z()   ) 
       < 2*std::sqrt(std::max(segIter1->localDirectionError().yy(),rhs->localDirectionError().xx())) &&
       fabs(segIter1->localDirection().y()/segIter1->localDirection().z()    - rhs->localDirection().y()/rhs->localDirection().z()   ) 
       < 2*std::sqrt(std::max(segIter1->localDirectionError().yy(),rhs->localDirectionError().yy())))
      result = true;
    */

  } // loop over duplicates

  return result;
}

bool MuonMesh::isDuplicateOf(const std::pair<CSCDetId,CSCSegmentRef>& rhs, 
                 const std::pair<CSCDetId,CSCSegmentRef>& lhs) const // this isDuplicateOf() deals with "overlapping chambers" duplicates
{
  bool result(false);
  
  // try to implement the simple case first just back-to-back segments without treatment of ME1/a ganging
  // ME1a should be a simple extension of this

  if(rhs.first.endcap() == lhs.first.endcap() &&
     rhs.first.station() == lhs.first.station() &&
     rhs.first.ring() == lhs.first.ring()) { // if same endcap,station,ring (minimal requirement for ovl candidate)
    /*
    std::cout << "Chamber 1: " << rhs.first << std::endl 
          << "Chamber 2: " << lhs.first << std::endl;

    std::cout << "Same endcap,ring,station." << std::endl;    
    */
    //create neighboring chamber labels, treat ring as (Z mod 36 or 18) + 1 number line: left, right defined accordingly.
    unsigned modulus = ((rhs.first.ring() != 1 || rhs.first.station() == 1) ? 36 : 18);
    int left_neighbor = (((rhs.first.chamber() - 1 + modulus)%modulus == 0 ) ? modulus : (rhs.first.chamber() - 1 + modulus)%modulus ); // + modulus to ensure positivity
    int right_neighbor = (((rhs.first.chamber() + 1)%modulus == 0 ) ? modulus : (rhs.first.chamber() + 1)%modulus );

    if(lhs.first.chamber() == left_neighbor || 
       lhs.first.chamber() == right_neighbor) { // if this is a neighboring chamber then it can be an overlap
      
      std::vector<CSCSegment> thesegments;
      thesegments.push_back(*(lhs.second));
      /*
      if(lhs.second->isME11a_duplicate())
    thesegments.insert(thesegments.begin(),
               lhs.second->duplicateSegments().begin(),
               lhs.second->duplicateSegments().end());
      */

      //std::cout << "lhs is in neighoring chamber of rhs." << std::endl;

      // rhs local direction info
      /*
      double rhs_dydz = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).y()/
                    geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).z();
      double rhs_dxdz = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).x()/
                    geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).z();
      double rhs_dydz_err = rhs.second->localDirectionError().yy();
      double rhs_dxdz_err = rhs.second->localDirectionError().xx();
      */
      
      //rhs global position info
      double rhs_phi = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).phi();
      double rhs_theta = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).theta();
      double rhs_z = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).z();

      for(std::vector<CSCSegment>::const_iterator ilhs = thesegments.begin(); ilhs != thesegments.end(); ++ilhs) {
    
    // lhs local direction info
    /*
    double lhs_dydz = geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).y()/
                      geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).z();
    double lhs_dxdz = geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).x()/
                      geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).z();
    double lhs_dydz_err = ilhs->localDirectionError().yy();
    double lhs_dxdz_err = ilhs->localDirectionError().xx();
    */
    
    //lhs global position info
    double lhs_phi = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).phi();
    double lhs_theta = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).theta();
    double lhs_z = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).z();
    /*
      std::cout << "RHS Segment Parameters:" << std::endl
      << "\t RHS Phi   : " << rhs_phi << std::endl
      << "\t RHS Theta : " << rhs_theta << std::endl
      << "\t RHS dx/dz : " << rhs_dxdz << " +- " << rhs_dxdz_err << std::endl
      << "\t RHS dy/dz : " << rhs_dydz << " +- " << rhs_dydz_err << std::endl; 
      
      std::cout << "LHS Segment Parameters:" << std::endl
      << "\t LHS Phi   : " << lhs_phi << std::endl
      << "\t LHS Theta : " << lhs_theta << std::endl
      << "\t LHS dx/dz : " << lhs_dxdz << " +- " << lhs_dxdz_err << std::endl
      << "\t LHS dy/dz : " << lhs_dydz << " +- " << lhs_dydz_err << std::endl; 
    */
    
    double phidiff = (fabs(rhs_phi - lhs_phi) > 2*M_PI ? fabs(rhs_phi - lhs_phi) - 2*M_PI : fabs(rhs_phi - lhs_phi));

    if(phidiff < OverlapDPhi && fabs(rhs_theta - lhs_theta) < OverlapDTheta && 
       fabs(rhs_z) < fabs(lhs_z) && rhs_z*lhs_z > 0) // phi overlap region is 3.5 degrees and rhs is infront of lhs
      result = true;
      } // loop over duplicate segments
    }// neighboring chamber
  } // same endcap,station,ring

  return result;
}

bool MuonMesh::isClusteredWith(const std::pair<CSCDetId,CSCSegmentRef>& lhs, 
                   const std::pair<CSCDetId,CSCSegmentRef>& rhs) const
{
  bool result(false);
  
  // try to implement the simple case first just back-to-back segments without treatment of ME1/a ganging
  // ME1a should be a simple extension of this

  //std::cout << lhs.first << ' ' << rhs.first << std::endl;

  if(rhs.first.endcap() == lhs.first.endcap() && lhs.first.station() < rhs.first.station()) {
          
      std::vector<CSCSegment> thesegments;
      thesegments.push_back(*(lhs.second));
      /*
      if(lhs.second->isME11a_duplicate())
    thesegments.insert(thesegments.begin(),
               lhs.second->duplicateSegments().begin(),
               lhs.second->duplicateSegments().end());
      */
      //std::cout << "lhs is in neighoring chamber of rhs." << std::endl;

      // rhs local direction info
      /*
      double rhs_dydz = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).y()/
                    geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).z();
      double rhs_dxdz = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).x()/
                    geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).z();
      double rhs_dydz_err = rhs.second->localDirectionError().yy();
      double rhs_dxdz_err = rhs.second->localDirectionError().xx();
      */
      
      //rhs global position info
      double rhs_phi = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).phi();
      double rhs_theta = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).theta();

      for(std::vector<CSCSegment>::const_iterator ilhs = thesegments.begin(); ilhs != thesegments.end(); ++ilhs) {
    
    // lhs local direction info
    /*
    double lhs_dydz = geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).y()/
                      geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).z();
    double lhs_dxdz = geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).x()/
                      geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).z();
    double lhs_dydz_err = ilhs->localDirectionError().yy();
    double lhs_dxdz_err = ilhs->localDirectionError().xx();
    */

    //lhs global position info
    double lhs_phi = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).phi();
    double lhs_theta = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).theta();
    /*
      std::cout << "RHS Segment Parameters:" << std::endl
      << "\t RHS Phi   : " << rhs_phi << std::endl
      << "\t RHS Theta : " << rhs_theta << std::endl
      << "\t RHS dx/dz : " << rhs_dxdz << " +- " << rhs_dxdz_err << std::endl
      << "\t RHS dy/dz : " << rhs_dydz << " +- " << rhs_dydz_err << std::endl; 
      
      std::cout << "LHS Segment Parameters:" << std::endl
      << "\t LHS Phi   : " << lhs_phi << std::endl
      << "\t LHS Theta : " << lhs_theta << std::endl
      << "\t LHS dx/dz : " << lhs_dxdz << " +- " << lhs_dxdz_err << std::endl
      << "\t LHS dy/dz : " << lhs_dydz << " +- " << lhs_dydz_err << std::endl; 
    */

    double phidiff = (fabs(rhs_phi - lhs_phi) > 2*M_PI ? fabs(rhs_phi - lhs_phi) - 2*M_PI : fabs(rhs_phi - lhs_phi));
    
    if(phidiff < ClusterDPhi && fabs(rhs_theta - lhs_theta) < ClusterDTheta) // phi overlap region is 37 degrees
      result = true;
      } // loop over duplicate segments    
  } // same endcap,station,ring

  return result;
}