CMS 3D CMS Logo

Classes | Public Types | Public Member Functions | Protected Member Functions | Private Member Functions | Private Attributes

MuonShowerInformationFiller Class Reference

#include <MuonShowerInformationFiller.h>

List of all members.

Classes

struct  AbsLessDPhi
struct  AbsLessDTheta
struct  LessAbsMag
struct  LessDPhi
struct  LessMag
struct  LessPerp
struct  LessPhi

Public Types

typedef
MuonTransientTrackingRecHit::ConstMuonRecHitPointer 
ConstMuonRecHitPointer
typedef
TransientTrackingRecHit::ConstRecHitContainer 
ConstRecHitContainer
typedef
MuonTransientTrackingRecHit::MuonRecHitContainer 
MuonRecHitContainer

Public Member Functions

void fillHitsByStation (const reco::Muon &)
reco::MuonShower fillShowerInformation (const reco::Muon &muon, const edm::Event &, const edm::EventSetup &)
 fill muon shower variables
 MuonShowerInformationFiller ()
 constructors
 MuonShowerInformationFiller (const edm::ParameterSet &)
virtual void setEvent (const edm::Event &)
 pass the Event to the algorithm at each event
void setServices (const edm::EventSetup &)
 set the services needed
 ~MuonShowerInformationFiller ()
 destructor

Protected Member Functions

const MuonServiceProxygetService () const

Private Member Functions

GlobalPoint crossingPoint (const GlobalPoint &, const GlobalPoint &, const BarrelDetLayer *) const
GlobalPoint crossingPoint (const GlobalPoint &, const GlobalPoint &, const Cylinder &) const
GlobalPoint crossingPoint (const GlobalPoint &, const GlobalPoint &, const BoundDisk &) const
GlobalPoint crossingPoint (const GlobalPoint &, const GlobalPoint &, const ForwardDetLayer *) const
std::vector< const GeomDet * > cscPositionToDets (const GlobalPoint &) const
std::vector< const GeomDet * > dtPositionToDets (const GlobalPoint &) const
MuonRecHitContainer findPerpCluster (MuonRecHitContainer &muonRecHits) const
MuonRecHitContainer findPhiCluster (MuonRecHitContainer &, const GlobalPoint &) const
TransientTrackingRecHit::ConstRecHitContainer findThetaCluster (TransientTrackingRecHit::ConstRecHitContainer &, const GlobalPoint &) const
std::vector< const GeomDet * > getCompatibleDets (const reco::Track &) const
TransientTrackingRecHit::ConstRecHitContainer hitsFromSegments (const GeomDet *, edm::Handle< DTRecSegment4DCollection >, edm::Handle< CSCSegmentCollection >) const

Private Attributes

std::string category_
std::vector< int > theAllStationHits
unsigned long long theCacheId_MT
unsigned long long theCacheId_TRH
std::vector< int > theCorrelatedStationHits
edm::ESHandle< CSCGeometrytheCSCGeometry
edm::InputTag theCSCRecHitLabel
edm::Handle
< CSCRecHit2DCollection
theCSCRecHits
edm::Handle< CSCSegmentCollectiontheCSCSegments
edm::InputTag theCSCSegmentsLabel
edm::InputTag theDT4DRecSegmentLabel
edm::Handle
< DTRecSegment4DCollection
theDT4DRecSegments
edm::ESHandle< DTGeometrytheDTGeometry
edm::InputTag theDTRecHitLabel
edm::Handle< DTRecHitCollectiontheDTRecHits
edm::ESHandle< MagneticFieldtheField
edm::ESHandle
< TransientTrackingRecHitBuilder
theMuonRecHitBuilder
std::string theMuonRecHitBuilderName
MuonServiceProxytheService
std::vector< float > theStationShowerDeltaR
std::vector< float > theStationShowerTSize
edm::ESHandle
< GeometricSearchTracker
theTracker
edm::ESHandle
< TransientTrackingRecHitBuilder
theTrackerRecHitBuilder
std::string theTrackerRecHitBuilderName
edm::ESHandle
< GlobalTrackingGeometry
theTrackingGeometry

Detailed Description

Description: class for muon shower identification

Date:
2011/03/21 22:41:25
Revision:
1.3
Author:
: A. Svyatkovskiy, Purdue University

Description: class for muon shower identification

Date:
2011/12/23 05:08:50
Revision:
1.7
Author:
: A. Svyatkovskiy, Purdue University

Definition at line 57 of file MuonShowerInformationFiller.h.


Member Typedef Documentation

Definition at line 63 of file MuonShowerInformationFiller.h.

Definition at line 61 of file MuonShowerInformationFiller.h.

Definition at line 62 of file MuonShowerInformationFiller.h.


Constructor & Destructor Documentation

MuonShowerInformationFiller::MuonShowerInformationFiller ( ) [inline]

constructors

Definition at line 68 of file MuonShowerInformationFiller.h.

{};
MuonShowerInformationFiller::MuonShowerInformationFiller ( const edm::ParameterSet par)

Definition at line 70 of file MuonShowerInformationFiller.cc.

