18 : g4InfoLabel(ps.getParameter<
std::
string>(
"moduleLabelG4")),
29 for (
int myStep = 0; myStep < 26; myStep++) {
35 ib.setCurrentFolder(
"EcalHitsV/EcalSimHitsValidation");
36 ib.setScope(MonitorElementData::Scope::RUN);
41 histo =
"EB crystals multiplicity";
44 histo =
"EB occupancy";
47 histo =
"EB longitudinal shower profile";
50 histo =
"EB hits energy spectrum";
53 histo =
"EB hits log10energy spectrum";
56 histo =
"EB hits log10energy spectrum vs normalized energy";
59 histo =
"EB hits log10energy spectrum vs normalized energy25";
62 histo =
"EB hits energy spectrum 2";
65 histo =
"EB crystal energy spectrum";
68 histo =
"EB crystal energy spectrum 2";
104 histo =
"EB E16oE25";
109 edm::LogInfo(
"EventInfo") <<
" Run = " <<
e.id().run() <<
" Event = " <<
e.id().event();
121 std::vector<PCaloHit> theEBCaloHits;
122 theEBCaloHits.insert(theEBCaloHits.end(), EcalHitsEB->begin(), EcalHitsEB->end());
126 double EBEnergy_ = 0.;
127 std::map<unsigned int, std::vector<PCaloHit *>, std::less<unsigned int>> CaloHitMap;
134 std::vector<double> econtr(140, 0.);
135 std::vector<double> econtr25(140, 0.);
138 uint32_t nEBHits = 0;
140 for (std::vector<PCaloHit>::iterator isim = theEBCaloHits.begin(); isim != theEBCaloHits.end(); ++isim) {
141 if (isim->time() > 500.) {
145 CaloHitMap[isim->id()].push_back(&(*isim));
149 LogDebug(
"HitInfo") <<
" CaloHit " << isim->getName() <<
"\n"
150 <<
" DetID = " << isim->id() <<
" EBDetId = " << ebid.ieta() <<
" " << ebid.iphi() <<
"\n"
151 <<
" Time = " << isim->time() <<
"\n"
152 <<
" Track Id = " << isim->geantTrackId() <<
"\n"
153 <<
" Energy = " << isim->energy();
157 uint32_t crystid = ebid.rawId();
158 ebmap[crystid] += isim->energy();
160 EBEnergy_ += isim->energy();
163 if (isim->energy() > 0) {
165 int log10i =
int((log10(isim->energy()) + 10.) * 10.);
166 if (log10i >= 0 && log10i < 140)
167 econtr[log10i] += isim->energy();
174 for (
std::map<uint32_t,
float, std::less<uint32_t>>::iterator it = ebmap.begin(); it != ebmap.end(); ++it)
178 for (
std::map<uint32_t,
float, std::less<uint32_t>>::iterator it = ebmap.begin(); it != ebmap.end(); ++it)
188 int by = myEBid.
iphi();
189 int bz = myEBid.
zside();
197 std::vector<uint32_t> ids25;
200 for (
unsigned i = 0;
i < 25;
i++) {
201 for (
unsigned int j = 0;
j < CaloHitMap[ids25[
i]].size();
j++) {
202 if (CaloHitMap[ids25[
i]][
j]->
energy() > 0) {
203 int log10i =
int((log10(CaloHitMap[ids25[
i]][
j]->
energy()) + 10.) * 10.);
204 if (log10i >= 0 && log10i < 140)
205 econtr25[log10i] += CaloHitMap[ids25[
i]][
j]->energy();
235 if (EBEnergy_ != 0) {
236 for (
int i = 0;
i < 140;
i++) {
242 for (
int i = 0;
i < 140;
i++) {
248 if (MyPEcalValidInfo.
isValid()) {
249 if (MyPEcalValidInfo->
eb1x1() > 0.) {
250 std::vector<float> BX0 = MyPEcalValidInfo->
bX0();
252 for (
int myStep = 0; myStep < 26; myStep++) {
261 int nCellInEta,
int nCellInPhi,
int centralEta,
int centralPhi,
int centralZ,
MapType &themap) {
263 float totalEnergy = 0.;
265 int goBackInEta = nCellInEta / 2;
266 int goBackInPhi = nCellInPhi / 2;
267 int startEta = centralZ * centralEta - goBackInEta;
268 int startPhi = centralPhi - goBackInPhi;
270 for (
int ieta = startEta;
ieta < startEta + nCellInEta;
ieta++) {
271 for (
int iphi = startPhi;
iphi < startPhi + nCellInPhi;
iphi++) {
278 }
else if (
iphi > 360) {
284 totalEnergy += themap[
index];
289 LogDebug(
"GeomInfo") << nCellInEta <<
" x " << nCellInPhi <<
" EB matrix energy = " << totalEnergy <<
" for "
290 << ncristals <<
" crystals";
295 int nCellInEta,
int nCellInPhi,
int centralEta,
int centralPhi,
int centralZ,
MapType &themap) {
297 std::vector<uint32_t> ids(nCellInEta * nCellInPhi);
299 int goBackInEta = nCellInEta / 2;
300 int goBackInPhi = nCellInPhi / 2;
301 int startEta = centralZ * centralEta - goBackInEta;
302 int startPhi = centralPhi - goBackInPhi;
304 for (
int ieta = startEta;
ieta < startEta + nCellInEta;
ieta++) {
305 for (
int iphi = startPhi;
iphi < startPhi + nCellInPhi;
iphi++) {
312 }
else if (
iphi > 360) {
317 ids[ncristals] =
index;
326 int nCellInEta,
int nCellInPhi,
int CentralEta,
int CentralPhi,
int CentralZ,
MapType &fillmap,
MapType &themap) {
327 int goBackInEta = nCellInEta / 2;
328 int goBackInPhi = nCellInPhi / 2;
330 int startEta = CentralZ * CentralEta - goBackInEta;
331 int startPhi = CentralPhi - goBackInPhi;
334 for (
int ieta = startEta;
ieta < startEta + nCellInEta;
ieta++) {
335 for (
int iphi = startPhi;
iphi < startPhi + nCellInPhi;
iphi++) {
342 }
else if (
iphi > 360) {
347 fillmap[
i++] = themap[
index];
353 if (fillmap[
i / 2] == themap[ebcenterid])
361 float e012 = themap[0] + themap[1] + themap[2];
362 float e036 = themap[0] + themap[3] + themap[6];
363 float e678 = themap[6] + themap[7] + themap[8];
364 float e258 = themap[2] + themap[5] + themap[8];
366 if ((e012 > e678 || e012 == e678) && (e036 > e258 || e036 == e258))
367 return E22 = themap[0] + themap[1] + themap[3] + themap[4];
368 else if ((e012 > e678 || e012 == e678) && (e036 < e258 || e036 == e258))
369 return E22 = themap[1] + themap[2] + themap[4] + themap[5];
370 else if ((e012 < e678 || e012 == e678) && (e036 > e258 || e036 == e258))
371 return E22 = themap[3] + themap[4] + themap[6] + themap[7];
372 else if ((e012 < e678 || e012 == e678) && (e036 < e258 || e036 == e258))
373 return E22 = themap[4] + themap[5] + themap[7] + themap[8];
381 float e0_4 = themap[0] + themap[1] + themap[2] + themap[3] + themap[4];
382 float e0_20 = themap[0] + themap[5] + themap[10] + themap[15] + themap[20];
383 float e4_24 = themap[4] + themap[9] + themap[14] + themap[19] + themap[24];
384 float e0_24 = themap[20] + themap[21] + themap[22] + themap[23] + themap[24];
386 if ((e0_4 > e0_24 || e0_4 == e0_24) && (e0_20 > e4_24 || e0_20 == e4_24))
387 return E44 = e33 + themap[0] + themap[1] + themap[2] + themap[3] + themap[5] + themap[10] + themap[15];
388 else if ((e0_4 > e0_24 || e0_4 == e0_24) && (e0_20 < e4_24 || e0_20 == e4_24))
389 return E44 = e33 + themap[1] + themap[2] + themap[3] + themap[4] + themap[9] + themap[14] + themap[19];
390 else if ((e0_4 < e0_24 || e0_4 == e0_24) && (e0_20 > e4_24 || e0_20 == e4_24))
391 return E44 = e33 + themap[5] + themap[10] + themap[15] + themap[20] + themap[21] + themap[22] + themap[23];
392 else if ((e0_4 < e0_24 || e0_4 == e0_24) && (e0_20 < e4_24 || e0_20 == e4_24))
393 return E44 = e33 + themap[21] + themap[22] + themap[23] + themap[24] + themap[9] + themap[14] + themap[19];
401 uint32_t unitWithMaxEnergy = 0;
404 MapType::iterator iter;
405 for (iter = themap.begin(); iter != themap.end(); iter++) {
408 unitWithMaxEnergy = (*iter).first;
412 LogDebug(
"GeomInfo") <<
" max energy of " <<
maxEnergy <<
" GeV in Unit id " << unitWithMaxEnergy;
413 return unitWithMaxEnergy;