drich-dev/src/test_pixel_gap_cuts.cpp

0001 // test if the pixel gap cuts are working
0002
0003 #include <spdlog/spdlog.h>
0004
0005 #include <TApplication.h>
0006 #include <TH2D.h>
0007 #include <TCanvas.h>
0008 #include <TStyle.h>
0009
0010 #include <podio/ROOTFrameReader.h>
0011 #include <podio/Frame.h>
0012 #include "DDRec/CellIDPositionConverter.h"
0013
0014 #include <edm4eic/MCRecoTrackerHitAssociationCollection.h>
0015
0016 #include <services/geometry/richgeo/ReadoutGeo.h>
0017
0018 int main(int argc, char** argv) {
0019
0020   // args
0021   if(argc<=1) {
0022     fmt::print("USAGE: {} [n/i] [file(default=out/rec.edm4hep.root)]\n", argv[0]);
0023     fmt::print(" [n/i]: n=non-interactive, i=interactive\n");
0024     return 2;
0025   }
0026   std::string interactiveOpt = std::string(argv[1]);
0027   std::string root_file_name = argc>2 ? std::string(argv[2]) : "out/rec.edm4hep.root";
0028
0029   // set interactive mode
0030   bool interactiveOn = interactiveOpt=="i";
0031   TApplication *app;
0032   if(interactiveOn)
0033     app = new TApplication("app", &argc, argv);
0034
0035   // logger
0036   auto logger = spdlog::default_logger()->clone("richgeo");
0037   logger->set_level(spdlog::level::trace);
0038
0039   // geometry from compact file
0040   std::string DETECTOR_PATH(getenv("DETECTOR_PATH"));
0041   std::string DETECTOR_CONFIG(getenv("DETECTOR_CONFIG"));
0042   if(DETECTOR_PATH.empty() || DETECTOR_CONFIG.empty()) {
0043     logger->error("cannot find default compact file, since env vars DETECTOR_PATH and DETECTOR_CONFIG are not set");
0044     return 1;
0045   }
0046   auto compactFile = DETECTOR_PATH + "/" + DETECTOR_CONFIG + ".xml";
0047   dd4hep::Detector & det = dd4hep::Detector::getInstance();
0048   det.fromXML(compactFile);
0049   dd4hep::rec::CellIDPositionConverter cellid_converter(det);
0050
0051   // ReadoutGeo
0052   richgeo::ReadoutGeo geo("DRICH", &det, &cellid_converter, logger);
0053
0054   // open input file
0055   auto reader = podio::ROOTFrameReader();
0056   reader.openFile(root_file_name);
0057   const std::string tree_name = "events";
0058
0059   // local hits histogram
0060   auto h = new TH2D("h","local MC SiPM hits",10000,-15,15,10000,-15,15);
0061
0062   // event loop
0063   for(unsigned e=0; e<reader.getEntries(tree_name); e++) {
0064     logger->trace("EVENT {}", e);
0065     auto frame = podio::Frame(reader.readNextEntry(tree_name));
0066     const auto& hit_assocs  = frame.get<edm4eic::MCRecoTrackerHitAssociationCollection>("DRICHRawHitsAssociations");
0067     if(!hit_assocs.isValid())
0068       throw std::runtime_error("cannot find hit associations");
0069
0070     for(const auto& hit_assoc : hit_assocs) {
0071       for(const auto& sim_hit : hit_assoc.getSimHits()) {
0072
0073         auto cellID = sim_hit.getCellID();
0074         auto pos    = sim_hit.getPosition();
0075         dd4hep::Position pos_global(pos.x*dd4hep::mm, pos.y*dd4hep::mm, pos.z*dd4hep::mm);
0076         auto pos_local = geo.GetSensorLocalPosition(cellID, pos_global);
0077         h->Fill(pos_local.y()/dd4hep::mm, pos_local.x()/dd4hep::mm);
0078
0079         if(logger->level() <= spdlog::level::trace) {
0080           logger->trace("pixel hit on cellID={:#018x}",cellID);
0081           auto print_pos = [&] (std::string name, dd4hep::Position p) {
0082             logger->trace("  {:>30} x={:.2f} y={:.2f} z={:.2f} [mm]: ", name, p.x()/dd4hep::mm,  p.y()/dd4hep::mm,  p.z()/dd4hep::mm);
0083           };
0084           print_pos("pos_local",  pos_local);
0085           print_pos("pos_global", pos_global);
0086           // auto dim = m_cellid_converter->cellDimensions(cellID);
0087           // for (std::size_t j = 0; j < std::size(dim); ++j)
0088           //   logger->trace("   - dimension {:<5} size: {:.2}",  j, dim[j]);
0089         }
0090       }
0091     }
0092   }
0093
0094   gStyle->SetOptStat(0);
0095   auto c = new TCanvas();
0096   h->Draw();
0097   fmt::print("NUMBER OF DIGITIZED PHOTONS: {}\n", h->GetEntries());
0098   if(interactiveOn) {
0099     fmt::print("\n\npress ^C to exit.\n\n");
0100     app->Run();
0101   } else {
0102     c->SaveAs("out/pixel_gaps.png");
0103   }
0104 }