References category_, edm::ParameterSet::getParameter(), MuonServiceProxy_cff::MuonServiceProxy, theAllStationHits, theCacheId_MT, theCacheId_TRH, theCorrelatedStationHits, theMuonRecHitBuilderName, theService, theStationShowerDeltaR, theStationShowerTSize, and theTrackerRecHitBuilderName.

                                                                                   :
  theService(0),
  theDTRecHitLabel(par.getParameter<InputTag>("DTRecSegmentLabel")),
  theCSCRecHitLabel(par.getParameter<InputTag>("CSCRecSegmentLabel")),
  theCSCSegmentsLabel(par.getParameter<InputTag>("CSCSegmentLabel")),
  theDT4DRecSegmentLabel(par.getParameter<InputTag>("DT4DRecSegmentLabel"))
{

  edm::ParameterSet serviceParameters = par.getParameter<edm::ParameterSet>("ServiceParameters");
  theService = new MuonServiceProxy(serviceParameters);

  theTrackerRecHitBuilderName = par.getParameter<string>("TrackerRecHitBuilder");
  theMuonRecHitBuilderName = par.getParameter<string>("MuonRecHitBuilder");

  theCacheId_TRH = 0;
  theCacheId_MT = 0;

  category_ = "MuonShowerInformationFiller";

  for (int istat = 0; istat < 4; istat++) {
      theStationShowerDeltaR.push_back(0.);
      theStationShowerTSize.push_back(0.);
      theAllStationHits.push_back(0);
      theCorrelatedStationHits.push_back(0);
  }

}
MuonShowerInformationFiller::~MuonShowerInformationFiller ( )

destructor

Definition at line 101 of file MuonShowerInformationFiller.cc.

References theService.

                                                          {
  if (theService) delete theService;
}

Member Function Documentation

GlobalPoint MuonShowerInformationFiller::crossingPoint ( const GlobalPoint p1,
const GlobalPoint p2,
const BarrelDetLayer dl 
) const [private]

Definition at line 439 of file MuonShowerInformationFiller.cc.

References BarrelDetLayer::specificSurface().

Referenced by crossingPoint(), and getCompatibleDets().

                                                                            {

  const BoundCylinder& bc = dl->specificSurface();
  return crossingPoint(p1, p2, bc);

}
GlobalPoint MuonShowerInformationFiller::crossingPoint ( const GlobalPoint p1,
const GlobalPoint p2,
const ForwardDetLayer dl 
) const [private]

Definition at line 490 of file MuonShowerInformationFiller.cc.

References crossingPoint(), and ForwardDetLayer::specificSurface().

                                                                             {

  const BoundDisk& bc = dl->specificSurface();
  return crossingPoint(p1, p2, bc);
   
}  
GlobalPoint MuonShowerInformationFiller::crossingPoint ( const GlobalPoint p1,
const GlobalPoint p2,
const BoundDisk disk 
) const [private]

Definition at line 499 of file MuonShowerInformationFiller.cc.

References BoundSurface::bounds(), Reference_intrackfit_cff::endcap, p2, GloballyPositioned< T >::position(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

                                                                         {

  float diskZ = disk.position().z();
  int endcap =  diskZ > 0 ? 1 : (diskZ < 0 ? -1 : 0);
  diskZ = diskZ + endcap*dynamic_cast<const SimpleDiskBounds&>(disk.bounds()).thickness()/2.;

  GlobalVector dp = p1 - p2;

  float slopeZ = dp.z()/dp.y();
  float y1 = diskZ / slopeZ;

  float slopeX = dp.z()/dp.x();
  float x1 = diskZ / slopeX;

  float z1 = diskZ;

  if (p2.z()*z1 > 0) {
    return GlobalPoint(x1, y1, z1);
  } else {
    return GlobalPoint(0, 0, 0);
  }

}
GlobalPoint MuonShowerInformationFiller::crossingPoint ( const GlobalPoint p1,
const GlobalPoint p2,
const Cylinder cyl 
) const [private]

Definition at line 448 of file MuonShowerInformationFiller.cc.

References a, n0, p2, CosmicsPD_Skims::radius, Cylinder::radius(), slope, mathSSE::sqrt(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

                                                                       {

  float radius = cyl.radius();

  GlobalVector dp = p1 - p2;
  float slope = dp.x()/dp.y();
  float a = p1.x() - slope * p1.y();

  float n2 = (1 + slope * slope);
  float n1 = 2*a*slope;
  float n0 = a*a - radius*radius;

  float y1 = 9999;
  float y2 = 9999;
  if ( n1*n1 - 4*n2*n0 > 0 ) {
    y1 = (-n1 + sqrt(n1*n1 - 4*n2*n0) ) / (2 * n2);
    y2 = (-n1 - sqrt(n1*n1 - 4*n2*n0) ) / (2 * n2);
  }

  float x1 = p1.x() + slope * (y1 - p1.y());
  float x2 = p1.x() + slope * (y2 - p1.y());

  float slopeZ = dp.z()/dp.y();

  float z1 = p1.z() + slopeZ * (y1 - p1.y());
  float z2 = p1.z() + slopeZ * (y2 - p1.y());

  // there are two crossing points, return the one that is in the same quadrant as point of extrapolation
  if ((p2.x()*x1 > 0) && (y1*p2.y() > 0) && (z1*p2.z() > 0)) {
    return GlobalPoint(x1, y1, z1); 
  } else { 
    return GlobalPoint(x2, y2, z2);
  }

}
vector< const GeomDet * > MuonShowerInformationFiller::cscPositionToDets ( const GlobalPoint gp) const [private]

Definition at line 598 of file MuonShowerInformationFiller.cc.

References category_, SiPixelRawToDigiRegional_cfi::deltaPhi, Reference_intrackfit_cff::endcap, i, LogTrace, M_PI, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), query::result, relativeConstraints::ring, mergeVDriftHistosByStation::sectors, relativeConstraints::station, theService, and PV3DBase< T, PVType, FrameType >::z().

Referenced by getCompatibleDets().

                                                                                                 {

  // determine the endcap side
  int endcap = 0;
  if (gp.z() > 0) {endcap = 1;} else {endcap = 2;} 

  // determine the csc station and range of rings
  int station = 5;

  // check all rings in a station
  if ( fabs(gp.z()) > 1000. && fabs(gp.z()) < 1055.0 ) {
    station = 4;
  }        
  else if ( fabs(gp.z()) > 910.0 && fabs(gp.z()) < 965.0) {
    station = 3;
  }
  else if ( fabs(gp.z()) > 800.0 && fabs(gp.z()) < 860.0) {
    station = 2;
  }
  else if ( fabs(gp.z()) > 570.0 && fabs(gp.z()) < 730.0) {
    station = 1;
  }

  vector<int> sectors;

  float phistep1 = M_PI/18.; //for all the rings except first rings for stations > 1
  float phistep2 = M_PI/9.;
  float phigp = (float)gp.phi();

  int ring = -1;

  // determine the ring
  if (station == 1) {

//FIX ME!!! 
      if (gp.perp() >  100 && gp.perp() < 270) ring = 1;
      else if (gp.perp() > 270 && gp.perp() < 450) ring = 2;
      else if (gp.perp() > 450 && gp.perp() < 695) ring = 3;
      else if (gp.perp() > 100 && gp.perp() < 270) ring = 4;

  } 
  else if (station == 2) {
      
      if (gp.perp() > 140 && gp.perp() < 350) ring = 1;
      else if (gp.perp() > 350 && gp.perp() < 700) ring = 2;

  }
  else if (station == 3) {

      if (gp.perp() > 160 && gp.perp() < 350) ring = 1;
      else if (gp.perp() > 350 && gp.perp() < 700) ring = 2;

  }
  else if (station == 4) {

      if (gp.perp() > 175 && gp.perp() < 350) ring = 1;
      else if (gp.perp() > 350 && gp.perp() < 700) ring = 2;

  }

  if (station > 1 && ring == 1) { 

   // we have 18 sectors in that case
   for (int i = 0; i < 18; i++) { 
      if ( fabs(deltaPhi(phigp, i*phistep2)) < phistep2 ) sectors.push_back(i+1);
    }

  } else {

   // we have 36 sectors in that case
   for (int i = 0; i < 36; i++) {
     if ( fabs(deltaPhi(phigp, i*phistep1)) < phistep1 ) sectors.push_back(i+1);
   }
}

   LogTrace(category_) << "CSC position to dets" << endl;
   LogTrace(category_) << "ring: " << ring << endl;
   LogTrace(category_) << "endcap: " << endcap << endl;
   LogTrace(category_) << "station: " << station << endl;
   LogTrace(category_) << "CSC number of sectors to consider: " << sectors.size() << endl;


 // check exceptional cases
   vector<const GeomDet*> result;
   if (station > 4 || station < 1) return result;
   if (endcap == 0) return result;
   if (ring == -1) return result;

   int minlayer = 1;
   int maxlayer = 6;

   for (vector<int>::const_iterator isector = sectors.begin(); isector != sectors.end(); ++isector) {
     for (int ilayer = minlayer; ilayer != maxlayer + 1; ++ ilayer) {
       CSCDetId cscid(endcap, station, ring, (*isector), ilayer);
       result.push_back(theService->trackingGeometry()->idToDet(cscid));
     }
   }

  return result;

}
vector< const GeomDet * > MuonShowerInformationFiller::dtPositionToDets ( const GlobalPoint gp) const [private]

Definition at line 529 of file MuonShowerInformationFiller.cc.

References category_, SiPixelRawToDigiRegional_cfi::deltaPhi, LogTrace, M_PI, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), query::result, mergeVDriftHistosByStation::sectors, relativeConstraints::station, theService, and PV3DBase< T, PVType, FrameType >::z().

