HcalRecProducers/src/HCALGPUAnalyzer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    ComparisonPlots/HCALGPUAnalyzer
0004 // Class:      HCALGPUAnalyzer
0005 //
0006 /**\class HCALGPUAnalyzer HCALGPUAnalyzer.cc ComparisonPlots/HCALGPUAnalyzer/plugins/HCALGPUAnalyzer.cc
0007
0008  Description: [one line class summary]
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Mariarosaria D'Alfonso
0015 //         Created:  Mon, 17 Dec 2018 16:22:58 GMT
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021 #include <string>
0022 #include <map>
0023 #include <iostream>
0024 using namespace std;
0025
0026 // user include files
0027 #include "FWCore/Framework/interface/Frameworkfwd.h"
0028 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0029
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032
0033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0034 #include "FWCore/Utilities/interface/InputTag.h"
0035 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0036 #include "FWCore/ServiceRegistry/interface/Service.h"
0037 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0038
0039 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
0040 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0041 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0042
0043 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0044 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0045
0046 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSimParameterMap.h"
0047
0048 #include "TH2F.h"
0049
0050 //
0051 // class declaration
0052 //
0053
0054 class HCALGPUAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0055 public:
0056   explicit HCALGPUAnalyzer(const edm::ParameterSet &);
0057   ~HCALGPUAnalyzer() override = default;
0058
0059   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0060
0061 private:
0062   void beginJob() override;
0063   void analyze(const edm::Event &, const edm::EventSetup &) override;
0064   void endJob() override;
0065
0066   // ----------member data ---------------------------
0067   //  void ClearVariables();
0068
0069   // some variables for storing information
0070   double Method0Energy, Method0EnergyGPU;
0071   double RecHitEnergy, RecHitEnergyGPU;
0072   double RecHitTime, RecHitTimeGPU;
0073   double iEta, iEtaGPU;
0074   double iPhi, iPhiGPU;
0075   int depth, depthGPU;
0076
0077   TH2F *hEnergy_2dMahi;
0078   TH2F *hEnergy_2dM0;
0079   TH2F *hTime_2dMahi;
0080
0081   TH2F *Unmatched;
0082   TH2F *Matched;
0083   TH1F *hEnergy_cpu;
0084   TH1F *hEnergy_gpu;
0085   TH1F *hEnergy_cpugpu;
0086   TH1F *hEnergy_cpugpu_rel;
0087   TH1F *hEnergyM0_cpu;
0088   TH1F *hEnergyM0_gpu;
0089   TH1F *hTime_cpu;
0090   TH1F *hTime_gpu;
0091
0092   // create the output file
0093   edm::Service<TFileService> FileService;
0094   // create the token to retrieve hit information
0095   edm::EDGetTokenT<HBHERecHitCollection> hRhToken;
0096   edm::EDGetTokenT<HBHERecHitCollection> hRhTokenGPU;
0097 };
0098
0099 //
0100 // constants, enums and typedefs
0101 //
0102
0103 //
0104 // static data member definitions
0105 //
0106
0107 //
0108 // constructors and destructor
0109 //
0110 HCALGPUAnalyzer::HCALGPUAnalyzer(const edm::ParameterSet &iConfig) {
0111   usesResource("TFileService");
0112
0113   hRhToken = consumes<HBHERecHitCollection>(iConfig.getUntrackedParameter<string>("HBHERecHits", "hbheprereco"));
0114   hRhTokenGPU = consumes<HBHERecHitCollection>(
0115       iConfig.getUntrackedParameter<string>("HBHERecHits", "hcalCPURecHitsProducer:recHitsLegacyHBHE"));
0116
0117   //
0118
0119   hEnergy_2dM0 = FileService->make<TH2F>("hEnergy_2dM0", "hEnergy_2dM0", 1000, 0., 100., 1000, 0., 100.);
0120   hEnergy_2dM0->GetXaxis()->SetTitle("Cpu M0 Energy");
0121   hEnergy_2dM0->GetYaxis()->SetTitle("GPU M0 Energy");
0122
0123   hEnergy_2dMahi = FileService->make<TH2F>("hEnergy_2dMahi", "hEnergy_2dMahi", 1000, 0., 100., 1000, 0., 100.);
0124   hEnergy_2dMahi->GetXaxis()->SetTitle("CPU Energy");
0125   hEnergy_2dMahi->GetYaxis()->SetTitle("GPU Energy");
0126
0127   hTime_2dMahi = FileService->make<TH2F>("hTime_2dMahi", "hTime_2dMahi", 250, -12.5, 12.5, 250, -12.5, 12.5);
0128   hTime_2dMahi->GetXaxis()->SetTitle("Mahi Time CPU");
0129   hTime_2dMahi->GetYaxis()->SetTitle("Mahi Time GPU");
0130
0131   //
0132
0133   hEnergyM0_cpu = FileService->make<TH1F>("hEnergyM0_cpu", "hEnergyM0_cpu", 100, 0., 100.);
0134   hEnergyM0_cpu->GetXaxis()->SetTitle("CPU Energy");
0135
0136   hEnergy_cpu = FileService->make<TH1F>("hEnergy_cpu", "hEnergy_cpu", 50, 0., 50.);
0137   hEnergy_cpu->GetXaxis()->SetTitle("CPU Energy");
0138
0139   hEnergy_gpu = FileService->make<TH1F>("hEnergy_gpu", "hEnergy_gpu", 50, 0., 50.);
0140   hEnergy_gpu->GetXaxis()->SetTitle("GPU Energy");
0141
0142   //
0143
0144   hEnergy_cpugpu = FileService->make<TH1F>("hEnergy_cpugpu", "hEnergy_cpugpu", 500, -2.5, 2.5);
0145   hEnergy_cpugpu->GetXaxis()->SetTitle("GPU Energy - CPU Energy [GeV]");
0146   hEnergy_cpugpu->GetYaxis()->SetTitle("# RecHits");
0147
0148   hEnergy_cpugpu_rel =
0149       FileService->make<TH1F>("hEnergy_cpugpu_rel", "hEnergy_cpugpu_rel ( E > 0.005 GeV)", 500, -2.5, 2.5);
0150   hEnergy_cpugpu_rel->GetXaxis()->SetTitle("(GPU Energy - CPU Energy) / CPU energy");
0151   hEnergy_cpugpu_rel->GetYaxis()->SetTitle("# RecHits");
0152
0153   //
0154
0155   hTime_cpu = FileService->make<TH1F>("hTime_cpu", "hTime_cpu", 50, -25., 25.);
0156   hTime_cpu->GetXaxis()->SetTitle("CPU Time");
0157
0158   hTime_gpu = FileService->make<TH1F>("hTime_gpu", "hTime_gpu", 50, -25., 25.);
0159   hTime_gpu->GetXaxis()->SetTitle("GPU Time");
0160
0161   Unmatched = FileService->make<TH2F>("Unmatched", "Unmatched (eta,phi)", 100, -50., 50., 85, 0., 85.);
0162   Matched = FileService->make<TH2F>("Matched", "Matched (eta,phi)", 100, -50., 50., 85, 0., 85.);
0163
0164   //now do what ever initialization is needed
0165 }
0166
0167 //
0168 // member functions
0169 //
0170
0171 // ------------ method called for each event  ------------
0172 void HCALGPUAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0173   using namespace edm;
0174
0175   // Read events
0176   Handle<HBHERecHitCollection> hRecHits;
0177   iEvent.getByToken(hRhToken, hRecHits);
0178
0179   Handle<HBHERecHitCollection> hRecHitsGPU;
0180   iEvent.getByToken(hRhTokenGPU, hRecHitsGPU);
0181
0182   // Loop over all rechits in one event
0183   for (int i = 0; i < (int)hRecHits->size(); i++) {
0184     // get ID information for the reconstructed hit
0185     HcalDetId detID_rh = (*hRecHits)[i].id().rawId();
0186
0187     // ID information can get us detector coordinates
0188     depth = (*hRecHits)[i].id().depth();
0189     iEta = detID_rh.ieta();
0190     iPhi = detID_rh.iphi();
0191
0192     // get some variables
0193     Method0Energy = (*hRecHits)[i].eraw();
0194     RecHitEnergy = (*hRecHits)[i].energy();
0195     RecHitTime = (*hRecHits)[i].time();
0196
0197     hEnergy_cpu->Fill(RecHitEnergy);
0198     hTime_cpu->Fill(RecHitTime);
0199
0200     /*
0201      cout << "Run " << i << ": ";
0202      cout << "Method0Energy: " << Method0Energy;
0203      cout << "RecHitEnergy: " << RecHitEnergy;
0204      cout << "depth: " << depth;
0205      cout << "iEta: " << iEta;
0206      cout << "iPhi: " << iPhi;
0207      cout << "RecHitTime" << RecHitTime;
0208      */
0209   }
0210
0211   for (int i = 0; i < (int)hRecHitsGPU->size(); i++) {
0212     // get ID information for the reconstructed hit
0213     HcalDetId detID_rh = (*hRecHitsGPU)[i].id().rawId();
0214
0215     // ID information can get us detector coordinates
0216     depthGPU = (*hRecHitsGPU)[i].id().depth();
0217     iEtaGPU = detID_rh.ieta();
0218     iPhiGPU = detID_rh.iphi();
0219
0220     // get some variables
0221     Method0EnergyGPU = (*hRecHitsGPU)[i].eraw();
0222     RecHitEnergyGPU = (*hRecHitsGPU)[i].energy();
0223     RecHitTimeGPU = (*hRecHitsGPU)[i].time();
0224
0225     hEnergy_gpu->Fill(RecHitEnergyGPU);
0226     hTime_gpu->Fill(RecHitTimeGPU);
0227
0228     /*
0229      cout << "Run " << i << ": ";
0230      cout << "Method0Energy: " << Method0EnergyGPU;
0231      cout << "RecHitEnergy: " << RecHitEnergyGPU;
0232      cout << "depth: " << depthGPU;
0233      cout << "iEta: " << iEtaGPU;
0234      cout << "iPhi: " << iPhiGPU;
0235      cout << "RecHitTime" << RecHitTimeGPU;
0236      */
0237   }
0238
0239   // Loop over all rechits in one event
0240   for (int i = 0; i < (int)hRecHits->size(); i++) {
0241     HcalDetId detID_rh = (*hRecHits)[i].id().rawId();
0242
0243     bool unmatched = true;
0244     //     cout << "--------------------------------------------------------" << endl;
0245
0246     for (int j = 0; j < (int)hRecHitsGPU->size(); j++) {
0247       HcalDetId detID_gpu = (*hRecHitsGPU)[j].id().rawId();
0248
0249       if ((detID_rh == detID_gpu)) {
0250         /*
0251      cout << "Mtime(cpu)" << (*hRecHits)[i].time() << endl;
0252      cout << "     Mtime(gpu)" << (*hRecHitsGPU)[j].time() << endl;
0253
0254      cout << "M0E(cpu)" << (*hRecHits)[i].eraw() << endl;
0255      cout << "     M0E(gpu)" << (*hRecHitsGPU)[j].eraw() << endl;
0256      */
0257
0258         auto relValue = ((*hRecHitsGPU)[j].energy() - (*hRecHits)[i].energy()) / (*hRecHits)[i].energy();
0259
0260         hEnergy_2dM0->Fill((*hRecHits)[i].eraw(), (*hRecHitsGPU)[j].eraw());
0261         hEnergy_2dMahi->Fill((*hRecHits)[i].energy(), (*hRecHitsGPU)[j].energy());
0262         hEnergy_cpugpu->Fill((*hRecHitsGPU)[j].energy() - (*hRecHits)[i].energy());
0263         if ((*hRecHits)[i].energy() > 0.005)
0264           hEnergy_cpugpu_rel->Fill(relValue);
0265         hTime_2dMahi->Fill((*hRecHits)[i].time(), (*hRecHitsGPU)[j].time());
0266
0267         /*
0268      if((relValue < - 0.9) and ((*hRecHits)[i].energy()>0.005)) {
0269        cout << "----------------------------------"<< endl;
0270        cout << " detID = " << detID_rh.rawId() << endl;
0271        cout << "ME(cpu)" << (*hRecHits)[i].energy() << endl;
0272        cout << "     ME(gpu)" << (*hRecHitsGPU)[j].energy() << endl;
0273      }
0274      */
0275
0276         Matched->Fill(detID_rh.ieta(), detID_rh.iphi());
0277
0278         unmatched = false;
0279       }
0280     }
0281
0282     ///
0283
0284     if (unmatched) {
0285       Unmatched->Fill(detID_rh.ieta(), detID_rh.iphi());
0286       //       cout << "   recHit not matched ="  << detID_rh << "  E(raw)=" << (*hRecHits)[i].eraw() << " E=" << (*hRecHits)[i].energy() << endl;
0287     }
0288   }
0289 }
0290
0291 // ------------ method called once each job just before starting event loop  ------------
0292 void HCALGPUAnalyzer::beginJob() {}
0293
0294 // ------------ method called once each job just after ending the event loop  ------------
0295 void HCALGPUAnalyzer::endJob() {}
0296
0297 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0298 void HCALGPUAnalyzer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0299   //The following says we do not know what parameters are allowed so do no validation
0300   // Please change this to state exactly what you do use, even if it is no parameters
0301   edm::ParameterSetDescription desc;
0302   desc.setUnknown();
0303   descriptions.addDefault(desc);
0304 }
0305
0306 //define this as a plug-in
0307 DEFINE_FWK_MODULE(HCALGPUAnalyzer);