EcalRecProducers/test/EcalLocalRecoTask.cc

0001 /*
0002  * \file EcalLocalReco.cc
0003  *
0004  *
0005 */
0006
0007 #include "RecoLocalCalo/EcalRecProducers/test/EcalLocalRecoTask.h"
0008 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0009 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0010 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0011
0012 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0013
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0017
0018 EcalLocalRecoTask::EcalLocalRecoTask(const edm::ParameterSet& ps)
0019     : verbose_(ps.getUntrackedParameter<bool>("verbose", false)),
0020       outputFile_(ps.getUntrackedParameter<std::string>("outputFile", "")),  // DQM ROOT output
0021       EBrecHitToken_(consumes<EBRecHitCollection>(ps.getParameter<edm::InputTag>("ebrechits"))),
0022       EErecHitToken_(consumes<EERecHitCollection>(ps.getParameter<edm::InputTag>("eerechits"))),
0023       ESrecHitToken_(consumes<ESRecHitCollection>(ps.getParameter<edm::InputTag>("esrechits"))),
0024       EBurecHitToken_(consumes<EBUncalibratedRecHitCollection>(ps.getParameter<edm::InputTag>("eburechits"))),
0025       EEurecHitToken_(consumes<EEUncalibratedRecHitCollection>(ps.getParameter<edm::InputTag>("eeurechits"))),
0026       EBdigiToken_(consumes<EBDigiCollection>(ps.getParameter<edm::InputTag>("ebdigis"))),
0027       EEdigiToken_(consumes<EEDigiCollection>(ps.getParameter<edm::InputTag>("eedigis"))),
0028       ESdigiToken_(consumes<ESDigiCollection>(ps.getParameter<edm::InputTag>("esdigis"))),
0029       cfToken_(consumes<CrossingFrame<PCaloHit>>(edm::InputTag("mix", "EcalHitsEB"))),
0030       pedestalToken_(esConsumes()) {
0031   if (outputFile_.size() != 0) {
0032     edm::LogInfo("EcalLocalRecoTaskInfo") << "histograms will be saved to '" << outputFile_.c_str() << "'";
0033   } else {
0034     edm::LogInfo("EcalLocalRecoTaskInfo") << "histograms will NOT be saved";
0035   }
0036
0037   if (verbose_) {
0038     edm::LogInfo("EcalLocalRecoTaskInfo") << "verbose switch is ON";
0039   } else {
0040     edm::LogInfo("EcalLocalRecoTaskInfo") << "verbose switch is OFF";
0041   }
0042
0043   dbe_ = 0;
0044
0045   // get hold of back-end interface
0046   dbe_ = edm::Service<DQMStore>().operator->();
0047
0048   meEBUncalibRecHitMaxSampleRatio_ = 0;
0049   meEBUncalibRecHitPedestal_ = 0;
0050   meEBUncalibRecHitOccupancy_ = 0;
0051   meEBRecHitSimHitRatio_ = 0;
0052
0053   Char_t histo[70];
0054
0055   if (dbe_) {
0056     dbe_->setCurrentFolder("EcalLocalRecoTask");
0057
0058     sprintf(histo, "EcalLocalRecoTask Barrel occupancy");
0059     meEBUncalibRecHitOccupancy_ = dbe_->book2D(histo, histo, 360, 0., 360., 170, -85., 85.);
0060
0061     sprintf(histo, "EcalLocalRecoTask Barrel Reco pedestals");
0062     meEBUncalibRecHitPedestal_ = dbe_->book1D(histo, histo, 1000, 0., 1000.);
0063
0064     sprintf(histo, "EcalLocalRecoTask RecHit Max Sample Ratio");
0065     meEBUncalibRecHitMaxSampleRatio_ = dbe_->book1D(histo, histo, 200, 0., 2.);
0066
0067     sprintf(histo, "EcalLocalRecoTask RecHit SimHit Ratio");
0068     meEBRecHitSimHitRatio_ = dbe_->book1D(histo, histo, 200, 0., 2.);
0069   }
0070 }
0071
0072 EcalLocalRecoTask::~EcalLocalRecoTask() {
0073   if (outputFile_.size() != 0 && dbe_)
0074     dbe_->save(outputFile_);
0075 }
0076
0077 void EcalLocalRecoTask::analyze(const edm::Event& e, const edm::EventSetup& c) {
0078   edm::LogInfo("EcalLocalRecoTaskInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
0079
0080   edm::Handle<EBDigiCollection> pEBDigis;
0081   edm::Handle<EEDigiCollection> pEEDigis;
0082   edm::Handle<ESDigiCollection> pESDigis;
0083
0084   e.getByToken(EBdigiToken_, pEBDigis);
0085   e.getByToken(EEdigiToken_, pEEDigis);
0086   e.getByToken(ESdigiToken_, pESDigis);
0087
0088   edm::Handle<EBUncalibratedRecHitCollection> pEBUncalibRecHit;
0089   edm::Handle<EEUncalibratedRecHitCollection> pEEUncalibRecHit;
0090
0091   e.getByToken(EBurecHitToken_, pEBUncalibRecHit);
0092   e.getByToken(EEurecHitToken_, pEEUncalibRecHit);
0093
0094   edm::Handle<EBRecHitCollection> pEBRecHit;
0095   edm::Handle<EERecHitCollection> pEERecHit;
0096   edm::Handle<ESRecHitCollection> pESRecHit;
0097
0098   e.getByToken(EBrecHitToken_, pEBRecHit);
0099   e.getByToken(EBrecHitToken_, pEERecHit);
0100   e.getByToken(ESrecHitToken_, pESRecHit);
0101
0102   edm::Handle<CrossingFrame<PCaloHit>> crossingFrame;
0103
0104   e.getByToken(cfToken_, crossingFrame);
0105
0106   auto barrelHits = std::make_unique<MixCollection<PCaloHit>>(crossingFrame.product());
0107
0108   MapType EBSimMap;
0109
0110   for (MixCollection<PCaloHit>::MixItr hitItr = barrelHits->begin(); hitItr != barrelHits->end(); ++hitItr) {
0111     EBDetId EBid = EBDetId(hitItr->id());
0112
0113     LogDebug("EcalLocalRecoTaskDebug") << " CaloHit " << hitItr->getName() << " DetID = " << hitItr->id() << "\n"
0114                                        << "Energy = " << hitItr->energy() << " Time = " << hitItr->time() << "\n"
0115                                        << "EBDetId = " << EBid.ieta() << " " << EBid.iphi();
0116
0117     uint32_t crystid = EBid.rawId();
0118     EBSimMap[crystid] += hitItr->energy();
0119   }
0120
0121   const EBDigiCollection* EBDigi = pEBDigis.product();
0122   const EBUncalibratedRecHitCollection* EBUncalibRecHit = pEBUncalibRecHit.product();
0123   const EBRecHitCollection* EBRecHit = pEBRecHit.product();
0124
0125   // loop over uncalibRecHit
0126   for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EBUncalibRecHit->begin();
0127        uncalibRecHit != EBUncalibRecHit->end();
0128        ++uncalibRecHit) {
0129     EBDetId EBid = EBDetId(uncalibRecHit->id());
0130     if (meEBUncalibRecHitOccupancy_)
0131       meEBUncalibRecHitOccupancy_->Fill(EBid.iphi(), EBid.ieta());
0132     if (meEBUncalibRecHitPedestal_)
0133       meEBUncalibRecHitPedestal_->Fill(uncalibRecHit->pedestal());
0134
0135     // Find corresponding recHit
0136     EcalRecHitCollection::const_iterator myRecHit = EBRecHit->find(EBid);
0137     // Find corresponding digi
0138     EBDigiCollection::const_iterator myDigi = EBDigi->find(EBid);
0139
0140     double eMax = 0.;
0141
0142     if (myDigi != EBDigi->end()) {
0143       for (unsigned int sample = 0; sample < myDigi->size(); ++sample) {
0144         double analogSample = EcalMGPASample((*myDigi)[sample]).adc();
0145         if (eMax < analogSample) {
0146           eMax = analogSample;
0147         }
0148       }
0149     } else
0150       continue;
0151
0152     const auto& myped = c.getData(pedestalToken_);
0153     EcalPedestals::const_iterator it = myped.getMap().find(EBid);
0154     if (it != myped.getMap().end()) {
0155       if (eMax > (*it).mean_x1 + 5 * (*it).rms_x1)  //only real signal RecHit
0156       {
0157         if (meEBUncalibRecHitMaxSampleRatio_)
0158           meEBUncalibRecHitMaxSampleRatio_->Fill((uncalibRecHit->amplitude() + uncalibRecHit->pedestal()) / eMax);
0159         edm::LogInfo("EcalLocalRecoTaskInfo")
0160             << " eMax = " << eMax << " Amplitude = " << uncalibRecHit->amplitude() + uncalibRecHit->pedestal();
0161       } else
0162         continue;
0163     } else
0164       continue;
0165
0166     if (myRecHit != EBRecHit->end()) {
0167       if (EBSimMap[EBid.rawId()] != 0.)
0168         if (meEBRecHitSimHitRatio_)
0169           meEBRecHitSimHitRatio_->Fill(myRecHit->energy() / EBSimMap[EBid.rawId()]);
0170     } else
0171       continue;
0172   }
0173 }
0174
0175 //define this as a plug-in
0176 DEFINE_FWK_MODULE(EcalLocalRecoTask);