Referenced by getCompatibleDets().

                                                                                                {

   int minwheel = -3;
   int maxwheel = -3;
   if ( gp.z() < -680.0 ) { minwheel = -3; maxwheel = -3;}
   else if ( gp.z() < -396.0 ) { minwheel = -2; maxwheel = -1;}
   else if (gp.z() < -126.8) { minwheel = -2; maxwheel = 0; }
   else if (gp.z() < 126.8) { minwheel = -1; maxwheel = 1; }
   else if (gp.z() < 396.0) { minwheel = 0; maxwheel = 2; }
   else if (gp.z() < 680.0) { minwheel = 1; maxwheel = 2; }
   else { minwheel = 3; maxwheel = 3; }

   int station = 5;
   if ( gp.perp() > 680.0 && gp.perp() < 755.0 ) station = 4;
   else if ( gp.perp() > 580.0 ) station = 3;
   else if ( gp.perp() > 480.0 ) station = 2;
   else if ( gp.perp() > 380.0 ) station = 1;
   else station = 0;

   vector<int> sectors;

   float phistep = M_PI/6;

   float phigp = (float)gp.phi();

   if ( fabs(deltaPhi(phigp, 0*phistep)) < phistep ) sectors.push_back(1);
   if ( fabs(deltaPhi(phigp, phistep))   < phistep ) sectors.push_back(2);
   if ( fabs(deltaPhi(phigp, 2*phistep)) < phistep ) sectors.push_back(3);
   if ( fabs(deltaPhi(phigp, 3*phistep)) < phistep ) {
       sectors.push_back(4);
       if (station == 4) sectors.push_back(13);
   }
   if ( fabs(deltaPhi(phigp, 4*phistep)) < phistep ) sectors.push_back(5);
   if ( fabs(deltaPhi(phigp, 5*phistep)) < phistep ) sectors.push_back(6);
   if ( fabs(deltaPhi(phigp, 6*phistep)) < phistep ) sectors.push_back(7);
   if ( fabs(deltaPhi(phigp, 7*phistep)) < phistep ) sectors.push_back(8);
   if ( fabs(deltaPhi(phigp, 8*phistep)) < phistep ) sectors.push_back(9);
   if ( fabs(deltaPhi(phigp, 9*phistep)) < phistep ) {
       sectors.push_back(10);
       if (station == 4) sectors.push_back(14);
   }
   if ( fabs(deltaPhi(phigp, 10*phistep)) < phistep ) sectors.push_back(11);
   if ( fabs(deltaPhi(phigp, 11*phistep)) < phistep ) sectors.push_back(12);

   LogTrace(category_) << "DT position to dets" << endl;
   LogTrace(category_) << "number of sectors to consider: " << sectors.size() << endl;   
   LogTrace(category_) << "station: " << station << " wheels: " << minwheel << " " << maxwheel << endl;

   vector<const GeomDet*> result;
   if (station > 4 || station < 1) return result;
   if (minwheel > 2 || maxwheel < -2) return result;

   for (vector<int>::const_iterator isector = sectors.begin(); isector != sectors.end(); ++isector ) {
     for (int iwheel = minwheel; iwheel != maxwheel + 1; ++iwheel) {
       DTChamberId chamberid(iwheel, station, (*isector));
       result.push_back(theService->trackingGeometry()->idToDet(chamberid));
     }
   }

   LogTrace(category_) << "number of GeomDets for this track: " << result.size() << endl;

   return result;

}
void MuonShowerInformationFiller::fillHitsByStation ( const reco::Muon muon)

