Multiplet.cpp 3.99 KB
Newer Older
Lennart Huth's avatar
Lennart Huth committed
1
2
3
4
5
6
7
8
9
10
11
/**
 * @file
 * @brief Implementation of Multiplet base object
 *
 * @copyright Copyright (c) 2017-2020 CERN and the Corryvreckan authors.
 * This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md".
 * In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an
 * Intergovernmental Organization or submit itself to any jurisdiction.
 */

#include "Multiplet.hpp"
12
13
#include "TMath.h"
#include "core/utils/unit.h"
Lennart Huth's avatar
Lennart Huth committed
14
15
16
#include "exceptions.h"

using namespace corryvreckan;
17
18
19
Multiplet::Multiplet(std::shared_ptr<Track> upstream, std::shared_ptr<Track> downstream) : Track() {
    m_upstream = upstream;
    m_downstream = downstream;
Lennart Huth's avatar
Lennart Huth committed
20

21
22
    // All clusters from up- and downstream should be referenced from this track:
    for(auto& cluster : m_upstream->getClusters()) {
23
24
        this->addCluster(cluster);
    }
25
    for(auto& cluster : m_downstream->getClusters()) {
26
27
28
29
        this->addCluster(cluster);
    }
}

Lennart Huth's avatar
Lennart Huth committed
30
ROOT::Math::XYPoint Multiplet::getKinkAt(const std::string&) const {
31
    return ROOT::Math::XYPoint(0, 0);
32
33
34
35
}

void Multiplet::calculateChi2() {

36
    chi2_ = m_upstream->getChi2() + m_downstream->getChi2() + sqrt(m_offsetAtScatterer.Dot(m_offsetAtScatterer));
37
    ndof_ = m_upstream->getNdof() + m_downstream->getNdof();
38
    chi2ndof_ = chi2_ / ndof_;
39
40
41
}

void Multiplet::calculateResiduals() {
42
    for(auto c : track_clusters_) {
43
        auto cluster = dynamic_cast<Cluster*>(c.GetObject());
44
        residual_global_[cluster->detectorID()] = cluster->global() - getIntercept(cluster->global().z());
45
        if(get_plane(cluster->detectorID()) != nullptr)
46
            residual_local_[cluster->detectorID()] =
47
                cluster->local() - get_plane(cluster->detectorID())->getToLocal() * getIntercept(cluster->global().z());
48
    }
49
50
51
52
53
}

void Multiplet::fit() {

    // tracks
54
55
    m_positionAtScatterer = ((m_downstream->getIntercept(m_scattererPosition) -
                              (ROOT::Math::XYZPoint(0, 0, 0) - m_upstream->getIntercept(m_scattererPosition))) /
56
                             2.);
57
    m_offsetAtScatterer = m_downstream->getIntercept(m_scattererPosition) - m_upstream->getIntercept(m_scattererPosition);
58
59

    // Calculate the angle
Lennart Huth's avatar
Lennart Huth committed
60
    ROOT::Math::XYZVector slopeUp = m_upstream->getDirection(m_scattererPosition);
Lennart Huth's avatar
Lennart Huth committed
61
    ROOT::Math::XYZVector slopeDown = m_downstream->getDirection(m_scattererPosition);
Lennart Huth's avatar
Lennart Huth committed
62
    //
Lennart Huth's avatar
Lennart Huth committed
63
    ROOT::Math::XYZVector kinks = (slopeDown /= slopeDown.z()) - (slopeUp /= slopeUp.z());
Lennart Huth's avatar
Lennart Huth committed
64
    m_kinkAtScatterer = ROOT::Math::XYVector(kinks.x(), kinks.y());
65
66

    this->calculateChi2();
67
    this->calculateResiduals();
68
    isFitted_ = true;
69
70
}

71
72
73
74
ROOT::Math::XYZPoint Multiplet::getIntercept(double z) const {
    return z == m_scattererPosition
               ? m_positionAtScatterer
               : (z < m_scattererPosition ? m_upstream->getIntercept(z) : m_downstream->getIntercept(z));
75
76
}

Lennart Huth's avatar
Lennart Huth committed
77
ROOT::Math::XYZPoint Multiplet::getState(const std::string& detectorID) const {
78
79
    return getClusterFromDetector(detectorID)->global().z() <= m_scattererPosition ? m_upstream->getState(detectorID)
                                                                                   : m_downstream->getState(detectorID);
80
81
}

Lennart Huth's avatar
Lennart Huth committed
82
ROOT::Math::XYZVector Multiplet::getDirection(const std::string& detectorID) const {
83
84
    return getClusterFromDetector(detectorID)->global().z() <= m_scattererPosition ? m_upstream->getDirection(detectorID)
                                                                                   : m_downstream->getDirection(detectorID);
85
86
}

87
88
89
90
XYZVector Multiplet::getDirection(const double& z) const {
    return (z <= m_scattererPosition ? m_upstream->getDirection(z) : m_downstream->getDirection(z));
}

91
void Multiplet::print(std::ostream& out) const {
92
    out << "Multiplet " << this->m_scattererPosition << ", " << this->m_positionAtScatterer << ", "
93
94
        << this->m_offsetAtScatterer << ", " << this->m_kinkAtScatterer << ", " << this->chi2_ << ", " << this->ndof_ << ", "
        << this->chi2ndof_ << ", " << this->timestamp();
Lennart Huth's avatar
Lennart Huth committed
95
}