#include <MuonShowerInformationFiller.h>

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 More...

	MuonShowerInformationFiller ()
	constructors More...

	MuonShowerInformationFiller (const edm::ParameterSet &, edm::ConsumesCollector &)

virtual void	setEvent (const edm::Event &)
	pass the Event to the algorithm at each event More...

void	setServices (const edm::EventSetup &)
	set the services needed More...

virtual	~MuonShowerInformationFiller ()
	destructor More...

Protected Member Functions
const MuonServiceProxy *	getService () 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 ForwardDetLayer *) const

GlobalPoint	crossingPoint (const GlobalPoint &, const GlobalPoint &, const Disk &) const

std::vector< const GeomDet * >	cscPositionToDets (const GlobalPoint &) const

std::vector< const GeomDet * >	dtPositionToDets (const GlobalPoint &) const

MuonRecHitContainer	findPerpCluster (MuonRecHitContainer &muonRecHits) 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< CSCGeometry >	theCSCGeometry

edm::InputTag	theCSCRecHitLabel

edm::Handle < CSCRecHit2DCollection >	theCSCRecHits

edm::EDGetTokenT < CSCRecHit2DCollection >	theCSCRecHitToken

edm::Handle< CSCSegmentCollection >	theCSCSegments

edm::InputTag	theCSCSegmentsLabel

edm::EDGetTokenT < CSCSegmentCollection >	theCSCSegmentsToken

edm::InputTag	theDT4DRecSegmentLabel

edm::Handle < DTRecSegment4DCollection >	theDT4DRecSegments

edm::EDGetTokenT < DTRecSegment4DCollection >	theDT4DRecSegmentToken

edm::ESHandle< DTGeometry >	theDTGeometry

edm::InputTag	theDTRecHitLabel

edm::Handle< DTRecHitCollection >	theDTRecHits

edm::EDGetTokenT < DTRecHitCollection >	theDTRecHitToken

edm::ESHandle< MagneticField >	theField

edm::ESHandle < TransientTrackingRecHitBuilder >	theMuonRecHitBuilder

std::string	theMuonRecHitBuilderName

MuonServiceProxy *	theService

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

Author: : A. Svyatkovskiy, Purdue University

Definition at line 57 of file MuonShowerInformationFiller.h.

Member Typedef Documentation

typedef MuonTransientTrackingRecHit::ConstMuonRecHitPointer MuonShowerInformationFiller::ConstMuonRecHitPointer

Definition at line 63 of file MuonShowerInformationFiller.h.

typedef TransientTrackingRecHit::ConstRecHitContainer MuonShowerInformationFiller::ConstRecHitContainer

Definition at line 61 of file MuonShowerInformationFiller.h.

typedef MuonTransientTrackingRecHit::MuonRecHitContainer MuonShowerInformationFiller::MuonRecHitContainer

Definition at line 62 of file MuonShowerInformationFiller.h.

Constructor & Destructor Documentation

MuonShowerInformationFiller::MuonShowerInformationFiller ( )

inline

constructors

Definition at line 68 of file MuonShowerInformationFiller.h.

68 {};

MuonShowerInformationFiller::MuonShowerInformationFiller	(	const edm::ParameterSet &	par,
		edm::ConsumesCollector &	iC
	)

Definition at line 68 of file MuonShowerInformationFiller.cc.

References category_, edm::ConsumesCollector::consumes(), edm::ParameterSet::getParameter(), MuonServiceProxy_cff::MuonServiceProxy, theAllStationHits, theCacheId_MT, theCacheId_TRH, theCorrelatedStationHits, theCSCRecHitLabel, theCSCRecHitToken, theCSCSegmentsLabel, theCSCSegmentsToken, theDT4DRecSegmentLabel, theDT4DRecSegmentToken, theDTRecHitLabel, theDTRecHitToken, 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"))
 {
 
   theDTRecHitToken =iC.consumes<DTRecHitCollection>(theDTRecHitLabel) ;
   theCSCRecHitToken =iC.consumes<CSCRecHit2DCollection>(theCSCRecHitLabel);
   theCSCSegmentsToken=iC.consumes<CSCSegmentCollection>(theCSCSegmentsLabel);
   theDT4DRecSegmentToken=iC.consumes<DTRecSegment4DCollection>(theDT4DRecSegmentLabel);
 
 
   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 ( )

virtual

destructor

Definition at line 105 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 414 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 Cylinder &	cyl
	)		const

private

Definition at line 423 of file MuonShowerInformationFiller.cc.

References a, reco::dp, n0, p2, Cylinder::radius(), CosmicsPD_Skims::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);
   }
 
 }

GlobalPoint MuonShowerInformationFiller::crossingPoint	(	const GlobalPoint &	p1,
		const GlobalPoint &	p2,
		const ForwardDetLayer *	dl
	)		const

private

Definition at line 465 of file MuonShowerInformationFiller.cc.

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

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