Definition at line 703 of file MuonShowerInformationFiller.cc.

References abs, begin, category_, CSC(), SiPixelRawToDigiRegional_cfi::deltaPhi, DetId::det(), cond::rpcobgas::detid, GeomDetEnumerators::DT, relativeConstraints::empty, end, findPerpCluster(), findPhiCluster(), findThetaCluster(), getCompatibleDets(), reco::Muon::globalTrack(), hitsFromSegments(), iseed, reco::Muon::isGlobalMuon(), reco::Muon::isStandAloneMuon(), LogTrace, PV3DBase< T, PVType, FrameType >::mag(), DTChamberId::maxLayerId, DTChamberId::maxSuperLayerId, DTChamberId::minLayerId, DTChamberId::minSuperLayerId, DetId::Muon, reco::Muon::outerTrack(), funct::pow(), DetId::rawId(), relativeConstraints::ring, dedefs::RPC, MuonTransientTrackingRecHit::specificBuild(), mathSSE::sqrt(), testRegression::stat, relativeConstraints::station, DetId::subdetId(), groupFilesInBlocks::temp, theAllStationHits, theCorrelatedStationHits, theCSCRecHits, theCSCSegments, theDT4DRecSegments, theDTRecHits, theStationShowerDeltaR, and theStationShowerTSize.

