CMS 3D CMS Logo

Public Types | Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes

TrackDetectorAssociator Class Reference

#include <TrackDetectorAssociator.h>

List of all members.

Public Types

typedef TrackAssociatorParameters AssociatorParameters
enum  Direction { Any, InsideOut, OutsideIn }

Public Member Functions

TrackDetMatchInfo associate (const edm::Event &, const edm::EventSetup &, const FreeTrajectoryState &, const AssociatorParameters &)
TrackDetMatchInfo associate (const edm::Event &iEvent, const edm::EventSetup &iSetup, const AssociatorParameters &parameters, const FreeTrajectoryState *innerState, const FreeTrajectoryState *outerState=0)
TrackDetMatchInfo associate (const edm::Event &, const edm::EventSetup &, const reco::Track &, const AssociatorParameters &, Direction direction=Any)
 associate using reco::Track
TrackDetMatchInfo associate (const edm::Event &, const edm::EventSetup &, const GlobalVector &, const GlobalPoint &, const int, const AssociatorParameters &)
 associate using 3-momentum, vertex and charge
TrackDetMatchInfo associate (const edm::Event &, const edm::EventSetup &, const SimTrack &, const SimVertex &, const AssociatorParameters &)
 associate using a simulated track
void setPropagator (const Propagator *)
 use a user configured propagator
 TrackDetectorAssociator ()
void useDefaultPropagator ()
 use the default propagator
 ~TrackDetectorAssociator ()

Static Public Member Functions

static bool crossedIP (const reco::Track &track)
static FreeTrajectoryState getFreeTrajectoryState (const edm::EventSetup &, const SimTrack &, const SimVertex &)
static FreeTrajectoryState getFreeTrajectoryState (const edm::EventSetup &, const GlobalVector &, const GlobalPoint &, const int)
static FreeTrajectoryState getFreeTrajectoryState (const edm::EventSetup &, const reco::Track &)
 get FreeTrajectoryState from different track representations

Private Member Functions

