Hello,
I would like to use a method of a Lorentz 4-vector on a RVec ROOT::Math::PtEtaPhiEVector
I can easily do this on the leading RVec elements, but how can I do this on all elements
I can not figure out how to do this.
I wrote a small example where I calculate the invariant mass of all photon pairs where I explain the
problem:
#!/usr/bin/env python
#
import ROOT
import os,sys
#include c++ in python
ROOT.gInterpreter.Declare(
"""
using Vec_t = ROOT::VecOps::RVec<float>;
using Vec4_t = ROOT::Math::PtEtaPhiEVector;
ROOT::VecOps::RVec< Vec4_t > ComputeInvariantMass(Vec_t& pt, Vec_t& eta, Vec_t& phi, Vec_t& e) {
ROOT::VecOps::RVec< Vec4_t> pair;
pair.reserve(pt.size());
for (int i=0; i<pt.size(); i++) {
for (int j=i+1; j<pt.size(); j++) {
ROOT::Math::PtEtaPhiEVector p1(pt[i], eta[i], phi[i], e[i]);
ROOT::Math::PtEtaPhiEVector p2(pt[j], eta[j], phi[j], e[j]);
pair.emplace_back(p1+p2);
}
}
return pair;
}
""")
#
if __name__ == "__main__":
path = "root://eospublic.cern.ch//eos/opendata/atlas/OutreachDatasets/2020-01-22/"
filename=path+'GamGam/Data/data_A.GamGam.root'
print ('filename= ',filename)
d = ROOT.RDataFrame("mini", filename)
d.Display()
# taking the leading element works
# d = d.Define("photonmass", 'ComputeInvariantMass(photon_pt, photon_eta, photon_phi, photon_E)[0].mass()')
#how to calculate the invariant mass of each element in the vector
d = d.Define("photonmass", 'ComputeInvariantMass(photon_pt, photon_eta, photon_phi, photon_E).mass()')
# this can not work like that since I need to access all elements of the RVEC and then execute the method
h1 = d.Histo1D("photonmass")
h1.Draw()
I am using root version 6.26/00 on lxplus7.
Thank you for your help.
Tancredi
Hi @tancredi ,
ComputeInvariantMass
returns an RVec
of 2 PtEtaPhiEVector
s. You have a few options to transform that into an RVec of 2 invariant masses (i.e. 2 floats).
-
you can do it directly in ComputeInvariantMass
, e.g. by calling pair.emplace_back((p1+p2).mass())
instead of pair.emplace_back(p1 + p2)
-
you can write a separate function that does the transformation, e.g. using the Map helper function from RVec, but you can also just do an explicit loop:
RVecD InvariantMasses(RVec<Vec4_t> &vecs) {
return Map(vecs, [](Vec4_t &vec) { return vec.mass(); });
}
- similarly, you can use
Map
or an explicit loop inside the Define
, e.g. (not using Map
here just to show how else this could be done):
d.Define("photonmass",\
'auto vecs = ComputeInvariantMass(photon_pt, photon_eta, photon_phi, photon_E); ROOT::RVecD masses(vecs.size()); for (int i = 0; i < vecs.size(); ++i) masses[i] = vecs.mass(); return masses;')
Personally I would go with option 2 and the Define
expression can become InvariantMasses(Make4Vectors(photon_pt, photon_eta, photon_phi, photon_E))
(where I have renamed ComputeInvariantMass
to Make4Vectors
.
Side notes
- another possible implementation of
ComputeInvariantMass
that saves some redundant constructions of PtEtaPhiEVector
s is:
auto vecs = Construct<Vec4_t>(pt, eta, phi, e); // creates an RVec<Vec4_t>
auto pair_idxs = Combinations(vecs, 2); // an RVec<RVec<size_t>> with a list of index pairs that span all pairs
RVec<Vec4_t> pairs;
pairs.reserve((pt.size() * pt.size() - 1) / 2);
for (auto idxs : pair_idxs)
pairs.emplace_back(vecs[idxs[0]] + vecs[idxs[1]]);
return pairs;
- I would like to point out RVec’s helper functions InvariantMass and InvariantMasses in case they can be useful, maybe together with Combinations to get the indexes that form all possible pairs.
Cheers,
Enrico