Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-02 03:39:40

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     // the direction of momentum of primary particles
0246     double pInit = 0;  // etaInit =0, phiInit =0, // UNUSED
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   // hit map for EB for matrices
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   //   EB Hit collection start
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   //   EE Hit collection start
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   // Hits from ES
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 }