test: debug missing mass

.gitignore

@@ -1,4 +1,4 @@
-/.cache
+.cache
 /install*
 *.lund
 *.dat

beamo/meson.build

@@ -0,0 +1,16 @@
+project('pythia_test', 'cpp')
+
+pythia_config = find_program('pythia8-config')
+pythia_link_args = []
+foreach arg : run_command(pythia_config, '--libs', check: true).stdout().strip().split()
+  if not arg.contains('-rpath')
+    pythia_link_args += arg
+  endif
+endforeach
+pythia_dep = declare_dependency(
+  include_directories: run_command(pythia_config, '--includedir', check: true).stdout().strip(),
+  link_args: pythia_link_args,
+)
+
+exe = executable('pythia_test', 'test.cpp', dependencies: [pythia_dep])
+test('pythia_test', exe)

beamo/test.cpp

@@ -0,0 +1,36 @@
+#include <Pythia8/Pythia.h>
+
+int main(int argc, char** argv)
+{
+  Pythia8::Pythia p;
+  auto& evt = p.event;
+  auto& proc = p.process;
+  auto& info = p.info;
+
+  p.readString("Beams:frameType = 2");
+  p.readString("WeakBosonExchange:ff2ff(t:gmZ) = on");
+  p.readString("PhaseSpace:Q2Min = 0.0");
+  p.readString("PhaseSpace:mHatMin = 0.0");
+  p.readString("SpaceShower:dipoleRecoil = off");
+
+  p.readString("Beams:idA = 11");
+  p.readString("Beams:idB = 2212");
+  p.readString("Beams:eA = 10.6");
+  p.readString("Beams:eB = 0.0");
+  p.readString("Random:setSeed = on");
+  p.readString("Random:seed = 82");
+
+  // p.readString("Next:numberShowInfo = 500");
+  // p.readString("Next:numberShowProcess = 0");
+  // p.readString("Next:numberShowEvent = 500");
+
+  p.init();
+  for(int i=0; i<100; i++) {
+    if(!p.next()) continue;
+    evt.list(false, false, 10);
+    // proc.list(false, false, 10);
+    // info.list();
+  }
+
+  return 0;
+}

draw_mx.C

@@ -0,0 +1,12 @@
+// ninja && clas-stringspinner --seed 29877 --trig 10000 --config beam_test --cut-inclusive 11,211,-211|grep M_X| awk '{print $2}' > output.dat
+void draw_mx(TString file_name = "build/output.dat") {
+  auto tr = new TTree();
+  tr->ReadFile(file_name, "mx/D");
+  new TCanvas();
+  tr->Draw("mx");
+  new TCanvas();
+  auto mx_dist = new TH1D("mx_dist","m_x",200,0,2);
+  tr->Project("mx_dist", "mx");
+  mx_dist->Draw();
+  mx_dist->Fit("gaus","","",0.8,0.95);
+}

meson.build

@@ -16,6 +16,7 @@ stringspinner_proj = subproject('stringspinner')
 # dependencies
 fmt_dep           = dependency('fmt', method: 'pkg-config', version: '>= 9.1.0')
 stringspinner_dep = dependency('stringspinner', fallback: ['stringspinner', 'stringspinner_dep'])
+pythia_dep = dependency('pythia_force_subproject', fallback: ['stringspinner', 'stringspinner_pythia_dep'])
 
 # main executable
 add_project_arguments('-DCLAS_STRINGSPINNER_VERSION="' + meson.project_version() + '"', language: ['cpp'])
@@ -24,7 +25,7 @@ main_exe = executable(
   'src' / meson.project_name() + '.cpp',
   'src/config/src/common.cc',
   include_directories: include_directories('src/config'),
-  dependencies: [ fmt_dep, stringspinner_dep ],
+  dependencies: [ fmt_dep, stringspinner_dep, pythia_dep ],
   install: true,
 )
 

pythia.env

@@ -0,0 +1,5 @@
+#!/usr/bin/env fish
+
+set -x PYTHIA8DIR /home/dilks/j/pythia/install/share/Pythia8/xmldoc
+set -xp PATH /home/dilks/j/pythia/install/bin
+set -xp LD_LIBRARY_PATH /home/dilks/j/pythia/install/lib

src/clas-stringspinner.cpp

@@ -8,8 +8,10 @@
 
 // configurations
 #include "config/clas12.h"
+#include "config/beam_test.h"
 static std::map<std::string, std::function<void(Pythia8::Pythia&)>> CONFIG_MAP = {
-  {"clas12", config_clas12}
+  {"clas12", config_clas12},
+  {"beam_test", config_beam_test},
 };
 
 // constants
