Tutorial june2020 part2

DataFormats/Math/src/classes_def.xml

@@ -4,12 +4,15 @@
   <class name="edm::Wrapper<edm::ValueMap<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > >"/>
   <class name="edm::Wrapper<ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<float>,ROOT::Math::DefaultCoordinateSystemTag> >"/>
   <class name="edm::Wrapper<std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > >"/>
+  <class name="edm::Wrapper<std::vector<ROOT::Math::DisplacementVector3D<ROOT::Math::CylindricalEta3D<float>,ROOT::Math::DefaultCoordinateSystemTag> > >"/>
+  <class name="edm::Wrapper<std::vector<ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<float>,ROOT::Math::DefaultCoordinateSystemTag> > >"/>
   <class name="std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > >"/>
 <selection>
   <class pattern="ROOT::Math::SMatrix<*>" />
   <class pattern="ROOT::Math::MatRepStd<*>" />
   <class pattern="ROOT::Math::RowOffsets<*>" />
   <class pattern="ROOT::Math::SVector<*>" />
+  <class pattern="std::vector<ROOT::Math::DisplacementVector3D<*>>" />
   <class pattern="std::vector<ROOT::Math::SMatrix<*>>" />
   <class pattern="std::vector<std::pair<ROOT::Math::PositionVector3D<*>,float> >"/>
 </selection>

Patatrack/Tutorial/plugins/BuildFile.xml

@@ -0,0 +1,8 @@
+<use name="cuda"/>
+<use name="DataFormats/Math"/>
+<use name="FWCore/Framework"/>
+<use name="FWCore/ParameterSet"/>
+<use name="FWCore/PluginManager"/>
+<use name="HeterogeneousCore/CUDAUtilities"/>
+<use name="HeterogeneousCore/CUDACore"/>
+<flags EDM_PLUGIN="1"/>

Patatrack/Tutorial/plugins/CompareCartesianVectors.cc

@@ -0,0 +1,84 @@
+// system include files
+#include <cmath>
+#include <iomanip>
+#include <iostream>
+#include <memory>
+#include <vector>
+
+// CMSSW include files
+#include "DataFormats/Math/interface/Vector3D.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/one/EDAnalyzer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+
+class CompareCartesianVectors : public edm::one::EDAnalyzer<> {
+public:
+  explicit CompareCartesianVectors(const edm::ParameterSet&);
+  ~CompareCartesianVectors() = default;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  using CartesianVectors = std::vector<math::XYZVectorF>;
+
+  virtual void analyze(const edm::Event&, const edm::EventSetup&) override;
+
+  bool compare(math::XYZVectorF const& first, math::XYZVectorF const& second) const;
+
+  edm::EDGetTokenT<CartesianVectors> first_;
+  edm::EDGetTokenT<CartesianVectors> second_;
+  const double precision_;
+};
+
+CompareCartesianVectors::CompareCartesianVectors(const edm::ParameterSet& config)
+    : first_(consumes<CartesianVectors>(config.getParameter<edm::InputTag>("first"))),
+      second_(consumes<CartesianVectors>(config.getParameter<edm::InputTag>("second"))),
+      precision_(config.getParameter<double>("precision")) {}
+
+void CompareCartesianVectors::analyze(const edm::Event& event, const edm::EventSetup& setup) {
+  auto const& first = event.get(first_);
+  auto const& second = event.get(second_);
+  if (first.size() != second.size()) {
+    std::cout << "The two collections have different size: " << first.size() << " and " << second.size() << "."
+              << std::endl;
+    return;
+  }
+  bool consistent = true;
+  for (unsigned int i = 0; i < first.size(); ++i) {
+    if (not compare(first[i], second[i])) {
+      if (consistent) {
+        std::cout << "Found inconsistent elements:" << std::endl;
+        consistent = false;
+      }
+      std::cout << std::setprecision(9) << "(" << first[i].x() << ", " << first[i].y() << ", " << first[i].z()
+                << ")  vs  (" << second[i].x() << ", " << second[i].y() << ", " << second[i].z() << ")" << std::endl;
+    }
+  }
+  if (consistent) {
+    std::cout << "All elements are consistent within " << precision_ << std::endl;
+  }
+}
+
+bool CompareCartesianVectors::compare(math::XYZVectorF const& first, math::XYZVectorF const& second) const {
+  if (std::abs(first.x() - second.x()) > std::abs(first.x() + second.x()) * precision_)
+    return false;
+  if (std::abs(first.y() - second.y()) > std::abs(first.y() + second.y()) * precision_)
+    return false;
+  if (std::abs(first.z() - second.z()) > std::abs(first.z() + second.z()) * precision_)
+    return false;
+  return true;
+}
+
+void CompareCartesianVectors::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag>("first", edm::InputTag("firstCartesianVectors"));
+  desc.add<edm::InputTag>("second", edm::InputTag("secondCartesianVectors"));
+  desc.add<double>("precision", 1.e-6);
+  descriptions.addWithDefaultLabel(desc);
+}
+
+// define this as a plug-in
+DEFINE_FWK_MODULE(CompareCartesianVectors);

