CMS 3D CMS Logo

ClusterTools.cc
Go to the documentation of this file.
2 
7 
9 
13 
14 #include "vdt/vdtMath.h"
15 
16 #include <iostream>
17 
18 using namespace hgcal;
20 }
21 
23  edm::ConsumesCollector& sumes) :
24  eetok( sumes.consumes<HGCRecHitCollection>(conf.getParameter<edm::InputTag>("HGCEEInput")) ),
25  fhtok( sumes.consumes<HGCRecHitCollection>(conf.getParameter<edm::InputTag>("HGCFHInput")) ),
26  bhtok( sumes.consumes<HGCRecHitCollection>(conf.getParameter<edm::InputTag>("HGCBHInput")) ) {
27  }
28 
30  rhtools_.getEvent(ev);
32  ev.getByToken(eetok, temp);
33  eerh_ = temp.product();
34  ev.getByToken(fhtok, temp);
35  fhrh_ = temp.product();
36  ev.getByToken(bhtok, temp);
37  bhrh_ = temp.product();
38 }
39 
42 }
43 
45  float energy=0.f, energyHad=0.f;
46  const auto& hits = clus.hitsAndFractions();
47  for( const auto& hit : hits ) {
48  const auto& id = hit.first;
49  const float fraction = hit.second;
50  if( id.det() == DetId::Forward ) {
51  switch( id.subdetId() ) {
52  case HGCEE:
53  energy += eerh_->find(id)->energy()*fraction;
54  break;
55  case HGCHEF:
56  {
57  const float temp = fhrh_->find(id)->energy();
58  energy += temp*fraction;
59  energyHad += temp*fraction;
60  }
61  break;
62  default:
63  throw cms::Exception("HGCalClusterTools")
64  << " Cluster contains hits that are not from HGCal! " << std::endl;
65  }
66  } else if ( id.det() == DetId::Hcal && id.subdetId() == HcalEndcap ) {
67  const float temp = bhrh_->find(id)->energy();
68  energy += temp*fraction;
69  energyHad += temp*fraction;
70  } else {
71  throw cms::Exception("HGCalClusterTools")
72  << " Cluster contains hits that are not from HGCal! " << std::endl;
73  }
74  }
75  float fraction = -1.f;
76  if( energy > 0.f ) {
77  fraction = energyHad/energy;
78  }
79  return fraction;
80 }
81 
82 
84  if( clu.clusters().size() == 0 ) return math::XYZPoint();
85 
86  double acc_x = 0.0;
87  double acc_y = 0.0;
88  double acc_z = 0.0;
89  double totweight = 0.;
90  double mcenergy = getMultiClusterEnergy(clu);
91  for( const auto& ptr : clu.clusters() ) {
92  if (mcenergy != 0) {
93  if(ptr->energy()<.01*mcenergy) continue; //cutoff < 1% layer contribution
94  }
95  const double weight = ptr->energy(); // weigth each corrdinate only by the total energy of the layer cluster
96  acc_x += ptr->x()*weight;
97  acc_y += ptr->y()*weight;
98  acc_z += ptr->z()*weight;
99  totweight += weight;
100  }
101  if (totweight != 0) {
102  acc_x /= totweight;
103  acc_y /= totweight;
104  acc_z /= totweight;
105  }
106  // return x/y/z in absolute coordinates
107  return math::XYZPoint(acc_x,acc_y,acc_z);
108 }
109 
110 int ClusterTools::getLayer(const DetId detid) const {
111  return rhtools_.getLayerWithOffset(detid);
112 }
113 
115  double acc = 0.0;
116  for(const auto& ptr : clu.clusters() ) {
117  acc += ptr->energy();
118  }
119  return acc;
120 }
121 
122 bool ClusterTools::getWidths(const reco::CaloCluster & clus,double & sigmaetaeta, double & sigmaphiphi, double & sigmaetaetal, double & sigmaphiphil ) const {
123 
124  if (getLayer(clus.hitsAndFractions()[0].first) > lastLayerEE) return false;
125  const math::XYZPoint & position(clus.position());
126  unsigned nhit=clus.hitsAndFractions().size();
127 
128  sigmaetaeta=0.;
129  sigmaphiphi=0.;
130  sigmaetaetal=0.;
131  sigmaphiphil=0.;
132 
133  double sumw=0.;
134  double sumlogw=0.;
135 
136  for (unsigned int ih=0; ih<nhit; ++ih) {
137  const DetId & id = (clus.hitsAndFractions())[ih].first;
138  if ((clus.hitsAndFractions())[ih].second==0.) continue;
139 
141  if (id.det()==DetId::Forward && id.subdetId()==HGCEE) {
142  const HGCRecHit * theHit = &(*eerh_->find(id));
143 
145  double weight = theHit->energy();
146  // take w0=2 To be optimized
147  double logweight = 0;
148  if (clus.energy() != 0) {
149  logweight = std::max(0.,2 + log(theHit->energy()/clus.energy()));
150  }
151  double deltaetaeta2 = (cellPos.eta()-position.eta())*(cellPos.eta()-position.eta());
152  double deltaphiphi2 = (cellPos.phi()-position.phi())*(cellPos.phi()-position.phi());
153  sigmaetaeta += deltaetaeta2* weight;
154  sigmaphiphi += deltaphiphi2 * weight;
155  sigmaetaetal += deltaetaeta2* logweight;
156  sigmaphiphil += deltaphiphi2 * logweight;
157  sumw += weight;
158  sumlogw += logweight;
159  }
160  }
161 
162 
163  if (sumw<=0.) return false;
164 
165  sigmaetaeta /= sumw;
166  sigmaetaeta = std::sqrt(sigmaetaeta);
167  sigmaphiphi /= sumw;
168  sigmaphiphi = std::sqrt(sigmaphiphi);
169 
170  if (sumlogw != 0) {
171  sigmaetaetal /= sumlogw;
172  sigmaetaetal = std::sqrt(sigmaetaetal);
173  sigmaphiphil /= sumlogw;
174  sigmaphiphil = std::sqrt(sigmaphiphil);
175  }
176 
177  return true;
178 }
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:129
void getEvent(const edm::Event &)
Definition: RecHitTools.cc:63
void getEvent(const edm::Event &)
Definition: ClusterTools.cc:29
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:460
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
std::vector< HGCRecHit >::const_iterator const_iterator
Definition: weight.py:1
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
Definition: CaloCluster.h:195
bool ev
math::XYZPoint getMultiClusterPosition(const reco::HGCalMultiCluster &) const
Definition: ClusterTools.cc:83
void getEventSetup(const edm::EventSetup &)
Definition: RecHitTools.cc:66
U second(std::pair< T, U > const &p)
RecHitTools rhtools_
Definition: ClusterTools.h:56
const edm::PtrVector< reco::BasicCluster > & clusters() const
float energy() const
Definition: CaloRecHit.h:17
T sqrt(T t)
Definition: SSEVec.h:18
double f[11][100]
double energy() const
cluster energy
Definition: CaloCluster.h:124
int getLayer(const DetId) const
const HGCRecHitCollection * eerh_
Definition: ClusterTools.h:58
Definition: DetId.h:18
double getMultiClusterEnergy(const reco::HGCalMultiCluster &) const
unsigned int getLayerWithOffset(const DetId &) const
Definition: RecHitTools.cc:139
T const * product() const
Definition: Handle.h:81
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
GlobalPoint getPosition(const DetId &id) const
Definition: RecHitTools.cc:73
float getClusterHadronFraction(const reco::CaloCluster &) const
Definition: ClusterTools.cc:44
const edm::EDGetTokenT< HGCRecHitCollection > fhtok
Definition: ClusterTools.h:57
T eta() const
Definition: PV3DBase.h:76
const edm::EDGetTokenT< HGCRecHitCollection > eetok
Definition: ClusterTools.h:57
iterator find(key_type k)
const HGCRecHitCollection * bhrh_
Definition: ClusterTools.h:58
HLT enums.
static const int lastLayerEE
Definition: ClusterTools.h:59
static int position[264][3]
Definition: ReadPGInfo.cc:509
const edm::EDGetTokenT< HGCRecHitCollection > bhtok
Definition: ClusterTools.h:57
bool getWidths(const reco::CaloCluster &clus, double &sigmaetaeta, double &sigmaphiphi, double &sigmaetaetalog, double &sigmaphiphilog) const
const HGCRecHitCollection * fhrh_
Definition: ClusterTools.h:58
void getEventSetup(const edm::EventSetup &)
Definition: ClusterTools.cc:40