GlobalPoint MuonShowerInformationFiller::crossingPoint	(	const GlobalPoint &	,
		const GlobalPoint &	,
		const Disk &
	)		const

private

vector< const GeomDet * > MuonShowerInformationFiller::cscPositionToDets ( const GlobalPoint & gp ) const

private

Definition at line 573 of file MuonShowerInformationFiller.cc.

References PV3DBase< T, PVType, FrameType >::barePhi(), category_, SiPixelRawToDigiRegional_cfi::deltaPhi, i, LogTrace, M_PI, PV3DBase< T, PVType, FrameType >::perp(), mps_fire::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.barePhi();
 
   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 504 of file MuonShowerInformationFiller.cc.

References PV3DBase< T, PVType, FrameType >::barePhi(), category_, SiPixelRawToDigiRegional_cfi::deltaPhi, LogTrace, M_PI, PV3DBase< T, PVType, FrameType >::perp(), mps_fire::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.barePhi();
 
    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 678 of file MuonShowerInformationFiller.cc.

References funct::abs(), begin, category_, filterCSVwithJSON::copy, MuonSubdetId::CSC, SiPixelRawToDigiRegional_cfi::deltaPhi, DetId::det(), MuonSubdetId::DT, end, findPerpCluster(), findThetaCluster(), getCompatibleDets(), reco::Muon::globalTrack(), hitsFromSegments(), iseed, reco::Muon::isGlobalMuon(), reco::Muon::isStandAloneMuon(), LogTrace, PV3DBase< T, PVType, FrameType >::mag(), bookConverter::max, DTChamberId::maxLayerId, DTChamberId::maxSuperLayerId, min(), DTChamberId::minLayerId, DTChamberId::minSuperLayerId, DetId::Muon, reco::Muon::outerTrack(), phi, funct::pow(), DetId::rawId(), relativeConstraints::ring, MuonSubdetId::RPC, MuonTransientTrackingRecHit::specificBuild(), mathSSE::sqrt(), relativeConstraints::station, DetId::subdetId(), groupFilesInBlocks::temp, theAllStationHits, theCorrelatedStationHits, theCSCRecHits, theCSCSegments, theDT4DRecSegments, theDTRecHits, theStationShowerDeltaR, theStationShowerTSize, and pileupCalc::upper.

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 muonRecHitsPhiBest;
   TransientTrackingRecHit::ConstRecHitContainer muonRecHitsThetaTemp, muonRecHitsThetaBest;
   // send station hits to the clustering algorithm
   for ( int stat = 0; stat != 4; stat++ ) {
     const size_t nhit = muonRecHits[stat].size();
     if (nhit < 2) continue; // Require at least 2 hits
     muonRecHitsPhiBest.clear();
     muonRecHitsPhiBest.reserve(nhit);
 
     // Cluster seeds by global position phi. Best cluster is chosen to give greatest dphi
     // Sort by phi (complexity = NLogN with enough memory, or = NLog^2N for insufficient mem)
     stable_sort(muonRecHits[stat].begin(), muonRecHits[stat].end(), LessPhi());
 
     // Search for gaps (complexity = N)
     std::vector<size_t> clUppers;
     for ( size_t ihit = 0; ihit < nhit; ++ihit ) {
       const size_t jhit = (ihit+1)%nhit;
       const double phi1 = muonRecHits[stat].at(ihit)->globalPosition().barePhi();
       const double phi2 = muonRecHits[stat].at(jhit)->globalPosition().barePhi();
 
       const double dphi = std::abs(deltaPhi(phi1, phi2));
       if ( dphi >= 0.05 ) clUppers.push_back(ihit);
     }
 
     //station shower sizes
     double dphimax = 0;
     if ( clUppers.empty() ) {
       // No gaps - there is only one cluster. Take all of them
       const double refPhi = muonRecHits[stat].at(0)->globalPosition().barePhi();
       double dphilo = 0, dphihi = 0;
       for ( auto& hit : muonRecHits[stat] ) {
         muonRecHitsPhiBest.push_back(hit);
         const double phi = hit->globalPosition().barePhi();
         dphilo = std::min(dphilo, deltaPhi(refPhi, phi));
         dphihi = std::max(dphihi, deltaPhi(refPhi, phi));
       }
       dphimax = std::abs(dphihi+dphilo);
     } else {
       // Loop over gaps to find the one with maximum dphi(begin, end)
       // By construction, number of gap must be greater than 1.
       size_t bestUpper = 0, bestLower = 0;
       for ( auto icl = clUppers.begin(); icl != clUppers.end(); ++icl ) {
         // upper bound of this cluster
         const size_t upper = *icl;
         // lower bound is +1 of preceeding upper bound
         const auto prevCl = (icl == clUppers.begin()) ? clUppers.end() : icl;
         const size_t lower = (*(prevCl-1)+1)%nhit;
 
         // At least two hit for a cluster
         if ( upper == lower ) continue;
 
         const double phi1 = muonRecHits[stat].at(upper)->globalPosition().barePhi();
         const double phi2 = muonRecHits[stat].at(lower)->globalPosition().barePhi();
 
         const double dphi = std::abs(deltaPhi(phi1, phi2));
         if ( dphimax < dphi ) {
           dphimax = dphi;
           bestUpper = upper;
           bestLower = lower;
         }
       }
 
       if ( bestUpper > bestLower ) {
         muonRecHitsPhiBest.reserve(bestUpper-bestLower+1);
         std::copy(muonRecHits[stat].begin()+bestLower,
                   muonRecHits[stat].begin()+bestUpper+1,
                   std::back_inserter(muonRecHitsPhiBest));
       } else if ( bestUpper < bestLower ) {
         muonRecHitsPhiBest.reserve(bestUpper+(nhit-bestLower)+1);
         std::copy(muonRecHits[stat].begin(),
                   muonRecHits[stat].begin()+bestUpper+1,
                   std::back_inserter(muonRecHitsPhiBest));
         std::copy(muonRecHits[stat].begin()+bestLower,
                   muonRecHits[stat].end(),
                   std::back_inserter(muonRecHitsPhiBest));
       }
     }
 
      //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(deltaPhi(muonRecHitsPhiBest.front()->globalPosition().barePhi(),muonRecHitsPhiBest.back()->globalPosition().barePhi()),2)+pow(muonRecHitsThetaBest.front()->globalPosition().theta()-muonRecHitsThetaBest.back()->globalPosition().theta(),2));
 
       }//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 112 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 308 of file MuonShowerInformationFiller.cc.

