#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 18 of file HFClusterAlgo.h.

Constructor & Destructor Documentation

HFClusterAlgo::HFClusterAlgo ( )

Definition at line 16 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 86 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, mps_fire::i, HFClusterAlgo::HFCompleteHit::id, indexByEta(), dqmiolumiharvest::j, isotrackApplyRegressor::k, DetachedQuadStep_cff::seeds, and EZArrayFL< T >::size().

Referenced by 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 31 of file HFClusterAlgo.h.

References isMC(), and m_isMC.

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

31 { m_isMC = isMC; }

HFClusterAlgo::isMC

void isMC(bool isMC)

Definition: HFClusterAlgo.h:31

HFClusterAlgo::m_isMC

bool m_isMC

Definition: HFClusterAlgo.h:48

bool HFClusterAlgo::isPMTHit ( const HFRecHit & hfr )

private

Definition at line 389 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 193 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, HcalDetId::ieta(), HcalDetId::ietaAbs(), indexByEta(), HcalDetId::iphi(), 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 400 of file HFClusterAlgo.cc.

References cuy::ii.

Referenced by HFEMClusterProducer::beginRun().

                                 {
   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 37 of file HFClusterAlgo.cc.

References cuy::ii, M_PI, and MCMaterialCorrections_3XX.

                                              {
   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 42 of file HFClusterAlgo.h.

friend class CompareHFCore

friend

Definition at line 43 of file HFClusterAlgo.h.

Member Data Documentation

std::vector<double> HFClusterAlgo::m_correctionByEta

private

Definition at line 51 of file HFClusterAlgo.h.

int HFClusterAlgo::m_correctionSet

private

Definition at line 49 of file HFClusterAlgo.h.

std::vector<double> HFClusterAlgo::m_cutByEta

private

Definition at line 50 of file HFClusterAlgo.h.

bool HFClusterAlgo::m_forcePulseFlagMC

private

Definition at line 47 of file HFClusterAlgo.h.

bool HFClusterAlgo::m_isMC

private

Definition at line 48 of file HFClusterAlgo.h.

Referenced by isMC().

double HFClusterAlgo::m_maximumRenergy

private

Definition at line 45 of file HFClusterAlgo.h.

double HFClusterAlgo::m_maximumSL

private

Definition at line 45 of file HFClusterAlgo.h.

double HFClusterAlgo::m_minTowerEnergy

private

Definition at line 45 of file HFClusterAlgo.h.

std::vector<double> HFClusterAlgo::m_seedmnEta

private

Definition at line 52 of file HFClusterAlgo.h.

std::vector<double> HFClusterAlgo::m_seedmnPhi

private

Definition at line 54 of file HFClusterAlgo.h.

std::vector<double> HFClusterAlgo::m_seedMXeta

private

Definition at line 53 of file HFClusterAlgo.h.

std::vector<double> HFClusterAlgo::m_seedMXphi

private

Definition at line 55 of file HFClusterAlgo.h.

double HFClusterAlgo::m_seedThreshold

private

Definition at line 45 of file HFClusterAlgo.h.

bool HFClusterAlgo::m_usePMTFlag

private

Definition at line 46 of file HFClusterAlgo.h.

bool HFClusterAlgo::m_usePulseFlag

private

Definition at line 47 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