pid/compact/scripts/cherenkov_xy_plot.cxx

Go to the documentation of this file.

R__LOAD_LIBRARY(libfmt.so)
#include "fmt/core.h"
R__LOAD_LIBRARY(libDDG4IO.so)
R__LOAD_LIBRARY(libGenDetectors.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 "TCanvas.h"
#include "TChain.h"
#include <random>
 
#include "npdet/PhotoMultiplierHit.h"
 
//#include "lcio2/TrackerRawDataData.h"
//#include "lcio2/TrackerRawData.h"
 
using PMHit     = npdet::PhotoMultiplierHit;
using PMHitRVec = ROOT::VecOps::RVec<npdet::PhotoMultiplierHit*>;
using PMTHitVector = std::vector<npdet::PhotoMultiplierHit*>;
 
void cherenkov_xy_plot(const char* fname = "derp.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("EVENT");
  t->Add(fname);
 
  ROOT::RDataFrame d0(*t);
  std::cout << t->GetBranch("ForwardRICHHits")->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("GenericRICH_example.xml");
  //dd4hep::rec::CellIDPositionConverter cellid_converter(detector);
 
  //dd4hep::rec::SurfaceManager& surfMan = *detector.extension<dd4hep::rec::SurfaceManager>() ;
  //auto surfMap = surfMan.map( "world" ) ;
  
  //auto nhits = [] (std::vector<npdet::sim::Geant4Tracker::Hit*>& hits){ return (int) hits.size(); };
 
  auto pmt_x = [&](const PMTHitVector& hits) {
    //return ROOT::VecOps::Map(hits, [](const PMHit& h) { return h.position.x(); });
    std::vector<double> res;
    for(auto h : hits) {
      if(h->position.z() < 400) {
        res.push_back(h->position.x());
      }
    }
    return res;
  };
  auto pmt_y = [&](const PMTHitVector& hits) {
    std::vector<double> res;
    for(auto h : hits) {
      if(h->position.z() < 400) {
        res.push_back(h->position.y());
      }
    }
    return res;
  };
 
  //auto pmt_x = [&](const PMTHitVector& hits) {
  //  return ROOT::VecOps::Map(hits, [](const PMHit& h) { return h.position.x(); });
  //};
 
  auto d1 = d0
  .Define("nhits", [](const PMTHitVector& hits) { return hits.size(); }, {"ForwardRICHHits"})
  .Define("pmt_x",pmt_x,{"ForwardRICHHits"})
  .Define("pmt_y",pmt_y,{"ForwardRICHHits"});
 
  auto h1    = d1.Histo1D("nhits");
  auto h1_x  = d1.Histo1D("pmt_x");
  auto h2 = d1.Histo2D({"cerxy", "cer_xy", 100u, -1000, -1000, 100u, -1000, -1000}, "pmt_x", "pmt_y");
  auto graph = d1.Graph("nhits", "nhits");
 
  TCanvas* c = new TCanvas();
 
  h2->DrawCopy("colz");
 
  //graph->DrawClone("ap");
  //graph->SetMarkerStyle(20);
  //c->Update();
 
}