Tracking4D.cpp 10.9 KB
Newer Older
1
#include "Tracking4D.h"
2
#include <TCanvas.h>
Simon Spannagel's avatar
Simon Spannagel committed
3
#include <TDirectory.h>
4
#include "objects/KDTree.hpp"
5

6
using namespace corryvreckan;
7
using namespace std;
8

9
Tracking4D::Tracking4D(Configuration config, std::vector<std::shared_ptr<Detector>> detectors)
10
    : Module(std::move(config), std::move(detectors)) {
11

Simon Spannagel's avatar
Simon Spannagel committed
12
    // Default values for cuts
13
14
15
    timingCut = m_config.get<double>("timingCut", static_cast<double>(Units::convert(200, "ns")));
    spatialCut = m_config.get<double>("spatialCut", static_cast<double>(Units::convert(0.2, "mm")));
    minHitsOnTrack = m_config.get<size_t>("minHitsOnTrack", 6);
16
    excludeDUT = m_config.get<bool>("excludeDUT", true);
17
18
}

19
void Tracking4D::initialise() {
Simon Spannagel's avatar
Simon Spannagel committed
20
21

    // Set up histograms
22
23
24
25
26
27
28
29
30
31
32
33
    std::string title = "Track #chi^{2};#chi^{2};events";
    trackChi2 = new TH1F("trackChi2", title.c_str(), 150, 0, 150);
    title = "Track #chi^{2}/ndof;#chi^{2}/ndof;events";
    trackChi2ndof = new TH1F("trackChi2ndof", title.c_str(), 100, 0, 50);
    title = "Clusters per track;clusters;tracks";
    clustersPerTrack = new TH1F("clustersPerTrack", title.c_str(), 10, 0, 10);
    title = "Track multiplicity;tracks;events";
    tracksPerEvent = new TH1F("tracksPerEvent", title.c_str(), 100, 0, 100);
    title = "Track angle X;angle_{x} [rad];events";
    trackAngleX = new TH1F("trackAngleX", title.c_str(), 2000, -0.01, 0.01);
    title = "Track angle Y;angle_{y} [rad];events";
    trackAngleY = new TH1F("trackAngleY", title.c_str(), 2000, -0.01, 0.01);
Simon Spannagel's avatar
Simon Spannagel committed
34

35
    // Loop over all planes
36
    for(auto& detector : get_detectors()) {
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
        auto detectorID = detector->name();

        TDirectory* directory = getROOTDirectory();
        TDirectory* local_directory = directory->mkdir(detectorID.c_str());
        if(local_directory == nullptr) {
            throw RuntimeError("Cannot create or access local ROOT directory for module " + this->getUniqueName());
        }
        local_directory->cd();

        title = detectorID + " Residual X;x_{track}-x [mm];events";
        residualsX[detectorID] = new TH1F("residualsX", title.c_str(), 500, -0.1, 0.1);
        title = detectorID + " Residual X, size 1;x_{track}-x [mm];events";
        residualsXwidth1[detectorID] = new TH1F("residualsXwidth1", title.c_str(), 500, -0.1, 0.1);
        title = detectorID + " Residual X, size 2;x_{track}-x [mm];events";
        residualsXwidth2[detectorID] = new TH1F("residualsXwidth2", title.c_str(), 500, -0.1, 0.1);
        title = detectorID + " Residual X, size 3;x_{track}-x [mm];events";
        residualsXwidth3[detectorID] = new TH1F("residualsXwidth3", title.c_str(), 500, -0.1, 0.1);
        title = detectorID + " Residual Y;y_{track}-y [mm];events";
        residualsY[detectorID] = new TH1F("residualsY", title.c_str(), 500, -0.1, 0.1);
        title = detectorID + " Residual Y, size 1;y_{track}-y [mm];events";
        residualsYwidth1[detectorID] = new TH1F("residualsYwidth1", title.c_str(), 500, -0.1, 0.1);
        title = detectorID + " Residual Y, size 2;y_{track}-y [mm];events";
        residualsYwidth2[detectorID] = new TH1F("residualsYwidth2", title.c_str(), 500, -0.1, 0.1);
        title = detectorID + " Residual Y, size 3;y_{track}-y [mm];events";
        residualsYwidth3[detectorID] = new TH1F("residualsYwidth3", title.c_str(), 500, -0.1, 0.1);
62

63
        directory->cd();
Simon Spannagel's avatar
Simon Spannagel committed
64
    }
65
66
}

