Commit dd5bbbf0 authored by Simon Spannagel's avatar Simon Spannagel
Browse files

Merge branch 'master' into clipboard

parents 77bdb090 b0d1b9a6
Pipeline #1079419 passed with stages
in 17 minutes and 58 seconds
......@@ -223,6 +223,9 @@ Configuration Detector::getConfiguration() const {
if(this->isReference()) {
roles.push_back("reference");
}
if(this->isAuxiliary()) {
roles.push_back("auxiliary");
}
if(!roles.empty()) {
config.setArray("role", roles);
......@@ -322,12 +325,10 @@ bool Detector::hitMasked(Track* track, int tolerance) const {
// Functions to get row and column from local position
double Detector::getRow(const PositionVector3D<Cartesian3D<double>> localPosition) const {
// (1-m_nPixelsY%2)/2. --> add 1/2 pixel pitch if even number of rows
double row = localPosition.Y() / m_pitch.Y() + static_cast<double>(m_nPixels.Y()) / 2. + (1 - m_nPixels.Y() % 2) / 2.;
return row;
}
double Detector::getColumn(const PositionVector3D<Cartesian3D<double>> localPosition) const {
// (1-m_nPixelsX%2)/2. --> add 1/2 pixel pitch if even number of columns
double column = localPosition.X() / m_pitch.X() + static_cast<double>(m_nPixels.X()) / 2. + (1 - m_nPixels.X() % 2) / 2.;
return column;
}
......@@ -336,14 +337,19 @@ double Detector::getColumn(const PositionVector3D<Cartesian3D<double>> localPosi
PositionVector3D<Cartesian3D<double>> Detector::getLocalPosition(double column, double row) const {
return PositionVector3D<Cartesian3D<double>>(
m_pitch.X() * (column - m_nPixels.X() / 2.), m_pitch.Y() * (row - m_nPixels.Y() / 2.), 0.);
m_pitch.X() * (column - (m_nPixels.X()) / 2.), m_pitch.Y() * (row - (m_nPixels.Y()) / 2.), 0.);
}
// Function to get in-pixel position
ROOT::Math::XYVector Detector::inPixel(const double column, const double row) const {
// a pixel ranges from (col-0.5) to (col+0.5)
return XYVector(m_pitch.X() * (column - floor(column) - 0.5), m_pitch.Y() * (row - floor(row) - 0.5));
}
ROOT::Math::XYVector Detector::inPixel(const PositionVector3D<Cartesian3D<double>> localPosition) const {
double column = getColumn(localPosition);
double row = getRow(localPosition);
return XYVector(m_pitch.X() * (column - floor(column)), m_pitch.Y() * (row - floor(row)));
return inPixel(column, row);
}
// Check if track position is within ROI:
......
......@@ -214,6 +214,14 @@ namespace corryvreckan {
// Function to get local position from column (x) and row (y) coordinates
PositionVector3D<Cartesian3D<double>> getLocalPosition(double column, double row) const;
/**
* Transformation from local (sensor) coordinates to in-pixel coordinates
* @param column Column address ranging from int_column-0.5*pitch to int_column+0.5*pitch
* @param row Row address ranging from int_column-0.5*pitch to int_column+0.5*pitch
* @return Position within a single pixel cell, given in units of length
*/
XYVector inPixel(const double column, const double row) const;
/**
* Transformation from local (sensor) coordinates to in-pixel coordinates
* @param localPosition Local position on the sensor
......
......@@ -21,7 +21,6 @@ std::shared_ptr<Detector> globalDetector;
AlignmentDUTResidual::AlignmentDUTResidual(Configuration config, std::shared_ptr<Detector> detector)
: Module(std::move(config), detector), m_detector(detector) {
m_numberOfTracksForAlignment = m_config.get<size_t>("number_of_tracks", 20000);
nIterations = m_config.get<size_t>("iterations", 3);
m_pruneTracks = m_config.get<bool>("prune_tracks", false);
......@@ -121,15 +120,6 @@ StatusCode AlignmentDUTResidual::run(std::shared_ptr<Clipboard> clipboard) {
}
}
// If we have enough tracks for the alignment, tell the event loop to finish
if(m_alignmenttracks.size() >= m_numberOfTracksForAlignment) {
LOG(STATUS) << "Accumulated " << m_alignmenttracks.size() << " tracks, interrupting processing.";
if(m_discardedtracks > 0) {
LOG(STATUS) << "Discarded " << m_discardedtracks << " input tracks.";
}
return StatusCode::EndRun;
}
// Otherwise keep going
return StatusCode::Success;
}
......@@ -198,6 +188,10 @@ void AlignmentDUTResidual::MinimiseResiduals(Int_t&, Double_t*, Double_t& result
void AlignmentDUTResidual::finalise() {
if(m_discardedtracks > 0) {
LOG(STATUS) << "Discarded " << m_discardedtracks << " input tracks.";
}
// Make the fitting object
TVirtualFitter* residualFitter = TVirtualFitter::Fitter(nullptr, 50);
residualFitter->SetFCN(MinimiseResiduals);
......
......@@ -59,7 +59,6 @@ namespace corryvreckan {
int m_discardedtracks{};
size_t nIterations;
size_t m_numberOfTracksForAlignment;
bool m_pruneTracks;
bool m_alignPosition;
bool m_alignOrientation;
......
......@@ -10,7 +10,6 @@ This module performs translational and rotational DUT alignment. The alignment i
This module uses tracks for alignment. The module moves the detector is is instantiated for and minimises the unbiased residuals calculated from the track intercepts with the plane.
### Parameters
* `number_of_tracks`: Number of tracks used in the alignment method chosen. Default value is `20000`.
* `iterations`: Number of times the chosen alignment method is to be iterated. Default value is `3`.
* `align_position`: Boolean to select whether to align the X and Y displacements of the detector or not. Note that the Z displacement is never aligned. The default value is `true`.
* `align_orientation`: Boolean to select whether to align the three rotations of the detector under consideration or not. The default value is `true`.
......@@ -31,7 +30,10 @@ For the detector under consideration, the following plots are produced:
### Usage
```toml
[Corryvreckan]
# The global track limit can be used to restrict the alignment:
number_of_tracks = 200000
[AlignmentDUTResidual]
number_of_tracks = 1000000
log_level = INFO
```
......@@ -16,7 +16,6 @@ AlignmentMillepede::AlignmentMillepede(Configuration config, std::vector<std::sh
: Module(std::move(config), std::move(detectors)) {
m_excludeDUT = m_config.get<bool>("exclude_dut", false);
m_numberOfTracksForAlignment = m_config.get<size_t>("number_of_tracks", 20000);
m_dofs = m_config.getArray<bool>("dofs", {});
m_nIterations = m_config.get<size_t>("iterations", 5);
......@@ -80,14 +79,6 @@ StatusCode AlignmentMillepede::run(std::shared_ptr<Clipboard> clipboard) {
Track* alignmentTrack = new Track(*track);
m_alignmenttracks.push_back(alignmentTrack);
}
// If we have enough tracks for the alignment, tell the event loop to finish
if(m_alignmenttracks.size() >= m_numberOfTracksForAlignment) {
LOG(STATUS) << "Accumulated " << m_alignmenttracks.size() << " tracks, interrupting processing.";
return StatusCode::EndRun;
}
// Otherwise keep going
return StatusCode::Success;
}
......
......@@ -95,7 +95,6 @@ namespace corryvreckan {
const unsigned int m);
TrackVector m_alignmenttracks;
size_t m_numberOfTracksForAlignment;
/// Number of global derivatives
unsigned int m_nagb;
......
......@@ -6,14 +6,13 @@
### Description
This implementation of the Millepede module has been taken from the [Kepler framework](https://gitlab.cern.ch/lhcb/Kepler) used for test beam data analysis within the LHCb collaboration. It has been written by Christoph Hombach and has seen contributions from Tim Evans and Heinrich Schindler. This version is only slightly adapted to the Corryvreckan framework.
The Millepede algorthm allows a simultaneous fit of both the tracks and the alignment constants.
The Millepede algorithm allows a simultaneous fit of both the tracks and the alignment constants.
The modules stops if the convergence, i.e. the absolute sum of all corrections over the total number of parameters, is smaller than the configured value.
### Parameters
* `exclude_dut` : Exclude the DUT from the alignment procedure. Default value
is `false`.
* `number_of_tracks`: Number of tracks used in the alignment method chosen. Default value is `20000`.
* `iterations`: Number of times the chosen alignment method is to be iterated. Default value is `3`.
* `dofs`: Degrees of freedom to be aligned. This parameter should be given as vector of six boolean values for the parameters "Translation X", "Translation Y", "Translation Z", "Rotation X", "Rotation Y" and "Rotation Z". The default setting is an alignment of all parameters except for "Translation Z", i.e. `dofs = true, true, false, true, true, true`.
* `residual_cut`: Residual cut to reject a track as an outlier. Default value is `0.05mm`;
......
......@@ -13,7 +13,6 @@ int detNum;
AlignmentTrackChi2::AlignmentTrackChi2(Configuration config, std::vector<std::shared_ptr<Detector>> detectors)
: Module(std::move(config), std::move(detectors)) {
m_numberOfTracksForAlignment = m_config.get<size_t>("number_of_tracks", 20000);
nIterations = m_config.get<size_t>("iterations", 3);
m_pruneTracks = m_config.get<bool>("prune_tracks", false);
......@@ -63,15 +62,6 @@ StatusCode AlignmentTrackChi2::run(std::shared_ptr<Clipboard> clipboard) {
m_alignmenttracks.push_back(alignmentTrack);
}
// If we have enough tracks for the alignment, tell the event loop to finish
if(m_alignmenttracks.size() >= m_numberOfTracksForAlignment) {
LOG(STATUS) << "Accumulated " << m_alignmenttracks.size() << " tracks, interrupting processing.";
if(m_discardedtracks > 0) {
LOG(INFO) << "Discarded " << m_discardedtracks << " input tracks.";
}
return StatusCode::EndRun;
}
// Otherwise keep going
return StatusCode::Success;
}
......@@ -130,8 +120,9 @@ void AlignmentTrackChi2::MinimiseTrackChi2(Int_t&, Double_t*, Double_t& result,
void AlignmentTrackChi2::finalise() {
// If not enough tracks were produced, do nothing
// if(m_alignmenttracks.size() < m_numberOfTracksForAlignment) return;
if(m_discardedtracks > 0) {
LOG(INFO) << "Discarded " << m_discardedtracks << " input tracks.";
}
// Make the fitting object
TVirtualFitter* residualFitter = TVirtualFitter::Fitter(nullptr, 50);
......
......@@ -36,7 +36,6 @@ namespace corryvreckan {
int m_discardedtracks{};
size_t nIterations;
size_t m_numberOfTracksForAlignment;
bool m_pruneTracks;
bool m_alignPosition;
bool m_alignOrientation;
......
......@@ -6,11 +6,10 @@
### Description
This module performs translational and rotational telescope plane alignment. The alignment is performed with respect to the reference plane set in the configuration file.
This module uses the tracks produced by the `BasicTracking` module to align the telescope planes. If fewer than half of the tracks have associated clusters, a warning is produced on terminal.
For each telescope detector except the reference plane, this method moves the detector, refits all of the tracks, and minimises the chi^2 of these new tracks. The parameters `detectorToAlign` and `DUT` are not used in this method as it automatically iterates through all telescope planes except the DUT (if a DUT is present).
This module uses tracks on the clipboard to align the telescope planes.
For each telescope detector except the reference plane, this method moves the detector, refits all of the tracks, and minimises the chi^2 of these new tracks. This method automatically iterates through all devices contributing to the track.
### Parameters
* `number_of_tracks`: Number of tracks used in the alignment method chosen. Default value is `20000`.
* `iterations`: Number of times the chosen alignment method is to be iterated. Default value is `3`.
* `align_position`: Boolean to select whether to align the X and Y displacements of the detector or not. Note that the Z displacement is never aligned. The default value is `true`.
* `align_orientation`: Boolean to select whether to align the three rotations of the detector under consideration or not. The default value is `true`.
......@@ -30,7 +29,10 @@ For each detector the following plots are produced:
### Usage
```toml
[Corryvreckan]
# The global track limit can be used to restrict the alignment:
number_of_tracks = 200000
[AlignmentTrackChi2]
number_of_tracks = 1000000
log_level = INFO
```
......@@ -107,7 +107,7 @@ void AnalysisDUT::initialise() {
// In-pixel studies:
auto pitch_x = static_cast<double>(Units::convert(m_detector->pitch().X(), "um"));
auto pitch_y = static_cast<double>(Units::convert(m_detector->pitch().Y(), "um"));
auto mod_axes = "x_{track} mod " + std::to_string(pitch_x) + "#mum;y_{track} mod " + std::to_string(pitch_y) + "#mum;";
std::string mod_axes = "in-pixel x_{track} [#mum];in-pixel y_{track} [#mum];";
// cut flow histogram
std::string title = m_detector->name() + ": number of tracks discarded by different cuts;cut type;tracks";
......@@ -118,60 +118,136 @@ void AnalysisDUT::initialise() {
hCutHisto->GetXaxis()->SetBinLabel(4, "Close to frame begin/end");
title = "DUT x resolution;" + mod_axes + "MAD(#Deltax) [#mum]";
rmsxvsxmym = new TProfile2D(
"rmsxvsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y);
rmsxvsxmym = new TProfile2D("rmsxvsxmym",
title.c_str(),
static_cast<int>(pitch_x),
-pitch_x / 2.,
pitch_x / 2.,
static_cast<int>(pitch_y),
-pitch_y / 2.,
pitch_y / 2.);
title = "DUT y resolution;" + mod_axes + "MAD(#Deltay) [#mum]";
rmsyvsxmym = new TProfile2D(
"rmsyvsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y);
rmsyvsxmym = new TProfile2D("rmsyvsxmym",
title.c_str(),
static_cast<int>(pitch_x),
-pitch_x / 2.,
pitch_x / 2.,
static_cast<int>(pitch_y),
-pitch_y / 2.,
pitch_y / 2.);
title = "DUT resolution;" + mod_axes + "MAD(#sqrt{#Deltax^{2}+#Deltay^{2}}) [#mum]";
rmsxyvsxmym = new TProfile2D(
"rmsxyvsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y);
rmsxyvsxmym = new TProfile2D("rmsxyvsxmym",
title.c_str(),
static_cast<int>(pitch_x),
-pitch_x / 2.,
pitch_x / 2.,
static_cast<int>(pitch_y),
-pitch_y / 2.,
pitch_y / 2.);
title = "DUT cluster charge map;" + mod_axes + "<cluster charge> [ke]";
qvsxmym = new TProfile2D(
"qvsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y, 0, 250);
qvsxmym = new TProfile2D("qvsxmym",
title.c_str(),
static_cast<int>(pitch_x),
-pitch_x / 2.,
pitch_x / 2.,
static_cast<int>(pitch_y),
-pitch_y / 2.,
pitch_y / 2.,
0,
250);
title = "DUT cluster charge map, Moyal approx;" + mod_axes + "cluster charge MPV [ke]";
qMoyalvsxmym = new TProfile2D(
"qMoyalvsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y, 0, 250);
qMoyalvsxmym = new TProfile2D("qMoyalvsxmym",
title.c_str(),
static_cast<int>(pitch_x),
-pitch_x / 2.,
pitch_x / 2.,
static_cast<int>(pitch_y),
-pitch_y / 2.,
pitch_y / 2.,
0,
250);
title = "DUT seed pixel charge map;" + mod_axes + "<seed pixel charge> [ke]";
pxqvsxmym = new TProfile2D(
"pxqvsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y, 0, 250);
pxqvsxmym = new TProfile2D("pxqvsxmym",
title.c_str(),
static_cast<int>(pitch_x),
-pitch_x / 2.,
pitch_x / 2.,
static_cast<int>(pitch_y),
-pitch_y / 2.,
pitch_y / 2.,
0,
250);
title = "DUT cluster size map;" + mod_axes + "<pixels/cluster>";
npxvsxmym = new TProfile2D(
"npxvsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y, 0, 4.5);
npxvsxmym = new TProfile2D("npxvsxmym",
title.c_str(),
static_cast<int>(pitch_x),
-pitch_x / 2.,
pitch_x / 2.,
static_cast<int>(pitch_y),
-pitch_y / 2.,
pitch_y / 2.,
0,
4.5);
title = "DUT 1-pixel cluster map;" + mod_axes + "clusters";
npx1vsxmym =
new TH2F("npx1vsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y);
npx1vsxmym = new TH2F("npx1vsxmym",
title.c_str(),
static_cast<int>(pitch_x),
-pitch_x / 2.,
pitch_x / 2.,
static_cast<int>(pitch_y),
-pitch_y / 2.,
pitch_y / 2.);
title = "DUT 2-pixel cluster map;" + mod_axes + "clusters";
npx2vsxmym =
new TH2F("npx2vsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y);
npx2vsxmym = new TH2F("npx2vsxmym",
title.c_str(),
static_cast<int>(pitch_x),
-pitch_x / 2.,
pitch_x / 2.,
static_cast<int>(pitch_y),
-pitch_y / 2.,
pitch_y / 2.);
title = "DUT 3-pixel cluster map;" + mod_axes + "clusters";
npx3vsxmym =
new TH2F("npx3vsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y);
npx3vsxmym = new TH2F("npx3vsxmym",
title.c_str(),
static_cast<int>(pitch_x),
-pitch_x / 2.,
pitch_x / 2.,
static_cast<int>(pitch_y),
-pitch_y / 2.,
pitch_y / 2.);
title = "DUT 4-pixel cluster map;" + mod_axes + "clusters";
npx4vsxmym =
new TH2F("npx4vsxmym", title.c_str(), static_cast<int>(pitch_x), 0, pitch_x, static_cast<int>(pitch_y), 0, pitch_y);
npx4vsxmym = new TH2F("npx4vsxmym",
title.c_str(),
static_cast<int>(pitch_x),
-pitch_x / 2.,
pitch_x / 2.,
static_cast<int>(pitch_y),
-pitch_y / 2.,
pitch_y / 2.);
// Efficiency maps
title = "hPixelEfficiencyMap" + mod_axes + "efficiency";
hPixelEfficiencyMap = new TProfile2D("hPixelEfficiencyMap",
"hPixelEfficiencyMap;column;row;efficiency",
static_cast<int>(pitch_x),
0,
pitch_x,
-pitch_x / 2.,
pitch_x / 2.,
static_cast<int>(pitch_y),
0,
pitch_y,
-pitch_y / 2.,
pitch_y / 2.,
0,
1);
title = "hChipEfficiencyMap;column; row; efficiency";
hChipEfficiencyMap = new TProfile2D("hChipEfficiencyMap",
"hChipEfficiencyMap;column;row;efficiency",
m_detector->nPixels().X(),
......@@ -182,6 +258,7 @@ void AnalysisDUT::initialise() {
m_detector->nPixels().Y(),
0,
1);
title = "hGlobalEfficiencyMap;efficiency";
hGlobalEfficiencyMap = new TProfile2D("hGlobalEfficiencyMap",
"hGlobalEfficiencyMap;column;row;efficiency",
300,
......
......@@ -37,10 +37,18 @@ void AnalysisEfficiency::initialise() {
if(nbins_x > 1e4 || nbins_y > 1e4) {
throw InvalidValueError(m_config, "inpixel_bin_size", "Too many bins for in-pixel histograms.");
}
std::string title = m_detector->name() + " Pixel efficiency map;x_{track} mod " + std::to_string(pitch_x) +
"#mum;y_{track} mod " + std::to_string(pitch_y) + "#mum;efficiency";
hPixelEfficiencyMap_trackPos =
new TProfile2D("pixelEfficiencyMap_trackPos", title.c_str(), nbins_x, 0, pitch_x, nbins_y, 0, pitch_y, 0, 1);
std::string title =
m_detector->name() + " Pixel efficiency map;in-pixel x_{track} [#mum];in-pixel y_{track} #mum;efficiency";
hPixelEfficiencyMap_trackPos = new TProfile2D("pixelEfficiencyMap_trackPos",
title.c_str(),
nbins_x,
-pitch_x / 2.,
pitch_x / 2.,
nbins_y,
-pitch_y / 2.,
pitch_y / 2.,
0,
1);
title = m_detector->name() + " Chip efficiency map;x [px];y [px];efficiency";
hChipEfficiencyMap_trackPos = new TProfile2D("chipEfficiencyMap_trackPos",
title.c_str(),
......
......@@ -228,22 +228,22 @@ void AnalysisTimingATLASpix::initialise() {
0,
m_detector->nPixels().Y());
hHitMapAssoc_inPixel = new TH2F("hitMapAssoc_inPixel",
"hitMapAssoc_inPixel; x_{track} mod 130 #mum; y_{track} mod 130 #mum",
"hitMapAssoc_inPixel; in-pixel x_{track} [#mum]; in-pixel y_{track} [#mum]",
static_cast<int>(pitch_x),
0,
pitch_x,
-pitch_x / 2.,
pitch_x / 2.,
static_cast<int>(pitch_y),
0,
pitch_y);
-pitch_y / 2.,
pitch_y / 2.);
hHitMapAssoc_inPixel_highCharge =
new TH2F("hitMapAssoc_inPixel_highCharge",
"hitMapAssoc_inPixel_highCharge; x_{track} mod 130 #mum; y_{track} mod 130 #mum",
"hitMapAssoc_inPixel_highCharge; in-pixel x_{track} [#mum]; in-pixel y_{track} [#mum]",
static_cast<int>(pitch_x),
0,
pitch_x,
-pitch_x / 2.,
pitch_x / 2.,
static_cast<int>(pitch_y),
0,
pitch_y);
-pitch_y / 2.,
pitch_y / 2.);
hClusterMapAssoc = new TH2F("hClusterMapAssoc",
"hClusterMapAssoc; x_{cluster} [px]; x_{cluster} [px]; # entries",
m_detector->nPixels().X(),
......
......@@ -28,6 +28,8 @@ void Clustering4D::initialise() {
clusterPositionGlobal = new TH2F("clusterPositionGlobal", title.c_str(), 400, -10., 10., 400, -10., 10.);
title = ";cluster timestamp [ns]; # events";
clusterTimes = new TH1F("clusterTimes", title.c_str(), 3e6, 0, 3e9);
title = m_detector->name() + " Cluster multiplicity;clusters;events";
clusterMultiplicity = new TH1F("clusterMultiplicity", title.c_str(), 50, 0, 50);
}
// Sort function for pixels from low to high times
......@@ -116,6 +118,8 @@ StatusCode Clustering4D::run(std::shared_ptr<Clipboard> clipboard) {
deviceClusters->push_back(cluster);
}
clusterMultiplicity->Fill(static_cast<double>(deviceClusters->size()));
// Put the clusters on the clipboard
clipboard->put(deviceClusters, m_detector->name());
LOG(DEBUG) << "Made " << deviceClusters->size() << " clusters for device " << m_detector->name();
......@@ -193,12 +197,11 @@ void Clustering4D::calculateClusterCentre(Cluster* cluster) {
return;
}
// Create object with local cluster position
PositionVector3D<Cartesian3D<double>> positionLocal(m_detector->pitch().X() * (column - m_detector->nPixels().X() / 2),
m_detector->pitch().Y() * (row - m_detector->nPixels().Y() / 2),
0);
// Calculate local cluster position
auto positionLocal = m_detector->getLocalPosition(column, row);
// Calculate global cluster position
PositionVector3D<Cartesian3D<double>> positionGlobal = m_detector->localToGlobal(positionLocal);
auto positionGlobal = m_detector->localToGlobal(positionLocal);
// Set the cluster parameters
cluster->setColumn(column);
......
......@@ -38,6 +38,7 @@ namespace corryvreckan {
TH1F* clusterCharge;
TH2F* clusterPositionGlobal;
TH1F* clusterTimes;
TH1F* clusterMultiplicity;
double timingCut;
int neighbour_radius_row;
......
......@@ -24,6 +24,7 @@ For each detector the following plots are produced:
* Cluster charge histogram
* 2D cluster positions in global coordinates
* Cluster times
* Cluster multiplicity
### Usage
```toml
......
......@@ -201,11 +201,10 @@ void ClusteringSpatial::calculateClusterCentre(Cluster* cluster) {
<< Units::display(cluster->timestamp(), "us");
// Create object with local cluster position
PositionVector3D<Cartesian3D<double>> positionLocal(m_detector->pitch().X() * (column - m_detector->nPixels().X() / 2.),
m_detector->pitch().Y() * (row - m_detector->nPixels().Y() / 2.),
0);
auto positionLocal = m_detector->getLocalPosition(column, row);
// Calculate global cluster position
PositionVector3D<Cartesian3D<double>> positionGlobal = m_detector->localToGlobal(positionLocal);
auto positionGlobal = m_detector->localToGlobal(positionLocal);
// Set the cluster parameters
cluster->setRow(row);
......
......@@ -16,17 +16,21 @@ void EtaCalculation::initialise() {
double pitchX = m_detector->pitch().X();
double pitchY = m_detector->pitch().Y();
std::string title = m_detector->name() + " #eta distribution X";
m_etaDistributionX = new TH2F("etaDistributionX", title.c_str(), 100., 0., pitchX, 100., 0., pitchY);
m_etaDistributionX =
new TH2F("etaDistributionX", title.c_str(), 100., -pitchX / 2., pitchX / 2., 100., -pitchY / 2., pitchY / 2.);
title = m_detector->name() + " #eta distribution Y";
m_etaDistributionY = new TH2F("etaDistributionY", title.c_str(), 100., 0., pitchX, 100., 0., pitchY);
m_etaDistributionY =
new TH2F("etaDistributionY", title.c_str(), 100., -pitchX / 2., pitchX / 2., 100., -pitchY / 2., pitchY / 2.);
title = m_detector->name() + " #eta profile X";
m_etaDistributionXprofile = new TProfile("etaDistributionXprofile", title.c_str(), 100., 0., pitchX, 0., pitchY);
m_etaDistributionXprofile =
new TProfile("etaDistributionXprofile", title.c_str(), 100., -pitchX / 2., pitchX / 2., -pitchY / 2., pitchY);
title = m_detector->name() + " #eta profile Y";
m_etaDistributionYprofile = new TProfile("etaDistributionYprofile", title.c_str(), 100., 0., pitchX, 0., pitchY);
m_etaDistributionYprofile =
new TProfile("etaDistributionYprofile", title.c_str(), 100., -pitchX / 2., pitchX / 2., -pitchY / 2., pitchY / 2.);
// Prepare fit functions - we need them for every detector as they might have different pitches
m_etaFitX = new TF1("etaFormulaX", m_etaFormulaX.c_str(), 0, pitchX);
m_etaFitY = new TF1("etaFormulaY", m_etaFormulaY.c_str(), 0, pitchY);
m_etaFitX = new TF1("etaFormulaX", m_etaFormulaX.c_str(), -pitchX / 2., pitchX / 2.);