SmearedFarForwardParticles.cpp 14.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include <algorithm>
#include <cmath>
#include <fmt/format.h>

#include "Gaudi/Algorithm.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiAlg/Producer.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiKernel/RndmGenerators.h"

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

// Event Model related classes
#include "dd4pod/Geant4ParticleCollection.h"
#include "eicd/ReconstructedParticleCollection.h"
#include "eicd/VectorPolar.h"

19
20
21
22
23
24
25
26
27
namespace {
  enum DetectorTags {
    kTagB0 = 1,
    kTagRP = 2,
    kTagOMD = 3,
    kTagZDC = 4
  };
}

28
namespace Jug::Fast {
29

30
class SmearedFarForwardParticles : public GaudiAlgorithm, AlgorithmIDMixin<> {
31
32
public:
  DataHandle<dd4pod::Geant4ParticleCollection> m_inputParticles{"inputMCParticles", Gaudi::DataHandle::Reader, this};
33
  DataHandle<eic::ReconstructedParticleCollection> m_outputParticles{"SmearedFarForwardParticles", Gaudi::DataHandle::Writer,
34
                                                                     this};
35

36
37
38
39
  Gaudi::Property<bool> m_enableZDC{this, "enableZDC", true};
  Gaudi::Property<bool> m_enableB0{this, "enableB0", true};
  Gaudi::Property<bool> m_enableRP{this, "enableRP", true};
  Gaudi::Property<bool> m_enableOMD{this, "enableOMD", true};
40
41
42

  // Beam energy, only used to determine the RP/OMD momentum ranges
  Gaudi::Property<double> m_ionBeamEnergy{this, "ionBeamEnergy", 100.};
43
  // RP default to 10-on-100 setting
44
  // Pz > 60% of beam energy (60% x 100GeV = 60GeV)
45
46
47
  // theta from 0.2mrad -> 5mrad
  Gaudi::Property<double> m_thetaMinRP{this, "thetaMinRP", 0.2e-3};
  Gaudi::Property<double> m_thetaMaxRP{this, "thetaMaxRP", 5e-3};
48
  Gaudi::Property<double> m_pMinRigidityRP{this, "pMinRigidityRP", 0.60};
49
50
51
52
53
54
55
56
57
58
59
  // B0
  Gaudi::Property<double> m_thetaMinB0{this, "thetaMinB0", 6.0e-3};
  Gaudi::Property<double> m_thetaMaxB0{this, "thetaMaxB0", 20.0e-3};
  // OMD default to 10-on-100 setting
  // 25% < P/Ebeam < 60% of beam energy (25% x 100GeV = 25GeV and 60% x 100GeV = 60GeV)
  // Angles both given for the small angle full-acceptance part,
  // and for the larger angle part where we only measure |phi| > rad
  Gaudi::Property<double> m_thetaMinFullOMD{this, "thetaMinFullOMD", 0.};
  Gaudi::Property<double> m_thetaMaxFullOMD{this, "thetaMaxFullOMD", 2e-3};
  Gaudi::Property<double> m_thetaMinPartialOMD{this, "thetaMinPartialOMD", 2.0e-3};
  Gaudi::Property<double> m_thetaMaxPartialOMD{this, "thetaMaxPartialOMD", 5.0e-3};
60
61
  Gaudi::Property<double> m_pMinRigidityOMD{this, "pMinRigidityOMD", 0.25};
  Gaudi::Property<double> m_pMaxRigidityOMD{this, "pMaxRigidityOMD", 0.60};
62

63
  // Crossing angle, set to -25mrad
64
  Gaudi::Property<double> m_crossingAngle{this, "crossingAngle", -0.025}; //-0.025}; -- causes double rotation with afterburner
65

66
67
  Rndm::Numbers m_gaussDist;

68
69
70
71
  // Monte Carlo particle source identifier
  const int32_t m_kMonteCarloSource{uniqueID<int32_t>("mcparticles")};

  private:
72
  using RecData = eic::ReconstructedParticle;
73
74

