TrackingSpatial.cpp 7.4 KB
Newer Older
1
2
#include "TrackingSpatial.h"
#include <TDirectory.h>
3
#include "objects/KDTree.hpp"
4

5
using namespace corryvreckan;
6
using namespace std;
7

8
TrackingSpatial::TrackingSpatial(Configuration config, std::vector<std::shared_ptr<Detector>> detectors)
9
    : Module(std::move(config), std::move(detectors)) {
10
11
12
    spatialCut = m_config.get<double>("spatial_cut", static_cast<double>(Units::convert(200, "um")));
    minHitsOnTrack = m_config.get<size_t>("min_hits_on_track", 6);
    excludeDUT = m_config.get<bool>("exclude_dut", true);
13
14
15
}

/*
16
17

 This algorithm performs the track finding using only spatial information
18
 (no timing). It is based on a linear extrapolation along the z axis, followed
19
 by a nearest neighbour search, and should be well adapted to testbeam
20
 reconstruction with a mostly colinear beam.
21

22
23
 */

24
void TrackingSpatial::initialise() {
Simon Spannagel's avatar
Simon Spannagel committed
25
26

    // Set up histograms
27
28
29
30
31
32
33
34
35
36
37
38
39
40
    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);

    // Loop over all planes
41
    for(auto& detector : get_detectors()) {
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
        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 Y;y_{track}-y [mm];events";
        residualsY[detectorID] = new TH1F("residualsY", title.c_str(), 500, -0.1, 0.1);

        directory->cd();
Simon Spannagel's avatar
Simon Spannagel committed
57
    }
58
59
}

