#include <EcalHaloAlgo.h>

Public Member Functions
reco::EcalHaloData	Calculate (const CaloGeometry &TheCaloGeometry, edm::Handle< reco::PhotonCollection > &ThePhotons, edm::Handle< reco::SuperClusterCollection > &TheSuperClusters, edm::Handle< EBRecHitCollection > &TheEBRecHits, edm::Handle< EERecHitCollection > &TheEERecHits, edm::Handle< ESRecHitCollection > &TheESRecHits)

	EcalHaloAlgo ()

float	GetAngleCut ()

float	GetEBRecHitEnergyThreshold ()

float	GetEERecHitEnergyThreshold ()

float	GetESRecHitEnergyThreshold ()

float	GetPhiWedgeEnergyThreshold ()

int	GetPhiWedgeNHitsThreshold ()

float	GetRoundnessCut ()

void	SetAngleCut (float a=4.)

void	SetPhiWedgeEnergyThreshold (float SumE)

void	SetPhiWedgeNHitsThreshold (int nhits)

void	SetPhiWedgeThresholds (float SumE, int nhits)

void	SetRecHitEnergyThresholds (float EB, float EE, float ES)

void	SetRoundnessCut (float r=100.)

	~EcalHaloAlgo ()

Private Attributes
float	AngleCut

float	EBRecHitEnergyThreshold

float	EERecHitEnergyThreshold

float	ESRecHitEnergyThreshold

int	NHitsThreshold

float	RoundnessCut

float	SumEnergyThreshold

Detailed Description

Definition at line 42 of file EcalHaloAlgo.h.

Constructor & Destructor Documentation

EcalHaloAlgo::EcalHaloAlgo ( )

Definition at line 17 of file EcalHaloAlgo.cc.

 {
   RoundnessCut = 0 ;
   AngleCut = 0 ;
   EBRecHitEnergyThreshold = 0.;
   EERecHitEnergyThreshold = 0.;
   ESRecHitEnergyThreshold = 0.;
   SumEnergyThreshold = 0.;
   NHitsThreshold =0;
 }

EcalHaloAlgo::~EcalHaloAlgo ( )

inline

Definition at line 48 of file EcalHaloAlgo.h.

48 {}

Member Function Documentation

EcalHaloData EcalHaloAlgo::Calculate	(	const CaloGeometry &	TheCaloGeometry,
		edm::Handle< reco::PhotonCollection > &	ThePhotons,
		edm::Handle< reco::SuperClusterCollection > &	TheSuperClusters,
		edm::Handle< EBRecHitCollection > &	TheEBRecHits,
		edm::Handle< EERecHitCollection > &	TheEERecHits,
		edm::Handle< ESRecHitCollection > &	TheESRecHits
	)

Definition at line 28 of file EcalHaloAlgo.cc.

