File indexing completed on 2024-04-06 12:29:49
0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Utilities/interface/InputTag.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/Utilities/interface/Exception.h"
0012 #include "FWCore/Utilities/interface/transform.h"
0013
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0016
0017 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0018 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0019 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0020 #include "DataFormats/Math/interface/Point3D.h"
0021 #include "DataFormats/Math/interface/Vector3D.h"
0022
0023 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0024 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0025 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0026 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0027 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0028 #include "SimG4CMS/Calo/interface/CaloHitID.h"
0029
0030 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0031 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0032 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0033
0034 #include <TH1F.h>
0035 #include <TH2F.h>
0036 #include <cmath>
0037 #include <iostream>
0038 #include <fstream>
0039 #include <map>
0040 #include <memory>
0041 #include <string>
0042 #include <vector>
0043
0044
0045
0046 class EcalSimHitStudy : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0047 public:
0048 EcalSimHitStudy(const edm::ParameterSet& ps);
0049 ~EcalSimHitStudy() override {}
0050 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0051
0052 protected:
0053 void beginJob() override;
0054 void analyze(edm::Event const&, edm::EventSetup const&) override;
0055 void beginRun(edm::Run const&, edm::EventSetup const&) override {}
0056 void endRun(edm::Run const&, edm::EventSetup const&) override {}
0057 void analyzeHits(std::vector<PCaloHit>&, int);
0058
0059 private:
0060 struct EcalHit {
0061 uint16_t id;
0062 double time, energy;
0063 EcalHit(uint16_t i = 0, double t = 0, double e = 0) : id(i), time(t), energy(e) {}
0064 };
0065 static const int ndets_ = 2;
0066 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tokGeom_;
0067 const std::string g4Label_;
0068 const std::vector<std::string> hitLab_;
0069 const double maxEnergy_, tmax_, w0_;
0070 const int selX_;
0071 const edm::EDGetTokenT<edm::HepMCProduct> tok_evt_;
0072 const std::vector<edm::EDGetTokenT<edm::PCaloHitContainer>> toks_calo_;
0073 const CaloGeometry* geometry_;
0074 TH1F *ptInc_, *etaInc_, *phiInc_, *eneInc_;
0075 TH1F *hit_[ndets_], *time_[ndets_], *timeAll_[ndets_];
0076 TH1F *edepEM_[ndets_], *edepHad_[ndets_], *edep_[ndets_];
0077 TH1F *etot_[ndets_], *etotg_[ndets_], *edepAll_[ndets_];
0078 TH1F *r1by9_[ndets_], *r1by25_[ndets_], *r9by25_[ndets_];
0079 TH1F *sEtaEta_[ndets_], *sEtaPhi_[ndets_], *sPhiPhi_[ndets_];
0080 TH2F *poszp_[ndets_], *poszn_[ndets_];
0081 };
0082
0083 EcalSimHitStudy::EcalSimHitStudy(const edm::ParameterSet& ps)
0084 : tokGeom_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0085 g4Label_(ps.getUntrackedParameter<std::string>("ModuleLabel", "g4SimHits")),
0086 hitLab_(ps.getUntrackedParameter<std::vector<std::string>>("CaloCollection")),
0087 maxEnergy_(ps.getUntrackedParameter<double>("MaxEnergy", 200.0)),
0088 tmax_(ps.getUntrackedParameter<double>("TimeCut", 100.0)),
0089 w0_(ps.getUntrackedParameter<double>("W0", 4.7)),
0090 selX_(ps.getUntrackedParameter<int>("SelectX", -1)),
0091 tok_evt_(consumes<edm::HepMCProduct>(
0092 edm::InputTag(ps.getUntrackedParameter<std::string>("SourceLabel", "VtxSmeared")))),
0093 toks_calo_{edm::vector_transform(hitLab_, [this](const std::string& name) {
0094 return consumes<edm::PCaloHitContainer>(edm::InputTag{g4Label_, name});
0095 })} {
0096 usesResource(TFileService::kSharedResource);
0097
0098 edm::LogVerbatim("HitStudy") << "Module Label: " << g4Label_ << " Hits: " << hitLab_[0] << ", " << hitLab_[1]
0099 << " MaxEnergy: " << maxEnergy_ << " Tmax: " << tmax_ << " Select " << selX_;
0100 }
0101
0102 void EcalSimHitStudy::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0103 edm::ParameterSetDescription desc;
0104 std::vector<std::string> calonames;
0105 desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
0106 desc.addUntracked<std::vector<std::string>>("CaloCollection", calonames);
0107 desc.addUntracked<std::string>("SourceLabel", "VtxSmeared");
0108 desc.addUntracked<double>("MaxEnergy", 200.0);
0109 desc.addUntracked<double>("TimeCut", 100.0);
0110 desc.addUntracked<int>("SelectX", -1);
0111 descriptions.add("EcalSimHitStudy", desc);
0112 }
0113
0114 void EcalSimHitStudy::beginJob() {
0115 edm::Service<TFileService> tfile;
0116 if (!tfile.isAvailable())
0117 throw cms::Exception("BadConfig") << "TFileService unavailable: "
0118 << "please add it to config file";
0119 char name[20], title[120];
0120 sprintf(title, "Incident PT (GeV)");
0121 ptInc_ = tfile->make<TH1F>("PtInc", title, 1000, 0., maxEnergy_);
0122 ptInc_->GetXaxis()->SetTitle(title);
0123 ptInc_->GetYaxis()->SetTitle("Events");
0124 sprintf(title, "Incident Energy (GeV)");
0125 eneInc_ = tfile->make<TH1F>("EneInc", title, 1000, 0., maxEnergy_);
0126 eneInc_->GetXaxis()->SetTitle(title);
0127 eneInc_->GetYaxis()->SetTitle("Events");
0128 sprintf(title, "Incident #eta");
0129 etaInc_ = tfile->make<TH1F>("EtaInc", title, 200, -5., 5.);
0130 etaInc_->GetXaxis()->SetTitle(title);
0131 etaInc_->GetYaxis()->SetTitle("Events");
0132 sprintf(title, "Incident #phi");
0133 phiInc_ = tfile->make<TH1F>("PhiInc", title, 200, -3.1415926, 3.1415926);
0134 phiInc_->GetXaxis()->SetTitle(title);
0135 phiInc_->GetYaxis()->SetTitle("Events");
0136 std::string dets[ndets_] = {"EB", "EE"};
0137 for (int i = 0; i < ndets_; i++) {
0138 sprintf(name, "Hit%d", i);
0139 sprintf(title, "Number of hits in %s", dets[i].c_str());
0140 hit_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 1000, 0., 20000.);
0141 hit_[i]->GetXaxis()->SetTitle(title);
0142 hit_[i]->GetYaxis()->SetTitle("Events");
0143 hit_[i]->Sumw2();
0144 sprintf(name, "Time%d", i);
0145 sprintf(title, "Time of the hit (ns) in %s", dets[i].c_str());
0146 time_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 1000, 0., 1000.);
0147 time_[i]->GetXaxis()->SetTitle(title);
0148 time_[i]->GetYaxis()->SetTitle("Hits");
0149 time_[i]->Sumw2();
0150 sprintf(name, "TimeAll%d", i);
0151 sprintf(title, "Hit time (ns) in %s (for first hit in the cell)", dets[i].c_str());
0152 timeAll_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 1000, 0., 1000.);
0153 timeAll_[i]->GetXaxis()->SetTitle(title);
0154 timeAll_[i]->GetYaxis()->SetTitle("Hits");
0155 timeAll_[i]->Sumw2();
0156 sprintf(name, "Edep%d", i);
0157 sprintf(title, "Energy deposit (GeV) in %s", dets[i].c_str());
0158 edep_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 5000, 0., maxEnergy_);
0159 edep_[i]->GetXaxis()->SetTitle(title);
0160 edep_[i]->GetYaxis()->SetTitle("Hits");
0161 edep_[i]->Sumw2();
0162 sprintf(name, "EdepAll%d", i);
0163 sprintf(title, "Total Energy deposit in the cell (GeV) in %s", dets[i].c_str());
0164 edepAll_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 5000, 0., maxEnergy_);
0165 edepAll_[i]->GetXaxis()->SetTitle(title);
0166 edepAll_[i]->GetYaxis()->SetTitle("Hits");
0167 edepAll_[i]->Sumw2();
0168 sprintf(name, "EdepEM%d", i);
0169 sprintf(title, "Energy deposit (GeV) by EM particles in %s", dets[i].c_str());
0170 edepEM_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 5000, 0., maxEnergy_);
0171 edepEM_[i]->GetXaxis()->SetTitle(title);
0172 edepEM_[i]->GetYaxis()->SetTitle("Hits");
0173 edepEM_[i]->Sumw2();
0174 sprintf(name, "EdepHad%d", i);
0175 sprintf(title, "Energy deposit (GeV) by hadrons in %s", dets[i].c_str());
0176 edepHad_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 5000, 0., maxEnergy_);
0177 edepHad_[i]->GetXaxis()->SetTitle(title);
0178 edepHad_[i]->GetYaxis()->SetTitle("Hits");
0179 edepHad_[i]->Sumw2();
0180 sprintf(name, "Etot%d", i);
0181 sprintf(title, "Total energy deposit (GeV) in %s", dets[i].c_str());
0182 etot_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 5000, 0., maxEnergy_);
0183 etot_[i]->GetXaxis()->SetTitle(title);
0184 etot_[i]->GetYaxis()->SetTitle("Events");
0185 etot_[i]->Sumw2();
0186 sprintf(name, "EtotG%d", i);
0187 sprintf(title, "Total energy deposit (GeV) in %s (t < 100 ns)", dets[i].c_str());
0188 etotg_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 5000, 0., maxEnergy_);
0189 etotg_[i]->GetXaxis()->SetTitle(title);
0190 etotg_[i]->GetYaxis()->SetTitle("Events");
0191 etotg_[i]->Sumw2();
0192 sprintf(name, "r1by9%d", i);
0193 sprintf(title, "E1/E9 in %s", dets[i].c_str());
0194 r1by9_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 100, 0.0, 1.0);
0195 r1by9_[i]->GetXaxis()->SetTitle(title);
0196 r1by9_[i]->GetYaxis()->SetTitle("Events");
0197 r1by9_[i]->Sumw2();
0198 sprintf(name, "r1by25%d", i);
0199 sprintf(title, "E1/E25 in %s", dets[i].c_str());
0200 r1by25_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 100, 0.0, 1.0);
0201 r1by25_[i]->GetXaxis()->SetTitle(title);
0202 r1by25_[i]->GetYaxis()->SetTitle("Events");
0203 r1by25_[i]->Sumw2();
0204 sprintf(name, "r9by25%d", i);
0205 sprintf(title, "E9/E25 in %s", dets[i].c_str());
0206 r9by25_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 100, 0.0, 1.0);
0207 r9by25_[i]->GetXaxis()->SetTitle(title);
0208 r9by25_[i]->GetYaxis()->SetTitle("Events");
0209 r9by25_[i]->Sumw2();
0210 double ymax = (i == 0) ? 0.0005 : 0.005;
0211 sprintf(name, "sEtaEta%d", i);
0212 sprintf(title, "Cov(#eta,#eta) in %s", dets[i].c_str());
0213 sEtaEta_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 1000, 0.0, ymax);
0214 sEtaEta_[i]->GetXaxis()->SetTitle(title);
0215 sEtaEta_[i]->GetYaxis()->SetTitle("Events");
0216 sEtaEta_[i]->Sumw2();
0217 sprintf(name, "sEtaPhi%d", i);
0218 sprintf(title, "Cov(#eta,#phi) in %s", dets[i].c_str());
0219 sEtaPhi_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 1000, 0.0, ymax);
0220 sEtaPhi_[i]->GetXaxis()->SetTitle(title);
0221 sEtaPhi_[i]->GetYaxis()->SetTitle("Events");
0222 sEtaPhi_[i]->Sumw2();
0223 ymax = (i == 0) ? 0.001 : 0.01;
0224 sprintf(name, "sPhiPhi%d", i);
0225 sprintf(title, "Cov(#phi,#phi) in %s", dets[i].c_str());
0226 sPhiPhi_[i] = tfile->make<TH1F>(name, dets[i].c_str(), 1000, 0.0, ymax);
0227 sPhiPhi_[i]->GetXaxis()->SetTitle(title);
0228 sPhiPhi_[i]->GetYaxis()->SetTitle("Events");
0229 sPhiPhi_[i]->Sumw2();
0230 if (i == 0) {
0231 sprintf(title, "%s+", dets[i].c_str());
0232 poszp_[i] = tfile->make<TH2F>("poszp0", title, 100, 0, 100, 360, 0, 360);
0233 poszp_[i]->GetXaxis()->SetTitle("i#eta");
0234 poszp_[i]->GetYaxis()->SetTitle("i#phi");
0235 sprintf(title, "%s-", dets[i].c_str());
0236 poszn_[i] = tfile->make<TH2F>("poszn0", title, 100, 0, 100, 360, 0, 360);
0237 poszn_[i]->GetXaxis()->SetTitle("i#eta");
0238 poszn_[i]->GetYaxis()->SetTitle("i#phi");
0239 } else {
0240 sprintf(title, "%s+", dets[i].c_str());
0241 poszp_[i] = tfile->make<TH2F>("poszp1", title, 100, -200, 200, 100, -200, 200);
0242 poszp_[i]->GetXaxis()->SetTitle("x (cm)");
0243 poszp_[i]->GetYaxis()->SetTitle("y (cm)");
0244 sprintf(title, "%s-", dets[i].c_str());
0245 poszn_[i] = tfile->make<TH2F>("poszn1", title, 100, -200, 200, 100, -200, 200);
0246 poszn_[i]->GetXaxis()->SetTitle("x (cm)");
0247 poszn_[i]->GetYaxis()->SetTitle("y (cm)");
0248 }
0249 poszp_[i]->GetYaxis()->SetTitleOffset(1.2);
0250 poszp_[i]->Sumw2();
0251 poszn_[i]->GetYaxis()->SetTitleOffset(1.2);
0252 poszn_[i]->Sumw2();
0253 }
0254 }
0255
0256 void EcalSimHitStudy::analyze(const edm::Event& e, const edm::EventSetup& iS) {
0257 #ifdef EDM_ML_DEBUG
0258 edm::LogVerbatim("HitStudy") << "Run = " << e.id().run() << " Event = " << e.id().event();
0259 #endif
0260
0261 geometry_ = &iS.getData(tokGeom_);
0262
0263 double eInc = 0, etaInc = 0, phiInc = 0;
0264 int type(-1);
0265 edm::Handle<edm::HepMCProduct> EvtHandle;
0266 e.getByToken(tok_evt_, EvtHandle);
0267 if (EvtHandle.isValid()) {
0268 const HepMC::GenEvent* myGenEvent = EvtHandle->GetEvent();
0269
0270 HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin();
0271 if (p != myGenEvent->particles_end()) {
0272 eInc = (*p)->momentum().e();
0273 etaInc = (*p)->momentum().eta();
0274 phiInc = (*p)->momentum().phi();
0275 }
0276 double ptInc = eInc / std::cosh(etaInc);
0277 ptInc_->Fill(ptInc);
0278 eneInc_->Fill(eInc);
0279 etaInc_->Fill(etaInc);
0280 phiInc_->Fill(phiInc);
0281
0282 if (std::abs(etaInc) < 1.46)
0283 type = 0;
0284 else if (std::abs(etaInc) > 1.49 && std::abs(etaInc) < 3.0)
0285 type = 1;
0286 }
0287
0288 int typeMin = (type < 0) ? 0 : type;
0289 int typeMax = (type < 0) ? 1 : type;
0290 for (int type = typeMin; type <= typeMax; ++type) {
0291 edm::Handle<edm::PCaloHitContainer> hitsCalo;
0292 e.getByToken(toks_calo_[type], hitsCalo);
0293 bool getHits = (hitsCalo.isValid());
0294 #ifdef EDM_ML_DEBUG
0295 edm::LogVerbatim("HitStudy") << "EcalSimHitStudy: Input flags Hits " << getHits << " with " << hitsCalo->size()
0296 << " hits";
0297 #endif
0298 if (getHits) {
0299 std::vector<PCaloHit> caloHits;
0300 caloHits.insert(caloHits.end(), hitsCalo->begin(), hitsCalo->end());
0301 if (!caloHits.empty())
0302 analyzeHits(caloHits, type);
0303 }
0304 }
0305 }
0306
0307 void EcalSimHitStudy::analyzeHits(std::vector<PCaloHit>& hits, int indx) {
0308 unsigned int nEC(0);
0309 std::map<unsigned int, EcalHit> hitMap;
0310 double etot(0), etotG(0);
0311 for (auto hit : hits) {
0312 double edep = hit.energy();
0313 double time = hit.time();
0314 unsigned int id_ = hit.id();
0315 double edepEM = hit.energyEM();
0316 double edepHad = hit.energyHad();
0317 if (indx == 0) {
0318 if ((hit.depth() == 1) || (hit.depth() == 2))
0319 continue;
0320 }
0321 if (time <= tmax_) {
0322 auto it = hitMap.find(id_);
0323 if (it == hitMap.end()) {
0324 uint16_t dep = hit.depth();
0325 hitMap[id_] = EcalHit(dep, time, edep);
0326 } else {
0327 (it->second).energy += edep;
0328 }
0329 etotG += edep;
0330 ++nEC;
0331 }
0332 time_[indx]->Fill(time);
0333 edep_[indx]->Fill(edep);
0334 edepEM_[indx]->Fill(edepEM);
0335 edepHad_[indx]->Fill(edepHad);
0336 etot += edep;
0337 }
0338
0339 double edepM(0);
0340 unsigned int idM(0);
0341 uint16_t depM(0);
0342 for (auto it : hitMap) {
0343 if (it.second.energy > edepM) {
0344 idM = it.first;
0345 edepM = it.second.energy;
0346 depM = it.second.id;
0347 }
0348 }
0349
0350 bool select(true);
0351 if (selX_ >= 0) {
0352 if ((depM & 0X4) != 0)
0353 select = (selX_ > 0);
0354 else
0355 select = (selX_ == 0);
0356 }
0357 #ifdef EDM_ML_DEBUG
0358 edm::LogVerbatim("HitStudy") << "EcalSimHitStudy::analyzeHits: Index " << indx << " Emax " << edepM << " IDMax "
0359 << std::hex << idM << ":" << depM << std::dec << " Select " << select << ":" << selX_
0360 << " Hits " << hits.size() << ":" << nEC << ":" << hitMap.size() << " ETotal " << etot
0361 << ":" << etotG;
0362 #endif
0363 if (select) {
0364 etot_[indx]->Fill(etot);
0365 etotg_[indx]->Fill(etotG);
0366 hit_[indx]->Fill(double(nEC));
0367 for (auto it : hitMap) {
0368 timeAll_[indx]->Fill((it.second).time);
0369 edepAll_[indx]->Fill((it.second).energy);
0370 DetId id(it.first);
0371 if (indx == 0) {
0372 if (EBDetId(id).zside() >= 0)
0373 poszp_[indx]->Fill(EBDetId(id).ietaAbs(), EBDetId(id).iphi());
0374 else
0375 poszn_[indx]->Fill(EBDetId(id).ietaAbs(), EBDetId(id).iphi());
0376 } else {
0377 GlobalPoint gpos = geometry_->getGeometry(id)->getPosition();
0378 if (EEDetId(id).zside() >= 0)
0379 poszp_[indx]->Fill(gpos.x(), gpos.y());
0380 else
0381 poszn_[indx]->Fill(gpos.x(), gpos.y());
0382 }
0383 }
0384
0385 math::XYZVector meanPosition(0.0, 0.0, 0.0);
0386 std::vector<math::XYZVector> position;
0387 std::vector<double> energy;
0388 double e9(0), e25(0);
0389 for (auto it : hitMap) {
0390 DetId id(it.first);
0391 int deta(99), dphi(99), dz(0);
0392 if (indx == 0) {
0393 deta = std::abs(EBDetId(id).ietaAbs() - EBDetId(idM).ietaAbs());
0394 dphi = std::abs(EBDetId(id).iphi() - EBDetId(idM).iphi());
0395 if (dphi > 180)
0396 dphi = std::abs(dphi - 360);
0397 dz = std::abs(EBDetId(id).zside() - EBDetId(idM).zside());
0398 } else {
0399 deta = std::abs(EEDetId(id).ix() - EEDetId(idM).ix());
0400 dphi = std::abs(EEDetId(id).iy() - EEDetId(idM).iy());
0401 dz = std::abs(EEDetId(id).zside() - EEDetId(idM).zside());
0402 }
0403 if (deta <= 1 && dphi <= 1 && dz < 1)
0404 e9 += (it.second).energy;
0405 if (deta <= 2 && dphi <= 2 && dz < 1) {
0406 e25 += (it.second).energy;
0407 GlobalPoint gpos = geometry_->getGeometry(id)->getPosition();
0408 math::XYZVector pos(gpos.x(), gpos.y(), gpos.z());
0409 meanPosition += (it.second).energy * pos;
0410 position.push_back(pos);
0411 energy.push_back((it.second).energy);
0412 }
0413 }
0414 double r1by9 = (e9 > 0) ? (edepM / e9) : -1;
0415 double r1by25 = (e25 > 0) ? (edepM / e25) : -1;
0416 double r9by25 = (e25 > 0) ? (e9 / e25) : -1;
0417
0418 meanPosition /= e25;
0419 double denom(0), numEtaEta(0), numEtaPhi(0), numPhiPhi(0);
0420 for (unsigned int k = 0; k < position.size(); ++k) {
0421 double dEta = position[k].eta() - meanPosition.eta();
0422 double dPhi = position[k].phi() - meanPosition.phi();
0423 if (dPhi > +M_PI) {
0424 dPhi = 2 * M_PI - dPhi;
0425 }
0426 if (dPhi < -M_PI) {
0427 dPhi = 2 * M_PI + dPhi;
0428 }
0429
0430 double w = std::max(0.0, (w0_ + std::log(energy[k] / e25)));
0431 denom += w;
0432 numEtaEta += std::abs(w * dEta * dEta);
0433 numEtaPhi += std::abs(w * dEta * dPhi);
0434 numPhiPhi += std::abs(w * dPhi * dPhi);
0435 #ifdef EDM_ML_DEBUG
0436 edm::LogVerbatim("HitStudy") << "[" << k << "] dEta " << dEta << " dPhi " << dPhi << " Wt " << energy[k] / e25
0437 << ":" << std::log(energy[k] / e25) << ":" << w;
0438 #endif
0439 }
0440 double sEtaEta = (denom > 0) ? (numEtaEta / denom) : -1.0;
0441 double sEtaPhi = (denom > 0) ? (numEtaPhi / denom) : -1.0;
0442 double sPhiPhi = (denom > 0) ? (numPhiPhi / denom) : -1.0;
0443
0444 #ifdef EDM_ML_DEBUG
0445 edm::LogVerbatim("HitStudy") << "EcalSimHitStudy::Ratios " << r1by9 << " : " << r1by25 << " : " << r9by25
0446 << " Covariances " << sEtaEta << " : " << sEtaPhi << " : " << sPhiPhi;
0447 #endif
0448 r1by9_[indx]->Fill(r1by9);
0449 r1by25_[indx]->Fill(r1by25);
0450 r9by25_[indx]->Fill(r9by25);
0451 sEtaEta_[indx]->Fill(sEtaEta);
0452 sEtaPhi_[indx]->Fill(sEtaPhi);
0453 sPhiPhi_[indx]->Fill(sPhiPhi);
0454 }
0455 }
0456
0457
0458 DEFINE_FWK_MODULE(EcalSimHitStudy);