Referenced by fillShowerInformation().

                                                                        {

  reco::TrackRef track;
  if ( muon.isGlobalMuon() )            track = muon.globalTrack();
  else if ( muon.isStandAloneMuon() )    track = muon.outerTrack();
  else return;

  // split 1D rechits by station
  vector<MuonRecHitContainer> muonRecHits(4);

  // split rechits from segs by station
  vector<TransientTrackingRecHit::ConstRecHitContainer> muonCorrelatedHits(4);  

  // get vector of GeomDets compatible with a track
  vector<const GeomDet*> compatibleLayers = getCompatibleDets(*track);

  // for special cases: CSC station 1
  MuonRecHitContainer tmpCSC1;
  bool dtOverlapToCheck = false;
  bool cscOverlapToCheck = false;

  for (vector<const GeomDet*>::const_iterator igd = compatibleLayers.begin(); igd != compatibleLayers.end(); igd++ )  {

    // get det id
    DetId geoId = (*igd)->geographicalId();

    // skip tracker hits
    if (geoId.det()!= DetId::Muon) continue;

    // DT 
    if ( geoId.subdetId() == MuonSubdetId::DT ) {

      // get station
      DTChamberId detid(geoId.rawId());
      int station = detid.station();
      int wheel = detid.wheel();

      // get rechits from segments per station
      TransientTrackingRecHit::ConstRecHitContainer muonCorrelatedHitsTmp = hitsFromSegments(*igd, theDT4DRecSegments, theCSCSegments);     
      TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_begin = muonCorrelatedHitsTmp.begin();
      TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_end = muonCorrelatedHitsTmp.end();

      for (; hits_begin!= hits_end;++hits_begin) {
        muonCorrelatedHits.at(station-1).push_back(*hits_begin);
      }

      //check overlap certain wheels and stations
      if (abs(wheel) == 2 && station != 4 &&  station != 1) dtOverlapToCheck = true;

      // loop over all superlayers of a DT chamber
      for (int isuperlayer = DTChamberId::minSuperLayerId; isuperlayer != DTChamberId::maxSuperLayerId + 1; ++isuperlayer) {
        // loop over all layers inside the superlayer
        for (int ilayer = DTChamberId::minLayerId; ilayer != DTChamberId::maxLayerId+1; ++ilayer) {
          DTLayerId lid(detid, isuperlayer, ilayer);
          DTRecHitCollection::range dRecHits = theDTRecHits->get(lid);
          for (DTRecHitCollection::const_iterator rechit = dRecHits.first; rechit != dRecHits.second;++rechit) {
            vector<const TrackingRecHit*> subrechits = (*rechit).recHits();
            for (vector<const TrackingRecHit*>::iterator irechit = subrechits.begin(); irechit != subrechits.end(); ++irechit) {
                     muonRecHits.at(station-1).push_back(MuonTransientTrackingRecHit::specificBuild((&**igd),&**irechit));
                }
             }
          }
       }
    }
    else if (geoId.subdetId() == MuonSubdetId::CSC) {

      // get station
      CSCDetId did(geoId.rawId());
      int station = did.station();
      int ring = did.ring();

      //get rechits from segments by station      
      TransientTrackingRecHit::ConstRecHitContainer muonCorrelatedHitsTmp = hitsFromSegments(*igd, theDT4DRecSegments, theCSCSegments);
      TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_begin = muonCorrelatedHitsTmp.begin();
      TransientTrackingRecHit::ConstRecHitContainer::const_iterator hits_end = muonCorrelatedHitsTmp.end();

      for (; hits_begin!= hits_end;++hits_begin) {
        muonCorrelatedHits.at(station-1).push_back(*hits_begin);
      }

      if ((station == 1 && ring == 3) && dtOverlapToCheck) cscOverlapToCheck = true;

      // split 1D rechits by station
      CSCRecHit2DCollection::range dRecHits = theCSCRecHits->get(did);
      for (CSCRecHit2DCollection::const_iterator rechit = dRecHits.first; rechit != dRecHits.second; ++rechit) {

        if (!cscOverlapToCheck) {
           muonRecHits.at(station-1).push_back(MuonTransientTrackingRecHit::specificBuild((&**igd),&*rechit));
           } else {
             tmpCSC1.push_back(MuonTransientTrackingRecHit::specificBuild((&**igd),&*rechit));

             //sort by perp, then insert to appropriate container
             MuonRecHitContainer temp = findPerpCluster(tmpCSC1);
             if (temp.empty()) continue;

             float center;
             if (temp.size() > 1) {
               center = (temp.front()->globalPosition().perp() + temp.back()->globalPosition().perp())/2.;
             } else {
               center = temp.front()->globalPosition().perp();
             }
             temp.clear();
             
             if (center > 550.) {
                muonRecHits.at(2).insert(muonRecHits.at(2).end(),tmpCSC1.begin(),tmpCSC1.end());
              } else {
              muonRecHits.at(1).insert(muonRecHits.at(1).end(),tmpCSC1.begin(),tmpCSC1.end());
           }
           tmpCSC1.clear();
         }
       }
     } else if (geoId.subdetId() == MuonSubdetId::RPC) {
       LogTrace(category_) << "Wrong subdet id" << endl;
     }
  }//loop over GeomDets compatible with a track

  // calculate number of all and correlated hits    
  for (int stat = 0; stat < 4; stat++) {
    theCorrelatedStationHits[stat] = muonCorrelatedHits.at(stat).size();     
    theAllStationHits[stat] = muonRecHits[stat].size();
  }
  LogTrace(category_) << "Hits used by the segments, by station "
       << theCorrelatedStationHits.at(0) << " "
       << theCorrelatedStationHits.at(1) << " "
       << theCorrelatedStationHits.at(2) << " "
       << theCorrelatedStationHits.at(3) << endl;

  LogTrace(category_) << "All DT 1D/CSC 2D  hits, by station "
       << theAllStationHits.at(0) << " "
       << theAllStationHits.at(1) << " "
       << theAllStationHits.at(2) << " "
       << theAllStationHits.at(3) << endl;

  //station shower sizes
  MuonTransientTrackingRecHit::MuonRecHitContainer muonRecHitsPhiTemp, muonRecHitsPhiBest;
  TransientTrackingRecHit::ConstRecHitContainer muonRecHitsThetaTemp, muonRecHitsThetaBest;

  // send station hits to the clustering algorithm
  for ( int stat = 0; stat != 4; stat++ ) {
    if (!muonRecHits[stat].empty()) {
      stable_sort(muonRecHits[stat].begin(), muonRecHits[stat].end(), LessPhi());

      float dphimax = 0;
      for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator iseed = muonRecHits[stat].begin(); iseed != muonRecHits[stat].end(); ++iseed) {
          if (!(*iseed)->isValid()) continue;
          GlobalPoint refpoint = (*iseed)->globalPosition(); //starting from the one with smallest value of phi
          muonRecHitsPhiTemp.clear();
          muonRecHitsPhiTemp = findPhiCluster(muonRecHits[stat], refpoint); //get clustered hits for this iseed
      if (muonRecHitsPhiTemp.size() > 1) {
         float dphi = fabs(deltaPhi((float)muonRecHitsPhiTemp.back()->globalPosition().phi(), (float)muonRecHitsPhiTemp.front()->globalPosition().phi()));
         if (dphi > dphimax) {
            dphimax = dphi;
            muonRecHitsPhiBest = muonRecHitsPhiTemp;
          }
        } //at least two hits
      }//loop over seeds 

     //fill showerTs
      if (!muonRecHitsPhiBest.empty()) {
        muonRecHits[stat] = muonRecHitsPhiBest;
        stable_sort(muonRecHits[stat].begin(), muonRecHits[stat].end(), LessAbsMag());
        muonRecHits[stat].front();
        GlobalPoint refpoint = muonRecHits[stat].front()->globalPosition();
        theStationShowerTSize.at(stat) = refpoint.mag() * dphimax;
      }

     //for theta
     if (!muonCorrelatedHits.at(stat).empty()) {

       float dthetamax = 0;
       for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator iseed = muonCorrelatedHits.at(stat).begin(); iseed != muonCorrelatedHits.at(stat).end(); ++iseed) {
           if (!(*iseed)->isValid()) continue;
           GlobalPoint refpoint = (*iseed)->globalPosition(); //starting from the one with smallest value of phi
           muonRecHitsThetaTemp.clear();
           muonRecHitsThetaTemp = findThetaCluster(muonCorrelatedHits.at(stat), refpoint);
       }//loop over seeds 
       if (muonRecHitsThetaTemp.size() > 1) {
         float dtheta = fabs((float)muonRecHitsThetaTemp.back()->globalPosition().theta() - (float)muonRecHitsThetaTemp.front()->globalPosition().theta());
         if (dtheta > dthetamax) {
           dthetamax = dtheta;
           muonRecHitsThetaBest = muonRecHitsThetaTemp;
         }
       } //at least two hits
     }//not empty container2

     //fill deltaRs
     if (muonRecHitsThetaBest.size() > 1 && muonRecHitsPhiBest.size() > 1)
       theStationShowerDeltaR.at(stat) = sqrt(pow(muonRecHitsPhiBest.front()->globalPosition().phi()-muonRecHitsPhiBest.back()->globalPosition().phi(),2)+pow(muonRecHitsThetaBest.front()->globalPosition().theta()-muonRecHitsThetaBest.back()->globalPosition().theta(),2));

        }//not empty container
      }//loop over station

       LogTrace(category_) << "deltaR around a track containing all the station hits, by station "
       << theStationShowerDeltaR.at(0) << " "
       << theStationShowerDeltaR.at(1) << " "
       << theStationShowerDeltaR.at(2) << " "
       << theStationShowerDeltaR.at(3) << endl;


       LogTrace(category_) << "Transverse cluster size, by station "
       << theStationShowerTSize.at(0) << " "
       << theStationShowerTSize.at(1) << " "
       << theStationShowerTSize.at(2) << " "
       << theStationShowerTSize.at(3) << endl;

  return;

}
reco::MuonShower MuonShowerInformationFiller::fillShowerInformation ( const reco::Muon muon,
const edm::Event iEvent,
const edm::EventSetup iSetup 
)

