22 prescaleFactor_ =
pset.getUntrackedParameter<
int>(
"prescaleFactor", 1);
24 barrelEcalHits_token_ = consumes<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit> > >(
26 endcapEcalHits_token_ = consumes<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit> > >(
32 seleXtalMinEnergy_ =
pset.getParameter<
double>(
"seleXtalMinEnergy");
33 clusSeedThr_ =
pset.getParameter<
double>(
"clusSeedThr");
34 clusEtaSize_ =
pset.getParameter<
int>(
"clusEtaSize");
35 clusPhiSize_ =
pset.getParameter<
int>(
"clusPhiSize");
36 ParameterLogWeighted_ =
pset.getParameter<
bool>(
"ParameterLogWeighted");
37 ParameterX0_ =
pset.getParameter<
double>(
"ParameterX0");
38 ParameterT0_barl_ =
pset.getParameter<
double>(
"ParameterT0_barl");
39 ParameterW0_ =
pset.getParameter<
double>(
"ParameterW0");
41 selePtGammaOne_ =
pset.getParameter<
double>(
"selePtGammaOne");
42 selePtGammaTwo_ =
pset.getParameter<
double>(
"selePtGammaTwo");
43 seleS4S9GammaOne_ =
pset.getParameter<
double>(
"seleS4S9GammaOne");
44 seleS4S9GammaTwo_ =
pset.getParameter<
double>(
"seleS4S9GammaTwo");
45 selePtPi0_ =
pset.getParameter<
double>(
"selePtPi0");
46 selePi0Iso_ =
pset.getParameter<
double>(
"selePi0Iso");
47 selePi0BeltDR_ =
pset.getParameter<
double>(
"selePi0BeltDR");
48 selePi0BeltDeta_ =
pset.getParameter<
double>(
"selePi0BeltDeta");
49 seleMinvMaxPi0_ =
pset.getParameter<
double>(
"seleMinvMaxPi0");
50 seleMinvMinPi0_ =
pset.getParameter<
double>(
"seleMinvMinPi0");
60 currentFolder_.str(
"");
61 currentFolder_ <<
"Egamma/PiZeroAnalyzer/";
64 hMinvPi0EB_ = ibooker.
book1D(
"Pi0InvmassEB",
"Pi0 Invariant Mass in EB", 100, 0., 0.5);
67 hPt1Pi0EB_ = ibooker.
book1D(
"Pt1Pi0EB",
"Pt 1st most energetic Pi0 photon in EB", 100, 0., 20.);
70 hPt2Pi0EB_ = ibooker.
book1D(
"Pt2Pi0EB",
"Pt 2nd most energetic Pi0 photon in EB", 100, 0., 20.);
73 hPtPi0EB_ = ibooker.
book1D(
"PtPi0EB",
"Pi0 Pt in EB", 100, 0., 20.);
76 hIsoPi0EB_ = ibooker.
book1D(
"IsoPi0EB",
"Pi0 Iso in EB", 50, 0., 1.);
83 if (nEvt_ % prescaleFactor_)
86 LogInfo(
"PiZeroAnalyzer") <<
"PiZeroAnalyzer Analyzing event number: " <<
e.id() <<
" Global Counter " << nEvt_
90 bool validEcalRecHits =
true;
93 e.getByToken(barrelEcalHits_token_, barrelHitHandle);
94 if (!barrelHitHandle.
isValid()) {
95 edm::LogError(
"PhotonProducer") <<
"Error! Can't get the product: barrelEcalHits_token_";
96 validEcalRecHits =
false;
100 e.getByToken(endcapEcalHits_token_, endcapHitHandle);
102 if (!endcapHitHandle.
isValid()) {
103 edm::LogError(
"PhotonProducer") <<
"Error! Can't get the product: endcapEcalHits_token_";
104 validEcalRecHits =
false;
107 if (validEcalRecHits)
108 makePizero(esup, barrelHitHandle, endcapHitHandle);
129 std::map<DetId, EcalRecHit> recHitsEB_map;
131 std::vector<EcalRecHit>
seeds;
135 vector<EBDetId> usedXtals;
140 static const int MAXCLUS = 2000;
143 vector<float> etClus;
144 vector<float> etaClus;
145 vector<float> phiClus;
146 vector<EBDetId> max_hit;
147 vector<vector<EcalRecHit> > RecHitsCluster;
148 vector<float> s4s9Clus;
151 for (itb = rhEB->
begin(); itb != rhEB->
end(); ++itb) {
153 double energy = itb->energy();
154 if (
energy > seleXtalMinEnergy_) {
155 std::pair<DetId, EcalRecHit> map_entry(itb->id(), *itb);
156 recHitsEB_map.insert(map_entry);
158 if (
energy > clusSeedThr_)
159 seeds.push_back(*itb);
162 sort(
seeds.begin(),
seeds.end(), [](
auto& x,
auto& y) {
return (x.energy() > y.energy()); });
163 for (std::vector<EcalRecHit>::iterator itseed =
seeds.begin(); itseed !=
seeds.end(); itseed++) {
164 EBDetId seed_id = itseed->id();
165 std::vector<EBDetId>::const_iterator usedIds;
167 bool seedAlreadyUsed =
false;
168 for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
169 if (*usedIds == seed_id) {
170 seedAlreadyUsed =
true;
178 std::vector<DetId> clus_v = topology_p->
getWindow(seed_id, clusEtaSize_, clusPhiSize_);
180 std::vector<std::pair<DetId, float> > clus_used;
182 vector<EcalRecHit> RecHitsInWindow;
184 double simple_energy = 0;
186 for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
189 bool HitAlreadyUsed =
false;
190 for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
191 if (*usedIds == *det) {
192 HitAlreadyUsed =
true;
198 if (recHitsEB_map.find(*det) != recHitsEB_map.end()) {
200 std::map<DetId, EcalRecHit>::iterator aHit;
201 aHit = recHitsEB_map.find(*det);
202 usedXtals.push_back(*det);
203 RecHitsInWindow.push_back(aHit->second);
204 clus_used.push_back(std::pair<DetId, float>(*det, 1.));
205 simple_energy = simple_energy + aHit->second.energy();
210 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
211 float p0x_s = simple_energy *
sin(theta_s) *
cos(clus_pos.phi());
212 float p0y_s = simple_energy *
sin(theta_s) *
sin(clus_pos.phi());
214 float et_s =
sqrt(p0x_s * p0x_s + p0y_s * p0y_s);
218 eClus.push_back(simple_energy);
219 etClus.push_back(et_s);
220 etaClus.push_back(clus_pos.eta());
221 phiClus.push_back(clus_pos.phi());
222 max_hit.push_back(seed_id);
223 RecHitsCluster.push_back(RecHitsInWindow);
227 for (
int i = 0;
i < 4;
i++)
228 s4s9_[
i] = itseed->energy();
229 for (
unsigned int j = 0;
j < RecHitsInWindow.size();
j++) {
235 s4s9_[0] += RecHitsInWindow[
j].energy();
238 s4s9_[0] += RecHitsInWindow[
j].energy();
239 s4s9_[1] += RecHitsInWindow[
j].energy();
243 s4s9_[1] += RecHitsInWindow[
j].energy();
249 if (((
EBDetId)RecHitsInWindow[
j].
id()).iphi() == seed_id.
iphi() - 1 ||
251 s4s9_[0] += RecHitsInWindow[
j].energy();
252 s4s9_[3] += RecHitsInWindow[
j].energy();
256 s4s9_[1] += RecHitsInWindow[
j].energy();
257 s4s9_[2] += RecHitsInWindow[
j].energy();
261 if ((((
EBDetId)RecHitsInWindow[
j].
id()).
ieta() == seed_id.
ieta() + 1 && seed_id.
ieta() != -1) ||
265 s4s9_[3] += RecHitsInWindow[
j].energy();
268 s4s9_[2] += RecHitsInWindow[
j].energy();
269 s4s9_[3] += RecHitsInWindow[
j].energy();
273 s4s9_[2] += RecHitsInWindow[
j].energy();
278 cout <<
" (EBDetId)RecHitsInWindow[j].id()).ieta() " << ((
EBDetId)RecHitsInWindow[
j].
id()).
ieta()
279 <<
" seed_id.ieta() " << seed_id.
ieta() << endl;
280 cout <<
" Problem with S4 calculation " << endl;
286 s4s9Clus.push_back(*max_element(s4s9_, s4s9_ + 4) / simple_energy);
290 if (nClus == MAXCLUS)
298 static const int MAXPI0S = 200;
301 vector<EBDetId> scXtals;
306 for (Int_t
i = 0;
i < nClus;
i++) {
307 for (Int_t
j =
i + 1;
j < nClus;
j++) {
309 if (etClus[
i] > selePtGammaOne_ && etClus[
j] > selePtGammaTwo_ && s4s9Clus[
i] > seleS4S9GammaOne_ &&
310 s4s9Clus[
j] > seleS4S9GammaTwo_) {
311 float theta_0 = 2. * atan(
exp(-etaClus[
i]));
312 float theta_1 = 2. * atan(
exp(-etaClus[
j]));
314 float p0x = eClus[
i] *
sin(theta_0) *
cos(phiClus[
i]);
315 float p1x = eClus[
j] *
sin(theta_1) *
cos(phiClus[
j]);
316 float p0y = eClus[
i] *
sin(theta_0) *
sin(phiClus[
i]);
317 float p1y = eClus[
j] *
sin(theta_1) *
sin(phiClus[
j]);
318 float p0z = eClus[
i] *
cos(theta_0);
319 float p1z = eClus[
j] *
cos(theta_1);
321 float pt_pi0 =
sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
323 if (pt_pi0 < selePtPi0_)
325 float m_inv =
sqrt((eClus[
i] + eClus[
j]) * (eClus[
i] + eClus[
j]) - (p0x + p1x) * (p0x + p1x) -
326 (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
327 if ((m_inv < seleMinvMaxPi0_) && (m_inv > seleMinvMinPi0_)) {
332 TVector3 pi0vect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
333 for (Int_t
k = 0;
k < nClus;
k++) {
334 if (
k ==
i ||
k ==
j)
336 TVector3 Clusvect = TVector3(eClus[
k] *
sin(2. * atan(
exp(-etaClus[
k]))) *
cos(phiClus[
k]),
337 eClus[
k] *
sin(2. * atan(
exp(-etaClus[
k]))) *
sin(phiClus[
k]),
338 eClus[
k] *
cos(2. * atan(
exp(-etaClus[
k]))));
339 float dretaclpi0 = fabs(etaClus[
k] - pi0vect.Eta());
340 float drclpi0 = Clusvect.DeltaR(pi0vect);
342 if ((drclpi0 < selePi0BeltDR_) && (dretaclpi0 < selePi0BeltDeta_)) {
343 Iso = Iso + etClus[
k];
344 IsoClus.push_back(
k);
348 if (Iso / pt_pi0 < selePi0Iso_) {
349 hMinvPi0EB_->Fill(m_inv);
350 hPt1Pi0EB_->Fill(etClus[
i]);
351 hPt2Pi0EB_->Fill(etClus[
j]);
352 hPtPi0EB_->Fill(pt_pi0);
353 hIsoPi0EB_->Fill(Iso / pt_pi0);
358 if (npi0_s == MAXPI0S)