Newer
Older
1
2
3
4
5
6
7
8
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
#include <iostream>
#include "ROOT/RDataFrame.hxx"
#include "dd4pod/Geant4ParticleCollection.h"
#include "eicd/ClusterCollection.h"
#include "eicd/ClusterData.h"
// If root cannot find a dd4pod header make sure to add the include dir to ROOT_INCLUDE_PATH:
// export ROOT_INCLUDE_PATH=$HOME/include:$ROOT_INCLUDE_PATH
using ROOT::RDataFrame;
using namespace ROOT::VecOps;
auto momentum = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
std::vector<double> result;
for (size_t i = 0; i < in.size(); ++i) {
result.push_back(in[i].E());
}
return result;
};
auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
std::vector<double> result;
for (size_t i = 0; i < in.size(); ++i) {
result.push_back(in[i].Theta()*180/M_PI);
}
return result;
};
auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
std::vector<ROOT::Math::PxPyPzMVector> result;
ROOT::Math::PxPyPzMVector lv;
for (size_t i = 0; i < in.size(); ++i) {
lv.SetCoordinates(in[i].psx, in[i].psy, in[i].psz, in[i].mass);
result.push_back(lv);
}
return result;
};
auto pt = [](std::vector<dd4pod::Geant4ParticleData> const& in) {
std::vector<float> result;
for (size_t i = 0; i < in.size(); ++i) {
result.push_back(std::sqrt(in[i].psx * in[i].psx + in[i].psy * in[i].psy));
}
return result;
};
auto eta = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
std::vector<float> result;
ROOT::Math::PxPyPzMVector lv;
for (size_t i = 0; i < in.size(); ++i) {
lv.SetCoordinates(in[i].psx, in[i].psy, in[i].psz, in[i].mass);
result.push_back(lv.Eta());
}
return result;
};
auto delta_E = [](std::vector<ROOT::Math::PxPyPzMVector> const& thrown, const std::vector<eic::ClusterData>& clusters) {
std::vector<double> result;
double best = 1000000.0;
for (const auto& p : thrown) {
for (const auto& c : clusters) {
double dE = (p.E() - c.energy);
result.push_back(dE );
if( dE < best) {
best = dE;
}
}
}
return best;
};
void barrel_clusters(const char* in_fname = "topside/rec_barrel_clusters.root")
{
ROOT::EnableImplicitMT();
ROOT::RDataFrame df("events", in_fname);
auto d0 = df.Define("isThrown", "mcparticles2.genStatus == 1")
.Define("thrownParticles", "mcparticles2[isThrown]")
.Define("thrownP", fourvec, {"thrownParticles"})
.Define("thrownEta", eta, {"thrownParticles"})
.Define("thrownTheta", theta, {"thrownP"})
.Define("thrownMomentum", momentum, {"thrownP"})
.Define("delta_E", delta_E, {"thrownP","SimpleClusters"})
.Define("nclusters", "SimpleClusters.size()")
.Define("Ecluster",
[](const std::vector<eic::ClusterData>& in) {
std::vector<double> res;
for (const auto& i : in)
res.push_back(i.energy);
return res;
},
{"SimpleClusters"});
auto d1 = d0.Filter("nclusters==1");
auto c_nclusters1 = d1.Count();
auto c_thrown = d0.Count();
auto h_eta_thrown = d0.Histo1D({"h_eta_thrown", " ; #eta ", 100, -5.0, 0.0}, "thrownEta");
auto h_theta_thrown = d0.Histo1D({"h_theta_thrown", "; #theta", 100, 30.0, 180.0}, "thrownTheta");
auto h_momentum_thrown = d0.Histo1D({"h_momentum_thrown", "; E [GeV]", 100, 0.0, 30.0}, "thrownMomentum");
auto h_nclusters = d0.Histo1D({"h_nclusters", "; N clusters", 6, 0, 6}, "nclusters");
auto h_Ecluster = d0.Histo1D({"h_Ecluster", "; cluster E [GeV]", 100, 0, 1}, "Ecluster");
auto h_Ecluster1 = d1.Histo1D({"h_Ecluster1", "One cluster events; cluster E [GeV]", 100, 0, 30}, "Ecluster");
auto h_momentum_thrown1 = d1.Histo1D({"h_momentum_thrown", "; E [GeV]", 100, 0.0, 30.0},"thrownMomentum");
auto h_delta_E = d0.Histo1D({"h_delta_E", "; #Delta E [GeV]", 100, -10, 10}, "delta_E");
//auto h_Ecluster2 = d1.Filter("thrownMomentum > 4").Histo1D({"h_Ecluster1", "One cluster events; cluster E [GeV]",100, 0,30},"Ecluster");
//auto h_momentum_thrown2 = d1.Filter("thrownMomentum > 4").Histo1D({"h_momentum_thrown", "; E [GeV]", 100, 0.0, 30.0},"thrownMomentum");
auto c = new TCanvas();
h_eta_thrown->DrawCopy();
c->SaveAs("results/clustering/barrel_clusters_etaThrown.png");
h_theta_thrown->DrawCopy();
c->SaveAs("results/clustering/barrel_clusters_thetaThrown.png");
h_nclusters->DrawCopy();
c->SaveAs("results/clustering/barrel_clusters_nclusters.png");
h_Ecluster->DrawCopy();
h_Ecluster1->SetLineColor(2);
h_Ecluster1->DrawCopy("same");
//h_Ecluster2->SetLineColor(4);
//h_Ecluster2->DrawCopy("same");
c->SaveAs("results/clustering/barrel_clusters_Ecluster.png");
h_momentum_thrown->DrawCopy();
c->SaveAs("results/clustering/barrel_clusters_momentum_thrown.png");
h_delta_E->DrawCopy();
c->SaveAs("results/clustering/barrel_clusters_delta_E.png");
std::cout << (*c_nclusters1) << " / " << (*c_thrown) << " = single cluster events/all \n";
//c->SaveAs("results/crystal_cal_electrons_Ecluster.png");
//std::string outfilename = "rdf_test.root";
//df2.Snapshot("events", outfilename, {"MCParticles_pt", "mcparticles"});
return 0;
}