Logo ROOT  
Reference Guide
RVec.hxx
Go to the documentation of this file.
1// Author: Enrico Guiraud, Enric Tejedor, Danilo Piparo CERN 01/2018
2
3/*************************************************************************
4 * Copyright (C) 1995-2018, Rene Brun and Fons Rademakers. *
5 * All rights reserved. *
6 * *
7 * For the licensing terms see $ROOTSYS/LICENSE. *
8 * For the list of contributors see $ROOTSYS/README/CREDITS. *
9 *************************************************************************/
10
11/**
12 \defgroup vecops VecOps
13*/
14
15#ifndef ROOT_RVEC
16#define ROOT_RVEC
17
18#ifdef _WIN32
19 #define _VECOPS_USE_EXTERN_TEMPLATES false
20#else
21 #define _VECOPS_USE_EXTERN_TEMPLATES true
22#endif
23
26#include <ROOT/RStringView.hxx>
27#include <TError.h> // R__ASSERT
28#include <ROOT/TypeTraits.hxx>
29
30#include <algorithm>
31#include <cmath>
32#include <numeric> // for inner_product
33#include <sstream>
34#include <stdexcept>
35#include <type_traits>
36#include <vector>
37#include <utility>
38#include <tuple>
39
40#define _USE_MATH_DEFINES // enable definition of M_PI
41#ifdef _WIN32
42// cmath does not expose M_PI on windows
43#include <math.h>
44#else
45#include <cmath>
46#endif
47
48#ifdef R__HAS_VDT
49#include <vdt/vdtMath.h>
50#endif
51
52
53namespace ROOT {
54namespace VecOps {
55template<typename T>
56class RVec;
57}
58
59namespace Detail {
60namespace VecOps {
61
62template<typename T>
64
65template <typename... T>
66std::size_t GetVectorsSize(std::string_view id, const RVec<T> &... vs)
67{
68 constexpr const auto nArgs = sizeof...(T);
69 const std::size_t sizes[] = {vs.size()...};
70 if (nArgs > 1) {
71 for (auto i = 1UL; i < nArgs; i++) {
72 if (sizes[0] == sizes[i])
73 continue;
74 std::string msg(id);
75 msg += ": input RVec instances have different lengths!";
76 throw std::runtime_error(msg);
77 }
78 }
79 return sizes[0];
80}
81
82template <typename F, typename... T>
83auto MapImpl(F &&f, const RVec<T> &... vs) -> RVec<decltype(f(vs[0]...))>
84{
85 const auto size = GetVectorsSize("Map", vs...);
86 RVec<decltype(f(vs[0]...))> ret(size);
87
88 for (auto i = 0UL; i < size; i++)
89 ret[i] = f(vs[i]...);
90
91 return ret;
92}
93
94template <typename Tuple_t, std::size_t... Is>
95auto MapFromTuple(Tuple_t &&t, std::index_sequence<Is...>)
96 -> decltype(MapImpl(std::get<std::tuple_size<Tuple_t>::value - 1>(t), std::get<Is>(t)...))
97{
98 constexpr const auto tupleSizeM1 = std::tuple_size<Tuple_t>::value - 1;
99 return MapImpl(std::get<tupleSizeM1>(t), std::get<Is>(t)...);
100}
101
102}
103}
104
105namespace Internal {
106namespace VecOps {
107
108// We use this helper to workaround a limitation of compilers such as
109// gcc 4.8 amd clang on osx 10.14 for which std::vector<bool>::emplace_back
110// is not defined.
111template <typename T, typename... Args>
112void EmplaceBack(T &v, Args &&... args)
113{
114 v.emplace_back(std::forward<Args>(args)...);
115}
116
117template <typename... Args>
118void EmplaceBack(std::vector<bool> &v, Args &&... args)
119{
120 v.push_back(std::forward<Args>(args)...);
121}
122
123} // End of VecOps NS
124} // End of Internal NS
125
126namespace VecOps {
127// clang-format off
128/**
129\class ROOT::VecOps::RVec
130\ingroup vecops
131\brief A "std::vector"-like collection of values implementing handy operation to analyse them
132\tparam T The type of the contained objects
133
134A RVec is a container designed to make analysis of values' collections fast and easy.
135Its storage is contiguous in memory and its interface is designed such to resemble to the one
136of the stl vector. In addition the interface features methods and external functions to ease
137the manipulation and analysis of the data in the RVec.
138
139\htmlonly
140<a href="https://doi.org/10.5281/zenodo.1253756"><img src="https://zenodo.org/badge/DOI/10.5281/zenodo.1253756.svg" alt="DOI"></a>
141\endhtmlonly
142
143## Table of Contents
144- [Example](#example)
145- [Owning and adopting memory](#owningandadoptingmemory)
146- [Sorting and manipulation of indices](#sorting)
147- [Usage in combination with RDataFrame](#usagetdataframe)
148- [Reference for the RVec class](#RVecdoxyref)
149
150Also see the [reference for RVec helper functions](https://root.cern/doc/master/namespaceROOT_1_1VecOps.html).
151
152## <a name="example"></a>Example
153Suppose to have an event featuring a collection of muons with a certain pseudorapidity,
154momentum and charge, e.g.:
155~~~{.cpp}
156std::vector<short> mu_charge {1, 1, -1, -1, -1, 1, 1, -1};
157std::vector<float> mu_pt {56, 45, 32, 24, 12, 8, 7, 6.2};
158std::vector<float> mu_eta {3.1, -.2, -1.1, 1, 4.1, 1.6, 2.4, -.5};
159~~~
160Suppose you want to extract the transverse momenta of the muons satisfying certain
161criteria, for example consider only negatively charged muons with a pseudorapidity
162smaller or equal to 2 and with a transverse momentum greater than 10 GeV.
163Such a selection would require, among the other things, the management of an explicit
164loop, for example:
165~~~{.cpp}
166std::vector<float> goodMuons_pt;
167const auto size = mu_charge.size();
168for (size_t i=0; i < size; ++i) {
169 if (mu_pt[i] > 10 && abs(mu_eta[i]) <= 2. && mu_charge[i] == -1) {
170 goodMuons_pt.emplace_back(mu_pt[i]);
171 }
172}
173~~~
174These operations become straightforward with RVec - we just need to *write what
175we mean*:
176~~~{.cpp}
177auto goodMuons_pt = mu_pt[ (mu_pt > 10.f && abs(mu_eta) <= 2.f && mu_charge == -1) ]
178~~~
179Now the clean collection of transverse momenta can be used within the rest of the data analysis, for
180example to fill a histogram.
181
182## <a name="owningandadoptingmemory"></a>Owning and adopting memory
183RVec has contiguous memory associated to it. It can own it or simply adopt it. In the latter case,
184it can be constructed with the address of the memory associated to it and its length. For example:
185~~~{.cpp}
186std::vector<int> myStlVec {1,2,3};
187RVec<int> myRVec(myStlVec.data(), myStlVec.size());
188~~~
189In this case, the memory associated to myStlVec and myRVec is the same, myRVec simply "adopted it".
190If any method which implies a re-allocation is called, e.g. *emplace_back* or *resize*, the adopted
191memory is released and new one is allocated. The previous content is copied in the new memory and
192preserved.
193
194## <a name="#sorting"></a>Sorting and manipulation of indices
195
196### Sorting
197RVec complies to the STL interfaces when it comes to iterations. As a result, standard algorithms
198can be used, for example sorting:
199~~~{.cpp}
200RVec<double> v{6., 4., 5.};
201std::sort(v.begin(), v.end());
202~~~
203
204For convinience, helpers are provided too:
205~~~{.cpp}
206auto sorted_v = Sort(v);
207auto reversed_v = Reverse(v);
208~~~
209
210### Manipulation of indices
211
212It is also possible to manipulated the RVecs acting on their indices. For example,
213the following syntax
214~~~{.cpp}
215RVec<double> v0 {9., 7., 8.};
216auto v1 = Take(v0, {1, 2, 0});
217~~~
218will yield a new RVec<double> the content of which is the first, second and zeroth element of
219v0, i.e. `{7., 8., 9.}`.
220
221The `Argsort` helper extracts the indices which order the content of a `RVec`. For example,
222this snippet accomplish in a more expressive way what we just achieved:
223~~~{.cpp}
224auto v1_indices = Argsort(v0); // The content of v1_indices is {1, 2, 0}.
225v1 = Take(v0, v1_indices);
226~~~
227
228The `Take` utility allows to extract portions of the `RVec`. The content to be *taken*
229can be specified with an `RVec` of indices or an integer. If the integer is negative,
230elements will be picked starting from the end of the container:
231~~~{.cpp}
232RVec<float> vf {1.f, 2.f, 3.f, 4.f};
233auto vf_1 = Take(vf, {1, 3}); // The content is {2.f, 4.f}
234auto vf_2 = Take(vf, 2); // The content is {1.f, 2.f}
235auto vf_3 = Take(vf, -3); // The content is {2.f, 3.f, 4.f}
236~~~
237
238## <a name="usagetdataframe"></a>Usage in combination with RDataFrame
239RDataFrame leverages internally RVecs. Suppose to have a dataset stored in a
240TTree which holds these columns (here we choose C arrays to represent the
241collections, they could be as well std::vector instances):
242~~~{.bash}
243 nPart "nPart/I" An integer representing the number of particles
244 px "px[nPart]/D" The C array of the particles' x component of the momentum
245 py "py[nPart]/D" The C array of the particles' y component of the momentum
246 E "E[nPart]/D" The C array of the particles' Energy
247~~~
248Suppose you'd like to plot in a histogram the transverse momenta of all particles
249for which the energy is greater than 200 MeV.
250The code required would just be:
251~~~{.cpp}
252RDataFrame d("mytree", "myfile.root");
253using doubles = RVec<double>;
254auto cutPt = [](doubles &pxs, doubles &pys, doubles &Es) {
255 auto all_pts = sqrt(pxs * pxs + pys * pys);
256 auto good_pts = all_pts[Es > 200.];
257 return good_pts;
258 };
259
260auto hpt = d.Define("pt", cutPt, {"px", "py", "E"})
261 .Histo1D("pt");
262hpt->Draw();
263~~~
264And if you'd like to express your selection as a string:
265~~~{.cpp}
266RDataFrame d("mytree", "myfile.root");
267auto hpt = d.Define("pt", "sqrt(pxs * pxs + pys * pys)[E>200]")
268 .Histo1D("pt");
269hpt->Draw();
270~~~
271<a name="RVecdoxyref"></a>
272**/
273// clang-format on
274template <typename T>
275class RVec {
276 // Here we check if T is a bool. This is done in order to decide what type
277 // to use as a storage. If T is anything but bool, we use a vector<T, RAdoptAllocator<T>>.
278 // If T is a bool, we opt for a plain vector<bool> otherwise we'll not be able
279 // to write the data type given the shortcomings of TCollectionProxy design.
280 static constexpr const auto IsVecBool = std::is_same<bool, T>::value;
281public:
282 using Impl_t = typename std::conditional<IsVecBool, std::vector<bool>, std::vector<T, ::ROOT::Detail::VecOps::RAdoptAllocator<T>>>::type;
283 using value_type = typename Impl_t::value_type;
284 using size_type = typename Impl_t::size_type;
285 using difference_type = typename Impl_t::difference_type;
286 using reference = typename Impl_t::reference;
287 using const_reference = typename Impl_t::const_reference;
288 using pointer = typename Impl_t::pointer;
289 using const_pointer = typename Impl_t::const_pointer;
290 // The data_t and const_data_t types are chosen to be void in case T is a bool.
291 // This way we can as elegantly as in the STL return void upon calling the data() method.
294 using iterator = typename Impl_t::iterator;
295 using const_iterator = typename Impl_t::const_iterator;
296 using reverse_iterator = typename Impl_t::reverse_iterator;
297 using const_reverse_iterator = typename Impl_t::const_reverse_iterator;
298
299private:
301
302public:
303 // constructors
304 RVec() {}
305
306 explicit RVec(size_type count) : fData(count) {}
307
308 RVec(size_type count, const T &value) : fData(count, value) {}
309
310 RVec(const RVec<T> &v) : fData(v.fData) {}
311
312 RVec(RVec<T> &&v) : fData(std::move(v.fData)) {}
313
314 RVec(const std::vector<T> &v) : fData(v.cbegin(), v.cend()) {}
315
316 RVec(pointer p, size_type n) : fData(n, T(), ROOT::Detail::VecOps::RAdoptAllocator<T>(p)) {}
317
318 template <class InputIt>
319 RVec(InputIt first, InputIt last) : fData(first, last) {}
320
321 RVec(std::initializer_list<T> init) : fData(init) {}
322
323 // assignment
325 {
326 fData = v.fData;
327 return *this;
328 }
329
331 {
332 std::swap(fData, v.fData);
333 return *this;
334 }
335
336 RVec<T> &operator=(std::initializer_list<T> ilist)
337 {
338 fData = ilist;
339 return *this;
340 }
341
342 // conversion
343 template <typename U, typename = std::enable_if<std::is_convertible<T, U>::value>>
344 operator RVec<U>() const
345 {
346 RVec<U> ret(size());
347 std::copy(begin(), end(), ret.begin());
348 return ret;
349 }
350
351 const Impl_t &AsVector() const { return fData; }
352 Impl_t &AsVector() { return fData; }
353
354 // accessors
355 reference at(size_type pos) { return fData.at(pos); }
356 const_reference at(size_type pos) const { return fData.at(pos); }
357 /// No exception thrown. The user specifies the desired value in case the RVec is shorter than `pos`.
358 value_type at(size_type pos, value_type fallback) { return pos < fData.size() ? fData[pos] : fallback; }
359 /// No exception thrown. The user specifies the desired value in case the RVec is shorter than `pos`.
360 value_type at(size_type pos, value_type fallback) const { return pos < fData.size() ? fData[pos] : fallback; }
361 reference operator[](size_type pos) { return fData[pos]; }
362 const_reference operator[](size_type pos) const { return fData[pos]; }
363
364 template <typename V, typename = std::enable_if<std::is_convertible<V, bool>::value>>
365 RVec operator[](const RVec<V> &conds) const
366 {
367 const size_type n = conds.size();
368
369 if (n != size())
370 throw std::runtime_error("Cannot index RVec with condition vector of different size");
371
372 RVec<T> ret;
373 ret.reserve(n);
374 for (size_type i = 0; i < n; ++i)
375 if (conds[i])
376 ret.emplace_back(fData[i]);
377 return ret;
378 }
379
380 reference front() { return fData.front(); }
381 const_reference front() const { return fData.front(); }
382 reference back() { return fData.back(); }
383 const_reference back() const { return fData.back(); }
384 data_t data() noexcept { return fData.data(); }
385 const_data_t data() const noexcept { return fData.data(); }
386 // iterators
387 iterator begin() noexcept { return fData.begin(); }
388 const_iterator begin() const noexcept { return fData.begin(); }
389 const_iterator cbegin() const noexcept { return fData.cbegin(); }
390 iterator end() noexcept { return fData.end(); }
391 const_iterator end() const noexcept { return fData.end(); }
392 const_iterator cend() const noexcept { return fData.cend(); }
393 reverse_iterator rbegin() noexcept { return fData.rbegin(); }
394 const_reverse_iterator rbegin() const noexcept { return fData.rbegin(); }
395 const_reverse_iterator crbegin() const noexcept { return fData.crbegin(); }
396 reverse_iterator rend() noexcept { return fData.rend(); }
397 const_reverse_iterator rend() const noexcept { return fData.rend(); }
398 const_reverse_iterator crend() const noexcept { return fData.crend(); }
399 // capacity
400 bool empty() const noexcept { return fData.empty(); }
401 size_type size() const noexcept { return fData.size(); }
402 size_type max_size() const noexcept { return fData.size(); }
403 void reserve(size_type new_cap) { fData.reserve(new_cap); }
404 size_type capacity() const noexcept { return fData.capacity(); }
405 void shrink_to_fit() { fData.shrink_to_fit(); };
406 // modifiers
407 void clear() noexcept { fData.clear(); }
408 iterator erase(iterator pos) { return fData.erase(pos); }
409 iterator erase(iterator first, iterator last) { return fData.erase(first, last); }
410 void push_back(T &&value) { fData.push_back(std::forward<T>(value)); }
411 void push_back(const value_type& value) { fData.push_back(value); };
412 template <class... Args>
413 reference emplace_back(Args &&... args)
414 {
415 ROOT::Internal::VecOps::EmplaceBack(fData, std::forward<Args>(args)...);
416 return fData.back();
417 }
418 /// This method is intended only for arithmetic types unlike the std::vector
419 /// corresponding one which is generic.
420 template<typename U = T, typename std::enable_if<std::is_arithmetic<U>::value, int>* = nullptr>
422 {
423 return fData.emplace(pos, value);
424 }
425 void pop_back() { fData.pop_back(); }
426 void resize(size_type count) { fData.resize(count); }
427 void resize(size_type count, const value_type &value) { fData.resize(count, value); }
428 void swap(RVec<T> &other) { std::swap(fData, other.fData); }
429};
430
431///@name RVec Unary Arithmetic Operators
432///@{
433
434#define RVEC_UNARY_OPERATOR(OP) \
435template <typename T> \
436RVec<T> operator OP(const RVec<T> &v) \
437{ \
438 RVec<T> ret(v); \
439 for (auto &x : ret) \
440 x = OP x; \
441return ret; \
442} \
443
448#undef RVEC_UNARY_OPERATOR
449
450///@}
451///@name RVec Binary Arithmetic Operators
452///@{
453
454#define ERROR_MESSAGE(OP) \
455 "Cannot call operator " #OP " on vectors of different sizes."
456
457#define RVEC_BINARY_OPERATOR(OP) \
458template <typename T0, typename T1> \
459auto operator OP(const RVec<T0> &v, const T1 &y) \
460 -> RVec<decltype(v[0] OP y)> \
461{ \
462 RVec<decltype(v[0] OP y)> ret(v.size()); \
463 auto op = [&y](const T0 &x) { return x OP y; }; \
464 std::transform(v.begin(), v.end(), ret.begin(), op); \
465 return ret; \
466} \
467 \
468template <typename T0, typename T1> \
469auto operator OP(const T0 &x, const RVec<T1> &v) \
470 -> RVec<decltype(x OP v[0])> \
471{ \
472 RVec<decltype(x OP v[0])> ret(v.size()); \
473 auto op = [&x](const T1 &y) { return x OP y; }; \
474 std::transform(v.begin(), v.end(), ret.begin(), op); \
475 return ret; \
476} \
477 \
478template <typename T0, typename T1> \
479auto operator OP(const RVec<T0> &v0, const RVec<T1> &v1) \
480 -> RVec<decltype(v0[0] OP v1[0])> \
481{ \
482 if (v0.size() != v1.size()) \
483 throw std::runtime_error(ERROR_MESSAGE(OP)); \
484 \
485 RVec<decltype(v0[0] OP v1[0])> ret(v0.size()); \
486 auto op = [](const T0 &x, const T1 &y) { return x OP y; }; \
487 std::transform(v0.begin(), v0.end(), v1.begin(), ret.begin(), op); \
488 return ret; \
489} \
490
499#undef RVEC_BINARY_OPERATOR
500
501///@}
502///@name RVec Assignment Arithmetic Operators
503///@{
504
505#define RVEC_ASSIGNMENT_OPERATOR(OP) \
506template <typename T0, typename T1> \
507RVec<T0>& operator OP(RVec<T0> &v, const T1 &y) \
508{ \
509 auto op = [&y](T0 &x) { return x OP y; }; \
510 std::transform(v.begin(), v.end(), v.begin(), op); \
511 return v; \
512} \
513 \
514template <typename T0, typename T1> \
515RVec<T0>& operator OP(RVec<T0> &v0, const RVec<T1> &v1) \
516{ \
517 if (v0.size() != v1.size()) \
518 throw std::runtime_error(ERROR_MESSAGE(OP)); \
519 \
520 auto op = [](T0 &x, const T1 &y) { return x OP y; }; \
521 std::transform(v0.begin(), v0.end(), v1.begin(), v0.begin(), op); \
522 return v0; \
523} \
524
535#undef RVEC_ASSIGNMENT_OPERATOR
536
537///@}
538///@name RVec Comparison and Logical Operators
539///@{
540
541#define RVEC_LOGICAL_OPERATOR(OP) \
542template <typename T0, typename T1> \
543auto operator OP(const RVec<T0> &v, const T1 &y) \
544 -> RVec<int> /* avoid std::vector<bool> */ \
545{ \
546 RVec<int> ret(v.size()); \
547 auto op = [y](const T0 &x) -> int { return x OP y; }; \
548 std::transform(v.begin(), v.end(), ret.begin(), op); \
549 return ret; \
550} \
551 \
552template <typename T0, typename T1> \
553auto operator OP(const T0 &x, const RVec<T1> &v) \
554 -> RVec<int> /* avoid std::vector<bool> */ \
555{ \
556 RVec<int> ret(v.size()); \
557 auto op = [x](const T1 &y) -> int { return x OP y; }; \
558 std::transform(v.begin(), v.end(), ret.begin(), op); \
559 return ret; \
560} \
561 \
562template <typename T0, typename T1> \
563auto operator OP(const RVec<T0> &v0, const RVec<T1> &v1) \
564 -> RVec<int> /* avoid std::vector<bool> */ \
565{ \
566 if (v0.size() != v1.size()) \
567 throw std::runtime_error(ERROR_MESSAGE(OP)); \
568 \
569 RVec<int> ret(v0.size()); \
570 auto op = [](const T0 &x, const T1 &y) -> int { return x OP y; }; \
571 std::transform(v0.begin(), v0.end(), v1.begin(), ret.begin(), op); \
572 return ret; \
573} \
574
583#undef RVEC_LOGICAL_OPERATOR
584
585///@}
586///@name RVec Standard Mathematical Functions
587///@{
588
589/// \cond
590template <typename T> struct PromoteTypeImpl;
591
592template <> struct PromoteTypeImpl<float> { using Type = float; };
593template <> struct PromoteTypeImpl<double> { using Type = double; };
594template <> struct PromoteTypeImpl<long double> { using Type = long double; };
595
596template <typename T> struct PromoteTypeImpl { using Type = double; };
597
598template <typename T>
599using PromoteType = typename PromoteTypeImpl<T>::Type;
600
601template <typename U, typename V>
602using PromoteTypes = decltype(PromoteType<U>() + PromoteType<V>());
603
604/// \endcond
605
606#define RVEC_UNARY_FUNCTION(NAME, FUNC) \
607 template <typename T> \
608 RVec<PromoteType<T>> NAME(const RVec<T> &v) \
609 { \
610 RVec<PromoteType<T>> ret(v.size()); \
611 auto f = [](const T &x) { return FUNC(x); }; \
612 std::transform(v.begin(), v.end(), ret.begin(), f); \
613 return ret; \
614 }
615
616#define RVEC_BINARY_FUNCTION(NAME, FUNC) \
617 template <typename T0, typename T1> \
618 RVec<PromoteTypes<T0, T1>> NAME(const T0 &x, const RVec<T1> &v) \
619 { \
620 RVec<PromoteTypes<T0, T1>> ret(v.size()); \
621 auto f = [&x](const T1 &y) { return FUNC(x, y); }; \
622 std::transform(v.begin(), v.end(), ret.begin(), f); \
623 return ret; \
624 } \
625 \
626 template <typename T0, typename T1> \
627 RVec<PromoteTypes<T0, T1>> NAME(const RVec<T0> &v, const T1 &y) \
628 { \
629 RVec<PromoteTypes<T0, T1>> ret(v.size()); \
630 auto f = [&y](const T1 &x) { return FUNC(x, y); }; \
631 std::transform(v.begin(), v.end(), ret.begin(), f); \
632 return ret; \
633 } \
634 \
635 template <typename T0, typename T1> \
636 RVec<PromoteTypes<T0, T1>> NAME(const RVec<T0> &v0, const RVec<T1> &v1) \
637 { \
638 if (v0.size() != v1.size()) \
639 throw std::runtime_error(ERROR_MESSAGE(NAME)); \
640 \
641 RVec<PromoteTypes<T0, T1>> ret(v0.size()); \
642 auto f = [](const T0 &x, const T1 &y) { return FUNC(x, y); }; \
643 std::transform(v0.begin(), v0.end(), v1.begin(), ret.begin(), f); \
644 return ret; \
645 } \
646
647#define RVEC_STD_UNARY_FUNCTION(F) RVEC_UNARY_FUNCTION(F, std::F)
648#define RVEC_STD_BINARY_FUNCTION(F) RVEC_BINARY_FUNCTION(F, std::F)
649
654
658
663
668
676
683
690
695#undef RVEC_STD_UNARY_FUNCTION
696
697///@}
698///@name RVec Fast Mathematical Functions with Vdt
699///@{
700
701#ifdef R__HAS_VDT
702#define RVEC_VDT_UNARY_FUNCTION(F) RVEC_UNARY_FUNCTION(F, vdt::F)
703
704RVEC_VDT_UNARY_FUNCTION(fast_expf)
705RVEC_VDT_UNARY_FUNCTION(fast_logf)
706RVEC_VDT_UNARY_FUNCTION(fast_sinf)
707RVEC_VDT_UNARY_FUNCTION(fast_cosf)
708RVEC_VDT_UNARY_FUNCTION(fast_tanf)
709RVEC_VDT_UNARY_FUNCTION(fast_asinf)
710RVEC_VDT_UNARY_FUNCTION(fast_acosf)
711RVEC_VDT_UNARY_FUNCTION(fast_atanf)
712
713RVEC_VDT_UNARY_FUNCTION(fast_exp)
714RVEC_VDT_UNARY_FUNCTION(fast_log)
715RVEC_VDT_UNARY_FUNCTION(fast_sin)
716RVEC_VDT_UNARY_FUNCTION(fast_cos)
717RVEC_VDT_UNARY_FUNCTION(fast_tan)
718RVEC_VDT_UNARY_FUNCTION(fast_asin)
719RVEC_VDT_UNARY_FUNCTION(fast_acos)
720RVEC_VDT_UNARY_FUNCTION(fast_atan)
721#undef RVEC_VDT_UNARY_FUNCTION
722
723#endif // R__HAS_VDT
724
725#undef RVEC_UNARY_FUNCTION
726
727///@}
728
729/// Inner product
730///
731/// Example code, at the ROOT prompt:
732/// ~~~{.cpp}
733/// using namespace ROOT::VecOps;
734/// RVec<float> v1 {1., 2., 3.};
735/// RVec<float> v2 {4., 5., 6.};
736/// auto v1_dot_v2 = Dot(v1, v2);
737/// v1_dot_v2
738/// // (float) 32.f
739/// ~~~
740template <typename T, typename V>
741auto Dot(const RVec<T> &v0, const RVec<V> &v1) -> decltype(v0[0] * v1[0])
742{
743 if (v0.size() != v1.size())
744 throw std::runtime_error("Cannot compute inner product of vectors of different sizes");
745 return std::inner_product(v0.begin(), v0.end(), v1.begin(), decltype(v0[0] * v1[0])(0));
746}
747
748/// Sum elements of an RVec
749///
750/// Example code, at the ROOT prompt:
751/// ~~~{.cpp}
752/// using namespace ROOT::VecOps;
753/// RVec<float> v {1.f, 2.f, 3.f};
754/// auto v_sum = Sum(v);
755/// v_sum
756/// // (float) 6.f
757/// ~~~
758template <typename T>
759T Sum(const RVec<T> &v)
760{
761 return std::accumulate(v.begin(), v.end(), T(0));
762}
763
764/// Get the mean of the elements of an RVec
765///
766/// The return type is a double precision floating point number.
767/// Example code, at the ROOT prompt:
768/// ~~~{.cpp}
769/// using namespace ROOT::VecOps;
770/// RVec<float> v {1.f, 2.f, 4.f};
771/// auto v_mean = Mean(v);
772/// v_mean
773/// // (double) 2.3333333
774/// ~~~
775template <typename T>
776double Mean(const RVec<T> &v)
777{
778 if (v.empty()) return 0.;
779 return double(Sum(v)) / v.size();
780}
781
782/// Get the greatest element of an RVec
783///
784/// Example code, at the ROOT prompt:
785/// ~~~~{.cpp}
786/// using namespace ROOT::VecOps;
787/// RVec<float> v {1.f, 2.f, 4.f};
788/// auto v_max = Max(v)
789/// v_max
790/// (float) 4.f
791/// ~~~~
792template <typename T>
793T Max(const RVec<T> &v)
794{
795 return *std::max_element(v.begin(), v.end());
796}
797
798/// Get the smallest element of an RVec
799///
800/// Example code, at the ROOT prompt:
801/// ~~~~{.cpp}
802/// using namespace ROOT::VecOps;
803/// RVec<float> v {1.f, 2.f, 4.f};
804/// auto v_min = Min(v)
805/// v_min
806/// (float) 1.f
807/// ~~~~
808template <typename T>
809T Min(const RVec<T> &v)
810{
811 return *std::min_element(v.begin(), v.end());
812}
813
814/// Get the index of the greatest element of an RVec
815/// In case of multiple occurrences of the maximum values,
816/// the index corresponding to the first occurrence is returned.
817///
818/// Example code, at the ROOT prompt:
819/// ~~~~{.cpp}
820/// using namespace ROOT::VecOps;
821/// RVec<float> v {1.f, 2.f, 4.f};
822/// auto v_argmax = ArgMax(v);
823/// v_argmax
824/// // (int) 2
825/// ~~~~
826template <typename T>
827std::size_t ArgMax(const RVec<T> &v)
828{
829 return std::distance(v.begin(), std::max_element(v.begin(), v.end()));
830}
831
832/// Get the index of the smallest element of an RVec
833/// In case of multiple occurrences of the minimum values,
834/// the index corresponding to the first occurrence is returned.
835///
836/// Example code, at the ROOT prompt:
837/// ~~~~{.cpp}
838/// using namespace ROOT::VecOps;
839/// RVec<float> v {1.f, 2.f, 4.f};
840/// auto v_argmin = ArgMin(v);
841/// v_argmin
842/// // (int) 0
843/// ~~~~
844template <typename T>
845std::size_t ArgMin(const RVec<T> &v)
846{
847 return std::distance(v.begin(), std::min_element(v.begin(), v.end()));
848}
849
850/// Get the variance of the elements of an RVec
851///
852/// The return type is a double precision floating point number.
853/// Example code, at the ROOT prompt:
854/// ~~~{.cpp}
855/// using namespace ROOT::VecOps;
856/// RVec<float> v {1.f, 2.f, 4.f};
857/// auto v_var = Var(v);
858/// v_var
859/// // (double) 2.3333333
860/// ~~~
861template <typename T>
862double Var(const RVec<T> &v)
863{
864 const std::size_t size = v.size();
865 if (size < std::size_t(2)) return 0.;
866 T sum_squares(0), squared_sum(0);
867 auto pred = [&sum_squares, &squared_sum](const T& x) {sum_squares+=x*x; squared_sum+=x;};
868 std::for_each(v.begin(), v.end(), pred);
869 squared_sum *= squared_sum;
870 const auto dsize = (double) size;
871 return 1. / (dsize - 1.) * (sum_squares - squared_sum / dsize );
872}
873
874/// Get the standard deviation of the elements of an RVec
875///
876/// The return type is a double precision floating point number.
877/// Example code, at the ROOT prompt:
878/// ~~~{.cpp}
879/// using namespace ROOT::VecOps;
880/// RVec<float> v {1.f, 2.f, 4.f};
881/// auto v_sd = StdDev(v);
882/// v_sd
883/// // (double) 1.5275252
884/// ~~~
885template <typename T>
886double StdDev(const RVec<T> &v)
887{
888 return std::sqrt(Var(v));
889}
890
891/// Create new collection applying a callable to the elements of the input collection
892///
893/// Example code, at the ROOT prompt:
894/// ~~~{.cpp}
895/// using namespace ROOT::VecOps;
896/// RVec<float> v {1.f, 2.f, 4.f};
897/// auto v_square = Map(v, [](float f){return f* 2.f;});
898/// v_square
899/// // (ROOT::VecOps::RVec<float> &) { 2.00000f, 4.00000f, 8.00000f }
900///
901/// RVec<float> x({1.f, 2.f, 3.f});
902/// RVec<float> y({4.f, 5.f, 6.f});
903/// RVec<float> z({7.f, 8.f, 9.f});
904/// auto mod = [](float x, float y, float z) { return sqrt(x * x + y * y + z * z); };
905/// auto v_mod = Map(x, y, z, mod);
906/// v_mod
907/// // (ROOT::VecOps::RVec<float> &) { 8.12404f, 9.64365f, 11.2250f }
908/// ~~~
909template <typename... Args>
910auto Map(Args &&... args)
911 -> decltype(ROOT::Detail::VecOps::MapFromTuple(std::forward_as_tuple(args...),
912 std::make_index_sequence<sizeof...(args) - 1>()))
913{
914 /*
915 Here the strategy in order to generalise the previous implementation of Map, i.e.
916 `RVec Map(RVec, F)`, here we need to move the last parameter of the pack in first
917 position in order to be able to invoke the Map function with automatic type deduction.
918 This is achieved in two steps:
919 1. Forward as tuple the pack to MapFromTuple
920 2. Invoke the MapImpl helper which has the signature `template<...T, F> RVec MapImpl(F &&f, RVec<T>...)`
921 NOTA BENE: the signature is very heavy but it is one of the lightest ways to manage in C++11
922 to build the return type based on the template args.
923 */
924 return ROOT::Detail::VecOps::MapFromTuple(std::forward_as_tuple(args...),
925 std::make_index_sequence<sizeof...(args) - 1>());
926}
927
928/// Create a new collection with the elements passing the filter expressed by the predicate
929///
930/// Example code, at the ROOT prompt:
931/// ~~~{.cpp}
932/// using namespace ROOT::VecOps;
933/// RVec<int> v {1, 2, 4};
934/// auto v_even = Filter(v, [](int i){return 0 == i%2;});
935/// v_even
936/// // (ROOT::VecOps::RVec<int> &) { 2, 4 }
937/// ~~~
938template <typename T, typename F>
940{
941 const auto thisSize = v.size();
942 RVec<T> w;
943 w.reserve(thisSize);
944 for (auto &&val : v) {
945 if (f(val))
946 w.emplace_back(val);
947 }
948 return w;
949}
950
951/// Return true if any of the elements equates to true, return false otherwise.
952///
953/// Example code, at the ROOT prompt:
954/// ~~~{.cpp}
955/// using namespace ROOT::VecOps;
956/// RVec<int> v {0, 1, 0};
957/// auto anyTrue = Any(v);
958/// anyTrue
959/// // (bool) true
960/// ~~~
961template <typename T>
962auto Any(const RVec<T> &v) -> decltype(v[0] == true)
963{
964 for (auto &&e : v)
965 if (e == true)
966 return true;
967 return false;
968}
969
970/// Return true if all of the elements equate to true, return false otherwise.
971///
972/// Example code, at the ROOT prompt:
973/// ~~~{.cpp}
974/// using namespace ROOT::VecOps;
975/// RVec<int> v {0, 1, 0};
976/// auto allTrue = All(v);
977/// allTrue
978/// // (bool) false
979/// ~~~
980template <typename T>
981auto All(const RVec<T> &v) -> decltype(v[0] == false)
982{
983 for (auto &&e : v)
984 if (e == false)
985 return false;
986 return true;
987}
988
989template <typename T>
990void swap(RVec<T> &lhs, RVec<T> &rhs)
991{
992 lhs.swap(rhs);
993}
994
995/// Return an RVec of indices that sort the input RVec
996///
997/// Example code, at the ROOT prompt:
998/// ~~~{.cpp}
999/// using namespace ROOT::VecOps;
1000/// RVec<double> v {2., 3., 1.};
1001/// auto sortIndices = Argsort(v);
1002/// sortIndices
1003/// // (ROOT::VecOps::RVec<unsigned long> &) { 2, 0, 1 }
1004/// ~~~
1005template <typename T>
1007{
1008 using size_type = typename RVec<T>::size_type;
1009 RVec<size_type> i(v.size());
1010 std::iota(i.begin(), i.end(), 0);
1011 std::sort(i.begin(), i.end(), [&v](size_type i1, size_type i2) { return v[i1] < v[i2]; });
1012 return i;
1013}
1014
1015/// Return elements of a vector at given indices
1016///
1017/// Example code, at the ROOT prompt:
1018/// ~~~{.cpp}
1019/// using namespace ROOT::VecOps;
1020/// RVec<double> v {2., 3., 1.};
1021/// auto vTaken = Take(v, {0,2});
1022/// vTaken
1023/// // (ROOT::VecOps::RVec<double>) { 2.0000000, 1.0000000 }
1024/// ~~~
1025template <typename T>
1026RVec<T> Take(const RVec<T> &v, const RVec<typename RVec<T>::size_type> &i)
1027{
1028 using size_type = typename RVec<T>::size_type;
1029 const size_type isize = i.size();
1030 RVec<T> r(isize);
1031 for (size_type k = 0; k < isize; k++)
1032 r[k] = v[i[k]];
1033 return r;
1034}
1035
1036/// Return first or last `n` elements of an RVec
1037///
1038/// if `n > 0` and last elements if `n < 0`.
1039///
1040/// Example code, at the ROOT prompt:
1041/// ~~~{.cpp}
1042/// using namespace ROOT::VecOps;
1043/// RVec<double> v {2., 3., 1.};
1044/// auto firstTwo = Take(v, 2);
1045/// firstTwo
1046/// // (ROOT::VecOps::RVec<double>) { 2.0000000, 3.0000000 }
1047/// auto lastOne = Take(v, -1);
1048/// lastOne
1049/// // (ROOT::VecOps::RVec<double>) { 1.0000000 }
1050/// ~~~
1051template <typename T>
1052RVec<T> Take(const RVec<T> &v, const int n)
1053{
1054 using size_type = typename RVec<T>::size_type;
1055 const size_type size = v.size();
1056 const size_type absn = std::abs(n);
1057 if (absn > size) {
1058 std::stringstream ss;
1059 ss << "Try to take " << absn << " elements but vector has only size " << size << ".";
1060 throw std::runtime_error(ss.str());
1061 }
1062 RVec<T> r(absn);
1063 if (n < 0) {
1064 for (size_type k = 0; k < absn; k++)
1065 r[k] = v[size - absn + k];
1066 } else {
1067 for (size_type k = 0; k < absn; k++)
1068 r[k] = v[k];
1069 }
1070 return r;
1071}
1072
1073/// Return copy of reversed vector
1074///
1075/// Example code, at the ROOT prompt:
1076/// ~~~{.cpp}
1077/// using namespace ROOT::VecOps;
1078/// RVec<double> v {2., 3., 1.};
1079/// auto v_reverse = Reverse(v);
1080/// v_reverse
1081/// // (ROOT::VecOps::RVec<double>) { 1.0000000, 3.0000000, 2.0000000 }
1082/// ~~~
1083template <typename T>
1085{
1086 RVec<T> r(v);
1087 std::reverse(r.begin(), r.end());
1088 return r;
1089}
1090
1091/// Return copy of RVec with elements sorted in ascending order
1092///
1093/// This helper is different from ArgSort since it does not return an RVec of indices,
1094/// but an RVec of values.
1095///
1096/// Example code, at the ROOT prompt:
1097/// ~~~{.cpp}
1098/// using namespace ROOT::VecOps;
1099/// RVec<double> v {2., 3., 1.};
1100/// auto v_sorted = Sort(v);
1101/// v_sorted
1102/// // (ROOT::VecOps::RVec<double>) { 1.0000000, 2.0000000, 3.0000000 }
1103/// ~~~
1104template <typename T>
1106{
1107 RVec<T> r(v);
1108 std::sort(r.begin(), r.end());
1109 return r;
1110}
1111
1112/// Return copy of RVec with elements sorted based on a comparison operator
1113///
1114/// The comparison operator has to fullfill the same requirements of the
1115/// predicate of by std::sort.
1116///
1117///
1118/// This helper is different from ArgSort since it does not return an RVec of indices,
1119/// but an RVec of values.
1120///
1121/// Example code, at the ROOT prompt:
1122/// ~~~{.cpp}
1123/// using namespace ROOT::VecOps;
1124/// RVec<double> v {2., 3., 1.};
1125/// auto v_sorted = Sort(v, [](double x, double y) {return 1/x < 1/y;});
1126/// v_sorted
1127/// // (ROOT::VecOps::RVec<double>) { 3.0000000, 2.0000000, 1.0000000 }
1128/// ~~~
1129template <typename T, typename Compare>
1131{
1132 RVec<T> r(v);
1133 std::sort(r.begin(), r.end(), std::forward<Compare>(c));
1134 return r;
1135}
1136
1137/// Return the indices that represent all combinations of the elements of two
1138/// RVecs.
1139///
1140/// The type of the return value is an RVec of two RVecs containing indices.
1141///
1142/// Example code, at the ROOT prompt:
1143/// ~~~{.cpp}
1144/// using namespace ROOT::VecOps;
1145/// auto comb_idx = Combinations(3, 2);
1146/// comb_idx
1147/// // (ROOT::VecOps::RVec<ROOT::VecOps::RVec<ROOT::VecOps::RVec<double>::size_type> >) { { 0, 0, 1, 1, 2, 2 }, { 0, 1, 0, 1, 0, 1 } }
1148/// ~~~
1149inline RVec<RVec<std::size_t>> Combinations(const std::size_t size1, const std::size_t size2)
1150{
1151 using size_type = std::size_t;
1153 r[0].resize(size1*size2);
1154 r[1].resize(size1*size2);
1155 size_type c = 0;
1156 for(size_type i=0; i<size1; i++) {
1157 for(size_type j=0; j<size2; j++) {
1158 r[0][c] = i;
1159 r[1][c] = j;
1160 c++;
1161 }
1162 }
1163 return r;
1164}
1165
1166/// Return the indices that represent all combinations of the elements of two
1167/// RVecs.
1168///
1169/// The type of the return value is an RVec of two RVecs containing indices.
1170///
1171/// Example code, at the ROOT prompt:
1172/// ~~~{.cpp}
1173/// using namespace ROOT::VecOps;
1174/// RVec<double> v1 {1., 2., 3.};
1175/// RVec<double> v2 {-4., -5.};
1176/// auto comb_idx = Combinations(v1, v2);
1177/// comb_idx
1178/// // (ROOT::VecOps::RVec<ROOT::VecOps::RVec<ROOT::VecOps::RVec<double>::size_type> >) { { 0, 0, 1, 1, 2, 2 }, { 0, 1,
1179/// 0, 1, 0, 1 } }
1180/// ~~~
1181template <typename T1, typename T2>
1183{
1184 return Combinations(v1.size(), v2.size());
1185}
1186
1187/// Return the indices that represent all unique combinations of the
1188/// elements of a given RVec.
1189///
1190/// ~~~{.cpp}
1191/// using namespace ROOT::VecOps;
1192/// RVec<double> v {1., 2., 3., 4.};
1193/// auto v_1 = Combinations(v, 1);
1194/// v_1
1195/// (ROOT::VecOps::RVec<ROOT::VecOps::RVec<ROOT::VecOps::RVec<double>::size_type> >) { { 0, 1, 2, 3 } }
1196/// auto v_2 = Combinations(v, 2);
1197/// auto v_2
1198/// (ROOT::VecOps::RVec<ROOT::VecOps::RVec<ROOT::VecOps::RVec<double>::size_type> >) { { 0, 0, 0, 1, 1, 2 }, { 1, 2, 3, 2, 3, 3 } }
1199/// auto v_3 = Combinations(v, 3);
1200/// v_3
1201/// (ROOT::VecOps::RVec<ROOT::VecOps::RVec<ROOT::VecOps::RVec<double>::size_type> >) { { 0, 0, 0, 1 }, { 1, 1, 2, 2 }, { 2, 3, 3, 3 } }
1202/// auto v_4 = Combinations(v, 4);
1203/// v_4
1204/// (ROOT::VecOps::RVec<ROOT::VecOps::RVec<ROOT::VecOps::RVec<double>::size_type> >) { { 0 }, { 1 }, { 2 }, { 3 } }
1205/// ~~~
1206template <typename T>
1208{
1209 using size_type = typename RVec<T>::size_type;
1210 const size_type s = v.size();
1211 if (n > s) {
1212 std::stringstream ss;
1213 ss << "Cannot make unique combinations of size " << n << " from vector of size " << s << ".";
1214 throw std::runtime_error(ss.str());
1215 }
1217 for(size_type k=0; k<s; k++)
1218 indices[k] = k;
1220 for(size_type k=0; k<n; k++)
1221 c[k].emplace_back(indices[k]);
1222 while (true) {
1223 bool run_through = true;
1224 long i = n - 1;
1225 for (; i>=0; i--) {
1226 if (indices[i] != i + s - n){
1227 run_through = false;
1228 break;
1229 }
1230 }
1231 if (run_through) {
1232 return c;
1233 }
1234 indices[i]++;
1235 for (long j=i+1; j<(long)n; j++)
1236 indices[j] = indices[j-1] + 1;
1237 for(size_type k=0; k<n; k++)
1238 c[k].emplace_back(indices[k]);
1239 }
1240}
1241
1242/// Return the indices of the elements which are not zero
1243///
1244/// Example code, at the ROOT prompt:
1245/// ~~~{.cpp}
1246/// using namespace ROOT::VecOps;
1247/// RVec<double> v {2., 0., 3., 0., 1.};
1248/// auto nonzero_idx = Nonzero(v);
1249/// nonzero_idx
1250/// // (ROOT::VecOps::RVec<ROOT::VecOps::RVec<double>::size_type>) { 0, 2, 4 }
1251/// ~~~
1252template <typename T>
1254{
1255 using size_type = typename RVec<T>::size_type;
1257 const auto size = v.size();
1258 r.reserve(size);
1259 for(size_type i=0; i<size; i++) {
1260 if(v[i] != 0) {
1261 r.emplace_back(i);
1262 }
1263 }
1264 return r;
1265}
1266
1267/// Return the intersection of elements of two RVecs.
1268///
1269/// Each element of v1 is looked up in v2 and added to the returned vector if
1270/// found. Following, the order of v1 is preserved. If v2 is already sorted, the
1271/// optional argument v2_is_sorted can be used to toggle of the internal sorting
1272/// step, therewith optimising runtime.
1273///
1274/// Example code, at the ROOT prompt:
1275/// ~~~{.cpp}
1276/// using namespace ROOT::VecOps;
1277/// RVec<double> v1 {1., 2., 3.};
1278/// RVec<double> v2 {-4., -5., 2., 1.};
1279/// auto v1_intersect_v2 = Intersect(v1, v2);
1280/// v1_intersect_v2
1281/// // (ROOT::VecOps::RVec<double>) { 1.0000000, 2.0000000 }
1282/// ~~~
1283template <typename T>
1284RVec<T> Intersect(const RVec<T>& v1, const RVec<T>& v2, bool v2_is_sorted = false)
1285{
1286 RVec<T> v2_sorted;
1287 if (!v2_is_sorted) v2_sorted = Sort(v2);
1288 const auto v2_begin = v2_is_sorted ? v2.begin() : v2_sorted.begin();
1289 const auto v2_end = v2_is_sorted ? v2.end() : v2_sorted.end();
1290 RVec<T> r;
1291 const auto size = v1.size();
1292 r.reserve(size);
1293 using size_type = typename RVec<T>::size_type;
1294 for(size_type i=0; i<size; i++) {
1295 if (std::binary_search(v2_begin, v2_end, v1[i])) {
1296 r.emplace_back(v1[i]);
1297 }
1298 }
1299 return r;
1300}
1301
1302/// Return the elements of v1 if the condition c is true and v2 if the
1303/// condition c is false.
1304///
1305/// Example code, at the ROOT prompt:
1306/// ~~~{.cpp}
1307/// using namespace ROOT::VecOps;
1308/// RVec<double> v1 {1., 2., 3.};
1309/// RVec<double> v2 {-1., -2., -3.};
1310/// auto c = v1 > 1;
1311/// c
1312/// // (ROOT::VecOps::RVec<int> &) { 0, 1, 1 }
1313/// auto if_c_v1_else_v2 = Where(c, v1, v2);
1314/// if_c_v1_else_v2
1315/// // (ROOT::VecOps::RVec<double> &) { -1.0000000, 2.0000000, 3.0000000 }
1316/// ~~~
1317template <typename T>
1318RVec<T> Where(const RVec<int>& c, const RVec<T>& v1, const RVec<T>& v2)
1319{
1320 using size_type = typename RVec<T>::size_type;
1321 const size_type size = c.size();
1322 RVec<T> r;
1323 r.reserve(size);
1324 for (size_type i=0; i<size; i++) {
1325 r.emplace_back(c[i] != 0 ? v1[i] : v2[i]);
1326 }
1327 return r;
1328}
1329
1330/// Return the elements of v1 if the condition c is true and sets the value v2
1331/// if the condition c is false.
1332///
1333/// Example code, at the ROOT prompt:
1334/// ~~~{.cpp}
1335/// using namespace ROOT::VecOps;
1336/// RVec<double> v1 {1., 2., 3.};
1337/// double v2 = 4.;
1338/// auto c = v1 > 1;
1339/// c
1340/// // (ROOT::VecOps::RVec<int> &) { 0, 1, 1 }
1341/// auto if_c_v1_else_v2 = Where(c, v1, v2);
1342/// if_c_v1_else_v2
1343/// // (ROOT::VecOps::RVec<double>) { 4.0000000, 2.0000000, 3.0000000 }
1344/// ~~~
1345template <typename T>
1347{
1348 using size_type = typename RVec<T>::size_type;
1349 const size_type size = c.size();
1350 RVec<T> r;
1351 r.reserve(size);
1352 for (size_type i=0; i<size; i++) {
1353 r.emplace_back(c[i] != 0 ? v1[i] : v2);
1354 }
1355 return r;
1356}
1357
1358/// Return the elements of v2 if the condition c is false and sets the value v1
1359/// if the condition c is true.
1360///
1361/// Example code, at the ROOT prompt:
1362/// ~~~{.cpp}
1363/// using namespace ROOT::VecOps;
1364/// double v1 = 4.;
1365/// RVec<double> v2 {1., 2., 3.};
1366/// auto c = v2 > 1;
1367/// c
1368/// // (ROOT::VecOps::RVec<int> &) { 0, 1, 1 }
1369/// auto if_c_v1_else_v2 = Where(c, v1, v2);
1370/// if_c_v1_else_v2
1371/// // (ROOT::VecOps::RVec<double>) { 1.0000000, 4.0000000, 4.0000000 }
1372/// ~~~
1373template <typename T>
1375{
1376 using size_type = typename RVec<T>::size_type;
1377 const size_type size = c.size();
1378 RVec<T> r;
1379 r.reserve(size);
1380 for (size_type i=0; i<size; i++) {
1381 r.emplace_back(c[i] != 0 ? v1 : v2[i]);
1382 }
1383 return r;
1384}
1385
1386/// Return a vector with the value v2 if the condition c is false and sets the
1387/// value v1 if the condition c is true.
1388///
1389/// Example code, at the ROOT prompt:
1390/// ~~~{.cpp}
1391/// using namespace ROOT::VecOps;
1392/// double v1 = 4.;
1393/// double v2 = 2.;
1394/// RVec<int> c {0, 1, 1};
1395/// auto if_c_v1_else_v2 = Where(c, v1, v2);
1396/// if_c_v1_else_v2
1397/// // (ROOT::VecOps::RVec<double>) { 2.0000000, 4.0000000, 4.0000000 }
1398/// ~~~
1399template <typename T>
1401{
1402 using size_type = typename RVec<T>::size_type;
1403 const size_type size = c.size();
1404 RVec<T> r;
1405 r.reserve(size);
1406 for (size_type i=0; i<size; i++) {
1407 r.emplace_back(c[i] != 0 ? v1 : v2);
1408 }
1409 return r;
1410}
1411
1412/// Return the concatenation of two RVecs.
1413///
1414/// Example code, at the ROOT prompt:
1415/// ~~~{.cpp}
1416/// using namespace ROOT::VecOps;
1417/// RVec<float> rvf {0.f, 1.f, 2.f};
1418/// RVec<int> rvi {7, 8, 9};
1419/// Concatenate(rvf, rvi);
1420/// // (ROOT::VecOps::RVec<float>) { 2.0000000, 4.0000000, 4.0000000 }
1421/// ~~~
1424{
1425 RVec<Common_t> res;
1426 res.reserve(v0.size() + v1.size());
1427 auto &resAsVect = res.AsVector();
1428 auto &v0AsVect = v0.AsVector();
1429 auto &v1AsVect = v1.AsVector();
1430 resAsVect.insert(resAsVect.begin(), v0AsVect.begin(), v0AsVect.end());
1431 resAsVect.insert(resAsVect.end(), v1AsVect.begin(), v1AsVect.end());
1432 return res;
1433}
1434
1435/// Return the angle difference \f$\Delta \phi\f$ of two scalars.
1436///
1437/// The function computes the closest angle from v1 to v2 with sign and is
1438/// therefore in the range \f$[-\pi, \pi]\f$.
1439/// The computation is done per default in radians \f$c = \pi\f$ but can be switched
1440/// to degrees \f$c = 180\f$.
1441template <typename T>
1442T DeltaPhi(T v1, T v2, const T c = M_PI)
1443{
1444 static_assert(std::is_floating_point<T>::value,
1445 "DeltaPhi must be called with floating point values.");
1446 auto r = std::fmod(v2 - v1, 2.0 * c);
1447 if (r < -c) {
1448 r += 2.0 * c;
1449 }
1450 else if (r > c) {
1451 r -= 2.0 * c;
1452 }
1453 return r;
1454}
1455
1456/// Return the angle difference \f$\Delta \phi\f$ in radians of two vectors.
1457///
1458/// The function computes the closest angle from v1 to v2 with sign and is
1459/// therefore in the range \f$[-\pi, \pi]\f$.
1460/// The computation is done per default in radians \f$c = \pi\f$ but can be switched
1461/// to degrees \f$c = 180\f$.
1462template <typename T>
1463RVec<T> DeltaPhi(const RVec<T>& v1, const RVec<T>& v2, const T c = M_PI)
1464{
1465 using size_type = typename RVec<T>::size_type;
1466 const size_type size = v1.size();
1467 auto r = RVec<T>(size);
1468 for (size_type i = 0; i < size; i++) {
1469 r[i] = DeltaPhi(v1[i], v2[i], c);
1470 }
1471 return r;
1472}
1473
1474/// Return the angle difference \f$\Delta \phi\f$ in radians of a vector and a scalar.
1475///
1476/// The function computes the closest angle from v1 to v2 with sign and is
1477/// therefore in the range \f$[-\pi, \pi]\f$.
1478/// The computation is done per default in radians \f$c = \pi\f$ but can be switched
1479/// to degrees \f$c = 180\f$.
1480template <typename T>
1481RVec<T> DeltaPhi(const RVec<T>& v1, T v2, const T c = M_PI)
1482{
1483 using size_type = typename RVec<T>::size_type;
1484 const size_type size = v1.size();
1485 auto r = RVec<T>(size);
1486 for (size_type i = 0; i < size; i++) {
1487 r[i] = DeltaPhi(v1[i], v2, c);
1488 }
1489 return r;
1490}
1491
1492/// Return the angle difference \f$\Delta \phi\f$ in radians of a scalar and a vector.
1493///
1494/// The function computes the closest angle from v1 to v2 with sign and is
1495/// therefore in the range \f$[-\pi, \pi]\f$.
1496/// The computation is done per default in radians \f$c = \pi\f$ but can be switched
1497/// to degrees \f$c = 180\f$.
1498template <typename T>
1499RVec<T> DeltaPhi(T v1, const RVec<T>& v2, const T c = M_PI)
1500{
1501 using size_type = typename RVec<T>::size_type;
1502 const size_type size = v2.size();
1503 auto r = RVec<T>(size);
1504 for (size_type i = 0; i < size; i++) {
1505 r[i] = DeltaPhi(v1, v2[i], c);
1506 }
1507 return r;
1508}
1509
1510/// Return the square of the distance on the \f$\eta\f$-\f$\phi\f$ plane (\f$\Delta R\f$) from
1511/// the collections eta1, eta2, phi1 and phi2.
1512///
1513/// The function computes \f$\Delta R^2 = (\eta_1 - \eta_2)^2 + (\phi_1 - \phi_2)^2\f$
1514/// of the given collections eta1, eta2, phi1 and phi2. The angle \f$\phi\f$ can
1515/// be set to radian or degrees using the optional argument c, see the documentation
1516/// of the DeltaPhi helper.
1517template <typename T>
1518RVec<T> DeltaR2(const RVec<T>& eta1, const RVec<T>& eta2, const RVec<T>& phi1, const RVec<T>& phi2, const T c = M_PI)
1519{
1520 const auto dphi = DeltaPhi(phi1, phi2, c);
1521 return (eta1 - eta2) * (eta1 - eta2) + dphi * dphi;
1522}
1523
1524/// Return the distance on the \f$\eta\f$-\f$\phi\f$ plane (\f$\Delta R\f$) from
1525/// the collections eta1, eta2, phi1 and phi2.
1526///
1527/// The function computes \f$\Delta R = \sqrt{(\eta_1 - \eta_2)^2 + (\phi_1 - \phi_2)^2}\f$
1528/// of the given collections eta1, eta2, phi1 and phi2. The angle \f$\phi\f$ can
1529/// be set to radian or degrees using the optional argument c, see the documentation
1530/// of the DeltaPhi helper.
1531template <typename T>
1532RVec<T> DeltaR(const RVec<T>& eta1, const RVec<T>& eta2, const RVec<T>& phi1, const RVec<T>& phi2, const T c = M_PI)
1533{
1534 return sqrt(DeltaR2(eta1, eta2, phi1, phi2, c));
1535}
1536
1537/// Return the distance on the \f$\eta\f$-\f$\phi\f$ plane (\f$\Delta R\f$) from
1538/// the scalars eta1, eta2, phi1 and phi2.
1539///
1540/// The function computes \f$\Delta R = \sqrt{(\eta_1 - \eta_2)^2 + (\phi_1 - \phi_2)^2}\f$
1541/// of the given scalars eta1, eta2, phi1 and phi2. The angle \f$\phi\f$ can
1542/// be set to radian or degrees using the optional argument c, see the documentation
1543/// of the DeltaPhi helper.
1544template <typename T>
1545T DeltaR(T eta1, T eta2, T phi1, T phi2, const T c = M_PI)
1546{
1547 const auto dphi = DeltaPhi(phi1, phi2, c);
1548 return std::sqrt((eta1 - eta2) * (eta1 - eta2) + dphi * dphi);
1549}
1550
1551/// Return the invariant mass of two particles given the collections of the quantities
1552/// transverse momentum (pt), rapidity (eta), azimuth (phi) and mass.
1553///
1554/// The function computes the invariant mass of two particles with the four-vectors
1555/// (pt1, eta2, phi1, mass1) and (pt2, eta2, phi2, mass2).
1556template <typename T>
1558 const RVec<T>& pt1, const RVec<T>& eta1, const RVec<T>& phi1, const RVec<T>& mass1,
1559 const RVec<T>& pt2, const RVec<T>& eta2, const RVec<T>& phi2, const RVec<T>& mass2)
1560{
1561 std::size_t size = pt1.size();
1562
1563 R__ASSERT(eta1.size() == size && phi1.size() == size && mass1.size() == size);
1564 R__ASSERT(pt2.size() == size && phi2.size() == size && mass2.size() == size);
1565
1566 RVec<T> inv_masses(size);
1567
1568 for (std::size_t i = 0u; i < size; ++i) {
1569 // Conversion from (pt, eta, phi, mass) to (x, y, z, e) coordinate system
1570 const auto x1 = pt1[i] * std::cos(phi1[i]);
1571 const auto y1 = pt1[i] * std::sin(phi1[i]);
1572 const auto z1 = pt1[i] * std::sinh(eta1[i]);
1573 const auto e1 = std::sqrt(x1 * x1 + y1 * y1 + z1 * z1 + mass1[i] * mass1[i]);
1574
1575 const auto x2 = pt2[i] * std::cos(phi2[i]);
1576 const auto y2 = pt2[i] * std::sin(phi2[i]);
1577 const auto z2 = pt2[i] * std::sinh(eta2[i]);
1578 const auto e2 = std::sqrt(x2 * x2 + y2 * y2 + z2 * z2 + mass2[i] * mass2[i]);
1579
1580 // Addition of particle four-vector elements
1581 const auto e = e1 + e2;
1582 const auto x = x1 + x2;
1583 const auto y = y1 + y2;
1584 const auto z = z1 + z2;
1585
1586 inv_masses[i] = std::sqrt(e * e - x * x - y * y - z * z);
1587 }
1588
1589 // Return invariant mass with (+, -, -, -) metric
1590 return inv_masses;
1591}
1592
1593/// Return the invariant mass of multiple particles given the collections of the
1594/// quantities transverse momentum (pt), rapidity (eta), azimuth (phi) and mass.
1595///
1596/// The function computes the invariant mass of multiple particles with the
1597/// four-vectors (pt, eta, phi, mass).
1598template <typename T>
1599T InvariantMass(const RVec<T>& pt, const RVec<T>& eta, const RVec<T>& phi, const RVec<T>& mass)
1600{
1601 const std::size_t size = pt.size();
1602
1603 R__ASSERT(eta.size() == size && phi.size() == size && mass.size() == size);
1604
1605 T x_sum = 0.;
1606 T y_sum = 0.;
1607 T z_sum = 0.;
1608 T e_sum = 0.;
1609
1610 for (std::size_t i = 0u; i < size; ++ i) {
1611 // Convert to (e, x, y, z) coordinate system and update sums
1612 const auto x = pt[i] * std::cos(phi[i]);
1613 x_sum += x;
1614 const auto y = pt[i] * std::sin(phi[i]);
1615 y_sum += y;
1616 const auto z = pt[i] * std::sinh(eta[i]);
1617 z_sum += z;
1618 const auto e = std::sqrt(x * x + y * y + z * z + mass[i] * mass[i]);
1619 e_sum += e;
1620 }
1621
1622 // Return invariant mass with (+, -, -, -) metric
1623 return std::sqrt(e_sum * e_sum - x_sum * x_sum - y_sum * y_sum - z_sum * z_sum);
1624}
1625
1626////////////////////////////////////////////////////////////////////////////
1627/// \brief Build an RVec of objects starting from RVecs of input to their constructors.
1628/// \tparam T Type of the objects contained in the created RVec.
1629/// \tparam Args_t Pack of types templating the input RVecs.
1630/// \param[in] args The RVecs containing the values used to initialise the output objects.
1631/// \return The RVec of objects initialised with the input parameters.
1632///
1633/// Example code, at the ROOT prompt:
1634/// ~~~{.cpp}
1635/// using namespace ROOT::VecOps;
1636/// RVec<float> pts = {15.5, 34.32, 12.95};
1637/// RVec<float> etas = {0.3, 2.2, 1.32};
1638/// RVec<float> phis = {0.1, 3.02, 2.2};
1639/// RVec<float> masses = {105.65, 105.65, 105.65};
1640/// auto fourVecs = Construct<ROOT::Math::PtEtaPhiMVector>(pts, etas, phis, masses);
1641/// cout << fourVecs << endl;
1642/// // { (15.5,0.3,0.1,105.65), (34.32,2.2,3.02,105.65), (12.95,1.32,2.2,105.65) }
1643/// ~~~
1644template <typename T, typename... Args_t>
1646{
1647 const auto size = ::ROOT::Detail::VecOps::GetVectorsSize("Construct", args...);
1648 RVec<T> ret;
1649 ret.reserve(size);
1650 for (auto i = 0UL; i < size; ++i) {
1651 ret.emplace_back(args[i]...);
1652 }
1653 return ret;
1654}
1655
1656////////////////////////////////////////////////////////////////////////////////
1657/// Print a RVec at the prompt:
1658template <class T>
1659std::ostream &operator<<(std::ostream &os, const RVec<T> &v)
1660{
1661 // In order to print properly, convert to 64 bit int if this is a char
1662 constexpr bool mustConvert = std::is_same<char, T>::value || std::is_same<signed char, T>::value ||
1663 std::is_same<unsigned char, T>::value || std::is_same<wchar_t, T>::value ||
1664 std::is_same<char16_t, T>::value || std::is_same<char32_t, T>::value;
1666 os << "{ ";
1667 auto size = v.size();
1668 if (size) {
1669 for (std::size_t i = 0; i < size - 1; ++i) {
1670 os << (Print_t)v[i] << ", ";
1671 }
1672 os << (Print_t)v[size - 1];
1673 }
1674 os << " }";
1675 return os;
1676}
1677
1678#if (_VECOPS_USE_EXTERN_TEMPLATES)
1679
1680#define RVEC_EXTERN_UNARY_OPERATOR(T, OP) \
1681 extern template RVec<T> operator OP<T>(const RVec<T> &);
1682
1683#define RVEC_EXTERN_BINARY_OPERATOR(T, OP) \
1684 extern template auto operator OP<T, T>(const T &x, const RVec<T> &v) \
1685 -> RVec<decltype(x OP v[0])>; \
1686 extern template auto operator OP<T, T>(const RVec<T> &v, const T &y) \
1687 -> RVec<decltype(v[0] OP y)>; \
1688 extern template auto operator OP<T, T>(const RVec<T> &v0, const RVec<T> &v1)\
1689 -> RVec<decltype(v0[0] OP v1[0])>;
1690
1691#define RVEC_EXTERN_ASSIGN_OPERATOR(T, OP) \
1692 extern template RVec<T> &operator OP<T, T>(RVec<T> &, const T &); \
1693 extern template RVec<T> &operator OP<T, T>(RVec<T> &, const RVec<T> &);
1694
1695#define RVEC_EXTERN_LOGICAL_OPERATOR(T, OP) \
1696 extern template RVec<int> operator OP<T, T>(const RVec<T> &, const T &); \
1697 extern template RVec<int> operator OP<T, T>(const T &, const RVec<T> &); \
1698 extern template RVec<int> operator OP<T, T>(const RVec<T> &, const RVec<T> &);
1699
1700#define RVEC_EXTERN_FLOAT_TEMPLATE(T) \
1701 extern template class RVec<T>; \
1702 RVEC_EXTERN_UNARY_OPERATOR(T, +) \
1703 RVEC_EXTERN_UNARY_OPERATOR(T, -) \
1704 RVEC_EXTERN_UNARY_OPERATOR(T, !) \
1705 RVEC_EXTERN_BINARY_OPERATOR(T, +) \
1706 RVEC_EXTERN_BINARY_OPERATOR(T, -) \
1707 RVEC_EXTERN_BINARY_OPERATOR(T, *) \
1708 RVEC_EXTERN_BINARY_OPERATOR(T, /) \
1709 RVEC_EXTERN_ASSIGN_OPERATOR(T, +=) \
1710 RVEC_EXTERN_ASSIGN_OPERATOR(T, -=) \
1711 RVEC_EXTERN_ASSIGN_OPERATOR(T, *=) \
1712 RVEC_EXTERN_ASSIGN_OPERATOR(T, /=) \
1713 RVEC_EXTERN_LOGICAL_OPERATOR(T, <) \
1714 RVEC_EXTERN_LOGICAL_OPERATOR(T, >) \
1715 RVEC_EXTERN_LOGICAL_OPERATOR(T, ==) \
1716 RVEC_EXTERN_LOGICAL_OPERATOR(T, !=) \
1717 RVEC_EXTERN_LOGICAL_OPERATOR(T, <=) \
1718 RVEC_EXTERN_LOGICAL_OPERATOR(T, >=) \
1719 RVEC_EXTERN_LOGICAL_OPERATOR(T, &&) \
1720 RVEC_EXTERN_LOGICAL_OPERATOR(T, ||)
1721
1722#define RVEC_EXTERN_INTEGER_TEMPLATE(T) \
1723 extern template class RVec<T>; \
1724 RVEC_EXTERN_UNARY_OPERATOR(T, +) \
1725 RVEC_EXTERN_UNARY_OPERATOR(T, -) \
1726 RVEC_EXTERN_UNARY_OPERATOR(T, ~) \
1727 RVEC_EXTERN_UNARY_OPERATOR(T, !) \
1728 RVEC_EXTERN_BINARY_OPERATOR(T, +) \
1729 RVEC_EXTERN_BINARY_OPERATOR(T, -) \
1730 RVEC_EXTERN_BINARY_OPERATOR(T, *) \
1731 RVEC_EXTERN_BINARY_OPERATOR(T, /) \
1732 RVEC_EXTERN_BINARY_OPERATOR(T, %) \
1733 RVEC_EXTERN_BINARY_OPERATOR(T, &) \
1734 RVEC_EXTERN_BINARY_OPERATOR(T, |) \
1735 RVEC_EXTERN_BINARY_OPERATOR(T, ^) \
1736 RVEC_EXTERN_ASSIGN_OPERATOR(T, +=) \
1737 RVEC_EXTERN_ASSIGN_OPERATOR(T, -=) \
1738 RVEC_EXTERN_ASSIGN_OPERATOR(T, *=) \
1739 RVEC_EXTERN_ASSIGN_OPERATOR(T, /=) \
1740 RVEC_EXTERN_ASSIGN_OPERATOR(T, %=) \
1741 RVEC_EXTERN_ASSIGN_OPERATOR(T, &=) \
1742 RVEC_EXTERN_ASSIGN_OPERATOR(T, |=) \
1743 RVEC_EXTERN_ASSIGN_OPERATOR(T, ^=) \
1744 RVEC_EXTERN_ASSIGN_OPERATOR(T, >>=) \
1745 RVEC_EXTERN_ASSIGN_OPERATOR(T, <<=) \
1746 RVEC_EXTERN_LOGICAL_OPERATOR(T, <) \
1747 RVEC_EXTERN_LOGICAL_OPERATOR(T, >) \
1748 RVEC_EXTERN_LOGICAL_OPERATOR(T, ==) \
1749 RVEC_EXTERN_LOGICAL_OPERATOR(T, !=) \
1750 RVEC_EXTERN_LOGICAL_OPERATOR(T, <=) \
1751 RVEC_EXTERN_LOGICAL_OPERATOR(T, >=) \
1752 RVEC_EXTERN_LOGICAL_OPERATOR(T, &&) \
1753 RVEC_EXTERN_LOGICAL_OPERATOR(T, ||)
1754
1755RVEC_EXTERN_INTEGER_TEMPLATE(char)
1756RVEC_EXTERN_INTEGER_TEMPLATE(short)
1757RVEC_EXTERN_INTEGER_TEMPLATE(int)
1758RVEC_EXTERN_INTEGER_TEMPLATE(long)
1759//RVEC_EXTERN_INTEGER_TEMPLATE(long long)
1760
1761RVEC_EXTERN_INTEGER_TEMPLATE(unsigned char)
1762RVEC_EXTERN_INTEGER_TEMPLATE(unsigned short)
1763RVEC_EXTERN_INTEGER_TEMPLATE(unsigned int)
1764RVEC_EXTERN_INTEGER_TEMPLATE(unsigned long)
1765//RVEC_EXTERN_INTEGER_TEMPLATE(unsigned long long)
1766
1767RVEC_EXTERN_FLOAT_TEMPLATE(float)
1768RVEC_EXTERN_FLOAT_TEMPLATE(double)
1769
1770#undef RVEC_EXTERN_UNARY_OPERATOR
1771#undef RVEC_EXTERN_BINARY_OPERATOR
1772#undef RVEC_EXTERN_ASSIGN_OPERATOR
1773#undef RVEC_EXTERN_LOGICAL_OPERATOR
1774#undef RVEC_EXTERN_INTEGER_TEMPLATE
1775#undef RVEC_EXTERN_FLOAT_TEMPLATE
1776
1777#define RVEC_EXTERN_UNARY_FUNCTION(T, NAME, FUNC) \
1778 extern template RVec<PromoteType<T>> NAME(const RVec<T> &);
1779
1780#define RVEC_EXTERN_STD_UNARY_FUNCTION(T, F) RVEC_EXTERN_UNARY_FUNCTION(T, F, std::F)
1781
1782#define RVEC_EXTERN_BINARY_FUNCTION(T0, T1, NAME, FUNC) \
1783 extern template RVec<PromoteTypes<T0, T1>> NAME(const RVec<T0> &, const T1 &); \
1784 extern template RVec<PromoteTypes<T0, T1>> NAME(const T0 &, const RVec<T1> &); \
1785 extern template RVec<PromoteTypes<T0, T1>> NAME(const RVec<T0> &, const RVec<T1> &);
1786
1787#define RVEC_EXTERN_STD_BINARY_FUNCTION(T, F) RVEC_EXTERN_BINARY_FUNCTION(T, T, F, std::F)
1788
1789#define RVEC_EXTERN_STD_FUNCTIONS(T) \
1790 RVEC_EXTERN_STD_UNARY_FUNCTION(T, abs) \
1791 RVEC_EXTERN_STD_BINARY_FUNCTION(T, fdim) \
1792 RVEC_EXTERN_STD_BINARY_FUNCTION(T, fmod) \
1793 RVEC_EXTERN_STD_BINARY_FUNCTION(T, remainder) \
1794 RVEC_EXTERN_STD_UNARY_FUNCTION(T, exp) \
1795 RVEC_EXTERN_STD_UNARY_FUNCTION(T, exp2) \
1796 RVEC_EXTERN_STD_UNARY_FUNCTION(T, expm1) \
1797 RVEC_EXTERN_STD_UNARY_FUNCTION(T, log) \
1798 RVEC_EXTERN_STD_UNARY_FUNCTION(T, log10) \
1799 RVEC_EXTERN_STD_UNARY_FUNCTION(T, log2) \
1800 RVEC_EXTERN_STD_UNARY_FUNCTION(T, log1p) \
1801 RVEC_EXTERN_STD_BINARY_FUNCTION(T, pow) \
1802 RVEC_EXTERN_STD_UNARY_FUNCTION(T, sqrt) \
1803 RVEC_EXTERN_STD_UNARY_FUNCTION(T, cbrt) \
1804 RVEC_EXTERN_STD_BINARY_FUNCTION(T, hypot) \
1805 RVEC_EXTERN_STD_UNARY_FUNCTION(T, sin) \
1806 RVEC_EXTERN_STD_UNARY_FUNCTION(T, cos) \
1807 RVEC_EXTERN_STD_UNARY_FUNCTION(T, tan) \
1808 RVEC_EXTERN_STD_UNARY_FUNCTION(T, asin) \
1809 RVEC_EXTERN_STD_UNARY_FUNCTION(T, acos) \
1810 RVEC_EXTERN_STD_UNARY_FUNCTION(T, atan) \
1811 RVEC_EXTERN_STD_BINARY_FUNCTION(T, atan2) \
1812 RVEC_EXTERN_STD_UNARY_FUNCTION(T, sinh) \
1813 RVEC_EXTERN_STD_UNARY_FUNCTION(T, cosh) \
1814 RVEC_EXTERN_STD_UNARY_FUNCTION(T, tanh) \
1815 RVEC_EXTERN_STD_UNARY_FUNCTION(T, asinh) \
1816 RVEC_EXTERN_STD_UNARY_FUNCTION(T, acosh) \
1817 RVEC_EXTERN_STD_UNARY_FUNCTION(T, atanh) \
1818 RVEC_EXTERN_STD_UNARY_FUNCTION(T, floor) \
1819 RVEC_EXTERN_STD_UNARY_FUNCTION(T, ceil) \
1820 RVEC_EXTERN_STD_UNARY_FUNCTION(T, trunc) \
1821 RVEC_EXTERN_STD_UNARY_FUNCTION(T, round) \
1822 RVEC_EXTERN_STD_UNARY_FUNCTION(T, erf) \
1823 RVEC_EXTERN_STD_UNARY_FUNCTION(T, erfc) \
1824 RVEC_EXTERN_STD_UNARY_FUNCTION(T, lgamma) \
1825 RVEC_EXTERN_STD_UNARY_FUNCTION(T, tgamma) \
1826
1827RVEC_EXTERN_STD_FUNCTIONS(float)
1828RVEC_EXTERN_STD_FUNCTIONS(double)
1829#undef RVEC_EXTERN_STD_UNARY_FUNCTION
1830#undef RVEC_EXTERN_STD_BINARY_FUNCTION
1831#undef RVEC_EXTERN_STD_UNARY_FUNCTIONS
1832
1833#ifdef R__HAS_VDT
1834
1835#define RVEC_EXTERN_VDT_UNARY_FUNCTION(T, F) RVEC_EXTERN_UNARY_FUNCTION(T, F, vdt::F)
1836
1837RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_expf)
1838RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_logf)
1839RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_sinf)
1840RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_cosf)
1841RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_tanf)
1842RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_asinf)
1843RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_acosf)
1844RVEC_EXTERN_VDT_UNARY_FUNCTION(float, fast_atanf)
1845
1846RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_exp)
1847RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_log)
1848RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_sin)
1849RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_cos)
1850RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_tan)
1851RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_asin)
1852RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_acos)
1853RVEC_EXTERN_VDT_UNARY_FUNCTION(double, fast_atan)
1854
1855#endif // R__HAS_VDT
1856
1857#endif // _VECOPS_USE_EXTERN_TEMPLATES
1858
1859} // End of VecOps NS
1860
1861// Allow to use RVec as ROOT::RVec
1862using ROOT::VecOps::RVec;
1863
1864} // End of ROOT NS
1865
1866#endif // ROOT_RVEC
double
Definition: Converters.cxx:921
long
Definition: Converters.cxx:858
ROOT::R::TRInterface & r
Definition: Object.C:4
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
#define e(i)
Definition: RSha256.hxx:103
#define RVEC_UNARY_OPERATOR(OP)
Definition: RVec.hxx:434
#define RVEC_ASSIGNMENT_OPERATOR(OP)
Definition: RVec.hxx:505
#define RVEC_STD_BINARY_FUNCTION(F)
Definition: RVec.hxx:648
#define RVEC_BINARY_OPERATOR(OP)
Definition: RVec.hxx:457
#define RVEC_LOGICAL_OPERATOR(OP)
Definition: RVec.hxx:541
#define RVEC_STD_UNARY_FUNCTION(F)
Definition: RVec.hxx:647
static const double x2[5]
static const double x1[5]
#define M_PI
Definition: Rotated.cxx:105
#define R__ASSERT(e)
Definition: TError.h:96
Int_t Compare(const void *item1, const void *item2)
int type
Definition: TGX11.cxx:120
double atan2(double, double)
double tanh(double)
double cosh(double)
double ceil(double)
double acos(double)
double sinh(double)
double cos(double)
double pow(double, double)
double log10(double)
double floor(double)
double atan(double)
double tan(double)
double sqrt(double)
double sin(double)
double asin(double)
double exp(double)
double log(double)
A "std::vector"-like collection of values implementing handy operation to analyse them.
Definition: RVec.hxx:275
typename std::conditional< IsVecBool, std::vector< bool >, std::vector< T, ::ROOT::Detail::VecOps::RAdoptAllocator< T > > >::type Impl_t
Definition: RVec.hxx:282
typename Impl_t::iterator iterator
Definition: RVec.hxx:294
static constexpr const auto IsVecBool
Definition: RVec.hxx:280
const_iterator begin() const noexcept
Definition: RVec.hxx:388
RVec(const RVec< T > &v)
Definition: RVec.hxx:310
reverse_iterator rbegin() noexcept
Definition: RVec.hxx:393
const_iterator cbegin() const noexcept
Definition: RVec.hxx:389
typename Impl_t::value_type value_type
Definition: RVec.hxx:283
reference operator[](size_type pos)
Definition: RVec.hxx:361
void clear() noexcept
Definition: RVec.hxx:407
iterator erase(iterator pos)
Definition: RVec.hxx:408
RVec< T > & operator=(std::initializer_list< T > ilist)
Definition: RVec.hxx:336
typename Impl_t::const_pointer const_pointer
Definition: RVec.hxx:289
RVec(std::initializer_list< T > init)
Definition: RVec.hxx:321
void resize(size_type count)
Definition: RVec.hxx:426
void push_back(T &&value)
Definition: RVec.hxx:410
const_reference back() const
Definition: RVec.hxx:383
typename Impl_t::difference_type difference_type
Definition: RVec.hxx:285
const_reverse_iterator crbegin() const noexcept
Definition: RVec.hxx:395
RVec operator[](const RVec< V > &conds) const
Definition: RVec.hxx:365
iterator erase(iterator first, iterator last)
Definition: RVec.hxx:409
RVec< T > & operator=(RVec< T > &&v)
Definition: RVec.hxx:330
const_reverse_iterator crend() const noexcept
Definition: RVec.hxx:398
const_reference at(size_type pos) const
Definition: RVec.hxx:356
typename Impl_t::const_iterator const_iterator
Definition: RVec.hxx:295
bool empty() const noexcept
Definition: RVec.hxx:400
reference back()
Definition: RVec.hxx:382
typename Impl_t::reference reference
Definition: RVec.hxx:286
void push_back(const value_type &value)
Definition: RVec.hxx:411
iterator end() noexcept
Definition: RVec.hxx:390
size_type size() const noexcept
Definition: RVec.hxx:401
RVec(size_type count, const T &value)
Definition: RVec.hxx:308
const_iterator end() const noexcept
Definition: RVec.hxx:391
iterator emplace(const_iterator pos, U value)
This method is intended only for arithmetic types unlike the std::vector corresponding one which is g...
Definition: RVec.hxx:421
typename Impl_t::const_reference const_reference
Definition: RVec.hxx:287
const_iterator cend() const noexcept
Definition: RVec.hxx:392
RVec(RVec< T > &&v)
Definition: RVec.hxx:312
RVec(size_type count)
Definition: RVec.hxx:306
typename Impl_t::const_reverse_iterator const_reverse_iterator
Definition: RVec.hxx:297
const_reference front() const
Definition: RVec.hxx:381
void reserve(size_type new_cap)
Definition: RVec.hxx:403
reference at(size_type pos)
Definition: RVec.hxx:355
const_reverse_iterator rend() const noexcept
Definition: RVec.hxx:397
typename std::conditional< IsVecBool, void, typename Impl_t::pointer >::type data_t
Definition: RVec.hxx:292
reference front()
Definition: RVec.hxx:380
typename Impl_t::reverse_iterator reverse_iterator
Definition: RVec.hxx:296
const_reference operator[](size_type pos) const
Definition: RVec.hxx:362
typename Impl_t::pointer pointer
Definition: RVec.hxx:288
size_type capacity() const noexcept
Definition: RVec.hxx:404
reverse_iterator rend() noexcept
Definition: RVec.hxx:396
RVec< T > & operator=(const RVec< T > &v)
Definition: RVec.hxx:324
typename std::conditional< IsVecBool, void, typename Impl_t::const_pointer >::type const_data_t
Definition: RVec.hxx:293
void shrink_to_fit()
Definition: RVec.hxx:405
const Impl_t & AsVector() const
Definition: RVec.hxx:351
value_type at(size_type pos, value_type fallback)
No exception thrown. The user specifies the desired value in case the RVec is shorter than pos.
Definition: RVec.hxx:358
iterator begin() noexcept
Definition: RVec.hxx:387
Impl_t & AsVector()
Definition: RVec.hxx:352
const_reverse_iterator rbegin() const noexcept
Definition: RVec.hxx:394
value_type at(size_type pos, value_type fallback) const
No exception thrown. The user specifies the desired value in case the RVec is shorter than pos.
Definition: RVec.hxx:360
typename Impl_t::size_type size_type
Definition: RVec.hxx:284
void resize(size_type count, const value_type &value)
Definition: RVec.hxx:427
reference emplace_back(Args &&... args)
Definition: RVec.hxx:413
const_data_t data() const noexcept
Definition: RVec.hxx:385
RVec(const std::vector< T > &v)
Definition: RVec.hxx:314
data_t data() noexcept
Definition: RVec.hxx:384
void swap(RVec< T > &other)
Definition: RVec.hxx:428
RVec(InputIt first, InputIt last)
Definition: RVec.hxx:319
RVec(pointer p, size_type n)
Definition: RVec.hxx:316
void pop_back()
Definition: RVec.hxx:425
size_type max_size() const noexcept
Definition: RVec.hxx:402
TPaveText * pt
Type
enumeration specifying the integration types.
double erfc(double x)
Complementary error function.
double tgamma(double x)
The gamma function is defined to be the extension of the factorial to real numbers.
double lgamma(double x)
Calculates the logarithm of the gamma function.
double erf(double x)
Error function encountered in integrating the normal distribution.
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
basic_string_view< char > string_view
#define F(x, y, z)
EvaluateInfo init(std::vector< RooRealProxy > parameters, std::vector< ArrayWrapper * > wrappers, std::vector< double * > arrays, size_t begin, size_t batchSize)
std::size_t GetVectorsSize(std::string_view id, const RVec< T > &... vs)
Definition: RVec.hxx:66
auto MapFromTuple(Tuple_t &&t, std::index_sequence< Is... >) -> decltype(MapImpl(std::get< std::tuple_size< Tuple_t >::value - 1 >(t), std::get< Is >(t)...))
Definition: RVec.hxx:95
ROOT::VecOps::RVec< T > RVec
Definition: RVec.hxx:63
auto MapImpl(F &&f, const RVec< T > &... vs) -> RVec< decltype(f(vs[0]...))>
Definition: RVec.hxx:83
void EmplaceBack(T &v, Args &&... args)
Definition: RVec.hxx:112
double T(double x)
Definition: ChebyshevPol.h:34
double log1p(double x)
declarations for functions which are not implemented by some compilers
Definition: Math.h:98
double expm1(double x)
exp(x) -1 with error cancellation when x is small
Definition: Math.h:110
double inner_product(const LAVector &, const LAVector &)
T Max(const RVec< T > &v)
Get the greatest element of an RVec.
Definition: RVec.hxx:793
RVec< T > Reverse(const RVec< T > &v)
Return copy of reversed vector.
Definition: RVec.hxx:1084
RVec< T > Intersect(const RVec< T > &v1, const RVec< T > &v2, bool v2_is_sorted=false)
Return the intersection of elements of two RVecs.
Definition: RVec.hxx:1284
auto All(const RVec< T > &v) -> decltype(v[0]==false)
Return true if all of the elements equate to true, return false otherwise.
Definition: RVec.hxx:981
T Sum(const RVec< T > &v)
Sum elements of an RVec.
Definition: RVec.hxx:759
RVec< typename RVec< T >::size_type > Nonzero(const RVec< T > &v)
Return the indices of the elements which are not zero.
Definition: RVec.hxx:1253
auto Map(Args &&... args) -> decltype(ROOT::Detail::VecOps::MapFromTuple(std::forward_as_tuple(args...), std::make_index_sequence< sizeof...(args) - 1 >()))
Create new collection applying a callable to the elements of the input collection.
Definition: RVec.hxx:910
RVec< T > Sort(const RVec< T > &v)
Return copy of RVec with elements sorted in ascending order.
Definition: RVec.hxx:1105
T InvariantMass(const RVec< T > &pt, const RVec< T > &eta, const RVec< T > &phi, const RVec< T > &mass)
Return the invariant mass of multiple particles given the collections of the quantities transverse mo...
Definition: RVec.hxx:1599
std::ostream & operator<<(std::ostream &os, const RVec< T > &v)
Print a RVec at the prompt:
Definition: RVec.hxx:1659
RVec< Common_t > Concatenate(const RVec< T0 > &v0, const RVec< T1 > &v1)
Return the concatenation of two RVecs.
Definition: RVec.hxx:1423
RVec< T > DeltaR(const RVec< T > &eta1, const RVec< T > &eta2, const RVec< T > &phi1, const RVec< T > &phi2, const T c=M_PI)
Return the distance on the - plane ( ) from the collections eta1, eta2, phi1 and phi2.
Definition: RVec.hxx:1532
RVec< T > InvariantMasses(const RVec< T > &pt1, const RVec< T > &eta1, const RVec< T > &phi1, const RVec< T > &mass1, const RVec< T > &pt2, const RVec< T > &eta2, const RVec< T > &phi2, const RVec< T > &mass2)
Return the invariant mass of two particles given the collections of the quantities transverse momentu...
Definition: RVec.hxx:1557
RVec< T > DeltaR2(const RVec< T > &eta1, const RVec< T > &eta2, const RVec< T > &phi1, const RVec< T > &phi2, const T c=M_PI)
Return the square of the distance on the - plane ( ) from the collections eta1, eta2,...
Definition: RVec.hxx:1518
double Mean(const RVec< T > &v)
Get the mean of the elements of an RVec.
Definition: RVec.hxx:776
RVec< T > Take(const RVec< T > &v, const RVec< typename RVec< T >::size_type > &i)
Return elements of a vector at given indices.
Definition: RVec.hxx:1026
void swap(RVec< T > &lhs, RVec< T > &rhs)
Definition: RVec.hxx:990
RVec< T > Construct(const RVec< Args_t > &... args)
Build an RVec of objects starting from RVecs of input to their constructors.
Definition: RVec.hxx:1645
T DeltaPhi(T v1, T v2, const T c=M_PI)
Return the angle difference of two scalars.
Definition: RVec.hxx:1442
auto Dot(const RVec< T > &v0, const RVec< V > &v1) -> decltype(v0[0] *v1[0])
Inner product.
Definition: RVec.hxx:741
T Min(const RVec< T > &v)
Get the smallest element of an RVec.
Definition: RVec.hxx:809
RVec< RVec< std::size_t > > Combinations(const std::size_t size1, const std::size_t size2)
Return the indices that represent all combinations of the elements of two RVecs.
Definition: RVec.hxx:1149
RVec< T > Where(const RVec< int > &c, const RVec< T > &v1, const RVec< T > &v2)
Return the elements of v1 if the condition c is true and v2 if the condition c is false.
Definition: RVec.hxx:1318
double StdDev(const RVec< T > &v)
Get the standard deviation of the elements of an RVec.
Definition: RVec.hxx:886
auto Any(const RVec< T > &v) -> decltype(v[0]==true)
Return true if any of the elements equates to true, return false otherwise.
Definition: RVec.hxx:962
RVec< typename RVec< T >::size_type > Argsort(const RVec< T > &v)
Return an RVec of indices that sort the input RVec.
Definition: RVec.hxx:1006
std::size_t ArgMin(const RVec< T > &v)
Get the index of the smallest element of an RVec In case of multiple occurrences of the minimum value...
Definition: RVec.hxx:845
double Var(const RVec< T > &v)
Get the variance of the elements of an RVec.
Definition: RVec.hxx:862
RVec< T > Filter(const RVec< T > &v, F &&f)
Create a new collection with the elements passing the filter expressed by the predicate.
Definition: RVec.hxx:939
std::size_t ArgMax(const RVec< T > &v)
Get the index of the greatest element of an RVec In case of multiple occurrences of the maximum value...
Definition: RVec.hxx:827
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
Definition: StringConv.hxx:21
static constexpr double s
Definition: first.py:1