File indexing completed on 2024-04-06 12:32:06
0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0003 #include "Validation/EcalHits/interface/EcalSimHitsValidProducer.h"
0004
0005 #include "DataFormats/Math/interface/Point3D.h"
0006 #include "SimDataFormats/ValidationFormats/interface/PValidationFormats.h"
0007 #include "SimG4CMS/Calo/interface/CaloG4HitCollection.h"
0008 #include "SimG4Core/Notification/interface/BeginOfEvent.h"
0009 #include "SimG4Core/Notification/interface/BeginOfTrack.h"
0010 #include "SimG4Core/Notification/interface/EndOfEvent.h"
0011
0012 #include "G4PrimaryParticle.hh"
0013 #include "G4PrimaryVertex.hh"
0014 #include "G4SDManager.hh"
0015 #include "G4Step.hh"
0016
0017 #include <iostream>
0018
0019 EcalSimHitsValidProducer::EcalSimHitsValidProducer(const edm::ParameterSet &iPSet)
0020 : ee1(0.0),
0021 ee4(0.0),
0022 ee9(0.0),
0023 ee16(0.0),
0024 ee25(0.0),
0025 eb1(0.0),
0026 eb4(0.0),
0027 eb9(0.0),
0028 eb16(0.0),
0029 eb25(0.0),
0030 totalEInEE(0.0),
0031 totalEInEB(0),
0032 totalEInES(0.0),
0033 totalEInEEzp(0.0),
0034 totalEInEEzm(0.0),
0035 totalEInESzp(0.0),
0036 totalEInESzm(0.0),
0037 totalHits(0),
0038 nHitsInEE(0),
0039 nHitsInEB(0),
0040 nHitsInES(0),
0041 nHitsIn1ES(0),
0042 nHitsIn2ES(0),
0043 nCrystalInEB(0),
0044 nCrystalInEEzp(0),
0045 nCrystalInEEzm(0),
0046 nHitsIn1ESzp(0),
0047 nHitsIn1ESzm(0),
0048 nHitsIn2ESzp(0),
0049 nHitsIn2ESzm(0),
0050 thePID(0),
0051 label(iPSet.getUntrackedParameter<std::string>("instanceLabel", "EcalValidInfo")) {
0052 produces<PEcalValidInfo>(label);
0053
0054 for (int i = 0; i < 26; i++) {
0055 eBX0[i] = 0.0;
0056 eEX0[i] = 0.0;
0057 }
0058 setMT(true);
0059 }
0060
0061 EcalSimHitsValidProducer::~EcalSimHitsValidProducer() {}
0062
0063 void EcalSimHitsValidProducer::produce(edm::Event &e, const edm::EventSetup &) {
0064 std::unique_ptr<PEcalValidInfo> product(new PEcalValidInfo);
0065 fillEventInfo(*product);
0066 e.put(std::move(product), label);
0067 }
0068
0069 void EcalSimHitsValidProducer::fillEventInfo(PEcalValidInfo &product) {
0070 if (ee1 != 0) {
0071 product.ee1 = ee1;
0072 product.ee4 = ee4;
0073 product.ee9 = ee9;
0074 product.ee16 = ee16;
0075 product.ee25 = ee25;
0076 for (int i = 0; i < 26; i++) {
0077 product.eEX0.push_back(eEX0[i]);
0078 }
0079 }
0080
0081 if (eb1 != 0) {
0082 product.eb1 = eb1;
0083 product.eb4 = eb4;
0084 product.eb9 = eb9;
0085 product.eb16 = eb16;
0086 product.eb25 = eb25;
0087 for (int i = 0; i < 26; i++) {
0088 product.eBX0.push_back(eBX0[i]);
0089 }
0090 }
0091
0092 product.totalEInEE = totalEInEE;
0093 product.totalEInEB = totalEInEB;
0094 product.totalEInES = totalEInES;
0095
0096 product.totalEInEEzp = totalEInEEzp;
0097 product.totalEInEEzm = totalEInEEzm;
0098
0099 product.totalEInESzp = totalEInESzp;
0100 product.totalEInESzm = totalEInESzm;
0101
0102 product.totalHits = totalHits;
0103 product.nHitsInEE = nHitsInEE;
0104 product.nHitsInEB = nHitsInEB;
0105 product.nHitsInES = nHitsInES;
0106 product.nHitsIn1ES = nHitsIn1ES;
0107 product.nHitsIn2ES = nHitsIn2ES;
0108 product.nCrystalInEB = nCrystalInEB;
0109 product.nCrystalInEEzp = nCrystalInEEzp;
0110 product.nCrystalInEEzm = nCrystalInEEzm;
0111
0112 product.nHitsIn1ESzp = nHitsIn1ESzp;
0113 product.nHitsIn1ESzm = nHitsIn1ESzm;
0114 product.nHitsIn2ESzp = nHitsIn2ESzp;
0115 product.nHitsIn2ESzm = nHitsIn2ESzm;
0116
0117 product.eOf1ES = eOf1ES;
0118 product.eOf2ES = eOf2ES;
0119 product.zOfES = zOfES;
0120
0121 product.eOf1ESzp = eOf1ESzp;
0122 product.eOf1ESzm = eOf1ESzm;
0123 product.eOf2ESzp = eOf2ESzp;
0124 product.eOf2ESzm = eOf2ESzm;
0125
0126 product.phiOfEECaloG4Hit = phiOfEECaloG4Hit;
0127 product.etaOfEECaloG4Hit = etaOfEECaloG4Hit;
0128 product.eOfEECaloG4Hit = eOfEECaloG4Hit;
0129 product.eOfEEPlusCaloG4Hit = eOfEEPlusCaloG4Hit;
0130 product.eOfEEMinusCaloG4Hit = eOfEEMinusCaloG4Hit;
0131 product.tOfEECaloG4Hit = tOfEECaloG4Hit;
0132
0133 product.phiOfESCaloG4Hit = phiOfESCaloG4Hit;
0134 product.etaOfESCaloG4Hit = etaOfESCaloG4Hit;
0135 product.eOfESCaloG4Hit = eOfESCaloG4Hit;
0136 product.tOfESCaloG4Hit = tOfESCaloG4Hit;
0137
0138 product.phiOfEBCaloG4Hit = phiOfEBCaloG4Hit;
0139 product.etaOfEBCaloG4Hit = etaOfEBCaloG4Hit;
0140 product.eOfEBCaloG4Hit = eOfEBCaloG4Hit;
0141 product.tOfEBCaloG4Hit = tOfEBCaloG4Hit;
0142
0143 product.theMomentum = theMomentum;
0144 product.theVertex = theVertex;
0145 product.thePID = thePID;
0146 }
0147
0148 void EcalSimHitsValidProducer::update(const BeginOfEvent *) {
0149 ee1 = 0.0;
0150 ee4 = 0.0;
0151 ee9 = 0.0;
0152 ee16 = 0.0;
0153 ee25 = 0.0;
0154
0155 eb1 = 0.0;
0156 eb4 = 0.0;
0157 eb9 = 0.0;
0158 eb16 = 0.0;
0159 eb25 = 0.0;
0160
0161 totalEInEE = 0.0;
0162 totalEInEB = 0.0;
0163 totalEInES = 0.0;
0164
0165 totalEInEEzp = 0.0;
0166 totalEInEEzm = 0.0;
0167 totalEInESzp = 0.0;
0168 totalEInESzm = 0.0;
0169
0170 totalHits = 0;
0171 nHitsInEE = 0;
0172 nHitsInEB = 0;
0173 nHitsInES = 0;
0174 nHitsIn1ES = 0;
0175 nHitsIn2ES = 0;
0176 nCrystalInEB = 0;
0177 nCrystalInEEzp = 0;
0178 nCrystalInEEzm = 0;
0179
0180 nHitsIn1ESzp = 0;
0181 nHitsIn1ESzm = 0;
0182 nHitsIn2ESzp = 0;
0183 nHitsIn2ESzm = 0;
0184
0185 for (int i = 0; i < 26; i++) {
0186 eBX0[i] = 0.0;
0187 eEX0[i] = 0.0;
0188 }
0189
0190 eOf1ES.clear();
0191 eOf2ES.clear();
0192 zOfES.clear();
0193
0194 eOf1ESzp.clear();
0195 eOf1ESzm.clear();
0196 eOf2ESzp.clear();
0197 eOf2ESzm.clear();
0198
0199 phiOfEECaloG4Hit.clear();
0200 etaOfEECaloG4Hit.clear();
0201 tOfEECaloG4Hit.clear();
0202 eOfEECaloG4Hit.clear();
0203 eOfEEPlusCaloG4Hit.clear();
0204 eOfEEMinusCaloG4Hit.clear();
0205
0206 phiOfESCaloG4Hit.clear();
0207 etaOfESCaloG4Hit.clear();
0208 tOfESCaloG4Hit.clear();
0209 eOfESCaloG4Hit.clear();
0210
0211 phiOfEBCaloG4Hit.clear();
0212 etaOfEBCaloG4Hit.clear();
0213 tOfEBCaloG4Hit.clear();
0214 eOfEBCaloG4Hit.clear();
0215 }
0216
0217 void EcalSimHitsValidProducer::update(const EndOfEvent *evt) {
0218 int trackID = 0;
0219 G4PrimaryParticle *thePrim = nullptr;
0220 int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
0221 if (nvertex <= 0) {
0222 edm::LogWarning("EcalSimHitsValidProducer") << " No Vertex in this Event!";
0223 } else {
0224 for (int i = 0; i < nvertex; i++) {
0225 G4PrimaryVertex *avertex = (*evt)()->GetPrimaryVertex(i);
0226 if (avertex == nullptr)
0227 edm::LogWarning("EcalSimHitsValidProducer") << " Pointer to vertex is NULL!";
0228 else {
0229 float x0 = avertex->GetX0();
0230 float y0 = avertex->GetY0();
0231 float z0 = avertex->GetZ0();
0232 float t0 = avertex->GetT0();
0233 theVertex.SetCoordinates(x0, y0, z0, t0);
0234
0235 int npart = avertex->GetNumberOfParticle();
0236 if (npart == 0)
0237 edm::LogWarning("EcalSimHitsValidProducer") << " No primary particle in this event";
0238 else {
0239 if (thePrim == nullptr)
0240 thePrim = avertex->GetPrimary(trackID);
0241 }
0242 }
0243 }
0244
0245
0246 double pInit = 0;
0247 if (thePrim != nullptr) {
0248 double px = thePrim->GetPx();
0249 double py = thePrim->GetPy();
0250 double pz = thePrim->GetPz();
0251 theMomentum.SetCoordinates(px, py, pz, 0.);
0252
0253 pInit = std::sqrt(px * px + py * py + pz * pz);
0254 if (pInit == 0)
0255 edm::LogWarning("EcalSimHitsValidProducer") << " Primary has p = 0 ; ";
0256 else {
0257 theMomentum.SetE(pInit);
0258 }
0259
0260 thePID = thePrim->GetPDGcode();
0261 } else {
0262 edm::LogWarning("EcalSimHitsValidProducer") << " Could not find the primary particle!!";
0263 }
0264 }
0265
0266 G4HCofThisEvent *allHC = (*evt)()->GetHCofThisEvent();
0267 int EBHCid = G4SDManager::GetSDMpointer()->GetCollectionID("EcalHitsEB");
0268 int EEHCid = G4SDManager::GetSDMpointer()->GetCollectionID("EcalHitsEE");
0269 int SEHCid = G4SDManager::GetSDMpointer()->GetCollectionID("EcalHitsES");
0270
0271 CaloG4HitCollection *theEBHC = (CaloG4HitCollection *)allHC->GetHC(EBHCid);
0272 CaloG4HitCollection *theEEHC = (CaloG4HitCollection *)allHC->GetHC(EEHCid);
0273 CaloG4HitCollection *theSEHC = (CaloG4HitCollection *)allHC->GetHC(SEHCid);
0274
0275 nHitsInEE = theEEHC->entries();
0276 nHitsInEB = theEBHC->entries();
0277 nHitsInES = theSEHC->entries();
0278 totalHits = nHitsInEE + nHitsInEB + nHitsInES;
0279
0280
0281 MapType ebmap;
0282 int theebhc_entries = theEBHC->entries();
0283 for (int j = 0; j < theebhc_entries; j++) {
0284 CaloG4Hit *aHit = (*theEBHC)[j];
0285 totalEInEB += aHit->getEnergyDeposit();
0286 float he = aHit->getEnergyDeposit();
0287 float htime = aHit->getTimeSlice();
0288
0289 math::XYZPoint hpos = aHit->getEntry();
0290 float htheta = hpos.theta();
0291 float heta = -log(tan(htheta * 0.5));
0292 float hphi = hpos.phi();
0293
0294 phiOfEBCaloG4Hit.push_back(hphi);
0295 etaOfEBCaloG4Hit.push_back(heta);
0296 tOfEBCaloG4Hit.push_back(htime);
0297 eOfEBCaloG4Hit.push_back(he);
0298 uint32_t crystid = aHit->getUnitID();
0299 ebmap[crystid] += aHit->getEnergyDeposit();
0300 }
0301
0302 nCrystalInEB = ebmap.size();
0303
0304
0305 MapType eemap, eezpmap, eezmmap;
0306 int theeehc_entries = theEEHC->entries();
0307 for (int j = 0; j < theeehc_entries; j++) {
0308 CaloG4Hit *aHit = (*theEEHC)[j];
0309 totalEInEE += aHit->getEnergyDeposit();
0310 float he = aHit->getEnergyDeposit();
0311 float htime = aHit->getTimeSlice();
0312
0313 math::XYZPoint hpos = aHit->getEntry();
0314 float htheta = hpos.theta();
0315 float heta = -log(tan(htheta * 0.5));
0316 float hphi = hpos.phi();
0317 phiOfEECaloG4Hit.push_back(hphi);
0318 etaOfEECaloG4Hit.push_back(heta);
0319 tOfEECaloG4Hit.push_back(htime);
0320 eOfEECaloG4Hit.push_back(he);
0321
0322 uint32_t crystid = aHit->getUnitID();
0323 EEDetId myEEid(crystid);
0324 if (myEEid.zside() == -1) {
0325 totalEInEEzm += aHit->getEnergyDeposit();
0326 eOfEEMinusCaloG4Hit.push_back(he);
0327 eezmmap[crystid] += aHit->getEnergyDeposit();
0328 }
0329 if (myEEid.zside() == 1) {
0330 totalEInEEzp += aHit->getEnergyDeposit();
0331 eOfEEPlusCaloG4Hit.push_back(he);
0332 eezpmap[crystid] += aHit->getEnergyDeposit();
0333 }
0334
0335 eemap[crystid] += aHit->getEnergyDeposit();
0336 }
0337
0338 nCrystalInEEzm = eezmmap.size();
0339 nCrystalInEEzp = eezpmap.size();
0340
0341
0342 int thesehc_entries = theSEHC->entries();
0343 for (int j = 0; j < thesehc_entries; j++) {
0344 CaloG4Hit *aHit = (*theSEHC)[j];
0345 totalEInES += aHit->getEnergyDeposit();
0346 ESDetId esid = ESDetId(aHit->getUnitID());
0347
0348 if (esid.zside() == -1) {
0349 totalEInESzm += aHit->getEnergyDeposit();
0350
0351 if (esid.plane() == 1) {
0352 nHitsIn1ESzm++;
0353 eOf1ESzm.push_back(aHit->getEnergyDeposit());
0354 } else if (esid.plane() == 2) {
0355 nHitsIn2ESzm++;
0356 eOf2ESzm.push_back(aHit->getEnergyDeposit());
0357 }
0358 }
0359 if (esid.zside() == 1) {
0360 totalEInESzp += aHit->getEnergyDeposit();
0361
0362 if (esid.plane() == 1) {
0363 nHitsIn1ESzp++;
0364 eOf1ESzp.push_back(aHit->getEnergyDeposit());
0365 } else if (esid.plane() == 2) {
0366 nHitsIn2ESzp++;
0367 eOf2ESzp.push_back(aHit->getEnergyDeposit());
0368 }
0369 }
0370 }
0371
0372 uint32_t eemaxid = getUnitWithMaxEnergy(eemap);
0373 uint32_t ebmaxid = getUnitWithMaxEnergy(ebmap);
0374 if (eemap[eemaxid] > ebmap[ebmaxid]) {
0375 uint32_t centerid = getUnitWithMaxEnergy(eemap);
0376 EEDetId myEEid(centerid);
0377 int ix = myEEid.ix();
0378 int iy = myEEid.iy();
0379 int iz = myEEid.zside();
0380
0381 ee1 = energyInEEMatrix(1, 1, ix, iy, iz, eemap);
0382 ee9 = energyInEEMatrix(3, 3, ix, iy, iz, eemap);
0383 ee25 = energyInEEMatrix(5, 5, ix, iy, iz, eemap);
0384 MapType neweemap;
0385 if (fillEEMatrix(3, 3, ix, iy, iz, neweemap, eemap)) {
0386 ee4 = eCluster2x2(neweemap);
0387 }
0388 if (fillEEMatrix(5, 5, ix, iy, iz, neweemap, eemap)) {
0389 ee16 = eCluster4x4(ee9, neweemap);
0390 }
0391 } else {
0392 uint32_t ebcenterid = getUnitWithMaxEnergy(ebmap);
0393 EBDetId myEBid(ebcenterid);
0394 int bx = myEBid.ietaAbs();
0395 int by = myEBid.iphi();
0396 int bz = myEBid.zside();
0397 eb1 = energyInEBMatrix(1, 1, bx, by, bz, ebmap);
0398 eb9 = energyInEBMatrix(3, 3, bx, by, bz, ebmap);
0399 eb25 = energyInEBMatrix(5, 5, bx, by, bz, ebmap);
0400
0401 MapType newebmap;
0402 if (fillEBMatrix(3, 3, bx, by, bz, newebmap, ebmap)) {
0403 eb4 = eCluster2x2(newebmap);
0404 }
0405 if (fillEBMatrix(5, 5, bx, by, bz, newebmap, ebmap)) {
0406 eb16 = eCluster4x4(eb9, newebmap);
0407 }
0408 }
0409 }
0410
0411 void EcalSimHitsValidProducer::update(const G4Step *aStep) {
0412 G4StepPoint *preStepPoint = aStep->GetPreStepPoint();
0413 const G4ThreeVector &hitPoint = preStepPoint->GetPosition();
0414 G4VPhysicalVolume *currentPV = preStepPoint->GetPhysicalVolume();
0415 const G4String &name = currentPV->GetName();
0416 std::string crystal;
0417 crystal.assign(name, 0, 4);
0418
0419 float Edeposit = aStep->GetTotalEnergyDeposit();
0420 if (crystal == "EFRY" && Edeposit > 0.0) {
0421 float z = hitPoint.z();
0422 float detz = fabs(fabs(z) - 3200);
0423 int x0 = (int)floor(detz / 8.9);
0424 if (x0 < 26) {
0425 eEX0[x0] += Edeposit;
0426 }
0427 }
0428 if (crystal == "EBRY" && Edeposit > 0.0) {
0429 float x = hitPoint.x();
0430 float y = hitPoint.y();
0431 float r = sqrt(x * x + y * y);
0432 float detr = r - 1290;
0433 int x0 = (int)floor(detr / 8.9);
0434 if (x0 < 26) {
0435 eBX0[x0] += Edeposit;
0436 }
0437 }
0438 }
0439
0440 bool EcalSimHitsValidProducer::fillEEMatrix(
0441 int nCellInEta, int nCellInPhi, int CentralEta, int CentralPhi, int CentralZ, MapType &fillmap, MapType &themap) {
0442 int goBackInEta = nCellInEta / 2;
0443 int goBackInPhi = nCellInPhi / 2;
0444
0445 int startEta = CentralEta - goBackInEta;
0446 int startPhi = CentralPhi - goBackInPhi;
0447
0448 int i = 0;
0449 for (int ieta = startEta; ieta < startEta + nCellInEta; ieta++) {
0450 for (int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++) {
0451 uint32_t index;
0452
0453 if (EEDetId::validDetId(ieta, iphi, CentralZ)) {
0454 index = EEDetId(ieta, iphi, CentralZ).rawId();
0455 } else {
0456 continue;
0457 }
0458
0459 fillmap[i++] = themap[index];
0460 }
0461 }
0462 uint32_t centerid = getUnitWithMaxEnergy(themap);
0463
0464 if (fillmap[i / 2] == themap[centerid])
0465 return true;
0466 else
0467 return false;
0468 }
0469
0470 bool EcalSimHitsValidProducer::fillEBMatrix(
0471 int nCellInEta, int nCellInPhi, int CentralEta, int CentralPhi, int CentralZ, MapType &fillmap, MapType &themap) {
0472 int goBackInEta = nCellInEta / 2;
0473 int goBackInPhi = nCellInPhi / 2;
0474
0475 int startEta = CentralZ * CentralEta - goBackInEta;
0476 int startPhi = CentralPhi - goBackInPhi;
0477
0478 int i = 0;
0479 for (int ieta = startEta; ieta < startEta + nCellInEta; ieta++) {
0480 for (int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++) {
0481 uint32_t index;
0482 if (abs(ieta) > 85 || abs(ieta) < 1) {
0483 continue;
0484 }
0485 if (iphi < 1) {
0486 index = EBDetId(ieta, iphi + 360).rawId();
0487 } else if (iphi > 360) {
0488 index = EBDetId(ieta, iphi - 360).rawId();
0489 } else {
0490 index = EBDetId(ieta, iphi).rawId();
0491 }
0492 fillmap[i++] = themap[index];
0493 }
0494 }
0495
0496 uint32_t ebcenterid = getUnitWithMaxEnergy(themap);
0497
0498 if (fillmap[i / 2] == themap[ebcenterid])
0499 return true;
0500 else
0501 return false;
0502 }
0503
0504 float EcalSimHitsValidProducer::eCluster2x2(MapType &themap) {
0505 float E22 = 0.;
0506 float e012 = themap[0] + themap[1] + themap[2];
0507 float e036 = themap[0] + themap[3] + themap[6];
0508 float e678 = themap[6] + themap[7] + themap[8];
0509 float e258 = themap[2] + themap[5] + themap[8];
0510
0511 if ((e012 > e678 || e012 == e678) && (e036 > e258 || e036 == e258))
0512 E22 = themap[0] + themap[1] + themap[3] + themap[4];
0513 else if ((e012 > e678 || e012 == e678) && (e036 < e258 || e036 == e258))
0514 E22 = themap[1] + themap[2] + themap[4] + themap[5];
0515 else if ((e012 < e678 || e012 == e678) && (e036 > e258 || e036 == e258))
0516 E22 = themap[3] + themap[4] + themap[6] + themap[7];
0517 else if ((e012 < e678 || e012 == e678) && (e036 < e258 || e036 == e258))
0518 E22 = themap[4] + themap[5] + themap[7] + themap[8];
0519
0520 return E22;
0521 }
0522
0523 float EcalSimHitsValidProducer::eCluster4x4(float e33, MapType &themap) {
0524 float E44 = 0.;
0525 float e0_4 = themap[0] + themap[1] + themap[2] + themap[3] + themap[4];
0526 float e0_20 = themap[0] + themap[5] + themap[10] + themap[15] + themap[20];
0527 float e4_24 = themap[4] + themap[9] + themap[14] + themap[19] + themap[24];
0528 float e0_24 = themap[20] + themap[21] + themap[22] + themap[23] + themap[24];
0529
0530 if ((e0_4 > e0_24 || e0_4 == e0_24) && (e0_20 > e4_24 || e0_20 == e4_24))
0531 E44 = e33 + themap[0] + themap[1] + themap[2] + themap[3] + themap[5] + themap[10] + themap[15];
0532 else if ((e0_4 > e0_24 || e0_4 == e0_24) && (e0_20 < e4_24 || e0_20 == e4_24))
0533 E44 = e33 + themap[1] + themap[2] + themap[3] + themap[4] + themap[9] + themap[14] + themap[19];
0534 else if ((e0_4 < e0_24 || e0_4 == e0_24) && (e0_20 > e4_24 || e0_20 == e4_24))
0535 E44 = e33 + themap[5] + themap[10] + themap[15] + themap[20] + themap[21] + themap[22] + themap[23];
0536 else if ((e0_4 < e0_24 || e0_4 == e0_24) && (e0_20 < e4_24 || e0_20 == e4_24))
0537 E44 = e33 + themap[21] + themap[22] + themap[23] + themap[24] + themap[9] + themap[14] + themap[19];
0538
0539 return E44;
0540 }
0541
0542 float EcalSimHitsValidProducer::energyInEEMatrix(
0543 int nCellInX, int nCellInY, int centralX, int centralY, int centralZ, MapType &themap) {
0544 int ncristals = 0;
0545 float totalEnergy = 0.;
0546
0547 int goBackInX = nCellInX / 2;
0548 int goBackInY = nCellInY / 2;
0549 int startX = centralX - goBackInX;
0550 int startY = centralY - goBackInY;
0551
0552 for (int ix = startX; ix < startX + nCellInX; ix++) {
0553 for (int iy = startY; iy < startY + nCellInY; iy++) {
0554 uint32_t index;
0555 if (EEDetId::validDetId(ix, iy, centralZ)) {
0556 index = EEDetId(ix, iy, centralZ).rawId();
0557 } else {
0558 continue;
0559 }
0560
0561 totalEnergy += themap[index];
0562 ncristals += 1;
0563
0564 LogDebug("EcalSimHitsValidProducer")
0565 << " EnergyInEEMatrix: ix - iy - E = " << ix << " " << iy << " " << themap[index];
0566 }
0567 }
0568
0569 LogDebug("EcalSimHitsValidProducer") << " EnergyInEEMatrix: energy in " << nCellInX << " cells in x times "
0570 << nCellInY << " cells in y matrix = " << totalEnergy << " for " << ncristals
0571 << " crystals";
0572
0573 return totalEnergy;
0574 }
0575
0576 float EcalSimHitsValidProducer::energyInEBMatrix(
0577 int nCellInEta, int nCellInPhi, int centralEta, int centralPhi, int centralZ, MapType &themap) {
0578 int ncristals = 0;
0579 float totalEnergy = 0.;
0580
0581 int goBackInEta = nCellInEta / 2;
0582 int goBackInPhi = nCellInPhi / 2;
0583 int startEta = centralZ * centralEta - goBackInEta;
0584 int startPhi = centralPhi - goBackInPhi;
0585
0586 for (int ieta = startEta; ieta < startEta + nCellInEta; ieta++) {
0587 for (int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++) {
0588 uint32_t index;
0589 if (abs(ieta) > 85 || abs(ieta) < 1) {
0590 continue;
0591 }
0592 if (iphi < 1) {
0593 index = EBDetId(ieta, iphi + 360).rawId();
0594 } else if (iphi > 360) {
0595 index = EBDetId(ieta, iphi - 360).rawId();
0596 } else {
0597 index = EBDetId(ieta, iphi).rawId();
0598 }
0599
0600 totalEnergy += themap[index];
0601 ncristals += 1;
0602
0603 LogDebug("EcalSimHitsValidProducer")
0604 << " EnergyInEBMatrix: ieta - iphi - E = " << ieta << " " << iphi << " " << themap[index];
0605 }
0606 }
0607
0608 LogDebug("EcalSimHitsValidProducer") << " EnergyInEBMatrix: energy in " << nCellInEta << " cells in eta times "
0609 << nCellInPhi << " cells in phi matrix = " << totalEnergy << " for " << ncristals
0610 << " crystals";
0611
0612 return totalEnergy;
0613 }
0614
0615 uint32_t EcalSimHitsValidProducer::getUnitWithMaxEnergy(MapType &themap) {
0616 uint32_t unitWithMaxEnergy = 0;
0617 float maxEnergy = 0.;
0618
0619 MapType::iterator iter;
0620 for (iter = themap.begin(); iter != themap.end(); iter++) {
0621 if (maxEnergy < (*iter).second) {
0622 maxEnergy = (*iter).second;
0623 unitWithMaxEnergy = (*iter).first;
0624 }
0625 }
0626
0627 LogDebug("EcalSimHitsValidProducer") << " Max energy of " << maxEnergy << " MeV was found in Unit id 0x" << std::hex
0628 << unitWithMaxEnergy << std::dec;
0629
0630 return unitWithMaxEnergy;
0631 }
