#include <HFClusterAlgo.h>

Classes
struct	HFCompleteHit

Public Member Functions
void	clusterize (const HFRecHitCollection &hf, const CaloGeometry &geom, reco::HFEMClusterShapeCollection &clusters, reco::SuperClusterCollection &SuperClusters)

	HFClusterAlgo ()

void	isMC (bool isMC)

void	setup (double minTowerEnergy, double seedThreshold, double maximumSL, double m_maximumRenergy, bool usePMTflag, bool usePulseflag, bool forcePulseFlagMC, int correctionSet)

Private Member Functions
bool	isPMTHit (const HFRecHit &hfr)

bool	makeCluster (const HcalDetId &seedid, const HFRecHitCollection &hf, const CaloGeometry &geom, reco::HFEMClusterShape &clusShp, reco::SuperCluster &SClus)

Private Attributes
std::vector< double >	m_correctionByEta

int	m_correctionSet

std::vector< double >	m_cutByEta

bool	m_forcePulseFlagMC

bool	m_isMC

double	m_maximumRenergy

double	m_maximumSL

double	m_minTowerEnergy

double	m_seedThreshold

bool	m_usePMTFlag

bool	m_usePulseFlag

Friends
class	CompareHFCompleteHitET

class	CompareHFCore

Detailed Description

Author: K. Klapoetke – Minnesota; Kevin Klapoetke (Minnesota)

$Id:version 1.2

Definition at line 22 of file HFClusterAlgo.h.

Constructor & Destructor Documentation

HFClusterAlgo::HFClusterAlgo ( )

Definition at line 19 of file HFClusterAlgo.cc.

 {
   m_isMC=true; // safest
   
 }

Member Function Documentation

void HFClusterAlgo::clusterize	(	const HFRecHitCollection &	hf,
		const CaloGeometry &	geom,
		reco::HFEMClusterShapeCollection &	clusters,
		reco::SuperClusterCollection &	SuperClusters
	)

Analyze the hits

Definition at line 75 of file HFClusterAlgo.cc.

References edm::SortedCollection< T, SORT >::begin(), edm::SortedCollection< T, SORT >::end(), HFClusterAlgo::HFCompleteHit::energy, HFClusterAlgo::HFCompleteHit::et, eta(), edm::SortedCollection< T, SORT >::find(), relativeConstraints::geom, CaloGeometry::getPosition(), HcalForward, i, HFClusterAlgo::HFCompleteHit::id, indexByEta(), j, gen::k, and python.multivaluedict::sort().