fill muon shower variables

Definition at line 108 of file MuonShowerInformationFiller.cc.

References fillHitsByStation(), reco::MuonShower::nStationCorrelatedHits, reco::MuonShower::nStationHits, setEvent(), setServices(), reco::MuonShower::stationShowerDeltaR, reco::MuonShower::stationShowerSizeT, theAllStationHits, theCorrelatedStationHits, theService, theStationShowerDeltaR, and theStationShowerTSize.

Referenced by MuonShowerInformationProducer::produce().

                                                                                                                                            {

  reco::MuonShower returnShower;
  
  // Update the services
  theService->update(iSetup);
  setEvent(iEvent);
  setServices(theService->eventSetup());

  fillHitsByStation(muon);
  std::vector<int> nStationHits = theAllStationHits;
  std::vector<int> nStationCorrelatedHits = theCorrelatedStationHits;
  std::vector<float> stationShowerSizeT = theStationShowerTSize;
  std::vector<float> stationShowerDeltaR = theStationShowerDeltaR;
  
  returnShower.nStationHits = nStationHits; 
  returnShower.nStationCorrelatedHits = nStationCorrelatedHits;
  returnShower.stationShowerSizeT = stationShowerSizeT;
  returnShower.stationShowerDeltaR = stationShowerDeltaR; 

  return returnShower;

}
MuonTransientTrackingRecHit::MuonRecHitContainer MuonShowerInformationFiller::findPerpCluster ( MuonTransientTrackingRecHit::MuonRecHitContainer muonRecHits) const [private]

Definition at line 333 of file MuonShowerInformationFiller.cc.

References perp(), query::result, and launcher::step.

Referenced by fillHitsByStation().

                                                                                                              {

  if ( muonRecHits.empty() ) return muonRecHits;

  stable_sort(muonRecHits.begin(), muonRecHits.end(), LessPerp());

  MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator
  seedhit = min_element(muonRecHits.begin(), muonRecHits.end(), LessPerp());

  MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ihigh = seedhit;
  MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ilow = seedhit;

  float step = 0.1;
  while (ihigh != muonRecHits.end()-1 && ( fabs((*(ihigh+1))->globalPosition().perp() - (*ihigh)->globalPosition().perp() ) < step)  ) {
    ihigh++;
  }
  while (ilow != muonRecHits.begin() && ( fabs((*ilow)->globalPosition().perp() - (*(ilow -1))->globalPosition().perp()) < step ) ) {
    ilow--;
  }

  MuonTransientTrackingRecHit::MuonRecHitContainer result(ilow, ihigh);

  return result;

}
MuonTransientTrackingRecHit::MuonRecHitContainer MuonShowerInformationFiller::findPhiCluster ( MuonTransientTrackingRecHit::MuonRecHitContainer muonRecHits,
const GlobalPoint refpoint 
) const [private]

Definition at line 276 of file MuonShowerInformationFiller.cc.

References category_, SiPixelRawToDigiRegional_cfi::deltaPhi, LogTrace, phi, query::result, and launcher::step.

Referenced by fillHitsByStation().

                                                                 {

  if ( muonRecHits.empty() ) return muonRecHits;

  //clustering step by phi
  float step = 0.05;
  MuonTransientTrackingRecHit::MuonRecHitContainer result;

  stable_sort(muonRecHits.begin(), muonRecHits.end(), AbsLessDPhi(refpoint));

  for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator ihit = muonRecHits.begin(); ihit != muonRecHits.end() - 1; ++ihit) {
      if (fabs(deltaPhi((*(ihit+1))->globalPosition().phi(), (*ihit)->globalPosition().phi() )) < step) {
          result.push_back(*ihit);  
        } else {
           break;
       }
  } 

  LogTrace(category_) <<  "phi front: " << muonRecHits.front()->globalPosition().phi() << endl;    
  LogTrace(category_) <<  "phi back: " << muonRecHits.back()->globalPosition().phi() << endl;

  return result;

}
TransientTrackingRecHit::ConstRecHitContainer MuonShowerInformationFiller::findThetaCluster ( TransientTrackingRecHit::ConstRecHitContainer muonRecHits,
const GlobalPoint refpoint 
) const [private]

Definition at line 306 of file MuonShowerInformationFiller.cc.

References query::result, launcher::step, and theta().

Referenced by fillHitsByStation().

                                                                 {

  if ( muonRecHits.empty() ) return muonRecHits;

  //clustering step by theta
  float step = 0.05;
  TransientTrackingRecHit::ConstRecHitContainer result;

  stable_sort(muonRecHits.begin(), muonRecHits.end(), AbsLessDTheta(refpoint));

  for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit = muonRecHits.begin(); ihit != muonRecHits.end() - 1; ++ihit) {
      if (fabs((*(ihit+1))->globalPosition().theta() - (*ihit)->globalPosition().theta() ) < step) {
          result.push_back(*ihit);
        } else {
           break;
       }
  }

  return result;

}
vector< const GeomDet * > MuonShowerInformationFiller::getCompatibleDets ( const reco::Track track) const [private]

