Added --exit-ca flag

README.md

@@ -59,6 +59,13 @@ abconv my.hepmc -f hepmc2 --ab-off --plot-off
 | -l,--limit UINT      | Limit number of events to process. (Shutdown after this number of parsed events)|
 | --ab-off             | No afterburner is applied|
 | --plot-off           | Don't produce validation plots|
+|--exit-ca             | Check existing crossing angle and exit if CA>1mrad (1) |
+
+> (1) How `--exit-ca` works exactly: when afterburner processes the first event it checks if 
+> it has 2 beam particles (fails with non zero code if not) and then calculates their crossing angle
+> with the `--exit-ca` is set and crossing angle absolute value is > 1mrad program ends with 0 code.
+> This method is very coarse as if in a source file the crossing angle is 0 but all beam effects do exist,
+> beam effects will be applied twice
 
 ## Validation
 

cpp/abconv/ArgumentProcessor.cc

@@ -16,13 +16,13 @@ UserArguments ArgumentProcessor::Process(int argc, char **argv)
     std::string config_file;
     bool ab_off = false;
     bool plot_off = false;
+    bool check_ca = false;
     ulong process_limit = 0;
     std::string input_format = "auto";
     std::string output_format = "hepmc3";
     long ev_start = 0;
     long ev_end = 0;
 
-
     app.add_option("-o,--output", output_base_name, "Base name for Output files ((!) no extension)");
     app.add_option("-c,--config", config_file, "Beams configuration file");
     app.add_option("-i, --in-format", input_format, "Input format: auto [default], hepmc2, hepmc3, hpe, lhef, gz, treeroot, root");
@@ -34,6 +34,7 @@ UserArguments ArgumentProcessor::Process(int argc, char **argv)
 
     app.add_flag("--ab-off", ab_off, "No afterburner is applied");
     app.add_flag("--plot-off", plot_off, "Don't produce validation plots");
+    app.add_flag("--exit-ca", check_ca, "Check existing crossing angle and exit if CA>1mrad");
 
     app.add_option("input file", optAllFiles, "Input file");
 
@@ -74,6 +75,7 @@ UserArguments ArgumentProcessor::Process(int argc, char **argv)
     result.EventProcessLimit = process_limit;
     result.StartEventIndex = ev_start;
     result.EndEventIndex = ev_end;
+    result.ExitOnCrossingAngle = check_ca;
     return result;
 }
 

cpp/abconv/ArgumentProcessor.hh

@@ -37,6 +37,7 @@ struct UserArguments
 
     long StartEventIndex;               /// Start event index (all skipped before that)
     long EndEventIndex;                 /// End event index (stop after that)
+    bool ExitOnCrossingAngle;           /// Check existing crossing angle and exit if CA>1mrad
 };
 
 

cpp/abconv/Converter.cc

@@ -35,7 +35,8 @@ void ab::abconv::Converter::convert() {
 
         if(0 == events_processed) {
             // The first event, determine beams configuration
-            get_beams_setup(evt);
+            bool must_exit = check_beams_setup(evt);
+            if(must_exit) return;
 
         }
 
@@ -130,7 +131,7 @@ void ab::abconv::Converter::print_processed_events(long count) {
     if(count % div == 0 ) printf("Events parsed: %li\n", count);
 }
 
-void ab::abconv::Converter::get_beams_setup(const HepMC3::GenEvent& event) {
+bool ab::abconv::Converter::check_beams_setup(const HepMC3::GenEvent& event) {
 
     std::vector<std::shared_ptr<const HepMC3::GenParticle>> beam_particles;
 
@@ -160,9 +161,13 @@ void ab::abconv::Converter::get_beams_setup(const HepMC3::GenEvent& event) {
     printf("   pdg: %d   p: %.1f, %.1f, %.1f   e:%.1f\n", pdg_one, mom_one.x(), mom_one.y(), mom_one.z(), mom_one.e());
     printf("   pdg: %d   p: %.1f, %.1f, %.1f   e:%.1f\n", pdg_two, mom_two.x(), mom_two.y(), mom_two.z(), mom_two.e());
     double one_dot_two = mom_one.x()*mom_two.x() + mom_one.y()*mom_two.y() + mom_one.z()*mom_two.z();
-    double cangle = std::acos(one_dot_two/(mom_one.length()*mom_two.length())) - M_PI;
-    printf("   crossing angle: %.0f [mrad]\n", cangle*1000);
+    double crossing_angle = std::acos(one_dot_two / (mom_one.length() * mom_two.length())) - M_PI;
+    printf("   crossing angle: %.0f [mrad]\n", crossing_angle * 1000);
     printf("=========================\n");
-}
-
 
+    if(_exit_on_ca && fabs(crossing_angle)>0.001) {
+        printf("(!) Existing crossing angle > 1mrad and --exit-ca flag is used. Exiting\n");
+        return true;
+    }
+    return false;
+}
\ No newline at end of file

cpp/abconv/Converter.hh

@@ -40,6 +40,11 @@ namespace ab{
                 this->_ap_callback = std::move(after_callback);
             }
 
+            void set_exit_on_ca(bool b) {
+                _exit_on_ca = b;
+            }
+
+
         private:
             std::shared_ptr<HepMC3::Reader> _reader;
             std::shared_ptr<HepMC3::Writer> _writer;
@@ -48,10 +53,13 @@ namespace ab{
             uint64_t _first_event_number=0;
             uint64_t _last_event_number=0;
             uint64_t _events_limit=0;
+            bool _exit_on_ca=false;
+
             std::function<void(HepMC3::GenEvent&)> _ap_callback=nullptr;
             static void print_processed_events(long count);
 
-            void get_beams_setup(const HepMC3::GenEvent& event);
+            bool check_beams_setup(const HepMC3::GenEvent& event);
+
         };
     }
 }

cpp/abconv/main.cc

@@ -59,6 +59,7 @@ int main(int argc, char** argv)
     conv.set_first_event_number(arg.StartEventIndex);
     conv.set_last_event_number(arg.EndEventIndex);
     conv.set_after_process_callback(plotting_callback);
+    conv.set_exit_on_ca(arg.ExitOnCrossingAngle);
 
     // CONVERT!!!
     conv.convert();

cpp/benchmark/Histogrammer.cc

@@ -11,8 +11,6 @@
 #include <HepMC3/GenParticle.h>
 #include <HepMC3/GenVertex.h>
 
-extern double vertex_x, vertex_y, vertex_z, vertex_t;
-
 
 Histogrammer::Histogrammer(std::string output_file):
         _output_file_path(std::move(output_file)),
@@ -140,10 +138,10 @@ void Histogrammer::process_event(HepMC3::GenEvent &event) {
 
     }
 
-    vtxX->Fill(vertex_x);
-    vtxY->Fill(vertex_y);
-    vtxZ->Fill(vertex_z);
-    vtxT->Fill(vertex_t);
+    vtxX->Fill(vtx_one->position().x());
+    vtxY->Fill(vtx_one->position().y());
+    vtxZ->Fill(vtx_one->position().z());
+    vtxT->Fill(vtx_one->position().t());
 
     vtx2X->Fill(vtx_two->position().x());
     vtx2Y->Fill(vtx_two->position().y());