rec_central_electrons.cxx

#include "ROOT/RDataFrame.hxx"
#include "TCanvas.h"
#include "TLegend.h"
#include "TH1D.h"
#include "TProfile.h"

#include <iostream>

R__LOAD_LIBRARY(libeicd.so)
R__LOAD_LIBRARY(libDD4pod.so)
#include "dd4pod/Geant4ParticleCollection.h"
#include "eicd/TrackParametersCollection.h"
#include "eicd/ClusterCollection.h"
#include "eicd/ClusterData.h"
#include "eicd/TrackerHitCollection.h"

using ROOT::RDataFrame;
using namespace ROOT::VecOps;

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;
}

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 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;
};

auto delta_p_over_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)/p1);
    }
  }
  return res;
};

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("theta_thrown", theta, {"thrownP"})
                 .Define("theta0", "theta_thrown[0]")
                 .Define("nTracks", "outputTrackParameters.size()")
                 .Define("p_track", p_track, {"outputTrackParameters"})
                 .Define("p_track1", p_track, {"outputTrackParameters1"})
                 .Define("p_track2", p_track, {"outputTrackParameters2"})
                 .Define("delta_p0",delta_p, {"p_track", "p_thrown"})
                 .Define("delta_p1",delta_p, {"p_track1", "p_thrown"})
                 .Define("delta_p2",delta_p, {"p_track2", "p_thrown"})
                 .Define("delta_p_over_p0",delta_p_over_p, {"p_track", "p_thrown"})
                 .Define("delta_p_over_p1",delta_p_over_p, {"p_track1", "p_thrown"})
                 .Define("delta_p_over_p2",delta_p_over_p, {"p_track2", "p_thrown"})
                 .Define("N_VtxBarrelHits",[](std::vector<eic::TrackerHitData> hits) { return hits.size();},{"VertexBarrelRecHits"})
                 .Define("N_SiBarrelHits", [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"TrackerBarrelRecHits"})
                 .Define("N_SiEndcapHits", [](std::vector<eic::TrackerHitData> hits) { return hits.size();}, {"TrackerEndcapRecHits"})
                 ;

  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_p0  = df0.Histo1D({"h_delta_p0", "Truth Track Init; GeV/c ",  100, -10, 10}, "delta_p0");
  auto h_delta_p1 = df0.Histo1D({"h_delta_p1", "Ecal Cluster Init; GeV/c ", 100, -10, 10}, "delta_p1");
  auto h_delta_p2 = df0.Histo1D({"h_delta_p2", "Ecal Cluster, innner vtx hit Init; GeV/c ", 100, -10, 10}, "delta_p2");

  auto h_delta_p0_over_p = df0.Histo1D({"h_delta_p0_over_p",  "Truth Track Init; delta p/p ",  100, -0.1, 0.1}, "delta_p_over_p0");
  auto h_delta_p1_over_p = df0.Histo1D({"h_delta_p1_over_p", "Ecal Cluster Init; delta p/p ", 100, -0.1, 0.1}, "delta_p_over_p1");
  auto h_delta_p2_over_p = df0.Histo1D({"h_delta_p2_over_p", "Ecal Cluster, innner vtx hit Init; delta p/p ", 100, -0.1, 0.1}, "delta_p_over_p2");

  auto hSiBarrel_N_vs_theta = df0.Histo1D({"hSiBarrel_N_vs_theta", "; #theta [deg.]",   20, 0, 180 }, "theta0", "N_SiBarrelHits");
  auto hSiEndcap_N_vs_theta = df0.Histo1D({"hSiEndcap_N_vs_theta", "; #theta [deg.]",   20, 0, 180 }, "theta0", "N_SiEndcapHits");
  auto hVtxBarrel_N_vs_theta = df0.Histo1D({"hVtxBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_VtxBarrelHits");

  auto hSiBarrel_Nhits  = df0.Histo1D({"hSiBarrel_Nhits", "; #theta [deg.]",   20, 0, 20 }, "theta0", "N_SiBarrelHits");
  auto hSiEndcap_Nhits  = df0.Histo1D({"hSiEndcap_Nhits", "; #theta [deg.]",   20, 0, 20 }, "theta0", "N_SiEndcapHits");
  auto hVtxBarrel_Nhits = df0.Histo1D({"hVtxBarrel_Nhits", "; #theta [deg.]", 20, 0, 20 }, "theta0", "N_VtxBarrelHits");

  auto hSiBarrel_Ntheta = df0.Histo1D({"hSiBarrel_Ntheta", "; #theta [deg.]",   20, 0, 180 }, "theta0");
  auto hSiEndcap_Ntheta = df0.Histo1D({"hSiEndcap_Ntheta", "; #theta [deg.]",   20, 0, 180 }, "theta0");
  auto hVtxBarrel_Ntheta = df0.Histo1D({"hVtxBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");

  auto c = new TCanvas();

  h_nTracks->DrawCopy();
  c->SaveAs("results/tracking/rec_central_electrons_nTracks.png");
  c->SaveAs("results/tracking/rec_central_electrons_nTracks.pdf");

  h_pTracks->DrawCopy();
  c->SaveAs("results/tracking/rec_central_electrons_pTracks.png");
  c->SaveAs("results/tracking/rec_central_electrons_pTracks.pdf");

  c = new TCanvas();
  THStack * hs = new THStack("hs_delta_p","; GeV/c ");
  TH1D* h1 = (TH1D*) h_delta_p0->Clone();
  hs->Add(h1);
  h1 = (TH1D*) h_delta_p1->Clone();
  h1->SetLineColor(2);
  hs->Add(h1);
  h1 = (TH1D*) h_delta_p2->Clone();
  h1->SetLineColor(4);
  h1->SetFillStyle(3001);
  h1->SetFillColor(4);
  hs->Add(h1);
  hs->Draw("nostack");
  c->BuildLegend();
  c->SaveAs("results/tracking/rec_central_electrons_delta_p.png");
  c->SaveAs("results/tracking/rec_central_electrons_delta_p.pdf");

  c  = new TCanvas();
  hs = new THStack("hs_delta_p_over_p","; delta p/p ");
  h1 = (TH1D*) h_delta_p0_over_p->Clone();
  hs->Add(h1);
  h1 = (TH1D*) h_delta_p1_over_p->Clone();
  h1->SetLineColor(2);
  hs->Add(h1);
  h1 = (TH1D*) h_delta_p2_over_p->Clone();
  h1->SetLineColor(4);
  h1->SetFillStyle(3001);
  h1->SetFillColor(4);
  hs->Add(h1);
  hs->Draw("nostack");
  c->BuildLegend();
  c->SaveAs("results/tracking/rec_central_electrons_delta_p_over_p.png");
  c->SaveAs("results/tracking/rec_central_electrons_delta_p_over_p.pdf");

  c  = new TCanvas();
  hs = new THStack("n_hits","; #theta  ");
  h1 = (TH1D*) hSiBarrel_N_vs_theta->Clone();
  auto h2 = (TH1D*) hSiBarrel_Ntheta->Clone();
  h1->Divide(h2);
  hs->Add(h1);
  h1 = (TH1D*) hSiEndcap_N_vs_theta->Clone();
  h2 = (TH1D*) hSiEndcap_Ntheta->Clone();
  h1->Divide(h2);
  h1->SetLineColor(2);
  hs->Add(h1);
  //h1 = (TH1D*) hVtxBarrel_vs_theta->Clone();
  //h1->SetLineColor(4);
  //h1->SetFillStyle(3001);
  //h1->SetFillColor(4);
  //hs->Add(h1);
  hs->Draw("nostack, hist");
  c->BuildLegend();
  c->SaveAs("results/tracking/rec_central_electrons_n_hits_vs_theta.png");
  c->SaveAs("results/tracking/rec_central_electrons_n_hits_vs_theta.pdf");

  c  = new TCanvas();
  hs = new THStack("theta","; #theta  ");
  h1 = (TH1D*) hSiBarrel_N_vs_theta->Clone();
  h2 = (TH1D*) hSiBarrel_Ntheta->Clone();
  //h1->Divide(h2);
  hs->Add(h2);
  h1 = (TH1D*) hSiEndcap_N_vs_theta->Clone();
  h2 = (TH1D*) hSiEndcap_Ntheta->Clone();
  //h1->Divide(h2);
  h1->SetLineColor(2);
  h2->SetLineColor(2);
  hs->Add(h2);
  //h1 = (TH1D*) hVtxBarrel_vs_theta->Clone();
  //h1->SetLineColor(4);
  //h1->SetFillStyle(3001);
  //h1->SetFillColor(4);
  //hs->Add(h1);
  hs->Draw("nostack hist");
  c->BuildLegend();
  c->SaveAs("results/tracking/rec_central_electrons_theta.png");
  c->SaveAs("results/tracking/rec_central_electrons_theta.pdf");

  c  = new TCanvas();
  hs = new THStack("hits","; hits  ");
  h1 = (TH1D*) hSiBarrel_Nhits->Clone();
  hs->Add(h1);
  h1 = (TH1D*) hSiEndcap_Nhits->Clone();
  h1->SetLineColor(2);
  h2->SetLineColor(2);
  hs->Add(h2);
  h1 = (TH1D*) hVtxBarrel_Nhits->Clone();
  h1->SetLineColor(4);
  h1->SetFillStyle(3001);
  h1->SetFillColor(4);
  hs->Add(h1);
  hs->Draw("nostack hist");
  c->BuildLegend();
  c->SaveAs("results/tracking/rec_central_electrons_nhits.png");
  c->SaveAs("results/tracking/rec_central_electrons_nhits.pdf");


  return 0;
}