Newer
Older
R__LOAD_LIBRARY(libfmt.so)
#include "fmt/core.h"
R__LOAD_LIBRARY(libDDG4IO.so)
#include "Math/Vector3D.h"
#include "Math/Vector4D.h"
#include "Math/VectorUtil.h"
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
#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"
#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 "lcio2/MCParticleData.h"
//#include "lcio2/ReconstructedParticleData.h"
/** Hit position example.
*
* Makes use of ROOT's TDataframe
* (https://root.cern.ch/doc/master/classROOT_1_1Experimental_1_1TDataFrame.html)
*
* There are two different ways to produce hits associated with 3d positions in DD4hep:
*
* 1. With the "segmentation" which is a logical construction that divides the volume into pieces.
* A segmentation hit associates a (unique) cellID to identify the position of the senstive element.
*
* 2. With the "surfaces" which produce hits associated with the actual volume's volumeID.
* A surface hit recordes the volumeID and other information for each hit.
*
* It is important to remember that these two identifiers (CellID and VolumeID) are different
* and mixing between them will cause problems. However, the VolumeID can be obtained from the CellID
* but not the other way around around. In this way they can work together to provide geometry
* related hit information.
*
* This example shows how to get both kinds of hit positions.
*
*/
void example_hit_position(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("gem_tracker_disc.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" ) ;
// Simple lambda to define nhits branch
auto nhits = [] (const std::vector<dd4hep::sim::Geant4Tracker::Hit*>& hits){ return hits.size(); };
auto hit_position = [&](const std::vector<dd4hep::sim::Geant4Tracker::Hit*>& hits){
std::vector<std::array<double,2>> result;
for(const auto& h: hits) {
// The actual hit position:
auto pos0 = (h->position/10.0);
// **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);
// 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({(pos1-pos0).x(),(pos1-pos0).y()});
}
return result;
};
auto deltax = [&](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 deltay = [&](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<dd4hep::sim::Geant4Tracker::Hit*>& hits) {
for(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("xy_hit_pos", hit_position, {"GEMTrackerHits"})
.Define("deltax", deltax, {"xy_hit_pos"})
.Define("deltay", deltay, {"xy_hit_pos"});
// .Define("derp",
// [](const std::vector<dd4hep::sim::Geant4Tracker::Hit*>& hits){
// std::vector<double> res;
// std::transform(hits.begin(),
// hits.end(),
// std::back_inserter(res),
// [&](auto h) { return h->position.x()/100.0 - 1.0; } );
// return res;
// }, {"GEMTrackerHits"})
//auto h0 = d1.Histo1D(TH1D("h0", "nhits; ", 100, -10.0,10.0), "deltay");
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");
}