#include <GlobalHaloAlgo.h>

Public Member Functions
reco::GlobalHaloData	Calculate (const CaloGeometry &TheCaloGeometry, const CSCGeometry &TheCSCGeometry, const reco::CaloMET &TheCaloMET, edm::Handle< edm::View< reco::Candidate > > &TheCaloTowers, edm::Handle< CSCSegmentCollection > &TheCSCSegments, edm::Handle< CSCRecHit2DCollection > &TheCSCRecHits, const reco::CSCHaloData &TheCSCHaloData, const reco::EcalHaloData &TheEcalHaloData, const reco::HcalHaloData &TheHcalHaloData)

	GlobalHaloAlgo ()

void	SetCaloTowerEtThreshold (float EtMin)

void	SetEcalMatchingRadius (float min, float max)

void	SetHcalMatchingRadius (float min, float max)

	~GlobalHaloAlgo ()

Private Attributes
float	Ecal_R_Max

float	Ecal_R_Min

float	Hcal_R_Max

float	Hcal_R_Min

float	TowerEtThreshold

Detailed Description

Definition at line 57 of file GlobalHaloAlgo.h.

Constructor & Destructor Documentation

GlobalHaloAlgo::GlobalHaloAlgo ( )

Definition at line 30 of file GlobalHaloAlgo.cc.

 {
   // 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.
   
 }

GlobalHaloAlgo::~GlobalHaloAlgo ( )

inline

Definition at line 62 of file GlobalHaloAlgo.h.

62 {}

Member Function Documentation

reco::GlobalHaloData GlobalHaloAlgo::Calculate	(	const CaloGeometry &	TheCaloGeometry,
		const CSCGeometry &	TheCSCGeometry,
		const reco::CaloMET &	TheCaloMET,
		edm::Handle< edm::View< reco::Candidate > > &	TheCaloTowers,
		edm::Handle< CSCSegmentCollection > &	TheCSCSegments,
		edm::Handle< CSCRecHit2DCollection > &	TheCSCRecHits,
		const reco::CSCHaloData &	TheCSCHaloData,
		const reco::EcalHaloData &	TheEcalHaloData,
		const reco::HcalHaloData &	TheHcalHaloData
	)

Definition at line 40 of file GlobalHaloAlgo.cc.

References funct::abs(), Abs(), EnergyCorrector::c, MuonSubdetId::CSC, CaloTower::emEt(), CaloTower::et(), reco::CSCHaloData::GetCSCTrackImpactPositions(), reco::GlobalHaloData::GetMatchedEcalPhiWedges(), reco::GlobalHaloData::GetMatchedHcalPhiWedges(), reco::HcalHaloData::GetPhiWedges(), reco::EcalHaloData::GetPhiWedges(), CaloTower::hadEt(), i, CSCGeometry::idToDetUnit(), CaloTower::ieta(), CaloTower::iphi(), edm::HandleBase::isValid(), DetId::Muon, PV3DBase< T, PVType, FrameType >::phi(), reco::LeafCandidate::phi(), Phi_To_EcaliPhi(), Phi_To_HcaliPhi(), reco::LeafCandidate::pt(), alignCSCRings::r, reco::GlobalHaloData::SetMETCorrections(), reco::GlobalHaloData::SetMETOverSumEt(), reco::PhiWedge::SetOverlappingCSCRecHits(), reco::PhiWedge::SetOverlappingCSCSegments(), reco::MET::sumEt(), GeomDet::surface(), x, PV3DBase< T, PVType, FrameType >::x(), detailsBasic3DVector::y, and PV3DBase< T, PVType, FrameType >::y().

Referenced by reco::GlobalHaloDataProducer::produce().

 {
   
   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;
 }