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	resetForRun ()
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
std::vector< double >	m_seedmnEta
std::vector< double >	m_seedmnPhi
std::vector< double >	m_seedMXeta
std::vector< double >	m_seedMXphi
double	m_seedThreshold
bool	m_usePMTFlag
bool	m_usePulseFlag

Friends
class	CompareHFCompleteHitET
class	CompareHFCore

Detailed Description

Author

Kevin Klapoetke (Minnesota)

$Id:version 1.2

Author

K. Klapoetke – Minnesota

Definition at line 22 of file HFClusterAlgo.h.

Constructor & Destructor Documentation

HFClusterAlgo::HFClusterAlgo

(

)

Definition at line 20 of file HFClusterAlgo.cc.

{
  m_isMC=true; // safest
  
}

Member Function Documentation

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

Analyze the hits

Definition at line 91 of file HFClusterAlgo.cc.

References edm::SortedCollection< T, SORT >::begin(), edm::SortedCollection< T, SORT >::end(), HFClusterAlgo::HFCompleteHit::energy, HFClusterAlgo::HFCompleteHit::et, PVValHelper::eta, edm::SortedCollection< T, SORT >::find(), relativeConstraints::geom, CaloGeometry::getGeometry(), CaloGeometry::getPosition(), HcalForward, photonIsolationHIProducer_cfi::hf, mps_fire::i, HFClusterAlgo::HFCompleteHit::id, indexByEta(), gen::k, SimDataFormats::CaloAnalysis::sc, and EZArrayFL< T >::size().

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

                                                          {
  
  const CaloGeometry* geom(&geomO);
  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
      auto ccg=geom->getGeometry(j->id()); 
      const CaloCellGeometry::CornersVec& CellCorners=ccg->getCorners();
      for(size_t sc=0;sc<CellCorners.size();sc++){
    if(fabs(CellCorners[sc].z())<1200){
      if(fabs(CellCorners[sc].eta())<m_seedmnEta[iz])m_seedmnEta[iz]=fabs(CellCorners[sc].eta());
      if(fabs(CellCorners[sc].eta())>m_seedMXeta[iz])m_seedMXeta[iz]=fabs(CellCorners[sc].eta());
    }
      }
    }
    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 clusterize(), fastPrimaryVertexProducer_cfi::clusters, relativeConstraints::geom, photonIsolationHIProducer_cfi::hf, isMC(), m_isMC, and resetForRun().

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

{ m_isMC=isMC; }

bool HFClusterAlgo::isPMTHit ( const HFRecHit &  hfr )

private

Definition at line 399 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 196 of file HFClusterAlgo.cc.

References funct::abs(), edm::SortedCollection< T, SORT >::end(), PVValHelper::eta, PV3DBase< T, PVType, FrameType >::eta(), edm::SortedCollection< T, SORT >::find(), CaloGeometry::getPosition(), HcalForward, training_settings::idx, HcalDetId::ieta(), HcalDetId::ietaAbs(), indexByEta(), HcalDetId::iphi(), cmsBatch::log, M_PI, SiStripPI::max, AlCaHLTBitMon_ParallelJobs::p, PV3DBase< T, PVType, FrameType >::phi(), 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;
  const GlobalPoint& sp=geom->getPosition(seedid);
  std::vector<double> coreCanid;
  std::vector<double>::const_iterator ci;
  HFRecHitCollection::const_iterator 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
      }
    




      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
      const 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;
  math::XYZPoint xyzclus(x_,y_,z_);
  //calcualte position, final
  double eta=xyzclus.eta();//w_e/w+sp.eta();
  
  double phi=xyzclus.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
  int idx= fabs(seedid.ieta())-29;
  int ipx=seedid.iphi();
  double Cphi =(phi-m_seedmnPhi[ipx])/(m_seedMXphi[ipx]-m_seedmnPhi[ipx]);
  double Ceta=(fabs(eta)- m_seedmnEta[idx])/(m_seedMXeta[idx]-m_seedmnEta[idx]);
    
  //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::resetForRun

(

)

Definition at line 409 of file HFClusterAlgo.cc.

References cuy::ii.

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

                                {
  edm::LogInfo("HFClusterAlgo")<<"Resetting for Run!";
  for(int ii=0;ii<13;ii++){
    m_cutByEta.push_back(-1);
    m_seedmnEta.push_back(99);
    m_seedMXeta.push_back(-1);
  }
}

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 hfClusterShapes_cfi::forcePulseFlagMC, cuy::ii, M_PI, hfClusterShapes_cfi::maximumRenergy, hfClusterShapes_cfi::maximumSL, MCMaterialCorrections_3XX, hfClusterShapes_cfi::minTowerEnergy, and UEAnalysisJets_cfi::seedThreshold.

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);
    m_seedmnEta.push_back(99);
    m_seedMXeta.push_back(-1);
  }
 for(int ii=0;ii<73;ii++){
    double minphi=0.0872664*(ii-1);
    double maxphi=0.0872664*(ii+1);
    while (minphi < -M_PI)
      minphi+=2*M_PI;
    while (minphi > M_PI)
      minphi-=2*M_PI;
    while (maxphi < -M_PI)
      maxphi+=2*M_PI;
    while (maxphi > M_PI)
      maxphi-=2*M_PI;
    if(ii==37)
      minphi=-3.1415904;
    m_seedmnPhi.push_back(minphi);
    m_seedMXphi.push_back(maxphi);
 }
 
 // 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_3XX[ii];
 } 
}

Friends And Related Function Documentation

friend class CompareHFCompleteHitET

friend

Definition at line 40 of file HFClusterAlgo.h.

friend class CompareHFCore

friend

Definition at line 41 of file HFClusterAlgo.h.

Member Data Documentation

std::vector<double> HFClusterAlgo::m_correctionByEta

private

Definition at line 49 of file HFClusterAlgo.h.

int HFClusterAlgo::m_correctionSet

private

Definition at line 47 of file HFClusterAlgo.h.

std::vector<double> HFClusterAlgo::m_cutByEta

private

Definition at line 48 of file HFClusterAlgo.h.

bool HFClusterAlgo::m_forcePulseFlagMC

private

Definition at line 45 of file HFClusterAlgo.h.

bool HFClusterAlgo::m_isMC

private

Definition at line 46 of file HFClusterAlgo.h.

Referenced by isMC().

double HFClusterAlgo::m_maximumRenergy

private

Definition at line 43 of file HFClusterAlgo.h.

double HFClusterAlgo::m_maximumSL

private

Definition at line 43 of file HFClusterAlgo.h.

double HFClusterAlgo::m_minTowerEnergy

private

Definition at line 43 of file HFClusterAlgo.h.

std::vector<double> HFClusterAlgo::m_seedmnEta

private

Definition at line 50 of file HFClusterAlgo.h.

std::vector<double> HFClusterAlgo::m_seedmnPhi

private

Definition at line 52 of file HFClusterAlgo.h.

std::vector<double> HFClusterAlgo::m_seedMXeta

private

Definition at line 51 of file HFClusterAlgo.h.

std::vector<double> HFClusterAlgo::m_seedMXphi

private

Definition at line 53 of file HFClusterAlgo.h.

double HFClusterAlgo::m_seedThreshold

private

Definition at line 43 of file HFClusterAlgo.h.

bool HFClusterAlgo::m_usePMTFlag

private

Definition at line 44 of file HFClusterAlgo.h.

bool HFClusterAlgo::m_usePulseFlag

private

Definition at line 45 of file HFClusterAlgo.h.