AnalysisTelescope.cpp 7.02 KB
Newer Older
1
#include "AnalysisTelescope.h"
2
#include "objects/Cluster.hpp"
3
#include "objects/MCParticle.hpp"
4

5
6
#include <TDirectory.h>

7
8
9
10
11
12
using namespace corryvreckan;
using namespace std;

AnalysisTelescope::AnalysisTelescope(Configuration config, std::vector<std::shared_ptr<Detector>> detectors)
    : Module(std::move(config), std::move(detectors)) {

13
    chi2ndofCut = m_config.get<double>("chi2ndof_cut", 3.);
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
}

void AnalysisTelescope::initialise() {

    // Initialise biased telescope residuals per telescope device and telescope resolution plots (using MCparticles) at the
    // position of the DUTs
    for(auto& detector : get_detectors()) {
        TDirectory* directory = getROOTDirectory();
        TDirectory* local_directory = directory->mkdir(detector->name().c_str());
        if(local_directory == nullptr) {
            throw RuntimeError("Cannot create or access local ROOT directory for module " + this->getUniqueName());
        }
        local_directory->cd();

        if(detector->isDUT()) {
            std::string title = detector->name() + " Telescope resolution X;x_{track}-x_{MC} [mm];events";
            telescopeResolutionX[detector->name()] = new TH1F("telescopeResolutionX", title.c_str(), 600, -0.2, 0.2);
            title = detector->name() + " Telescope resolution Y;y_{track}-y_{MC} [mm];events";
            telescopeResolutionY[detector->name()] = new TH1F("telescopeResolutionY", title.c_str(), 600, -0.2, 0.2);
        } else {
            std::string title = detector->name() + " Biased residual X (local);x_{track}-x_{cluster} [mm];events";
            telescopeResidualsLocalX[detector->name()] = new TH1F("residualX_local", title.c_str(), 400, -0.2, 0.2);
            title = detector->name() + " Biased residual Y (local);y_{track}-y_{cluster} [mm];events";
            telescopeResidualsLocalY[detector->name()] = new TH1F("residualY_local", title.c_str(), 400, -0.2, 0.2);
            title = detector->name() + " Biased residual X (global);x_{track}-x_{cluster} [mm];events";
            telescopeResidualsX[detector->name()] = new TH1F("residualX_global", title.c_str(), 400, -0.2, 0.2);
            title = detector->name() + " Biased residual Y (global);y_{track}-y_{cluster} [mm];events";
            telescopeResidualsY[detector->name()] = new TH1F("residualY_global", title.c_str(), 400, -0.2, 0.2);

            title = detector->name() + " Biased MC residual X (local);x_{track}-x_{MC} [mm];events";
            telescopeMCresidualsLocalX[detector->name()] = new TH1F("residualX_MC_local", title.c_str(), 400, -0.2, 0.2);
            title = detector->name() + " Biased MC residual Y (local);y_{track}-y_{MC} [mm];events";
            telescopeMCresidualsLocalY[detector->name()] = new TH1F("residualY_MC_local", title.c_str(), 400, -0.2, 0.2);
            title = detector->name() + " Biased MC residual X (global);x_{track}-x_{MC} [mm];events";
            telescopeMCresidualsX[detector->name()] = new TH1F("residualX_MC_global", title.c_str(), 400, -0.2, 0.2);
            title = detector->name() + " Biased MC residual Y (global);y_{track}-y_{MC} [mm];events";
            telescopeMCresidualsY[detector->name()] = new TH1F("residualY_MC_global", title.c_str(), 400, -0.2, 0.2);
        }

        directory->cd();
    }
}

ROOT::Math::XYZPoint AnalysisTelescope::closestApproach(ROOT::Math::XYZPoint position, MCParticles* particles) {
    // Find the closest MC particle
    double smallestDistance(DBL_MAX);
    ROOT::Math::XYZPoint particlePosition;
    for(auto& particle : (*particles)) {
        ROOT::Math::XYZPoint entry = particle->getLocalStart();
        ROOT::Math::XYZPoint exit = particle->getLocalEnd();
        ROOT::Math::XYZPoint centre((entry.X() + exit.X()) / 2., (entry.Y() + exit.Y()) / 2., (entry.Z() + exit.Z()) / 2.);
        double distance = sqrt((centre.X() - position.X()) * (centre.X() - position.X()) +
                               (centre.Y() - position.Y()) * (centre.Y() - position.Y()));
        if(distance < smallestDistance) {
            particlePosition.SetXYZ(centre.X(), centre.Y(), centre.Z());
        }
    }
    return particlePosition;
}

74
StatusCode AnalysisTelescope::run(std::shared_ptr<Clipboard> clipboard) {
75
76
77
78
79

    // Get the tracks from the clipboard
    Tracks* tracks = reinterpret_cast<Tracks*>(clipboard->get("tracks"));
    if(tracks == nullptr) {
        LOG(DEBUG) << "No tracks on the clipboard";
80
        return StatusCode::Success;
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
    }

    for(auto& track : (*tracks)) {

        // Cut on the chi2/ndof
        if(track->chi2ndof() > chi2ndofCut) {
            continue;
        }

        // Loop over clusters of the track:
        for(auto& cluster : track->clusters()) {
            auto detector = get_detector(cluster->detectorID());
            if(detector == nullptr || detector->isDUT()) {
                continue;
            }

            auto name = detector->name();
98
            ROOT::Math::XYZPoint intercept = track->intercept(cluster->global().z());
99
            auto interceptLocal = detector->globalToLocal(intercept);
100
101
102
103
            telescopeResidualsLocalX[name]->Fill(cluster->local().x() - interceptLocal.X());
            telescopeResidualsLocalY[name]->Fill(cluster->local().y() - interceptLocal.Y());
            telescopeResidualsX[name]->Fill(cluster->global().x() - intercept.X());
            telescopeResidualsY[name]->Fill(cluster->global().y() - intercept.Y());
104
105
106
107
108
109
110
111

            // Get the MC particles from the clipboard
            MCParticles* mcParticles = reinterpret_cast<MCParticles*>(clipboard->get(name, "mcparticles"));
            if(mcParticles == nullptr) {
                continue;
            }

            ROOT::Math::XYZPoint particlePosition = closestApproach(cluster->local(), mcParticles);
112
113
            telescopeMCresidualsLocalX[name]->Fill(cluster->local().x() + detector->size().X() / 2 - particlePosition.X());
            telescopeMCresidualsLocalY[name]->Fill(cluster->local().y() + detector->size().Y() / 2 - particlePosition.Y());
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
            telescopeMCresidualsX[name]->Fill(interceptLocal.X() + detector->size().X() / 2 - particlePosition.X());
            telescopeMCresidualsY[name]->Fill(interceptLocal.Y() + detector->size().Y() / 2 - particlePosition.Y());
        }

        // Calculate telescope resolution at DUT
        for(auto& detector : get_detectors()) {
            if(!detector->isDUT()) {
                continue;
            }

            // Get the MC particles from the clipboard
            MCParticles* mcParticles = reinterpret_cast<MCParticles*>(clipboard->get(detector->name(), "mcparticles"));
            if(mcParticles == nullptr) {
                continue;
            }

            auto intercept = detector->getIntercept(track);
            auto interceptLocal = detector->globalToLocal(intercept);
            auto particlePosition = closestApproach(interceptLocal, mcParticles);

            telescopeResolutionX[detector->name()]->Fill(interceptLocal.X() + detector->size().X() / 2 -
                                                         particlePosition.X());
        }
    }

    // Return value telling analysis to keep running
140
    return StatusCode::Success;
141
}