CMS 3D CMS Logo

Public Member Functions | Public Attributes | Private Attributes

CSCHaloAlgo Class Reference

#include <CSCHaloAlgo.h>

List of all members.

Public Member Functions

reco::CSCHaloData Calculate (const CSCGeometry &TheCSCGeometry, edm::Handle< reco::TrackCollection > &TheCSCTracks, edm::Handle< reco::MuonCollection > &TheMuons, edm::Handle< CSCSegmentCollection > &TheCSCSegments, edm::Handle< CSCRecHit2DCollection > &TheCSCRecHits, edm::Handle< L1MuGMTReadoutCollection > &TheL1GMTReadout, edm::Handle< edm::TriggerResults > &TheHLTResults, const edm::TriggerNames *triggerNames, const edm::Handle< CSCALCTDigiCollection > &TheALCTs)
 CSCHaloAlgo ()
void SetDetaThreshold (float x)
void SetDphiThreshold (float x)
void SetExpectedBX (int x)
void SetMatchingDEtaThreshold (float x)
void SetMatchingDPhiThreshold (float x)
void SetMatchingDWireThreshold (int x)
void SetMinMaxInnerRadius (float min, float max)
void SetMinMaxOuterMomentumTheta (float min, float max)
void SetMinMaxOuterRadius (float min, float max)
void SetNormChi2Threshold (float x)
void SetRecHitTime0 (float x)
void SetRecHitTimeWindow (float x)
 ~CSCHaloAlgo ()

Public Attributes

std::vector< edm::InputTagvIT_HLTBit

Private Attributes

float deta_threshold
float dphi_threshold
int expected_BX
float matching_deta_threshold
float matching_dphi_threshold
int matching_dwire_threshold
float max_inner_radius
float max_outer_radius
float max_outer_theta
float min_inner_radius
float min_outer_radius
float min_outer_theta
float norm_chi2_threshold
float recHit_t0
float recHit_twindow

Detailed Description

Definition at line 80 of file CSCHaloAlgo.h.


Constructor & Destructor Documentation

CSCHaloAlgo::CSCHaloAlgo ( )
CSCHaloAlgo::~CSCHaloAlgo ( ) [inline]

Definition at line 84 of file CSCHaloAlgo.h.

{}

Member Function Documentation

reco::CSCHaloData CSCHaloAlgo::Calculate ( const CSCGeometry TheCSCGeometry,
edm::Handle< reco::TrackCollection > &  TheCSCTracks,
edm::Handle< reco::MuonCollection > &  TheMuons,
edm::Handle< CSCSegmentCollection > &  TheCSCSegments,
edm::Handle< CSCRecHit2DCollection > &  TheCSCRecHits,
edm::Handle< L1MuGMTReadoutCollection > &  TheL1GMTReadout,
edm::Handle< edm::TriggerResults > &  TheHLTResults,
const edm::TriggerNames triggerNames,
const edm::Handle< CSCALCTDigiCollection > &  TheALCTs 
)

Definition at line 35 of file CSCHaloAlgo.cc.