Definition at line 362 of file MuonShowerInformationFiller.cc.

References begin, category_, crossingPoint(), cscPositionToDets(), dtPositionToDets(), end, reco::TrackBase::eta(), TrajectoryStateOnSurface::globalPosition(), trajectoryStateTransform::innerStateOnSurface(), LogTrace, trajectoryStateTransform::outerStateOnSurface(), reco::TrackBase::p(), reco::TrackBase::phi(), theService, pileupDistInMC::total, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by fillHitsByStation().

                                                                                                  {

  vector<const GeomDet*> total;
  total.reserve(1000);

  LogTrace(category_)  << "Consider a track " << track.p() << " eta: " << track.eta() << " phi " << track.phi() << endl;

  
  TrajectoryStateOnSurface innerTsos = trajectoryStateTransform::innerStateOnSurface(track, *theService->trackingGeometry(), &*theService->magneticField());
  TrajectoryStateOnSurface outerTsos = trajectoryStateTransform::outerStateOnSurface(track, *theService->trackingGeometry(), &*theService->magneticField());

  GlobalPoint innerPos = innerTsos.globalPosition();
  GlobalPoint outerPos = outerTsos.globalPosition();

  vector<GlobalPoint> allCrossingPoints;

  const vector<DetLayer*>& dtlayers = theService->detLayerGeometry()->allDTLayers();

  for (vector<DetLayer*>::const_iterator iLayer = dtlayers.begin(); iLayer != dtlayers.end(); ++iLayer) {

    // crossing points of track with cylinder
    GlobalPoint xPoint = crossingPoint(innerPos, outerPos, dynamic_cast<const BarrelDetLayer*>(*iLayer));
      
    // check if point is inside the detector
    if ((fabs(xPoint.y()) < 1000.0) && (fabs(xPoint.z()) < 1500 ) && 
             (!(xPoint.y() == 0 && xPoint.x() == 0 && xPoint.z() == 0))) allCrossingPoints.push_back(xPoint);
  }

  stable_sort(allCrossingPoints.begin(), allCrossingPoints.end(), LessMag(innerPos) );

  vector<const GeomDet*> tempDT;

  for (vector<GlobalPoint>::const_iterator ipos = allCrossingPoints.begin(); ipos != allCrossingPoints.end(); ++ipos) {

    tempDT = dtPositionToDets(*ipos);
    vector<const GeomDet*>::const_iterator begin = tempDT.begin();
    vector<const GeomDet*>::const_iterator end = tempDT.end();
     
    for (; begin!=end;++begin) {
      total.push_back(*begin);
    } 

  }
  allCrossingPoints.clear();

  const vector<DetLayer*>& csclayers = theService->detLayerGeometry()->allCSCLayers();
  for (vector<DetLayer*>::const_iterator iLayer = csclayers.begin(); iLayer != csclayers.end(); ++iLayer) {

    GlobalPoint xPoint = crossingPoint(innerPos, outerPos, dynamic_cast<const ForwardDetLayer*>(*iLayer));

    // check if point is inside the detector
    if ((fabs(xPoint.y()) < 1000.0) && (fabs(xPoint.z()) < 1500.0) 
           && (!(xPoint.y() == 0 && xPoint.x() == 0 && xPoint.z() == 0))) allCrossingPoints.push_back(xPoint);
   }
   stable_sort(allCrossingPoints.begin(), allCrossingPoints.end(), LessMag(innerPos) );

   vector<const GeomDet*> tempCSC;
   for (vector<GlobalPoint>::const_iterator ipos = allCrossingPoints.begin(); ipos != allCrossingPoints.end(); ++ipos) {

     tempCSC = cscPositionToDets(*ipos);
     vector<const GeomDet*>::const_iterator begin = tempCSC.begin();
     vector<const GeomDet*>::const_iterator end = tempCSC.end();

     for (; begin!=end;++begin) {
       total.push_back(*begin);
     }

  }

  return total;

}
const MuonServiceProxy* MuonShowerInformationFiller::getService ( ) const [inline, protected]

Definition at line 88 of file MuonShowerInformationFiller.h.

References theService.

{ return theService; }
TransientTrackingRecHit::ConstRecHitContainer MuonShowerInformationFiller::hitsFromSegments ( const GeomDet geomDet,
edm::Handle< DTRecSegment4DCollection dtSegments,
edm::Handle< CSCSegmentCollection cscSegments 
) const [private]

Definition at line 177 of file MuonShowerInformationFiller.cc.

References MuonTransientTrackingRecHitBreaker::breakInSubRecHits(), category_, CSC(), GeomDetEnumerators::DT, GeomDet::geographicalId(), LogTrace, mag(), DetId::rawId(), MuonTransientTrackingRecHit::specificBuild(), and DetId::subdetId().