60
StatusCode TrackingSpatial::run(std::shared_ptr<Clipboard> clipboard) {
Simon Spannagel's avatar
Simon Spannagel committed
61
62
63

    // Container for all clusters, and detectors in tracking
    map<string, KDTree*> trees;
64
    vector<std::shared_ptr<Detector>> detectors;
65
    Clusters* referenceClusters = nullptr;
Simon Spannagel's avatar
Simon Spannagel committed
66

67
    // Loop over all detectors and get clusters
Simon Spannagel's avatar
Simon Spannagel committed
68
    double minZ = 1000.;
69
    std::string seedPlane;
70
    for(auto& detector : get_detectors()) {
71
        string detectorID = detector->name();
Simon Spannagel's avatar
Simon Spannagel committed
72
73

        // Get the clusters
74
        Clusters* tempClusters = reinterpret_cast<Clusters*>(clipboard->get(detectorID, "clusters"));
75
        if(tempClusters != nullptr) {
Simon Spannagel's avatar
Simon Spannagel committed
76
            // Store the clusters of the first plane in Z as the reference
77
            if(detector->displacement().Z() < minZ) {
Simon Spannagel's avatar
Simon Spannagel committed
78
                referenceClusters = tempClusters;
79
                seedPlane = detector->name();
80
                minZ = detector->displacement().Z();
Simon Spannagel's avatar
Simon Spannagel committed
81
82
83
84
85
86
87
88
            }
            if(tempClusters->size() == 0)
                continue;

            // Make trees of the clusters on each plane
            KDTree* clusterTree = new KDTree();
            clusterTree->buildSpatialTree(*tempClusters);
            trees[detectorID] = clusterTree;
89
            detectors.push_back(detector);
Simon Spannagel's avatar
Simon Spannagel committed
90
91
            LOG(DEBUG) << "Picked up " << tempClusters->size() << " clusters on device " << detectorID;
        }
92
    }
93

Simon Spannagel's avatar
Simon Spannagel committed
94
    // If there are no detectors then stop trying to track
Simon Spannagel's avatar
Simon Spannagel committed
95
    if(detectors.empty() || referenceClusters == nullptr) {
96
97
98
99
100
101
        // Clean up tree objects
        for(auto tree = trees.cbegin(); tree != trees.cend();) {
            delete tree->second;
            tree = trees.erase(tree);
        }

Simon Spannagel's avatar
Simon Spannagel committed
102
        return StatusCode::NoData;
103
    }
Simon Spannagel's avatar
Simon Spannagel committed
104

105
106
107
    // Output track container
    Tracks* tracks = new Tracks();

Simon Spannagel's avatar
Simon Spannagel committed
108
    // Loop over all clusters
109
    for(auto& cluster : (*referenceClusters)) {
Simon Spannagel's avatar
Simon Spannagel committed
110
111
112
113
114
115
116
117
118
119
120

        // Make a new track
        Track* track = new Track();

        // Add the cluster to the track
        track->addCluster(cluster);

        // Loop over each subsequent planes. For each plane, if extrapolate
        // the hit from the previous plane along the z axis, and look for
        // a neighbour on the new plane. We started on the most upstream
        // plane, so first detector is 1 (not 0)
121
122
        for(auto& detector : detectors) {
            auto detectorID = detector->name();
123
            if(trees.count(detectorID) == 0) {
Simon Spannagel's avatar
Simon Spannagel committed
124
                continue;
125
            }
126
127
            if(detectorID == seedPlane)
                continue;
128
129

            // Check if the DUT should be excluded and obey:
130
            if(excludeDUT && detector->isDUT()) {
131
                continue;
132
            }
Simon Spannagel's avatar
Simon Spannagel committed
133

134
135
136
            // Get the closest neighbour
            LOG(DEBUG) << "- looking for nearest cluster on device " << detectorID;
            Cluster* closestCluster = trees[detectorID]->getClosestNeighbour(cluster);
Simon Spannagel's avatar
Simon Spannagel committed
137

138
            // Check if it is within the spatial window
139
140
141
142
            double distance = sqrt((cluster->global().x() - closestCluster->global().x()) *
                                       (cluster->global().x() - closestCluster->global().x()) +
                                   (cluster->global().y() - closestCluster->global().y()) *
                                       (cluster->global().y() - closestCluster->global().y()));
Simon Spannagel's avatar
Simon Spannagel committed
143

144
            if(distance > spatialCut)
145
                continue;
146

147
148
149
150
            // Add the cluster to the track
            track->addCluster(closestCluster);
            LOG(DEBUG) << "- added cluster to track. Distance is " << distance;
        }
151

152
153
        // Now should have a track with one cluster from each plane
        if(track->nClusters() < minHitsOnTrack) {
154
155
            LOG(DEBUG) << "Not enough clusters on the track, found " << track->nClusters() << " but " << minHitsOnTrack
                       << " required.";
156
157
            delete track;
            continue;
158
        }
159

160
161
162
163
164
165
166
167
        // Fit the track
        track->fit();

        // Save the track
        tracks->push_back(track);

        // Fill histograms
        trackChi2->Fill(track->chi2());
168
        clustersPerTrack->Fill(static_cast<double>(track->nClusters()));
169
        trackChi2ndof->Fill(track->chi2ndof());
170
171
        trackAngleX->Fill(atan(track->direction().X()));
        trackAngleY->Fill(atan(track->direction().Y()));
172
173

        // Make residuals
174
        for(auto& trackCluster : track->clusters()) {
175
            string detectorID = trackCluster->detectorID();
176
177
178
            ROOT::Math::XYZPoint intercept = track->intercept(trackCluster->global().z());
            residualsX[detectorID]->Fill(intercept.X() - trackCluster->global().x());
            residualsY[detectorID]->Fill(intercept.Y() - trackCluster->global().y());
179
        }
180
    }
181

182
    // Save the tracks on the clipboard
183
    tracksPerEvent->Fill(static_cast<double>(tracks->size()));
184
    if(tracks->size() > 0) {
185
        clipboard->put("tracks", reinterpret_cast<Objects*>(tracks));
186
    }
187

188
    // Clean up tree objects
189
190
191
    for(auto tree = trees.cbegin(); tree != trees.cend();) {
        delete tree->second;
        tree = trees.erase(tree);
192
193
    }

194
    LOG(DEBUG) << "End of event";
195
    return StatusCode::Success;
196
}