@@ -435,7 +437,7 @@ int main(int argc, char** argv)
   set_config(pyth, fmt::format("Beams:idA = {}", BEAM_PDG));
   set_config(pyth, fmt::format("Beams:idB = {}", target_pdg));
   set_config(pyth, fmt::format("Beams:eA = {}", beam_energy));
-  set_config(pyth, fmt::format("Beams:eB = {}", target_mass));
+  set_config(pyth, fmt::format("Beams:eB = {}", 0.0));
   //// seed
   set_config(pyth, "Random:setSeed = on");
   set_config(pyth, fmt::format("Random:seed = {}", seed));
@@ -492,6 +494,86 @@ int main(int argc, char** argv)
     }
     else cut_inclusive_passed = true;
 
+    ///////////////////////////////////////////////////////////////////////////////
+
+    // find the beam and target
+    int beam_idx = -1;
+    int target_idx = -1;
+    for(int idx = 0; idx < evt.size(); ++idx) {
+      auto const& par = evt[idx];
+      if(par.status() == -12) {
+        switch(par.id()) {
+          case 11: beam_idx = idx; break;
+          case 2212: target_idx = idx; break;
+        }
+      }
+    }
+    if(beam_idx < 0 || target_idx < 0)
+      continue;
+
+    // pythia's frame beam and target momenta
+    auto beam_E   = evt[beam_idx].e(); // pythia's frame
+    auto beam_p   = evt[beam_idx].pz();
+    auto target_E = evt[target_idx].e();
+    auto target_p = evt[target_idx].pz();
+    fmt::print("ROOT {:5.12} {:5.12} {:5.12}\n", beam_p, target_p, (evt[beam_idx].p()+evt[target_idx].p()).mCalc());
+
+    // lab frame beam and target momenta (H is energy, q is pz)
+    double beam_H   = beam_energy;
+    double beam_q   = std::sqrt(std::pow(beam_energy,2) - std::pow(pdt.constituentMass(11),2));
+    double target_H = target_mass;
+    double target_q = 0.0;
+
+    // solve for beta and gamma of the lorentz transformation from pythia frame to lab frame
+    auto get_boost = [](double E, double p, double H, double q, double& beta, double& gamma, double &a, double &b) {
+      double denom = (std::pow(E,2) - std::pow(p,2));
+      a = (E*H - p*q) / denom;
+      b = (H*p - E*q) / denom;
+      // fmt::print("a={} b={}\n", a, b);
+      gamma = a;
+      beta = b / a;
+    };
+    double beam_gamma,   beam_beta, beam_a, beam_b;
+    double target_gamma, target_beta, target_a, target_b;
+    get_boost(beam_E,   beam_p,   beam_H,   beam_q,   beam_gamma,   beam_beta, beam_a, beam_b);
+    get_boost(target_E, target_p, target_H, target_q, target_gamma, target_beta, target_a, target_b);
+
+    // transform the event from pythia frame to lab frame
+    // evt.bst(0., 0., -target_beta, target_gamma);
+    // evt.bst(0., 0., -beam_beta, beam_gamma);
+
+    // fmt::print("target_beta={} target_gamma={}\n", target_beta, target_gamma);
+    // fmt::print("beam_beta={} beam_gamma={}\n", beam_beta, beam_gamma);
+
+    // auto apply_boost = [beta=target_beta, gamma=target_gamma](Pythia8::Vec4 v) -> Pythia8::Vec4 {
+    //   return Pythia8::Vec4(
+    //       v.px(),
+    //       v.py(),
+    //       beta * gamma * v.e() + gamma * v.pz(),
+    //       gamma * v.e() + beta * gamma * v.pz()
+    //       );
+    // };
+    auto apply_boost = [a=beam_a, b=target_b](Pythia8::Vec4 v) -> Pythia8::Vec4 {
+      return Pythia8::Vec4(
+          v.px(),
+          v.py(),
+          -b*v.e() + a*v.pz(),
+          a*v.e() - b*v.pz()
+          );
+    };
+
+    // auto print_vec = [](std::string const& s, Pythia8::Vec4 v) {
+    //   fmt::print("{}: {:5.12}  {:5.12}  {:5.12}  {:5.12}\n", s, v.px(), v.py(), v.pz(), v.e());
+    // };
+    // print_vec("  beam", evt[beam_idx].p());
+    // print_vec("  beam", apply_boost(evt[beam_idx].p()));
+    // print_vec("target", evt[target_idx].p());
+    // print_vec("target", apply_boost(evt[target_idx].p()));
+
+    evt.list(false, false, 12);
+
+    ///////////////////////////////////////////////////////////////////////////////
+
     // loop over particles
     std::vector<LundParticle> lund_particles;
     Verbose("Particles:");