References funct::abs(), Abs(), angle(), CompareTime(), DetId::Ecal, reco::EcalHaloData::GetPhiWedges(), reco::EcalHaloData::GetShowerShapesAngle(), reco::EcalHaloData::GetShowerShapesRoundness(), CaloGeometry::getSubdetectorGeometry(), reco::EcalHaloData::GetSuperClusters(), i, hit::id, EBDetId::ieta(), edm::helper::Filler< Map >::insert(), EBDetId::iphi(), edm::HandleBase::isValid(), j, edm::Handle< T >::product(), edm::RefVector< C, T, F >::push_back(), DetId::rawId(), and reco::PhiWedge::SetPlusZOriginConfidence().

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

 {  
   EcalHaloData TheEcalHaloData;
 
   // Store energy sum of rechits as a function of iPhi (iphi goes from 1 to 72)       
   float SumE[361];           
   // Store number of rechits as a function of iPhi
   int NumHits[361];               
   // Store minimum time of rechit as a function of iPhi
   float MinTimeHits[361];                             
   // Store maximum time of rechit as a function of iPhi 
   float MaxTimeHits[361];                      
 
   // initialize
   for(int i = 0 ; i < 361 ; i++ )
     {
       SumE[i] = 0.;
       NumHits[i] = 0 ;
       MinTimeHits[i] = 9999.;
       MaxTimeHits[i] = -9999.;
     }
 
   // Loop over EB RecHits
   for(EBRecHitCollection::const_iterator hit = TheEBRecHits->begin() ; hit != TheEBRecHits->end() ; hit++ )
     {
       // Arbitrary threshold to kill noise (needs to be optimized with data)
       if (hit->energy() < EBRecHitEnergyThreshold ) continue;
       
       // Get Det Id of the rechit
       DetId id = DetId(hit->id()); 
       const CaloSubdetectorGeometry* TheSubGeometry = 0;                                                                         
       const CaloCellGeometry* cell = 0 ;                                                                                    
 
       // Get EB geometry 
       TheSubGeometry = TheCaloGeometry.getSubdetectorGeometry(DetId::Ecal, 1);                                                         
       EBDetId EcalID(id.rawId());
       if( TheSubGeometry )
     cell = TheSubGeometry->getGeometry(id);
   
       if(cell)
     {
       // GlobalPoint globalpos = cell->getPosition();
       //      float r = TMath::Sqrt ( globalpos.y()*globalpos.y() + globalpos.x()*globalpos.x());
       int iPhi = EcalID.iphi();
 
       if( iPhi < 361 )    // just to be safe
         {
           //iPhi = (iPhi-1)/5 +1;  // convert ecal iphi to phiwedge iphi  (e.g. there are 5 crystal per phi wedge, as in calotowers )
           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++ )
   for( int iPhi = 1 ; iPhi < 361; iPhi++ )
     {
       if( SumE[iPhi] >= SumEnergyThreshold && NumHits[iPhi] > NHitsThreshold )
     {
       // Build PhiWedge and store to EcalHaloData 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
       
       // Loop over EB RecHits
       std::vector<const EcalRecHit*>  Hits;
       for(EBRecHitCollection::const_iterator hit = TheEBRecHits->begin() ; hit != TheEBRecHits->end() ; hit++ )
         {
           if (hit->energy() < EBRecHitEnergyThreshold ) continue;
           
           // Get Det Id of the rechit
           DetId id = DetId(hit->id()); 
           EBDetId EcalID(id.rawId());
           int Hit_iPhi = EcalID.iphi();
           //Hit_iPhi = (Hit_iPhi-1)/5 +1; // convert ecal iphi to phiwedge iphi
           if( Hit_iPhi != iPhi ) 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++ )
         {
           DetId id_i = DetId(Hits[i]->id()); 
           EBDetId EcalID_i(id_i.rawId());
           int ieta_i = EcalID_i.ieta();
           for( unsigned int j = (i+1) ; j < Hits.size() ; j++ )
         {
           DetId id_j = DetId(Hits[j]->id() );
           EBDetId EcalID_j(id_j.rawId());
           int ieta_j = EcalID_j.ieta();
           if( ieta_i > ieta_j ) PlusToMinus += TMath::Abs(ieta_j - ieta_i );
           else MinusToPlus += TMath::Abs(ieta_j -ieta_i) ;
         }
         }
       
       float PlusZOriginConfidence = (PlusToMinus+MinusToPlus) ? PlusToMinus / (PlusToMinus+MinusToPlus) : -1.;
       wedge.SetPlusZOriginConfidence(PlusZOriginConfidence);
       TheEcalHaloData.GetPhiWedges().push_back(wedge);
     }
     }
 
   std::vector<float> vShowerShapes_Roundness;
   std::vector<float> vShowerShapes_Angle ;
   for(reco::SuperClusterCollection::const_iterator cluster = TheSuperClusters->begin() ; cluster != TheSuperClusters->end() ; cluster++ )
     {
       if( abs(cluster->eta()) <= 1.48  )
     { 
       vector<float> shapes = EcalClusterTools::roundnessBarrelSuperClusters( *cluster, (*TheEBRecHits.product()));
       float roundness = shapes[0];
       float angle = shapes[1];
       
       // Check if supercluster belongs to photon and passes the cuts on roundness and angle, if so store the reference to it
       if( (roundness >=0 && roundness < GetRoundnessCut()) &&  angle >= 0 && angle < GetAngleCut() )
         {
           edm::Ref<SuperClusterCollection> TheClusterRef( TheSuperClusters, cluster - TheSuperClusters->begin());
           bool BelongsToPhoton = false;
           if( ThePhotons.isValid() )
         {
           for(reco::PhotonCollection::const_iterator iPhoton = ThePhotons->begin() ; iPhoton != ThePhotons->end() ; iPhoton++ )
             {
               if(iPhoton->isEB())
             if ( TheClusterRef ==  iPhoton->superCluster() )
               {
                 BelongsToPhoton = true;
                 break;
               }
             }
         }
           //Only store refs to suspicious EB SuperClusters which belong to Photons
           //Showershape variables are more discriminating for these cases  
           if( BelongsToPhoton )
         { 
           TheEcalHaloData.GetSuperClusters().push_back( TheClusterRef ) ; 
         }
         }
       vShowerShapes_Roundness.push_back(shapes[0]);
       vShowerShapes_Angle.push_back(shapes[1]);
     }
       else
     { 
       vShowerShapes_Roundness.push_back(-1.);
       vShowerShapes_Angle.push_back(-1.);
     }
     }
   
   edm::ValueMap<float>::Filler TheRoundnessFiller( TheEcalHaloData.GetShowerShapesRoundness() );
   TheRoundnessFiller.insert( TheSuperClusters, vShowerShapes_Roundness.begin(), vShowerShapes_Roundness.end() );
   TheRoundnessFiller.fill();  
 
   edm::ValueMap<float>::Filler TheAngleFiller( TheEcalHaloData.GetShowerShapesAngle() );
   TheAngleFiller.insert( TheSuperClusters, vShowerShapes_Angle.begin() , vShowerShapes_Angle.end() );
   TheAngleFiller.fill();
 
   return TheEcalHaloData;
 }