diff --git a/scripts/tutorial2_cell_size.cxx b/scripts/tutorial2_cell_size.cxx new file mode 100644 index 0000000000000000000000000000000000000000..90d248d2b9a80c3f4071fb31d7d9dff8d63cfb43 --- /dev/null +++ b/scripts/tutorial2_cell_size.cxx @@ -0,0 +1,118 @@ +R__LOAD_LIBRARY(libfmt.so) +#include "fmt/core.h" + +#include "Math/Vector3D.h" +#include "Math/Vector4D.h" +#include "Math/VectorUtil.h" + +#include "ROOT/RDataFrame.hxx" +#include "TCanvas.h" +#include "TChain.h" +#include "TF1.h" +#include "TFile.h" +#include "TH1D.h" +#include "TH1F.h" +#include "TLegend.h" +#include "TMath.h" +#include "TRandom3.h" +#include "TTree.h" + +#include <algorithm> +#include <iterator> +#include <tuple> +#include <vector> +#include <cmath> + +#include "DD4hep/Detector.h" +#include "DDRec/CellIDPositionConverter.h" + +#include "dd4pod/Geant4ParticleCollection.h" +#include "dd4pod/TrackerHitCollection.h" + +/** Hit position example. + * + */ +void tutorial2_cell_size(const char* fname = "gem_tracker_sim.root") { + using namespace ROOT::Math; + + TChain* t = new TChain("events"); + t->Add(fname); + + ROOT::RDataFrame d0(*t, {"GEMTrackerHits", "MCParticles"}); + //std::cout << d0.DescribeDataset() << std::endl; + + // ------------------------- + // Get the DD4hep instance + // Load the compact XML file + // Initialize the position converter tool + dd4hep::Detector& detector = dd4hep::Detector::getInstance(); + detector.fromCompact("gem_tracker.xml"); + dd4hep::rec::CellIDPositionConverter cellid_converter(detector); + + // Simple lambda to define nhits branch + auto nhits = [](const std::vector<dd4pod::TrackerHitData>& evt) { return (int)evt.size(); }; + + auto local_position = [&](const std::vector<dd4pod::TrackerHitData>& hits) { + std::vector<std::array<double, 2>> result; + for (const auto& h : hits) { + // The actual hit position: + auto pos0 = (h.position); + + // **This is the most important part.** It provides a way of getting the + // geometry information in a flexible way. + // The segmentation hit postion: + auto pos1 = cellid_converter.position(h.cellID); + fmt::print(" Hit Position : {},{},{}\n", pos0.x / 10.0, pos0.y / 10.0, pos0.z / 10.0); + fmt::print("Segmentation-Cell Position : {},{},{}\n", pos1.x(), pos1.y(), pos1.z()); + + auto context = cellid_converter.findContext( h.cellID ) ; + dd4hep::Readout r = cellid_converter.findReadout( context->element ) ; + dd4hep::Segmentation seg = r.segmentation() ; + auto cell_dim = seg.cellDimensions(h.cellID); + std::cout << " dim "; + for(const auto& dim : cell_dim) { + std::cout << dim << ", "; + } + std::cout << "\n"; + + result.push_back({pos0.x / 10.0 - pos1.x(), pos0.y / 10.0 - pos1.y()}); + } + return result; + }; + + auto x_pos = [&](const std::vector<std::array<double, 2>>& xypos) { + std::vector<double> result; + for (const auto& h : xypos) { + result.push_back(h.at(0)); + } + return result; + }; + auto y_pos = [&](const std::vector<std::array<double, 2>>& xypos) { + std::vector<double> result; + for (auto h : xypos) { + result.push_back(h.at(1)); + } + return result; + }; + + auto d1 = d0.Define("nhits", nhits, {"GEMTrackerHits"}) + .Filter([=](const std::vector<dd4pod::TrackerHitData>& hits) { + for (auto h : hits) { + auto pos = ROOT::Math::XYZVector(h.position.x,h.position.y,h.position.z); + if ((pos.r() > 100.0) && (std::abs(pos.phi()-M_PI/2.0)< M_PI/6)) { + return true; + } + } + return false; + }, + {"GEMTrackerHits"}) + .Define("xy_hit_pos", local_position, {"GEMTrackerHits"}) + .Define("x_pos", x_pos, {"xy_hit_pos"}) + .Define("y_pos", y_pos, {"xy_hit_pos"}); + + auto h_local_pos = d1.Histo2D(TH2D("h_local_pos", ";x [cm]; y [cm] ", 100, -5.0, 5.0, 100, -5.0, 5.0), "x_pos", "y_pos"); + + TCanvas* c = new TCanvas(); + h_local_pos->DrawClone("colz"); + c->SaveAs("results/hit_position.png"); +}