  public:
75
  SmearedFarForwardParticles(const std::string& name, ISvcLocator* svcLoc)
76
      : GaudiAlgorithm(name, svcLoc), AlgorithmIDMixin(name, info()) {
77
78
    declareProperty("inputMCParticles", m_inputParticles, "mcparticles");
    declareProperty("outputParticles", m_outputParticles, "ReconstructedParticles");
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
  }
  StatusCode initialize() override {
    if (GaudiAlgorithm::initialize().isFailure())
      return StatusCode::FAILURE;

    IRndmGenSvc* randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
    // use 0 for mean and 1 for standard deviation. Can rescale appropriately for the
    // different subsystems
    StatusCode sc = m_gaussDist.initialize(randSvc, Rndm::Gauss(0.0, 1.0));
    if (!sc.isSuccess()) {
      return StatusCode::FAILURE;
    }
    return StatusCode::SUCCESS;
  }
  StatusCode execute() override {
    const auto& mc = *(m_inputParticles.get());
    auto& rc       = *(m_outputParticles.createAndPut());

97
    std::vector<std::vector<RecData>> rc_parts;
98
99
100
101
102
103
104
105
106
107
108
109
110
111
    if (m_enableZDC) {
      rc_parts.push_back(zdc(mc));
    }
    if (m_enableRP) {
      rc_parts.push_back(rp(mc));
    }
    if (m_enableB0) {
      rc_parts.push_back(b0(mc));
    }
    if (m_enableOMD) {
      rc_parts.push_back(omd(mc));
    }
    for (const auto& det : rc_parts) {
      for (const auto& part : det) {
112
        rc.push_back(part);
113
114
115
116
117
118
119
120
      }
    }
    return StatusCode::SUCCESS;
  }

private:
  // ZDC smearing as in eic_smear
  // https://github.com/eic/eicsmeardetectors/blob/9a1831dd97bf517b80a06043b9ee4bfb96b483d8/SmearMatrixDetector_0_1_FF.cxx#L224
121
122
  std::vector<RecData> zdc(const dd4pod::Geant4ParticleCollection& mc) {
    std::vector<RecData> rc;
123
    for (const auto& part : mc) {
124
      if (part.genStatus() > 1) {
Sylvester Joosten's avatar
Sylvester Joosten committed
125
126
127
        if (msgLevel(MSG::DEBUG)) {
          debug() << "ignoring particle with genStatus = " << part.genStatus() << endmsg;
        }
128
129
130
        continue;
      }
      // only detect neutrons and photons
131
      const auto mom_ion = rotateLabToIonDirection(part.ps());
132
133
134
135
      if (part.pdgID() != 2112 && part.pdgID() != 22) {
        continue;
      }
      // only 0-->4.5 mrad
136
      if (mom_ion.theta() > 4.5 / 1000.) {
137
138
139
        continue;
      }
      const double E    = std::hypot(part.ps().mag(), part.mass());
140
141
142
143
144
145
146
147
148
149
150
      double conTerm = 0.05; //default 5%
      double stoTerm = 0.5;  //default 50%
      double angTerm = 0.003; //3mrad

      if(part.pdgID() == 2112){
        conTerm = 0.05; //default 5%
        stoTerm = 0.5;  //default 50%
        angTerm = 0.003; //3mrad
      }
      else if(part.pdgID() == 22){  //EMCAL expected to have slightly better performance
        conTerm = 0.03; //default 3%
151
152
        stoTerm = 0.10;  //default 10% for WSciFi
        angTerm = 0.001; //1mrad is the detault for the block size
153
154
155
      }

      const double dE   = sqrt((conTerm * E) * (conTerm * E) + stoTerm * stoTerm * E) * m_gaussDist(); //50%/SqrtE + 5%
156
      const double Es   = E + dE;
157
      const double th   = mom_ion.theta();
158
      const double dth  = (angTerm / sqrt(E)) * m_gaussDist();
159
      const double ths  = th + dth;
160
      const double phi  = mom_ion.phi();
161
162
163
      const double dphi = 0;
      const double phis = phi + dphi;
      const double moms = sqrt(Es * Es - part.mass() * part.mass());
164
165
      const eic::VectorPolar mom3s_ion{moms, ths, phis};
      const auto mom3s = rotateIonToLabDirection(mom3s_ion);
166
167
168
169
170
171
      eic::ReconstructedParticle rec_part;
      rec_part.ID({part.ID(), algorithmID()});
      rec_part.p(mom3s);
      rec_part.v({part.vs().x, part.vs().y, part.vs().z});
      rec_part.time(static_cast<float>(part.time()));
      rec_part.pid(part.pdgID());
172
      rec_part.status(kTagZDC);
173
174
175
176
177
178
179
180
      rec_part.charge(static_cast<int16_t>(part.charge()));
      rec_part.weight(1.);
      rec_part.direction({mom3s.theta(), mom3s.phi()});
      rec_part.momentum(static_cast<float>(moms));
      rec_part.energy(static_cast<float>(Es));
      rec_part.mass(static_cast<float>(part.mass()));
      rec_part.mcID({part.ID(), m_kMonteCarloSource});
      rc.push_back(rec_part);
181

182
183
184
185
186
187
188
189
190
191
192
193
      if (msgLevel(MSG::DEBUG)) {
        debug()
            << fmt::format(
                   "Found ZDC particle: {}, Etrue: {}, Emeas: {}, ptrue: {}, pmeas: {}, theta_true: {}, theta_meas: {}",
                   part.pdgID(), E, rec_part.energy(), part.ps().mag(), rec_part.p().mag(), th, rec_part.p().theta())
            << endmsg;
      }
    }
    return rc;
  }
  // Fast B0 as in
  // https://github.com/eic/eicsmeardetectors/blob/9a1831dd97bf517b80a06043b9ee4bfb96b483d8/SmearMatrixDetector_0_1_FF.cxx#L254
194
195
  std::vector<RecData> b0(const dd4pod::Geant4ParticleCollection& mc) {
    std::vector<RecData> rc;
196
    for (const auto& part : mc) {
Sylvester Joosten's avatar
Sylvester Joosten committed
197
198
199
200
      if (part.genStatus() > 1) {
        if (msgLevel(MSG::DEBUG)) {
          debug() << "ignoring particle with genStatus = " << part.genStatus() << endmsg;
        }
201
202
203
204
205
206
207
208
        continue;
      }
      // only detect charged hadrons and photons
      if (part.pdgID() != 2212 && part.pdgID() != -2212 && part.pdgID() != 211 && part.pdgID() != -211 &&
          part.pdgID() != 321 && part.pdgID() != -321 && part.pdgID() != 22) {
        continue;
      }
      // only 6-->20 mrad
209
      const auto mom_ion = removeCrossingAngle(part.ps()); //rotateLabToIonDirection(part.ps());
210
      if (mom_ion.theta() < m_thetaMinB0 || mom_ion.theta() > m_thetaMaxB0) {
211
212
        continue;
      }
213
214
215
216
217
218
219
220
      auto rc_part = smearMomentum(part);
      // we don't detect photon energy, just its angles and presence
      if (part.pdgID() == 22) {
        rc_part.p({0,0,0});
        rc_part.energy(0);
      }
      rc_part.status(kTagB0);
      rc.push_back(rc_part);
221
      if (msgLevel(MSG::DEBUG)) {
222
        auto& rec_part = rc.back();
223
224
225
226
227
228
229
230
231
232
233
        debug() << fmt::format("Found B0 particle: {}, ptrue: {}, pmeas: {}, pttrue: {}, ptmeas: {}, theta_true: {}, "
                               "theta_meas: {}",
                               part.pdgID(), part.ps().mag(), rec_part.momentum(), std::hypot(part.ps().x, part.ps().y),
                               std::hypot(rec_part.p().x, rec_part.p().y), part.ps().theta(), rec_part.p().theta())
                << endmsg;
      }
    }

    return rc;
  }

234
235
  std::vector<RecData> rp(const dd4pod::Geant4ParticleCollection& mc) {
    std::vector<RecData> rc;
236
    for (const auto& part : mc) {
Sylvester Joosten's avatar
Sylvester Joosten committed
237
238
239
240
      if (part.genStatus() > 1) {
        if (msgLevel(MSG::DEBUG)) {
          debug() << "ignoring particle with genStatus = " << part.genStatus() << endmsg;
        }
241
242
243
244
245
246
        continue;
      }
      // only detect protons
      if (part.pdgID() != 2212) {
        continue;
      }
247
      const auto mom_ion = removeCrossingAngle(part.ps()); //rotateLabToIonDirection(part.ps());
248
249
      if (mom_ion.theta() < m_thetaMinRP || mom_ion.theta() > m_thetaMaxRP ||
          mom_ion.z < m_pMinRigidityRP * m_ionBeamEnergy) {
250
251
        continue;
      }
252
253
254
      auto rc_part = smearMomentum(part);
      rc_part.status(kTagRP);
      rc.push_back(rc_part);
255
      if (msgLevel(MSG::DEBUG)) {
256
        auto& rec_part = rc.back();
257
258
259
260
261
262
263
264
265
266
        debug() << fmt::format("Found RP particle: {}, ptrue: {}, pmeas: {}, pttrue: {}, ptmeas: {}, theta_true: {}, "
                               "theta_meas: {}",
                               part.pdgID(), part.ps().mag(), rec_part.momentum(), std::hypot(part.ps().x, part.ps().y),
                               std::hypot(rec_part.p().x, rec_part.p().y), part.ps().theta(), rec_part.p().theta())
                << endmsg;
      }
    }
    return rc;
  }

267
268
  std::vector<RecData> omd(const dd4pod::Geant4ParticleCollection& mc) {
    std::vector<RecData> rc;
269
    for (const auto& part : mc) {
Sylvester Joosten's avatar
Sylvester Joosten committed
270
271
272
273
      if (part.genStatus() > 1) {
        if (msgLevel(MSG::DEBUG)) {
          debug() << "ignoring particle with genStatus = " << part.genStatus() << endmsg;
        }
274
275
276
277
278
279
        continue;
      }
      // only detect protons
      if (part.pdgID() != 2212) {
        continue;
      }
280
      const auto mom_ion = removeCrossingAngle(part.ps()); //rotateLabToIonDirection(part.ps());
281
      if (mom_ion.z < m_pMinRigidityOMD * m_ionBeamEnergy || mom_ion.z > m_pMaxRigidityOMD * m_ionBeamEnergy) {
282
283
284
        continue;
      }
      // angle cut
285
286
287
288
289
290
      //const double phi          = (mom_ion.phi() < M_PI) ? mom_ion.phi() : mom_ion.phi() - 2 * M_PI;
      //const bool in_small_angle = (mom_ion.theta() > m_thetaMinFullOMD && mom_ion.theta() < m_thetaMaxFullOMD);
      //const bool in_large_angle = (mom_ion.theta() > m_thetaMinPartialOMD && mom_ion.theta() < m_thetaMaxPartialOMD);
      //if (!in_small_angle || (std::abs(phi) > 1 && !in_large_angle)) {
      //  continue;
      //}
291
292
293
      auto rc_part = smearMomentum(part);
      rc_part.status(kTagOMD);
      rc.push_back(rc_part);
294
      if (msgLevel(MSG::DEBUG)) {
295
        auto& rec_part = rc.back();
296
297
298
299
300
301
302
303
304
305
306
307
        debug() << fmt::format("Found OMD particle: {}, ptrue: {}, pmeas: {}, pttrue: {}, ptmeas: {}, theta_true: {}, "
                               "theta_meas: {}",
                               part.pdgID(), part.ps().mag(), rec_part.momentum(), std::hypot(part.ps().x, part.ps().y),
                               std::hypot(rec_part.p().x, rec_part.p().y), part.ps().theta(), rec_part.p().theta())
                << endmsg;
      }
    }
    return rc;
  }

