Skip to content
Snippets Groups Projects
Commit c5fdeb17 authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

new file: scripts/tutorial2_cell_size.cxx

parent 754cc8a4
No related branches found
No related tags found
No related merge requests found
Pipeline #10099 failed
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");
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment