benchmarks/zdc/simple_info_plot_histograms.cxx

0001 R__LOAD_LIBRARY(libGenDetectors.so)
0002 R__LOAD_LIBRARY(libfmt.so)
0003 #include "fmt/core.h"
0004 R__LOAD_LIBRARY(libDDG4IO.so)
0005 #include "DD4hep/Detector.h"
0006 #include "DDG4/Geant4Data.h"
0007 #include "DDRec/CellIDPositionConverter.h"
0008 #include "DDRec/SurfaceManager.h"
0009 #include "DD4hep/VolumeManager.h"
0010 #include "DD4hep/detail/Handle.inl"
0011 #include "DD4hep/detail/ObjectsInterna.h"
0012 #include "DD4hep/detail/DetectorInterna.h"
0013 #include "DD4hep/detail/VolumeManagerInterna.h"
0014 #include "DDRec/Surface.h"
0015 #include "DD4hep/Volumes.h"
0016 #include "DD4hep/DetElement.h"
0017 #include "DD4hep/NamedObject.h"
0018 #include "DD4hep/IDDescriptor.h"
0019 #include "DD4hep/ConditionsMap.h"
0020 #include "TGeoMatrix.h"
0021 #include "ROOT/RDataFrame.hxx"
0022 #include "Math/Vector3D.h"
0023 #include "Math/Vector4D.h"
0024 #include "Math/VectorUtil.h"
0025 #include "TCanvas.h"
0026 #include "TLegend.h"
0027 #include "TMath.h"
0028 #include "TRandom3.h"
0029 #include "TFile.h"
0030 #include "TH1F.h"
0031 #include "TH1D.h"
0032 #include "TTree.h"
0033 #include "TChain.h"
0034 #include <random>
0035 #include <iostream>
0036 #include "TStyle.h"
0037
0038 void simple_info_plot_histograms(const char* fname = "sim_output/output_zdc_photons.edm4hep.root"){
0039
0040   // Setting for graphs
0041   gROOT->SetStyle("Plain");
0042   gStyle->SetOptFit(1);
0043   gStyle->SetLineWidth(2);
0044   gStyle->SetPadTickX(1);
0045   gStyle->SetPadTickY(1);
0046   gStyle->SetPadGridX(1);
0047   gStyle->SetPadGridY(1);
0048   gStyle->SetPadLeftMargin(0.14);
0049
0050   TChain* t = new TChain("events");
0051   t->Add(fname);
0052
0053   ROOT::RDataFrame d0(*t);//, {"ZDCHits","MCParticles"});
0054
0055   // Detector
0056   dd4hep::Detector& detector = dd4hep::Detector::getInstance();
0057   detector.fromCompact("benchmarks/zdc/ZDC_example.xml");
0058   // Volume
0059   dd4hep::VolumeManager volman = dd4hep::VolumeManager::getVolumeManager(detector);
0060   // CellID Coverter
0061   dd4hep::rec::CellIDPositionConverter cellid_converter(detector);
0062
0063   // Number of hits
0064   auto nhits = [] (std::vector<dd4hep::sim::Geant4Calorimeter::Hit*>& hits){ return (int) hits.size(); };
0065
0066   // Cell ID
0067   auto cellID = [&] (const std::vector<dd4hep::sim::Geant4Calorimeter::Hit*>& hits) {
0068         std::vector<double> result;
0069         for(const auto& h: hits)
0070                 result.push_back(h->cellID);
0071   return result;
0072   };
0073
0074   // Volume ID
0075   auto volID = [&] (const std::vector<dd4hep::sim::Geant4Calorimeter::Hit*>& hits) {
0076         std::vector<double> result;
0077         for(const auto& h: hits) {
0078         // method 1: use cell ID to get volume ID
0079                 auto volcontext = cellid_converter.findContext(h->cellID);
0080         //auto volid = volcontext->identifier;
0081
0082         // method 2: use detector element, readout, segmentation, then volume ID
0083         //dd4hep::Readout r = cellid_converter.findReadout(volcontext->element);
0084         //dd4hep::Segmentation seg = r.segmentation();
0085         //auto volid = seg.volumeID(h->cellID);
0086
0087          result.push_back(volcontext->identifier);
0088         }
0089   return result;
0090   };
0091
0092   // Detector ID
0093   auto detID = [&] (const std::vector<dd4hep::sim::Geant4Calorimeter::Hit*>& hits) {
0094         std::vector<double> result;
0095         for(const auto& h: hits) {
0096         auto detelement = volman.lookupDetector(h->cellID);
0097         result.push_back(detelement.volumeID());
0098     }
0099   return result;
0100   };
0101
0102   // Hit position X
0103   auto hit_x_position = [&] (const std::vector<dd4hep::sim::Geant4Calorimeter::Hit*>& hits) {
0104     std::vector<double> result;
0105     for(const auto& h: hits)
0106         result.push_back(h->position.x()); //mm
0107   return result;
0108   };
0109
0110   // Hit position Y
0111   auto hit_y_position = [&] (const std::vector<dd4hep::sim::Geant4Calorimeter::Hit*>& hits) {
0112         std::vector<double> result;
0113         for(const auto& h: hits)
0114         result.push_back(h->position.y()); //mm
0115   return result;
0116   };
0117
0118   // Hit position Z
0119   auto hit_z_position = [&] (const std::vector<dd4hep::sim::Geant4Calorimeter::Hit*>& hits) {
0120         std::vector<double> result;
0121         for(const auto& h: hits)
0122         result.push_back(h->position.z()); //mm
0123   return result;
0124   };
0125
0126   // Energy deposition
0127   auto e_dep = [&] (const std::vector<dd4hep::sim::Geant4Calorimeter::Hit*>& hits) {
0128         std::vector<double> result;
0129         for(const auto& h: hits)
0130         result.push_back(h->energyDeposit); //GeV
0131   return result;
0132   };
0133
0134   auto d1 = d0.Define("nhits", nhits, {"ZDCHits"})
0135           .Define("cellID", cellID, {"ZDCHits"})
0136           .Define("volID", volID, {"ZDCHits"})
0137           .Define("detID", detID, {"ZDCHits"})
0138               .Define("hit_x_position", hit_x_position, {"ZDCHits"})
0139           .Define("hit_y_position", hit_y_position, {"ZDCHits"})
0140               .Define("hit_z_position", hit_z_position, {"ZDCHits"})
0141           .Define("e_dep", e_dep, {"ZDCHits"})
0142           ;
0143
0144   // Define Histograms
0145   auto h0 = d1.Histo1D({"h0", "nhits histogram; nhits; Events", 100, 0,5000}, "nhits");
0146   auto h1 = d1.Histo1D({"h1", "hit position X histogram; hit position X [mm]; Events", 60,-30,30}, "hit_x_position");
0147   auto h2 = d1.Histo1D({"h2", "hit position Y histogram; hit position Y [mm]; Events", 100,-30,80}, "hit_y_position");
0148   auto h3 = d1.Histo1D({"h3", "hit position Z histogram; hit position Z [mm]; Events", 100,1000,1300}, "hit_z_position");
0149   auto h4 = d1.Histo1D({"h4", "energy deposition histogram; energy deposition [GeV]; Events", 100,0,300}, "e_dep");
0150   auto h5 = d1.Histo1D({"h5", "detector ID; detector ID; Events", 3,-0.5,2.5}, "detID");
0151   auto h6 = d1.Histo1D({"h6", "volume ID; volume ID; Events", 100,0,50000000}, "volID");
0152
0153   auto h7 = d1.Histo2D({"h7", "hit position Y vs. X histogram; hit position X [mm]; hit position Y [mm]", 100,-30,30, 100, -30, 30}, "hit_x_position", "hit_y_position");
0154
0155   auto n0 = d1.Filter([](int n){ return (n>0); },{"nhits"}).Count();
0156
0157   d1.Snapshot("info_EVENT","sim_output/info_zdc_photons.edm4hep.root");
0158   std::cout << *n0 << " events with nonzero hits\n";
0159
0160   TCanvas *c1 = new TCanvas("c1","c1",800,600);
0161   c1->SetLogy(1);
0162   h0->GetYaxis()->SetTitleOffset(1.4);
0163   h0->SetLineWidth(2);
0164   h0->SetLineColor(kBlack);
0165   h0->DrawClone();
0166   c1->SaveAs("sim_output/nhits_histo_zdc_photons.png");
0167
0168   TCanvas *c2 = new TCanvas("c2","c2",1000,1000);
0169   c2->Divide(2,2);
0170   c2->cd(1);
0171   h1->GetYaxis()->SetTitleOffset(1.7);
0172   h1->SetLineWidth(2);
0173   h1->SetLineColor(kBlack);
0174   h1->DrawClone();
0175
0176   c2->cd(2);
0177   h2->GetYaxis()->SetTitleOffset(1.7);
0178   h2->SetLineWidth(2);
0179   h2->SetLineColor(kBlack);
0180   h2->DrawClone();
0181
0182   c2->cd(3);
0183   h3->GetYaxis()->SetTitleOffset(1.7);
0184   h3->SetLineWidth(2);
0185   h3->SetLineColor(kBlack);
0186   h3->DrawClone();
0187   c2->SaveAs("sim_output/hit_postion_histo_zdc_photons.png");
0188
0189   TCanvas *c3 = new TCanvas("c3","c3",600,600);
0190   c3->cd();
0191   c3->SetLogy(1);
0192   h4->GetYaxis()->SetTitleOffset(1.4);
0193   h4->SetLineWidth(2);
0194   h4->SetLineColor(kBlack);
0195   h4->DrawClone();
0196   c3->SaveAs("sim_output/edep_histo_zdc_photons.png");
0197
0198   TCanvas *c4 = new TCanvas("c4","c4",1000,600);
0199   c4->Divide(2,1);
0200   c4->SetLogy(0);
0201   c4->cd(1);
0202   h5->GetYaxis()->SetTitleOffset(2.0);
0203   h5->SetLineWidth(2);
0204   h5->SetLineColor(kBlack);
0205   h5->DrawClone();
0206
0207   c4->cd(2);
0208   h6->GetYaxis()->SetTitleOffset(2.0);
0209   h6->SetLineWidth(2);
0210   h6->SetLineColor(kBlack);
0211   h6->DrawClone();
0212   c4->SaveAs("sim_output/detID_volID_histo_zdc_photons.png");
0213
0214   TCanvas *c5 = new TCanvas("c5","c5", 500, 500);
0215   h7->Draw("COLZ");
0216   c5->SaveAs("results/far_forward/zdc/zdc_x_y_hit_map.png");
0217
0218   if(*n0<1) {
0219     std::quick_exit(1);
0220   }
0221 }
0222