|
|
|
#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 |
$Id:version 1.2
Definition at line 22 of file HFClusterAlgo.h.
| HFClusterAlgo::HFClusterAlgo | ( | ) |
Definition at line 19 of file HFClusterAlgo.cc.
{
m_isMC=true; // safest
}
| 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().
| 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(), eta(), PV3DBase< T, PVType, FrameType >::eta(), edm::SortedCollection< T, SORT >::find(), CaloGeometry::getPosition(), HcalForward, i, HcalDetId::ieta(), HcalDetId::ietaAbs(), indexByEta(), HcalDetId::iphi(), funct::log(), M_PI, max(), AlCaHLTBitMon_ParallelJobs::p, PV3DBase< T, PVType, FrameType >::phi(), 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];
}
}
friend class CompareHFCompleteHitET [friend] |
Definition at line 38 of file HFClusterAlgo.h.
friend class CompareHFCore [friend] |
Definition at line 39 of file HFClusterAlgo.h.
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.
1.7.3