HFClusterAlgo Class Reference

#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 82 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().

{ m_isMC=isMC; }

bool HFClusterAlgo::isPMTHit ( const HFRecHit &  hfr )

private

Definition at line 409 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 178 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(), create_public_lumi_plots::log, M_PI, max(), AlCaHLTBitMon_ParallelJobs::p, phi, PV3DBase< T, PVType, FrameType >::phi(), evf::utils::sid, python.multivaluedict::sort(), w(), 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 51 of file HFClusterAlgo.cc.

References MCMaterialCorrections, and ZplusMC2010Corrections.

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];
  }
  if (m_correctionSet==2) { // corrections for material from MC + 2010 Z-based ccalibration
    for (int ii=0; ii<13*2; ii++) 
      m_correctionByEta[ii]=ZplusMC2010Corrections[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.