InclusiveKinematicsDA.cpp 8.16 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
#include "Gaudi/Algorithm.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiAlg/Producer.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/PhysicalConstants.h"
#include <algorithm>
#include <cmath>

#include "JugBase/IParticleSvc.h"
#include "JugBase/DataHandle.h"

13
14
#include "JugBase/Utilities/Beam.h"
#include "JugBase/Utilities/Boost.h"
15

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

19
// Event Model related classes
20
#include "edm4hep/MCParticleCollection.h"
21
#include "eicd/MCRecoParticleAssociationCollection.h"
22
23
24
25
26
#include "eicd/ReconstructedParticleCollection.h"
#include "eicd/InclusiveKinematicsCollection.h"

namespace Jug::Reco {

Sylvester Joosten's avatar
Sylvester Joosten committed
27
class InclusiveKinematicsDA : public GaudiAlgorithm {
Wouter Deconinck's avatar
Wouter Deconinck committed
28
private:
29
30
  DataHandle<edm4hep::MCParticleCollection> m_inputMCParticleCollection{
    "MCParticles",
31
32
    Gaudi::DataHandle::Reader,
    this};
Sylvester Joosten's avatar
Sylvester Joosten committed
33
  DataHandle<eicd::ReconstructedParticleCollection> m_inputParticleCollection{
34
35
36
    "ReconstructedParticles",
    Gaudi::DataHandle::Reader,
    this};
37
38
39
40
  DataHandle<eicd::MCRecoParticleAssociationCollection> m_inputParticleAssociation{
    "MCRecoParticleAssociation",
    Gaudi::DataHandle::Reader,
    this};
Sylvester Joosten's avatar
Sylvester Joosten committed
41
  DataHandle<eicd::InclusiveKinematicsCollection> m_outputInclusiveKinematicsCollection{
42
43
44
45
46
47
48
    "InclusiveKinematicsDA",
    Gaudi::DataHandle::Writer,
    this};

  Gaudi::Property<double> m_crossingAngle{this, "crossingAngle", -0.025 * Gaudi::Units::radian};

  SmartIF<IParticleSvc> m_pidSvc;
49
  double m_proton{0}, m_neutron{0}, m_electron{0};
50

Wouter Deconinck's avatar
Wouter Deconinck committed
51
public:
52
  InclusiveKinematicsDA(const std::string& name, ISvcLocator* svcLoc)
Sylvester Joosten's avatar
Sylvester Joosten committed
53
      : GaudiAlgorithm(name, svcLoc) {
54
    declareProperty("inputMCParticles", m_inputMCParticleCollection, "MCParticles");
55
    declareProperty("inputParticles", m_inputParticleCollection, "ReconstructedParticles");
56
    declareProperty("inputAssociations", m_inputParticleAssociation, "MCRecoParticleAssociation");
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
    declareProperty("outputData", m_outputInclusiveKinematicsCollection, "InclusiveKinematicsDA");
  }

  StatusCode initialize() override {
    if (GaudiAlgorithm::initialize().isFailure())
      return StatusCode::FAILURE;

    m_pidSvc = service("ParticleSvc");
    if (!m_pidSvc) {
      error() << "Unable to locate Particle Service. "
              << "Make sure you have ParticleSvc in the configuration."
              << endmsg;
      return StatusCode::FAILURE;
    }
    m_proton = m_pidSvc->particle(2212).mass;
72
    m_neutron = m_pidSvc->particle(2112).mass;
73
74
75
76
77
78
79
    m_electron = m_pidSvc->particle(11).mass;

    return StatusCode::SUCCESS;
  }

