benchmarks/zdc/simple_checking.cxx

0001 //R__LOAD_LIBRARY(libfmt.so)
0002 //#include "fmt/core.h"
0003 R__LOAD_LIBRARY(libDDG4IO.so)
0004 //
0005 //#include "DD4hep/Detector.h"
0006 #include "DDG4/Geant4Data.h"
0007 //#include "DDRec/CellIDPositionConverter.h"
0008 //#include "DDRec/SurfaceManager.h"
0009 //#include "DDRec/Surface.h"
0010 #include "ROOT/RDataFrame.hxx"
0011 //
0012 //#include "lcio2/MCParticleData.h"
0013 //#include "lcio2/ReconstructedParticleData.h"
0014
0015 //#include "Math/Vector3D.h"
0016 //#include "Math/Vector4D.h"
0017 //#include "Math/VectorUtil.h"
0018 #include "TCanvas.h"
0019 //#include "TLegend.h"
0020 //#include "TMath.h"
0021 //#include "TRandom3.h"
0022 //#include "TFile.h"
0023 //#include "TH1F.h"
0024 //#include "TH1D.h"
0025 //#include "TTree.h"
0026 #include "TChain.h"
0027 //#include "TF1.h"
0028 #include <random>
0029 //#include "lcio2/TrackerRawDataData.h"
0030 //#include "lcio2/TrackerRawData.h"
0031
0032 void simple_checking(const char* fname = "sim_output/output_zdc_photons.edm4hep.root"){
0033 std::cout << "testing 1\n";
0034   ROOT::EnableImplicitMT(); // Tell ROOT you want to go parallel
0035   //using namespace lcio2;
0036   double degree = TMath::Pi()/180.0;
0037
0038   TChain* t = new TChain("events");
0039   t->Add(fname);
0040
0041   ROOT::RDataFrame d0(*t);//, {"GEMTrackerHintits","MCParticles"});
0042 std::cout << "testing 2\n";
0043   //std::cout << t->GetBranch("GEMTrackerHits")->GetClassName() << std::endl;
0044   //std::vector<dd4hep::sim::Geant4Tracker::Hit*>
0045
0046   // -------------------------
0047   // Get the DD4hep instance
0048   // Load the compact XML file
0049   // Initialize the position converter tool
0050   //dd4hep::Detector& detector = dd4hep::Detector::getInstance();
0051   //detector.fromCompact("gem_tracker_disc.xml");
0052   //dd4hep::rec::CellIDPositionConverter cellid_converter(detector);
0053
0054   //// -------------------------
0055   //// Get the surfaces map
0056   //dd4hep::rec::SurfaceManager& surfMan = *detector.extension<dd4hep::rec::SurfaceManager>() ;
0057   //auto surfMap = surfMan.map( "world" ) ;
0058
0059   auto nhits = [] (std::vector<dd4hep::sim::Geant4Calorimeter::Hit*>& hits){ return (int) hits.size(); };
0060 std::cout << "testing 3\n";
0061   //for(const auto& h: hits){
0062   //  //std::cout << (h->position/10.0) << std::endl;
0063   //  //std::cout << cellid_converter.position(h->cellID) << std::endl;
0064   //  //dd4hep::rec::SurfaceMap::const_iterator
0065   //  const auto si = _surfMap.find( cellid_converter.findContext(h->cellID)->identifier ); //identifier=volumeID
0066   //  dd4hep::rec::ISurface* surf = (si != _surfMap.end() ?  si->second  : 0);
0067   //  dd4hep::rec::Vector3D pos =  surf->origin();//fit_global(pivot[0],pivot[1],pivot[2]);
0068   //  //std::cout << pos.x() << ", " << pos.y() << ", " << pos.z()<< std::endl;
0069   //  // transform lcio units to dd4hep units, see documentation for other functions
0070   //  //DDSurfaces::Vector2D fit_local = surf->globalToLocal( dd4hep::mm * fit_global );
0071   //}
0072   //  return hits.size(); };
0073
0074   //auto digitize_gem_hits =
0075   //  [&](const std::vector<dd4hep::sim::Geant4Tracker::Hit*>& hits) {
0076   //    std::vector<lcio2::TrackerRawDataData> digi_hits;
0077   //    std::normal_distribution<> time_dist(0,2.0);
0078   //    std::normal_distribution<> adc_dist(5.0,3.0);
0079
0080   //    std::map<int64_t,lcio2::TrackerRawDataData> hits_by_id;
0081   //    for(const auto& h: hits) {
0082   //      //lcio2::TrackerRawDataData ahit;
0083   //      auto& ahit = hits_by_id[(int64_t)h->cellID];
0084   //      auto pos = h->position/10.0; //cm
0085
0086   //      ahit.cellID0   = h->cellID;
0087   //      ahit.cellID1   = h->cellID;
0088   //      ahit.channelID = h->cellID;
0089   //      //fmt::print("{} vs {} vs {}\n", id1, id2,id3);
0090   //      fmt::print("{} vs {}\n", h->cellID, ahit.cellID0);
0091   //      // time is not kept from dd4hep hit, instead using z position as crude substitute
0092   //      ahit.time = pos.z() + time_dist(gen);
0093   //      ahit.adc  = adc_dist(gen);
0094   //      //digi_hits.push_back(ahit);
0095   //    }
0096   //    for (auto& [cell, hit] : hits_by_id) {
0097   //      //fmt::print("{} vs {}\n", cell, hit.cellID0);
0098   //      digi_hits.push_back(hit);
0099   //    }
0100   //    return digi_hits;
0101   //  };
0102
0103   auto d1 = d0.Define("nhits", nhits, {"ZDCHits"})
0104               //.Filter([](int n){ return (n>4); },{"nhits"})
0105               //.Define("delta",hit_position, {"GEMTrackerHits"})
0106               //.Define("RawTrackerHits", digitize_gem_hits, {"GEMTrackerHits"})
0107               ;
0108 std::cout << "testing 4\n";
0109   auto h0 = d1.Histo1D(TH1D("h0", "nhits; ", 20, 0,20), "nhits");
0110
0111   auto n0 = d1.Filter([](int n){ return (n>0); },{"nhits"}).Count();
0112
0113   TCanvas* c = new TCanvas();
0114
0115   //d1.Snapshot("digitized_EVENT","test_gem_tracker_digi.root");
0116 std::cout << "testing 5\n";
0117   std::cout << *n0 << " events with nonzero hits\n";
0118 std::cout << "testing 6\n";
0119   if(*n0<5) {
0120     std::quick_exit(1);
0121   }
0122
0123 }
0124