Referenced by HFEMClusterProducer::produce().

                                                           {
   
   std::vector<HFCompleteHit> protoseeds, seeds;
   HFRecHitCollection::const_iterator j,j2;
   std::vector<HFCompleteHit>::iterator i;
   std::vector<HFCompleteHit>::iterator k;
   int dP, dE, PWrap;
   bool isok=true;
   HFEMClusterShape clusShp;
  
   SuperCluster Sclus;
   bool doCluster=false;
   
   for (j=hf.begin(); j!= hf.end(); j++)  {
     const int aieta=j->id().ietaAbs();
     int iz=(aieta-29);
     // only long fibers and not 29,40,41 allowed to be considered as seeds
     if (j->id().depth()!=1) continue;
     if (aieta==40 || aieta==41 || aieta==29) continue;
    
 
     if (iz<0 || iz>12) {
       edm::LogWarning("HFClusterAlgo") << "Strange invalid HF hit: " << j->id();
       continue;
     }
 
     if (m_cutByEta[iz]<0) {
       double eta=geom.getPosition(j->id()).eta();
       m_cutByEta[iz]=m_seedThreshold*cosh(eta); // convert ET to E for this ring
     }
     double elong=j->energy()*m_correctionByEta[indexByEta(j->id())];
     if (elong>m_cutByEta[iz]) {
       j2=hf.find(HcalDetId(HcalForward,j->id().ieta(),j->id().iphi(),2));
       double eshort=(j2==hf.end())?(0):(j2->energy());
       if (j2!=hf.end())
          eshort*=m_correctionByEta[indexByEta(j2->id())];
       if (((elong-eshort)/(elong+eshort))>m_maximumSL) continue;
       //if ((m_usePMTFlag)&&(j->flagField(4,1))) continue;
       //if ((m_usePulseFlag)&&(j->flagField(1,1))) continue;
       if(isPMTHit(*j)) continue;
 
       HFCompleteHit ahit;
       double eta=geom.getPosition(j->id()).eta();
       ahit.id=j->id();
       ahit.energy=elong;
       ahit.et=ahit.energy/cosh(eta);
       protoseeds.push_back(ahit);
     }
   }
 
   if(!protoseeds.empty()){   
     std::sort(protoseeds.begin(), protoseeds.end(), CompareHFCompleteHitET());
     for (i=protoseeds.begin(); i!= protoseeds.end(); i++)  {
       isok=true;
       doCluster=false;
 
       if ( (i==protoseeds.begin()) && (isok) ) {
     doCluster=true;
       }else {
     // check for overlap with existing clusters 
     for (k=seeds.begin(); isok && k!=seeds.end(); k++) { //i->hits, k->seeds
       
       for (dE=-2; dE<=2; dE++)
         for (dP=-4;dP<=4; dP+=2) {
           PWrap=k->id.iphi()+dP;    
           if (PWrap<0) 
         PWrap+=72;
           if (PWrap>72)
         PWrap-=72;
           
           if ( (i->id.iphi()==PWrap) && (i->id.ieta()==k->id.ieta()+dE))
         isok = false;
         }
     }
     if (isok) {
       doCluster=true;
     }
       }
       if (doCluster) { 
     seeds.push_back(*i);
 
     bool clusterOk=makeCluster( i->id(),hf, geom,clusShp,Sclus);
     if (clusterOk) { // cluster is _not_ ok if seed is rejected due to other cuts
       clusterShapes.push_back(clusShp);
       SuperClusters.push_back(Sclus);
     }
 
       }
     }//end protoseed loop
   }//end if seeCount
 }

void HFClusterAlgo::isMC ( bool isMC )

inline

Definition at line 28 of file HFClusterAlgo.h.

References isMC(), and m_isMC.

Referenced by isMC(), and HFEMClusterProducer::produce().

28 { m_isMC=isMC; }

HFClusterAlgo::isMC

void isMC(bool isMC)

Definition: HFClusterAlgo.h:28

HFClusterAlgo::m_isMC

bool m_isMC

Definition: HFClusterAlgo.h:44

bool HFClusterAlgo::isPMTHit ( const HFRecHit & hfr )

private

Definition at line 402 of file HFClusterAlgo.cc.

References CaloRecHit::flagField(), HcalCaloFlagLabels::HFDigiTime, and HcalCaloFlagLabels::HFLongShort.

                                                {
 
   bool pmthit=false;
 
   if((hfr.flagField(HcalCaloFlagLabels::HFLongShort))&&(m_usePMTFlag)) pmthit=true;
   if (!(m_isMC && !m_forcePulseFlagMC))
     if((hfr.flagField(HcalCaloFlagLabels::HFDigiTime))&&(m_usePulseFlag)) pmthit=true;
  
   return pmthit;
 
 }

bool HFClusterAlgo::makeCluster	(	const HcalDetId &	seedid,
		const HFRecHitCollection &	hf,
		const CaloGeometry &	geom,
		reco::HFEMClusterShape &	clusShp,
		reco::SuperCluster &	SClus
	)

private

Definition at line 171 of file HFClusterAlgo.cc.

