TDataFrame dealing with heppy+podio

heppyanalysis.cxx

#include <ROOT/TDataFrame.hxx>
#include "TVector3.h"
#include <TCanvas.h>
#include <TApplication.h>
#include <TSystem.h>
#include <ROOT/TVec.hxx>

/*
 * Running from lxplus:
 *   source /cvmfs/sft.cern.ch/lcg/views/dev3/latest/x86_64-slc6-gcc62-opt/setup.sh
 *   g++ -g -o analysis heppyAnalysis.cxx `root-config --cflags --glibs` -lROOTVecOps -lROOTDataFrame
 */

// Test to reproduce Heppy analysis
int main(){

   TApplication app("app", nullptr, nullptr);

   auto fname = "FCCDelphesOutput.root";
   TFile f1(fname);

   //f1.MakeProject("dictsDelphes", "*", "RECREATE+");
   gSystem->Load("dictsDelphes/dictsDelphes.so");

   using ints   = ROOT::Experimental::VecOps::TVec<int>;
   using floats = ROOT::Experimental::VecOps::TVec<float>;
   using u_ints = ROOT::Experimental::VecOps::TVec<unsigned int>;
   
   ROOT::Experimental::TDataFrame df("events", fname);

   auto getPts = [](floats &pxs, floats &pys){
       auto all_pts = sqrt(pxs * pxs + pys * pys);
       return all_pts;
   };

   // Apply simple selectors
   auto selectors = df.Define("j_pts", getPts, {"jets.core.p4.px", "jets.core.p4.py"})
                      .Define("m_pts", getPts, {"muons.core.p4.px", "muons.core.p4.py"})
                      .Define("e_pts", getPts, {"electrons.core.p4.px", "electrons.core.p4.py"})
                      .Define("j_pts_30", "j_pts[j_pts > 30]") // selector: jets_30
                      .Define("e_pts_30", "e_pts[e_pts > 30]") // selector: sel_electrons
                      .Define("m_pts_30", "m_pts[m_pts > 30]") // selector: sel_muons
                    ;

   auto get_sumpt = [](floats &pts){
       float sum = 0;
       for(auto pt : pts)
           sum += pt;
       return sum;
   };

   // Define total pt for the particles
   // Heppy: ptc.iso.sumpt
   // iso-> Isolation: '''Holds the results of an isolation calculation.'''
   auto ptsums =  selectors.Define("m_sumpt", get_sumpt, {"m_pts"})
                           .Define("e_sumpt", get_sumpt, {"e_pts"})
                   ;
   
   auto getVects = [](floats &pxs, floats &pys, floats &pzs){
      ROOT::Experimental::VecOps::TVec<TVector3> all_vects;
      for(int i=0; i<pxs.size(); ++i){
         all_vects.emplace_back(TVector3(pxs[i], pys[i], pzs[i]));
      }
      return all_vects;
   };

   // New df to quickly identify problems
   auto vectors = df.Define("j_vects", getVects, {"jets.core.p4.px", "jets.core.p4.py", "jets.core.p4.pz"})
                    .Define("e_vects", getVects, {"electrons.core.p4.px", "electrons.core.p4.py", "electrons.core.p4.pz"})
                    .Define("m_vects", getVects, {"muons.core.p4.px", "muons.core.p4.py", "muons.core.p4.pz"})
                 ;

   //using TVector3s ROOT::Experimental::VecOps::TVec<TVector3>; 
   //auto getDeltas = [](TVector3s ptcs_a, TVector3s ptcs_b){
   //   for(auto v : ptcs_a)
   //}

   // Not working (see example.root)
   // Different vector sizes
   //   m_sumpt : 100 elements
   //   m_pts : 86 elements
   auto muons_h = ptsums.Define("good_muons", "m_sumpt[(m_sumpt / m_pts) < 0.1]").Histo1D("good_muons");

   TCanvas c12;
   muons_h->Draw();
   c12.Draw();

   app.Run();
   return 0;
}