load_geo.cxx

Go to the documentation of this file.

R__LOAD_LIBRARY(libDDG4IO.so)
 
#include "DD4hep/Detector.h"
#include "DDG4/Geant4Data.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
#include "ROOT/RDataFrame.hxx"
//
//#include <vector>
//#include <tuple>
//#include <algorithm>
//#include <iterator>
//
#include "DD4hep/Detector.h"
#include "DDG4/Geant4Data.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
#include "ROOT/RDataFrame.hxx"
 
//#include "lcio2/MCParticleData.h"
//#include "lcio2/ReconstructedParticleData.h"
 
//#include "Math/Vector3D.h"
//#include "Math/Vector4D.h"
//#include "Math/VectorUtil.h"
#include "TCanvas.h"
//#include "TLegend.h"
//#include "TMath.h"
//#include "TRandom3.h"
//#include "TFile.h"
//#include "TH1F.h"
//#include "TH1D.h"
//#include "TTree.h"
#include "TChain.h"
//#include "TF1.h"
void load_geo(const char* fname = "test_tracker_disc.root"){
 
  using namespace ROOT::Math;
 
  ROOT::EnableImplicitMT(4);
  TChain* t = new TChain("EVENT");
  t->Add(fname);
 
  ROOT::RDataFrame d0(*t, {"GEMTrackerHits","MCParticles"});
 
  //How to get the type of the initial branch: (vector<dd4hep::sim::Geant4Tracker::Hit*>)
  //std::cout << t->GetBranch("GEMTrackerHits")->GetClassName() << std::endl;
  
  // -------------------------
  // Get the DD4hep instance
  // Load the compact XML file
  // Initialize the position converter tool
  dd4hep::Detector& detector = dd4hep::Detector::getInstance();
  detector.fromCompact("IRChamber_example.xml");
  dd4hep::rec::CellIDPositionConverter cellid_converter(detector);
 
  // -------------------------
  // Get the surfaces map
  //dd4hep::rec::SurfaceManager& surfMan = *detector.extension<dd4hep::rec::SurfaceManager>() ;
  //auto surfMap = surfMan.map( "world" ) ;
 
  //auto nhits = [] (const std::vector<dd4hep::sim::Geant4Tracker::Hit*>& hits){ return hits.size(); };
 
 
  //auto deltax = [&](const std::vector<dd4hep::sim::Geant4Tracker::Hit*>& hits){
  //  std::vector<double> result;
  //  for(const auto& h: hits) {
  //    auto cell = h->cellID;
  //    // The actual hit position:
  //    XYZVector pos0 = (h->position/10.0);
  //    // The segmentation hit postion:
  //    XYZVector pos1 = cellid_converter.position(cell);
 
  //    auto seg_pos = cellid_converter.position(cell);
 
  //    //auto                 context  = cellid_converter.findContext( cell ) ;
  //    //dd4hep::Readout      r        = cellid_converter.findReadout( context->element ) ;
  //    //dd4hep::Segmentation seg      = r.segmentation() ;
  //    //auto                 cell_dim = seg.cellDimensions(cell);
  //    auto                 cell_dim = cellid_converter.cellDimensions(cell);
  //    std::cout << " dim ";
  //    for(const auto& dim : cell_dim) {
  //      std::cout << dim << ", ";
  //    }
  //    //XYZVector pos1 = pos11;
  //    // Use the cellID to get the VolumeID (VolumeManagerContext.identifier = VolumeID)
  //    const auto si = surfMap->find( cellid_converter.findContext(h->cellID)->identifier );
  //    // Get the Surface (http://test-dd4hep.web.cern.ch/test-dd4hep/doxygen/html/classdd4hep_1_1rec_1_1_i_surface.html)
  //    dd4hep::rec::ISurface* surf = (si != surfMap->end() ?  si->second  : nullptr);
  //    dd4hep::rec::Vector3D  pos2;
  //    // Get the origin of the surface volume
  //    if(!surf) pos2 = surf->origin();
  //    std::cout << "              Hit Position : " << pos0 << std::endl;
  //    std::cout << "Segmentation-Cell Position : " << pos1 << std::endl;
  //    std::cout << "   Surface Origin Position : " << pos2 << std::endl;
  //    result.push_back( VectorUtil::RotateZ((pos1-pos0), -1.0*pos1.Phi()+TMath::Pi()/2.0).x() );
  //  }
  //  return result;
  //};
 
  //auto deltay = [&](const std::vector<dd4hep::sim::Geant4Tracker::Hit*>& hits){
  //  std::vector<double> result;
  //  for(const auto& h: hits) {
  //    // The actual hit position:
  //    XYZVector pos0 = (h->position/10.0);
  //    // The segmentation hit postion:
  //    XYZVector pos1 = cellid_converter.position(h->cellID);
  //    // Use the cellID to get the VolumeID (VolumeManagerContext.identifier = VolumeID)
  //    const auto si = surfMap->find( cellid_converter.findContext(h->cellID)->identifier );
  //    // Get the Surface (http://test-dd4hep.web.cern.ch/test-dd4hep/doxygen/html/classdd4hep_1_1rec_1_1_i_surface.html)
  //    dd4hep::rec::ISurface* surf = (si != surfMap->end() ?  si->second  : nullptr);
  //    dd4hep::rec::Vector3D  pos2;
  //    // Get the origin of the surface volume
  //    if(!surf) pos2 = surf->origin();
  //    //std::cout << "              Hit Position : " << pos0 << std::endl;
  //    //std::cout << "Segmentation-Cell Position : " << pos1 << std::endl;
  //    //std::cout << "   Surface Origin Position : " << pos2 << std::endl;
  //    result.push_back( VectorUtil::RotateZ((pos1-pos0), -1.0*pos1.Phi()+TMath::Pi()/2.0).y() );
  //  }
  //  return result;
  //};
 
  //auto d1 = d0
  //  .Define("nhits",nhits, {"GEMTrackerHits"})
  //  .Filter([=](const std::vector<dd4hep::sim::Geant4Tracker::Hit*>& hits) {
  //    for(const auto& h: hits) {
  //      auto pos = h->position/10;
  //      if( TMath::Sqrt(pos.x()*pos.x()+ pos.y()*pos.y()) > 30.0 ){
  //        return true;
  //      }
  //    }
  //    return false;}, {"GEMTrackerHits"})
  //  .Define("deltax", deltax, {"GEMTrackerHits"})
  //  .Define("deltay", deltay, {"GEMTrackerHits"});
 
 
  //auto h0 = d1.Histo2D( TH2D("h0", "nhits; ", 100, -2.0, 2.0, 100, -2.0, 2.0), "deltax", "deltay");
 
  //TCanvas* c = new TCanvas();
  //h0->DrawClone("colz");
 
  //d1.Snapshot("EVENT","derp.root");
}