@@ -600,6 +682,69 @@ int main(int argc, char** argv)
         break;
     }
 
+
+    //////////////////////////////////////////////////////////////////////////////////
+
+    /*
+    for(int idx = 0; idx < evt.size(); ++idx) {
+      auto const& par = evt[idx];
+      if(par.isFinal()) {
+        auto p = std::hypot(par.px(),par.py(),par.pz());
+        auto mass = std::sqrt( std::pow(par.e(),2) - std::pow(p,2) );
+        fmt::print("TEST: pdg={} p={} e()={} m()={} mass={} diff={}\n", par.id(), p, par.e(), par.m(), mass, par.m() - mass);
+      }
+    }
+    */
+
+    // find the scattered electron
+    int ele_idx = -1;
+    double ele_p = -10000;
+    for(int idx = 0; idx < evt.size(); ++idx) {
+      auto const& par = evt[idx];
+      if(par.isFinal() && par.id() == 11) {
+        auto this_p = std::hypot(par.px(), par.py(), par.pz());
+        if(this_p > ele_p) {
+          ele_p = par.e();
+          ele_idx = idx;
+        }
+      }
+    }
+    // fmt::print("ele_idx = {}\n", ele_idx);
+    if(ele_idx < 0) continue;
+
+    // loop over pi+ pi- dihadrons
+    for(int idxA = 0; idxA < evt.size(); ++idxA) {
+      auto const& parA = evt[idxA];
+      if(parA.isFinal() && parA.id() == 211) {
+        for(int idxB = 0; idxB < evt.size(); ++idxB) {
+          auto const& parB = evt[idxB];
+          if(parB.isFinal() && parB.id() == -211) {
+
+            // auto p_beam = evt[beam_idx].p();
+            // auto p_target = evt[target_idx].p();
+            // auto p_ele = evt[ele_idx].p();
+            // auto p_pip = parA.p();
+            // auto p_pim = parB.p();
+
+            auto p_beam = apply_boost(evt[beam_idx].p());
+            auto p_target = apply_boost(evt[target_idx].p());
+            auto p_ele = apply_boost(evt[ele_idx].p());
+            auto p_pip = apply_boost(parA.p());
+            auto p_pim = apply_boost(parB.p());
+
+            // calculate missing mass
+            auto vecW = p_beam + p_target - p_ele;
+            auto vecPh = p_pip + p_pim;
+            auto M_X = (vecW - vecPh).mCalc();
+            fmt::print("M_X {}\n", M_X);
+          }
+        }
+      }
+    }
+
+    //////////////////////////////////////////////////////////////////////////////////
+
+
   } // end EVENT LOOP
 
   fmt::print("GENERATED LUND FILE: {}\n", out_file);

src/config/beam_test.h

@@ -0,0 +1,22 @@
+#include "src/common.h"
+static void config_beam_test(Pythia8::Pythia& pyth) {
+
+  // the beams are back-to-back, but with different energies; this is the recommended setting for fixed target
+  set_config(pyth, "Beams:frameType = 2");
+
+  // interaction mechanism: $\gamma*/Z^0$ $t$-channel exchange, with full interference between $\gamma*$ and $Z^0$
+  set_config(pyth, "WeakBosonExchange:ff2ff(t:gmZ) = on");
+
+  // phase-space cuts
+  set_config(pyth, "PhaseSpace:Q2Min = 0.0");   // minimum Q^2
+  set_config(pyth, "PhaseSpace:mHatMin = 0.0"); // minimum invariant mass (for low-x)
+
+  // set dipole recoil; turning this on doesn't appear to change the kinematic distributions noticeably
+  set_config(pyth, "SpaceShower:dipoleRecoil = off");
+
+  // handle event printouts
+  set_config(pyth, "Next:numberShowInfo = 0");
+  set_config(pyth, "Next:numberShowProcess = 0");
+  set_config(pyth, "Next:numberShowEvent = 0");
+
+}

subprojects/stringspinner/meson.build

@@ -65,7 +65,7 @@ if get_option('install_example')
     'stringspinner-example',
     sources: ['stringspinner/dis.cc'],
     link_with: stringspinner_libs,
-    dependencies: [ pythia_dep ],
+    dependencies: [ stringspinner_pythia_dep ],
     install: true,
     cpp_args: ['-Wno-sign-compare', '-Wno-unused-variable'],
   )