/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoMET/METAlgorithms/src/HcalHaloAlgo.cc

Go to the documentation of this file.

#include "RecoMET/METAlgorithms/interface/HcalHaloAlgo.h"

/*
  [class]:  HcalHaloAlgo
  [authors]: R. Remington, The University of Florida
  [description]: See HcalHaloAlgo.h
  [date]: October 15, 2009
*/

using namespace std;
using namespace edm;
using namespace reco;

#include <iomanip>
bool CompareTime(const HBHERecHit* x, const HBHERecHit* y ){ return x->time() < y->time() ;}

HcalHaloAlgo::HcalHaloAlgo()
{
  HBRecHitEnergyThreshold = 0.;
  HERecHitEnergyThreshold = 0.;
  SumEnergyThreshold = 0.;
  NHitsThreshold = 0;
}

HcalHaloData HcalHaloAlgo::Calculate(const CaloGeometry& TheCaloGeometry, edm::Handle<HBHERecHitCollection>& TheHBHERecHits)
{
  HcalHaloData TheHcalHaloData;
  
  // Store Energy sum of rechits as a function of iPhi (iPhi goes from 1 to 72)
  float SumE[73];
  // Store Number of rechits as a function of iPhi 
  int NumHits[73];
  // Store minimum time of rechit as a function of iPhi
  float MinTimeHits[73];
  // Store maximum time of rechit as a function of iPhi
  float MaxTimeHits[73];
  for(unsigned int i = 0 ; i < 73 ; i++ ) 
    {
      SumE[i] = 0;
      NumHits[i]= 0;
      MinTimeHits[i] = 0.;
      MaxTimeHits[i] = 0.;
    }
  
  for( HBHERecHitCollection::const_iterator hit = TheHBHERecHits->begin() ; hit != TheHBHERecHits->end() ; hit++ )
    {
      HcalDetId id = HcalDetId(hit->id());                                                                                                    
      switch ( id.subdet() )                                                                                         
        {      
        case HcalBarrel:                                                                           
          if(hit->energy() < HBRecHitEnergyThreshold )continue;
          break;                                                                                                                  
        case HcalEndcap:                                                                                          
          if(hit->energy() < HERecHitEnergyThreshold ) continue;
          break;
        default:
          continue;
        }
      
      int iEta = id.ieta();
      int iPhi = id.iphi();
      if(iPhi < 73 && TMath::Abs(iEta) < 23 )
        { 
          SumE[iPhi]+= hit->energy();
          NumHits[iPhi] ++;
          
          float time = hit->time();
          MinTimeHits[iPhi] = time < MinTimeHits[iPhi] ? time : MinTimeHits[iPhi];
          MaxTimeHits[iPhi] = time > MaxTimeHits[iPhi] ? time : MaxTimeHits[iPhi];
        }
    }
  
  for( int iPhi = 1 ; iPhi < 73 ; iPhi++ )
    {
      if( SumE[iPhi] >= SumEnergyThreshold && NumHits[iPhi] > NHitsThreshold )
        {
          // Build PhiWedge and store to HcalHaloData if energy or #hits pass thresholds
          PhiWedge wedge(SumE[iPhi], iPhi, NumHits[iPhi], MinTimeHits[iPhi], MaxTimeHits[iPhi]);
          
          // Loop over rechits again to calculate direction based on timing info
          std::vector<const HBHERecHit*> Hits;
          for( HBHERecHitCollection::const_iterator hit = TheHBHERecHits->begin() ; hit != TheHBHERecHits->end() ; hit++ )
            {

              HcalDetId id = HcalDetId(hit->id());
              if( id.iphi() != iPhi ) continue;
              if( TMath::Abs(id.ieta() ) > 22 ) continue;  // has to overlap geometrically w/ HB
              switch ( id.subdet() )
                {
                case HcalBarrel:
                  if(hit->energy() < HBRecHitEnergyThreshold )continue;
                  break;
                case HcalEndcap:
                  if(hit->energy() < HERecHitEnergyThreshold ) continue;
                  break;
                default:
                  continue;
                }
              Hits.push_back(&(*hit));
            }
              
          std::sort( Hits.begin() , Hits.end() , CompareTime);
          float MinusToPlus = 0.;
          float PlusToMinus = 0.;
          for( unsigned int i = 0 ; i < Hits.size() ; i++ )
            {
              HcalDetId id_i = HcalDetId(Hits[i]->id() );
              int ieta_i = id_i.ieta();
              for( unsigned int j = (i+1) ; j < Hits.size() ; j++ )
                {
                  HcalDetId id_j = HcalDetId(Hits[j]->id() );
                  int ieta_j = id_j.ieta();
                  if( ieta_i > ieta_j ) PlusToMinus += TMath::Abs(ieta_i - ieta_j ) ;
                  else MinusToPlus += TMath::Abs(ieta_i - ieta_j);
                }
            }
          float PlusZOriginConfidence = (PlusToMinus + MinusToPlus )? PlusToMinus / ( PlusToMinus + MinusToPlus ) : -1. ;
          wedge.SetPlusZOriginConfidence( PlusZOriginConfidence );
          TheHcalHaloData.GetPhiWedges().push_back( wedge );
        }
    }
  return TheHcalHaloData;
  
}