  // all momentum smearing in EIC-smear for the far-forward region uses
  // the same 2 relations for P and Pt smearing (B0, RP, OMD)
308
  RecData smearMomentum(const dd4pod::ConstGeant4Particle& part) {
309
    const auto mom_ion = rotateLabToIonDirection(part.ps());
310
    const double p     = mom_ion.mag();
311
    const double dp    = (0.025 * p) * m_gaussDist();
312
    const double ps    = p + dp;
313

314
    const double pt  = std::hypot(mom_ion.x, mom_ion.y);
315
    const double dpt = (0.03 * pt) * m_gaussDist();
316
    // just apply relative smearing on px and py
317
318
319
320
321
322
    const double dpxs = (0.03 * mom_ion.x) * m_gaussDist(); //+ (1 + dpt / pt);
    const double dpys = (0.03 * mom_ion.y) * m_gaussDist(); //+ (1 + dpt / pt);

    const double pxs = mom_ion.x + dpxs; 
    const double pys = mom_ion.y + dpys; 

323
    // now get pz
324
    const double pzs = sqrt(ps * ps - pxs * pxs - pys * pys);
325

326
327
    // And build our 3-vector
    const eic::VectorXYZ psmear_ion = {pxs, pys, pzs};
328
    const auto psmear               = rotateIonToLabDirection(psmear_ion);
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
    eic::ReconstructedParticle rec_part;
    rec_part.ID({part.ID(), algorithmID()});
    rec_part.p(psmear);
    rec_part.v({part.vs().x, part.vs().y, part.vs().z});
    rec_part.time(static_cast<float>(part.time()));
    rec_part.pid(part.pdgID());
    rec_part.status(0);
    rec_part.charge(static_cast<int16_t>(part.charge()));
    rec_part.weight(1.);
    rec_part.direction({psmear.theta(), psmear.phi()});
    rec_part.momentum(static_cast<float>(ps));
    rec_part.energy(std::hypot(ps, part.mass()));
    rec_part.mass(static_cast<float>(part.mass()));
    rec_part.mcID({part.ID(), m_kMonteCarloSource});
    return rec_part;
344
345
  }

346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
  // Rotate 25mrad about the y-axis
  eic::VectorXYZ rotateLabToIonDirection(const eic::VectorXYZ& vec) const {
    const double sth = sin(-m_crossingAngle);
    const double cth = cos(-m_crossingAngle);
    return {cth * vec.x + sth * vec.z, vec.y, -sth * vec.x + cth * vec.z};
  }
  eic::VectorXYZ rotateLabToIonDirection(const dd4pod::VectorXYZ& vec) const {
    return rotateLabToIonDirection(eic::VectorXYZ{vec.x, vec.y, vec.z});
  }