67
StatusCode Tracking4D::run(std::shared_ptr<Clipboard> clipboard) {
Simon Spannagel's avatar
Simon Spannagel committed
68
69
70
71
72

    LOG(DEBUG) << "Start of event";
    // Container for all clusters, and detectors in tracking
    map<string, KDTree*> trees;
    vector<string> detectors;
73
    Clusters* referenceClusters = nullptr;
Simon Spannagel's avatar
Simon Spannagel committed
74

75
    // Loop over all planes and get clusters
Simon Spannagel's avatar
Simon Spannagel committed
76
    bool firstDetector = true;
77
    std::string seedPlane;
78
    for(auto& detector : get_detectors()) {
79
        string detectorID = detector->name();
Simon Spannagel's avatar
Simon Spannagel committed
80
81

        // Get the clusters
82
83
        Clusters* tempClusters = reinterpret_cast<Clusters*>(clipboard->get(detectorID, "clusters"));
        if(tempClusters == nullptr || tempClusters->size() == 0) {
Simon Spannagel's avatar
Simon Spannagel committed
84
85
86
87
            LOG(DEBUG) << "Detector " << detectorID << " does not have any clusters on the clipboard";
        } else {
            // Store them
            LOG(DEBUG) << "Picked up " << tempClusters->size() << " clusters from " << detectorID;
88
            if(firstDetector && !detector->isDUT()) {
Simon Spannagel's avatar
Simon Spannagel committed
89
                referenceClusters = tempClusters;
90
91
                seedPlane = detector->name();
                LOG(DEBUG) << "Seed plane is " << seedPlane;
92
                firstDetector = false;
Simon Spannagel's avatar
Simon Spannagel committed
93
            }
94

Simon Spannagel's avatar
Simon Spannagel committed
95
96
97
98
99
            KDTree* clusterTree = new KDTree();
            clusterTree->buildTimeTree(*tempClusters);
            trees[detectorID] = clusterTree;
            detectors.push_back(detectorID);
        }
100
    }
101

Simon Spannagel's avatar
Simon Spannagel committed
102
    // If there are no detectors then stop trying to track
103
    if(detectors.size() == 0 || referenceClusters == nullptr) {
Simon Spannagel's avatar
Simon Spannagel committed
104
105
106
107
108
109
        // Clean up tree objects
        for(auto tree = trees.cbegin(); tree != trees.cend();) {
            delete tree->second;
            tree = trees.erase(tree);
        }

110
        return StatusCode::Success;
111
112
113
114
    }

    // Output track container
    Tracks* tracks = new Tracks();
Simon Spannagel's avatar
Simon Spannagel committed
115
116
117

    // Loop over all clusters
    map<Cluster*, bool> used;
118
    for(auto& cluster : (*referenceClusters)) {
Simon Spannagel's avatar
Simon Spannagel committed
119
120

        // Make a new track
Simon Spannagel's avatar
Simon Spannagel committed
121
        LOG(DEBUG) << "Looking at next seed cluster";
122

Simon Spannagel's avatar
Simon Spannagel committed
123
124
125
126
127
128
        Track* track = new Track();
        // Add the cluster to the track
        track->addCluster(cluster);
        track->setTimestamp(cluster->timestamp());
        used[cluster] = true;

129
        // Loop over each subsequent plane and look for a cluster within the timing cuts
130
        for(auto& detectorID : detectors) {
131
132
            // Get the detector
            auto det = get_detector(detectorID);
133

134
            // Check if the DUT should be excluded and obey:
135
            if(excludeDUT && det->isDUT()) {
136
                LOG(DEBUG) << "Skipping DUT plane.";
Simon Spannagel's avatar
Simon Spannagel committed
137
                continue;
138
            }
139

140
            if(detectorID == seedPlane)
Simon Spannagel's avatar
Simon Spannagel committed
141
                continue;
142
            if(trees.count(detectorID) == 0)
Simon Spannagel's avatar
Simon Spannagel committed
143
144
                continue;

145
            // Get all neighbours within the timing cut
146
            LOG(DEBUG) << "Searching for neighbouring cluster on " << detectorID;
147
            LOG(DEBUG) << "- cluster time is " << Units::display(cluster->timestamp(), {"ns", "us", "s"});
148
            Cluster* closestCluster = nullptr;
Simon Spannagel's avatar
Simon Spannagel committed
149
            double closestClusterDistance = spatialCut;
150
            Clusters neighbours = trees[detectorID]->getAllClustersInTimeWindow(cluster, timingCut);
Simon Spannagel's avatar
Simon Spannagel committed
151
152
153
154
155
156
157

            LOG(DEBUG) << "- found " << neighbours.size() << " neighbours";

            // Now look for the spatially closest cluster on the next plane
            double interceptX, interceptY;
            if(track->nClusters() > 1) {
                track->fit();
158

159
                PositionVector3D<Cartesian3D<double>> interceptPoint = det->getIntercept(track);
Simon Spannagel's avatar
Simon Spannagel committed
160
161
162
                interceptX = interceptPoint.X();
                interceptY = interceptPoint.Y();
            } else {
163
164
                interceptX = cluster->global().x();
                interceptY = cluster->global().y();
Simon Spannagel's avatar
Simon Spannagel committed
165
166
167
            }

            // Loop over each neighbour in time
168
            for(size_t ne = 0; ne < neighbours.size(); ne++) {
Simon Spannagel's avatar
Simon Spannagel committed
169
170
171
                Cluster* newCluster = neighbours[ne];

                // Calculate the distance to the previous plane's cluster/intercept
172
173
                double distance = sqrt((interceptX - newCluster->global().x()) * (interceptX - newCluster->global().x()) +
                                       (interceptY - newCluster->global().y()) * (interceptY - newCluster->global().y()));
Simon Spannagel's avatar
Simon Spannagel committed
174
175
176
177
178
179
180
181

                // If this is the closest keep it
                if(distance < closestClusterDistance) {
                    closestClusterDistance = distance;
                    closestCluster = newCluster;
                }
            }

182
            if(closestCluster == nullptr) {
183
                LOG(DEBUG) << "No cluster within spatial cut.";
Simon Spannagel's avatar
Simon Spannagel committed
184
                continue;
185
            }
Simon Spannagel's avatar
Simon Spannagel committed
186
187
188
189
190
191
192
193

            // Add the cluster to the track
            LOG(DEBUG) << "- added cluster to track";
            track->addCluster(closestCluster);
        } //*/

        // Now should have a track with one cluster from each plane
        if(track->nClusters() < minHitsOnTrack) {
Simon Spannagel's avatar
Simon Spannagel committed
194
195
            LOG(DEBUG) << "Not enough clusters on the track, found " << track->nClusters() << " but " << minHitsOnTrack
                       << " required.";
Simon Spannagel's avatar
Simon Spannagel committed
196
197
            delete track;
            continue;
198
        }
199

Simon Spannagel's avatar
Simon Spannagel committed
200
201
202
203
204
205
        // Fit the track and save it
        track->fit();
        tracks->push_back(track);

        // Fill histograms
        trackChi2->Fill(track->chi2());
206
        clustersPerTrack->Fill(static_cast<double>(track->nClusters()));
Simon Spannagel's avatar
Simon Spannagel committed
207
        trackChi2ndof->Fill(track->chi2ndof());
208
209
        trackAngleX->Fill(atan(track->direction().X()));
        trackAngleY->Fill(atan(track->direction().Y()));
Simon Spannagel's avatar
Simon Spannagel committed
210
211
212

        // Make residuals
        Clusters trackClusters = track->clusters();
213
        for(auto& trackCluster : trackClusters) {
Simon Spannagel's avatar
Simon Spannagel committed
214
            string detectorID = trackCluster->detectorID();
215
216
            ROOT::Math::XYZPoint intercept = track->intercept(trackCluster->global().z());
            residualsX[detectorID]->Fill(intercept.X() - trackCluster->global().x());
217
            if(trackCluster->columnWidth() == 1)
218
                residualsXwidth1[detectorID]->Fill(intercept.X() - trackCluster->global().x());
219
            if(trackCluster->columnWidth() == 2)
220
                residualsXwidth2[detectorID]->Fill(intercept.X() - trackCluster->global().x());
221
            if(trackCluster->columnWidth() == 3)
222
223
                residualsXwidth3[detectorID]->Fill(intercept.X() - trackCluster->global().x());
            residualsY[detectorID]->Fill(intercept.Y() - trackCluster->global().y());
224
            if(trackCluster->rowWidth() == 1)
225
                residualsYwidth1[detectorID]->Fill(intercept.Y() - trackCluster->global().y());
226
            if(trackCluster->rowWidth() == 2)
227
                residualsYwidth2[detectorID]->Fill(intercept.Y() - trackCluster->global().y());
228
            if(trackCluster->rowWidth() == 3)
229
                residualsYwidth3[detectorID]->Fill(intercept.Y() - trackCluster->global().y());
Simon Spannagel's avatar
Simon Spannagel committed
230
        }
231
232
233
234

        // Improve the track timestamp by taking the average of all planes
        double avg_track_time = 0;
        for(auto& trackCluster : trackClusters) {
Simon Spannagel's avatar
Simon Spannagel committed
235
            avg_track_time += static_cast<double>(Units::convert(trackCluster->timestamp(), "ns"));
236
            avg_track_time -= static_cast<double>(Units::convert(trackCluster->global().z(), "mm") / (299.792458));
237
        }
Simon Spannagel's avatar
Simon Spannagel committed
238
        track->setTimestamp(avg_track_time / static_cast<double>(track->nClusters()));
239
    }
240

Simon Spannagel's avatar
Simon Spannagel committed
241
242
    // Save the tracks on the clipboard
    if(tracks->size() > 0) {
243
        clipboard->put("tracks", reinterpret_cast<Objects*>(tracks));
244
        tracksPerEvent->Fill(static_cast<double>(tracks->size()));
Jens Kroeger's avatar
Jens Kroeger committed
245
246
    } else {
        delete tracks;
247
    }
248

Simon Spannagel's avatar
Simon Spannagel committed
249
    // Clean up tree objects
Simon Spannagel's avatar
Simon Spannagel committed
250
251
252
    for(auto tree = trees.cbegin(); tree != trees.cend();) {
        delete tree->second;
        tree = trees.erase(tree);
Simon Spannagel's avatar
Simon Spannagel committed
253
    }
254

Simon Spannagel's avatar
Simon Spannagel committed
255
    LOG(DEBUG) << "End of event";
256
    return StatusCode::Success;
257
}