HcalRawToDigi/bin/makeHcalRaw2DigiGpuValidationPlots.cpp

0001 #include <iomanip>
0002 #include <iostream>
0003 #include <string>
0004 #include <vector>
0005
0006 #include <TCanvas.h>
0007 #include <TFile.h>
0008 #include <TH1D.h>
0009 #include <TH2D.h>
0010 #include <TPaveStats.h>
0011 #include <TTree.h>
0012
0013 #include "CUDADataFormats/HcalDigi/interface/DigiCollection.h"
0014 #include "DataFormats/Common/interface/Wrapper.h"
0015 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0016
0017 #define CREATE_HIST_1D(varname, nbins, first, last) auto varname = new TH1D(#varname, #varname, nbins, first, last)
0018
0019 #define CREATE_HIST_2D(varname, nbins, first, last) \
0020   auto varname = new TH2D(#varname, #varname, nbins, first, last, nbins, first, last)
0021
0022 QIE11DigiCollection filterQIE11(QIE11DigiCollection const& coll) {
0023   QIE11DigiCollection out;
0024   out.reserve(coll.size());
0025
0026   for (uint32_t i = 0; i < coll.size(); i++) {
0027     auto const df = coll[i];
0028     auto const id = HcalDetId{df.id()};
0029     if (id.subdetId() != HcalEndcap)
0030       continue;
0031
0032     out.push_back(QIE11DataFrame{df});
0033   }
0034
0035   return out;
0036 }
0037
0038 int main(int argc, char* argv[]) {
0039   if (argc < 3) {
0040     std::cout << "run with: ./<exe> <path to input file> <path to output file>\n";
0041     exit(0);
0042   }
0043
0044   auto filterf01HE = [](QIE11DigiCollection const& coll) {
0045     QIE11DigiCollection out{coll.samples(), coll.subdetId()};
0046     out.reserve(coll.size());
0047
0048     for (uint32_t i = 0; i < coll.size(); i++) {
0049       auto const df = QIE11DataFrame{coll[i]};
0050       auto const id = HcalDetId{df.id()};
0051       if ((df.flavor() == 0 or df.flavor() == 1) and id.subdetId() == HcalEndcap)
0052         out.push_back(df);
0053     }
0054
0055     return out;
0056   };
0057
0058   auto filterf3HB = [](QIE11DigiCollection const& coll) {
0059     QIE11DigiCollection out{coll.samples(), coll.subdetId()};
0060     out.reserve(coll.size());
0061
0062     for (uint32_t i = 0; i < coll.size(); i++) {
0063       auto const df = QIE11DataFrame{coll[i]};
0064       auto const did = HcalDetId{df.id()};
0065       if (df.flavor() == 3 and did.subdetId() == HcalBarrel)
0066         out.push_back(df);
0067     }
0068
0069     return out;
0070   };
0071
0072   // branches to use
0073   using Collectionf01 =
0074       hcal::DigiCollection<hcal::Flavor1, calo::common::VecStoragePolicy<calo::common::CUDAHostAllocatorAlias>>;
0075   using Collectionf5 =
0076       hcal::DigiCollection<hcal::Flavor5, calo::common::VecStoragePolicy<calo::common::CUDAHostAllocatorAlias>>;
0077   using Collectionf3 =
0078       hcal::DigiCollection<hcal::Flavor3, calo::common::VecStoragePolicy<calo::common::CUDAHostAllocatorAlias>>;
0079   edm::Wrapper<Collectionf01>* wgpuf01he = nullptr;
0080   edm::Wrapper<Collectionf5>* wgpuf5hb = nullptr;
0081   edm::Wrapper<Collectionf3>* wgpuf3hb = nullptr;
0082   edm::Wrapper<QIE11DigiCollection>* wcpuf01he = nullptr;
0083   edm::Wrapper<HBHEDigiCollection>* wcpuf5hb = nullptr;
0084
0085   std::string inFileName{argv[1]};
0086   std::string outFileName{argv[2]};
0087
0088   // prep output
0089   TFile rfout{outFileName.c_str(), "recreate"};
0090
0091   CREATE_HIST_1D(hADCf01HEGPU, 256, 0, 256);
0092   CREATE_HIST_1D(hADCf01HECPU, 256, 0, 256);
0093   CREATE_HIST_1D(hADCf5HBGPU, 128, 0, 128);
0094   CREATE_HIST_1D(hADCf5HBCPU, 128, 0, 128);
0095   CREATE_HIST_1D(hADCf3HBGPU, 256, 0, 256);
0096   CREATE_HIST_1D(hADCf3HBCPU, 256, 0, 256);
0097   CREATE_HIST_1D(hTDCf01HEGPU, 64, 0, 64);
0098   CREATE_HIST_1D(hTDCf01HECPU, 64, 0, 64);
0099
0100   CREATE_HIST_2D(hADCf01HEGPUvsCPU, 256, 0, 256);
0101   CREATE_HIST_2D(hADCf3HBGPUvsCPU, 256, 0, 256);
0102   CREATE_HIST_2D(hADCf5HBGPUvsCPU, 128, 0, 128);
0103   CREATE_HIST_2D(hTDCf01HEGPUvsCPU, 64, 0, 64);
0104   CREATE_HIST_2D(hTDCf3HBGPUvsCPU, 4, 0, 4);
0105
0106   // prep input
0107   TFile rfin{inFileName.c_str()};
0108   TTree* rt = (TTree*)rfin.Get("Events");
0109   rt->SetBranchAddress("QIE11DataFrameHcalDataFrameContainer_hcalDigis__RECO.", &wcpuf01he);
0110   rt->SetBranchAddress("HBHEDataFramesSorted_hcalDigis__RECO.", &wcpuf5hb);
0111   rt->SetBranchAddress(
0112       "hcalFlavor5calocommonCUDAHostAllocatorAliascalocommonVecStoragePolicyhcalDigiCollection_hcalCPUDigisProducer_"
0113       "f5HBDigis_RECO.",
0114       &wgpuf5hb);
0115   rt->SetBranchAddress(
0116       "hcalFlavor1calocommonCUDAHostAllocatorAliascalocommonVecStoragePolicyhcalDigiCollection_hcalCPUDigisProducer_"
0117       "f01HEDigis_RECO.",
0118       &wgpuf01he);
0119   rt->SetBranchAddress(
0120       "hcalFlavor3calocommonCUDAHostAllocatorAliascalocommonVecStoragePolicyhcalDigiCollection_hcalCPUDigisProducer_"
0121       "f3HBDigis_RECO.",
0122       &wgpuf3hb);
0123
0124   // accumulate
0125   auto const nentries = rt->GetEntries();
0126   std::cout << ">>> nentries = " << nentries << std::endl;
0127   for (int ie = 0; ie < nentries; ++ie) {
0128     rt->GetEntry(ie);
0129
0130     auto const& f01HEProduct = wgpuf01he->bareProduct();
0131     auto const& f5HBProduct = wgpuf5hb->bareProduct();
0132     auto const& f3HBProduct = wgpuf3hb->bareProduct();
0133     auto const& qie11Product = wcpuf01he->bareProduct();
0134     auto const qie11Filteredf01 = filterf01HE(qie11Product);
0135     auto const qie11Filteredf3 = filterf3HB(qie11Product);
0136     auto const& qie8Product = wcpuf5hb->bareProduct();
0137
0138     auto const ngpuf01he = f01HEProduct.ids.size();
0139     auto const ngpuf5hb = f5HBProduct.ids.size();
0140     auto const ngpuf3hb = f3HBProduct.ids.size();
0141     auto const ncpuf01he = qie11Filteredf01.size();
0142     auto const ncpuf5hb = qie8Product.size();
0143     auto const ncpuf3hb = qie11Filteredf3.size();
0144
0145     /*
0146         printf("ngpuf01he = %u nqie11 = %u ncpuf01he = %u ngpuf5hb = %u ncpuf5hb = %u\n",
0147             f01HEProduct.size(), qie11Product.size(), qie11Filtered.size(),
0148             f5HBProduct.size(),
0149             static_cast<uint32_t>(qie8Product.size()));
0150             */
0151
0152     if (ngpuf01he != ncpuf01he) {
0153       std::cerr << "*** mismatch in number of flavor 01 digis for event " << ie << std::endl
0154                 << ">>> ngpuf01he = " << ngpuf01he << std::endl
0155                 << ">>> ncpuf01he = " << ncpuf01he << std::endl;
0156     }
0157
0158     {
0159       auto const& idsgpu = f01HEProduct.ids;
0160       auto const& datagpu = f01HEProduct.data;
0161
0162       for (uint32_t ich = 0; ich < ncpuf01he; ich++) {
0163         auto const cpudf = QIE11DataFrame{qie11Filteredf01[ich]};
0164         auto const cpuid = cpudf.id();
0165         auto iter2idgpu = std::find(idsgpu.begin(), idsgpu.end(), cpuid);
0166
0167         if (iter2idgpu == idsgpu.end()) {
0168           std::cerr << "missing " << HcalDetId{cpuid} << std::endl;
0169           continue;
0170         }
0171
0172         // FIXME: cna fail...
0173         assert(*iter2idgpu == cpuid);
0174
0175         auto const ptrdiff = iter2idgpu - idsgpu.begin();
0176         auto const nsamples_gpu = hcal::compute_nsamples<hcal::Flavor1>(f01HEProduct.stride);
0177         auto const nsamples_cpu = qie11Filteredf01.samples();
0178         assert(static_cast<uint32_t>(nsamples_cpu) == nsamples_gpu);
0179
0180         uint32_t ichgpu = ptrdiff;
0181         uint32_t offset = ichgpu * f01HEProduct.stride;
0182         uint16_t const* df_start = datagpu.data() + offset;
0183         for (uint32_t sample = 0u; sample < nsamples_gpu; sample++) {
0184           auto const cpuadc = cpudf[sample].adc();
0185           auto const gpuadc = hcal::adc_for_sample<hcal::Flavor1>(df_start, sample);
0186           auto const cputdc = cpudf[sample].tdc();
0187           auto const gputdc = hcal::tdc_for_sample<hcal::Flavor1>(df_start, sample);
0188           auto const cpucapid = cpudf[sample].capid();
0189           auto const gpucapid = hcal::capid_for_sample<hcal::Flavor1>(df_start, sample);
0190
0191           hADCf01HEGPU->Fill(gpuadc);
0192           hADCf01HECPU->Fill(cpuadc);
0193           hTDCf01HEGPU->Fill(gputdc);
0194           hTDCf01HECPU->Fill(cputdc);
0195           hADCf01HEGPUvsCPU->Fill(cpuadc, gpuadc);
0196           hTDCf01HEGPUvsCPU->Fill(cputdc, gputdc);
0197
0198           // At RAW Decoding level there must not be any mistmatches
0199           // in the adc values at all!
0200           assert(static_cast<uint8_t>(cpuadc) == gpuadc);
0201           assert(static_cast<uint8_t>(cputdc) == gputdc);
0202           assert(static_cast<uint8_t>(cpucapid) == gpucapid);
0203         }
0204       }
0205     }
0206
0207     if (ngpuf3hb != ncpuf3hb) {
0208       std::cerr << "*** mismatch in number of flavor 3 digis for event " << ie << std::endl
0209                 << ">>> ngpuf01he = " << ngpuf3hb << std::endl
0210                 << ">>> ncpuf01he = " << ncpuf3hb << std::endl;
0211     }
0212
0213     {
0214       auto const& idsgpu = f3HBProduct.ids;
0215       auto const& datagpu = f3HBProduct.data;
0216
0217       for (uint32_t ich = 0; ich < ncpuf3hb; ich++) {
0218         auto const cpudf = QIE11DataFrame{qie11Filteredf3[ich]};
0219         auto const cpuid = cpudf.id();
0220         auto iter2idgpu = std::find(idsgpu.begin(), idsgpu.end(), cpuid);
0221
0222         if (iter2idgpu == idsgpu.end()) {
0223           std::cerr << "missing " << HcalDetId{cpuid} << std::endl;
0224           continue;
0225         }
0226
0227         // FIXME: cna fail...
0228         assert(*iter2idgpu == cpuid);
0229
0230         auto const ptrdiff = iter2idgpu - idsgpu.begin();
0231         auto const nsamples_gpu = hcal::compute_nsamples<hcal::Flavor3>(f3HBProduct.stride);
0232         auto const nsamples_cpu = qie11Filteredf3.samples();
0233         assert(static_cast<uint32_t>(nsamples_cpu) == nsamples_gpu);
0234
0235         uint32_t ichgpu = ptrdiff;
0236         uint32_t offset = ichgpu * f3HBProduct.stride;
0237         uint16_t const* df_start = datagpu.data() + offset;
0238         for (uint32_t sample = 0u; sample < nsamples_gpu; sample++) {
0239           auto const cpuadc = cpudf[sample].adc();
0240           auto const gpuadc = hcal::adc_for_sample<hcal::Flavor3>(df_start, sample);
0241           auto const cputdc = cpudf[sample].tdc();
0242           auto const gputdc = hcal::tdc_for_sample<hcal::Flavor3>(df_start, sample);
0243
0244           hADCf3HBGPU->Fill(gpuadc);
0245           hADCf3HBCPU->Fill(cpuadc);
0246           hADCf3HBGPUvsCPU->Fill(cpuadc, gpuadc);
0247           hTDCf3HBGPUvsCPU->Fill(cputdc, gputdc);
0248
0249           // At RAW Decoding level there must not be any mistmatches
0250           // in the adc values at all!
0251           assert(static_cast<uint8_t>(cpuadc) == gpuadc);
0252           assert(static_cast<uint8_t>(cputdc) == gputdc);
0253         }
0254       }
0255     }
0256
0257     if (ngpuf5hb != ncpuf5hb) {
0258       std::cerr << "*** mismatch in number of flavor 5 digis for event " << ie << std::endl
0259                 << ">>> ngpuf5hb = " << ngpuf5hb << std::endl
0260                 << ">>> ncpuf5hb = " << ncpuf5hb << std::endl;
0261     }
0262
0263     {
0264       auto const& idsgpu = f5HBProduct.ids;
0265       auto const& datagpu = f5HBProduct.data;
0266       for (uint32_t i = 0; i < ncpuf5hb; i++) {
0267         auto const cpudf = qie8Product[i];
0268         auto const cpuid = cpudf.id().rawId();
0269         auto iter2idgpu = std::find(idsgpu.begin(), idsgpu.end(), cpuid);
0270         if (iter2idgpu == idsgpu.end()) {
0271           std::cerr << "missing " << HcalDetId{cpuid} << std::endl;
0272           continue;
0273         }
0274
0275         assert(*iter2idgpu == cpuid);
0276
0277         auto const ptrdiff = iter2idgpu - idsgpu.begin();
0278         auto const nsamples_gpu = hcal::compute_nsamples<hcal::Flavor5>(f5HBProduct.stride);
0279         auto const nsamples_cpu = qie8Product[0].size();
0280         assert(static_cast<uint32_t>(nsamples_cpu) == nsamples_gpu);
0281
0282         uint32_t offset = ptrdiff * f5HBProduct.stride;
0283         uint16_t const* df_start = datagpu.data() + offset;
0284         for (uint32_t sample = 0u; sample < nsamples_gpu; sample++) {
0285           auto const cpuadc = cpudf.sample(sample).adc();
0286           auto const gpuadc = hcal::adc_for_sample<hcal::Flavor5>(df_start, sample);
0287           auto const cpucapid = cpudf.sample(sample).capid();
0288           auto const gpucapid = hcal::capid_for_sample<hcal::Flavor1>(df_start, sample);
0289
0290           hADCf5HBGPU->Fill(gpuadc);
0291           hADCf5HBCPU->Fill(cpuadc);
0292           hADCf5HBGPUvsCPU->Fill(cpuadc, gpuadc);
0293
0294           // the must for us at RAW Decoding stage
0295           assert(static_cast<uint8_t>(cpuadc) == gpuadc);
0296           assert(static_cast<uint8_t>(cpucapid) == gpucapid);
0297         }
0298       }
0299     }
0300   }
0301
0302   {
0303     TCanvas c{"plots", "plots", 4200, 6200};
0304     c.Divide(3, 3);
0305     c.cd(1);
0306     {
0307       gPad->SetLogy();
0308       hADCf01HECPU->SetLineColor(kBlack);
0309       hADCf01HECPU->SetLineWidth(1.);
0310       hADCf01HECPU->Draw("");
0311       hADCf01HEGPU->SetLineColor(kBlue);
0312       hADCf01HEGPU->SetLineWidth(1.);
0313       hADCf01HEGPU->Draw("sames");
0314       gPad->Update();
0315       auto stats = (TPaveStats*)hADCf01HEGPU->FindObject("stats");
0316       auto y2 = stats->GetY2NDC();
0317       auto y1 = stats->GetY1NDC();
0318       stats->SetY2NDC(y1);
0319       stats->SetY1NDC(y1 - (y2 - y1));
0320     }
0321     c.cd(2);
0322     {
0323       gPad->SetLogy();
0324       hADCf5HBCPU->SetLineColor(kBlack);
0325       hADCf5HBCPU->SetLineWidth(1.);
0326       hADCf5HBCPU->Draw("");
0327       hADCf5HBGPU->SetLineColor(kBlue);
0328       hADCf5HBGPU->SetLineWidth(1.);
0329       hADCf5HBGPU->Draw("sames");
0330       gPad->Update();
0331       auto stats = (TPaveStats*)hADCf5HBGPU->FindObject("stats");
0332       auto y2 = stats->GetY2NDC();
0333       auto y1 = stats->GetY1NDC();
0334       stats->SetY2NDC(y1);
0335       stats->SetY1NDC(y1 - (y2 - y1));
0336     }
0337     c.cd(3);
0338     {
0339       gPad->SetLogy();
0340       hADCf3HBCPU->SetLineColor(kBlack);
0341       hADCf3HBCPU->SetLineWidth(1.);
0342       hADCf3HBCPU->Draw("");
0343       hADCf3HBGPU->SetLineColor(kBlue);
0344       hADCf3HBGPU->SetLineWidth(1.);
0345       hADCf3HBGPU->Draw("sames");
0346       gPad->Update();
0347       auto stats = (TPaveStats*)hADCf3HBGPU->FindObject("stats");
0348       auto y2 = stats->GetY2NDC();
0349       auto y1 = stats->GetY1NDC();
0350       stats->SetY2NDC(y1);
0351       stats->SetY1NDC(y1 - (y2 - y1));
0352     }
0353     c.cd(4);
0354     hADCf01HEGPUvsCPU->Draw("colz");
0355     c.cd(5);
0356     hADCf5HBGPUvsCPU->Draw("colz");
0357     c.cd(6);
0358     hADCf3HBGPUvsCPU->Draw("colz");
0359     c.cd(7);
0360     {
0361       gPad->SetLogy();
0362       hTDCf01HECPU->SetLineColor(kBlack);
0363       hTDCf01HECPU->SetLineWidth(1.);
0364       hTDCf01HECPU->Draw("");
0365       hTDCf01HEGPU->SetLineColor(kBlue);
0366       hTDCf01HEGPU->SetLineWidth(1.);
0367       hTDCf01HEGPU->Draw("sames");
0368       gPad->Update();
0369       auto stats = (TPaveStats*)hTDCf01HEGPU->FindObject("stats");
0370       auto y2 = stats->GetY2NDC();
0371       auto y1 = stats->GetY1NDC();
0372       stats->SetY2NDC(y1);
0373       stats->SetY1NDC(y1 - (y2 - y1));
0374     }
0375     c.cd(8);
0376     hTDCf01HEGPUvsCPU->Draw("colz");
0377     c.cd(9);
0378     hTDCf3HBGPUvsCPU->Draw("colz");
0379
0380     c.SaveAs("plots.pdf");
0381   }
0382
0383   rfin.Close();
0384   rfout.Write();
0385   rfout.Close();
0386 }