  eic::VectorXYZ rotateIonToLabDirection(const eic::VectorXYZ& vec) const {
    const double sth = sin(m_crossingAngle);
    const double cth = cos(m_crossingAngle);
    return {cth * vec.x + sth * vec.z, vec.y, -sth * vec.x + cth * vec.z};
  }
  eic::VectorXYZ rotateIonToLabDirection(const dd4pod::VectorXYZ& vec) const {
    return rotateIonToLabDirection(eic::VectorXYZ{vec.x, vec.y, vec.z});
  }
364
365
366
367
368
369
370
371
372
373
374
375
376
377
  eic::VectorXYZ removeCrossingAngle(const dd4pod::VectorXYZ& vec) const {

    Double_t px = vec.x;
    Double_t py = vec.y;
    Double_t pz = vec.z;
    Double_t s = std::sin(-m_crossingAngle);
    Double_t c = std::cos(-m_crossingAngle);
    Double_t zz = pz;
    Double_t xx = px;
    pz = c*zz - s*xx;
    px = s*zz + c*xx;
    return {px, py, pz};
    //momenta.push_back(ROOT::Math::PxPyPzMVector{px, py, pz, i1.mass});
  }
378

379
}; 
380

381
DECLARE_COMPONENT(SmearedFarForwardParticles)
382

383
} // namespace Jug::Fast
384