From b589590f55ba4bd7149137953f395d5596f9ac05 Mon Sep 17 00:00:00 2001
From: Whitney Armstrong <warmstrong@anl.gov>
Date: Mon, 2 Nov 2020 20:56:22 -0600
Subject: [PATCH] Added new tracking analysis script. new file:
scripts/rec_central_electrons.cxx
---
tracking/scripts/rec_central_electrons.cxx | 125 +++++++++++++++++++++
1 file changed, 125 insertions(+)
create mode 100644 tracking/scripts/rec_central_electrons.cxx
diff --git a/tracking/scripts/rec_central_electrons.cxx b/tracking/scripts/rec_central_electrons.cxx
new file mode 100644
index 00000000..accea3f0
--- /dev/null
+++ b/tracking/scripts/rec_central_electrons.cxx
@@ -0,0 +1,125 @@
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "eicd/TrackParametersCollection.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterData.h"
+
+//namespace edm4hep {
+//
+//std::vector<float> pt (std::vector<MCParticleData> const& in){
+// std::vector<float> result;
+// for (size_t i = 0; i < in.size(); ++i) {
+// result.push_back(std::sqrt(in[i].momentum.x * in[i].momentum.x + in[i].momentum.y * in[i].momentum.y));
+// }
+// return result;
+//}
+//
+//std::vector<float> eta(std::vector<MCParticleData> const& in){
+// std::vector<float> result;
+// ROOT::Math::PxPyPzMVector lv;
+// for (size_t i = 0; i < in.size(); ++i) {
+// lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
+// result.push_back(lv.Eta());
+// }
+// return result;
+//}
+//
+//std::vector<float> cos_theta(std::vector<MCParticleData> const& in){
+// std::vector<float> result;
+// ROOT::Math::PxPyPzMVector lv;
+// for (size_t i = 0; i < in.size(); ++i) {
+// lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
+// result.push_back(cos(lv.Theta()));
+// }
+// return result;
+//}
+//
+//}
+auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
+ std::vector<double> result;
+ for (size_t i = 0; i < in.size(); ++i) {
+ result.push_back(std::abs(1.0/(in[i].qOverP)));
+ }
+ return result;
+};
+
+
+std::vector<float> 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;
+}
+
+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;
+};
+
+
+int rec_central_electrons(const char* fname = "topside/rec_central_electrons.root") {
+
+ ROOT::EnableImplicitMT();
+ ROOT::RDataFrame df("events", fname);
+
+ auto df0 = df.Define("isThrown", "mcparticles2.genStatus == 1")
+ .Define("thrownParticles", "mcparticles2[isThrown]")
+ .Define("thrownP", fourvec, {"thrownParticles"})
+ .Define("p_thrown", momentum, {"thrownP"})
+ .Define("nTracks", "outputTrackParameters.size()")
+ .Define("p_track", p_track, {"outputTrackParameters"})
+ .Define("delta_p",
+ [](const std::vector<double> &tracks,
+ const std::vector<double> &thrown) {
+ std::vector<double> res;
+ for (const auto &p1 : thrown) {
+ for (const auto &p2 : tracks) {
+ res.push_back(p1 - p2);
+ }
+ }
+ return res;
+ },
+ {"p_track", "p_thrown"});
+
+ auto h_nTracks = df0.Histo1D({"h_nTracks", "; N tracks ", 10, 0, 10}, "nTracks");
+ auto h_pTracks = df0.Histo1D({"h_pTracks", "; GeV/c ", 100, 0, 10}, "p_track");
+ auto h_delta_p = df0.Histo1D({"h_delta_p", "; GeV/c ", 100, -10, 10}, "delta_p");
+
+ auto c = new TCanvas();
+
+ h_nTracks->DrawCopy();
+ c->SaveAs("results/rec_central_electrons_nTracks.png");
+
+ h_pTracks->DrawCopy();
+ c->SaveAs("results/rec_central_electrons_pTracks.png");
+
+ h_delta_p->DrawCopy();
+ c->SaveAs("results/rec_central_electrons_delta_p.png");
+
+ return 0;
+}
--
GitLab