Working afterburner tests

cpp/afterburner/Afterburner.cc

@@ -112,6 +112,13 @@ CLHEP::Hep3Vector ab::Afterburner::smear_beam_divergence(const CLHEP::Hep3Vector
 
     // Resulting
     return y_smear_out_accelerator_plane * x_smear_in_accelerator_plane * beam_dir;
+
+
+}
+
+
+ab::AfterburnerEventResult ab::Afterburner::process_event() {
+    return process_event(CLHEP::HepLorentzVector(0,0,0,0));
 }
 
 //! move vertex in translation,boost,rotation according to vertex settings
@@ -229,7 +236,7 @@ ab::AfterburnerEventResult ab::Afterburner::process_event(const CLHEP::HepLorent
     } else {
         // need a rotation
         CLHEP::Hep3Vector rotation_axis = (beamA_vec - beamB_vec).cross(z_axis);
-        const double rotation_angle_to_z = -acos(cos_rotation_angle_to_z);
+        const double rotation_angle_to_z = acos(cos_rotation_angle_to_z);   // TODO back to -acos? 
 
         result.rotation = CLHEP::HepRotation(rotation_axis, rotation_angle_to_z);
 
@@ -334,3 +341,4 @@ void ab::Afterburner::print() const {
          << ", " << _cfg.beam_two.z_shift_ver << endl;
     cout << "=========================\n";
 }
+

cpp/afterburner/Afterburner.hh

@@ -30,6 +30,8 @@ namespace ab{
     public:
         Afterburner();
 
+        AfterburnerEventResult process_event();
+
         AfterburnerEventResult process_event(const CLHEP::HepLorentzVector &init_vtx);
 
         void print() const;

cpp/afterburner/CMakeLists.txt

@@ -2,13 +2,18 @@ cmake_minimum_required( VERSION 3.5 )
 set(CMAKE_CXX_STANDARD 17)
 project(EICAfterburnerLibrary LANGUAGES CXX)
 
-set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/../cmake/Modules ${CMAKE_MODULE_PATH} )
+
 
 #include(../cmake/FindCLHEP.cmake)
 
 message(STATUS "EIC Afterburner library" )
 
+# CLHEP Package might be built by user and has clhepConfig.cmake or might be built by system
+find_package(CLHEP)
+if(NOT CLHEP_FOUND)
+    set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/../cmake/Modules ${CMAKE_MODULE_PATH} )
     find_package(CLHEP REQUIRED)
+endif()
 find_package(GSL REQUIRED)
 
 add_library(afterburner

cpp/benchmark/ArgumentProcessor.cc

@@ -9,30 +9,39 @@
 UserArguments ArgumentProcessor::Process(int argc, char **argv)
 {
     UserArguments result;                                          // This function result
-    CLI::App app{mAppDescription};
+    CLI::App app{get_description()};
 
-    std::string optOutputName("g4e_output");
+    std::string optOutputName("ab_output.hepmc");
     std::vector<std::string> optAllFiles;
+    std::string benchmarkName("default");
+
     app.add_option("-o,--output", optOutputName, "Base name for Output files");
-    app.add_option("files", optAllFiles, "Input files. ");
+    app.add_option("-b,--benchmark", benchmarkName, "Benchmark name: default, crossing");
+    app.add_option("input file", optAllFiles, "Input file");
 
     // Parse everything
     try {
         app.parse(argc, argv);
+        if(optAllFiles.empty()) {
+            throw std::runtime_error("Please, provide an input file");
+        }
     } catch(const CLI::ParseError &e) {
         app.exit(e);
         throw;
     }
 
     // Input files (macros and data files)
-    result.AllFileNames = optAllFiles;
+    result.InputFileName = optAllFiles[0];
 
      // Output file name:
-    result.OutputBaseName = optOutputName;
+    result.OutputFileName = optOutputName;
 
+    result.BenchmarkName = benchmarkName;
     return result;
 }
 
-std::string ArgumentProcessor::mAppDescription =
-        "EIC Afterburner benchmarks";
+std::string ArgumentProcessor::get_description() {
+    return "EIC Afterburner benchmarks";
+}
+
 

cpp/benchmark/ArgumentProcessor.hh

@@ -25,9 +25,9 @@
 /// environment variables provided by users
 struct UserArguments
 {
-    std::vector<std::string> AllFileNames;          /// List of all arguments which are not flags=values (macro file names, input file names, etc.)
-    std::string SourceGenerator;                    /// The generator to use
-    std::string OutputBaseName;                     /// Desired output file(s) name (without suffix and extension)
+    std::string InputFileName;          /// List of all arguments which are not flags=values (macro file names, input file names, etc.)
+    std::string OutputFileName;                     /// Desired output file(s) name (without suffix and extension)
+    std::string BenchmarkName;                      /// Benchmark name
 };
 
 
@@ -40,7 +40,8 @@ public:
 private:
 
     /// Processes both program arguments and environment variables to build the right ProgramArguments
-    static std::string mAppDescription;
+
+    static std::string get_description();
 };
 
 