bool addTAMuonSegmentMatch (TAMuonChamberMatch &, const RecSegment *, const AssociatorParameters &)
void fillCaloTowers (const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
void fillCaloTruth (const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
void fillEcal (const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
void fillHcal (const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
void fillHO (const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
void fillMuon (const edm::Event &, TrackDetMatchInfo &, const AssociatorParameters &)
void fillPreshower (const edm::Event &iEvent, TrackDetMatchInfo &info, const AssociatorParameters &)
DetIdAssociator::MapRange getMapRange (const std::pair< float, float > &delta, const float dR)
math::XYZPoint getPoint (const GlobalPoint &point)
math::XYZPoint getPoint (const LocalPoint &point)
void getTAMuonChamberMatches (std::vector< TAMuonChamberMatch > &matches, const AssociatorParameters &parameters)
math::XYZVector getVector (const GlobalVector &vec)
math::XYZVector getVector (const LocalVector &vec)
void init (const edm::EventSetup &)

Private Attributes

CachedTrajectory cachedTrajectory_
edm::ESHandle< DetIdAssociatorcaloDetIdAssociator_
PropagatordefProp_
edm::ESHandle< DetIdAssociatorecalDetIdAssociator_
edm::ESHandle< DetIdAssociatorhcalDetIdAssociator_
edm::ESHandle< DetIdAssociatorhoDetIdAssociator_
const PropagatorivProp_
edm::ESHandle< DetIdAssociatormuonDetIdAssociator_
edm::ESHandle< DetIdAssociatorpreshowerDetIdAssociator_
edm::ESHandle< CaloGeometrytheCaloGeometry_
edm::ESWatcher
< IdealMagneticFieldRecord
theMagneticFeildWatcher_
edm::ESHandle
< GlobalTrackingGeometry
theTrackingGeometry_
bool useDefaultPropagator_

Detailed Description

Definition at line 53 of file TrackDetectorAssociator.h.


Member Typedef Documentation

Definition at line 58 of file TrackDetectorAssociator.h.


Member Enumeration Documentation

Enumerator:
Any 
InsideOut 
OutsideIn 

Definition at line 59 of file TrackDetectorAssociator.h.


Constructor & Destructor Documentation

TrackDetectorAssociator::TrackDetectorAssociator ( )

Definition at line 100 of file TrackDetectorAssociator.cc.

{
   ivProp_ = 0;
   defProp_ = 0;
   useDefaultPropagator_ = false;
}
TrackDetectorAssociator::~TrackDetectorAssociator ( )

Definition at line 107 of file TrackDetectorAssociator.cc.

{
   if (defProp_) delete defProp_;
}

Member Function Documentation

bool TrackDetectorAssociator::addTAMuonSegmentMatch ( TAMuonChamberMatch matchedChamber,
const RecSegment segment,
const AssociatorParameters parameters 
) [private]

Definition at line 842 of file TrackDetectorAssociator.cc.

References reco::deltaPhi(), TrackAssociatorParameters::dRMuon, PV3DBase< T, PVType, FrameType >::eta(), TrajectoryStateOnSurface::freeState(), TrackingRecHit::geographicalId(), DTRecSegment4D::hasPhi(), TAMuonSegmentMatch::hasPhi, TAMuonSegmentMatch::hasZed, DTRecSegment4D::hasZed(), TAMuonChamberMatch::id, DTRecSegment2D::ist0Valid(), RecSegment::localDirection(), RecSegment::localDirectionError(), TrackingRecHit::localPosition(), TrackingRecHit::localPositionError(), LogTrace, M_PI, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), DTRecSegment4D::phiSegment(), FreeTrajectoryState::position(), funct::pow(), DetId::rawId(), TAMuonSegmentMatch::segmentGlobalPosition, TAMuonSegmentMatch::segmentLocalDirection, TAMuonSegmentMatch::segmentLocalErrorDxDz, TAMuonSegmentMatch::segmentLocalErrorDyDz, TAMuonSegmentMatch::segmentLocalErrorXDxDz, TAMuonSegmentMatch::segmentLocalErrorXX, TAMuonSegmentMatch::segmentLocalErrorXY, TAMuonSegmentMatch::segmentLocalErrorYDyDz, TAMuonSegmentMatch::segmentLocalErrorYY, TAMuonSegmentMatch::segmentLocalPosition, TAMuonChamberMatch::segments, DTRecSegment2D::specificRecHits(), mathSSE::sqrt(), DTRecSegment2D::t0(), TAMuonSegmentMatch::t0, GeomDet::toGlobal(), TAMuonChamberMatch::tState, LocalError::xx(), LocalError::xy(), and LocalError::yy().

{
   LogTrace("TrackAssociator")
     << "Segment local position: " << segment->localPosition() << "\n"
     << std::hex << segment->geographicalId().rawId() << "\n";
   
   const GeomDet* chamber = muonDetIdAssociator_->getGeomDet(matchedChamber.id);
   TrajectoryStateOnSurface trajectoryStateOnSurface = matchedChamber.tState;
   GlobalPoint segmentGlobalPosition = chamber->toGlobal(segment->localPosition());

   LogTrace("TrackAssociator")
     << "Segment global position: " << segmentGlobalPosition << " \t (R_xy,eta,phi): "
     << segmentGlobalPosition.perp() << "," << segmentGlobalPosition.eta() << "," << segmentGlobalPosition.phi() << "\n";

   LogTrace("TrackAssociator")
     << "\teta hit: " << segmentGlobalPosition.eta() << " \tpropagator: " << trajectoryStateOnSurface.freeState()->position().eta() << "\n"
     << "\tphi hit: " << segmentGlobalPosition.phi() << " \tpropagator: " << trajectoryStateOnSurface.freeState()->position().phi() << std::endl;
   
   bool isGood = false;
   bool isDTWithoutY = false;
   const DTRecSegment4D* dtseg = dynamic_cast<const DTRecSegment4D*>(segment);
   if ( dtseg && (! dtseg->hasZed()) )
      isDTWithoutY = true;

   double deltaPhi(fabs(segmentGlobalPosition.phi()-trajectoryStateOnSurface.freeState()->position().phi()));
   if(deltaPhi>M_PI) deltaPhi = fabs(deltaPhi-M_PI*2.);

   if( isDTWithoutY )
     {
        isGood = deltaPhi < parameters.dRMuon;
        // Be in chamber
        isGood &= fabs(segmentGlobalPosition.eta()-trajectoryStateOnSurface.freeState()->position().eta()) < .3;
     } else isGood = sqrt( pow(segmentGlobalPosition.eta()-trajectoryStateOnSurface.freeState()->position().eta(),2) + 
                           deltaPhi*deltaPhi) < parameters.dRMuon;

   if(isGood) {
      TAMuonSegmentMatch muonSegment;
      muonSegment.segmentGlobalPosition = getPoint(segmentGlobalPosition);
      muonSegment.segmentLocalPosition = getPoint( segment->localPosition() );
      muonSegment.segmentLocalDirection = getVector( segment->localDirection() );
      muonSegment.segmentLocalErrorXX = segment->localPositionError().xx();
      muonSegment.segmentLocalErrorYY = segment->localPositionError().yy();
      muonSegment.segmentLocalErrorXY = segment->localPositionError().xy();
      muonSegment.segmentLocalErrorDxDz = segment->localDirectionError().xx();
      muonSegment.segmentLocalErrorDyDz = segment->localDirectionError().yy();
      
      // DANGEROUS - compiler cannot guaranty parameters ordering
      // AlgebraicSymMatrix segmentCovMatrix = segment->parametersError();
      // muonSegment.segmentLocalErrorXDxDz = segmentCovMatrix[2][0];
      // muonSegment.segmentLocalErrorYDyDz = segmentCovMatrix[3][1];
      muonSegment.segmentLocalErrorXDxDz = -999;
      muonSegment.segmentLocalErrorYDyDz = -999;
      muonSegment.hasZed = true;
      muonSegment.hasPhi = true;
      
      // timing information
      muonSegment.t0 = 0;
      if ( dtseg ) {
        if ( (dtseg->hasPhi()) && (! isDTWithoutY) ) {
          int phiHits = dtseg->phiSegment()->specificRecHits().size();
          //      int zHits = dtseg->zSegment()->specificRecHits().size();
          int hits=0;
          double t0=0.;
          // TODO: cuts on hit numbers not optimized in any way yet...
          if (phiHits>5 && dtseg->phiSegment()->ist0Valid()) {
            t0+=dtseg->phiSegment()->t0()*phiHits;
            hits+=phiHits;
            LogTrace("TrackAssociator") << " Phi t0: " << dtseg->phiSegment()->t0() << " hits: " << phiHits;
          }
          // the z segments seem to contain little useful information...
//        if (zHits>3) {
//          t0+=s->zSegment()->t0()*zHits;
//          hits+=zHits;
//          std::cout << "   Z t0: " << s->zSegment()->t0() << " hits: " << zHits << std::endl;
//        }
          if (hits) muonSegment.t0 = t0/hits;
//        std::cout << " --- t0: " << muonSegment.t0 << std::endl;
        } else {
           // check and set dimensionality
           if (isDTWithoutY) muonSegment.hasZed = false;
           if (! dtseg->hasPhi()) muonSegment.hasPhi = false;
        }
      }
      matchedChamber.segments.push_back(muonSegment);
   }

   return isGood;
}
TrackDetMatchInfo TrackDetectorAssociator::associate ( const edm::Event iEvent,
const edm::EventSetup iSetup,
const FreeTrajectoryState fts,
const AssociatorParameters parameters 
)

propagate a track across the whole detector and find associated objects. Association is done in two modes 1) an object is associated to a track only if crossed by track 2) an object is associated to a track if it is withing an eta-phi cone of some radius with respect to a track. (the cone origin is at (0,0,0)) Trajectory bending in eta-phi is taking into account when matching is performed

associate using FreeTrajectoryState

Definition at line 157 of file TrackDetectorAssociator.cc.

Referenced by BetaCalculatorTK::addInfoToCandidate(), EcalCosmicsHists::analyze(), IsolatedTracksCone::analyze(), spr::chargeIsolation(), spr::chargeIsolationEcal(), spr::chargeIsolationHcal(), spr::coneChargeIsolation(), muonisolation::JetExtractor::deposit(), muonisolation::CaloExtractorByAssociator::deposits(), MuonIdProducer::fillMuonId(), SoftElectronProducer::produce(), HighPtTrackEcalDetIdProducer::produce(), InterestingTrackEcalDetIdProducer::produce(), cms::MuonMETValueMapProducer::produce(), AlCaIsoTracksProducer::produce(), and ReduceHcalRecHitCollectionProducer::produce().

{
   return associate(iEvent,iSetup,parameters,&fts);
}
TrackDetMatchInfo TrackDetectorAssociator::associate ( const edm::Event iEvent,
const edm::EventSetup iSetup,
const SimTrack track,
const SimVertex vertex,
const AssociatorParameters parameters 
)

associate using a simulated track

Definition at line 1087 of file TrackDetectorAssociator.cc.

{
   return associate(iEvent, iSetup, getFreeTrajectoryState(iSetup, track, vertex), parameters);
}
TrackDetMatchInfo TrackDetectorAssociator::associate ( const edm::Event iEvent,
const edm::EventSetup iSetup,
const AssociatorParameters parameters,
const FreeTrajectoryState innerState,
const FreeTrajectoryState outerState = 0 
)

associate using inner and outer most states of a track in the silicon tracker.

Definition at line 165 of file TrackDetectorAssociator.cc.

References Exception, ExpressReco_HICollisions_FallBack::fillCaloTowers, info, init, LogTrace, ExpressReco_HICollisions_FallBack::maxZ, ExpressReco_HICollisions_FallBack::minZ, PV3DBase< T, PVType, FrameType >::perp(), FreeTrajectoryState::position(), TrackDetMatchInfo::setCaloGeometry(), TrackDetMatchInfo::stateAtIP, TrackDetMatchInfo::trkGlobPosAtEcal, TrackDetMatchInfo::trkGlobPosAtHcal, TrackDetMatchInfo::trkGlobPosAtHO, TrackDetMatchInfo::trkMomAtEcal, TrackDetMatchInfo::trkMomAtHcal, TrackDetMatchInfo::trkMomAtHO, TrackAssociatorParameters::truthMatch, TrackAssociatorParameters::useCalo, TrackAssociatorParameters::useEcal, TrackAssociatorParameters::useHcal, TrackAssociatorParameters::useHO, TrackAssociatorParameters::useMuon, TrackAssociatorParameters::usePreshower, and PV3DBase< T, PVType, FrameType >::z().

{
   TrackDetMatchInfo info;
   if (! parameters.useEcal && ! parameters.useCalo && ! parameters.useHcal &&
       ! parameters.useHO && ! parameters.useMuon && !parameters.usePreshower)
     throw cms::Exception("ConfigurationError") << 
     "Configuration error! No subdetector was selected for the track association.";
   // TimerStack timers;
   // timers.push("TrackDetectorAssociator::associate",TimerStack::DetailedMonitoring);
   
   SteppingHelixStateInfo trackOrigin(*innerState);
   info.stateAtIP = *innerState;
   cachedTrajectory_.setStateAtIP(trackOrigin);
   
   init( iSetup );
   // get track trajectory
   // timers.push("TrackDetectorAssociator::fillEcal::propagation");
   // ECAL points (EB+EE)
   // If the phi angle between a track entrance and exit points is more
   // than 2 crystals, it is possible that the track will cross 3 crystals
   // and therefore one has to check at least 3 points along the track
   // trajectory inside ECAL. In order to have a chance to cross 4 crystalls
   // in the barrel, a track should have P_t as low as 3 GeV or smaller
   // If it's necessary, number of points along trajectory can be increased
   
   info.setCaloGeometry(theCaloGeometry_);
   
   // timers.push("TrackDetectorAssociator::associate::getTrajectories");
   cachedTrajectory_.reset_trajectory();
   // estimate propagation outer boundaries based on 
   // requested sub-detector information. For now limit
   // propagation region only if muon matching is not 
   // requested.
   double HOmaxR = hoDetIdAssociator_->volume().maxR();
   double HOmaxZ = hoDetIdAssociator_->volume().maxZ();
   double minR = ecalDetIdAssociator_->volume().minR();
   double minZ = preshowerDetIdAssociator_->volume().minZ();
   cachedTrajectory_.setMaxHORadius(HOmaxR);
   cachedTrajectory_.setMaxHOLength(HOmaxZ*2.);
   cachedTrajectory_.setMinDetectorRadius(minR);
   cachedTrajectory_.setMinDetectorLength(minZ*2.);

   double maxR(0);
   double maxZ(0);

   if (parameters.useMuon) {
     maxR = muonDetIdAssociator_->volume().maxR();
     maxZ = muonDetIdAssociator_->volume().maxZ();
     cachedTrajectory_.setMaxDetectorRadius(maxR);
     cachedTrajectory_.setMaxDetectorLength(maxZ*2.);
   }
   else {
     maxR = HOmaxR;
     maxZ = HOmaxZ;
     cachedTrajectory_.setMaxDetectorRadius(HOmaxR);
     cachedTrajectory_.setMaxDetectorLength(HOmaxZ*2.);
   }
   
   // If track extras exist and outerState is before HO maximum, then use outerState
   if (outerState) {
     if (outerState->position().perp()<HOmaxR && fabs(outerState->position().z())<HOmaxZ) {
       LogTrace("TrackAssociator") << "Using outerState as trackOrigin at Rho=" << outerState->position().perp()
             << "  Z=" << outerState->position().z() << "\n";
       trackOrigin = SteppingHelixStateInfo(*outerState);
     }
     else if(innerState) {
       LogTrace("TrackAssociator") << "Using innerState as trackOrigin at Rho=" << innerState->position().perp()
             << "  Z=" << innerState->position().z() << "\n";
       trackOrigin = SteppingHelixStateInfo(*innerState);
     }
   }

   if ( ! cachedTrajectory_.propagateAll(trackOrigin) ) return info;
   
   // get trajectory in calorimeters
   cachedTrajectory_.findEcalTrajectory( ecalDetIdAssociator_->volume() );
   cachedTrajectory_.findHcalTrajectory( hcalDetIdAssociator_->volume() );
   cachedTrajectory_.findHOTrajectory( hoDetIdAssociator_->volume() );
   cachedTrajectory_.findPreshowerTrajectory( preshowerDetIdAssociator_->volume() );

   info.trkGlobPosAtEcal = getPoint( cachedTrajectory_.getStateAtEcal().position() );
   info.trkGlobPosAtHcal = getPoint( cachedTrajectory_.getStateAtHcal().position() );
   info.trkGlobPosAtHO  = getPoint( cachedTrajectory_.getStateAtHO().position() );
   
   info.trkMomAtEcal = cachedTrajectory_.getStateAtEcal().momentum();
   info.trkMomAtHcal = cachedTrajectory_.getStateAtHcal().momentum();
   info.trkMomAtHO   = cachedTrajectory_.getStateAtHO().momentum();
   
   // timers.pop_and_push("TrackDetectorAssociator::associate::fillInfo");
     
   if (parameters.useEcal) fillEcal( iEvent, info, parameters);
   if (parameters.useCalo) fillCaloTowers( iEvent, info, parameters);
   if (parameters.useHcal) fillHcal( iEvent, info, parameters);
   if (parameters.useHO)   fillHO( iEvent, info, parameters);
   if (parameters.usePreshower) fillPreshower( iEvent, info, parameters);
   if (parameters.useMuon) fillMuon( iEvent, info, parameters);
   if (parameters.truthMatch) fillCaloTruth( iEvent, info, parameters);
   
   return info;
}
TrackDetMatchInfo TrackDetectorAssociator::associate ( const edm::Event iEvent,
const edm::EventSetup iSetup,
const GlobalVector momentum,
const GlobalPoint vertex,
const int  charge,
const AssociatorParameters parameters 
)

associate using 3-momentum, vertex and charge

Definition at line 91 of file PFBlockAlgo.cc.

References PFBlockAlgo::associate(), PFBlockAlgo::blocks_, gather_cfg::cout, PFBlockAlgo::debug_, PFBlockAlgo::elements_, PFBlockAlgo::link(), and NONE.

                                                     {

    
#ifdef PFLOW_DEBUG
  if(debug_ ) cout<<"PFBlockAlgo::associate start ----"<<endl;
#endif

  if( last!= elements_.end() ) {
    PFBlockLink::Type linktype = PFBlockLink::NONE;
    double dist = -1;
    PFBlock::LinkTest linktest = PFBlock::LINKTEST_RECHIT;
    link( *last, *next, linktype, linktest, dist ); 

   
    if(dist<-0.5) {
#ifdef PFLOW_DEBUG
      if(debug_ ) cout<<"link failed"<<endl;
#endif
      return ++next; // association failed
    }
    else {
      // add next element to the current pflowblock
      blocks_->back().addElement( *next );  

      // (*next)->setIndex( blocks_->back()->indexToLastElement() );
      
      // this is not necessary? 
      // next->setPFBlock(this);
      
      // create a link between next and last
      links.push_back( PFBlockLink(linktype, 
                                   linktest,
                                   dist,
                                   (*last)->index(), 
                                   (*next)->index() ) );
      // not necessary ?
      //       next->connect( links_.size()-1 );
      //       last->connect( links_.size()-1 );      
    }
  }
  else {
    // add next element to this eflowblock
#ifdef PFLOW_DEBUG
    if(debug_ ) cout<<"adding to block element "<<(**next)<<endl;
#endif
    blocks_->back().addElement( *next );
    // (*next)->setIndex( blocks_->back()->indexToLastElement() );   
    // next->setPFBlock(this);
  }

  // recursive call: associate next and other unused elements 
  
  //   IE afterNext = next;
  //   ++afterNext;
  //  cout<<"last "<<**last<<" next "<<**next<<endl;
  
  for(IE ie = elements_.begin(); 
      ie != elements_.end();) {
    
    if( ie == last || ie == next ) {
      ++ie;
      continue;
    } 
    
    // *ie already included to a block
    if( (*ie)->locked() ) {
#ifdef PFLOW_DEBUG
      if(debug_ ) cout<<"element "<<(**ie)<<"already used"<<endl;
#endif
      ++ie;
      continue;
    }    
    
    
#ifdef PFLOW_DEBUG
    if(debug_ ) cout<<"calling associate "<<(**next)<<" & "<<(**ie)<<endl;
#endif
    ie = associate(next, ie, links);
  }       

#ifdef PFLOW_DEBUG
  if(debug_ ) {
    cout<<"**** deleting element "<<endl;
    cout<<**next<<endl;
  }
#endif
  delete *next;

#ifdef PFLOW_DEBUG
  if(debug_ ) {
    cout<<"**** removing element "<<endl;
  }
#endif

  IE iteratorToNextFreeElement = elements_.erase( next );

#ifdef PFLOW_DEBUG
  if(debug_ ) cout<<"PFBlockAlgo::associate stop ----"<<endl;
#endif

  return iteratorToNextFreeElement;
}
TrackDetMatchInfo TrackDetectorAssociator::associate ( const edm::Event iEvent,
const edm::EventSetup iSetup,
const reco::Track track,
const AssociatorParameters parameters,
Direction  direction = Any 
)

associate using reco::Track

Definition at line 995 of file TrackDetectorAssociator.cc.

References ecalTB2006H4_GenSimDigiReco_cfg::bField, Exception, reco::Track::extra(), edm::EventSetup::get(), TrajectoryStateTransform::initialFreeState(), TrajectoryStateTransform::innerFreeState(), reco::Track::innerMomentum(), reco::Track::innerPosition(), edm::Ref< C, T, F >::isNull(), LogTrace, TrajectoryStateTransform::outerFreeState(), reco::Track::outerMomentum(), reco::Track::outerPosition(), PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), FreeTrajectoryState::position(), query::result, and PV3DBase< T, PVType, FrameType >::z().

{
   double currentStepSize = cachedTrajectory_.getPropagationStep();
   TrajectoryStateTransform tsTransform;
   edm::ESHandle<MagneticField> bField;
   iSetup.get<IdealMagneticFieldRecord>().get(bField);

   if(track.extra().isNull()) {
      if ( direction != InsideOut ) 
        throw cms::Exception("FatalError") << 
        "No TrackExtra information is available and association is done with something else than InsideOut track.\n" <<
        "Either change the parameter or provide needed data!\n";
     LogTrace("TrackAssociator") << "Track Extras not found\n";
     FreeTrajectoryState initialState = tsTransform.initialFreeState(track,&*bField);
     return associate(iEvent, iSetup, parameters, &initialState); // 5th argument is null pointer
   }
   
   LogTrace("TrackAssociator") << "Track Extras found\n";
   FreeTrajectoryState innerState = tsTransform.innerFreeState(track,&*bField);
   FreeTrajectoryState outerState = tsTransform.outerFreeState(track,&*bField);
   FreeTrajectoryState referenceState = tsTransform.initialFreeState(track,&*bField);
   
   LogTrace("TrackAssociator") << "inner track state (rho, z, phi):" << 
     track.innerPosition().Rho() << ", " << track.innerPosition().z() <<
     ", " << track.innerPosition().phi() << "\n";
   LogTrace("TrackAssociator") << "innerFreeState (rho, z, phi):" << 
     innerState.position().perp() << ", " << innerState.position().z() <<
     ", " << innerState.position().phi() << "\n";
   
   LogTrace("TrackAssociator") << "outer track state (rho, z, phi):" << 
     track.outerPosition().Rho() << ", " << track.outerPosition().z() <<
     ", " << track.outerPosition().phi() << "\n";
   LogTrace("TrackAssociator") << "outerFreeState (rho, z, phi):" << 
     outerState.position().perp() << ", " << outerState.position().z() <<
     ", " << outerState.position().phi() << "\n";
   
   // InsideOut first
   if ( crossedIP( track ) ) {
      switch ( direction ) {
       case InsideOut:
       case Any:
         return associate(iEvent, iSetup, parameters, &referenceState, &outerState);
         break;
       case OutsideIn:
           {
              cachedTrajectory_.setPropagationStep( -fabs(currentStepSize) );
              TrackDetMatchInfo result = associate(iEvent, iSetup, parameters, &innerState, &referenceState);
              cachedTrajectory_.setPropagationStep( currentStepSize );
              return result;
              break;
           }
      }
   } else {
      switch ( direction ) {
       case InsideOut:
         return associate(iEvent, iSetup, parameters, &innerState, &outerState);
         break;
       case OutsideIn:
           {
              cachedTrajectory_.setPropagationStep( -fabs(currentStepSize) );
              TrackDetMatchInfo result = associate(iEvent, iSetup, parameters, &outerState, &innerState);
              cachedTrajectory_.setPropagationStep( currentStepSize );
              return result;
              break;
           }
       case Any:
           {
              // check if we deal with clear outside-in case
              if ( track.innerPosition().Dot( track.innerMomentum() ) < 0 &&
                   track.outerPosition().Dot( track.outerMomentum() ) < 0 )
                {
                   cachedTrajectory_.setPropagationStep( -fabs(currentStepSize) );
                   TrackDetMatchInfo result;
                   if ( track.innerPosition().R() < track.outerPosition().R() )
                     result = associate(iEvent, iSetup, parameters, &innerState, &outerState);
                   else
                     result = associate(iEvent, iSetup, parameters, &outerState, &innerState);
                   cachedTrajectory_.setPropagationStep( currentStepSize );
                   return result;
                }
           }
      }
   }
        
   // all other cases  
   return associate(iEvent, iSetup, parameters, &innerState, &outerState);
}
bool TrackDetectorAssociator::crossedIP ( const reco::Track track) [static]

Definition at line 1106 of file TrackDetectorAssociator.cc.

References reco::Track::innerMomentum(), reco::Track::innerPosition(), reco::Track::outerMomentum(), and reco::Track::outerPosition().

Referenced by MuonIdProducer::produce().

{
   bool crossed = true;
   crossed &= (track.innerPosition().rho() > 3 ); // something close to active volume
   crossed &= (track.outerPosition().rho() > 3 ); // something close to active volume
   crossed &=  ( ( track.innerPosition().x()*track.innerMomentum().x() +
                   track.innerPosition().y()*track.innerMomentum().y() < 0 ) !=
                 ( track.outerPosition().x()*track.outerMomentum().x() + 
                   track.outerPosition().y()*track.outerMomentum().y() < 0 ) );
   return crossed;
}
void TrackDetectorAssociator::fillCaloTowers ( const edm::Event iEvent,
TrackDetMatchInfo info,
const AssociatorParameters parameters 
) [private]

Definition at line 353 of file TrackDetectorAssociator.cc.

References TrackAssociatorParameters::accountForTrajectoryChangeCalo, ExpressReco_HICollisions_FallBack::caloTowers, TrackDetMatchInfo::crossedTowerIds, TrackDetMatchInfo::crossedTowers, TrackAssociatorParameters::dRHcal, TrackAssociatorParameters::dRHcalPreselection, edm::Event::getByLabel(), CachedTrajectory::IpToHcal, TrackDetMatchInfo::isGoodCalo, edm::HandleBase::isValid(), LogTrace, TrackAssociatorParameters::theCaloTowerCollectionLabel, and TrackDetMatchInfo::towers.

{
   // TimerStack timers;
   // timers.push("TrackDetectorAssociator::fillCaloTowers");

   // use ECAL and HCAL trajectories to match a tower. (HO isn't used for matching).
   std::vector<GlobalPoint> trajectory;
   const std::vector<SteppingHelixStateInfo>& ecalTrajectoryStates = cachedTrajectory_.getEcalTrajectory();
   const std::vector<SteppingHelixStateInfo>& hcalTrajectoryStates = cachedTrajectory_.getHcalTrajectory();
   for(std::vector<SteppingHelixStateInfo>::const_iterator itr = ecalTrajectoryStates.begin();
       itr != ecalTrajectoryStates.end(); itr++) trajectory.push_back(itr->position());
   for(std::vector<SteppingHelixStateInfo>::const_iterator itr = hcalTrajectoryStates.begin();
       itr != hcalTrajectoryStates.end(); itr++) trajectory.push_back(itr->position());
   
   if(trajectory.empty()) {
      LogTrace("TrackAssociator") << "HCAL trajectory is empty; moving on\n";
      info.isGoodCalo = 0;
      return;
   }
   info.isGoodCalo = 1;
   
   // find crossed CaloTowers
   // timers.push("TrackDetectorAssociator::fillCaloTowers::access::CaloTowers");
   edm::Handle<CaloTowerCollection> caloTowers;

   iEvent.getByLabel( parameters.theCaloTowerCollectionLabel, caloTowers );
   if (!caloTowers.isValid())  throw cms::Exception("FatalError") << "Unable to find CaloTowers in event!\n";
   
   // timers.pop_and_push("TrackDetectorAssociator::fillCaloTowers::matching");
   std::set<DetId> caloTowerIdsInRegion;
   if (parameters.accountForTrajectoryChangeCalo){
      // get trajectory change with respect to initial state
      DetIdAssociator::MapRange mapRange = getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHcal),
                                                       parameters.dRHcalPreselection);
      caloTowerIdsInRegion = caloDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0],mapRange);
   } else caloTowerIdsInRegion = caloDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0], parameters.dRHcalPreselection);
   
   LogTrace("TrackAssociator") << "Towers in the region: " << caloTowerIdsInRegion.size();
   std::set<DetId> caloTowerIdsInACone = caloDetIdAssociator_->getDetIdsInACone(caloTowerIdsInRegion, trajectory, parameters.dRHcal);
   LogTrace("TrackAssociator") << "Towers in the cone: " << caloTowerIdsInACone.size();
   std::vector<DetId> crossedCaloTowerIds = caloDetIdAssociator_->getCrossedDetIds(caloTowerIdsInRegion, trajectory);
   LogTrace("TrackAssociator") << "Towers crossed: " << crossedCaloTowerIds.size();
   
   info.crossedTowerIds = crossedCaloTowerIds;
   
   // add CaloTowers
   // timers.pop_and_push("TrackDetectorAssociator::fillCaloTowers::addCrossedTowers");
   for(std::vector<DetId>::const_iterator itr=crossedCaloTowerIds.begin(); itr!=crossedCaloTowerIds.end();itr++)
     {
        CaloTowerCollection::const_iterator tower = (*caloTowers).find(*itr);
        if(tower != (*caloTowers).end()) 
          info.crossedTowers.push_back(&*tower);
        else
          LogTrace("TrackAssociator") << "Crossed CaloTower is not found for DetId: " << (*itr).rawId();
     }
   
   // timers.pop_and_push("TrackDetectorAssociator::fillCaloTowers::addTowersInTheRegion");
   for(std::set<DetId>::const_iterator itr=caloTowerIdsInACone.begin(); itr!=caloTowerIdsInACone.end();itr++)
     {
        CaloTowerCollection::const_iterator tower = (*caloTowers).find(*itr);
        if(tower != (*caloTowers).end()) 
          info.towers.push_back(&*tower);
        else 
          LogTrace("TrackAssociator") << "CaloTower from the cone is not found for DetId: " << (*itr).rawId();
     }
   
}
void TrackDetectorAssociator::fillCaloTruth ( const edm::Event iEvent,
TrackDetMatchInfo info,
const AssociatorParameters parameters 
) [private]