Patatrack/Tutorial/plugins/ConvertToCartesianVectors.cc

@@ -0,0 +1,65 @@
+// system include files
+#include <cmath>
+#include <memory>
+#include <vector>
+
+// CMSSW include files
+#include "DataFormats/Math/interface/Vector3D.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/StreamID.h"
+
+class ConvertToCartesianVectors : public edm::stream::EDProducer<> {
+public:
+  explicit ConvertToCartesianVectors(const edm::ParameterSet&);
+  ~ConvertToCartesianVectors() = default;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  using CartesianVectors = std::vector<math::XYZVectorF>;
+  using CylindricalVectors = std::vector<math::RhoEtaPhiVectorF>;
+
+  static void convert(math::RhoEtaPhiVectorF const& cilindrical, math::XYZVectorF & cartesian);
+
+  virtual void produce(edm::Event&, const edm::EventSetup&) override;
+
+  edm::EDGetTokenT<CylindricalVectors> input_;
+  edm::EDPutTokenT<CartesianVectors> output_;
+};
+
+ConvertToCartesianVectors::ConvertToCartesianVectors(const edm::ParameterSet& config)
+    : input_(consumes<CylindricalVectors>(config.getParameter<edm::InputTag>("input"))) {
+  output_ = produces<CartesianVectors>();
+}
+
+void ConvertToCartesianVectors::convert(math::RhoEtaPhiVectorF const& cilindrical, math::XYZVectorF & cartesian) {
+    cartesian.SetCoordinates(cilindrical.rho() * std::cos(cilindrical.phi()),
+                             cilindrical.rho() * std::sin(cilindrical.phi()),
+                             cilindrical.rho() * std::sinh(cilindrical.eta()));
+}
+
+void ConvertToCartesianVectors::produce(edm::Event& event, const edm::EventSetup& setup) {
+  auto const& input = event.get(input_);
+  auto elements = input.size();
+  auto product = std::make_unique<CartesianVectors>(elements);
+
+  // convert the vectors from cylindrical to cartesian coordinates
+  for (unsigned int i = 0; i < elements; ++i) {
+    convert(input[i], (*product)[i]);
+  }
+
+  event.put(output_, std::move(product));
+}
+
+void ConvertToCartesianVectors::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag>("input", edm::InputTag("cylindricalVectors"));
+  descriptions.addWithDefaultLabel(desc);
+}
+
+// define this as a plug-in
+DEFINE_FWK_MODULE(ConvertToCartesianVectors);

Patatrack/Tutorial/plugins/ConvertToCartesianVectorsCUDA.cc

