Belle II Software  release-06-02-00
ECLCalDigitVariables.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 // analysis
10 #include <analysis/VariableManager/Manager.h>
11 
12 // framework
13 #include <framework/core/Module.h>
14 #include <framework/datastore/StoreObjPtr.h>
15 #include <framework/datastore/StoreArray.h>
16 
17 // for crystal geometry
18 #include <ecl/geometry/ECLGeometryPar.h>
19 #include <framework/geometry/B2Vector3.h>
20 
21 // dataobjects
22 #include <mdst/dataobjects/ECLCluster.h>
23 #include <mdst/dataobjects/Track.h>
24 #include <mdst/dataobjects/MCParticle.h>
25 
26 #include <analysis/dataobjects/Particle.h>
27 #include <ecl/dataobjects/ECLCalDigit.h>
28 #include <ecl/dataobjects/ECLShower.h>
29 #include <ecl/dataobjects/ECLCellIdMapping.h>
30 #include <tracking/dataobjects/ExtHit.h>
31 
32 #include <ecl/dataobjects/ECLDsp.h>
33 #include <ecl/geometry/ECLGeometryPar.h>
34 
35 using namespace std;
36 
37 namespace Belle2 {
43  namespace ECLCalDigitVariable {
44 
45  // enum with available data types
46  enum varType {
47  energy = 1,
48  time = 2,
49  timeResolution = 3,
50  twoComponentChi2 = 10,
51  twoComponentTotalEnergy = 11,
52  twoComponentHadronEnergy = 12,
53  twoComponentDiodeEnergy = 13,
54  twoComponentSavedChi2_PhotonHadron = 14,
55  twoComponentSavedChi2_PileUpPhoton = 15,
56  twoComponentSavedChi2_PhotonDiode = 16,
57  twoComponentFitType = 17,
58  weight = 18,
59  phi = 20,
60  theta = 21,
61  phiId = 22,
62  thetaId = 23,
63  cellId = 24,
64  mcenergy = 25,
65  usedforenergy = 26,
66  R_geom = 27, // Requires a geometry environment
67  phiOffset = 30,
68  thetaOffset = 31,
69  phiPointing = 32,
70  thetaPointing = 33
71  };
72 
73  // enum with available center types
74  enum centerType {
75  maxCell = 0,
76  extCell = 1
77  };
78 
79 
81  int getCenterCell(const Particle* particle)
82  {
83  // get maximum cell id for this cluster (ignore weights)
84  const ECLCluster* cluster = particle->getECLCluster();
85  if (cluster) {
86  int maxCellId = -1;
87  double maxEnergy = -1.;
88 
89  auto clusterDigitRelations = cluster->getRelationsTo<ECLCalDigit>();
90  for (unsigned int ir = 0; ir < clusterDigitRelations.size(); ++ir) {
91  const auto calDigit = clusterDigitRelations.object(ir);
92 
93  if (calDigit->getEnergy() > maxEnergy) {
94  maxEnergy = calDigit->getEnergy();
95  maxCellId = calDigit->getCellId();
96  }
97  }
98 
99  return maxCellId;
100  }
101 
102  return -1;
103  }
104 
106  int getExtCell(const Particle* particle)
107  {
108  Const::EDetector myDetID = Const::EDetector::ECL;
109  Const::ChargedStable hypothesis = Const::pion;
110  int pdgCode = abs(hypothesis.getPDGCode());
111 
112  const Track* track = particle->getTrack();
113  if (track) {
114  for (const auto& extHit : track->getRelationsTo<ExtHit>()) {
115  if (abs(extHit.getPdgCode()) != pdgCode) continue;
116  if ((extHit.getDetectorID() != myDetID)) continue;
117  if (extHit.getStatus() != EXT_ENTER) continue;
118 
119  int copyid = extHit.getCopyID();
120 
121  if (copyid == -1) continue;
122  const int cellid = copyid + 1;
123  return cellid;
124  }
125  }
126 
127  return -1;
128  }
129 
131  // This function requires a geometry environment
132  double getCellIdMagnitude(int cellid)
133  {
134  Belle2::ECL::ECLGeometryPar* geom = Belle2::ECL::ECLGeometryPar::Instance();
135  Belle2::B2Vector3D position = geom->GetCrystalPos(cellid - 1);
136  return position.Mag();
137  }
138 
140  double getCalDigitExpertByEnergyRank(const Particle* particle, const std::vector<double>& vars)
141  {
142 
143  if (vars.size() != 2) {
144  B2FATAL("Need exactly two parameters (energy rank, variable id).");
145  }
146 
147  if (int(std::lround(vars[0])) < 0) B2FATAL("Index cannot be negative.");
148 
149  const unsigned int indexIn = int(std::lround(vars[0]));
150 
151  const int varid = int(std::lround(vars[1]));
152 
153  //EnergyToSort vector is used for sorting digits by digit energy measured by FPGAs
154  std::vector<std::tuple<double, unsigned int>> energyToSort;
155 
156  const ECLCluster* cluster = particle->getECLCluster();
157 
158  if (cluster) {
159 
160  auto relatedDigits = cluster->getRelationsTo<ECLCalDigit>();
161 
162  if (indexIn < relatedDigits.size()) {
163 
164  for (unsigned int iRel = 0; iRel < relatedDigits.size(); iRel++) {
165 
166  const auto caldigit = relatedDigits.object(iRel);
167 
168  energyToSort.emplace_back(caldigit->getEnergy(), iRel);
169 
170  }
171 
172  } else {
173  return std::numeric_limits<double>::quiet_NaN();
174  }
175 
176  std::sort(energyToSort.begin(), energyToSort.end(), std::greater<>());
177 
178  const auto [digitEnergy, caldigitIndex] = energyToSort[indexIn];
179 
180  const auto caldigitSelected = relatedDigits.object(caldigitIndex);
181 
182  // Mapping object for phi & theta
183  StoreObjPtr<ECLCellIdMapping> mapping;
184  if (!mapping) {
185  B2ERROR("Mapping not found, did you forget to run the eclFillCellIdMapping module?");
186  return std::numeric_limits<double>::quiet_NaN();
187  }
188 
189  if (varid == varType::energy) {
190  return caldigitSelected->getEnergy();
191  } else if (varid == varType::time) {
192  return caldigitSelected->getTime();
193  } else if (varid == varType::twoComponentChi2) {
194  return caldigitSelected->getTwoComponentChi2();
195  } else if (varid == varType::twoComponentTotalEnergy) {
196  return caldigitSelected->getTwoComponentTotalEnergy();
197  } else if (varid == varType::twoComponentHadronEnergy) {
198  return caldigitSelected->getTwoComponentHadronEnergy();
199  } else if (varid == varType::twoComponentSavedChi2_PhotonHadron) {
200  return caldigitSelected->getTwoComponentSavedChi2(ECLDsp::photonHadron);
201  } else if (varid == varType::twoComponentSavedChi2_PileUpPhoton) {
202  return caldigitSelected->getTwoComponentSavedChi2(ECLDsp::photonHadronBackgroundPhoton);
203  } else if (varid == varType::twoComponentSavedChi2_PhotonDiode) {
204  return caldigitSelected->getTwoComponentSavedChi2(ECLDsp::photonDiodeCrossing);
205  } else if (varid == varType::twoComponentDiodeEnergy) {
206  return caldigitSelected->getTwoComponentDiodeEnergy();
207  } else if (varid == varType::twoComponentFitType) {
208  return int(caldigitSelected->getTwoComponentFitType());
209  } else if (varid == varType::cellId) {
210  return caldigitSelected->getCellId();
211  } else if (varid == varType::weight) {
212  const auto weight = relatedDigits.weight(caldigitIndex);
213  return weight;
214  } else if (varid == varType::phi) {
215  return mapping->getCellIdToPhi(caldigitSelected->getCellId());
216  } else if (varid == varType::theta) {
217  return mapping->getCellIdToTheta(caldigitSelected->getCellId());
218  } else if (varid == varType::R_geom) {
219  return getCellIdMagnitude(caldigitSelected->getCellId());
220  } else {
221  B2FATAL("variable id not found.");
222  }
223 
224  }
225 
226  return std::numeric_limits<double>::quiet_NaN();
227 
228  }
229 
231  double getCalDigitExpert(const Particle* particle, const std::vector<double>& vars)
232  {
233  if (vars.size() != 4) {
234  B2FATAL("Need exactly four parameters (cellid, neighbour area size, variable id, and cluster center (0) or ext (1)).");
235  }
236 
237  StoreObjPtr<ECLCellIdMapping> mapping;
238  const unsigned int posid = int(std::lround(vars[0]));
239  const int nneighbours = int(std::lround(vars[1]));
240  const int varid = int(std::lround(vars[2]));
241  const int extid = int(std::lround(vars[3]));
242 
243  if (!mapping) {
244  B2ERROR("Mapping not found, did you forget to run the eclFillCellIdMapping module?");
245  return std::numeric_limits<double>::quiet_NaN();
246  }
247  if (nneighbours != 5 and nneighbours != 7) {
248  B2FATAL("Please request 5 or 7 neighbour area.");
249  return std::numeric_limits<double>::quiet_NaN();
250  }
251 
252  // get maximum cell id for this cluster (ignore weights)
253  // TODO: if eclshowers exist we could skip that step.
254  int maxCellId = -1;
255  if (extid == centerType::extCell) {
256  maxCellId = getExtCell(particle);
257  } else {
258  maxCellId = getCenterCell(particle);
259  }
260 
261 
262  if (maxCellId < 0) return std::numeric_limits<double>::quiet_NaN();
263 
264  // get the requested neighbourid
265  int neighbourid = -1;
266  std::vector<short int> neighbours;
267 
268  if (nneighbours == 5) {
269  neighbours = mapping->getCellIdToNeighbour5(maxCellId);
270  } else if (nneighbours == 7) {
271  neighbours = mapping->getCellIdToNeighbour7(maxCellId);
272  }
273 
274  if (posid < neighbours.size()) {
275  neighbourid = neighbours[posid];
276  } else {
277  B2WARNING("This position id is not contained in the requested neighbours.");
278  return std::numeric_limits<double>::quiet_NaN();
279  }
280 
281  // some variables can be returned even if no ECLCalDigit is present
282  if (varid == varType::phi) {
283  return mapping->getCellIdToPhi(neighbourid);
284  } else if (varid == varType::theta) {
285  return mapping->getCellIdToTheta(neighbourid);
286  } else if (varid == varType::R_geom) {
287  return getCellIdMagnitude(neighbourid);
288  } else if (varid == varType::phiId) {
289  return mapping->getCellIdToPhiId(neighbourid);
290  } else if (varid == varType::thetaId) {
291  return mapping->getCellIdToThetaId(neighbourid);
292  } else if (varid == varType::cellId) {
293  return neighbourid;
294  }
295  //... and some really need a ECLCalDigit being present
296  else {
297  const int storearraypos = mapping->getCellIdToStoreArray(neighbourid);
298  StoreArray<ECLCalDigit> eclCalDigits;
299 
300  if (storearraypos >= 0) {
301  if (varid == varType::energy) {
302  return eclCalDigits[storearraypos]->getEnergy();
303  } else if (varid == varType::time) {
304  return eclCalDigits[storearraypos]->getTime();
305  } else if (varid == varType::timeResolution) {
306  return eclCalDigits[storearraypos]->getTimeResolution();
307  } else if (varid == varType::twoComponentChi2) {
308  return eclCalDigits[storearraypos]->getTwoComponentChi2();
309  } else if (varid == varType::twoComponentTotalEnergy) {
310  return eclCalDigits[storearraypos]->getTwoComponentTotalEnergy();
311  } else if (varid == varType::twoComponentHadronEnergy) {
312  return eclCalDigits[storearraypos]->getTwoComponentHadronEnergy();
313  } else if (varid == varType::twoComponentSavedChi2_PhotonHadron) {
314  return eclCalDigits[storearraypos]->getTwoComponentSavedChi2(ECLDsp::photonHadron);
315  } else if (varid == varType::twoComponentSavedChi2_PileUpPhoton) {
316  return eclCalDigits[storearraypos]->getTwoComponentSavedChi2(ECLDsp::photonHadronBackgroundPhoton);
317  } else if (varid == varType::twoComponentSavedChi2_PhotonDiode) {
318  return eclCalDigits[storearraypos]->getTwoComponentSavedChi2(ECLDsp::photonDiodeCrossing);
319  } else if (varid == varType::twoComponentDiodeEnergy) {
320  return eclCalDigits[storearraypos]->getTwoComponentDiodeEnergy();
321  } else if (varid == varType::twoComponentFitType) {
322  return int(eclCalDigits[storearraypos]->getTwoComponentFitType());
323  } else if (varid == varType::mcenergy) {
324  // loop over all related MCParticles
325  auto digitMCRelations = eclCalDigits[storearraypos]->getRelationsTo<MCParticle>();
326  double edep = 0.0;
327  for (unsigned int i = 0; i < digitMCRelations.size(); ++i) {
328  edep += digitMCRelations.weight(i);
329  }
330  return edep;
331  } else if (varid == varType::usedforenergy) {
332  const ECLCluster* cluster = particle->getECLCluster();
333  if (cluster) {
334  unsigned int cellid = eclCalDigits[storearraypos]->getCellId();
335 
336  std::vector<unsigned int> listCellIds;
337  auto clusterShowerRelations = cluster->getRelationsWith<ECLShower>(); // should never happen to have more than one shower
338  for (unsigned int ir = 0; ir < clusterShowerRelations.size(); ++ir) {
339  const auto shower = clusterShowerRelations.object(ir);
340  listCellIds = shower->getListOfCrystalsForEnergy();
341  }
342 
343  if (std::find(listCellIds.begin(), listCellIds.end(), cellid) != listCellIds.end()) { //find is faster than count
344  return 1;
345  }
346 
347  return 0;
348  }
349  return std::numeric_limits<double>::quiet_NaN();
350 
351  }
352  } else {
353  return std::numeric_limits<double>::quiet_NaN();
354  }
355  }
356 
357  return std::numeric_limits<double>::quiet_NaN();
358  }
359 
360  double getExtCellExpert(const Particle* particle, int varid, bool front)
361  {
362  ECL::ECLGeometryPar* geometry = ECL::ECLGeometryPar::Instance();
363  if (!geometry) {
364  B2ERROR("Geometry not found!");
365  return std::numeric_limits<double>::quiet_NaN();
366  }
367  const Track* track = particle->getTrack();
368  if (track) {
369  ExtHit* edgeExtHit = nullptr;
370  if (front) {
371  for (const auto& extHit : track->getRelationsTo<ExtHit>()) {
372  if (extHit.getDetectorID() != Const::EDetector::ECL) continue;
373  if (extHit.getStatus() != EXT_ENTER) continue;
374  int crystalID = extHit.getCopyID() - 1;
375  if (crystalID == -1) continue;
376  edgeExtHit = new ExtHit(extHit);
377  break;
378  }
379  } else {
380  auto extHits = track->getRelationsTo<ExtHit>();
381  for (unsigned int iextHit(extHits.size() - 1); iextHit > 0; --iextHit) {
382  const auto extHit = extHits[iextHit];
383  if (extHit->getDetectorID() != Const::EDetector::ECL) continue;
384  if (extHit->getStatus() != EXT_EXIT) continue;
385  int crystalID = extHit->getCopyID() - 1;
386  if (crystalID == -1) break;
387  edgeExtHit = new ExtHit(*extHit);
388  break;
389  }
390  }
391 
392  if (!edgeExtHit) return std::numeric_limits<double>::quiet_NaN();
393  const TVector3& extHitPosition = edgeExtHit->getPosition();
394  const TVector3& trackPointing = edgeExtHit->getMomentum();
395 
396  geometry->Mapping(edgeExtHit->getCopyID() - 1);
397  const int thetaID = geometry->GetThetaID();
398  const int phiID = geometry->GetPhiID();
399 
400  const TVector3& crystalCenterPosition = geometry->GetCrystalPos(geometry->GetCellID(thetaID, phiID));
401  const TVector3& crystalOrientation = geometry->GetCrystalVec(geometry->GetCellID(thetaID, phiID));
402  const TVector3& crystalPositionOnSurface = crystalCenterPosition - (crystalCenterPosition - extHitPosition).Dot(
403  crystalOrientation.Unit()) * crystalOrientation.Unit();
404 
405  if (varid == varType::phiOffset) {
406  return extHitPosition.DeltaPhi(crystalPositionOnSurface);
407  } else if (varid == varType::thetaOffset) {
408  return extHitPosition.Theta() - crystalPositionOnSurface.Theta();
409  } else if (varid == varType::phiPointing) {
410  return trackPointing.DeltaPhi(crystalOrientation);
411  } else if (varid == varType::thetaPointing) {
412  return trackPointing.Theta() - crystalOrientation.Theta();
413  }
414  }
415 
416  return std::numeric_limits<double>::quiet_NaN();
417  }
418  }
419 
420  namespace Variable {
421 
423  double getECLCalDigitEnergyByEnergyRank(const Particle* particle, const std::vector<double>& vars)
424  {
425  if (vars.size() != 1) {
426  B2FATAL("Need exactly one parameters (energy index).");
427  }
428  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::energy};
429  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
430  }
431 
433  double getECLCalDigitTimeByEnergyRank(const Particle* particle, const std::vector<double>& vars)
434  {
435  if (vars.size() != 1) {
436  B2FATAL("Need exactly one parameters (energy index).");
437  }
438  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::time};
439  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
440  }
441 
443  double getCellIdByEnergyRank(const Particle* particle, const std::vector<double>& vars)
444  {
445  if (vars.size() != 1) {
446  B2FATAL("Need exactly one parameters (energy index).");
447  }
448  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::cellId};
449  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
450  }
451 
453  double getTwoComponentFitTypeByEnergyRank(const Particle* particle, const std::vector<double>& vars)
454  {
455  if (vars.size() != 1) {
456  B2FATAL("Need exactly one parameters (energy index).");
457  }
458  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::twoComponentFitType};
459  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
460  }
461 
463  double getTwoComponentChi2ByEnergyRank(const Particle* particle, const std::vector<double>& vars)
464  {
465  if (vars.size() != 1) {
466  B2FATAL("Need exactly one parameters (energy index).");
467  }
468  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::twoComponentChi2};
469  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
470  }
471 
473  double getTwoComponentTotalEnergyByEnergyRank(const Particle* particle, const std::vector<double>& vars)
474  {
475  if (vars.size() != 1) {
476  B2FATAL("Need exactly one parameters (energy index).");
477  }
478  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::twoComponentTotalEnergy};
479  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
480  }
481 
483  double getTwoComponentHadronEnergyByEnergyRank(const Particle* particle, const std::vector<double>& vars)
484  {
485  if (vars.size() != 1) {
486  B2FATAL("Need exactly one parameters (energy index).");
487  }
488  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::twoComponentHadronEnergy};
489  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
490  }
491 
493  double getTwoComponentDiodeEnergyByEnergyRank(const Particle* particle, const std::vector<double>& vars)
494  {
495  if (vars.size() != 1) {
496  B2FATAL("Need exactly one parameters (energy index).");
497  }
498  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::twoComponentDiodeEnergy};
499  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
500  }
501 
503  double getTwoComponentChi2SavedByEnergyRank_PhotonHadron(const Particle* particle, const std::vector<double>& vars)
504  {
505  if (vars.size() != 1) {
506  B2FATAL("Need exactly one parameters (energy index).");
507  }
508  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::twoComponentSavedChi2_PhotonHadron};
509  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
510  }
511 
513  double getTwoComponentChi2SavedByEnergyRank_PileUpPhoton(const Particle* particle, const std::vector<double>& vars)
514  {
515  if (vars.size() != 1) {
516  B2FATAL("Need exactly one parameters (energy index).");
517  }
518  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::twoComponentSavedChi2_PileUpPhoton};
519  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
520  }
521 
523  double getTwoComponentChi2SavedByEnergyRank_PhotonDiode(const Particle* particle, const std::vector<double>& vars)
524  {
525  if (vars.size() != 1) {
526  B2FATAL("Need exactly one parameters (energy index).");
527  }
528  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::twoComponentSavedChi2_PhotonDiode};
529  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
530  }
531 
533  double getWeightByEnergyRank(const Particle* particle, const std::vector<double>& vars)
534  {
535  if (vars.size() != 1) {
536  B2FATAL("Need exactly one parameters (energy index).");
537  }
538  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::weight};
539  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
540  }
541 
543  double getECLCalDigitEnergy(const Particle* particle, const std::vector<double>& vars)
544  {
545  if (vars.size() != 2) {
546  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
547  }
548 
549  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::energy, ECLCalDigitVariable::centerType::maxCell};
550  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
551  }
552 
554  double getMCEnergy(const Particle* particle, const std::vector<double>& vars)
555  {
556  if (vars.size() != 2) {
557  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
558  }
559 
560  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::mcenergy, ECLCalDigitVariable::centerType::maxCell};
561  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
562  }
563 
564 
566  double getExtECLCalDigitEnergy(const Particle* particle, const std::vector<double>& vars)
567  {
568  if (vars.size() != 2) {
569  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
570  }
571 
572  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::energy, ECLCalDigitVariable::centerType::extCell};
573  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
574  }
575 
576 
578  double getECLCalDigitTime(const Particle* particle, const std::vector<double>& vars)
579  {
580  if (vars.size() != 2) {
581  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
582  }
583 
584  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::time, ECLCalDigitVariable::centerType::maxCell};
585  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
586  }
587 
589  double getExtECLCalDigitTime(const Particle* particle, const std::vector<double>& vars)
590  {
591  if (vars.size() != 2) {
592  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
593  }
594 
595  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::time, ECLCalDigitVariable::centerType::extCell};
596  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
597  }
598 
600  double getExtECLCalDigitTimeResolution(const Particle* particle, const std::vector<double>& vars)
601  {
602  if (vars.size() != 2) {
603  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
604  }
605 
606  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::timeResolution, ECLCalDigitVariable::centerType::extCell};
607  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
608  }
609 
611  double getECLCalDigitTimeResolution(const Particle* particle, const std::vector<double>& vars)
612  {
613  if (vars.size() != 2) {
614  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
615  }
616 
617  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::timeResolution, ECLCalDigitVariable::centerType::maxCell};
618  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
619  }
620 
622  double getTwoComponentChi2(const Particle* particle, const std::vector<double>& vars)
623  {
624  if (vars.size() != 2) {
625  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
626  }
627 
628  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentChi2, ECLCalDigitVariable::centerType::maxCell};
629  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
630  }
631 
633  double getExtECLCalDigitTwoComponentChi2(const Particle* particle, const std::vector<double>& vars)
634  {
635  if (vars.size() != 2) {
636  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
637  }
638 
639  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentChi2, ECLCalDigitVariable::centerType::extCell};
640  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
641  }
642 
644  double getTwoComponentTotalEnergy(const Particle* particle, const std::vector<double>& vars)
645  {
646  if (vars.size() != 2) {
647  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
648  }
649 
650  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentTotalEnergy, ECLCalDigitVariable::centerType::maxCell};
651  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
652  }
653 
655  double getExtECLCalDigitTwoComponentTotalEnergy(const Particle* particle, const std::vector<double>& vars)
656  {
657  if (vars.size() != 2) {
658  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
659  }
660 
661  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentTotalEnergy, ECLCalDigitVariable::centerType::extCell};
662  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
663  }
664 
666  double getTwoComponentHadronEnergy(const Particle* particle, const std::vector<double>& vars)
667  {
668  if (vars.size() != 2) {
669  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
670  }
671 
672  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentHadronEnergy, ECLCalDigitVariable::centerType::maxCell};
673  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
674  }
675 
677  double getUsedForClusterEnergy(const Particle* particle, const std::vector<double>& vars)
678  {
679  if (vars.size() != 2) {
680  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
681  }
682 
683  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::usedforenergy, ECLCalDigitVariable::centerType::maxCell};
684  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
685  }
686 
688  double getExtECLCalDigitTwoComponentHadronEnergy(const Particle* particle, const std::vector<double>& vars)
689  {
690  if (vars.size() != 2) {
691  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
692  }
693 
694  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentHadronEnergy, ECLCalDigitVariable::centerType::extCell};
695  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
696  }
697 
699  double getECLCalDigitTwoComponentDiodeEnergy(const Particle* particle, const std::vector<double>& vars)
700  {
701  if (vars.size() != 2) {
702  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
703  }
704 
705  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentDiodeEnergy, ECLCalDigitVariable::centerType::maxCell};
706  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
707  }
708 
710  double getExtECLCalDigitTwoComponentDiodeEnergy(const Particle* particle, const std::vector<double>& vars)
711  {
712  if (vars.size() != 2) {
713  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
714  }
715 
716  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentDiodeEnergy, ECLCalDigitVariable::centerType::extCell};
717  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
718  }
719 
721  double getECLCalDigitTwoComponentFitType(const Particle* particle, const std::vector<double>& vars)
722  {
723  if (vars.size() != 2) {
724  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
725  }
726 
727  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentFitType, ECLCalDigitVariable::centerType::maxCell};
728  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
729  }
730 
732  double getExtECLCalDigitTwoComponentFitType(const Particle* particle, const std::vector<double>& vars)
733  {
734  if (vars.size() != 2) {
735  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
736  }
737 
738  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentFitType, ECLCalDigitVariable::centerType::extCell};
739  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
740  }
741 
743  double getECLCalDigitTwoComponentChi2Saved_PhotonHadron(const Particle* particle, const std::vector<double>& vars)
744  {
745  if (vars.size() != 2) {
746  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
747  }
748 
749  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentSavedChi2_PhotonHadron, ECLCalDigitVariable::centerType::maxCell};
750  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
751  }
752 
754  double getExtECLCalDigitTwoComponentChi2Saved_PhotonHadron(const Particle* particle, const std::vector<double>& vars)
755  {
756  if (vars.size() != 2) {
757  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
758  }
759 
760  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentSavedChi2_PhotonHadron, ECLCalDigitVariable::centerType::extCell};
761  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
762  }
763 
765  double getECLCalDigitTwoComponentChi2Saved_PileUpPhoton(const Particle* particle, const std::vector<double>& vars)
766  {
767  if (vars.size() != 2) {
768  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
769  }
770 
771  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentSavedChi2_PileUpPhoton, ECLCalDigitVariable::centerType::maxCell};
772  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
773  }
774 
776  double getExtECLCalDigitTwoComponentChi2Saved_PileUpPhoton(const Particle* particle, const std::vector<double>& vars)
777  {
778  if (vars.size() != 2) {
779  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
780  }
781 
782  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentSavedChi2_PileUpPhoton, ECLCalDigitVariable::centerType::extCell};
783  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
784  }
785 
787  double getECLCalDigitTwoComponentChi2Saved_PhotonDiode(const Particle* particle, const std::vector<double>& vars)
788  {
789  if (vars.size() != 2) {
790  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
791  }
792 
793  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentSavedChi2_PhotonDiode, ECLCalDigitVariable::centerType::maxCell};
794  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
795  }
796 
798  double getExtECLCalDigitTwoComponentChi2Saved_PhotonDiode(const Particle* particle, const std::vector<double>& vars)
799  {
800  if (vars.size() != 2) {
801  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
802  }
803 
804  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::twoComponentSavedChi2_PhotonDiode, ECLCalDigitVariable::centerType::extCell};
805  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
806  }
807 
809  double getPhi(const Particle* particle, const std::vector<double>& vars)
810  {
811  if (vars.size() != 2) {
812  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
813  }
814 
815  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::phi, ECLCalDigitVariable::centerType::maxCell};
816  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
817  }
818 
820  double getExtPhi(const Particle* particle, const std::vector<double>& vars)
821  {
822  if (vars.size() != 2) {
823  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
824  }
825 
826  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::phi, ECLCalDigitVariable::centerType::extCell};
827  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
828  }
829 
831  double getTheta(const Particle* particle, const std::vector<double>& vars)
832  {
833  if (vars.size() != 2) {
834  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
835  }
836 
837  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::theta, ECLCalDigitVariable::centerType::maxCell};
838  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
839  }
840 
842  double getExtTheta(const Particle* particle, const std::vector<double>& vars)
843  {
844  if (vars.size() != 2) {
845  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
846  }
847 
848  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::theta, ECLCalDigitVariable::centerType::extCell};
849  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
850  }
851 
853  double getPhiId(const Particle* particle, const std::vector<double>& vars)
854  {
855  if (vars.size() != 2) {
856  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
857  }
858 
859  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::phiId, ECLCalDigitVariable::centerType::maxCell};
860  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
861  }
862 
864  double getExtPhiId(const Particle* particle, const std::vector<double>& vars)
865  {
866  if (vars.size() != 2) {
867  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
868  }
869 
870  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::phiId, ECLCalDigitVariable::centerType::extCell};
871  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
872  }
873 
875  double getThetaId(const Particle* particle, const std::vector<double>& vars)
876  {
877  if (vars.size() != 2) {
878  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
879  }
880 
881  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::thetaId, ECLCalDigitVariable::centerType::maxCell};
882  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
883  }
884 
886  double getExtThetaId(const Particle* particle, const std::vector<double>& vars)
887  {
888  if (vars.size() != 2) {
889  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
890  }
891 
892  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::thetaId, ECLCalDigitVariable::centerType::extCell};
893  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
894  }
895 
897  double getCellId(const Particle* particle, const std::vector<double>& vars)
898  {
899  if (vars.size() != 2) {
900  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
901  }
902 
903  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::cellId, ECLCalDigitVariable::centerType::maxCell};
904  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
905  }
906 
908  double getCenterCellId(const Particle* particle)
909  {
910  const int centercellid = ECLCalDigitVariable::getCenterCell(particle);
911 
912  if (centercellid < 0) return std::numeric_limits<double>::quiet_NaN();
913  return centercellid;
914  }
915 
917  double getCenterCellCrystalTheta(const Particle* particle)
918  {
919  const int centercellid = ECLCalDigitVariable::getCenterCell(particle);
920  StoreObjPtr<ECLCellIdMapping> mapping;
921 
922  if (!mapping) {
923  B2ERROR("Mapping not found, did you forget to run the eclFillCellIdMapping module?");
924  return std::numeric_limits<double>::quiet_NaN();
925  }
926 
927  if (centercellid < 0) return std::numeric_limits<double>::quiet_NaN();
928  return mapping->getCellIdToTheta(centercellid);
929  }
930 
932  double getCenterCellCrystalPhi(const Particle* particle)
933  {
934  const int centercellid = ECLCalDigitVariable::getCenterCell(particle);
935  StoreObjPtr<ECLCellIdMapping> mapping;
936 
937  if (!mapping) {
938  B2ERROR("Mapping not found, did you forget to run the eclFillCellIdMapping module?");
939  return std::numeric_limits<double>::quiet_NaN();
940  }
941 
942  if (centercellid < 0) return std::numeric_limits<double>::quiet_NaN();
943  return mapping->getCellIdToPhi(centercellid);
944  }
945 
947  double getExtCellId(const Particle* particle)
948  {
949  const int extcellid = ECLCalDigitVariable::getExtCell(particle);
950 
951  if (extcellid < 0) return std::numeric_limits<double>::quiet_NaN();
952  return extcellid;
953  }
954 
956  double getCenterCellThetaId(const Particle* particle)
957  {
958  const int centercellid = ECLCalDigitVariable::getCenterCell(particle);
959  StoreObjPtr<ECLCellIdMapping> mapping;
960 
961  if (!mapping) {
962  B2ERROR("Mapping not found, did you forget to run the eclFillCellIdMapping module?");
963  return std::numeric_limits<double>::quiet_NaN();
964  }
965 
966  if (centercellid < 0) return std::numeric_limits<double>::quiet_NaN();
967  return mapping->getCellIdToThetaId(centercellid);
968  }
969 
971  double getCenterCellPhiId(const Particle* particle)
972  {
973  const int centercellid = ECLCalDigitVariable::getCenterCell(particle);
974  StoreObjPtr<ECLCellIdMapping> mapping;
975 
976  if (!mapping) {
977  B2ERROR("Mapping not found, did you forget to run the eclFillCellIdMapping module?");
978  return std::numeric_limits<double>::quiet_NaN();
979  }
980 
981  if (centercellid < 0) return std::numeric_limits<double>::quiet_NaN();
982  return mapping->getCellIdToPhiId(centercellid);
983  }
984 
986  double getExtCellThetaId(const Particle* particle)
987  {
988  const int extcellid = ECLCalDigitVariable::getExtCell(particle);
989  StoreObjPtr<ECLCellIdMapping> mapping;
990 
991  if (!mapping) {
992  B2ERROR("Mapping not found, did you forget to run the eclFillCellIdMapping module?");
993  return std::numeric_limits<double>::quiet_NaN();
994  }
995 
996  if (extcellid < 0) return std::numeric_limits<double>::quiet_NaN();
997  return mapping->getCellIdToThetaId(extcellid);
998  }
999 
1001  double getExtCellPhiId(const Particle* particle)
1002  {
1003  const int extcellid = ECLCalDigitVariable::getExtCell(particle);
1004  StoreObjPtr<ECLCellIdMapping> mapping;
1005 
1006  if (!mapping) {
1007  B2ERROR("Mapping not found, did you forget to run the eclFillCellIdMapping module?");
1008  return std::numeric_limits<double>::quiet_NaN();
1009  }
1010 
1011  if (extcellid < 0) return std::numeric_limits<double>::quiet_NaN();
1012  return mapping->getCellIdToPhiId(extcellid);
1013  }
1014 
1016  double getExtCellCrystalTheta(const Particle* particle)
1017  {
1018  const int extcellid = ECLCalDigitVariable::getExtCell(particle);
1019  StoreObjPtr<ECLCellIdMapping> mapping;
1020 
1021  if (!mapping) {
1022  B2ERROR("Mapping not found, did you forget to run the eclFillCellIdMapping module?");
1023  return std::numeric_limits<double>::quiet_NaN();
1024  }
1025 
1026  if (extcellid < 0) return std::numeric_limits<double>::quiet_NaN();
1027  return mapping->getCellIdToTheta(extcellid);
1028  }
1029 
1031  double getExtCellCrystalPhi(const Particle* particle)
1032  {
1033  const int extcellid = ECLCalDigitVariable::getExtCell(particle);
1034  StoreObjPtr<ECLCellIdMapping> mapping;
1035 
1036  if (!mapping) {
1037  B2ERROR("Mapping not found, did you forget to run the eclFillCellIdMapping module?");
1038  return std::numeric_limits<double>::quiet_NaN();
1039  }
1040 
1041  if (extcellid < 0) return std::numeric_limits<double>::quiet_NaN();
1042  return mapping->getCellIdToPhi(extcellid);
1043  }
1044 
1046  double getCenterCellIndex(const Particle* particle, const std::vector<double>& vars)
1047  {
1048  if (vars.size() != 1) {
1049  B2FATAL("Need exactly one parameters (neighbour area size).");
1050  }
1051 
1052  StoreObjPtr<ECLCellIdMapping> mapping;
1053  const int nneighbours = int(std::lround(vars[0]));
1054 
1055  if (!mapping) {
1056  B2ERROR("Mapping not found, did you forget to run the eclFillCellIdMapping module?");
1057  return std::numeric_limits<double>::quiet_NaN();
1058  }
1059  if (nneighbours != 5 and nneighbours != 7) {
1060  B2FATAL("Please request 5 or 7 neighbour area.");
1061  return std::numeric_limits<double>::quiet_NaN();
1062  }
1063 
1064  const int centercellid = ECLCalDigitVariable::getCenterCell(particle);
1065  if (centercellid < 0) return std::numeric_limits<double>::quiet_NaN();
1066 
1067  std::vector<short int> neighbours;
1068 
1069  if (nneighbours == 5) {
1070  neighbours = mapping->getCellIdToNeighbour5(centercellid);
1071  } else if (nneighbours == 7) {
1072  neighbours = mapping->getCellIdToNeighbour7(centercellid);
1073  }
1074 
1075  for (unsigned int idx = 0; idx < neighbours.size(); idx++) {
1076  if (neighbours[idx] == centercellid) return idx;
1077  }
1078 
1079  return std::numeric_limits<double>::quiet_NaN();
1080  }
1081 
1082 
1084  double getExtCenterCellIndex(const Particle* particle, const std::vector<double>& vars)
1085  {
1086  if (vars.size() != 1) {
1087  B2FATAL("Need exactly one parameters (neighbour area size).");
1088  }
1089 
1090  StoreObjPtr<ECLCellIdMapping> mapping;
1091  const int nneighbours = int(std::lround(vars[0]));
1092 
1093  if (!mapping) {
1094  B2ERROR("Mapping not found, did you forget to run the eclFillCellIdMapping module?");
1095  return std::numeric_limits<double>::quiet_NaN();
1096  }
1097  if (nneighbours != 5 and nneighbours != 7) {
1098  B2FATAL("Please request 5 or 7 neighbour area.");
1099  return std::numeric_limits<double>::quiet_NaN();
1100  }
1101 
1102  const int centercellid = ECLCalDigitVariable::getExtCell(particle);
1103  if (centercellid < 0) return std::numeric_limits<double>::quiet_NaN();
1104 
1105  std::vector<short int> neighbours;
1106 
1107  if (nneighbours == 5) {
1108  neighbours = mapping->getCellIdToNeighbour5(centercellid);
1109  } else if (nneighbours == 7) {
1110  neighbours = mapping->getCellIdToNeighbour7(centercellid);
1111  }
1112 
1113  for (unsigned int idx = 0; idx < neighbours.size(); idx++) {
1114  if (neighbours[idx] == centercellid) return idx;
1115  }
1116 
1117  return std::numeric_limits<double>::quiet_NaN();
1118  }
1119 
1120  double getTotalECLCalDigitMCEnergy(const Particle* particle)
1121  {
1122  const MCParticle* mcparticle = particle->getRelatedTo<MCParticle>();
1123  if (mcparticle == nullptr)
1124  return std::numeric_limits<double>::quiet_NaN();
1125 
1126  double sum = 0.0;
1127  auto mcDigitRelations = mcparticle->getRelationsWith<ECLCalDigit>();
1128  for (unsigned int ir = 0; ir < mcDigitRelations.size(); ++ir) {
1129  sum += mcDigitRelations.weight(ir);
1130  }
1131 
1132  return sum;
1133 
1134  }
1135 
1136  double getClusterECLCalDigitMCEnergy(const Particle* particle)
1137  {
1138  // get MCParticle (return if there is none)
1139  const MCParticle* mcparticle = particle->getRelatedTo<MCParticle>();
1140  if (mcparticle == nullptr)
1141  return std::numeric_limits<double>::quiet_NaN();
1142 
1143  const ECLCluster* cluster = particle->getECLCluster();
1144  if (cluster) {
1145  std::vector<unsigned int> listCellIds;
1146  auto clusterShowerRelations = cluster->getRelationsWith<ECLShower>(); // should never happen to have more than one shower
1147  for (unsigned int ir = 0; ir < clusterShowerRelations.size(); ++ir) {
1148  const auto shower = clusterShowerRelations.object(ir);
1149  listCellIds = shower->getListOfCrystalsForEnergy();
1150  }
1151 
1152  double sum = 0.0;
1153  auto clusterDigitRelations = mcparticle->getRelationsWith<ECLCalDigit>();
1154  for (unsigned int ir = 0; ir < clusterDigitRelations.size(); ++ir) {
1155 
1156  // check if this digit is in the current cluster
1157  unsigned int cellid = clusterDigitRelations.object(ir)->getCellId();
1158  if (std::find(listCellIds.begin(), listCellIds.end(), cellid) != listCellIds.end()) { //find is faster than count
1159  sum += clusterDigitRelations.weight(ir);
1160  }
1161  }
1162 
1163  return sum;
1164  }
1165 
1166  return std::numeric_limits<float>::quiet_NaN();
1167 
1168  }
1169 
1170  double getExtFrontPositionPhiOffset(const Particle* particle)
1171  {
1172  return ECLCalDigitVariable::getExtCellExpert(particle, ECLCalDigitVariable::varType::phiOffset, true);
1173  }
1174 
1175  double getExtFrontPositionThetaOffset(const Particle* particle)
1176  {
1177  return ECLCalDigitVariable::getExtCellExpert(particle, ECLCalDigitVariable::varType::thetaOffset, true);
1178  }
1179 
1180  double getExtFrontPositionPhiPointing(const Particle* particle)
1181  {
1182  return ECLCalDigitVariable::getExtCellExpert(particle, ECLCalDigitVariable::varType::phiPointing, true);
1183  }
1184 
1185  double getExtFrontPositionThetaPointing(const Particle* particle)
1186  {
1187  return ECLCalDigitVariable::getExtCellExpert(particle, ECLCalDigitVariable::varType::thetaPointing, true);
1188  }
1189 
1190  double getExtBackPositionPhiOffset(const Particle* particle)
1191  {
1192  return ECLCalDigitVariable::getExtCellExpert(particle, ECLCalDigitVariable::varType::phiOffset, false);
1193  }
1194 
1195  double getExtBackPositionThetaOffset(const Particle* particle)
1196  {
1197  return ECLCalDigitVariable::getExtCellExpert(particle, ECLCalDigitVariable::varType::thetaOffset, false);
1198  }
1199 
1200  double getExtBackPositionPhiPointing(const Particle* particle)
1201  {
1202  return ECLCalDigitVariable::getExtCellExpert(particle, ECLCalDigitVariable::varType::phiPointing, false);
1203  }
1204 
1205  double getExtBackPositionThetaPointing(const Particle* particle)
1206  {
1207  return ECLCalDigitVariable::getExtCellExpert(particle, ECLCalDigitVariable::varType::thetaPointing, false);
1208  }
1209 
1211  double getPhiByEnergyRank(const Particle* particle, const std::vector<double>& vars)
1212  {
1213  if (vars.size() != 1) {
1214  B2FATAL("Need exactly one parameters (energy index).");
1215  }
1216  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::phi};
1217  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
1218  }
1219 
1221  double getThetaByEnergyRank(const Particle* particle, const std::vector<double>& vars)
1222  {
1223  if (vars.size() != 1) {
1224  B2FATAL("Need exactly one parameters (energy index).");
1225  }
1226  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::theta};
1227  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
1228  }
1229 
1231  double getR(const Particle* particle, const std::vector<double>& vars)
1232  {
1233  if (vars.size() != 2) {
1234  B2FATAL("Need exactly two parameters (cellid and neighbour area size).");
1235  }
1236 
1237  std::vector<double> parameters {vars[0], vars[1], ECLCalDigitVariable::varType::R_geom, ECLCalDigitVariable::centerType::maxCell};
1238  return ECLCalDigitVariable::getCalDigitExpert(particle, parameters);
1239  }
1240 
1242  double getRByEnergyRank(const Particle* particle, const std::vector<double>& vars)
1243  {
1244  if (vars.size() != 1) {
1245  B2FATAL("Need exactly one parameters (energy index).");
1246  }
1247  std::vector<double> parameters {vars[0], ECLCalDigitVariable::varType::R_geom};
1248  return ECLCalDigitVariable::getCalDigitExpertByEnergyRank(particle, parameters);
1249  }
1250 
1251  double getClusterNHitsThreshold(const Particle* particle, const std::vector<double>& vars)
1252  {
1253  if (vars.size() != 1) {
1254  B2FATAL("Need exactly one parameter (energy threshold in GeV).");
1255  }
1256  const double threshold = vars[0];
1257 
1258  const ECLCluster* cluster = particle->getECLCluster();
1259  if (cluster) {
1260  double nhits = 0.;
1261 
1262  auto clusterDigitRelations = cluster->getRelationsTo<ECLCalDigit>();
1263  for (unsigned int ir = 0; ir < clusterDigitRelations.size(); ++ir) {
1264  const auto calDigit = clusterDigitRelations.object(ir);
1265  const auto weight = clusterDigitRelations.weight(ir);
1266 
1267  // take the unweighted eclcaldigit energy for this check (closer to real hardware threshold)
1268  if (calDigit->getEnergy() > threshold) {
1269  nhits += weight;
1270  }
1271  }
1272 
1273  return nhits;
1274  }
1275  return std::numeric_limits<float>::quiet_NaN();
1276  }
1277 
1278 
1279  VARIABLE_GROUP("ECL Calibration (cDST)");
1280  REGISTER_VARIABLE("eclcaldigitEnergy(i, j)", getECLCalDigitEnergy,
1281  "[calibration] Returns the energy of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1282  REGISTER_VARIABLE("eclcaldigitTime(i, j)", getECLCalDigitTime,
1283  "[calibration] Returns the time of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1284  REGISTER_VARIABLE("eclcaldigitTimeResolution(i, j)", getECLCalDigitTimeResolution,
1285  "[calibration] Returns the time resolution (dt99) of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1286  REGISTER_VARIABLE("eclcaldigitTwoComponentChi2(i, j)", getTwoComponentChi2,
1287  "[calibration] Returns the two component fit chi2 of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1288  REGISTER_VARIABLE("eclcaldigitTwoComponentTotalEnergy(i, j)", getTwoComponentTotalEnergy,
1289  "[calibration] Returns the two component total energy of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1290  REGISTER_VARIABLE("eclcaldigitTwoComponentHadronEnergy(i, j)", getTwoComponentHadronEnergy,
1291  "[calibration] Returns the two component hadron energy of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1292  REGISTER_VARIABLE("eclcaldigitPhi(i, j)", getPhi,
1293  "[calibration] Returns phi of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1294  REGISTER_VARIABLE("eclcaldigitTheta(i, j)", getTheta,
1295  "[calibration] Returns theta of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1296  REGISTER_VARIABLE("eclcaldigitR(i, j)", getR,
1297  "Returns R (from a geometry object) of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1298  REGISTER_VARIABLE("eclcaldigitPhiId(i, j)", getPhiId,
1299  "[calibration] Returns the phi Id of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1300  REGISTER_VARIABLE("eclcaldigitThetaId(i, j)", getThetaId,
1301  "[calibration] Returns the theta Id of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1302  REGISTER_VARIABLE("eclcaldigitCellId(i, j)", getCellId,
1303  "[calibration] Returns the cell id of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours (1-based)");
1304  REGISTER_VARIABLE("eclcaldigitUsedForClusterEnergy(i, j)", getUsedForClusterEnergy,
1305  " [calibration] Returns the 0 (not used) 1 (used) of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours (1-based)");
1306 
1307  REGISTER_VARIABLE("eclcaldigitCenterCellId", getCenterCellId, "[calibration] Returns the center cell id");
1308  REGISTER_VARIABLE("eclcaldigitCenterCellThetaId", getCenterCellThetaId, "[calibration] Returns the center cell theta id");
1309  REGISTER_VARIABLE("eclcaldigitCenterCellPhiId", getCenterCellPhiId, "[calibration] Returns the center cell phi id");
1310  REGISTER_VARIABLE("eclcaldigitCenterCellCrystalTheta", getCenterCellCrystalTheta,
1311  "[calibration] Returns the center cell crystal theta");
1312  REGISTER_VARIABLE("eclcaldigitCenterCellCrystalPhi", getCenterCellCrystalPhi,
1313  "[calibration] Returns the center cell crystal phi");
1314  REGISTER_VARIABLE("eclcaldigitCenterCellIndex(i)", getCenterCellIndex,
1315  "[calibration] Returns the center cell index (within its 5x5 (i=5) or 7x7 (i=7) neighbours)");
1316  REGISTER_VARIABLE("eclcaldigitMCEnergy(i, j)", getMCEnergy,
1317  "[calibration] Returns the true deposited energy of all particles of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours (1-based)");
1318  REGISTER_VARIABLE("clusterNHitsThreshold(i)", getClusterNHitsThreshold,
1319  "[calibration] Returns sum of crystal weights sum(w_i) with w_i<=1 associated to this cluster above threshold (in GeV)");
1320 
1321  VARIABLE_GROUP("ECL Calibration (based on extrapolated tracks) (cDST)");
1322  REGISTER_VARIABLE("eclcaldigitExtEnergy(i, j)", getExtECLCalDigitEnergy,
1323  "[calibration] Returns the energy of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an ext track");
1324  REGISTER_VARIABLE("eclcaldigitExtTime(i, j)", getExtECLCalDigitTime,
1325  "[calibration] Returns the time of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an ext track");
1326  REGISTER_VARIABLE("eclcaldigitExtTimeResolution(i, j)", getExtECLCalDigitTimeResolution,
1327  "[calibration] Returns the time resolution (dt99) of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an ext track");
1328  REGISTER_VARIABLE("eclcaldigitExtTwoComponentTotalEnergy(i, j)", getExtECLCalDigitTwoComponentTotalEnergy,
1329  "[calibration] Returns the TwoComponentTotalEnergy of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an ext track");
1330  REGISTER_VARIABLE("eclcaldigitExtTwoComponentHadronEnergy(i, j)", getExtECLCalDigitTwoComponentHadronEnergy,
1331  "[calibration] Returns the TwoComponentHadronEnergy of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an ext track");
1332  REGISTER_VARIABLE("eclcaldigitExtTwoComponentChi2(i, j)", getExtECLCalDigitTwoComponentChi2,
1333  "[calibration] Returns the TwoComponentchi2 of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an ext track");
1334  REGISTER_VARIABLE("eclcaldigitExtPhi(i, j)", getExtPhi,
1335  "[calibration] Returns phi of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an extrapolated track");
1336  REGISTER_VARIABLE("eclcaldigitExtTheta(i, j)", getExtTheta,
1337  "[calibration] Returns theta of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an extrapolated track");
1338  REGISTER_VARIABLE("eclcaldigitExtPhiId(i, j)", getExtPhiId,
1339  "[calibration] Returns the phi Id of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an extrapolated track");
1340  REGISTER_VARIABLE("eclcaldigitExtThetaId(i, j)", getExtThetaId,
1341  "[calibration] Returns the theta Id of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an extrapolated track");
1342  REGISTER_VARIABLE("eclcaldigitExtCellId", getExtCellId, "[calibration] Returns the extrapolated cell id");
1343  REGISTER_VARIABLE("eclcaldigitExtCellThetaId", getExtCellThetaId, "[calibration] Returns the ext cell theta id");
1344  REGISTER_VARIABLE("eclcaldigitExtCellPhiId", getExtCellPhiId, "[calibration] Returns the ext cell phi id");
1345  REGISTER_VARIABLE("eclcaldigitExtCellCrystalTheta", getExtCellCrystalTheta, "[calibration] Returns the ext cell crystal theta");
1346  REGISTER_VARIABLE("eclcaldigitExtCellCrystalPhi", getExtCellCrystalPhi, "[calibration] Returns the ext cell crystal phi");
1347  REGISTER_VARIABLE("eclcaldigitExtCenterCellIndex(i)", getExtCenterCellIndex,
1348  "[calibration] Returns the center cell index (within its 5x5 (i=5) or 7x7 (i=7) neighbours) for an ext track");
1349 
1350  REGISTER_VARIABLE("eclcaldigitExtFrontPositionPhiOffset", getExtFrontPositionPhiOffset,
1351  "[calibration] Returns the difference in the azimuthal angle (in radians)"
1352  "between the position where the track hit the front face of the ECL and the"
1353  "center of the struck crystal projected onto the front surface.");
1354  REGISTER_VARIABLE("eclcaldigitExtFrontPositionThetaOffset", getExtFrontPositionThetaOffset,
1355  "[calibration] Returns the difference in the polar angle (in radians)"
1356  "between the position where the track hit the front face of the ECL and the"
1357  "center of the struck crystal projected onto the front surface.");
1358  REGISTER_VARIABLE("eclcaldigitExtFrontPositionPhiPointing", getExtFrontPositionPhiPointing,
1359  "[calibration] Returns the difference in the azimuthal angle (in radians)"
1360  "between the momentum direction when the track hit the front face of the ECL and the"
1361  "orientation of the struck crystal.");
1362  REGISTER_VARIABLE("eclcaldigitExtFrontPositionThetaPointing", getExtFrontPositionThetaPointing,
1363  "[calibration] Returns the difference in the polar angle (in radians)"
1364  "between the momentum direction when the track hit the front face of the ECL and the"
1365  "orientation of the struck crystal.");
1366 
1367  REGISTER_VARIABLE("eclcaldigitExtTwoComponentFitType(i, j)", getExtECLCalDigitTwoComponentFitType,
1368  "[calibration] Returns the TwoComponentFitType of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an ext track");
1369  REGISTER_VARIABLE("eclcaldigitExtTwoComponentDiodeEnergy(i, j)", getExtECLCalDigitTwoComponentDiodeEnergy,
1370  "[calibration] Returns the TwoComponentDiodeEnergy of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an ext track");
1371  REGISTER_VARIABLE("eclcaldigitExtTwoComponentChi2Saved_PhotonHadron(i, j)", getExtECLCalDigitTwoComponentChi2Saved_PhotonHadron,
1372  "[calibration] Returns the TwoComponentChi2Saved_PhotonHadron of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an ext track");
1373  REGISTER_VARIABLE("eclcaldigitExtTwoComponentChi2Saved_PileUpPhoton(i, j)", getExtECLCalDigitTwoComponentChi2Saved_PileUpPhoton,
1374  "[calibration] Returns the TwoComponentChi2Saved_PileUpPhoton of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an ext track");
1375  REGISTER_VARIABLE("eclcaldigitExtTwoComponentChi2Saved_PhotonDiode(i, j)", getExtECLCalDigitTwoComponentChi2Saved_PhotonDiode,
1376  "[calibration] Returns the TwoComponentChi2Saved_PhotonDiode of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours for an ext track");
1377 
1378  REGISTER_VARIABLE("eclcaldigitTwoComponentFitType(i, j)", getECLCalDigitTwoComponentFitType,
1379  "[calibration] Returns the TwoComponentFitType of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1380  REGISTER_VARIABLE("eclcaldigitTwoComponentDiodeEnergy(i, j)", getECLCalDigitTwoComponentDiodeEnergy,
1381  "[calibration] Returns the TwoComponentDiodeEnergy of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1382  REGISTER_VARIABLE("eclcaldigitTwoComponentChi2Saved_PhotonHadron(i, j)", getECLCalDigitTwoComponentChi2Saved_PhotonHadron,
1383  "[calibration] Returns the TwoComponentChi2Saved_PhotonHadron of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1384  REGISTER_VARIABLE("eclcaldigitTwoComponentChi2Saved_PileUpPhoton(i, j)", getECLCalDigitTwoComponentChi2Saved_PileUpPhoton,
1385  "[calibration] Returns the TwoComponentChi2Saved_PileUpPhoton of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1386  REGISTER_VARIABLE("eclcaldigitTwoComponentChi2Saved_PhotonDiode(i, j)", getECLCalDigitTwoComponentChi2Saved_PhotonDiode,
1387  "[calibration] Returns the TwoComponentChi2Saved_PhotonDiode of the i-th caldigit for 5x5 (j=5) or 7x7 (j=7) neighbours");
1388 
1389  REGISTER_VARIABLE("eclcaldigitEnergyByEnergyRank(i)", getECLCalDigitEnergyByEnergyRank,
1390  "[calibration] Returns the caldigit energy of the i-th highest energy caldigit in the cluster (i>=0)");
1391 
1392  REGISTER_VARIABLE("eclcaldigitTimeByEnergyRank(i)", getECLCalDigitTimeByEnergyRank,
1393  "[calibration] Returns the caldigit time of the i-th highest energy caldigit in the cluster (i>=0)");
1394 
1395  REGISTER_VARIABLE("eclcaldigitTwoComponentFitTypeByEnergyRank(i)", getTwoComponentFitTypeByEnergyRank,
1396  "[calibration] Returns the offline fit type of the i-th highest energy caldigit in the cluster (i>=0)");
1397 
1398  REGISTER_VARIABLE("eclcaldigitTwoComponentChi2ByEnergyRank(i)", getTwoComponentChi2ByEnergyRank,
1399  "[calibration] Returns the two component chi2 of the i-th highest energy caldigit in the cluster (i>=0)");
1400 
1401  REGISTER_VARIABLE("eclcaldigitTwoComponentEnergyByEnergyRank(i)", getTwoComponentTotalEnergyByEnergyRank,
1402  "[calibration] Returns the two component total energy of the i-th highest energy caldigit in the cluster (i>=0)");
1403 
1404  REGISTER_VARIABLE("eclcaldigitTwoComponentHadronEnergyByEnergyRank(i)", getTwoComponentHadronEnergyByEnergyRank,
1405  "[calibration] Returns the two component fit Hadron Energy of the i-th highest energy caldigit in the cluster (i>=0)");
1406 
1407  REGISTER_VARIABLE("eclcaldigitTwoComponentDiodeEnergyByEnergyRank(i)", getTwoComponentDiodeEnergyByEnergyRank,
1408  "[calibration] Returns the two component fit Diode Energy of the i-th highest energy caldigit in the cluster (i>=0)");
1409 
1410  REGISTER_VARIABLE("eclcaldigitTwoComponentChi2SavedByEnergyRank_PhotonHadron(i)", getTwoComponentChi2SavedByEnergyRank_PhotonHadron,
1411  "[calibration] Returns the chi2 for the photo+hadron fit type of the i-th highest energy caldigit in the cluster (i>=0)");
1412 
1413  REGISTER_VARIABLE("eclcaldigitTwoComponentChi2SavedByEnergyRank_PileUpPhoton(i)", getTwoComponentChi2SavedByEnergyRank_PileUpPhoton,
1414  "[calibration] Returns the chi2 for the photo+hadron+pile-up photon fit type of the i-th highest energy caldigit in the cluster (i>=0)");
1415 
1416  REGISTER_VARIABLE("eclcaldigitTwoComponentChi2SavedByEnergyRank_PhotonDiode(i)", getTwoComponentChi2SavedByEnergyRank_PhotonDiode,
1417  "[calibration] Returns the chi2 for the photo+diode fit type of the i-th highest energy caldigit in the cluster (i>=0)");
1418 
1419  REGISTER_VARIABLE("eclcaldigitWeightByEnergyRank(i)", getWeightByEnergyRank,
1420  "[calibration] Returns the weight of the i-th highest energy caldigit in the cluster (i>=0)");
1421 
1422  REGISTER_VARIABLE("eclcaldigitCellIdByEnergyRank(i)", getCellIdByEnergyRank,
1423  "[calibration] Returns the cell id of the i-th highest energy caldigit in the cluster (i>=0)");
1424 
1425  REGISTER_VARIABLE("totalECLCalDigitMCEnergy", getTotalECLCalDigitMCEnergy,
1426  "[calibration] Returns total deposited MC energy in all ECLCalDigits for the MC particle");
1427 
1428  REGISTER_VARIABLE("clusterECLCalDigitMCEnergy", getClusterECLCalDigitMCEnergy,
1429  "[calibration] Returns total deposited MC energy in all ECLCalDigits for the MC particle that are used to calculate the cluster energy");
1430 
1431  REGISTER_VARIABLE("eclcaldigitPhiByEnergyRank(i)", getPhiByEnergyRank,
1432  "Returns phi of the i-th highest energy caldigit in the cluster (i>=0)");
1433 
1434  REGISTER_VARIABLE("eclcaldigitThetaByEnergyRank(i)", getThetaByEnergyRank,
1435  "Returns theta of the i-th highest energy caldigit in the cluster (i>=0)");
1436 
1437  REGISTER_VARIABLE("eclcaldigitRByEnergyRank(i)", getRByEnergyRank,
1438  "Returns R of the i-th highest energy caldigit in the cluster (i>=0)");
1439  }
1440 
1441  // Create an empty module which allows basf2 to easily find the library and load it from the steering file
1442  class EnableECLCalDigitVariablesModule: public Module {}; // Register this module to create a .map lookup file.
1443  REG_MODULE(EnableECLCalDigitVariables);
1445 }
Belle2::B2Vector3::Mag
DataType Mag() const
The magnitude (rho in spherical coordinate system).
Definition: B2Vector3.h:144
Belle2::ECL::ECLGeometryPar
The Class for ECL Geometry Parameters.
Definition: ECLGeometryPar.h:39
Belle2::ECL::ECLGeometryPar::Instance
static ECLGeometryPar * Instance()
Static method to get a reference to the ECLGeometryPar instance.
Definition: ECLGeometryPar.cc:149
Belle2::Module
Base class for Modules.
Definition: Module.h:72
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17