CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoMuon/MuonIdentification/src/MuonShowerInformationFiller.cc

Go to the documentation of this file.
00001 
00014 #include "RecoMuon/MuonIdentification/interface/MuonShowerInformationFiller.h"
00015 
00016 // system include files
00017 #include <memory>
00018 #include <algorithm>
00019 #include <iostream>
00020 
00021 // user include files
00022 #include "FWCore/Framework/interface/Event.h"
00023 #include "FWCore/PluginManager/interface/PluginManager.h"
00024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00025 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00026 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
00027 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
00028 #include "DataFormats/DetId/interface/DetId.h"
00029 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00030 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
00031 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
00032 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
00033 #include "DataFormats/GeometrySurface/interface/Bounds.h"
00034 #include "RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h"
00035 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00036 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
00037 #include "MagneticField/Engine/interface/MagneticField.h"
00038 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00039 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00040 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00041 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
00042 #include "TrackingTools/DetLayers/interface/DetLayer.h"
00043 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
00044 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00045 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
00046 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
00047 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
00048 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00049 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBreaker.h"
00050 #include "RecoMuon/TrackingTools/interface/MuonTrajectoryBuilder.h"
00051 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
00052 #include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h"
00053 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00054 
00055 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
00056 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00057 
00058 #include "DataFormats/MuonReco/interface/MuonShower.h"
00059 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00060 #include "DataFormats/MuonReco/interface/Muon.h"
00061 #include "DataFormats/TrackReco/interface/Track.h"
00062 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00063 
00064 using namespace std;
00065 using namespace edm;
00066 
00067 //
00068 // Constructor
00069 //
00070 MuonShowerInformationFiller::MuonShowerInformationFiller(const edm::ParameterSet& par) :
00071   theService(0),
00072   theDTRecHitLabel(par.getParameter<InputTag>("DTRecSegmentLabel")),
00073   theCSCRecHitLabel(par.getParameter<InputTag>("CSCRecSegmentLabel")),
00074   theCSCSegmentsLabel(par.getParameter<InputTag>("CSCSegmentLabel")),
00075   theDT4DRecSegmentLabel(par.getParameter<InputTag>("DT4DRecSegmentLabel"))
00076 {
00077 
00078   edm::ParameterSet serviceParameters = par.getParameter<edm::ParameterSet>("ServiceParameters");
00079   theService = new MuonServiceProxy(serviceParameters);
00080 
00081   theTrackerRecHitBuilderName = par.getParameter<string>("TrackerRecHitBuilder");
00082   theMuonRecHitBuilderName = par.getParameter<string>("MuonRecHitBuilder");
00083 
00084   theCacheId_TRH = 0;
00085   theCacheId_MT = 0;
00086 
00087   category_ = "MuonShowerInformationFiller";
00088 
00089   for (int istat = 0; istat < 4; istat++) {
00090       theStationShowerDeltaR.push_back(0.);
00091       theStationShowerTSize.push_back(0.);
00092       theAllStationHits.push_back(0);
00093       theCorrelatedStationHits.push_back(0);
00094   }
00095 
00096 }
00097 
00098 //
00099 // Destructor
00100 //
00101 MuonShowerInformationFiller::~MuonShowerInformationFiller() {
00102   if (theService) delete theService;
00103 }
00104 
00105 //
00106 //Fill the MuonShower struct
00107 //
00108 reco::MuonShower MuonShowerInformationFiller::fillShowerInformation( const reco::Muon& muon, const edm::Event& iEvent, const edm::EventSetup& iSetup) {
00109 
00110   reco::MuonShower returnShower;
00111   
00112   // Update the services
00113   theService->update(iSetup);
00114   setEvent(iEvent);
00115   setServices(theService->eventSetup());
00116 
00117   fillHitsByStation(muon);
00118   std::vector<int> nStationHits = theAllStationHits;
00119   std::vector<int> nStationCorrelatedHits = theCorrelatedStationHits;
00120   std::vector<float> stationShowerSizeT = theStationShowerTSize;
00121   std::vector<float> stationShowerDeltaR = theStationShowerDeltaR;
00122   
00123   returnShower.nStationHits = nStationHits; 
00124   returnShower.nStationCorrelatedHits = nStationCorrelatedHits;
00125   returnShower.stationShowerSizeT = stationShowerSizeT;
00126   returnShower.stationShowerDeltaR = stationShowerDeltaR; 
00127 
00128   return returnShower;
00129 
00130 }
00131 
00132 //
00133 // Set Event
00134 //
00135 void MuonShowerInformationFiller::setEvent(const edm::Event& event) {
00136 
00137   // get all the necesary products
00138   event.getByLabel(theDTRecHitLabel, theDTRecHits);
00139   event.getByLabel(theCSCRecHitLabel, theCSCRecHits);
00140   event.getByLabel(theCSCSegmentsLabel, theCSCSegments);
00141   event.getByLabel(theDT4DRecSegmentLabel, theDT4DRecSegments);
00142 
00143   for (int istat = 0; istat < 4; istat++) {
00144       theStationShowerDeltaR.at(istat) = 0.;
00145       theStationShowerTSize.at(istat) = 0.;
00146       theAllStationHits.at(istat) = 0;
00147       theCorrelatedStationHits.at(istat) = 0;
00148   }
00149 }
00150 
00151 
00152 //
00153 // Set services
00154 //
00155 void MuonShowerInformationFiller::setServices(const EventSetup& setup) {
00156 
00157   // DetLayer Geometry
00158   setup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry);
00159   setup.get<IdealMagneticFieldRecord>().get(theField);
00160   setup.get<TrackerRecoGeometryRecord>().get(theTracker);
00161   setup.get<MuonGeometryRecord>().get(theCSCGeometry);
00162   setup.get<MuonGeometryRecord>().get(theDTGeometry);
00163 
00164   // Transient Rechit Builders
00165   unsigned long long newCacheId_TRH = setup.get<TransientRecHitRecord>().cacheIdentifier();
00166   if ( newCacheId_TRH != theCacheId_TRH ) {
00167     setup.get<TransientRecHitRecord>().get(theTrackerRecHitBuilderName,theTrackerRecHitBuilder);
00168     setup.get<TransientRecHitRecord>().get(theMuonRecHitBuilderName,theMuonRecHitBuilder);
00169   }
00170 
00171 }
00172 
00173 //
00174 // Get hits owned by segments
00175 //
00176 TransientTrackingRecHit::ConstRecHitContainer
00177 MuonShowerInformationFiller::hitsFromSegments(const GeomDet* geomDet, 
00178                             edm::Handle<DTRecSegment4DCollection> dtSegments, 
00179                             edm::Handle<CSCSegmentCollection> cscSegments) const {
00180 
00181   MuonTransientTrackingRecHit::MuonRecHitContainer segments;
00182 
00183   DetId geoId = geomDet->geographicalId();
00184 
00185   if (geoId.subdetId() == MuonSubdetId::DT) {
00186 
00187     DTChamberId chamberId(geoId.rawId());
00188 
00189     // loop on segments 4D                 
00190     DTRecSegment4DCollection::id_iterator chamberIdIt;
00191     for (chamberIdIt = dtSegments->id_begin();
00192          chamberIdIt != dtSegments->id_end();
00193          ++chamberIdIt){
00194 
00195       if (*chamberIdIt != chamberId) continue;
00196 
00197       // Get the range for the corresponding ChamberId              
00198       DTRecSegment4DCollection::range  range = dtSegments->get((*chamberIdIt));
00199   
00200       for (DTRecSegment4DCollection::const_iterator iseg = range.first;
00201         iseg!=range.second;++iseg) {
00202         if (iseg->dimension() != 4) continue;
00203         segments.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet,&*iseg));
00204       }
00205     } 
00206   }
00207   else if (geoId.subdetId() == MuonSubdetId::CSC) {
00208 
00209       CSCDetId did(geoId.rawId());
00210 
00211       for (  CSCSegmentCollection::id_iterator chamberId = cscSegments->id_begin();
00212          chamberId != cscSegments->id_end(); ++chamberId) {
00213 
00214       if ((*chamberId).chamber() != did.chamber()) continue;
00215 
00216       // Get the range for the corresponding ChamberId                                                                                                                     
00217       CSCSegmentCollection::range  range = cscSegments->get((*chamberId));
00218 
00219       for (CSCSegmentCollection::const_iterator iseg = range.first;
00220         iseg!=range.second;++iseg) {
00221         if (iseg->dimension() != 3) continue;
00222         segments.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet,&*iseg));
00223        }
00224     }
00225   }
00226   else {
00227     LogTrace(category_) << "Segments are not built in RPCs" << endl;
00228   }
00229 
00230   TransientTrackingRecHit::ConstRecHitContainer allhitscorrelated;
00231 
00232   if (segments.empty()) return allhitscorrelated;  
00233 
00234   TransientTrackingRecHit::ConstRecHitPointer muonRecHit(segments.front().get());
00235   allhitscorrelated = MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit,2);
00236 
00237   if (segments.size() == 1) return allhitscorrelated;
00238 
00239   for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator iseg = segments.begin() + 1;
00240        iseg != segments.end(); ++iseg) {
00241 
00242     TransientTrackingRecHit::ConstRecHitPointer muonRecHit((*iseg).get());
00243     TransientTrackingRecHit::ConstRecHitContainer hits1 = MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit,2);
00244 
00245     for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit1 = hits1.begin();
00246          ihit1 != hits1.end(); ++ihit1 ) {
00247 
00248       bool usedbefore = false;       
00249       //unused      DetId thisID = (*ihit1)->geographicalId();
00250       //LocalPoint lp1dinsegHit = (*ihit1)->localPosition();
00251       GlobalPoint gp1dinsegHit = (*ihit1)->globalPosition();
00252 
00253       for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit2 = allhitscorrelated.begin();
00254            ihit2 != allhitscorrelated.end(); ++ihit2 ) {
00255 
00256         //unused        DetId thisID2 = (*ihit2)->geographicalId();
00257         //LocalPoint lp1dinsegHit2 = (*ihit2)->localPosition();
00258         GlobalPoint gp1dinsegHit2 = (*ihit2)->globalPosition();
00259 
00260         if ( (gp1dinsegHit2 - gp1dinsegHit).mag() < 1.0 ) usedbefore = true;
00261 
00262       }
00263       if ( !usedbefore ) allhitscorrelated.push_back(*ihit1);
00264     }
00265   }
00266 
00267   return allhitscorrelated;
00268 
00269 }
00270 
00271 
00272 //
00273 // Find cluster
00274 //
00275 MuonTransientTrackingRecHit::MuonRecHitContainer 
00276 MuonShowerInformationFiller::findPhiCluster(MuonTransientTrackingRecHit::MuonRecHitContainer& muonRecHits, 
00277                               const GlobalPoint& refpoint) const {
00278 
00279   if ( muonRecHits.empty() ) return muonRecHits;
00280 
00281   //clustering step by phi
00282   float step = 0.05;
00283   MuonTransientTrackingRecHit::MuonRecHitContainer result;
00284 
00285   stable_sort(muonRecHits.begin(), muonRecHits.end(), AbsLessDPhi(refpoint));
00286 
00287   for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ihit = muonRecHits.begin(); ihit != muonRecHits.end() - 1; ++ihit) {
00288       if (fabs(deltaPhi((*(ihit+1))->globalPosition().phi(), (*ihit)->globalPosition().phi() )) < step) {
00289           result.push_back(*ihit);  
00290         } else {
00291            break;
00292        }
00293   } 
00294 
00295   LogTrace(category_) <<  "phi front: " << muonRecHits.front()->globalPosition().phi() << endl;    
00296   LogTrace(category_) <<  "phi back: " << muonRecHits.back()->globalPosition().phi() << endl;
00297 
00298   return result;
00299 
00300 }
00301 
00302 //
00303 //
00304 //
00305 TransientTrackingRecHit::ConstRecHitContainer
00306 MuonShowerInformationFiller::findThetaCluster(TransientTrackingRecHit::ConstRecHitContainer& muonRecHits,
00307                               const GlobalPoint& refpoint) const {
00308 
00309   if ( muonRecHits.empty() ) return muonRecHits;
00310 
00311   //clustering step by theta
00312   float step = 0.05;
00313   TransientTrackingRecHit::ConstRecHitContainer result;
00314 
00315   stable_sort(muonRecHits.begin(), muonRecHits.end(), AbsLessDTheta(refpoint));
00316 
00317   for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit = muonRecHits.begin(); ihit != muonRecHits.end() - 1; ++ihit) {
00318       if (fabs((*(ihit+1))->globalPosition().theta() - (*ihit)->globalPosition().theta() ) < step) {
00319           result.push_back(*ihit);
00320         } else {
00321            break;
00322        }
00323   }
00324 
00325   return result;
00326 
00327 }
00328 
00329 //
00330 //Used to treat overlap region
00331 //
00332 MuonTransientTrackingRecHit::MuonRecHitContainer
00333 MuonShowerInformationFiller::findPerpCluster(MuonTransientTrackingRecHit::MuonRecHitContainer& muonRecHits) const {
00334 
00335   if ( muonRecHits.empty() ) return muonRecHits;
00336 
00337   stable_sort(muonRecHits.begin(), muonRecHits.end(), LessPerp());
00338 
00339   MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator
00340   seedhit = min_element(muonRecHits.begin(), muonRecHits.end(), LessPerp());
00341 
00342   MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ihigh = seedhit;
00343   MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ilow = seedhit;
00344 
00345   float step = 0.1;
00346   while (ihigh != muonRecHits.end()-1 && ( fabs((*(ihigh+1))->globalPosition().perp() - (*ihigh)->globalPosition().perp() ) < step)  ) {
00347     ihigh++;
00348   }
00349   while (ilow != muonRecHits.begin() && ( fabs((*ilow)->globalPosition().perp() - (*(ilow -1))->globalPosition().perp()) < step ) ) {
00350     ilow--;
00351   }
00352 
00353   MuonTransientTrackingRecHit::MuonRecHitContainer result(ilow, ihigh);
00354 
00355   return result;
00356 
00357 }
00358 
00359 //
00360 // Get compatible dets
00361 //
00362 vector<const GeomDet*> MuonShowerInformationFiller::getCompatibleDets(const reco::Track& track) const {
00363 
00364   vector<const GeomDet*> total;
00365   total.reserve(1000);
00366 
00367   LogTrace(category_)  << "Consider a track " << track.p() << " eta: " << track.eta() << " phi " << track.phi() << endl;
00368 
00369   
00370   TrajectoryStateOnSurface innerTsos = trajectoryStateTransform::innerStateOnSurface(track, *theService->trackingGeometry(), &*theService->magneticField());
00371   TrajectoryStateOnSurface outerTsos = trajectoryStateTransform::outerStateOnSurface(track, *theService->trackingGeometry(), &*theService->magneticField());
00372 
00373   GlobalPoint innerPos = innerTsos.globalPosition();
00374   GlobalPoint outerPos = outerTsos.globalPosition();
00375 
00376   vector<GlobalPoint> allCrossingPoints;
00377 
00378   const vector<DetLayer*>& dtlayers = theService->detLayerGeometry()->allDTLayers();
00379 
00380   for (vector<DetLayer*>::const_iterator iLayer = dtlayers.begin(); iLayer != dtlayers.end(); ++iLayer) {
00381 
00382     // crossing points of track with cylinder
00383     GlobalPoint xPoint = crossingPoint(innerPos, outerPos, dynamic_cast<const BarrelDetLayer*>(*iLayer));
00384       
00385     // check if point is inside the detector
00386     if ((fabs(xPoint.y()) < 1000.0) && (fabs(xPoint.z()) < 1500 ) && 
00387              (!(xPoint.y() == 0 && xPoint.x() == 0 && xPoint.z() == 0))) allCrossingPoints.push_back(xPoint);
00388   }
00389 
00390   stable_sort(allCrossingPoints.begin(), allCrossingPoints.end(), LessMag(innerPos) );
00391 
00392   vector<const GeomDet*> tempDT;
00393 
00394   for (vector<GlobalPoint>::const_iterator ipos = allCrossingPoints.begin(); ipos != allCrossingPoints.end(); ++ipos) {
00395 
00396     tempDT = dtPositionToDets(*ipos);
00397     vector<const GeomDet*>::const_iterator begin = tempDT.begin();
00398     vector<const GeomDet*>::const_iterator end = tempDT.end();
00399      
00400     for (; begin!=end;++begin) {
00401       total.push_back(*begin);
00402     } 
00403 
00404   }
00405   allCrossingPoints.clear();
00406 
00407   const vector<DetLayer*>& csclayers = theService->detLayerGeometry()->allCSCLayers();
00408   for (vector<DetLayer*>::const_iterator iLayer = csclayers.begin(); iLayer != csclayers.end(); ++iLayer) {
00409 
00410     GlobalPoint xPoint = crossingPoint(innerPos, outerPos, dynamic_cast<const ForwardDetLayer*>(*iLayer));
00411 
00412     // check if point is inside the detector
00413     if ((fabs(xPoint.y()) < 1000.0) && (fabs(xPoint.z()) < 1500.0) 
00414            && (!(xPoint.y() == 0 && xPoint.x() == 0 && xPoint.z() == 0))) allCrossingPoints.push_back(xPoint);
00415    }
00416    stable_sort(allCrossingPoints.begin(), allCrossingPoints.end(), LessMag(innerPos) );
00417 
00418    vector<const GeomDet*> tempCSC;
00419    for (vector<GlobalPoint>::const_iterator ipos = allCrossingPoints.begin(); ipos != allCrossingPoints.end(); ++ipos) {
00420 
00421      tempCSC = cscPositionToDets(*ipos);
00422      vector<const GeomDet*>::const_iterator begin = tempCSC.begin();
00423      vector<const GeomDet*>::const_iterator end = tempCSC.end();
00424 
00425      for (; begin!=end;++begin) {
00426        total.push_back(*begin);
00427      }
00428 
00429   }
00430 
00431   return total;
00432 
00433 }
00434 
00435 
00436 //
00437 // Intersection point of track with barrel layer
00438 //
00439 GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1, 
00440                                             const GlobalPoint& p2,
00441                                             const BarrelDetLayer* dl) const {
00442 
00443   const BoundCylinder& bc = dl->specificSurface();
00444   return crossingPoint(p1, p2, bc);
00445 
00446 }
00447 
00448 GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1,
00449                                             const GlobalPoint& p2, 
00450                                             const Cylinder& cyl) const {
00451 
00452   float radius = cyl.radius();
00453 
00454   GlobalVector dp = p1 - p2;
00455   float slope = dp.x()/dp.y();
00456   float a = p1.x() - slope * p1.y();
00457 
00458   float n2 = (1 + slope * slope);
00459   float n1 = 2*a*slope;
00460   float n0 = a*a - radius*radius;
00461 
00462   float y1 = 9999;
00463   float y2 = 9999;
00464   if ( n1*n1 - 4*n2*n0 > 0 ) {
00465     y1 = (-n1 + sqrt(n1*n1 - 4*n2*n0) ) / (2 * n2);
00466     y2 = (-n1 - sqrt(n1*n1 - 4*n2*n0) ) / (2 * n2);
00467   }
00468 
00469   float x1 = p1.x() + slope * (y1 - p1.y());
00470   float x2 = p1.x() + slope * (y2 - p1.y());
00471 
00472   float slopeZ = dp.z()/dp.y();
00473 
00474   float z1 = p1.z() + slopeZ * (y1 - p1.y());
00475   float z2 = p1.z() + slopeZ * (y2 - p1.y());
00476 
00477   // there are two crossing points, return the one that is in the same quadrant as point of extrapolation
00478   if ((p2.x()*x1 > 0) && (y1*p2.y() > 0) && (z1*p2.z() > 0)) {
00479     return GlobalPoint(x1, y1, z1); 
00480   } else { 
00481     return GlobalPoint(x2, y2, z2);
00482   }
00483 
00484 }
00485 
00486 
00487 //
00488 // Intersection point of track with a forward layer
00489 //
00490 GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1, 
00491                                             const GlobalPoint& p2, 
00492                                             const ForwardDetLayer* dl) const {
00493 
00494   const BoundDisk& bc = dl->specificSurface();
00495   return crossingPoint(p1, p2, bc);
00496    
00497 }  
00498 
00499 GlobalPoint MuonShowerInformationFiller::crossingPoint(const GlobalPoint& p1, 
00500                                             const GlobalPoint& p2, 
00501                                             const BoundDisk& disk) const {
00502 
00503   float diskZ = disk.position().z();
00504   int endcap =  diskZ > 0 ? 1 : (diskZ < 0 ? -1 : 0);
00505   diskZ = diskZ + endcap*dynamic_cast<const SimpleDiskBounds&>(disk.bounds()).thickness()/2.;
00506 
00507   GlobalVector dp = p1 - p2;
00508 
00509   float slopeZ = dp.z()/dp.y();
00510   float y1 = diskZ / slopeZ;
00511 
00512   float slopeX = dp.z()/dp.x();
00513   float x1 = diskZ / slopeX;
00514 
00515   float z1 = diskZ;
00516 
00517   if (p2.z()*z1 > 0) {
00518     return GlobalPoint(x1, y1, z1);
00519   } else {
00520     return GlobalPoint(0, 0, 0);
00521   }
00522 
00523 }
00524 
00525 
00526 //
00527 // GeomDets along the track in DT
00528 //
00529 vector<const GeomDet*> MuonShowerInformationFiller::dtPositionToDets(const GlobalPoint& gp) const {
00530 
00531    int minwheel = -3;
00532    int maxwheel = -3;
00533    if ( gp.z() < -680.0 ) { minwheel = -3; maxwheel = -3;}
00534    else if ( gp.z() < -396.0 ) { minwheel = -2; maxwheel = -1;}
00535    else if (gp.z() < -126.8) { minwheel = -2; maxwheel = 0; }
00536    else if (gp.z() < 126.8) { minwheel = -1; maxwheel = 1; }
00537    else if (gp.z() < 396.0) { minwheel = 0; maxwheel = 2; }
00538    else if (gp.z() < 680.0) { minwheel = 1; maxwheel = 2; }
00539    else { minwheel = 3; maxwheel = 3; }
00540 
00541    int station = 5;
00542    if ( gp.perp() > 680.0 && gp.perp() < 755.0 ) station = 4;
00543    else if ( gp.perp() > 580.0 ) station = 3;
00544    else if ( gp.perp() > 480.0 ) station = 2;
00545    else if ( gp.perp() > 380.0 ) station = 1;
00546    else station = 0;
00547 
00548    vector<int> sectors;
00549 
00550    float phistep = M_PI/6;
00551 
00552    float phigp = (float)gp.phi();
00553 
00554    if ( fabs(deltaPhi(phigp, 0*phistep)) < phistep ) sectors.push_back(1);
00555    if ( fabs(deltaPhi(phigp, phistep))   < phistep ) sectors.push_back(2);
00556    if ( fabs(deltaPhi(phigp, 2*phistep)) < phistep ) sectors.push_back(3);
00557    if ( fabs(deltaPhi(phigp, 3*phistep)) < phistep ) {
00558        sectors.push_back(4);
00559        if (station == 4) sectors.push_back(13);
00560    }
00561    if ( fabs(deltaPhi(phigp, 4*phistep)) < phistep ) sectors.push_back(5);
00562    if ( fabs(deltaPhi(phigp, 5*phistep)) < phistep ) sectors.push_back(6);
00563    if ( fabs(deltaPhi(phigp, 6*phistep)) < phistep ) sectors.push_back(7);
00564    if ( fabs(deltaPhi(phigp, 7*phistep)) < phistep ) sectors.push_back(8);
00565    if ( fabs(deltaPhi(phigp, 8*phistep)) < phistep ) sectors.push_back(9);
00566    if ( fabs(deltaPhi(phigp, 9*phistep)) < phistep ) {
00567        sectors.push_back(10);
00568        if (station == 4) sectors.push_back(14);
00569    }
00570    if ( fabs(deltaPhi(phigp, 10*phistep)) < phistep ) sectors.push_back(11);
00571    if ( fabs(deltaPhi(phigp, 11*phistep)) < phistep ) sectors.push_back(12);
00572 
00573    LogTrace(category_) << "DT position to dets" << endl;
00574    LogTrace(category_) << "number of sectors to consider: " << sectors.size() << endl;   
00575    LogTrace(category_) << "station: " << station << " wheels: " << minwheel << " " << maxwheel << endl;
00576 
00577    vector<const GeomDet*> result;
00578    if (station > 4 || station < 1) return result;
00579    if (minwheel > 2 || maxwheel < -2) return result;
00580 
00581    for (vector<int>::const_iterator isector = sectors.begin(); isector != sectors.end(); ++isector ) {
00582      for (int iwheel = minwheel; iwheel != maxwheel + 1; ++iwheel) {
00583        DTChamberId chamberid(iwheel, station, (*isector));
00584        result.push_back(theService->trackingGeometry()->idToDet(chamberid));
00585      }
00586    }
00587 
00588    LogTrace(category_) << "number of GeomDets for this track: " << result.size() << endl;
00589 
00590    return result;
00591 
00592 }
00593 
00594 
00595 //
00596 // GeomDets along the track in CSC
00597 //
00598 vector<const GeomDet*> MuonShowerInformationFiller::cscPositionToDets(const GlobalPoint& gp) const {
00599 
00600   // determine the endcap side
00601   int endcap = 0;
00602   if (gp.z() > 0) {endcap = 1;} else {endcap = 2;} 
00603 
00604   // determine the csc station and range of rings
00605   int station = 5;
00606 
00607   // check all rings in a station
00608   if ( fabs(gp.z()) > 1000. && fabs(gp.z()) < 1055.0 ) {
00609     station = 4;
00610   }        
00611   else if ( fabs(gp.z()) > 910.0 && fabs(gp.z()) < 965.0) {
00612     station = 3;
00613   }
00614   else if ( fabs(gp.z()) > 800.0 && fabs(gp.z()) < 860.0) {
00615     station = 2;
00616   }
00617   else if ( fabs(gp.z()) > 570.0 && fabs(gp.z()) < 730.0) {
00618     station = 1;
00619   }
00620 
00621   vector<int> sectors;
00622 
00623   float phistep1 = M_PI/18.; //for all the rings except first rings for stations > 1
00624   float phistep2 = M_PI/9.;
00625   float phigp = (float)gp.phi();
00626 
00627   int ring = -1;
00628 
00629   // determine the ring
00630   if (station == 1) {
00631 
00632 //FIX ME!!! 
00633       if (gp.perp() >  100 && gp.perp() < 270) ring = 1;
00634       else if (gp.perp() > 270 && gp.perp() < 450) ring = 2;
00635       else if (gp.perp() > 450 && gp.perp() < 695) ring = 3;
00636       else if (gp.perp() > 100 && gp.perp() < 270) ring = 4;
00637 
00638   } 
00639   else if (station == 2) {
00640       
00641       if (gp.perp() > 140 && gp.perp() < 350) ring = 1;
00642       else if (gp.perp() > 350 && gp.perp() < 700) ring = 2;
00643 
00644   }
00645   else if (station == 3) {
00646 
00647       if (gp.perp() > 160 && gp.perp() < 350) ring = 1;
00648       else if (gp.perp() > 350 && gp.perp() < 700) ring = 2;
00649 
00650   }
00651   else if (station == 4) {
00652 
00653       if (gp.perp() > 175 && gp.perp() < 350) ring = 1;
00654       else if (gp.perp() > 350 && gp.perp() < 700) ring = 2;
00655 
00656   }
00657 
00658   if (station > 1 && ring == 1) { 
00659 
00660    // we have 18 sectors in that case
00661    for (int i = 0; i < 18; i++) { 
00662       if ( fabs(deltaPhi(phigp, i*phistep2)) < phistep2 ) sectors.push_back(i+1);
00663     }
00664 
00665   } else {
00666 
00667    // we have 36 sectors in that case
00668    for (int i = 0; i < 36; i++) {
00669      if ( fabs(deltaPhi(phigp, i*phistep1)) < phistep1 ) sectors.push_back(i+1);
00670    }
00671 }
00672 
00673    LogTrace(category_) << "CSC position to dets" << endl;
00674    LogTrace(category_) << "ring: " << ring << endl;
00675    LogTrace(category_) << "endcap: " << endcap << endl;
00676    LogTrace(category_) << "station: " << station << endl;
00677    LogTrace(category_) << "CSC number of sectors to consider: " << sectors.size() << endl;
00678 
00679 
00680  // check exceptional cases
00681    vector<const GeomDet*> result;
00682    if (station > 4 || station < 1) return result;
00683    if (endcap == 0) return result;
00684    if (ring == -1) return result;
00685 
00686    int minlayer = 1;
00687    int maxlayer = 6;
00688 
00689    for (vector<int>::const_iterator isector = sectors.begin(); isector != sectors.end(); ++isector) {
00690      for (int ilayer = minlayer; ilayer != maxlayer + 1; ++ ilayer) {
00691        CSCDetId cscid(endcap, station, ring, (*isector), ilayer);
00692        result.push_back(theService->trackingGeometry()->idToDet(cscid));
00693      }
00694    }
00695 
00696   return result;
00697 
00698 }
00699 
00700 //
00701 //Set class members
00702 //
00703 void MuonShowerInformationFiller::fillHitsByStation(const reco::Muon& muon) {
00704 
00705   reco::TrackRef track;
00706   if ( muon.isGlobalMuon() )            track = muon.globalTrack();
00707   else if ( muon.isStandAloneMuon() )    track = muon.outerTrack();
00708   else return;
00709 
00710   // split 1D rechits by station
00711   vector<MuonRecHitContainer> muonRecHits(4);
00712 
00713   // split rechits from segs by station
00714   vector<TransientTrackingRecHit::ConstRecHitContainer> muonCorrelatedHits(4);  
00715 
00716   // get vector of GeomDets compatible with a track
00717   vector<const GeomDet*> compatibleLayers = getCompatibleDets(*track);
00718 
00719   // for special cases: CSC station 1
00720   MuonRecHitContainer tmpCSC1;
00721   bool dtOverlapToCheck = false;
00722   bool cscOverlapToCheck = false;
00723 
00724   for (vector<const GeomDet*>::const_iterator igd = compatibleLayers.begin(); igd != compatibleLayers.end(); igd++ )  {
00725 
00726     // get det id
00727     DetId geoId = (*igd)->geographicalId();
00728 
00729     // skip tracker hits
00730     if (geoId.det()!= DetId::Muon) continue;
00731 
00732     // DT 
00733     if ( geoId.subdetId() == MuonSubdetId::DT ) {
00734 
00735       // get station
00736       DTChamberId detid(geoId.rawId());
00737       int station = detid.station();
00738       int wheel = detid.wheel();
00739 
00740       // get rechits from segments per station
00741       TransientTrackingRecHit::ConstRecHitContainer muonCorrelatedHitsTmp = hitsFromSegments(*igd, theDT4DRecSegments, theCSCSegments);     
00742       TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_begin = muonCorrelatedHitsTmp.begin();
00743       TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_end = muonCorrelatedHitsTmp.end();
00744 
00745       for (; hits_begin!= hits_end;++hits_begin) {
00746         muonCorrelatedHits.at(station-1).push_back(*hits_begin);
00747       }
00748 
00749       //check overlap certain wheels and stations
00750       if (abs(wheel) == 2 && station != 4 &&  station != 1) dtOverlapToCheck = true;
00751 
00752       // loop over all superlayers of a DT chamber
00753       for (int isuperlayer = DTChamberId::minSuperLayerId; isuperlayer != DTChamberId::maxSuperLayerId + 1; ++isuperlayer) {
00754         // loop over all layers inside the superlayer
00755         for (int ilayer = DTChamberId::minLayerId; ilayer != DTChamberId::maxLayerId+1; ++ilayer) {
00756           DTLayerId lid(detid, isuperlayer, ilayer);
00757           DTRecHitCollection::range dRecHits = theDTRecHits->get(lid);
00758           for (DTRecHitCollection::const_iterator rechit = dRecHits.first; rechit != dRecHits.second;++rechit) {
00759             vector<const TrackingRecHit*> subrechits = (*rechit).recHits();
00760             for (vector<const TrackingRecHit*>::iterator irechit = subrechits.begin(); irechit != subrechits.end(); ++irechit) {
00761                      muonRecHits.at(station-1).push_back(MuonTransientTrackingRecHit::specificBuild((&**igd),&**irechit));
00762                 }
00763              }
00764           }
00765        }
00766     }
00767     else if (geoId.subdetId() == MuonSubdetId::CSC) {
00768 
00769       // get station
00770       CSCDetId did(geoId.rawId());
00771       int station = did.station();
00772       int ring = did.ring();
00773 
00774       //get rechits from segments by station      
00775       TransientTrackingRecHit::ConstRecHitContainer muonCorrelatedHitsTmp = hitsFromSegments(*igd, theDT4DRecSegments, theCSCSegments);
00776       TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_begin = muonCorrelatedHitsTmp.begin();
00777       TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_end = muonCorrelatedHitsTmp.end();
00778 
00779       for (; hits_begin!= hits_end;++hits_begin) {
00780         muonCorrelatedHits.at(station-1).push_back(*hits_begin);
00781       }
00782 
00783       if ((station == 1 && ring == 3) && dtOverlapToCheck) cscOverlapToCheck = true;
00784 
00785       // split 1D rechits by station
00786       CSCRecHit2DCollection::range dRecHits = theCSCRecHits->get(did);
00787       for (CSCRecHit2DCollection::const_iterator rechit = dRecHits.first; rechit != dRecHits.second; ++rechit) {
00788 
00789         if (!cscOverlapToCheck) {
00790            muonRecHits.at(station-1).push_back(MuonTransientTrackingRecHit::specificBuild((&**igd),&*rechit));
00791            } else {
00792              tmpCSC1.push_back(MuonTransientTrackingRecHit::specificBuild((&**igd),&*rechit));
00793 
00794              //sort by perp, then insert to appropriate container
00795              MuonRecHitContainer temp = findPerpCluster(tmpCSC1);
00796              if (temp.empty()) continue;
00797 
00798              float center;
00799              if (temp.size() > 1) {
00800                center = (temp.front()->globalPosition().perp() + temp.back()->globalPosition().perp())/2.;
00801              } else {
00802                center = temp.front()->globalPosition().perp();
00803              }
00804              temp.clear();
00805              
00806              if (center > 550.) {
00807                 muonRecHits.at(2).insert(muonRecHits.at(2).end(),tmpCSC1.begin(),tmpCSC1.end());
00808               } else {
00809               muonRecHits.at(1).insert(muonRecHits.at(1).end(),tmpCSC1.begin(),tmpCSC1.end());
00810            }
00811            tmpCSC1.clear();
00812          }
00813        }
00814      } else if (geoId.subdetId() == MuonSubdetId::RPC) {
00815        LogTrace(category_) << "Wrong subdet id" << endl;
00816      }
00817   }//loop over GeomDets compatible with a track
00818 
00819   // calculate number of all and correlated hits    
00820   for (int stat = 0; stat < 4; stat++) {
00821     theCorrelatedStationHits[stat] = muonCorrelatedHits.at(stat).size();     
00822     theAllStationHits[stat] = muonRecHits[stat].size();
00823   }
00824   LogTrace(category_) << "Hits used by the segments, by station "
00825        << theCorrelatedStationHits.at(0) << " "
00826        << theCorrelatedStationHits.at(1) << " "
00827        << theCorrelatedStationHits.at(2) << " "
00828        << theCorrelatedStationHits.at(3) << endl;
00829 
00830   LogTrace(category_) << "All DT 1D/CSC 2D  hits, by station "
00831        << theAllStationHits.at(0) << " "
00832        << theAllStationHits.at(1) << " "
00833        << theAllStationHits.at(2) << " "
00834        << theAllStationHits.at(3) << endl;
00835 
00836   //station shower sizes
00837   MuonTransientTrackingRecHit::MuonRecHitContainer muonRecHitsPhiTemp, muonRecHitsPhiBest;
00838   TransientTrackingRecHit::ConstRecHitContainer muonRecHitsThetaTemp, muonRecHitsThetaBest;
00839 
00840   // send station hits to the clustering algorithm
00841   for ( int stat = 0; stat != 4; stat++ ) {
00842     if (!muonRecHits[stat].empty()) {
00843       stable_sort(muonRecHits[stat].begin(), muonRecHits[stat].end(), LessPhi());
00844 
00845       float dphimax = 0;
00846       for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator iseed = muonRecHits[stat].begin(); iseed != muonRecHits[stat].end(); ++iseed) {
00847           if (!(*iseed)->isValid()) continue;
00848           GlobalPoint refpoint = (*iseed)->globalPosition(); //starting from the one with smallest value of phi
00849           muonRecHitsPhiTemp.clear();
00850           muonRecHitsPhiTemp = findPhiCluster(muonRecHits[stat], refpoint); //get clustered hits for this iseed
00851       if (muonRecHitsPhiTemp.size() > 1) {
00852          float dphi = fabs(deltaPhi((float)muonRecHitsPhiTemp.back()->globalPosition().phi(), (float)muonRecHitsPhiTemp.front()->globalPosition().phi()));
00853          if (dphi > dphimax) {
00854             dphimax = dphi;
00855             muonRecHitsPhiBest = muonRecHitsPhiTemp;
00856           }
00857         } //at least two hits
00858       }//loop over seeds 
00859 
00860      //fill showerTs
00861       if (!muonRecHitsPhiBest.empty()) {
00862         muonRecHits[stat] = muonRecHitsPhiBest;
00863         stable_sort(muonRecHits[stat].begin(), muonRecHits[stat].end(), LessAbsMag());
00864         muonRecHits[stat].front();
00865         GlobalPoint refpoint = muonRecHits[stat].front()->globalPosition();
00866         theStationShowerTSize.at(stat) = refpoint.mag() * dphimax;
00867       }
00868 
00869      //for theta
00870      if (!muonCorrelatedHits.at(stat).empty()) {
00871 
00872        float dthetamax = 0;
00873        for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator iseed = muonCorrelatedHits.at(stat).begin(); iseed != muonCorrelatedHits.at(stat).end(); ++iseed) {
00874            if (!(*iseed)->isValid()) continue;
00875            GlobalPoint refpoint = (*iseed)->globalPosition(); //starting from the one with smallest value of phi
00876            muonRecHitsThetaTemp.clear();
00877            muonRecHitsThetaTemp = findThetaCluster(muonCorrelatedHits.at(stat), refpoint);
00878        }//loop over seeds 
00879        if (muonRecHitsThetaTemp.size() > 1) {
00880          float dtheta = fabs((float)muonRecHitsThetaTemp.back()->globalPosition().theta() - (float)muonRecHitsThetaTemp.front()->globalPosition().theta());
00881          if (dtheta > dthetamax) {
00882            dthetamax = dtheta;
00883            muonRecHitsThetaBest = muonRecHitsThetaTemp;
00884          }
00885        } //at least two hits
00886      }//not empty container2
00887 
00888      //fill deltaRs
00889      if (muonRecHitsThetaBest.size() > 1 && muonRecHitsPhiBest.size() > 1)
00890        theStationShowerDeltaR.at(stat) = sqrt(pow(muonRecHitsPhiBest.front()->globalPosition().phi()-muonRecHitsPhiBest.back()->globalPosition().phi(),2)+pow(muonRecHitsThetaBest.front()->globalPosition().theta()-muonRecHitsThetaBest.back()->globalPosition().theta(),2));
00891 
00892         }//not empty container
00893       }//loop over station
00894 
00895        LogTrace(category_) << "deltaR around a track containing all the station hits, by station "
00896        << theStationShowerDeltaR.at(0) << " "
00897        << theStationShowerDeltaR.at(1) << " "
00898        << theStationShowerDeltaR.at(2) << " "
00899        << theStationShowerDeltaR.at(3) << endl;
00900 
00901 
00902        LogTrace(category_) << "Transverse cluster size, by station "
00903        << theStationShowerTSize.at(0) << " "
00904        << theStationShowerTSize.at(1) << " "
00905        << theStationShowerTSize.at(2) << " "
00906        << theStationShowerTSize.at(3) << endl;
00907 
00908   return;
00909 
00910 }