Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
/*
* Reconstruct the cluster with Center of Gravity method
* Logarithmic weighting is used for mimicing energy deposit in transverse direction
*
* Author: Chao Peng (ANL), 09/27/2020
*/
#include <algorithm>
#include "Gaudi/Property.h"
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiKernel/ToolHandle.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
// FCCSW
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
// Event Model related classes
#include "eicd/ClusterCollection.h"
using namespace Gaudi::Units;
namespace Jug::Reco {
class ClusterRecoCoG : public GaudiAlgorithm
{
public:
Gaudi::Property<double> m_logWeightThres{this, "logWeightThres", 4.2};
DataHandle<eic::ClusterCollection>
m_clusterCollection{"clusterCollection", Gaudi::DataHandle::Reader, this};
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
ClusterRecoCoG(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
{
declareProperty("clusterCollection", m_clusterCollection, "");
}
StatusCode initialize() override
{
if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
return StatusCode::SUCCESS;
}
StatusCode execute() override
{
// input collections
const auto &clusters = *m_clusterCollection.get();
// reconstruct hit position for the cluster
for (auto &cl : clusters) {
reconstruct(cl);
info() << cl.energy()/GeV << " GeV, (" << cl.position()[0]/mm << ", "
<< cl.position()[1]/mm << ", " << cl.position()[2]/mm << ")" << endmsg;
}
return StatusCode::SUCCESS;
}
private:
void reconstruct(eic::Cluster cl) {
float totalE = 0.;
for (auto &hit : cl.hits()) {
totalE += hit.energy();
}
cl.energy(totalE);
// center of gravity with logarithmic weighting
float totalW = 0., x = 0., y = 0., z = 0.;
for (auto &hit : cl.hits()) {
float weight = m_logWeightThres + std::log(hit.energy()/totalE);
totalW += weight;
x += hit.position().x * weight;
y += hit.position().y * weight;
z += hit.position().z * weight;
}
cl.position() = std::array<float, 3>{x/totalW, y/totalW, z/totalW};
}
};
DECLARE_COMPONENT(ClusterRecoCoG)
} // namespace Jug::Reco