cpp/benchmark/CMakeLists.txt

@@ -7,6 +7,13 @@ find_package(HepMC3 REQUIRED)
 message(STATUS "BENCHMARKS HEPMC3_LIBRARIES ${HEPMC3_LIBRARIES}")
 message(STATUS "BENCHMARKS HEPMC3_INCLUDE_DIR ${HEPMC3_INCLUDE_DIR}")
 
+find_package(CLHEP)
+if(NOT CLHEP_FOUND)
+    set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/../cmake/Modules ${CMAKE_MODULE_PATH} )
+    find_package(CLHEP REQUIRED)
+endif()
+
+
 add_executable(benchmark ArgumentProcessor.hh ArgumentProcessor.cc main.cc)
 target_compile_features(benchmark
         PUBLIC

cpp/benchmark/main.cc

@@ -14,12 +14,18 @@
 #include "HepMC3/WriterHEPEVT.h"
 #include "HepMC3/ReaderFactory.h"
 
+#include <CLHEP/Vector/LorentzVector.h>
+#include <CLHEP/Vector/Boost.h>
+#include <CLHEP/Vector/Rotation.h>
+#include <CLHEP/Vector/EulerAngles.h>
+
+
 #include "ArgumentProcessor.hh"
 
 #include <afterburner/Afterburner.hh>
 #include <afterburner/AfterburnerConfig.hh>
 
-void convert_hepmc3_file(const std::string &input_file_name, const std::string &output_file_name, const ab::Afterburner &afterburner)
+void convert_hepmc3_file(const std::string &input_file_name, const std::string &output_file_name, ab::Afterburner& afterburner)
 {
     using namespace HepMC3;
 
@@ -46,6 +52,53 @@ void convert_hepmc3_file(const std::string &input_file_name, const std::string &
         if (last_event_number && evt.event_number()>last_event_number) continue;
         evt.set_run_info(input_file->run_info());
 
+        cout<<fixed;
+        std::cout<<"************\nORIGIN EVENT\n************\n";
+        Print::content(cout, evt);
+
+
+        // Run afterburner calculation
+        auto ab_result = afterburner.process_event();
+
+        // Rotate
+        auto axis = ab_result.rotation.axis();
+        auto delta = ab_result.rotation.delta();
+        auto euler = ab_result.rotation.eulerAngles();
+        auto theta = ab_result.rotation.theta();
+        auto psi = ab_result.rotation.psi();
+        auto phi = ab_result.rotation.phi();
+
+        // Instead of using HepMC rotate, we just rotate particle momentums (vertex rotation is done in the AB)
+        for ( auto p: evt.particles())
+        {
+            FourVector mom=p->momentum();
+            long double tempX=mom.x();
+            long double tempY=mom.y();
+            long double tempZ=mom.z();
+
+            CLHEP::Hep3Vector tmp(tempX, tempY, tempZ);
+            tmp = tmp.rotate(ab_result.rotation.eulerAngles());
+
+            FourVector temp(tmp.x(),tmp.y(),tmp.z(),mom.e());
+            p->set_momentum(temp);
+        }
+        std::cout<<"************\nAFTER ROTATION\n************\n";
+        Print::content(evt);
+
+        // Boost
+        auto boost = ab_result.boost.boostVector();
+        evt.boost(FourVector(boost.x(), boost.y(), boost.z(), 0));
+        std::cout<<"************\nAFTER BOOST\n************\n";
+        Print::content(evt);
+
+
+        // Translate
+        auto vtx = ab_result.vertex;
+        evt.shift_position_to(FourVector(vtx.x(), vtx.y(), vtx.z(), vtx.t()));
+        std::cout<<"************\nAFTER TRANSLATION\n************\n";
+        Print::content(evt);
+
+
         //Note the difference between ROOT and Ascii readers.
         // The former read GenRunInfo before first event and the later at the same time as first event.
         output_file->write_event(evt);
@@ -57,6 +110,7 @@ void convert_hepmc3_file(const std::string &input_file_name, const std::string &
             printf("Event limit reached:->events_parsed(%li) >= events_limit(%li)<-. Exit.\n",events_parsed , events_limit);
             break;
         }
+
     }
 
     if (input_file)   input_file->close();
@@ -74,9 +128,44 @@ int main(int argc, char** argv)
 
     // Afterburner instance
     ab::Afterburner afterburner;
-    afterburner.print();
 
-    convert_hepmc3_file(arguments.AllFileNames[0], arguments.OutputBaseName + ".hepmc", afterburner);
+    if(arguments.BenchmarkName == "crossing") {
+        // Only crossing angle benchmark
+
+        ab::AfterburnerConfig cfg;
+        //25mrad x-ing as in EIC CDR
+        const double EIC_hadron_crossing_angle = 25e-3;
+
+        cfg.beam_one.direction_theta = EIC_hadron_crossing_angle;   // beamA_theta
+        cfg.beam_one.direction_phi = 0;                             // beamA_phi
+        cfg.beam_two.direction_theta = M_PI;                        // beamB_theta
+        cfg.beam_two.direction_phi = 0;                             // beamB_phi
+
+        cfg.beam_one.divergence_hor = 0;
+        cfg.beam_one.divergence_ver = 0;
+        cfg.beam_two.divergence_hor = 0;
+        cfg.beam_two.divergence_ver = 0;
+
+        cfg.beam_one.z_shift_hor = 0;
+        cfg.beam_one.z_shift_ver = 0;
+        cfg.beam_two.z_shift_hor = 0;
+        cfg.beam_two.z_shift_ver = 0;
+
+        cfg.vertex_smear_width_x = 0;
+        cfg.vertex_smear_width_y = 0;
+        cfg.vertex_smear_width_z = 0;
+        cfg.vertex_smear_width_t = 0;
+
+        afterburner.set_config(cfg);
+    }
+    else
+    {
+        // Default configuration benchmark
+        // Nothing is here as everything is setup by default
+    }
+
+    afterburner.print();
+    convert_hepmc3_file(arguments.InputFileName, arguments.OutputFileName, afterburner);
 
     return EXIT_SUCCESS;
 }

python/gen_hepmc_event.py

@@ -13,10 +13,10 @@ def create_zonly_event(file_name):
     v1 = hm.GenVertex(hm.FourVector(1, 2, 3, 4))
     evt.add_vertex(v1)
 
-    p1 = hm.GenParticle(hm.FourVector(0, 0, 90, 90), 2212, 1)
+    p1 = hm.GenParticle(hm.FourVector(0, 0, 1000, 1000), 2212, 1)
     evt.add_particle(p1)
 
-    p2 = hm.GenParticle(hm.FourVector(0, 0, 7, 7), 11, 1)
+    p2 = hm.GenParticle(hm.FourVector(0, 0, 100, 100), 11, 1)
     evt.add_particle(p2)
 
     v1.add_particle_in(p1)
@@ -33,7 +33,7 @@ def read_events(file_name):
         raise "test_create_save_load.failed()"
 
     while True:
-        evt = hm.GenEvent()
+        evt = hm.GenEvent(hm.Units.GEV, hm.Units.CM)
         input_file.read_event(evt)
         if input_file.failed():
             break

python/zonly_event.hepmc

@@ -2,8 +2,8 @@ HepMC::Version 3.02.03
 HepMC::Asciiv3-START_EVENT_LISTING
 E 1 1 2
 U GEV MM
-P 1 0 2212 0.0000000000000000e+00 0.0000000000000000e+00 9.0000000000000000e+01 9.0000000000000000e+01 0.0000000000000000e+00 1
+P 1 0 2212 0.0000000000000000e+00 0.0000000000000000e+00 1.0000000000000000e+03 1.0000000000000000e+03 0.0000000000000000e+00 1
 V -1 0 [1] @ 1.0000000000000000e+00 2.0000000000000000e+00 3.0000000000000000e+00 4.0000000000000000e+00
-P 2 -1 11 0.0000000000000000e+00 0.0000000000000000e+00 7.0000000000000000e+00 7.0000000000000000e+00 0.0000000000000000e+00 1
+P 2 -1 11 0.0000000000000000e+00 0.0000000000000000e+00 1.0000000000000000e+02 1.0000000000000000e+02 0.0000000000000000e+00 1
 HepMC::Asciiv3-END_EVENT_LISTING