Definition at line 935 of file TrackDetectorAssociator.cc.

References TrackDetMatchInfo::ecalTrueEnergy, SimTrack::genpartIndex(), edm::Event::getByLabel(), edm::Event::getByType(), TrackDetMatchInfo::hcalTrueEnergy, TrackDetMatchInfo::hcalTrueEnergyCorrected, FreeTrajectoryState::momentum(), TrackDetMatchInfo::simTrack, TrackDetMatchInfo::stateAtIP, TrackDetMatchInfo::trkGlobPosAtHcal, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

{
   // get list of simulated tracks and their vertices
   using namespace edm;
   Handle<SimTrackContainer> simTracks;
   iEvent.getByType<SimTrackContainer>(simTracks);
   if (! simTracks.isValid() ) throw cms::Exception("FatalError") << "No simulated tracks found\n";
   
   Handle<SimVertexContainer> simVertices;
   iEvent.getByType<SimVertexContainer>(simVertices);
   if (! simVertices.isValid() ) throw cms::Exception("FatalError") << "No simulated vertices found\n";
   
   // get sim calo hits
   Handle<PCaloHitContainer> simEcalHitsEB;
   iEvent.getByLabel("g4SimHits","EcalHitsEB",simEcalHitsEB);
   if (! simEcalHitsEB.isValid() ) throw cms::Exception("FatalError") << "No simulated ECAL EB hits found\n";

   Handle<PCaloHitContainer> simEcalHitsEE;
   iEvent.getByLabel("g4SimHits","EcalHitsEE",simEcalHitsEE);
   if (! simEcalHitsEE.isValid() ) throw cms::Exception("FatalError") << "No simulated ECAL EE hits found\n";

   Handle<PCaloHitContainer> simHcalHits;
   iEvent.getByLabel("g4SimHits","HcalHits",simHcalHits);
   if (! simHcalHits.isValid() ) throw cms::Exception("FatalError") << "No simulated HCAL hits found\n";

   // find truth partner
   SimTrackContainer::const_iterator simTrack = simTracks->begin();
   for( ; simTrack != simTracks->end(); ++simTrack){
      math::XYZVector simP3( simTrack->momentum().x(), simTrack->momentum().y(), simTrack->momentum().z() );
      math::XYZVector recoP3( info.stateAtIP.momentum().x(), info.stateAtIP.momentum().y(), info.stateAtIP.momentum().z() );
      if ( ROOT::Math::VectorUtil::DeltaR(recoP3, simP3) < 0.1 ) break;
   }
   if ( simTrack != simTracks->end() ) {
      info.simTrack = &(*simTrack);
      double ecalTrueEnergy(0);
      double hcalTrueEnergy(0);
      
      // loop over calo hits
      for( PCaloHitContainer::const_iterator hit = simEcalHitsEB->begin(); hit != simEcalHitsEB->end(); ++hit )
        if ( hit->geantTrackId() == info.simTrack->genpartIndex() ) ecalTrueEnergy += hit->energy();
      
      for( PCaloHitContainer::const_iterator hit = simEcalHitsEE->begin(); hit != simEcalHitsEE->end(); ++hit )
        if ( hit->geantTrackId() == info.simTrack->genpartIndex() ) ecalTrueEnergy += hit->energy();
      
      for( PCaloHitContainer::const_iterator hit = simHcalHits->begin(); hit != simHcalHits->end(); ++hit )
        if ( hit->geantTrackId() == info.simTrack->genpartIndex() ) hcalTrueEnergy += hit->energy();
      
      info.ecalTrueEnergy = ecalTrueEnergy;
      info.hcalTrueEnergy = hcalTrueEnergy;
      info.hcalTrueEnergyCorrected = hcalTrueEnergy;
      if ( fabs(info.trkGlobPosAtHcal.eta()) < 1.3 )
        info.hcalTrueEnergyCorrected = hcalTrueEnergy*113.2;
      else 
        if ( fabs(info.trkGlobPosAtHcal.eta()) < 3.0 )
          info.hcalTrueEnergyCorrected = hcalTrueEnergy*167.2;
   }
}
void TrackDetectorAssociator::fillEcal ( const edm::Event iEvent,
TrackDetMatchInfo info,
const AssociatorParameters parameters 
) [private]

