/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoMET/METAlgorithms/src/GlobalHaloAlgo.cc

Go to the documentation of this file.

#include "RecoMET/METAlgorithms/interface/GlobalHaloAlgo.h"
/*
  [class]:  GlobalHaloAlgo
  [authors]: R. Remington, The University of Florida
  [description]: See GlobalHaloAlgo.h
  [date]: October 15, 2009
*/
using namespace std;
using namespace edm;
using namespace reco;

int Phi_To_HcaliPhi(float phi) 
{
  phi = phi < 0 ? phi + 2.*TMath::Pi() : phi ;
  float phi_degrees = phi * (360.) / ( 2. * TMath::Pi() ) ;
  int iPhi = (int) ( ( phi_degrees/5. ) + 1.);
   
  return iPhi < 73 ? iPhi : 73 ;
}

int Phi_To_EcaliPhi(float phi) 
{
  phi = phi < 0 ? phi + 2.*TMath::Pi() : phi ;
  float phi_degrees = phi * (360.) / ( 2. * TMath::Pi() ) ;
  int iPhi = (int) (  phi_degrees  + 1.);
   
  return iPhi < 361 ? iPhi : 360 ;
}

GlobalHaloAlgo::GlobalHaloAlgo()
{
  // Defaults are "loose"
  Ecal_R_Min = 110.;   // Tight: 200.
  Ecal_R_Max = 330.;   // Tight: 250. 
  Hcal_R_Min = 110.;   // Tight: 220.
  Hcal_R_Max = 490.;   // Tight: 350.
  
}