References CSCGeometry::chamber(), chambers, CSC(), CSCDetId, CSCDetId::endcap(), PV3DBase< T, PVType, FrameType >::eta(), reco::CSCHaloData::GetCSCTrackImpactPositions(), L1MuGMTReadoutCollection::getRecords(), reco::CSCHaloData::GetTracks(), CSCGeometry::idToDetUnit(), getHLTprescales::index, edm::HandleBase::isValid(), j, label, DetId::Muon, PV3DBase< T, PVType, FrameType >::phi(), Pi, edm::Handle< T >::product(), edm::RefVector< C, T, F >::push_back(), reco::CSCHaloData::SetHLTBit(), reco::CSCHaloData::SetNOutOfTimeHits(), reco::CSCHaloData::SetNOutOfTimeTriggers(), reco::CSCHaloData::SetNumberOfHaloTriggers(), findQualityFiles::size, GeomDet::surface(), theta(), Surface::toGlobal(), GeomDet::toGlobal(), edm::TriggerNames::triggerIndex(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::z(), and z.

{
  reco::CSCHaloData TheCSCHaloData;
  if( TheCSCTracks.isValid() )
    {
      for( reco::TrackCollection::const_iterator iTrack = TheCSCTracks->begin() ; iTrack != TheCSCTracks->end() ; iTrack++ )
        {
          bool StoreTrack = false;
          // Calculate global phi coordinate for central most rechit in the track
          float innermost_global_z = 1500.;
          float outermost_global_z = 0.;
          GlobalPoint InnerMostGlobalPosition(0.,0.,0.);  // smallest abs(z)
          GlobalPoint OuterMostGlobalPosition(0.,0.,0.);  // largest abs(z)
          int nCSCHits = 0;
          for(unsigned int j = 0 ; j < iTrack->extra()->recHits().size(); j++ )
            {
              edm::Ref<TrackingRecHitCollection> hit( iTrack->extra()->recHits(), j );
              if( !hit->isValid() ) continue;
              DetId TheDetUnitId(hit->geographicalId());
              if( TheDetUnitId.det() != DetId::Muon ) continue;
              if( TheDetUnitId.subdetId() != MuonSubdetId::CSC ) continue;

              const GeomDetUnit *TheUnit = TheCSCGeometry.idToDetUnit(TheDetUnitId);
              LocalPoint TheLocalPosition = hit->localPosition();  
              const BoundPlane& TheSurface = TheUnit->surface();
              const GlobalPoint TheGlobalPosition = TheSurface.toGlobal(TheLocalPosition);

              float z = TheGlobalPosition.z();
              if( TMath::Abs(z) < innermost_global_z )
                {
                  innermost_global_z = TMath::Abs(z);
                  InnerMostGlobalPosition = GlobalPoint( TheGlobalPosition);
                }
              if( TMath::Abs(z) > outermost_global_z )
                {
                  outermost_global_z = TMath::Abs(z);
                  OuterMostGlobalPosition = GlobalPoint( TheGlobalPosition );
                }
              nCSCHits ++;
            }

          if( nCSCHits < 3 ) continue; // This needs to be optimized, but is the minimum 
 
          if( OuterMostGlobalPosition.x() == 0. || OuterMostGlobalPosition.y() == 0. || OuterMostGlobalPosition.z() == 0. ) 
            continue;
          if( InnerMostGlobalPosition.x() == 0. || InnerMostGlobalPosition.y() == 0. || InnerMostGlobalPosition.z() == 0. )
            continue;
          
          //Its a CSC Track,store it if it passes halo selection 
          StoreTrack = true;      

          float deta = TMath::Abs( OuterMostGlobalPosition.eta() - InnerMostGlobalPosition.eta() );
          float dphi = TMath::ACos( TMath::Cos( OuterMostGlobalPosition.phi() - InnerMostGlobalPosition.phi() ) ) ;
          float theta = iTrack->outerMomentum().theta();
          float innermost_x = InnerMostGlobalPosition.x() ;
          float innermost_y = InnerMostGlobalPosition.y();
          float outermost_x = OuterMostGlobalPosition.x();
          float outermost_y = OuterMostGlobalPosition.y();
          float innermost_r = TMath::Sqrt(innermost_x *innermost_x + innermost_y * innermost_y );
          float outermost_r = TMath::Sqrt(outermost_x *outermost_x + outermost_y * outermost_y );
          
          if( deta < deta_threshold )
            StoreTrack = false;
          if( theta > min_outer_theta && theta < max_outer_theta )
            StoreTrack = false;
          if( dphi > dphi_threshold )
            StoreTrack = false;
          if( innermost_r < min_inner_radius )
            StoreTrack = false;
          if( innermost_r > max_inner_radius )
            StoreTrack = false;
          if( outermost_r < min_outer_radius )
            StoreTrack = false;
          if( outermost_r > max_outer_radius )
            StoreTrack  = false;
          if( iTrack->normalizedChi2() > norm_chi2_threshold )
            StoreTrack = false;



          if( StoreTrack )
            {
              TheCSCHaloData.GetCSCTrackImpactPositions().push_back( InnerMostGlobalPosition );

              edm::Ref<TrackCollection> TheTrackRef( TheCSCTracks, iTrack - TheCSCTracks->begin() ) ;
              TheCSCHaloData.GetTracks().push_back( TheTrackRef );
            }
        }
    }

   if( TheHLTResults.isValid() )
     {
       bool EventPasses = false;
       for( unsigned int index = 0 ; index < vIT_HLTBit.size(); index++)
         {
           if( vIT_HLTBit[index].label().size() )
             {
               //Get the HLT bit and check to make sure it is valid 
               unsigned int bit = triggerNames->triggerIndex( vIT_HLTBit[index].label().c_str());
               if( bit < TheHLTResults->size() )
                 {
                   //If any of the HLT names given by the user accept, then the event passes
                   if( TheHLTResults->accept( bit ) && !TheHLTResults->error( bit ) )
                     {
                       EventPasses = true;
                     }
                 }
             }
         }
       if( EventPasses )
         TheCSCHaloData.SetHLTBit(true);
       else
         TheCSCHaloData.SetHLTBit(false);
     }

   if( TheL1GMTReadout.isValid() )
     {
       L1MuGMTReadoutCollection const *gmtrc = TheL1GMTReadout.product ();
       std::vector < L1MuGMTReadoutRecord > gmt_records = gmtrc->getRecords ();
       std::vector < L1MuGMTReadoutRecord >::const_iterator igmtrr;
       
       int icsc = 0;
       int PlusZ = 0 ;
       int MinusZ = 0 ;
       // Check to see if CSC BeamHalo trigger is tripped
       for (igmtrr = gmt_records.begin (); igmtrr != gmt_records.end (); igmtrr++)
         {
           std::vector < L1MuRegionalCand >::const_iterator iter1;
           std::vector < L1MuRegionalCand > rmc;
           rmc = igmtrr->getCSCCands ();
           for (iter1 = rmc.begin (); iter1 != rmc.end (); iter1++)
             {
              if (!(*iter1).empty ())
                {
                  if ((*iter1).isFineHalo ())
                    {
                      float halophi = iter1->phiValue();
                      halophi = halophi > TMath::Pi() ? halophi - 2.*TMath::Pi() : halophi;
                      float haloeta = iter1->etaValue();
                      bool HaloIsGood = true;
                      // Check if halo trigger is faked by any collision muons
                      if( TheMuons.isValid() )
                        {
                          float dphi = 9999.;
                          float deta = 9999.;
                          for( reco::MuonCollection::const_iterator mu = TheMuons->begin(); mu != TheMuons->end() && HaloIsGood ; mu++ )
                            {
                              // Don't match with SA-only muons
                              if( mu->isStandAloneMuon() && !mu->isTrackerMuon() && !mu->isGlobalMuon() )  continue;
                              
                              /*
                              if(!mu->isTrackerMuon())
                                {
                                  if( mu->isStandAloneMuon() )
                                    {
                                      //make sure that this SA muon is not actually a halo-like muon
                                      float theta =  mu->outerTrack()->outerMomentum().theta();
                                      float deta = TMath::Abs(mu->outerTrack()->outerPosition().eta() - mu->outerTrack()->innerPosition().eta());
                                      if( theta < min_outer_theta || theta > max_outer_theta )  //halo-like
                                        continue;
                                      else if ( deta > deta_threshold ) //halo-like
                                        continue;
                                    }
                                }
                              */
                            
                              const std::vector<MuonChamberMatch> chambers = mu->matches();
                              for(std::vector<MuonChamberMatch>::const_iterator iChamber = chambers.begin();
                                  iChamber != chambers.end() ; iChamber ++ )
                                {
                                  if( iChamber->detector() != MuonSubdetId::CSC ) continue;
                                  for( std::vector<reco::MuonSegmentMatch>::const_iterator iSegment = iChamber->segmentMatches.begin() ; 
                                       iSegment != iChamber->segmentMatches.end(); ++iSegment )
                                    {
                                      edm::Ref<CSCSegmentCollection> cscSegment = iSegment->cscSegmentRef;
                                      std::vector<CSCRecHit2D> hits = cscSegment -> specificRecHits();
                                      for( std::vector<CSCRecHit2D>::iterator iHit = hits.begin();
                                           iHit != hits.end() ; iHit++ )
                                        {
                                          DetId TheDetUnitId(iHit->cscDetId());
                                          const GeomDetUnit *TheUnit = TheCSCGeometry.idToDetUnit(TheDetUnitId);
                                          LocalPoint TheLocalPosition = iHit->localPosition();
                                          const BoundPlane& TheSurface = TheUnit->surface();
                                          GlobalPoint TheGlobalPosition = TheSurface.toGlobal(TheLocalPosition);
                                          
                                          float phi_ = TheGlobalPosition.phi();
                                          float eta_ = TheGlobalPosition.eta();
                                          deta = deta < TMath::Abs( eta_ - haloeta ) ? deta : TMath::Abs( eta_ - haloeta );
                                          dphi = dphi < TMath::Abs( phi_ - halophi ) ? dphi : TMath::Abs( phi_ - halophi );
                                        }
                                    }
                                }
                              if ( dphi < matching_dphi_threshold && deta < matching_deta_threshold) 
                                HaloIsGood = false; // i.e., collision muon likely faked halo trigger
                            }
                        }
                      if( !HaloIsGood ) 
                        continue;
                      if( (*iter1).etaValue() > 0 )
                        PlusZ++;
                      else
                        MinusZ++;
                    }
                  else
                    icsc++;
                }
             }
         }
       TheCSCHaloData.SetNumberOfHaloTriggers(PlusZ, MinusZ);
     }

   // Loop over CSCALCTDigi collection to look for out-of-time chamber triggers 
   // A collision muon in real data should only have ALCTDigi::getBX() = 3 ( in MC, it will be 6 )
   // Note that there could be two ALCTs per chamber 
   short int n_alctsP=0;
   short int n_alctsM=0;
   if(TheALCTs.isValid())
     {
       for (CSCALCTDigiCollection::DigiRangeIterator j=TheALCTs->begin(); j!=TheALCTs->end(); j++) 
         {
           const CSCALCTDigiCollection::Range& range =(*j).second;
           CSCDetId detId((*j).first.rawId());
           for (CSCALCTDigiCollection::const_iterator digiIt = range.first; digiIt!=range.second; ++digiIt)
             {
               if( (*digiIt).isValid() && ( (*digiIt).getBX() < expected_BX ) )
                 {
                   int digi_endcap  = detId.endcap();
                   int digi_station = detId.station();
                   int digi_ring    = detId.ring();
                   int digi_chamber = detId.chamber();
                   int digi_wire    = digiIt->getKeyWG();
                   if( digi_station == 1 && digi_ring == 4 )   //hack
                     digi_ring = 1;
                   
                   bool DigiIsGood = true;
                   int dwire = 999.;
                   if( TheMuons.isValid() ) 
                     {
                       //Check if there are any collision muons with hits in the vicinity of the digi
                       for(reco::MuonCollection::const_iterator mu = TheMuons->begin(); mu!= TheMuons->end() && DigiIsGood ; mu++ )
                         {
                           if( !mu->isTrackerMuon() && !mu->isGlobalMuon() && mu->isStandAloneMuon() ) continue;
                           /*
                             if(!mu->isTrackerMuon())
                                {
                                  if( mu->isStandAloneMuon() )
                                    {
                                      //make sure that this SA muon is not actually a halo-like muon
                                      float theta =  mu->outerTrack()->outerMomentum().theta();
                                      float deta = TMath::Abs(mu->outerTrack()->outerPosition().eta() - mu->outerTrack()->innerPosition().eta());
                                      if( theta < min_outer_theta || theta > max_outer_theta )  //halo-like
                                        continue;
                                      else if ( deta > deta_threshold ) //halo-like
                                        continue;
                                    }
                                    }
                           */
                          
                           const std::vector<MuonChamberMatch> chambers = mu->matches();
                           for(std::vector<MuonChamberMatch>::const_iterator iChamber = chambers.begin();
                               iChamber != chambers.end(); iChamber ++ )
                             {
                               if( iChamber->detector() != MuonSubdetId::CSC ) continue;
                               for( std::vector<reco::MuonSegmentMatch>::const_iterator iSegment = iChamber->segmentMatches.begin();
                                    iSegment != iChamber->segmentMatches.end(); iSegment++ )
                                 {
                                   edm::Ref<CSCSegmentCollection> cscSegRef = iSegment->cscSegmentRef;
                                   std::vector<CSCRecHit2D> hits = cscSegRef->specificRecHits();
                                   for( std::vector<CSCRecHit2D>::iterator iHit = hits.begin();
                                        iHit != hits.end(); iHit++ )
                                     {
                                       if( iHit->cscDetId().endcap() != digi_endcap ) continue;
                                       if( iHit->cscDetId().station() != digi_station ) continue;
                                       if( iHit->cscDetId().ring() != digi_ring ) continue;
                                       if( iHit->cscDetId().chamber() != digi_chamber ) continue;
                                       CSCRecHit2D::ChannelContainer hitwires = iHit->wgroups();
                                       int nwires = hitwires.size();
                                       int center_id = nwires/2 + 1;
                                       int hit_wire = hitwires[center_id -1 ];
                                       dwire = dwire < TMath::Abs(hit_wire - digi_wire)? dwire : TMath::Abs(hit_wire - digi_wire );
                                     }
                                 }
                             }
                           if( dwire <= matching_dwire_threshold ) 
                             DigiIsGood = false;  // collision-like muon is close to this digi
                         }
                     }
                   // only count out of time digis if they are not matched to collision muons
                   if( DigiIsGood ) 
                     {
                       if( detId.endcap() == 1 ) 
                         n_alctsP++;
                       else if ( detId.endcap() ==  2) 
                         n_alctsM++;
                     }
                 }
             }
         }
     }
   TheCSCHaloData.SetNOutOfTimeTriggers(n_alctsP,n_alctsM);

   // Loop over the CSCRecHit2D collection to look for out-of-time recHits
   // Out-of-time is defined as tpeak outside [t_0 + TOF - t_window, t_0 + TOF + t_window]
   // where t_0 and t_window are configurable parameters
   short int n_recHitsP = 0;
   short int n_recHitsM = 0;
   if( TheCSCRecHits.isValid() )
     {
       CSCRecHit2DCollection::const_iterator dRHIter;
       for (dRHIter = TheCSCRecHits->begin(); dRHIter != TheCSCRecHits->end(); dRHIter++) 
         {
           if ( !((*dRHIter).isValid()) ) continue;  // only interested in valid hits
           CSCDetId idrec = (CSCDetId)(*dRHIter).cscDetId();
           float RHTime = (*dRHIter).tpeak();
           LocalPoint rhitlocal = (*dRHIter).localPosition();
           const CSCChamber* chamber = TheCSCGeometry.chamber(idrec);
           GlobalPoint globalPosition = chamber->toGlobal(rhitlocal);
           float globZ = globalPosition.z();
           if ( RHTime < (recHit_t0 - recHit_twindow) )
             {
               if( globZ > 0 )
                 n_recHitsP++;
               else
                 n_recHitsM++;
             }
           
           /*

           float globX = globalPosition.x();
           float globY = globalPosition.y();
           float globZ = globalPosition.z();
           float TOF = (sqrt(globX*globX+ globY*globY + globZ*globZ))/29.9792458 ; //cm -> ns
           if ( (RHTime < (recHit_t0 + TOF - recHit_twindow)) || (RHTime > (recHit_t0 + TOF + recHit_twindow)) )
             {
               if( globZ > 0 ) 
                 n_recHitsP++;
               else
                 n_recHitsM++;
             }
           */
         }
     }
   TheCSCHaloData.SetNOutOfTimeHits(n_recHitsP+n_recHitsM);

   return TheCSCHaloData;
}
void CSCHaloAlgo::SetDetaThreshold ( float  x) [inline]

Definition at line 93 of file CSCHaloAlgo.h.

References deta_threshold, and ExpressReco_HICollisions_FallBack::x.

void CSCHaloAlgo::SetDphiThreshold ( float  x) [inline]

Definition at line 96 of file CSCHaloAlgo.h.

References dphi_threshold, and ExpressReco_HICollisions_FallBack::x.

void CSCHaloAlgo::SetExpectedBX ( int  x) [inline]

Definition at line 100 of file CSCHaloAlgo.h.

References expected_BX, and ExpressReco_HICollisions_FallBack::x.

{ expected_BX = x ;}
void CSCHaloAlgo::SetMatchingDEtaThreshold ( float  x) [inline]
void CSCHaloAlgo::SetMatchingDPhiThreshold ( float  x) [inline]
void CSCHaloAlgo::SetMatchingDWireThreshold ( int  x) [inline]
void CSCHaloAlgo::SetMinMaxInnerRadius ( float  min,
float  max 
) [inline]

Definition at line 94 of file CSCHaloAlgo.h.

References max(), max_inner_radius, min, and min_inner_radius.

void CSCHaloAlgo::SetMinMaxOuterMomentumTheta ( float  min,
float  max 
) [inline]

Definition at line 101 of file CSCHaloAlgo.h.

References max(), max_outer_theta, min, and min_outer_theta.

void CSCHaloAlgo::SetMinMaxOuterRadius ( float  min,
float  max 
) [inline]

Definition at line 95 of file CSCHaloAlgo.h.

References max(), max_outer_radius, min, and min_outer_radius.

void CSCHaloAlgo::SetNormChi2Threshold ( float  x) [inline]
void CSCHaloAlgo::SetRecHitTime0 ( float  x) [inline]

Definition at line 98 of file CSCHaloAlgo.h.

References recHit_t0, and ExpressReco_HICollisions_FallBack::x.

{ recHit_t0 = x;}
void CSCHaloAlgo::SetRecHitTimeWindow ( float  x) [inline]

Definition at line 99 of file CSCHaloAlgo.h.

References recHit_twindow, and ExpressReco_HICollisions_FallBack::x.


Member Data Documentation

float CSCHaloAlgo::deta_threshold [private]

Definition at line 107 of file CSCHaloAlgo.h.

Referenced by SetDetaThreshold().

float CSCHaloAlgo::dphi_threshold [private]

Definition at line 114 of file CSCHaloAlgo.h.

Referenced by SetDphiThreshold().

int CSCHaloAlgo::expected_BX [private]

Definition at line 118 of file CSCHaloAlgo.h.

Referenced by SetExpectedBX().

Definition at line 120 of file CSCHaloAlgo.h.

Referenced by SetMatchingDEtaThreshold().

Definition at line 119 of file CSCHaloAlgo.h.

Referenced by SetMatchingDPhiThreshold().

Definition at line 121 of file CSCHaloAlgo.h.

Referenced by SetMatchingDWireThreshold().

Definition at line 111 of file CSCHaloAlgo.h.

Referenced by SetMinMaxInnerRadius().

Definition at line 113 of file CSCHaloAlgo.h.

Referenced by SetMinMaxOuterRadius().

Definition at line 108 of file CSCHaloAlgo.h.

Referenced by SetMinMaxOuterMomentumTheta().

Definition at line 110 of file CSCHaloAlgo.h.

Referenced by SetMinMaxInnerRadius().

Definition at line 112 of file CSCHaloAlgo.h.

Referenced by SetMinMaxOuterRadius().

Definition at line 109 of file CSCHaloAlgo.h.

Referenced by SetMinMaxOuterMomentumTheta().

Definition at line 115 of file CSCHaloAlgo.h.

Referenced by SetNormChi2Threshold().

float CSCHaloAlgo::recHit_t0 [private]

Definition at line 116 of file CSCHaloAlgo.h.

Referenced by SetRecHitTime0().

float CSCHaloAlgo::recHit_twindow [private]

Definition at line 117 of file CSCHaloAlgo.h.

Referenced by SetRecHitTimeWindow().

Definition at line 91 of file CSCHaloAlgo.h.