Definition at line 270 of file TrackDetectorAssociator.cc.

References TrackAssociatorParameters::accountForTrajectoryChangeCalo, TrackDetMatchInfo::crossedEcalIds, TrackDetMatchInfo::crossedEcalRecHits, TrackAssociatorParameters::dREcal, TrackAssociatorParameters::dREcalPreselection, egHLT::errCodes::EBRecHits, TrackDetMatchInfo::ecalRecHits, egHLT::errCodes::EERecHits, edm::Event::getByLabel(), CachedTrajectory::IpToEcal, TrackDetMatchInfo::isGoodEcal, edm::HandleBase::isValid(), LogTrace, TrackAssociatorParameters::theEBRecHitCollectionLabel, and TrackAssociatorParameters::theEERecHitCollectionLabel.

{
   // TimerStack timers;
   // timers.push("TrackDetectorAssociator::fillEcal");
   
   const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getEcalTrajectory();
   
   for(std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
       itr != trajectoryStates.end(); itr++) 
     LogTrace("TrackAssociator") << "ECAL trajectory point (rho, z, phi): " << itr->position().perp() <<
     ", " << itr->position().z() << ", " << itr->position().phi();

   std::vector<GlobalPoint> coreTrajectory;
   for(std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
       itr != trajectoryStates.end(); itr++) coreTrajectory.push_back(itr->position());
   
   if(coreTrajectory.empty()) {
      LogTrace("TrackAssociator") << "ECAL track trajectory is empty; moving on\n";
      info.isGoodEcal = 0;
      return;
   }
   info.isGoodEcal = 1;

   // Find ECAL crystals
   // timers.push("TrackDetectorAssociator::fillEcal::access::EcalBarrel");
   edm::Handle<EBRecHitCollection> EBRecHits;
   iEvent.getByLabel( parameters.theEBRecHitCollectionLabel, EBRecHits );
   if (!EBRecHits.isValid()) throw cms::Exception("FatalError") << "Unable to find EBRecHitCollection in the event!\n";

   // timers.pop_and_push("TrackDetectorAssociator::fillEcal::access::EcalEndcaps");
   edm::Handle<EERecHitCollection> EERecHits;
   iEvent.getByLabel( parameters.theEERecHitCollectionLabel, EERecHits );
   if (!EERecHits.isValid()) throw cms::Exception("FatalError") << "Unable to find EERecHitCollection in event!\n";

   // timers.pop_and_push("TrackDetectorAssociator::fillEcal::matching");
   // timers.push("TrackDetectorAssociator::fillEcal::matching::region");
   std::set<DetId> ecalIdsInRegion;
   if (parameters.accountForTrajectoryChangeCalo){
      // get trajectory change with respect to initial state
      DetIdAssociator::MapRange mapRange = getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToEcal),
                                                       parameters.dREcalPreselection);
      ecalIdsInRegion = ecalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0],mapRange);
   } else ecalIdsInRegion = ecalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dREcalPreselection);
   // timers.pop_and_push("TrackDetectorAssociator::fillEcal::matching::cone");
   LogTrace("TrackAssociator") << "ECAL hits in the region: " << ecalIdsInRegion.size();
   if (parameters.dREcalPreselection > parameters.dREcal)
     ecalIdsInRegion =  ecalDetIdAssociator_->getDetIdsInACone(ecalIdsInRegion, coreTrajectory, parameters.dREcal);
   LogTrace("TrackAssociator") << "ECAL hits in the cone: " << ecalIdsInRegion.size();
   std::vector<DetId> crossedEcalIds =  
     ecalDetIdAssociator_->getCrossedDetIds(ecalIdsInRegion, coreTrajectory);
   LogTrace("TrackAssociator") << "ECAL crossed hits " << crossedEcalIds.size();
   
   info.crossedEcalIds = crossedEcalIds;
   
   // add EcalRecHits
   // timers.pop_and_push("TrackDetectorAssociator::fillEcal::addCrossedHits");
   for(std::vector<DetId>::const_iterator itr=crossedEcalIds.begin(); itr!=crossedEcalIds.end();itr++)
   {
      std::vector<EcalRecHit>::const_iterator ebHit = (*EBRecHits).find(*itr);
      std::vector<EcalRecHit>::const_iterator eeHit = (*EERecHits).find(*itr);
      if(ebHit != (*EBRecHits).end()) 
         info.crossedEcalRecHits.push_back(&*ebHit);
      else if(eeHit != (*EERecHits).end()) 
         info.crossedEcalRecHits.push_back(&*eeHit);
      else  
         LogTrace("TrackAssociator") << "Crossed EcalRecHit is not found for DetId: " << itr->rawId();
   }
   // timers.pop_and_push("TrackDetectorAssociator::fillEcal::addHitsInTheRegion");
   for(std::set<DetId>::const_iterator itr=ecalIdsInRegion.begin(); itr!=ecalIdsInRegion.end();itr++)
   {
      std::vector<EcalRecHit>::const_iterator ebHit = (*EBRecHits).find(*itr);
      std::vector<EcalRecHit>::const_iterator eeHit = (*EERecHits).find(*itr);
      if(ebHit != (*EBRecHits).end()) 
         info.ecalRecHits.push_back(&*ebHit);
      else if(eeHit != (*EERecHits).end()) 
         info.ecalRecHits.push_back(&*eeHit);
      else 
         LogTrace("TrackAssociator") << "EcalRecHit from the cone is not found for DetId: " << itr->rawId();
   }
}
void TrackDetectorAssociator::fillHcal ( const edm::Event iEvent,
TrackDetMatchInfo info,
const AssociatorParameters parameters 
) [private]