@@ -0,0 +1,96 @@
+// system include files
+#include <cmath>
+#include <memory>
+#include <vector>
+
+// CMSSW include files
+#include "DataFormats/Math/interface/Vector3D.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/StreamID.h"
+#include "HeterogeneousCore/CUDAUtilities/interface/cudaCheck.h"
+
+#include "HeterogeneousCore/CUDAUtilities/interface/host_unique_ptr.h"
+#include "HeterogeneousCore/CUDAUtilities/interface/device_unique_ptr.h"
+#include "HeterogeneousCore/CUDAUtilities/interface/host_noncached_unique_ptr.h"
+#include "HeterogeneousCore/CUDACore/interface/ScopedContext.h"
+#include "cudavectors.h"
+
+class ConvertToCartesianVectorsCUDA : public edm::stream::EDProducer<edm::ExternalWork> {
+public:
+  explicit ConvertToCartesianVectorsCUDA(const edm::ParameterSet&);
+  ~ConvertToCartesianVectorsCUDA() = default;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  using CartesianVectors = std::vector<math::XYZVectorF>;
+  using CylindricalVectors = std::vector<math::RhoEtaPhiVectorF>;
+
+  void acquire(edm::Event const& event, edm::EventSetup const& setup, edm::WaitingTaskWithArenaHolder waitingTaskHolder) override;
+  virtual void produce(edm::Event&, const edm::EventSetup&) override;
+  cms::cuda::host::unique_ptr<cudavectors::CartesianVector[]> output_buffer_;
+  edm::EDGetTokenT<CylindricalVectors> input_;
+  edm::EDPutTokenT<CartesianVectors> output_;
+};
+
+ConvertToCartesianVectorsCUDA::ConvertToCartesianVectorsCUDA(const edm::ParameterSet& config)
+    : input_(consumes<CylindricalVectors>(config.getParameter<edm::InputTag>("input"))) {
+  output_ = produces<CartesianVectors>();
+}
+
+
+void ConvertToCartesianVectorsCUDA::acquire(const edm::Event& event, const edm::EventSetup& setup,
+                                                           edm::WaitingTaskWithArenaHolder waitingTaskHolder) {
+
+  cms::cuda::ScopedContextAcquire ctx{event.streamID(), std::move(waitingTaskHolder)};
+
+  auto const& input = event.get(input_);
+  auto elements = input.size();
+
+  // allocate memory on the GPU for the cylindrical and cartesian vectors
+
+  auto gpu_input = cms::cuda::make_device_unique<cudavectors::CylindricalVector[]>(elements, ctx.stream());
+  auto gpu_product = cms::cuda::make_device_unique<cudavectors::CartesianVector[]>(elements, ctx.stream());
+  auto cpu_input = cms::cuda::make_host_noncached_unique<cudavectors::CylindricalVector[]>(elements, cudaHostAllocWriteCombined);
+  output_buffer_ = cms::cuda::make_host_unique<cudavectors::CartesianVector[]>(elements, ctx.stream());
+
+  // copy the input data to the GPU
+
+  std::memcpy(cpu_input.get(), input.data(), sizeof(cudavectors::CylindricalVector) * elements);
+  cudaCheck(cudaMemcpyAsync(gpu_input.get(), cpu_input.get(), sizeof(cudavectors::CylindricalVector) * elements, cudaMemcpyHostToDevice, ctx.stream()));
+
+   // convert the vectors from cylindrical to cartesian coordinates, on the GPU
+
+   cudavectors::convertWrapper(gpu_input.get(), gpu_product.get(), elements, ctx.stream());
+
+  // copy the result from the GPU
+
+ cudaCheck(cudaMemcpyAsync(output_buffer_.get(), gpu_product.get(), sizeof(cudavectors::CartesianVector) * elements, cudaMemcpyDeviceToHost, ctx.stream()));
+
+ // free the GPU memory
+   // no need of explicit free operation 
+} 
+void ConvertToCartesianVectorsCUDA::produce(edm::Event& event, const edm::EventSetup& setup) {
+  //no need for a CUDA context here, because there are no CUDA operations
+
+   auto const& input = event.get(input_);
+   auto elements = input.size();
+
+  // instantiate the event product, copy the results from the output buffer, and free it
+   auto product = std::make_unique<CartesianVectors>(elements);
+   std::memcpy((void*) product->data(), output_buffer_.get(), sizeof(cudavectors::CartesianVector) * elements);
+   output_buffer_.reset();
+   event.put(output_, std::move(product));
+}
+
+void ConvertToCartesianVectorsCUDA::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag>("input", edm::InputTag("cylindricalVectors"));
+  descriptions.addWithDefaultLabel(desc);
+}
+// define this as a plug-in
+DEFINE_FWK_MODULE(ConvertToCartesianVectorsCUDA);

