Reference Manual

Go to the documentation of this file.
00001 #include "RecoMuon/MuonIdentification/interface/MuonMesh.h"
00002 #include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
00003 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
00004 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
00005 
00006 #include <utility>
00007 #include <algorithm>
00008 
00009 MuonMesh::MuonMesh(const edm::ParameterSet& parm) : doME1a(parm.getParameter<bool>("ME1a")),
00010                                                     doOverlaps(parm.getParameter<bool>("Overlap")),
00011                                                     doClustering(parm.getParameter<bool>("Clustering")),
00012                                                     OverlapDPhi(parm.getParameter<double>("OverlapDPhi")), 
00013                                                     OverlapDTheta(parm.getParameter<double>("OverlapDTheta")), 
00014                                                     ClusterDPhi(parm.getParameter<double>("ClusterDPhi")), 
00015                                                     ClusterDTheta(parm.getParameter<double>("ClusterDTheta")) {
00016 }
00017 
00018 void MuonMesh::fillMesh(std::vector<reco::Muon>* inputMuons) {
00019 
00020   for(std::vector<reco::Muon>::iterator muonIter1 = inputMuons->begin();
00021       muonIter1 != inputMuons->end();
00022       ++muonIter1) {
00023     if(muonIter1->isTrackerMuon()){
00024 
00025       mesh_[&*muonIter1]; // create this entry if it's a tracker muon
00026       for(std::vector<reco::Muon>::iterator muonIter2 = inputMuons->begin();
00027           muonIter2 != inputMuons->end();
00028           ++muonIter2) {
00029         if(muonIter2->isTrackerMuon()) {
00030           if(muonIter2 != muonIter1) {
00031             
00032             // now fill all the edges for muon1 based on overlaps with muon2
00033             for(std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = muonIter1->matches().begin();
00034                 chamberIter1 != muonIter1->matches().end();
00035                 ++chamberIter1) {             
00036               for(std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
00037                   segmentIter1 != chamberIter1->segmentMatches.end();
00038                   ++segmentIter1) {
00039                 for(std::vector<reco::MuonChamberMatch>::iterator chamberIter2 = muonIter2->matches().begin();
00040                     chamberIter2 != muonIter2->matches().end();
00041                     ++chamberIter2) {             
00042                   for(std::vector<reco::MuonSegmentMatch>::iterator segmentIter2 = chamberIter2->segmentMatches.begin();
00043                       segmentIter2 != chamberIter2->segmentMatches.end();
00044                       ++segmentIter2) {
00045                                       
00046 
00047                     //if(segmentIter1->mask & 0x1e0000 && segmentIter2->mask &0x1e0000) {
00048 
00049                       bool addsegment(false);
00050                       
00051                       if( segmentIter1->cscSegmentRef.isNonnull() && segmentIter2->cscSegmentRef.isNonnull() ) {
00052                                                 
00053                         if( doME1a && 
00054                             isDuplicateOf(segmentIter1->cscSegmentRef,segmentIter2->cscSegmentRef) &&
00055                             CSCDetId(chamberIter1->id).ring() == 4 && CSCDetId(chamberIter2->id).ring() == 4 &&
00056                             chamberIter1->id == chamberIter2->id ) {
00057                           addsegment = true;
00058                           //std::cout << "\tME1/a sharing detected." << std::endl;
00059                         }
00060                         
00061                         if( doOverlaps &&
00062                             isDuplicateOf(std::make_pair(CSCDetId(chamberIter1->id),segmentIter1->cscSegmentRef),
00063                                          std::make_pair(CSCDetId(chamberIter2->id),segmentIter2->cscSegmentRef)) ) {
00064                           addsegment = true;
00065                           //std::cout << "\tChamber Overlap sharing detected." << std::endl;
00066                         }
00067                         
00068                         if( doClustering &&
00069                             isClusteredWith(std::make_pair(CSCDetId(chamberIter1->id),segmentIter1->cscSegmentRef),
00070                                            std::make_pair(CSCDetId(chamberIter2->id),segmentIter2->cscSegmentRef)) ) {
00071                           addsegment = true;
00072                           //std::cout << "\tCluster sharing detected." << std::endl;
00073                         }
00074                         //std::cout << std::endl;
00075                       } // has valid csc segment ref
00076                       
00077                       if(addsegment) { // add segment if clusters/overlaps/replicant and doesn't already exist
00078                         
00079                         if(find(mesh_[&*muonIter1].begin(),
00080                                 mesh_[&*muonIter1].end(),
00081                                 std::make_pair(&*muonIter2,
00082                                                std::make_pair(&*chamberIter2,
00083                                                               &*segmentIter2)
00084                                                )
00085                                 ) == mesh_[&*muonIter1].end()) {
00086                           mesh_[&*muonIter1].push_back(std::make_pair(&*muonIter2,
00087                                                                       std::make_pair(&*chamberIter2,
00088                                                                                      &*segmentIter2)
00089                                                                       )
00090                                                        );
00091                         } // find
00092                       } // add segment?
00093                       //} // both segments won arbitration
00094                   } // segmentIter 2                                              
00095                 } // chamberIter2
00096               } //segmentIter1
00097             } // chamberIter1   
00098             
00099           } // if different muon        
00100         } // is tracker muon  
00101       } // muonIter2
00102       
00103     } // is tracker muon
00104   } // muonIter1
00105 
00106   // special cases
00107 
00108   // one muon: mark all segments belonging to a muon as cleaned as there are no other muons to fight with
00109   if(mesh_.size() == 1) {
00110     for(std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = mesh_.begin()->first->matches().begin();
00111         chamberIter1 != mesh_.begin()->first->matches().end();
00112         ++chamberIter1) {             
00113       for(std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
00114           segmentIter1 != chamberIter1->segmentMatches.end();
00115           ++segmentIter1) {
00116         segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);        
00117       } // segmentIter1
00118     } // chamberIter1
00119   } // if only one tracker muon set winner bit boosted arbitration
00120 
00121   // segments that are not shared amongst muons and the have won all segment arbitration flags need to be promoted
00122   // also promote DT segments
00123   if(mesh_.size() > 1) {
00124     for( MeshType::iterator i = mesh_.begin(); i != mesh_.end(); ++i ) {
00125       for( std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = i->first->matches().begin();
00126            chamberIter1 != i->first->matches().end();
00127            ++chamberIter1 ) {   
00128         for(std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
00129             segmentIter1 != chamberIter1->segmentMatches.end();
00130             ++segmentIter1) {
00131 
00132           bool shared(false);
00133           
00134           for( AssociationType::iterator j = i->second.begin();
00135                j != i->second.end();
00136                ++j ) {
00137 
00138             if( segmentIter1->cscSegmentRef.isNonnull() && 
00139                 j->second.second->cscSegmentRef.isNonnull() ) {               
00140               if(chamberIter1->id.subdetId() == MuonSubdetId::CSC &&
00141                  j->second.first->id.subdetId() == MuonSubdetId::CSC ) {
00142                 CSCDetId segIterId(chamberIter1->id), shareId(j->second.first->id);
00143 
00144                 if( doOverlaps &&
00145                     isDuplicateOf(std::make_pair(segIterId,segmentIter1->cscSegmentRef),
00146                                   std::make_pair(shareId,j->second.second->cscSegmentRef)) )
00147                   shared = true;
00148                 
00149                 if( doME1a && 
00150                     isDuplicateOf(segmentIter1->cscSegmentRef,j->second.second->cscSegmentRef) &&
00151                     segIterId.ring() == 4 && shareId.ring() == 4 &&
00152                     segIterId == segIterId)
00153                   shared = true;
00154 
00155                 if( doClustering && 
00156                     isClusteredWith(std::make_pair(CSCDetId(chamberIter1->id),segmentIter1->cscSegmentRef),
00157                                     std::make_pair(CSCDetId(j->second.first->id),j->second.second->cscSegmentRef)) )
00158                   shared = true;
00159               } // in CSCs?
00160             } // cscSegmentRef non null?
00161           } // j
00162           
00163           // Promote segments which have won all arbitration and are not shared or are DT segments
00164           if( ((segmentIter1->mask&0x1e0000) == 0x1e0000 && !shared) ||
00165               (chamberIter1->id.subdetId() == MuonSubdetId::DT && (segmentIter1->mask&0x1e0000)) ) 
00166               segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
00167 
00168         } // segmentIter1
00169       } // chamberIter1
00170     }// i
00171   } // if non-trivial case
00172 
00173 }
00174 
00175 void MuonMesh::pruneMesh() {
00176   
00177   for( MeshType::iterator i = mesh_.begin(); i != mesh_.end(); ++i ) {
00178 
00179         for( AssociationType::iterator j = i->second.begin();
00180          j != i->second.end();
00181          ++j ) {
00182 
00183       for( std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = i->first->matches().begin();
00184            chamberIter1 != i->first->matches().end();
00185            ++chamberIter1 ) {   
00186         for(std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
00187             segmentIter1 != chamberIter1->segmentMatches.end();
00188             ++segmentIter1) {
00189 
00190           if(j->second.second->cscSegmentRef.isNonnull() && segmentIter1->cscSegmentRef.isNonnull()) {
00191             
00192             bool me1a(false), overlap(false), cluster(false);
00193 
00194             // remove physical overlap duplicates first
00195             if( doOverlaps &&
00196                 isDuplicateOf(std::make_pair(CSCDetId(chamberIter1->id),segmentIter1->cscSegmentRef),
00197                               std::make_pair(CSCDetId(j->second.first->id),j->second.second->cscSegmentRef)) ) {
00198 
00199               if ( i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) >
00200                    j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ) {
00201                 
00202                 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByOvlClean);
00203                 segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
00204                 
00205                 overlap = true;
00206               } else if ( i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ==
00207                           j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ) { // muon with more matched stations wins
00208                 
00209                 if((segmentIter1->mask & 0x1e0000) > (j->second.second->mask & 0x1e0000)) { // segment with better match wins 
00210                   
00211                   segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByOvlClean);
00212                   segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
00213                   
00214                   overlap = true;                                  
00215                 } else { // ??
00216                   // leave this available for later
00217                 }                
00218               } // overlap duplicate resolution
00219             } // is overlap duplicate
00220             
00221             // do ME1/a arbitration second since the tie breaker depends on other stations
00222             // Unlike the other cleanings this one removes the bits from segments associated to tracks which
00223             // fail cleaning. (Instead of setting bits for the segments which win.)
00224             if( doME1a && 
00225                 isDuplicateOf(segmentIter1->cscSegmentRef,j->second.second->cscSegmentRef) &&
00226                 CSCDetId(chamberIter1->id).ring() == 4 && CSCDetId(j->second.first->id).ring() == 4 &&
00227                 chamberIter1->id ==  j->second.first->id ) {          
00228 
00229               if( j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) < 
00230                   i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ) {
00231                                 
00232                 for(AssociationType::iterator AsscIter1 = i->second.begin();
00233                     AsscIter1 != i->second.end();
00234                     ++AsscIter1) {
00235                   if(AsscIter1->second.second->cscSegmentRef.isNonnull())
00236                     if(j->first == AsscIter1->first &&
00237                        j->second.first == AsscIter1->second.first &&
00238                        isDuplicateOf(segmentIter1->cscSegmentRef,AsscIter1->second.second->cscSegmentRef)) {
00239                       AsscIter1->second.second->mask &= ~reco::MuonSegmentMatch::BelongsToTrackByME1aClean;
00240                     }           
00241                 }  
00242 
00243                 me1a = true;
00244               } else if ( j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) == 
00245                           i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ) { // muon with best arbitration wins
00246                                                 
00247                 bool bestArb(true);
00248                 
00249                 for(AssociationType::iterator AsscIter1 = i->second.begin();
00250                     AsscIter1 != i->second.end();
00251                     ++AsscIter1) {
00252                   if(AsscIter1->second.second->cscSegmentRef.isNonnull())
00253                     if(j->first == AsscIter1->first &&
00254                        j->second.first == AsscIter1->second.first &&
00255                        isDuplicateOf(segmentIter1->cscSegmentRef,AsscIter1->second.second->cscSegmentRef) && 
00256                        (segmentIter1->mask & 0x1e0000) < (AsscIter1->second.second->mask & 0x1e0000))
00257                       bestArb = false;
00258                 }
00259                 
00260                 if(bestArb) {
00261                   
00262                   for(AssociationType::iterator AsscIter1 = i->second.begin();
00263                       AsscIter1 != i->second.end();
00264                       ++AsscIter1) {
00265                     if(AsscIter1->second.second->cscSegmentRef.isNonnull())
00266                       if(j->first == AsscIter1->first &&
00267                          j->second.first == AsscIter1->second.first &&
00268                          isDuplicateOf(segmentIter1->cscSegmentRef,AsscIter1->second.second->cscSegmentRef)) {
00269                         AsscIter1->second.second->mask &= ~reco::MuonSegmentMatch::BelongsToTrackByME1aClean;
00270                       }
00271                   }
00272                   
00273                 }        
00274                 me1a = true;                            
00275 
00276               } // ME1/a duplicate resolution                 
00277             } // is ME1/aduplicate?
00278             
00279             if(doClustering && 
00280                isClusteredWith(std::make_pair(CSCDetId(chamberIter1->id),segmentIter1->cscSegmentRef),
00281                                std::make_pair(CSCDetId(j->second.first->id),j->second.second->cscSegmentRef))) {
00282 
00283               if (i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) >
00284                   j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) ) {
00285 
00286                  segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
00287                  segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
00288 
00289                  cluster = true;
00290               } else if (i->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000) <
00291                          j->first->numberOfMatches((reco::Muon::ArbitrationType)0x1e0000)) {
00292 
00293                 j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
00294                 j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
00295                 
00296                 cluster = true;
00297               } else { // muon with more matched stations wins
00298                  
00299                  if((segmentIter1->mask & 0x1e0000) > (j->second.second->mask & 0x1e0000)) { // segment with better match wins 
00300                    
00301                    segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
00302                    segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
00303                    
00304                    cluster = true;                                 
00305                  } else if ((segmentIter1->mask & 0x1e0000) < (j->second.second->mask & 0x1e0000)){ //
00306                    
00307                    j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByClusClean);
00308                    j->second.second->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);
00309                    
00310                    cluster = true;                                 
00311                  } else {
00312                  }
00313                } // cluster sharing resolution
00314               
00315             } // is clustered with?
00316             
00317           } // csc ref nonnull    
00318         } // segmentIter1
00319       } // chamberIter1       
00320     } // j, associated segments iterator
00321   } // i, map iterator
00322 
00323   // final step: make sure everything that's won a cleaning flag has the "BelongsToTrackByCleaning" flag
00324 
00325    for( MeshType::iterator i = mesh_.begin(); i != mesh_.end(); ++i ) {
00326      for( std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = i->first->matches().begin();
00327           chamberIter1 != i->first->matches().end();
00328           ++chamberIter1 ) {    
00329        for(std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
00330            segmentIter1 != chamberIter1->segmentMatches.end();
00331            ++segmentIter1) {
00332          // set cleaning bit if initial no cleaning bit but there are cleaning algorithm bits set.
00333          if( !segmentIter1->isMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning) &&
00334              segmentIter1->isMask(0xe00000) )
00335            segmentIter1->setMask(reco::MuonSegmentMatch::BelongsToTrackByCleaning);      
00336        }// segmentIter1
00337      } // chamberIter1
00338    } // i
00339   
00340 }
00341 
00342 bool MuonMesh::isDuplicateOf(const CSCSegmentRef& lhs, const CSCSegmentRef& rhs) const // this isDuplicateOf() deals with duplicate segments in ME1/a
00343 {
00344   bool result(false);
00345  
00346   if(!lhs->isME11a_duplicate())
00347     return result;
00348 
00349   std::vector<CSCSegment> lhs_duplicates = lhs->duplicateSegments();
00350       
00351   if(fabs(lhs->localPosition().x()        - rhs->localPosition().x()      ) < 1E-3 &&
00352      fabs(lhs->localPosition().y()        - rhs->localPosition().y()      ) < 1E-3 &&
00353      fabs(lhs->localDirection().x()/lhs->localDirection().z()    - rhs->localDirection().x()/rhs->localDirection().z()   ) < 1E-3 &&
00354      fabs(lhs->localDirection().y()/lhs->localDirection().z()    - rhs->localDirection().y()/rhs->localDirection().z()   ) < 1E-3 &&
00355      fabs(lhs->localPositionError().xx()  - rhs->localPositionError().xx()   ) < 1E-3 &&
00356      fabs(lhs->localPositionError().yy()  - rhs->localPositionError().yy()   ) < 1E-3 &&
00357      fabs(lhs->localDirectionError().xx() - rhs->localDirectionError().xx()) < 1E-3 &&
00358      fabs(lhs->localDirectionError().yy() - rhs->localDirectionError().yy()) < 1E-3)
00359     result = true;
00360 
00361   for( std::vector<CSCSegment>::const_iterator segIter1 = lhs_duplicates.begin();
00362        segIter1 != lhs_duplicates.end();
00363        ++segIter1 ) { // loop over lhs duplicates
00364 
00365     if(fabs(segIter1->localPosition().x()        - rhs->localPosition().x()      ) < 1E-3 &&
00366        fabs(segIter1->localPosition().y()        - rhs->localPosition().y()      ) < 1E-3 &&
00367        fabs(segIter1->localDirection().x()/segIter1->localDirection().z()    - rhs->localDirection().x()/rhs->localDirection().z()   ) < 1E-3 &&
00368        fabs(segIter1->localDirection().y()/segIter1->localDirection().z()    - rhs->localDirection().y()/rhs->localDirection().z()   ) < 1E-3 &&
00369        fabs(segIter1->localPositionError().xx()  - rhs->localPositionError().xx()   ) < 1E-3 &&
00370        fabs(segIter1->localPositionError().yy()  - rhs->localPositionError().yy()   ) < 1E-3 &&
00371        fabs(segIter1->localDirectionError().xx() - rhs->localDirectionError().xx()) < 1E-3 &&
00372        fabs(segIter1->localDirectionError().yy() - rhs->localDirectionError().yy()) < 1E-3)
00373       result = true;
00374     /*
00375     if(fabs(segIter1->localPosition().x()        - rhs->localPosition().x()      ) < 2*sqrt(segIter1->localPositionError().xx()) &&
00376        fabs(segIter1->localPosition().y()        - rhs->localPosition().y()      ) < 2*sqrt(segIter1->localPositionError().yy()) &&
00377        fabs(segIter1->localDirection().x()/segIter1->localDirection().z()    - rhs->localDirection().x()/rhs->localDirection().z()   ) 
00378        < 2*std::sqrt(std::max(segIter1->localDirectionError().yy(),rhs->localDirectionError().xx())) &&
00379        fabs(segIter1->localDirection().y()/segIter1->localDirection().z()    - rhs->localDirection().y()/rhs->localDirection().z()   ) 
00380        < 2*std::sqrt(std::max(segIter1->localDirectionError().yy(),rhs->localDirectionError().yy())))
00381       result = true;
00382     */
00383 
00384   } // loop over duplicates
00385 
00386   return result;
00387 }
00388 
00389 bool MuonMesh::isDuplicateOf(const std::pair<CSCDetId,CSCSegmentRef>& rhs, 
00390                              const std::pair<CSCDetId,CSCSegmentRef>& lhs) const // this isDuplicateOf() deals with "overlapping chambers" duplicates
00391 {
00392   bool result(false);
00393   
00394   // try to implement the simple case first just back-to-back segments without treatment of ME1/a ganging
00395   // ME1a should be a simple extension of this
00396 
00397   if(rhs.first.endcap() == lhs.first.endcap() &&
00398      rhs.first.station() == lhs.first.station() &&
00399      rhs.first.ring() == lhs.first.ring()) { // if same endcap,station,ring (minimal requirement for ovl candidate)
00400     /*
00401     std::cout << "Chamber 1: " << rhs.first << std::endl 
00402               << "Chamber 2: " << lhs.first << std::endl;
00403 
00404     std::cout << "Same endcap,ring,station." << std::endl;    
00405     */
00406     //create neighboring chamber labels, treat ring as (Z mod 36 or 18) + 1 number line: left, right defined accordingly.
00407     unsigned modulus = ((rhs.first.ring() != 1 || rhs.first.station() == 1) ? 36 : 18);
00408     int left_neighbor = (((rhs.first.chamber() - 1 + modulus)%modulus == 0 ) ? modulus : (rhs.first.chamber() - 1 + modulus)%modulus ); // + modulus to ensure positivity
00409     int right_neighbor = (((rhs.first.chamber() + 1)%modulus == 0 ) ? modulus : (rhs.first.chamber() + 1)%modulus );
00410 
00411     if(lhs.first.chamber() == left_neighbor || 
00412        lhs.first.chamber() == right_neighbor) { // if this is a neighboring chamber then it can be an overlap
00413       
00414       std::vector<CSCSegment> thesegments;
00415       thesegments.push_back(*(lhs.second));
00416       /*
00417       if(lhs.second->isME11a_duplicate())
00418         thesegments.insert(thesegments.begin(),
00419                            lhs.second->duplicateSegments().begin(),
00420                            lhs.second->duplicateSegments().end());
00421       */
00422 
00423       //std::cout << "lhs is in neighoring chamber of rhs." << std::endl;
00424 
00425       // rhs local direction info
00426       /*
00427       double rhs_dydz = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).y()/
00428                         geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).z();
00429       double rhs_dxdz = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).x()/
00430                         geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).z();
00431       double rhs_dydz_err = rhs.second->localDirectionError().yy();
00432       double rhs_dxdz_err = rhs.second->localDirectionError().xx();
00433       */
00434       
00435       //rhs global position info
00436       double rhs_phi = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).phi();
00437       double rhs_theta = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).theta();
00438       double rhs_z = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).z();
00439 
00440       for(std::vector<CSCSegment>::const_iterator ilhs = thesegments.begin(); ilhs != thesegments.end(); ++ilhs) {
00441         
00442         // lhs local direction info
00443         /*
00444         double lhs_dydz = geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).y()/
00445                           geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).z();
00446         double lhs_dxdz = geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).x()/
00447                           geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).z();
00448         double lhs_dydz_err = ilhs->localDirectionError().yy();
00449         double lhs_dxdz_err = ilhs->localDirectionError().xx();
00450         */
00451         
00452         //lhs global position info
00453         double lhs_phi = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).phi();
00454         double lhs_theta = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).theta();
00455         double lhs_z = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).z();
00456         /*
00457           std::cout << "RHS Segment Parameters:" << std::endl
00458           << "\t RHS Phi   : " << rhs_phi << std::endl
00459           << "\t RHS Theta : " << rhs_theta << std::endl
00460           << "\t RHS dx/dz : " << rhs_dxdz << " +- " << rhs_dxdz_err << std::endl
00461           << "\t RHS dy/dz : " << rhs_dydz << " +- " << rhs_dydz_err << std::endl; 
00462           
00463           std::cout << "LHS Segment Parameters:" << std::endl
00464           << "\t LHS Phi   : " << lhs_phi << std::endl
00465           << "\t LHS Theta : " << lhs_theta << std::endl
00466           << "\t LHS dx/dz : " << lhs_dxdz << " +- " << lhs_dxdz_err << std::endl
00467           << "\t LHS dy/dz : " << lhs_dydz << " +- " << lhs_dydz_err << std::endl; 
00468         */
00469         
00470         double phidiff = (fabs(rhs_phi - lhs_phi) > 2*M_PI ? fabs(rhs_phi - lhs_phi) - 2*M_PI : fabs(rhs_phi - lhs_phi));
00471 
00472         if(phidiff < OverlapDPhi && fabs(rhs_theta - lhs_theta) < OverlapDTheta && 
00473            fabs(rhs_z) < fabs(lhs_z) && rhs_z*lhs_z > 0) // phi overlap region is 3.5 degrees and rhs is infront of lhs
00474           result = true;
00475       } // loop over duplicate segments
00476     }// neighboring chamber
00477   } // same endcap,station,ring
00478 
00479   return result;
00480 }
00481 
00482 bool MuonMesh::isClusteredWith(const std::pair<CSCDetId,CSCSegmentRef>& lhs, 
00483                                const std::pair<CSCDetId,CSCSegmentRef>& rhs) const
00484 {
00485   bool result(false);
00486   
00487   // try to implement the simple case first just back-to-back segments without treatment of ME1/a ganging
00488   // ME1a should be a simple extension of this
00489 
00490   //std::cout << lhs.first << ' ' << rhs.first << std::endl;
00491 
00492   if(rhs.first.endcap() == lhs.first.endcap() && lhs.first.station() < rhs.first.station()) {
00493           
00494       std::vector<CSCSegment> thesegments;
00495       thesegments.push_back(*(lhs.second));
00496       /*
00497       if(lhs.second->isME11a_duplicate())
00498         thesegments.insert(thesegments.begin(),
00499                            lhs.second->duplicateSegments().begin(),
00500                            lhs.second->duplicateSegments().end());
00501       */
00502       //std::cout << "lhs is in neighoring chamber of rhs." << std::endl;
00503 
00504       // rhs local direction info
00505       /*
00506       double rhs_dydz = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).y()/
00507                         geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).z();
00508       double rhs_dxdz = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).x()/
00509                         geometry_->chamber(rhs.first)->toGlobal(rhs.second->localDirection()).z();
00510       double rhs_dydz_err = rhs.second->localDirectionError().yy();
00511       double rhs_dxdz_err = rhs.second->localDirectionError().xx();
00512       */
00513       
00514       //rhs global position info
00515       double rhs_phi = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).phi();
00516       double rhs_theta = geometry_->chamber(rhs.first)->toGlobal(rhs.second->localPosition()).theta();
00517 
00518       for(std::vector<CSCSegment>::const_iterator ilhs = thesegments.begin(); ilhs != thesegments.end(); ++ilhs) {
00519         
00520         // lhs local direction info
00521         /*
00522         double lhs_dydz = geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).y()/
00523                           geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).z();
00524         double lhs_dxdz = geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).x()/
00525                           geometry_->chamber(lhs.first)->toGlobal(ilhs->localDirection()).z();
00526         double lhs_dydz_err = ilhs->localDirectionError().yy();
00527         double lhs_dxdz_err = ilhs->localDirectionError().xx();
00528         */
00529 
00530         //lhs global position info
00531         double lhs_phi = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).phi();
00532         double lhs_theta = geometry_->chamber(lhs.first)->toGlobal(ilhs->localPosition()).theta();
00533         /*
00534           std::cout << "RHS Segment Parameters:" << std::endl
00535           << "\t RHS Phi   : " << rhs_phi << std::endl
00536           << "\t RHS Theta : " << rhs_theta << std::endl
00537           << "\t RHS dx/dz : " << rhs_dxdz << " +- " << rhs_dxdz_err << std::endl
00538           << "\t RHS dy/dz : " << rhs_dydz << " +- " << rhs_dydz_err << std::endl; 
00539           
00540           std::cout << "LHS Segment Parameters:" << std::endl
00541           << "\t LHS Phi   : " << lhs_phi << std::endl
00542           << "\t LHS Theta : " << lhs_theta << std::endl
00543           << "\t LHS dx/dz : " << lhs_dxdz << " +- " << lhs_dxdz_err << std::endl
00544           << "\t LHS dy/dz : " << lhs_dydz << " +- " << lhs_dydz_err << std::endl; 
00545         */
00546 
00547         double phidiff = (fabs(rhs_phi - lhs_phi) > 2*M_PI ? fabs(rhs_phi - lhs_phi) - 2*M_PI : fabs(rhs_phi - lhs_phi));
00548         
00549         if(phidiff < ClusterDPhi && fabs(rhs_theta - lhs_theta) < ClusterDTheta) // phi overlap region is 37 degrees
00550           result = true;
00551       } // loop over duplicate segments    
00552   } // same endcap,station,ring
00553 
00554   return result;
00555 }
CMSSW_4_4_3_patch1/src/RecoMuon/MuonIdentification/src/MuonMesh.cc