Definition at line 448 of file TrackDetectorAssociator.cc.

References TrackAssociatorParameters::accountForTrajectoryChangeCalo, runEdmFileComparison::collection, TrackDetMatchInfo::crossedHcalIds, TrackDetMatchInfo::crossedHcalRecHits, TrackAssociatorParameters::dRHcal, TrackAssociatorParameters::dRHcalPreselection, edm::Event::getByLabel(), TrackDetMatchInfo::hcalRecHits, info, CachedTrajectory::IpToHcal, TrackDetMatchInfo::isGoodHcal, edm::HandleBase::isValid(), LogTrace, and TrackAssociatorParameters::theHBHERecHitCollectionLabel.

{
   // TimerStack timers;
   // timers.push("TrackDetectorAssociator::fillHcal");

   const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getHcalTrajectory();

   std::vector<GlobalPoint> coreTrajectory;
   for(std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
       itr != trajectoryStates.end(); itr++) coreTrajectory.push_back(itr->position());   
   
   if(coreTrajectory.empty()) {
      LogTrace("TrackAssociator") << "HCAL trajectory is empty; moving on\n";
      info.isGoodHcal = 0;
      return;
   }
   info.isGoodHcal = 1;
   
   // find crossed Hcals
   // timers.push("TrackDetectorAssociator::fillHcal::access::Hcal");
   edm::Handle<HBHERecHitCollection> collection;

   iEvent.getByLabel( parameters.theHBHERecHitCollectionLabel, collection );
   if ( ! collection.isValid() ) throw cms::Exception("FatalError") << "Unable to find HBHERecHits in event!\n";
   
   // timers.pop_and_push("TrackDetectorAssociator::fillHcal::matching");
   std::set<DetId> idsInRegion;
   if (parameters.accountForTrajectoryChangeCalo){
      // get trajectory change with respect to initial state
      DetIdAssociator::MapRange mapRange = getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHcal),
                                                       parameters.dRHcalPreselection);
      idsInRegion = hcalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], mapRange);
   } else idsInRegion = hcalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dRHcalPreselection);
   
   LogTrace("TrackAssociator") << "HCAL hits in the region: " << idsInRegion.size() << "\n" << DetIdInfo::info(idsInRegion);
   std::set<DetId> idsInACone = hcalDetIdAssociator_->getDetIdsInACone(idsInRegion, coreTrajectory, parameters.dRHcal);
   LogTrace("TrackAssociator") << "HCAL hits in the cone: " << idsInACone.size() << "\n" << DetIdInfo::info(idsInACone);
   std::vector<DetId> crossedIds = 
     hcalDetIdAssociator_->getCrossedDetIds(idsInRegion, coreTrajectory);
   LogTrace("TrackAssociator") << "HCAL hits crossed: " << crossedIds.size() << "\n" << DetIdInfo::info(crossedIds);
   
   info.crossedHcalIds = crossedIds;
   // timers.pop_and_push("TrackDetectorAssociator::fillHcal::addCrossedHits");
   // add Hcal
   // timers.push("TrackDetectorAssociator::fillHcal::addHcal");
   for(std::vector<DetId>::const_iterator itr=crossedIds.begin(); itr!=crossedIds.end();itr++)
     {
        HBHERecHitCollection::const_iterator hit = (*collection).find(*itr);
        if( hit != (*collection).end() ) 
          info.crossedHcalRecHits.push_back(&*hit);
        else
          LogTrace("TrackAssociator") << "Crossed HBHERecHit is not found for DetId: " << itr->rawId();
     }
   // timers.pop_and_push("TrackDetectorAssociator::fillHcal::addHitsInTheRegion");
   for(std::set<DetId>::const_iterator itr=idsInACone.begin(); itr!=idsInACone.end();itr++)
     {
        HBHERecHitCollection::const_iterator hit = (*collection).find(*itr);
        if( hit != (*collection).end() ) 
          info.hcalRecHits.push_back(&*hit);
        else 
          LogTrace("TrackAssociator") << "HBHERecHit from the cone is not found for DetId: " << itr->rawId();
     }
}
void TrackDetectorAssociator::fillHO ( const edm::Event iEvent,
TrackDetMatchInfo info,
const AssociatorParameters parameters 
) [private]