reco::GlobalHaloData GlobalHaloAlgo::Calculate(const CaloGeometry& TheCaloGeometry, const CSCGeometry& TheCSCGeometry, const reco::CaloMET& TheCaloMET, edm::Handle< edm::View<Candidate> >& TheCaloTowers, edm::Handle<CSCSegmentCollection>& TheCSCSegments, edm::Handle<CSCRecHit2DCollection>& TheCSCRecHits,  const CSCHaloData& TheCSCHaloData, const EcalHaloData& TheEcalHaloData, const HcalHaloData& TheHcalHaloData)
{
  
  GlobalHaloData TheGlobalHaloData;
  float METOverSumEt = TheCaloMET.sumEt() ? TheCaloMET.pt() / TheCaloMET.sumEt() : 0 ;
  TheGlobalHaloData.SetMETOverSumEt(METOverSumEt);

  //int EcalOverlapping_CSCRecHits[73];
  //int EcalOverlapping_CSCSegments[73];

  int EcalOverlapping_CSCRecHits[361];
  int EcalOverlapping_CSCSegments[361];
  int HcalOverlapping_CSCRecHits[73];
  int HcalOverlapping_CSCSegments[73];

  if( TheCSCSegments.isValid() )
    {
      for(CSCSegmentCollection::const_iterator iSegment = TheCSCSegments->begin(); iSegment != TheCSCSegments->end(); iSegment++) 
        {
          bool EcalOverlap[361];
          bool HcalOverlap[73];
          for( int i = 0 ; i < 361 ; i++ ) 
            {
              EcalOverlap[i] = false;
              if( i < 73 ) HcalOverlap[i] = false;
            }
          
          std::vector<CSCRecHit2D> Hits = iSegment->specificRecHits() ;
          for(std::vector<CSCRecHit2D>::iterator iHit = Hits.begin() ; iHit != Hits.end(); iHit++ )
            {
              DetId TheDetUnitId(iHit->geographicalId());
              if( TheDetUnitId.det() != DetId::Muon ) continue;
              if( TheDetUnitId.subdetId() != MuonSubdetId::CSC ) continue;

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

              int Hcal_iphi = Phi_To_HcaliPhi( TheGlobalPosition.phi() ) ;
              int Ecal_iphi = Phi_To_EcaliPhi( TheGlobalPosition.phi() ) ;
              float x = TheGlobalPosition.x(); 
              float y = TheGlobalPosition.y();
              
              float r = TMath::Sqrt( x*x + y*y);
              
              if( r < Ecal_R_Max && r > Ecal_R_Min )
                EcalOverlap[Ecal_iphi] = true;
              if( r < Hcal_R_Max && r > Hcal_R_Max ) 
                HcalOverlap[Hcal_iphi] = true;
            }
          for( int i = 0 ; i < 361 ; i++ ) 
            {
              if( EcalOverlap[i] )  EcalOverlapping_CSCSegments[i]++;
              if( i < 73 && HcalOverlap[i] )
                HcalOverlapping_CSCSegments[i]++;
            }
        } 
    }
  if( TheCSCRecHits.isValid() )
    {
      for(CSCRecHit2DCollection::const_iterator iCSCRecHit = TheCSCRecHits->begin();   iCSCRecHit != TheCSCRecHits->end(); iCSCRecHit++ )
        {
          
          DetId TheDetUnitId(iCSCRecHit->geographicalId());
          if( TheDetUnitId.det() != DetId::Muon ) continue;
          if( TheDetUnitId.subdetId() != MuonSubdetId::CSC ) continue;
          
          const GeomDetUnit *TheUnit = TheCSCGeometry.idToDetUnit(TheDetUnitId);
          LocalPoint TheLocalPosition = iCSCRecHit->localPosition();  
          const BoundPlane& TheSurface = TheUnit->surface();
          const GlobalPoint TheGlobalPosition = TheSurface.toGlobal(TheLocalPosition);
          
          int Hcaliphi = Phi_To_HcaliPhi( TheGlobalPosition.phi() ) ;
          int Ecaliphi = Phi_To_EcaliPhi( TheGlobalPosition.phi() ) ;
          float x = TheGlobalPosition.x(); 
          float y = TheGlobalPosition.y();
          
          float r = TMath::Sqrt(x*x + y*y);
          
          if( r < Ecal_R_Max && r > Ecal_R_Min )
            EcalOverlapping_CSCRecHits[Ecaliphi] ++;
          if( r < Hcal_R_Max && r > Hcal_R_Max ) 
            HcalOverlapping_CSCRecHits[Hcaliphi] ++ ;
        }
    }  

  // In development....
  // Get Ecal Wedges
  std::vector<PhiWedge> EcalWedges = TheEcalHaloData.GetPhiWedges();
  
  // Get Hcal Wedges
  std::vector<PhiWedge> HcalWedges = TheHcalHaloData.GetPhiWedges();
  
  //Get Ref to CSC Tracks
  //edm::RefVector<reco::TrackCollection> TheCSCTracks = TheCSCHaloData.GetTracks();
  //for(unsigned int i = 0 ; i < TheCSCTracks.size() ; i++ )
  //edm::Ref<reco::TrackCollection> iTrack( TheCSCTracks, i );
  
  // Get global positions of central most rechit of CSC Halo tracks
  std::vector<GlobalPoint> TheGlobalPositions = TheCSCHaloData.GetCSCTrackImpactPositions();

  // Container to store Ecal/Hcal iPhi values matched to impact point of CSC tracks
  std::vector<int> vEcaliPhi, vHcaliPhi;

  // Keep track of number of calo pointing CSC halo tracks that do not match to Phi wedges
  int N_Unmatched_Tracks = 0;  
  
  for( std::vector<GlobalPoint>::iterator Pos = TheGlobalPositions.begin() ; Pos != TheGlobalPositions.end() ; Pos ++ ) 
    {
      // Calculate global phi coordinate for central most rechit in the track
      float global_phi = Pos->phi();
      float global_r = TMath::Sqrt(Pos->x()*Pos->x() + Pos->y()*Pos->y());
      
      // Convert global phi to iPhi
      int global_EcaliPhi = Phi_To_EcaliPhi( global_phi );
      int global_HcaliPhi = Phi_To_HcaliPhi( global_phi );
      bool MATCHED = false;
      
      //Loop over Ecal Phi Wedges 
      for( std::vector<PhiWedge>::iterator iWedge = EcalWedges.begin() ; iWedge != EcalWedges.end() ; iWedge++ )
        {
          if( (TMath::Abs( global_EcaliPhi - iWedge->iPhi() ) <= 5 ) && (global_r >  Ecal_R_Min && global_r < Ecal_R_Max ) )
            {
              bool StoreWedge = true;
              for( unsigned int i = 0 ; i< vEcaliPhi.size() ; i++ ) if ( vEcaliPhi[i] == iWedge->iPhi() ) StoreWedge = false;
              
              if( StoreWedge ) 
                {
                  PhiWedge NewWedge(*iWedge);
                  NewWedge.SetOverlappingCSCSegments( EcalOverlapping_CSCSegments[iWedge->iPhi()] );
                  NewWedge.SetOverlappingCSCRecHits( EcalOverlapping_CSCRecHits[iWedge->iPhi()] );
                  vEcaliPhi.push_back( iWedge->iPhi() );
                  TheGlobalHaloData.GetMatchedEcalPhiWedges().push_back( NewWedge );
                }
              MATCHED = true;
            }
        }
      //Loop over Hcal Phi Wedges 
      for( std::vector<PhiWedge>::iterator iWedge = HcalWedges.begin() ; iWedge != HcalWedges.end() ; iWedge++ )
        {
          if(  (TMath::Abs( global_HcaliPhi - iWedge->iPhi() ) <=  2 ) && (global_r > Hcal_R_Min && global_r < Hcal_R_Max ) )
            {
              bool StoreWedge  = true;
              for( unsigned int i = 0 ; i < vHcaliPhi.size() ; i++ ) if(  vHcaliPhi[i] == iWedge->iPhi() ) StoreWedge = false;
              
              if( StoreWedge ) 
                {
                  vHcaliPhi.push_back( iWedge->iPhi() ) ;
                  PhiWedge NewWedge(*iWedge);
                  NewWedge.SetOverlappingCSCSegments( HcalOverlapping_CSCSegments[iWedge->iPhi()] );
                  NewWedge.SetOverlappingCSCRecHits(  HcalOverlapping_CSCRecHits[iWedge->iPhi()] );               
                  PhiWedge wedge(*iWedge);
                  TheGlobalHaloData.GetMatchedHcalPhiWedges().push_back( NewWedge ) ; 
                }
              MATCHED = true;
            }
        }
      if( !MATCHED ) N_Unmatched_Tracks ++;
    }
  
  // Corrections to MEx, MEy
  float dMEx = 0.; 
  float dMEy = 0.;
  // Loop over calotowers and correct the MET for the towers that lie in the trajectory of the CSC Halo Tracks
  for( edm::View<Candidate>::const_iterator iCandidate = TheCaloTowers->begin() ; iCandidate != TheCaloTowers->end() ; iCandidate++ )
    {
      const Candidate* c = &(*iCandidate);
      if ( c )
        {
          const CaloTower* iTower = dynamic_cast<const CaloTower*> (c);
          if( iTower->et() < TowerEtThreshold ) continue;
          if( abs(iTower->ieta()) > 24 )  continue;   // not in barrel/endcap
          int iphi = iTower->iphi();
          for( unsigned int x = 0 ; x < vEcaliPhi.size() ; x++ )
            {
              if( iphi == vEcaliPhi[x] ) 
                {
                  dMEx += ( TMath::Cos(iTower->phi())*iTower->emEt() );
                  dMEy += ( TMath::Sin(iTower->phi())*iTower->emEt() );
                }
            }
          for( unsigned int x = 0 ; x < vHcaliPhi.size() ; x++ )
            {
              if( iphi == vHcaliPhi[x] ) 
                {
                  dMEx += ( TMath::Cos(iTower->phi() )*iTower->hadEt() ) ;
                  dMEy += ( TMath::Sin(iTower->phi() )*iTower->hadEt() ) ;
                }
            }
        }
    }
  
  TheGlobalHaloData.SetMETCorrections(dMEx, dMEy);
  return TheGlobalHaloData;
}