ClusterRecoCoG.cpp 4.03 KB
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
/*
 *  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:
34
    Gaudi::Property<double> m_logWeightBase{this, "logWeightBase", 3.6};
35
36
    Gaudi::Property<double> m_depthCorrection{this, "depthCorrection", 0.0};
    Gaudi::Property<std::string> m_moduleDimZName{this, "moduleDimZName", ""};
37
38
    DataHandle<eic::ClusterCollection>
        m_clusterCollection{"clusterCollection", Gaudi::DataHandle::Reader, this};
39
40
    // Pointer to the geometry service
    SmartIF<IGeoSvc> m_geoSvc;
41
42
43
44
45
46
47
48
49
50
51
52
53

    // 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;
        }
54
55
56
57
58
59
        m_geoSvc = service("GeoSvc");
        if (!m_geoSvc) {
            error() << "Unable to locate Geometry Service. "
                    << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
            return StatusCode::FAILURE;
        }
60
61
62
63
        // update depth correction if a name is provided
        if (!m_moduleDimZName.value().empty()) {
            m_depthCorrection = m_geoSvc->detector()->constantAsDouble(m_moduleDimZName);
        }
Jihee Kim's avatar
Jihee Kim committed
64
        //info() << "z_length " << depth << endmsg;
65
66
67
68
69
70
        return StatusCode::SUCCESS;
    }

    StatusCode execute() override
    {
        // input collections
Chao Peng's avatar
Chao Peng committed
71
        auto &clusters = *m_clusterCollection.get();
72
73
74
        // reconstruct hit position for the cluster
        for (auto &cl : clusters) {
            reconstruct(cl);
Chao Peng's avatar
Chao Peng committed
75
76
            // info() << cl.energy()/GeV << " GeV, (" << cl.position().x/mm << ", "
            //        << cl.position().y/mm << ", " << cl.position().z/mm << ")" << endmsg;
77
78
79
80
81
82
        }

        return StatusCode::SUCCESS;
    }

private:
83
84
85
86
87
88
89
90
91
    void reconstruct(eic::Cluster cl)
    {
        // no hits
        if (cl.hits_size() == 0) {
            return;
        }

        // calculate total energy, find the cell with the maximum energy deposit
        float totalE = 0., maxE = 0.;
Chao Peng's avatar
Chao Peng committed
92
        auto centerID = cl.hits_begin()->cellID();
93
        for (auto &hit : cl.hits()) {
94
95
96
97
            auto energy = hit.energy();
            totalE += energy;
            if (energy > maxE) {
                maxE = energy;
Chao Peng's avatar
Chao Peng committed
98
                centerID = hit.cellID();
99
            }
100
101
102
103
        }
        cl.energy(totalE);

        // center of gravity with logarithmic weighting
104
        float tw = 0., x = 0., y = 0., z = 0.;
105
        for (auto &hit : cl.hits()) {
106
            // suppress low energy contributions
107
108
            float w = std::max(0., m_logWeightBase + std::log(hit.energy()/totalE));
            tw += w;
Chao Peng's avatar
Chao Peng committed
109
110
111
            x += hit.local_x() * w;
            y += hit.local_y() * w;
            z += hit.local_z() * w;
112
        }
113
114
115
116

        // convert local position to global position, use the cell with max edep as a reference
        auto volman = m_geoSvc->detector()->volumeManager();
        auto alignment = volman.lookupDetector(centerID).nominal();
117
        auto gpos = alignment.localToWorld(dd4hep::Position(x/tw, y/tw, z/tw + m_depthCorrection));
118
119

        cl.position({gpos.x(), gpos.y(), gpos.z()});
120
121
122
123
124
125
126
    }
};

DECLARE_COMPONENT(ClusterRecoCoG)

} // namespace Jug::Reco