Definition at line 514 of file TrackDetectorAssociator.cc.

References TrackAssociatorParameters::accountForTrajectoryChangeCalo, runEdmFileComparison::collection, TrackDetMatchInfo::crossedHOIds, TrackDetMatchInfo::crossedHORecHits, TrackAssociatorParameters::dRHcal, TrackAssociatorParameters::dRHcalPreselection, edm::Event::getByLabel(), TrackDetMatchInfo::hoRecHits, CachedTrajectory::IpToHO, TrackDetMatchInfo::isGoodHO, edm::HandleBase::isValid(), LogTrace, and TrackAssociatorParameters::theHORecHitCollectionLabel.

{
   // TimerStack timers;
   // timers.push("TrackDetectorAssociator::fillHO");

   const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getHOTrajectory();

   std::vector<GlobalPoint> coreTrajectory;
   for(std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
       itr != trajectoryStates.end(); itr++) coreTrajectory.push_back(itr->position());

   if(coreTrajectory.empty()) {
      LogTrace("TrackAssociator") << "HO trajectory is empty; moving on\n";
      info.isGoodHO = 0;
      return;
   }
   info.isGoodHO = 1;
   
   // find crossed HOs
   // timers.pop_and_push("TrackDetectorAssociator::fillHO::access::HO");
   edm::Handle<HORecHitCollection> collection;

   iEvent.getByLabel( parameters.theHORecHitCollectionLabel, collection );
   if ( ! collection.isValid() ) throw cms::Exception("FatalError") << "Unable to find HORecHits in event!\n";
   
   // timers.pop_and_push("TrackDetectorAssociator::fillHO::matching");
   std::set<DetId> idsInRegion;
   if (parameters.accountForTrajectoryChangeCalo){
      // get trajectory change with respect to initial state
      DetIdAssociator::MapRange mapRange = getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHO),
                                                       parameters.dRHcalPreselection);
      idsInRegion = hoDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], mapRange);
   } else idsInRegion = hoDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dRHcalPreselection);
   
   LogTrace("TrackAssociator") << "idsInRegion.size(): " << idsInRegion.size();
   std::set<DetId> idsInACone = hoDetIdAssociator_->getDetIdsInACone(idsInRegion, coreTrajectory, parameters.dRHcal);
   LogTrace("TrackAssociator") << "idsInACone.size(): " << idsInACone.size();
   std::vector<DetId> crossedIds =
     hoDetIdAssociator_->getCrossedDetIds(idsInRegion, coreTrajectory);
   
   info.crossedHOIds = crossedIds;
   
   // add HO
   // timers.pop_and_push("TrackDetectorAssociator::fillHO::addCrossedHits");
   for(std::vector<DetId>::const_iterator itr=crossedIds.begin(); itr!=crossedIds.end();itr++)
     {
        HORecHitCollection::const_iterator hit = (*collection).find(*itr);
        if( hit != (*collection).end() ) 
          info.crossedHORecHits.push_back(&*hit);
        else
          LogTrace("TrackAssociator") << "Crossed HORecHit is not found for DetId: " << itr->rawId();
     }

   // timers.pop_and_push("TrackDetectorAssociator::fillHO::addHitsInTheRegion");
   for(std::set<DetId>::const_iterator itr=idsInACone.begin(); itr!=idsInACone.end();itr++)
     {
        HORecHitCollection::const_iterator hit = (*collection).find(*itr);
        if( hit != (*collection).end() ) 
          info.hoRecHits.push_back(&*hit);
        else 
          LogTrace("TrackAssociator") << "HORecHit from the cone is not found for DetId: " << itr->rawId();
     }
}
void TrackDetectorAssociator::fillMuon ( const edm::Event iEvent,
TrackDetMatchInfo info,
const AssociatorParameters parameters 
) [private]

Definition at line 768 of file TrackDetectorAssociator.cc.

References TrackDetMatchInfo::chambers, ExpressReco_HICollisions_FallBack::cscSegments, Exception, edm::Event::getByLabel(), edm::HandleBase::isValid(), LogTrace, TrackAssociatorParameters::theCSCSegmentCollectionLabel, and TrackAssociatorParameters::theDTRecSegment4DCollectionLabel.

{
   // TimerStack timers;
   // timers.push("TrackDetectorAssociator::fillMuon");

   // Get the segments from the event
   // timers.push("TrackDetectorAssociator::fillMuon::access");
   edm::Handle<DTRecSegment4DCollection> dtSegments;
   iEvent.getByLabel( parameters.theDTRecSegment4DCollectionLabel, dtSegments );
   if (! dtSegments.isValid()) 
     throw cms::Exception("FatalError") << "Unable to find DTRecSegment4DCollection in event!\n";
   
   edm::Handle<CSCSegmentCollection> cscSegments;
   iEvent.getByLabel( parameters.theCSCSegmentCollectionLabel, cscSegments );
   if (! cscSegments.isValid()) 
     throw cms::Exception("FatalError") << "Unable to find CSCSegmentCollection in event!\n";

   
   // check the map of available segments
   // if there is no segments in a given direction at all,
   // then there is no point to fly there.
   // 
   // MISSING
   // Possible solution: quick search for presence of segments 
   // for the set of DetIds

   // timers.pop_and_push("TrackDetectorAssociator::fillMuon::matchChembers");
   
   // get a set of matches corresponding to muon chambers
   std::vector<TAMuonChamberMatch> matchedChambers;
   LogTrace("TrackAssociator") << "Trying to Get ChamberMatches" << std::endl;
   getTAMuonChamberMatches(matchedChambers, parameters);
   LogTrace("TrackAssociator") << "Chambers matched: " << matchedChambers.size() << "\n";
   
   // Iterate over all chamber matches and fill segment matching 
   // info if it's available
   // timers.pop_and_push("TrackDetectorAssociator::fillMuon::findSemgents");
   for(std::vector<TAMuonChamberMatch>::iterator matchedChamber = matchedChambers.begin(); 
       matchedChamber != matchedChambers.end(); matchedChamber++)
     {
        const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet((*matchedChamber).id);
        // DT chamber
        if(const DTChamber* chamber = dynamic_cast<const DTChamber*>(geomDet) ) {
           // Get the range for the corresponding segments
           DTRecSegment4DCollection::range  range = dtSegments->get(chamber->id());
           // Loop over the segments of this chamber
           for (DTRecSegment4DCollection::const_iterator segment = range.first; segment!=range.second; segment++) {
             if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) {
                matchedChamber->segments.back().dtSegmentRef = DTRecSegment4DRef(dtSegments, segment - dtSegments->begin());
             }
           }
        }else{
           // CSC Chamber
           if(const CSCChamber* chamber = dynamic_cast<const CSCChamber*>(geomDet) ) {
              // Get the range for the corresponding segments
              CSCSegmentCollection::range  range = cscSegments->get(chamber->id());
              // Loop over the segments
              for (CSCSegmentCollection::const_iterator segment = range.first; segment!=range.second; segment++) {
                 if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) {
                     matchedChamber->segments.back().cscSegmentRef = CSCSegmentRef(cscSegments, segment - cscSegments->begin());
                 }
              }
           }else{
              throw cms::Exception("FatalError") << "Failed to cast GeomDet object to either DTChamber or CSCChamber. Who is this guy anyway?\n";
           }
        }
        info.chambers.push_back(*matchedChamber);
     }
}
void TrackDetectorAssociator::fillPreshower ( const edm::Event iEvent,
TrackDetMatchInfo info,
const AssociatorParameters parameters 
) [private]