  StatusCode execute() override {
    // input collections
80
81
82
    const auto& mcparts = *(m_inputMCParticleCollection.get());
    const auto& rcparts = *(m_inputParticleCollection.get());
    const auto& rcassoc = *(m_inputParticleAssociation.get());
83
84
85
    // output collection
    auto& out_kinematics = *(m_outputInclusiveKinematicsCollection.createAndPut());

86
    // Get incoming electron beam
87
    const auto ei_coll = Jug::Base::Beam::find_first_beam_electron(mcparts);
88
89
90
    if (ei_coll.size() == 0) {
      if (msgLevel(MSG::DEBUG)) {
        debug() << "No beam electron found" << endmsg;
91
      }
92
      return StatusCode::SUCCESS;
93
    }
94
    const PxPyPzEVector ei(
95
      Jug::Base::Beam::round_beam_four_momentum(
96
97
98
99
100
101
102
        ei_coll[0].getMomentum(),
        m_electron,
        {-5.0, -10.0, -18.0},
        0.0)
      );

    // Get incoming hadron beam
103
    const auto pi_coll = Jug::Base::Beam::find_first_beam_hadron(mcparts);
104
    if (pi_coll.size() == 0) {
105
      if (msgLevel(MSG::DEBUG)) {
106
        debug() << "No beam hadron found" << endmsg;
107
108
109
      }
      return StatusCode::SUCCESS;
    }
110
    const PxPyPzEVector pi(
111
      Jug::Base::Beam::round_beam_four_momentum(
112
113
114
115
116
117
118
        pi_coll[0].getMomentum(),
        pi_coll[0].getPDG() == 2212 ? m_proton : m_neutron,
        {41.0, 100.0, 275.0},
        m_crossingAngle)
      );

    // Get first scattered electron
119
    const auto ef_coll = Jug::Base::Beam::find_first_scattered_electron(mcparts);
120
    if (ef_coll.size() == 0) {
121
      if (msgLevel(MSG::DEBUG)) {
122
        debug() << "No truth scattered electron found" << endmsg;
123
124
125
      }
      return StatusCode::SUCCESS;
    }
126
127
128
129
130
131
132
133
134
135
136
137
    // Associate first scattered electron with reconstructed electrons
    //const auto ef_assoc = std::find_if(
    //  rcassoc.begin(),
    //  rcassoc.end(),
    //  [&ef_coll](const auto& a){ return a.getSimID() == ef_coll[0].id(); });
    auto ef_assoc = rcassoc.begin();
    for (; ef_assoc != rcassoc.end(); ++ef_assoc) {
      if (ef_assoc->getSimID() == ef_coll[0].id()) {
        break;
      }
    }
    if (!(ef_assoc != rcassoc.end())) {
Wouter Deconinck's avatar
Wouter Deconinck committed
138
      if (msgLevel(MSG::DEBUG)) {
139
        debug() << "Truth scattered electron not in reconstructed particles" << endmsg;
Wouter Deconinck's avatar
Wouter Deconinck committed
140
      }
141
      return StatusCode::SUCCESS;
Wouter Deconinck's avatar
Wouter Deconinck committed
142
    }
143
144
    const auto ef_rc{ef_assoc->getRec()};
    const auto ef_rc_id{ef_rc.id()};
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165

    // Loop over reconstructed particles to get all outgoing particles
    // ----------------------------------------------------------------- 
    // Right now, everything is taken from Reconstructed particles branches.
    // 
    // This means the tracking detector is used for charged particles to caculate the momentum,
    // and the magnitude of this momentum plus the true PID to calculate the energy.
    // No requirement is made that these particles produce a hit in any other detector
    // 
    // Using the Reconstructed particles branches also means that the reconstruction for neutrals is done using the
    // calorimeter(s) information for the energy and angles, and then using this energy and the true PID to get the
    // magnitude of the momentum.
    // -----------------------------------------------------------------

    //Sums in colinear frame
    double pxsum = 0;
    double pysum = 0;
    double pzsum = 0;
    double Esum = 0;
    double theta_e = 0;

166
    // Get boost to colinear frame
167
    auto boost = Jug::Base::Boost::determine_boost(ei, pi);
168
169
170
171
172
173
174

    for(const auto& p : rcparts) {
      // Get the scattered electron index and angle
      if (p.id() == ef_rc_id) {
        // Lorentz vector in lab frame
        PxPyPzEVector e_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
        // Boost to colinear frame
175
        PxPyPzEVector e_boosted = Jug::Base::Boost::apply_boost(boost, e_lab);
176
177
178
179
180
181
182
183

        theta_e = e_boosted.Theta();

      // Sum over all particles other than scattered electron
      } else {
        // Lorentz vector in lab frame
        PxPyPzEVector hf_lab(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy());
        // Boost to colinear frame
184
        PxPyPzEVector hf_boosted = Jug::Base::Boost::apply_boost(boost, hf_lab);
185
186
187
188
189
190

        pxsum += hf_boosted.Px();
        pysum += hf_boosted.Py();
        pzsum += hf_boosted.Pz();
        Esum += hf_boosted.E();
      }
191
192
193
194
195
196
197
    }

    // DIS kinematics calculations
    auto sigma_h = Esum - pzsum;
    auto ptsum = sqrt(pxsum*pxsum + pysum*pysum);
    auto theta_h = 2.*atan(sigma_h/ptsum);

198
199
    // If no scattered electron was found
    if (sigma_h <= 0) {
200
      if (msgLevel(MSG::DEBUG)) {
201
        debug() << "No scattered electron found" << endmsg;
202
      }
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
      return StatusCode::SUCCESS;
    }

    // Calculate kinematic variables
    const auto y_da = tan(theta_h/2.) / ( tan(theta_e/2.) + tan(theta_h/2.) );
    const auto Q2_da = 4.*ei.energy()*ei.energy() * ( 1. / tan(theta_e/2.) ) * ( 1. / (tan(theta_e/2.) + tan(theta_h/2.)) );
    const auto x_da = Q2_da / (4.*ei.energy()*pi.energy()*y_da);
    const auto nu_da = Q2_da / (2.*m_proton*x_da);
    const auto W_da = sqrt(m_proton*m_proton + 2*m_proton*nu_da - Q2_da);
    auto kin = out_kinematics.create(x_da, Q2_da, W_da, y_da, nu_da);
    kin.setScat(ef_rc);

    // Debugging output
    if (msgLevel(MSG::DEBUG)) {
      debug() << "pi = " << pi << endmsg;
      debug() << "ei = " << ei << endmsg;
      debug() << "x,Q2,W,y,nu = "
              << kin.getX() << ","
              << kin.getQ2() << ","
              << kin.getW() << ","
              << kin.getY() << ","
              << kin.getNu()
              << endmsg;
226
227
228
229
230
231
    }

    return StatusCode::SUCCESS;
  }
};

232
// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
233
234
235
DECLARE_COMPONENT(InclusiveKinematicsDA)

} // namespace Jug::Reco