References abs, edm::SortedCollection< T, SORT >::end(), PV3DBase< T, PVType, FrameType >::eta(), eta(), edm::SortedCollection< T, SORT >::find(), CaloGeometry::getPosition(), HcalForward, i, HcalDetId::ieta(), HcalDetId::ietaAbs(), indexByEta(), HcalDetId::iphi(), funct::log(), M_PI, max(), L1TEmulatorMonitor_cff::p, phi, PV3DBase< T, PVType, FrameType >::phi(), evf::utils::sid, python.multivaluedict::sort(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::z(), and HcalDetId::zside().

                                       {
             
 
   double w=0;//sum over all log E's
   double wgt=0;
   double w_e=0;//sum over ieat*energy
   double w_x=0;
   double w_y=0;
   double w_z=0;
   double wp_e=0;//sum over iphi*energy
   double e_e=0;//nonwieghted eta sum
   double e_ep=0; //nonweighted phi sum
  
   double l_3=0;//sum for enenergy in 3x3 long fibers etc.
   double s_3=0;
   double l_5=0;
   double s_5=0;
   double l_1=0;
   double s_1=0;
   int de, dp, phiWrap;
   double l_1e=0;
   GlobalPoint sp=geom.getPosition(seedid);
   std::vector<double> coreCanid;
   std::vector<double>::const_iterator ci;
   HFRecHitCollection::const_iterator i,is,il;
   std::vector<DetId> usedHits; 
  
   HFRecHitCollection::const_iterator si;
   HcalDetId sid(HcalForward,seedid.ieta(),seedid.iphi(),1);
   si=hf.find(sid);  
 
   bool clusterOk=true; // assume the best to start...
 
   // lots happens here
   // edge type 1 has 40/41 in 3x3 and 5x5
   bool edge_type1=seedid.ietaAbs()==39 && (seedid.iphi()%4)==3;
 
   double e_seed=si->energy()*m_correctionByEta[indexByEta(si->id())];
   
   for (de=-2; de<=2; de++)
     for (dp=-4;dp<=4; dp+=2) {
       phiWrap=seedid.iphi()+dp; 
       if (phiWrap<0) 
         phiWrap+=72;
       if (phiWrap>72)
         phiWrap-=72;
 
   
         /* Handling of phi-width change problems */
         if (edge_type1 && de==seedid.zside()) {
           if (dp==-2) { // we want it in the 3x3
             phiWrap-=2;
             if (phiWrap<0) 
               phiWrap+=72;
           } 
           else if (dp==-4) {
             continue; // but not double counted in 5x5
           }
         }
 
     HcalDetId idl(HcalForward,seedid.ieta()+de,phiWrap,1);
     HcalDetId ids(HcalForward,seedid.ieta()+de,phiWrap,2);
 
     
     il=hf.find(idl);
     is=hf.find(ids);        
 
 
 
 
     double e_long=1.0; 
     double e_short=0.0; 
     if (il!=hf.end()) e_long=il->energy()*m_correctionByEta[indexByEta(il->id())];
     if (e_long <= m_minTowerEnergy) e_long=0.0;
     if (is!=hf.end()) e_short=is->energy()*m_correctionByEta[indexByEta(is->id())];
     if (e_short <= m_minTowerEnergy) e_short=0.0;
     double eRatio=(e_long-e_short)/std::max(1.0,(e_long+e_short));
     
     // require S/L > a minimum amount for inclusion
     if ((abs(eRatio) > m_maximumSL)&&(std::max(e_long,e_short) > m_maximumRenergy)) {
       if (dp==0 && de==0) clusterOk=false; // somehow, the seed is hosed
       continue;
     }
      
     if((il!=hf.end())&&(isPMTHit(*il))){
       if (dp==0 && de==0) clusterOk=false; // somehow, the seed is hosed
       continue;//continue to next hit, do not include this one in cluster
     }
     
 
      /*
 
      old pmt code
      // cut on "PMT HIT" flag
      if ((il!=hf.end())&&(il->flagField(4,1))&&(m_usePMTFlag)) {//HFPET flag for lone/short doil->flagField(0,1)
      if (dp==0 && de==0) clusterOk=false; // somehow, the seed is hosed
      continue;//continue to next hit, do not include this one in cluster
      }
      
      // cut on "Pulse shape HIT" flag
      if ((il!=hf.end())&&(il->flagField(1,1))&&(m_usePulseFlag)) {//HF DIGI TIME flag
      if (dp==0 && de==0) clusterOk=false; // somehow, the seed is hosed
      continue;//continue to next hit, do not include this one in cluster
      }
      */
  
 
 
 
     if (e_long > m_minTowerEnergy && il!=hf.end()) {
 
       // record usage
       usedHits.push_back(idl.rawId());
       // always in the 5x5
       l_5+=e_long;
       // maybe in the 3x3
           if ((de>-2)&&(de<2)&&(dp>-4)&&(dp<4)) {
             l_3+=e_long;
           }
       // sometimes in the 1x1
       if ((dp==0)&&(de==0)) {
             l_1=e_long;
           }
 
       // maybe in the core?
       if ((de>-2)&&(de<2)&&(dp>-4)&&(dp<4)&&(e_long>(.5*e_seed))) {
             coreCanid.push_back(e_long);
           }
       
       // position calculation
       GlobalPoint p=geom.getPosition(idl);
           
           double d_p = p.phi()-sp.phi();
           while (d_p < -M_PI)
             d_p+=2*M_PI;
           while (d_p > M_PI)
             d_p-=2*M_PI;
           double d_e = p.eta()-sp.eta();
       
       wgt=log((e_long));
       if (wgt>0){
         w+=wgt;
         w_e+=(d_e)*wgt;
         wp_e+=(d_p)*wgt;
         e_e+=d_e;
         e_ep+=d_p;
        
         w_x+=(p.x())*wgt;//(p.x()-sp.x())*wgt;
         w_y+=(p.y())*wgt;
               w_z+=(p.z())*wgt;
       }
     } else {
       if (dp==0 && de==0) clusterOk=false; // somehow, the seed is hosed
     }
     
     if (e_short > m_minTowerEnergy && is!=hf.end()) {
       // record usage
       usedHits.push_back(ids.rawId());
       // always in the 5x5
       s_5+=e_short;
       // maybe in the 3x3
       if ((de>-2)&&(de<2)&&(dp>-4)&&(dp<4)) {
             s_3+=e_short;
           }
       // sometimes in the 1x1
       if ((dp==0)&&(de==0)) {
             s_1=e_short;
           }
     }
     }
 
 
   if (!clusterOk) return false;
 
   //Core sorting done here
   std::sort(coreCanid.begin(), coreCanid.end(), CompareHFCore());
   for (ci=coreCanid.begin();ci!=coreCanid.end();ci++){
     if(ci==coreCanid.begin()){
       l_1e=*ci;
     }else if (*ci>.5*l_1e){
       l_1e+=*ci;
     }
   }//core sorting end 
   
   double z_=w_z/w;    
   double x_=w_x/w;
   double y_=w_y/w;
   
   //calcualte position, final
   double eta=w_e/w+sp.eta();
   
   double phi=(wp_e/w)+sp.phi();
   
   while (phi < -M_PI)
     phi+=2*M_PI;
   while (phi > M_PI)
     phi-=2*M_PI;
   
   //calculate cell phi and cell eta
   static const double HFEtaBounds[14] = {2.853, 2.964, 3.139, 3.314, 3.489, 3.664, 3.839, 4.013, 4.191, 4.363, 4.538, 4.716, 4.889, 5.191};
   double RcellEta = fabs(eta);
   double Cphi = (phi>0.)?(fmod((phi),0.087*2)/(0.087*2)):((fmod((phi),0.087*2)/(0.087*2))+1.0);
   double Rbin = -1.0;
   for (int icell = 0; icell < 12; icell++ ){
     if ( (RcellEta>HFEtaBounds[icell]) && (RcellEta<HFEtaBounds[icell+1]) )
       Rbin = (RcellEta - HFEtaBounds[icell])/(HFEtaBounds[icell+1] - HFEtaBounds[icell]);
   }
   double Ceta=Rbin;
   
   while (phi< -M_PI)
     phi+=2*M_PI;
   while (phi > M_PI)
     phi-=2*M_PI;
   
   
   math::XYZPoint xyzclus(x_,y_,z_);
   
   //return  HFEMClusterShape, SuperCluster
   HFEMClusterShape myClusShp(l_1, s_1, l_3, s_3, l_5,s_5, l_1e,Ceta, Cphi,seedid);
   clusShp = myClusShp;
       
   SuperCluster MySclus(l_3,xyzclus);
   Sclus=MySclus;
 
   return clusterOk;
   
 }

void HFClusterAlgo::setup	(	double	minTowerEnergy,
		double	seedThreshold,
		double	maximumSL,
		double	m_maximumRenergy,
		bool	usePMTflag,
		bool	usePulseflag,
		bool	forcePulseFlagMC,
		int	correctionSet
	)

Definition at line 48 of file HFClusterAlgo.cc.

References MCMaterialCorrections.

Referenced by HFEMClusterProducer::HFEMClusterProducer().

                                                                                           {
   m_seedThreshold=seedThreshold;
   m_minTowerEnergy=minTowerEnergy;
   m_maximumSL=maximumSL;
   m_usePMTFlag=usePMTflag;
   m_usePulseFlag=usePulseflag;
   m_forcePulseFlagMC=forcePulseFlagMC;
   m_maximumRenergy=maximumRenergy;
   m_correctionSet = correctionSet;
 
   for(int ii=0;ii<13;ii++){
     m_cutByEta.push_back(-1);
   }
 
   // always set all the corrections to one...
   for (int ii=0; ii<13*2; ii++) 
     m_correctionByEta.push_back(1.0);
 
   if (m_correctionSet==1) { // corrections for material from MC
     for (int ii=0; ii<13*2; ii++) 
       m_correctionByEta[ii]=MCMaterialCorrections[ii];
   }
 
 }

Friends And Related Function Documentation

friend class CompareHFCompleteHitET

friend

Definition at line 38 of file HFClusterAlgo.h.

friend class CompareHFCore

friend

Definition at line 39 of file HFClusterAlgo.h.

Member Data Documentation

std::vector<double> HFClusterAlgo::m_correctionByEta

private

Definition at line 47 of file HFClusterAlgo.h.

int HFClusterAlgo::m_correctionSet

private

Definition at line 45 of file HFClusterAlgo.h.

std::vector<double> HFClusterAlgo::m_cutByEta

private

Definition at line 46 of file HFClusterAlgo.h.

bool HFClusterAlgo::m_forcePulseFlagMC

private

Definition at line 43 of file HFClusterAlgo.h.

bool HFClusterAlgo::m_isMC

private

Definition at line 44 of file HFClusterAlgo.h.

Referenced by isMC().

double HFClusterAlgo::m_maximumRenergy

private

Definition at line 41 of file HFClusterAlgo.h.

double HFClusterAlgo::m_maximumSL

private

Definition at line 41 of file HFClusterAlgo.h.

double HFClusterAlgo::m_minTowerEnergy

private

Definition at line 41 of file HFClusterAlgo.h.

double HFClusterAlgo::m_seedThreshold

private

Definition at line 41 of file HFClusterAlgo.h.

bool HFClusterAlgo::m_usePMTFlag

private

Definition at line 42 of file HFClusterAlgo.h.

bool HFClusterAlgo::m_usePulseFlag

private

Definition at line 43 of file HFClusterAlgo.h.

Classes

Public Member Functions

Private Member Functions

Private Attributes

Friends

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Friends And Related Function Documentation

Member Data Documentation