Definition at line 423 of file TrackDetectorAssociator.cc.

References TrackDetMatchInfo::crossedPreshowerIds, TrackAssociatorParameters::dRPreshowerPreselection, and LogTrace.

{
   std::vector<GlobalPoint> trajectory;
   const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getPreshowerTrajectory();
   for(std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
       itr != trajectoryStates.end(); itr++) trajectory.push_back(itr->position());
   
   if(trajectory.empty()) {
      LogTrace("TrackAssociator") << "Preshower trajectory is empty; moving on\n";
      return;
   }
   
   std::set<DetId> idsInRegion = 
     preshowerDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0], 
                                                       parameters.dRPreshowerPreselection);
   
   LogTrace("TrackAssociator") << "Number of Preshower Ids in the region: " << idsInRegion.size();
   std::vector<DetId> crossedIds = preshowerDetIdAssociator_->getCrossedDetIds(idsInRegion, trajectory);
   LogTrace("TrackAssociator") << "Number of Preshower Ids in crossed: " << crossedIds.size();
   info.crossedPreshowerIds = crossedIds;
}
FreeTrajectoryState TrackDetectorAssociator::getFreeTrajectoryState ( const edm::EventSetup iSetup,
const reco::Track track 
) [static]

get FreeTrajectoryState from different track representations

Definition at line 612 of file TrackDetectorAssociator.cc.

References ecalTB2006H4_GenSimDigiReco_cfg::bField, reco::TrackBase::charge(), ExpressReco_HICollisions_FallBack::e, edm::EventSetup::get(), reco::TrackBase::momentum(), point, and reco::TrackBase::vertex().

Referenced by BetaCalculatorTK::addInfoToCandidate(), IsolatedTracksCone::analyze(), spr::chargeIsolation(), spr::chargeIsolationEcal(), spr::chargeIsolationHcal(), spr::coneChargeIsolation(), SoftElectronProducer::produce(), InterestingTrackEcalDetIdProducer::produce(), cms::MuonMETValueMapProducer::produce(), and AlCaIsoTracksProducer::produce().

{
   edm::ESHandle<MagneticField> bField;
   iSetup.get<IdealMagneticFieldRecord>().get(bField);
   
   GlobalVector vector( track.momentum().x(), track.momentum().y(), track.momentum().z() );

   GlobalPoint point( track.vertex().x(), track.vertex().y(),  track.vertex().z() );

   GlobalTrajectoryParameters tPars(point, vector, track.charge(), &*bField);
   
   // FIX THIS !!!
   // need to convert from perigee to global or helix (curvilinear) frame
   // for now just an arbitrary matrix.
   CLHEP::HepSymMatrix covT(6,1); covT *= 1e-6; // initialize to sigma=1e-3
   CartesianTrajectoryError tCov(covT);
   
   return FreeTrajectoryState(tPars, tCov);
}
FreeTrajectoryState TrackDetectorAssociator::getFreeTrajectoryState ( const edm::EventSetup iSetup,
const SimTrack track,
const SimVertex vertex 
) [static]

Definition at line 580 of file TrackDetectorAssociator.cc.

References abs, DeDxDiscriminatorTools::charge(), CoreSimTrack::momentum(), point, CoreSimVertex::position(), and CoreSimTrack::type().

{
   GlobalVector vector( track.momentum().x(), track.momentum().y(), track.momentum().z() );
   GlobalPoint point( vertex.position().x(), vertex.position().y(), vertex.position().z() );

   int charge = track.type( )> 0 ? -1 : 1; // lepton convention
   if ( abs(track.type( )) == 211 || // pion
        abs(track.type( )) == 321 || // kaon
        abs(track.type( )) == 2212 )
     charge = track.type( )< 0 ? -1 : 1;
   return getFreeTrajectoryState(iSetup, vector, point, charge);
}
FreeTrajectoryState TrackDetectorAssociator::getFreeTrajectoryState ( const edm::EventSetup iSetup,
const GlobalVector momentum,
const GlobalPoint vertex,
const int  charge 
) [static]

Definition at line 595 of file TrackDetectorAssociator.cc.

References ecalTB2006H4_GenSimDigiReco_cfg::bField, ExpressReco_HICollisions_FallBack::e, and edm::EventSetup::get().

{
   edm::ESHandle<MagneticField> bField;
   iSetup.get<IdealMagneticFieldRecord>().get(bField);
   
   GlobalTrajectoryParameters tPars(vertex, momentum, charge, &*bField);
   
   CLHEP::HepSymMatrix covT(6,1); covT *= 1e-6; // initialize to sigma=1e-3
   CartesianTrajectoryError tCov(covT);
   
   return FreeTrajectoryState(tPars, tCov);
}
DetIdAssociator::MapRange TrackDetectorAssociator::getMapRange ( const std::pair< float, float > &  delta,
const float  dR 
) [private]

Definition at line 633 of file TrackDetectorAssociator.cc.

References DetIdAssociator::MapRange::dPhiMinus, DetIdAssociator::MapRange::dPhiPlus, DetIdAssociator::MapRange::dThetaMinus, DetIdAssociator::MapRange::dThetaPlus, and LogTrace.

{
   DetIdAssociator::MapRange mapRange;
   mapRange.dThetaPlus = dR;
   mapRange.dThetaMinus = dR;
   mapRange.dPhiPlus = dR;
   mapRange.dPhiMinus = dR;
   if ( delta.first > 0 ) 
     mapRange.dThetaPlus += delta.first;
   else
     mapRange.dThetaMinus += fabs(delta.first);
   if ( delta.second > 0 ) 
     mapRange.dPhiPlus += delta.second;
   else
     mapRange.dPhiMinus += fabs(delta.second);
   LogTrace("TrackAssociator") << "Selection range: (dThetaPlus, dThetaMinus, dPhiPlus, dPhiMinus, dRPreselection): " <<
     mapRange.dThetaPlus << ", " << mapRange.dThetaMinus << ", " << 
     mapRange.dPhiPlus << ", " << mapRange.dPhiMinus << ", " << dR;
   return mapRange;
}
math::XYZPoint TrackDetectorAssociator::getPoint ( const GlobalPoint point) [inline, private]
math::XYZPoint TrackDetectorAssociator::getPoint ( const LocalPoint point) [inline, private]
void TrackDetectorAssociator::getTAMuonChamberMatches ( std::vector< TAMuonChamberMatch > &  matches,
const AssociatorParameters parameters 
) [private]

Definition at line 655 of file TrackDetectorAssociator.cc.

References BoundSurface::bounds(), TrackAssociatorParameters::dRMuonPreselection, TrajectoryStateOnSurface::freeState(), CachedTrajectory::FullTrajectory, TAMuonChamberMatch::id, SteppingHelixStateInfo::isValid(), TrajectoryStateOnSurface::isValid(), Bounds::length(), CSCWireTopology::lengthOfPlane(), TAMuonChamberMatch::localDistanceX, TAMuonChamberMatch::localDistanceY, TrajectoryStateOnSurface::localError(), LogTrace, match(), SteppingHelixStateInfo::momentum(), TrackAssociatorParameters::muonMaxDistanceSigmaX, TrackAssociatorParameters::muonMaxDistanceSigmaY, TrackAssociatorParameters::muonMaxDistanceX, TrackAssociatorParameters::muonMaxDistanceY, CSCWireTopology::narrowWidthOfPlane(), CSCChamberSpecs::oddLayerGeometry(), SteppingHelixStateInfo::position(), FreeTrajectoryState::position(), LocalTrajectoryError::positionError(), mathSSE::sqrt(), GeomDet::surface(), GloballyPositioned< T >::toLocal(), TAMuonChamberMatch::tState, Vector3DBase< T, FrameTag >::unit(), CSCWireTopology::wideWidthOfPlane(), Bounds::width(), CSCWireTopology::wireAngle(), CSCLayerGeometry::wireTopology(), PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), PV3DBase< T, PVType, FrameType >::y(), CSCWireTopology::yOfWire(), and LocalError::yy().

{
   // Strategy:
   //    Propagate through the whole detector, estimate change in eta and phi 
   //    along the trajectory, add this to dRMuon and find DetIds around this 
   //    direction using the map. Then propagate fast to each surface and apply 
   //    final matching criteria.

   // TimerStack timers(TimerStack::Disableable);
   // timers.push("MuonDetIdAssociator::getTrajectoryInMuonDetector");
   // timers.push("MuonDetIdAssociator::getTrajectoryInMuonDetector::propagation",TimerStack::FastMonitoring);
   // timers.pop();
   // get the direction first
   SteppingHelixStateInfo trajectoryPoint = cachedTrajectory_.getStateAtHcal();
   // If trajectory point at HCAL is not valid, try to use the outer most state of the
   // trajectory instead.
   if(! trajectoryPoint.isValid() ) trajectoryPoint = cachedTrajectory_.getOuterState();
   if(! trajectoryPoint.isValid() ) {
      LogTrace("TrackAssociator") << 
        "trajectory position at HCAL is not valid. Assume the track cannot reach muon detectors and skip it";
      return;
   }

   GlobalVector direction = trajectoryPoint.momentum().unit();
   LogTrace("TrackAssociator") << "muon direction: " << direction << "\n\t and corresponding point: " <<
     trajectoryPoint.position() <<"\n";
   
   DetIdAssociator::MapRange mapRange = getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::FullTrajectory),
                                                    parameters.dRMuonPreselection);
     
   // and find chamber DetIds

   // timers.push("MuonDetIdAssociator::getTrajectoryInMuonDetector::getDetIdsCloseToAPoint",TimerStack::FastMonitoring);
   std::set<DetId> muonIdsInRegion = 
     muonDetIdAssociator_->getDetIdsCloseToAPoint(trajectoryPoint.position(), mapRange);
   // timers.pop_and_push("MuonDetIdAssociator::getTrajectoryInMuonDetector::matching",TimerStack::FastMonitoring);
   LogTrace("TrackAssociator") << "Number of chambers to check: " << muonIdsInRegion.size();
   for(std::set<DetId>::const_iterator detId = muonIdsInRegion.begin(); detId != muonIdsInRegion.end(); detId++)
   {
      const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet(*detId);
      // timers.push("MuonDetIdAssociator::getTrajectoryInMuonDetector::matching::localPropagation",TimerStack::FastMonitoring);
      TrajectoryStateOnSurface stateOnSurface = cachedTrajectory_.propagate( &geomDet->surface() );
      if (! stateOnSurface.isValid()) {
         LogTrace("TrackAssociator") << "Failed to propagate the track; moving on\n\t"<<
            detId->rawId() << " not crossed\n";
         continue;
      }
      // timers.pop_and_push("MuonDetIdAssociator::getTrajectoryInMuonDetector::matching::geometryAccess",TimerStack::FastMonitoring);
      LocalPoint localPoint = geomDet->surface().toLocal(stateOnSurface.freeState()->position());
      LocalError localError = stateOnSurface.localError().positionError();
      float distanceX = 0;
      float distanceY = 0;
      float sigmaX = 0.0;
      float sigmaY = 0.0;
      if(const CSCChamber* cscChamber = dynamic_cast<const CSCChamber*>(geomDet) ) {
         const CSCChamberSpecs* chamberSpecs = cscChamber->specs();
         if(! chamberSpecs) {
            LogTrace("TrackAssociator") << "Failed to get CSCChamberSpecs from CSCChamber; moving on\n";
            continue;
         }
         const CSCLayerGeometry* layerGeometry = chamberSpecs->oddLayerGeometry(1);
         if(! layerGeometry) {
            LogTrace("TrackAssociator") << "Failed to get CSCLayerGeometry from CSCChamberSpecs; moving on\n";
            continue;
         }
         const CSCWireTopology* wireTopology = layerGeometry->wireTopology();
         if(! wireTopology) {
            LogTrace("TrackAssociator") << "Failed to get CSCWireTopology from CSCLayerGeometry; moving on\n";
            continue;
         }

         float wideWidth      = wireTopology->wideWidthOfPlane();
         float narrowWidth    = wireTopology->narrowWidthOfPlane();
         float length         = wireTopology->lengthOfPlane();
         // If slanted, there is no y offset between local origin and symmetry center of wire plane
         float yOfFirstWire   = fabs(wireTopology->wireAngle())>1.E-06 ? -0.5*length : wireTopology->yOfWire(1);
         // y offset between local origin and symmetry center of wire plane
         float yCOWPOffset    = yOfFirstWire+0.5*length;

         // tangent of the incline angle from inside the trapezoid
         float tangent = (wideWidth-narrowWidth)/(2.*length);
         // y position wrt bottom of trapezoid
         float yPrime  = localPoint.y()+fabs(yOfFirstWire);
         // half trapezoid width at y' is 0.5 * narrowWidth + x side of triangle with the above tangent and side y'
         float halfWidthAtYPrime = 0.5*narrowWidth+yPrime*tangent;
         distanceX = fabs(localPoint.x()) - halfWidthAtYPrime;
         distanceY = fabs(localPoint.y()-yCOWPOffset) - 0.5*length;
         sigmaX = distanceX/sqrt(localError.xx());
         sigmaY = distanceY/sqrt(localError.yy());
      } else {
         distanceX = fabs(localPoint.x()) - geomDet->surface().bounds().width()/2.;
         distanceY = fabs(localPoint.y()) - geomDet->surface().bounds().length()/2.;
         sigmaX = distanceX/sqrt(localError.xx());
         sigmaY = distanceY/sqrt(localError.yy());
      }
      // timers.pop_and_push("MuonDetIdAssociator::getTrajectoryInMuonDetector::matching::checking",TimerStack::FastMonitoring);
      if ( (distanceX < parameters.muonMaxDistanceX && distanceY < parameters.muonMaxDistanceY) ||
           (sigmaX < parameters.muonMaxDistanceSigmaX && sigmaY < parameters.muonMaxDistanceSigmaY) ) {
         LogTrace("TrackAssociator") << "found a match, DetId: " << detId->rawId();
         TAMuonChamberMatch match;
         match.tState = stateOnSurface;
         match.localDistanceX = distanceX;
         match.localDistanceY = distanceY;
         match.id = *detId;
         matches.push_back(match);
      }
      //timers.pop();
   }
   //timers.pop();
   
}
math::XYZVector TrackDetectorAssociator::getVector ( const GlobalVector vec) [inline, private]
math::XYZVector TrackDetectorAssociator::getVector ( const LocalVector vec) [inline, private]
void TrackDetectorAssociator::init ( const edm::EventSetup iSetup) [private]

Definition at line 124 of file TrackDetectorAssociator.cc.

References anyDirection, SteppingHelixPropagator::applyRadX0Correction(), ecalTB2006H4_GenSimDigiReco_cfg::bField, edm::EventSetup::get(), SteppingHelixPropagator::setMaterialMode(), and SteppingHelixPropagator_cfi::SteppingHelixPropagator.

{
   // access the calorimeter geometry
   iSetup.get<CaloGeometryRecord>().get(theCaloGeometry_);
   if (!theCaloGeometry_.isValid()) 
     throw cms::Exception("FatalError") << "Unable to find CaloGeometryRecord in event!\n";
   
   // get the tracking Geometry
   iSetup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry_);
   if (!theTrackingGeometry_.isValid()) 
     throw cms::Exception("FatalError") << "Unable to find GlobalTrackingGeometryRecord in event!\n";
   
   if (useDefaultPropagator_ && (! defProp_ || theMagneticFeildWatcher_.check(iSetup) ) ) {
      // setup propagator
      edm::ESHandle<MagneticField> bField;
      iSetup.get<IdealMagneticFieldRecord>().get(bField);
      
      SteppingHelixPropagator* prop  = new SteppingHelixPropagator(&*bField,anyDirection);
      prop->setMaterialMode(false);
      prop->applyRadX0Correction(true);
      // prop->setDebug(true); // tmp
      defProp_ = prop;
      setPropagator(defProp_);
   }

   iSetup.get<DetIdAssociatorRecord>().get("EcalDetIdAssociator", ecalDetIdAssociator_);
   iSetup.get<DetIdAssociatorRecord>().get("HcalDetIdAssociator", hcalDetIdAssociator_);
   iSetup.get<DetIdAssociatorRecord>().get("HODetIdAssociator", hoDetIdAssociator_);
   iSetup.get<DetIdAssociatorRecord>().get("CaloDetIdAssociator", caloDetIdAssociator_);
   iSetup.get<DetIdAssociatorRecord>().get("MuonDetIdAssociator", muonDetIdAssociator_);
   iSetup.get<DetIdAssociatorRecord>().get("PreshowerDetIdAssociator", preshowerDetIdAssociator_);
}
void TrackDetectorAssociator::setPropagator ( const Propagator ptr)
void TrackDetectorAssociator::useDefaultPropagator ( )

Member Data Documentation

Definition at line 187 of file TrackDetectorAssociator.h.

Definition at line 193 of file TrackDetectorAssociator.h.

Definition at line 186 of file TrackDetectorAssociator.h.

Definition at line 190 of file TrackDetectorAssociator.h.

Definition at line 191 of file TrackDetectorAssociator.h.

Definition at line 192 of file TrackDetectorAssociator.h.

Definition at line 185 of file TrackDetectorAssociator.h.

Definition at line 194 of file TrackDetectorAssociator.h.

Definition at line 195 of file TrackDetectorAssociator.h.

Definition at line 197 of file TrackDetectorAssociator.h.

Definition at line 200 of file TrackDetectorAssociator.h.

Definition at line 198 of file TrackDetectorAssociator.h.

Definition at line 188 of file TrackDetectorAssociator.h.