Skip to content
Snippets Groups Projects
example_hit_recon.cxx 2.56 KiB
Newer Older
  • Learn to ignore specific revisions
  • R__LOAD_LIBRARY(libDDG4IO.so)
    
    
    #include <random>
    
    #include "Math/GenVector/VectorUtil.h"
    
    #include "Math/Vector3D.h"
    #include "Math/Vector4D.h"
    
    #include "ROOT/RDataFrame.hxx"
    
    #include "TCanvas.h"
    
    #include "TChain.h"
    #include "TF1.h"
    
    #include "TFile.h"
    #include "TH1D.h"
    
    #include "TH1F.h"
    #include "TMath.h"
    #include "TRandom3.h"
    
    #include "TTree.h"
    
    // DD4hep
    #include "DD4hep/Detector.h"
    #include "DDG4/Geant4Data.h"
    #include "DDRec/CellIDPositionConverter.h"
    #include "DDRec/Surface.h"
    
    #include "DDRec/SurfaceManager.h"
    
    
    #include "lcio2/TrackerHit.h"
    #include "lcio2/TrackerHitData.h"
    
    #include "lcio2/TrackerRawData.h"
    #include "lcio2/TrackerRawDataData.h"
    
    void example_hit_recon(const char* fname = "gem_tracker_digi.root"){
    
    
      ROOT::EnableImplicitMT(); // Tell ROOT you want to go parallel
      using namespace lcio2;
      double degree = TMath::Pi()/180.0;
    
      std::random_device rd;
      std::mt19937 gen(rd());
    
      TChain* t = new TChain("digitized_EVENT");
      t->Add(fname);
    
    
      ROOT::RDataFrame d0(*t);
    
      //std::cout << t->GetBranch("GEMTrackerHits")->GetClassName() << std::endl;
    
      
      auto nhits = [] (const std::vector<dd4hep::sim::Geant4Tracker::Hit*>& hits){ return hits.size(); };
    
      dd4hep::Detector& detector = dd4hep::Detector::getInstance();
    
      detector.fromCompact("gem_tracker.xml");
    
      dd4hep::rec::CellIDPositionConverter cellid_converter(detector);
    
      dd4hep::rec::SurfaceManager& surfMan = *detector.extension<dd4hep::rec::SurfaceManager>() ;
      //const SurfaceMap& _surfMap = *surfMan.map( "world" ) ;
      const auto& _surfMap = *surfMan.map( "world" ) ;
    
      auto recon_gem_hits = 
        [&](const std::vector<lcio2::TrackerRawDataData>&        digi_hits)
        {
          std::vector<lcio2::TrackerHitData> recon_hits;
          int i_hit = 0;
          for(const auto& h: digi_hits) {
            //std::cout << hits.at(i_hit)->cellID  << " vs " << h.cellID0 << "\n";
            //auto pos = cellid_converter.position(hits.at(i_hit)->cellID );
            auto pos = cellid_converter.position(h.cellID0);//hits.at(i_hit)->cellID );
            std::cout << pos << std::endl;
            lcio2::TrackerHitData ahit;
            ahit.cellID0   = h.cellID0;
            ahit.cellID1   = h.cellID1;
            ahit.EDep      = h.adc;
            ahit.EDepError = 0.01;
            ahit.time      = h.time;
            ahit.position  = {pos.x(), pos.y(), pos.z()};
            recon_hits.push_back(ahit);
            i_hit++;
          }
          return recon_hits;
        };
    
      auto d1 = d0
      .Filter("nhits>3")
      .Define("TrackerHits", recon_gem_hits, {"RawTrackerHits"});
    
      //TCanvas* c = new TCanvas();
      //h0->DrawClone();
    
    
      d1.Snapshot("recon_EVENT","gem_tracker_recon.root");