References perp(), mps_fire::result, and relval_parameters_module::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;
 
 }

TransientTrackingRecHit::ConstRecHitContainer MuonShowerInformationFiller::findThetaCluster	(	TransientTrackingRecHit::ConstRecHitContainer &	muonRecHits,
		const GlobalPoint &	refpoint
	)		const

private

Definition at line 281 of file MuonShowerInformationFiller.cc.

References mps_fire::result, relval_parameters_module::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 337 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<const DetLayer*>& dtlayers = theService->detLayerGeometry()->allDTLayers();
 
   for (auto 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<const DetLayer*>& csclayers = theService->detLayerGeometry()->allCSCLayers();
   for (auto 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

inlineprotected

Definition at line 88 of file MuonShowerInformationFiller.h.

References theService.

88 { return theService; }

MuonShowerInformationFiller::theService

MuonServiceProxy * theService

Definition: MuonShowerInformationFiller.h:97

TransientTrackingRecHit::ConstRecHitContainer MuonShowerInformationFiller::hitsFromSegments	(	const GeomDet *	geomDet,
		edm::Handle< DTRecSegment4DCollection >	dtSegments,
		edm::Handle< CSCSegmentCollection >	cscSegments
	)		const

private

Definition at line 181 of file MuonShowerInformationFiller.cc.

References MuonTransientTrackingRecHitBreaker::breakInSubRecHits(), category_, MuonSubdetId::CSC, MuonSubdetId::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());
   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));
     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 139 of file MuonShowerInformationFiller.cc.

References theAllStationHits, theCorrelatedStationHits, theCSCRecHits, theCSCRecHitToken, theCSCSegments, theCSCSegmentsToken, theDT4DRecSegments, theDT4DRecSegmentToken, theDTRecHits, theDTRecHitToken, theStationShowerDeltaR, and theStationShowerTSize.

Referenced by fillShowerInformation().

                                                                 {
 
   // get all the necesary products
   event.getByToken(theDTRecHitToken, theDTRecHits);
   event.getByToken(theCSCRecHitToken, theCSCRecHits);
   event.getByToken(theCSCSegmentsToken, theCSCSegments);
   event.getByToken(theDT4DRecSegmentToken, 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 )

set the services needed

Definition at line 159 of file MuonShowerInformationFiller.cc.

References edm::EventSetup::get(), theCacheId_TRH, theCSCGeometry, theDTGeometry, theField, theMuonRecHitBuilder, theMuonRecHitBuilderName, theTracker, theTrackerRecHitBuilder, theTrackerRecHitBuilderName, and theTrackingGeometry.

Referenced by fillShowerInformation().

                                                                      {
 
   // DetLayer Geometry
   setup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry);
   setup.get<IdealMagneticFieldRecord>().get(theField);
   setup.get<TrackerRecoGeometryRecord>().get(theTracker);
   setup.get<MuonGeometryRecord>().get(theCSCGeometry);
   setup.get<MuonGeometryRecord>().get(theDTGeometry);
 
   // Transient Rechit Builders
   unsigned long long newCacheId_TRH = setup.get<TransientRecHitRecord>().cacheIdentifier();
   if ( newCacheId_TRH != theCacheId_TRH ) {
     setup.get<TransientRecHitRecord>().get(theTrackerRecHitBuilderName,theTrackerRecHitBuilder);
     setup.get<TransientRecHitRecord>().get(theMuonRecHitBuilderName,theMuonRecHitBuilder);
   }
 
 }