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

Merge remote-tracking branch 'daniel/master'

parents c4d43463 41cef3ba
Pipeline #259294 passed with stages
in 3 minutes and 38 seconds
......@@ -21,6 +21,7 @@ void DUTAnalysis::initialise() {
tracksVersusTime = new TH1F("tracksVersusTime", "tracksVersusTime", 300000, 0, 300);
associatedTracksVersusTime = new TH1F("associatedTracksVersusTime", "associatedTracksVersusTime", 300000, 0, 300);
residualsX = new TH1F("residualsX", "residualsX", 400, -0.2, 0.2);
residualsX1pix = new TH1F("residualsX1pix", "residualsX1pix", 400, -0.2, 0.2);
residualsY = new TH1F("residualsY", "residualsY", 400, -0.2, 0.2);
clusterTotAssociated = new TH1F("clusterTotAssociated", "clusterTotAssociated", 20000, 0, 100000);
......@@ -46,6 +47,7 @@ void DUTAnalysis::initialise() {
if(m_useMCtruth) {
residualsXMCtruth = new TH1F("residualsXMCtruth", "residualsXMCtruth", 400, -0.2, 0.2);
telescopeResolution = new TH1F("telescopeResolution", "telescopeResolution", 400, -0.2, 0.2);
}
// Initialise member variables
......@@ -239,6 +241,8 @@ StatusCode DUTAnalysis::run(Clipboard* clipboard) {
LOG(TRACE) << "Found associated cluster";
associatedTracksVersusTime->Fill(Units::convert(track->timestamp(), "s"));
residualsX->Fill(xdistance);
if(cluster->size() == 1)
residualsX1pix->Fill(xdistance);
residualsY->Fill(ydistance);
clusterTotAssociated->Fill(cluster->tot());
clusterSizeAssociated->Fill(cluster->size());
......@@ -264,7 +268,9 @@ StatusCode DUTAnalysis::run(Clipboard* clipboard) {
particlePosition.SetXYZ(centre.X(), centre.Y(), centre.Z());
}
}
residualsXMCtruth->Fill(cluster->localX() - particlePosition.X());
residualsXMCtruth->Fill(cluster->localX() + 7.04 - particlePosition.X());
auto interceptLocal = *(detector->globalToLocal()) * intercept;
telescopeResolution->Fill(interceptLocal.X() + 7.04 - particlePosition.X());
}
// Fill power pulsing response
......
......@@ -24,6 +24,7 @@ namespace corryvreckan {
TH1F* tracksVersusTime;
TH1F* associatedTracksVersusTime;
TH1F* residualsX;
TH1F* residualsX1pix;
TH1F* residualsY;
TH1F* clusterTotAssociated;
TH1F* clusterSizeAssociated;
......@@ -36,6 +37,7 @@ namespace corryvreckan {
TH2F* residualsTimeVsSignal;
TH1F* residualsXMCtruth;
TH1F* telescopeResolution;
TH2F* hAssociatedTracksGlobalPosition;
TH2F* hUnassociatedTracksGlobalPosition;
......
......@@ -4,63 +4,103 @@ using namespace corryvreckan;
using namespace std;
EtaCorrection::EtaCorrection(Configuration config, std::vector<Detector*> detectors)
: Algorithm(std::move(config), std::move(detectors)) {}
: Algorithm(std::move(config), std::move(detectors)) {
m_DUT = m_config.get<std::string>("DUT");
m_chi2ndofCut = m_config.get<double>("chi2ndofCut", 100.);
m_etaFormulaX = m_config.get<std::string>("EtaFormulaX", "[0] + [1]*x + [2]*x^2 + [3]*x^3 + [4]*x^4 + [5]*x^5");
m_etaConstantsX = m_config.getArray<double>("EtaConstantsX", {});
m_etaFormulaY = m_config.get<std::string>("EtaFormulaY", "[0] + [1]*x + [2]*x^2 + [3]*x^3 + [4]*x^4 + [5]*x^5");
m_etaConstantsY = m_config.getArray<double>("EtaConstantsY", {});
}
void EtaCorrection::initialise() {
// Initialise histograms per device
for(auto& detector : get_detectors()) {
// Check if they are a Timepix3
if(detector->type() != "Timepix3")
continue;
// Simple histogram per device
string name = "plotForDevice_" + detector->name();
plotPerDevice[detector->name()] = new TH2F(name.c_str(), name.c_str(), 256, 0, 256, 256, 0, 256);
}
// Initialise single histograms
string name = "singlePlot";
singlePlot = new TH1F(name.c_str(), name.c_str(), 1000, 0, 1000);
m_detector = get_detector(m_DUT);
double pitchX = m_detector->pitchX();
double pitchY = m_detector->pitchY();
string name = "etaDistributionX";
m_etaDistributionX = new TH2F(name.c_str(), name.c_str(), 100., 0., pitchX, 100., 0., pitchY);
name = "etaDistributionY";
m_etaDistributionY = new TH2F(name.c_str(), name.c_str(), 100., 0., pitchX, 100., 0., pitchY);
name = "etaDistributionXprofile";
m_etaDistributionXprofile = new TProfile(name.c_str(), name.c_str(), 100., 0., pitchX, 0., pitchY);
name = "etaDistributionYprofile";
m_etaDistributionYprofile = new TProfile(name.c_str(), name.c_str(), 100., 0., pitchX, 0., pitchY);
name = "etaDistributionXcorrected";
m_etaDistributionXcorrected = new TH2F(name.c_str(), name.c_str(), 100., 0., pitchX, 100., 0., pitchY);
name = "etaDistributionYcorrected";
m_etaDistributionYcorrected = new TH2F(name.c_str(), name.c_str(), 100., 0., pitchX, 100., 0., pitchY);
// Initialise member variables
m_etaCorrectorX = new TF1("etaCorrectorX", m_etaFormulaX.c_str(), 0, m_detector->pitchX());
for(int x = 0; x < m_etaConstantsX.size(); x++)
m_etaCorrectorX->SetParameter(x, m_etaConstantsX[x]);
m_etaCorrectorY = new TF1("etaCorrectorY", m_etaFormulaY.c_str(), 0, m_detector->pitchY());
for(int y = 0; y < m_etaConstantsY.size(); y++)
m_etaCorrectorY->SetParameter(y, m_etaConstantsY[y]);
m_eventNumber = 0;
}
StatusCode EtaCorrection::run(Clipboard* clipboard) {
// Loop over all Timepix3 and make plots
for(auto& detector : get_detectors()) {
// Check if they are a Timepix3
if(detector->type() != "Timepix3")
continue;
// Get the tracks from the clipboard
Tracks* tracks = (Tracks*)clipboard->get("tracks");
if(tracks == NULL) {
LOG(DEBUG) << "No tracks on the clipboard";
return Success;
}
// Get the pixels
Pixels* pixels = (Pixels*)clipboard->get(detector->name(), "pixels");
if(pixels == NULL) {
LOG(DEBUG) << "Detector " << detector->name() << " does not have any pixels on the clipboard";
continue;
}
// Loop over all tracks and look at the associated clusters to plot the eta distribution
for(auto& track : (*tracks)) {
// Get the clusters
Clusters* clusters = (Clusters*)clipboard->get(detector->name(), "clusters");
if(clusters == NULL) {
LOG(DEBUG) << "Detector " << detector->name() << " does not have any clusters on the clipboard";
// Cut on the chi2/ndof
if(track->chi2ndof() > m_chi2ndofCut) {
continue;
}
// Loop over all pixels and make hitmaps
for(auto& pixel : (*pixels)) {
// Fill the plots for this device
plotPerDevice[detector->name()]->Fill(pixel->m_column, pixel->m_row);
// Get the in-pixel track intercept
PositionVector3D<Cartesian3D<double>> trackIntercept = m_detector->getIntercept(track);
PositionVector3D<Cartesian3D<double>> trackInterceptLocal = *(m_detector->globalToLocal()) * trackIntercept;
double pixelInterceptX = m_detector->inPixelX(trackInterceptLocal) / 1000.;
double pixelInterceptY = m_detector->inPixelY(trackInterceptLocal) / 1000.;
(pixelInterceptX > m_detector->pitchX() / 2. ? pixelInterceptX -= m_detector->pitchX() / 2.
: pixelInterceptX += m_detector->pitchX() / 2.);
(pixelInterceptY > m_detector->pitchY() / 2. ? pixelInterceptY -= m_detector->pitchY() / 2.
: pixelInterceptY += m_detector->pitchY() / 2.);
// Look at the associated clusters and plot the eta function
for(auto& dutCluster : track->associatedClusters()) {
// Only look at clusters from the DUT
if(dutCluster->detectorID() != m_DUT)
continue;
// Ignore single pixel clusters
if(dutCluster->size() == 1)
continue;
// Get the fraction along the pixel
double inPixelX = m_detector->pitchX() * (dutCluster->column() - floor(dutCluster->column()));
double inPixelY = m_detector->pitchY() * (dutCluster->row() - floor(dutCluster->row()));
if(dutCluster->columnWidth() == 2) {
m_etaDistributionX->Fill(inPixelX, pixelInterceptX);
m_etaDistributionXprofile->Fill(inPixelX, pixelInterceptX);
// Apply the eta correction
double newX = floor(dutCluster->localX() / m_detector->pitchX()) * m_detector->pitchX() +
m_etaCorrectorX->Eval(inPixelX);
PositionVector3D<Cartesian3D<double>> positionLocal(newX, dutCluster->localY(), 0);
PositionVector3D<Cartesian3D<double>> positionGlobal = *(m_detector->localToGlobal()) * positionLocal;
dutCluster->setClusterCentre(positionGlobal.X(), positionGlobal.Y(), positionGlobal.Z());
dutCluster->setClusterCentreLocal(positionLocal.X(), positionLocal.Y(), positionLocal.Z());
}
if(dutCluster->rowWidth() == 2) {
m_etaDistributionY->Fill(inPixelY, pixelInterceptY);
m_etaDistributionYprofile->Fill(inPixelY, pixelInterceptY);
}
}
}
// Fill single histogram
singlePlot->Fill(m_eventNumber);
// Increment event counter
m_eventNumber++;
......
......@@ -3,8 +3,12 @@
#include <iostream>
#include "TCanvas.h"
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "core/Detector.h"
#include "core/algorithm/Algorithm.h"
#include "objects/Cluster.h"
#include "objects/Pixel.h"
......@@ -23,14 +27,25 @@ namespace corryvreckan {
StatusCode run(Clipboard* clipboard);
void finalise();
// Histograms for several devices
std::map<std::string, TH2F*> plotPerDevice;
// Single histograms
TH1F* singlePlot;
// Member variables
std::string m_DUT;
int m_eventNumber;
double m_chi2ndofCut;
Detector* m_detector;
std::string m_etaFormulaX;
std::vector<double> m_etaConstantsX;
TF1* m_etaCorrectorX;
std::string m_etaFormulaY;
std::vector<double> m_etaConstantsY;
TF1* m_etaCorrectorY;
// Histograms
TH2F* m_etaDistributionX;
TH2F* m_etaDistributionY;
TProfile* m_etaDistributionXprofile;
TProfile* m_etaDistributionYprofile;
TH2F* m_etaDistributionXcorrected;
TH2F* m_etaDistributionYcorrected;
};
} // namespace corryvreckan
#endif // EtaCorrection_H
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment