|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
|
#include "Calibration/IsolatedParticles/interface/CaloConstants.h"
#include "Calibration/IsolatedParticles/interface/FindDistCone.h"
#include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
namespace spr {
// Cone clustering core
double getDistInPlaneTrackDir(const GlobalPoint& caloPoint,
const GlobalVector& caloVector,
const GlobalPoint& rechitPoint,
bool debug) {
const GlobalVector caloIntersectVector(caloPoint.x(), caloPoint.y(),
caloPoint.z()); //p
const GlobalVector caloUnitVector = caloVector.unit();
const GlobalVector rechitVector(rechitPoint.x(), rechitPoint.y(), rechitPoint.z());
const GlobalVector rechitUnitVector = rechitVector.unit();
double dotprod_denominator = caloUnitVector.dot(rechitUnitVector);
double dotprod_numerator = caloUnitVector.dot(caloIntersectVector);
double rechitdist = dotprod_numerator / dotprod_denominator;
const GlobalVector effectiveRechitVector = rechitdist * rechitUnitVector;
const GlobalPoint effectiveRechitPoint(
effectiveRechitVector.x(), effectiveRechitVector.y(), effectiveRechitVector.z());
GlobalVector distance_vector = effectiveRechitPoint - caloPoint;
if (debug) {
edm::LogVerbatim("IsoTrack") << "getDistInPlaneTrackDir: point " << caloPoint << " dirn " << caloVector
<< " numerator " << dotprod_numerator << " denominator " << dotprod_denominator
<< " distance " << distance_vector.mag();
}
if (dotprod_denominator > 0. && dotprod_numerator > 0.) {
return distance_vector.mag();
} else {
return 999999.;
}
}
// Not used, but here for reference
double getDistInCMatEcal(double eta1, double phi1, double eta2, double phi2, bool debug) {
double dR, Rec;
if (fabs(eta1) < spr::etaBEEcal)
Rec = spr::rFrontEB;
else
Rec = spr::zFrontEE;
double ce1 = cosh(eta1);
double ce2 = cosh(eta2);
double te1 = tanh(eta1);
double te2 = tanh(eta2);
double z = cos(phi1 - phi2) / ce1 / ce2 + te1 * te2;
if (z != 0)
dR = fabs(Rec * ce1 * sqrt(1. / z / z - 1.));
else
dR = 999999.;
if (debug)
edm::LogVerbatim("IsoTrack") << "getDistInCMatEcal: between (" << eta1 << ", " << phi1 << ") and (" << eta2
<< ", " << phi2 << " is " << dR;
return dR;
}
// Not used, but here for reference
double getDistInCMatHcal(double eta1, double phi1, double eta2, double phi2, bool debug) {
// Radii and eta from Geometry/HcalCommonData/data/hcalendcapalgo.xml
// and Geometry/HcalCommonData/data/hcalbarrelalgo.xml
double dR, Rec;
if (fabs(eta1) < spr::etaBEHcal)
Rec = spr::rFrontHB;
else
Rec = spr::zFrontHE;
double ce1 = cosh(eta1);
double ce2 = cosh(eta2);
double te1 = tanh(eta1);
double te2 = tanh(eta2);
double z = cos(phi1 - phi2) / ce1 / ce2 + te1 * te2;
if (z != 0)
dR = fabs(Rec * ce1 * sqrt(1. / z / z - 1.));
else
dR = 999999.;
if (debug)
edm::LogVerbatim("IsoTrack") << "getDistInCMatHcal: between (" << eta1 << ", " << phi1 << ") and (" << eta2
<< ", " << phi2 << " is " << dR;
return dR;
}
void getEtaPhi(HBHERecHitCollection::const_iterator hit,
std::vector<int>& RH_ieta,
std::vector<int>& RH_iphi,
std::vector<double>& RH_ene,
bool) {
RH_ieta.push_back(hit->id().ieta());
RH_iphi.push_back(hit->id().iphi());
RH_ene.push_back(hit->energy());
}
void getEtaPhi(edm::PCaloHitContainer::const_iterator hit,
std::vector<int>& RH_ieta,
std::vector<int>& RH_iphi,
std::vector<double>& RH_ene,
bool) {
// SimHit function not yet implemented.
RH_ieta.push_back(-9);
RH_iphi.push_back(-9);
RH_ene.push_back(-9.);
}
void getEtaPhi(EcalRecHitCollection::const_iterator hit,
std::vector<int>& RH_ieta,
std::vector<int>& RH_iphi,
std::vector<double>& RH_ene,
bool) {
// Ecal function not yet implemented.
RH_ieta.push_back(-9);
RH_iphi.push_back(-9);
RH_ene.push_back(-9.);
}
void getEtaPhi(HBHERecHitCollection::const_iterator hit, int& ieta, int& iphi, bool) {
ieta = hit->id().ieta();
iphi = hit->id().iphi();
}
void getEtaPhi(edm::PCaloHitContainer::const_iterator hit, int& ieta, int& iphi, bool) {
DetId id = DetId(hit->id());
if (id.det() == DetId::Hcal) {
ieta = ((HcalDetId)(hit->id())).ieta();
iphi = ((HcalDetId)(hit->id())).iphi();
} else if (id.det() == DetId::Ecal && id.subdetId() == EcalBarrel) {
ieta = ((EBDetId)(id)).ieta();
iphi = ((EBDetId)(id)).iphi();
} else {
ieta = 999;
iphi = 999;
}
}
void getEtaPhi(EcalRecHitCollection::const_iterator hit, int& ieta, int& iphi, bool) {
DetId id = hit->id();
if (id.subdetId() == EcalBarrel) {
ieta = ((EBDetId)(id)).ieta();
iphi = ((EBDetId)(id)).iphi();
} else {
ieta = 999;
iphi = 999;
}
}
double getEnergy(HBHERecHitCollection::const_iterator hit, int useRaw, bool) {
double energy = ((useRaw == 1) ? hit->eraw() : ((useRaw == 2) ? hit->eaux() : hit->energy()));
return energy;
}
double getEnergy(EcalRecHitCollection::const_iterator hit, int, bool) { return hit->energy(); }
double getEnergy(edm::PCaloHitContainer::const_iterator hit, int, bool) {
// This will not yet handle Ecal CaloHits!!
double samplingWeight = 1.;
// Hard coded sampling weights from JFH analysis of iso tracks
// Sept 2009.
DetId id = hit->id();
if (id.det() == DetId::Hcal) {
HcalDetId detId(id);
if (detId.subdet() == HcalBarrel)
samplingWeight = 114.1;
else if (detId.subdet() == HcalEndcap)
samplingWeight = 167.3;
else {
// ONLY protection against summing HO, HF simhits
return 0.;
}
}
return samplingWeight * hit->energy();
}
GlobalPoint getGpos(const CaloGeometry* geo, HBHERecHitCollection::const_iterator hit, bool) {
DetId detId(hit->id());
return (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(detId)))->getPosition(detId);
}
GlobalPoint getGpos(const CaloGeometry* geo, edm::PCaloHitContainer::const_iterator hit, bool) {
DetId detId(hit->id());
GlobalPoint point = (detId.det() == DetId::Hcal)
? (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(detId)))->getPosition(detId)
: geo->getPosition(detId);
return point;
}
GlobalPoint getGpos(const CaloGeometry* geo, EcalRecHitCollection::const_iterator hit, bool) {
// Not tested for EcalRecHits!!
if (hit->id().subdetId() == EcalEndcap) {
EEDetId EEid = EEDetId(hit->id());
return geo->getPosition(EEid);
} else { // (hit->id().subdetId() == EcalBarrel)
EBDetId EBid = EBDetId(hit->id());
return geo->getPosition(EBid);
}
}
double getRawEnergy(HBHERecHitCollection::const_iterator hit, int useRaw) {
double energy = ((useRaw == 1) ? hit->eraw() : ((useRaw == 2) ? hit->eaux() : hit->energy()));
return energy;
}
double getRawEnergy(EcalRecHitCollection::const_iterator hit, int) { return hit->energy(); }
double getRawEnergy(edm::PCaloHitContainer::const_iterator hit, int) { return hit->energy(); }
} // namespace spr
|