Beam.h 2 KB
Newer Older
1
2
3
4
5
6
7
8
#pragma once

#include "Math/Vector4D.h"
using ROOT::Math::PxPyPzEVector;

#include "edm4hep/MCParticleCollection.h"
#include "eicd/ReconstructedParticleCollection.h"

9
namespace Jug::Base::Beam {
10
11

  template<class collection>
12
  auto find_first_with_pdg(
13
14
15
16
17
18
19
20
21
22
23
24
25
      const collection& parts,
      const std::set<int32_t>& pdg) {
    collection c;
    for (const auto& p: parts) {
      if (pdg.count(p.getPDG()) > 0) {
        c.push_back(p);
        break;
      }
    }
    return c;
  }

  template<class collection>
26
  auto find_first_with_status_pdg(
27
28
29
30
31
32
33
34
35
36
37
38
39
      const collection& parts,
      const std::set<int32_t>& status,
      const std::set<int32_t>& pdg) {
    collection c;
    for (const auto& p: parts) {
      if (status.count(p.getGeneratorStatus()) > 0 && pdg.count(p.getPDG()) > 0) {
        c.push_back(p);
        break;
      }
    }
    return c;
  }

40
  inline auto find_first_beam_electron(const edm4hep::MCParticleCollection& mcparts) {
41
42
43
    return find_first_with_status_pdg(mcparts, {4}, {11});
  }

44
  inline auto find_first_beam_hadron(const edm4hep::MCParticleCollection& mcparts) {
45
46
47
    return find_first_with_status_pdg(mcparts, {4}, {2212, 2112});
  }

48
  inline auto find_first_scattered_electron(const edm4hep::MCParticleCollection& mcparts) {
49
50
51
    return find_first_with_status_pdg(mcparts, {1}, {11});
  }

52
  inline auto find_first_scattered_electron(const eicd::ReconstructedParticleCollection& rcparts) {
53
54
55
    return find_first_with_pdg(rcparts, {11});
  }

56
  inline
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
  PxPyPzEVector
  round_beam_four_momentum(
      const edm4hep::Vector3f& p_in,
      const float mass,
      const std::vector<float>& pz_set,
      const float crossing_angle = 0.0) {
    PxPyPzEVector p_out;
    for (const auto& pz : pz_set) {
      if (fabs(p_in.z / pz - 1) < 0.1) {
        p_out.SetPz(pz);
        break;
      }
    }
    p_out.SetPx(p_out.Pz() * sin(crossing_angle));
    p_out.SetPz(p_out.Pz() * cos(crossing_angle));
    p_out.SetE(std::hypot(p_out.Px(), p_out.Pz(), mass));
    return p_out;
  }

76
} // namespace Jug::Base::Beam