barrel_clusters.cxx

#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_over_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)/p.E();
      result.push_back(dE );
      if( dE < best) {
        best = dE;
      }
    }
  }
  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("delta_E_over_E", delta_E_over_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, -3, 3}, "delta_E");
  auto h_delta_E_over_E =
      d0.Histo1D({"h_delta_E_over_E", ";  #frac{E_{thrown}-E_{cluster}}{E_{thrown}} ", 100, -1, 1}, "delta_E_over_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");

  h_delta_E_over_E->DrawCopy();
  c->SaveAs("results/clustering/barrel_clusters_delta_E_over_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;
}