Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-18 08:23:52

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 //#define EDM_ML_DEBUG
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   // get handles to calogeometry
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 //define this as a plug-in
0458 DEFINE_FWK_MODULE(EcalSimHitStudy);