Patatrack/Tutorial/plugins/GenerateCylindricalVectors.cc

@@ -0,0 +1,63 @@
+// system include files
+#include <memory>
+#include <random>
+#include <vector>
+
+// CMSSW include files
+#include "DataFormats/Math/interface/Vector3D.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/StreamID.h"
+
+class GenerateCylindricalVectors : public edm::stream::EDProducer<> {
+public:
+  explicit GenerateCylindricalVectors(const edm::ParameterSet&);
+  ~GenerateCylindricalVectors() = default;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  using CylindricalVectors = std::vector<math::RhoEtaPhiVectorF>;
+
+  void beginStream(edm::StreamID) override {}
+  void produce(edm::Event&, const edm::EventSetup&) override;
+  void endStream() override {}
+
+  std::mt19937 generator_;
+  std::lognormal_distribution<float> genPt_;
+  std::uniform_real_distribution<float> genEta_;
+  std::uniform_real_distribution<float> genPhi_;
+  const uint32_t size_;
+
+  edm::EDPutTokenT<CylindricalVectors> output_;
+};
+
+GenerateCylindricalVectors::GenerateCylindricalVectors(const edm::ParameterSet& config)
+    : generator_(std::random_device()()),
+      genPt_(3, 0.6),
+      genEta_(-5., +5.),
+      genPhi_(0., 2 * M_PI),
+      size_(config.getParameter<uint32_t>("size"))  // number of CylindricalVectors to generate
+{
+  output_ = produces<CylindricalVectors>();
+}
+
+void GenerateCylindricalVectors::produce(edm::Event& event, const edm::EventSetup& setup) {
+  auto product = std::make_unique<CylindricalVectors>(size_);
+  for (auto& p : *product) {
+    p.SetCoordinates(genPt_(generator_), genEta_(generator_), genPhi_(generator_));
+  }
+  event.put(output_, std::move(product));
+}
+
+void GenerateCylindricalVectors::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<uint32_t>("size", 1000)->setComment("number of generated elements");
+  descriptions.addWithDefaultLabel(desc);
+}
+
+// define this as a plug-in
+DEFINE_FWK_MODULE(GenerateCylindricalVectors);

Patatrack/Tutorial/plugins/PrintCartesianVectors.cc

@@ -0,0 +1,47 @@
+// system include files
+#include <memory>
+#include <vector>
+
+// CMSSW include files
+#include "DataFormats/Math/interface/Vector3D.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/one/EDAnalyzer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+
+class PrintCartesianVectors : public edm::one::EDAnalyzer<> {
+public:
+  explicit PrintCartesianVectors(const edm::ParameterSet&);
+  ~PrintCartesianVectors() = default;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  using CartesianVectors = std::vector<math::XYZVectorF>;
+
+  virtual void analyze(const edm::Event&, const edm::EventSetup&) override;
+
+  edm::EDGetTokenT<CartesianVectors> input_;
+};
+
+PrintCartesianVectors::PrintCartesianVectors(const edm::ParameterSet& config)
+    : input_(consumes<CartesianVectors>(config.getParameter<edm::InputTag>("input"))) {}
+
+void PrintCartesianVectors::analyze(const edm::Event& event, const edm::EventSetup& setup) {
+  for (auto const& v : event.get(input_)) {
+    std::cout << std::fixed << "x: " << std::setw(6) << std::setprecision(2) << v.x() << ", y: " << std::setw(6)
+              << std::setprecision(2) << v.y() << ", z: " << std::setw(8) << std::setprecision(2) << v.z() << std::endl;
+  }
+  std::cout << std::endl;
+}
+
+void PrintCartesianVectors::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag>("input", edm::InputTag("cartesianVectors"));
+  descriptions.addWithDefaultLabel(desc);
+}
+
+// define this as a plug-in
+DEFINE_FWK_MODULE(PrintCartesianVectors);