Referenced by fillHitsByStation().

                                                                               {

  MuonTransientTrackingRecHit::MuonRecHitContainer segments;

  DetId geoId = geomDet->geographicalId();

  if (geoId.subdetId() == MuonSubdetId::DT) {

    DTChamberId chamberId(geoId.rawId());

    // loop on segments 4D                 
    DTRecSegment4DCollection::id_iterator chamberIdIt;
    for (chamberIdIt = dtSegments->id_begin();
         chamberIdIt != dtSegments->id_end();
         ++chamberIdIt){

      if (*chamberIdIt != chamberId) continue;

      // Get the range for the corresponding ChamberId              
      DTRecSegment4DCollection::range  range = dtSegments->get((*chamberIdIt));
  
      for (DTRecSegment4DCollection::const_iterator iseg = range.first;
        iseg!=range.second;++iseg) {
        if (iseg->dimension() != 4) continue;
        segments.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet,&*iseg));
      }
    } 
  }
  else if (geoId.subdetId() == MuonSubdetId::CSC) {

      CSCDetId did(geoId.rawId());

      for (  CSCSegmentCollection::id_iterator chamberId = cscSegments->id_begin();
         chamberId != cscSegments->id_end(); ++chamberId) {

      if ((*chamberId).chamber() != did.chamber()) continue;

      // Get the range for the corresponding ChamberId                                                                                                                     
      CSCSegmentCollection::range  range = cscSegments->get((*chamberId));

      for (CSCSegmentCollection::const_iterator iseg = range.first;
        iseg!=range.second;++iseg) {
        if (iseg->dimension() != 3) continue;
        segments.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet,&*iseg));
       }
    }
  }
  else {
    LogTrace(category_) << "Segments are not built in RPCs" << endl;
  }

  TransientTrackingRecHit::ConstRecHitContainer allhitscorrelated;

  if (segments.empty()) return allhitscorrelated;  

  TransientTrackingRecHit::ConstRecHitPointer muonRecHit(segments.front().get());
  allhitscorrelated = MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit,2);

  if (segments.size() == 1) return allhitscorrelated;

  for (MuonTransientTrackingRecHit::MuonRecHitContainer::const_iterator iseg = segments.begin() + 1;
       iseg != segments.end(); ++iseg) {

    TransientTrackingRecHit::ConstRecHitPointer muonRecHit((*iseg).get());
    TransientTrackingRecHit::ConstRecHitContainer hits1 = MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit,2);

    for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit1 = hits1.begin();
         ihit1 != hits1.end(); ++ihit1 ) {

      bool usedbefore = false;       
      //unused      DetId thisID = (*ihit1)->geographicalId();
      //LocalPoint lp1dinsegHit = (*ihit1)->localPosition();
      GlobalPoint gp1dinsegHit = (*ihit1)->globalPosition();

      for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ihit2 = allhitscorrelated.begin();
           ihit2 != allhitscorrelated.end(); ++ihit2 ) {

        //unused        DetId thisID2 = (*ihit2)->geographicalId();
        //LocalPoint lp1dinsegHit2 = (*ihit2)->localPosition();
        GlobalPoint gp1dinsegHit2 = (*ihit2)->globalPosition();

        if ( (gp1dinsegHit2 - gp1dinsegHit).mag() < 1.0 ) usedbefore = true;

      }
      if ( !usedbefore ) allhitscorrelated.push_back(*ihit1);
    }
  }

  return allhitscorrelated;

}
void MuonShowerInformationFiller::setEvent ( const edm::Event event) [virtual]

pass the Event to the algorithm at each event

Definition at line 135 of file MuonShowerInformationFiller.cc.

References theAllStationHits, theCorrelatedStationHits, theCSCRecHitLabel, theCSCRecHits, theCSCSegments, theCSCSegmentsLabel, theDT4DRecSegmentLabel, theDT4DRecSegments, theDTRecHitLabel, theDTRecHits, theStationShowerDeltaR, and theStationShowerTSize.

Referenced by fillShowerInformation().

                                                                {

  // get all the necesary products
  event.getByLabel(theDTRecHitLabel, theDTRecHits);
  event.getByLabel(theCSCRecHitLabel, theCSCRecHits);
  event.getByLabel(theCSCSegmentsLabel, theCSCSegments);
  event.getByLabel(theDT4DRecSegmentLabel, theDT4DRecSegments);

  for (int istat = 0; istat < 4; istat++) {
      theStationShowerDeltaR.at(istat) = 0.;
      theStationShowerTSize.at(istat) = 0.;
      theAllStationHits.at(istat) = 0;
      theCorrelatedStationHits.at(istat) = 0;
  }
}
void MuonShowerInformationFiller::setServices ( const edm::EventSetup setup)

Member Data Documentation

unsigned long long MuonShowerInformationFiller::theCacheId_MT [private]

Definition at line 181 of file MuonShowerInformationFiller.h.

Referenced by MuonShowerInformationFiller().

unsigned long long MuonShowerInformationFiller::theCacheId_TRH [private]

Definition at line 180 of file MuonShowerInformationFiller.h.

Referenced by MuonShowerInformationFiller(), and setServices().

Definition at line 202 of file MuonShowerInformationFiller.h.

Referenced by setServices().

Definition at line 190 of file MuonShowerInformationFiller.h.

Referenced by setEvent().

Definition at line 194 of file MuonShowerInformationFiller.h.

Referenced by fillHitsByStation(), and setEvent().

Definition at line 195 of file MuonShowerInformationFiller.h.

Referenced by fillHitsByStation(), and setEvent().

Definition at line 191 of file MuonShowerInformationFiller.h.

Referenced by setEvent().

Definition at line 192 of file MuonShowerInformationFiller.h.

Referenced by setEvent().

Definition at line 196 of file MuonShowerInformationFiller.h.

Referenced by fillHitsByStation(), and setEvent().

Definition at line 203 of file MuonShowerInformationFiller.h.

Referenced by setServices().

Definition at line 189 of file MuonShowerInformationFiller.h.

Referenced by setEvent().

Definition at line 193 of file MuonShowerInformationFiller.h.

Referenced by fillHitsByStation(), and setEvent().

Definition at line 201 of file MuonShowerInformationFiller.h.

Referenced by setServices().

Definition at line 187 of file MuonShowerInformationFiller.h.

Referenced by setServices().

Definition at line 186 of file MuonShowerInformationFiller.h.

Referenced by MuonShowerInformationFiller(), and setServices().

Definition at line 199 of file MuonShowerInformationFiller.h.

Referenced by setServices().

Definition at line 184 of file MuonShowerInformationFiller.h.

Referenced by setServices().

Definition at line 183 of file MuonShowerInformationFiller.h.

Referenced by MuonShowerInformationFiller(), and setServices().

Definition at line 200 of file MuonShowerInformationFiller.h.

Referenced by setServices().