132 std::map<DetId, EcalRecHit> recHitsEB_map;
134 std::vector<EcalRecHit>
seeds;
138 vector<EBDetId> usedXtals;
143 static const int MAXCLUS = 2000;
146 vector<float> etClus;
147 vector<float> etaClus;
148 vector<float> phiClus;
149 vector<EBDetId> max_hit;
150 vector<vector<EcalRecHit> > RecHitsCluster;
151 vector<float> s4s9Clus;
154 for (itb = rhEB->begin(); itb != rhEB->end(); ++itb) {
156 double energy = itb->energy();
158 std::pair<DetId, EcalRecHit> map_entry(itb->id(), *itb);
159 recHitsEB_map.insert(map_entry);
162 seeds.push_back(*itb);
165 sort(seeds.begin(), seeds.end(), [](
auto&
x,
auto&
y) {
return (
x.energy() >
y.energy()); });
166 for (std::vector<EcalRecHit>::iterator itseed = seeds.begin(); itseed != seeds.end(); itseed++) {
167 EBDetId seed_id = itseed->id();
168 std::vector<EBDetId>::const_iterator usedIds;
170 bool seedAlreadyUsed =
false;
171 for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
172 if (*usedIds == seed_id) {
173 seedAlreadyUsed =
true;
183 std::vector<std::pair<DetId, float> > clus_used;
185 vector<EcalRecHit> RecHitsInWindow;
187 double simple_energy = 0;
189 for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
192 bool HitAlreadyUsed =
false;
193 for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
194 if (*usedIds == *det) {
195 HitAlreadyUsed =
true;
201 if (recHitsEB_map.find(*det) != recHitsEB_map.end()) {
203 std::map<DetId, EcalRecHit>::iterator aHit;
204 aHit = recHitsEB_map.find(*det);
205 usedXtals.push_back(*det);
206 RecHitsInWindow.push_back(aHit->second);
207 clus_used.push_back(std::pair<DetId, float>(*det, 1.));
208 simple_energy = simple_energy + aHit->second.energy();
213 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
214 float p0x_s = simple_energy *
sin(theta_s) *
cos(clus_pos.phi());
215 float p0y_s = simple_energy *
sin(theta_s) *
sin(clus_pos.phi());
217 float et_s =
sqrt(p0x_s * p0x_s + p0y_s * p0y_s);
221 eClus.push_back(simple_energy);
222 etClus.push_back(et_s);
223 etaClus.push_back(clus_pos.eta());
224 phiClus.push_back(clus_pos.phi());
225 max_hit.push_back(seed_id);
226 RecHitsCluster.push_back(RecHitsInWindow);
230 for (
int i = 0;
i < 4;
i++)
231 s4s9_[
i] = itseed->energy();
232 for (
unsigned int j = 0;
j < RecHitsInWindow.size();
j++) {
234 if ((((
EBDetId)RecHitsInWindow[
j].
id()).ieta() == seed_id.
ieta() - 1 && seed_id.
ieta() != 1) ||
235 (seed_id.
ieta() == 1 && (((
EBDetId)RecHitsInWindow[
j].
id()).ieta() == seed_id.
ieta() - 2))) {
236 if (((
EBDetId)RecHitsInWindow[
j].id()).iphi() == seed_id.
iphi() - 1 ||
237 ((
EBDetId)RecHitsInWindow[
j].
id()).iphi() - 360 == seed_id.
iphi() - 1) {
238 s4s9_[0] += RecHitsInWindow[
j].energy();
240 if (((
EBDetId)RecHitsInWindow[
j].
id()).iphi() == seed_id.
iphi()) {
241 s4s9_[0] += RecHitsInWindow[
j].energy();
242 s4s9_[1] += RecHitsInWindow[
j].energy();
244 if (((
EBDetId)RecHitsInWindow[
j].
id()).iphi() == seed_id.
iphi() + 1 ||
245 ((
EBDetId)RecHitsInWindow[
j].
id()).iphi() - 360 == seed_id.
iphi() + 1) {
246 s4s9_[1] += RecHitsInWindow[
j].energy();
251 if (((
EBDetId)RecHitsInWindow[
j].
id()).ieta() == seed_id.
ieta()) {
252 if (((
EBDetId)RecHitsInWindow[
j].
id()).iphi() == seed_id.
iphi() - 1 ||
253 ((
EBDetId)RecHitsInWindow[
j].
id()).iphi() - 360 == seed_id.
iphi() - 1) {
254 s4s9_[0] += RecHitsInWindow[
j].energy();
255 s4s9_[3] += RecHitsInWindow[
j].energy();
257 if (((
EBDetId)RecHitsInWindow[
j].
id()).iphi() == seed_id.
iphi() + 1 ||
258 ((
EBDetId)RecHitsInWindow[
j].
id()).iphi() - 360 == seed_id.
iphi() + 1) {
259 s4s9_[1] += RecHitsInWindow[
j].energy();
260 s4s9_[2] += RecHitsInWindow[
j].energy();
264 if ((((
EBDetId)RecHitsInWindow[
j].
id()).ieta() == seed_id.
ieta() + 1 && seed_id.
ieta() != -1) ||
265 (seed_id.
ieta() == -1 && (((
EBDetId)RecHitsInWindow[
j].
id()).ieta() == seed_id.
ieta() + 2))) {
266 if (((
EBDetId)RecHitsInWindow[
j].id()).iphi() == seed_id.
iphi() - 1 ||
267 ((
EBDetId)RecHitsInWindow[
j].
id()).iphi() - 360 == seed_id.
iphi() - 1) {
268 s4s9_[3] += RecHitsInWindow[
j].energy();
270 if (((
EBDetId)RecHitsInWindow[
j].
id()).iphi() == seed_id.
iphi()) {
271 s4s9_[2] += RecHitsInWindow[
j].energy();
272 s4s9_[3] += RecHitsInWindow[
j].energy();
274 if (((
EBDetId)RecHitsInWindow[
j].
id()).iphi() == seed_id.
iphi() + 1 ||
275 ((
EBDetId)RecHitsInWindow[
j].
id()).iphi() - 360 == seed_id.
iphi() + 1) {
276 s4s9_[2] += RecHitsInWindow[
j].energy();
281 cout <<
" (EBDetId)RecHitsInWindow[j].id()).ieta() " << ((
EBDetId)RecHitsInWindow[
j].
id()).ieta()
282 <<
" seed_id.ieta() " << seed_id.
ieta() << endl;
283 cout <<
" Problem with S4 calculation " << endl;
289 s4s9Clus.push_back(*max_element(s4s9_, s4s9_ + 4) / simple_energy);
293 if (nClus == MAXCLUS)
301 static const int MAXPI0S = 200;
304 vector<EBDetId> scXtals;
309 for (Int_t
i = 0;
i < nClus;
i++) {
310 for (Int_t
j =
i + 1;
j < nClus;
j++) {
314 float theta_0 = 2. * atan(
exp(-etaClus[
i]));
315 float theta_1 = 2. * atan(
exp(-etaClus[
j]));
317 float p0x = eClus[
i] *
sin(theta_0) *
cos(phiClus[
i]);
318 float p1x = eClus[
j] *
sin(theta_1) *
cos(phiClus[
j]);
319 float p0y = eClus[
i] *
sin(theta_0) *
sin(phiClus[
i]);
320 float p1y = eClus[
j] *
sin(theta_1) *
sin(phiClus[
j]);
321 float p0z = eClus[
i] *
cos(theta_0);
322 float p1z = eClus[
j] *
cos(theta_1);
324 float pt_pi0 =
sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
328 float m_inv =
sqrt((eClus[
i] + eClus[
j]) * (eClus[
i] + eClus[j]) - (p0x + p1x) * (p0x + p1x) -
329 (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
335 TVector3 pi0vect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
336 for (Int_t
k = 0;
k < nClus;
k++) {
337 if (
k ==
i ||
k == j)
339 TVector3 Clusvect = TVector3(eClus[
k] *
sin(2. * atan(
exp(-etaClus[
k]))) *
cos(phiClus[k]),
340 eClus[k] *
sin(2. * atan(
exp(-etaClus[k]))) *
sin(phiClus[k]),
341 eClus[k] *
cos(2. * atan(
exp(-etaClus[k]))));
342 float dretaclpi0 = fabs(etaClus[k] - pi0vect.Eta());
343 float drclpi0 = Clusvect.DeltaR(pi0vect);
346 Iso = Iso + etClus[
k];
347 IsoClus.push_back(k);
361 if (npi0_s == MAXPI0S)
double clusSeedThr_
parameters needed for pizero finding
uint16_t *__restrict__ id
MonitorElement * hIsoPi0EB_
Sin< T >::type sin(const T &t)
MonitorElement * hPtPi0EB_
MonitorElement * hMinvPi0EB_
std::vector< EcalRecHit >::const_iterator const_iterator
edm::ParameterSet posCalcParameters_
Exp< T >::type exp(const T &t)
int iphi() const
get the crystal iphi
const edm::ESGetToken< CaloTopology, CaloTopologyRecord > caloTopologyToken_
const edm::ESGetToken< CaloGeometry, CaloGeometryRecord > caloGeometryToken_
Cos< T >::type cos(const T &t)
math::XYZPoint Calculate_Location(const HitsAndFractions &iDetIds, const edm::SortedCollection< HitType > *iRecHits, const CaloSubdetectorGeometry *iSubGeom, const CaloSubdetectorGeometry *iESGeom=nullptr)
int ieta() const
get the crystal ieta
virtual std::vector< DetId > getWindow(const DetId &id, const int &northSouthSize, const int &eastWestSize) const
XYZPointD XYZPoint
point in space with cartesian internal representation
MonitorElement * hPt2Pi0EB_
MonitorElement * hPt1Pi0EB_
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
double seleXtalMinEnergy_