Skip to content
GitLab
Commit 36f8f5d2 authored by Simon Spannagel's avatar Simon Spannagel
Browse files

Rework algorithms: replace Tee with LOG, put into corryvreckan namespace, adapt constructor

parent a71258a8
Hide whitespace changes
Inline Side-by-side
src/algorithms/Alignment.C
View file @ 36f8f5d2
......@@ -2,9 +2,10 @@
#include <TVirtualFitter.h>
Alignment::Alignment(bool debugging)
: Algorithm("Alignment"){
debug = debugging;
using namespace corryvreckan;
Alignment::Alignment(Configuration config, Clipboard* clipboard)
: Algorithm(std::move(config), clipboard){
m_numberOfTracksForAlignment = 20000;
nIterations = 3;
}
......@@ -22,23 +23,23 @@ void Alignment::initialise(Parameters* par){
// During run, just pick up tracks and save them till the end
StatusCode Alignment::run(Clipboard* clipboard){
// Get the tracks
Tracks* tracks = (Tracks*)clipboard->get("tracks");
if(tracks == NULL){
return Success;
}
// Make a local copy and store it
// Make a local copy and store it
for(int iTrack=0; iTrack<tracks->size(); iTrack++){
Track* track = (*tracks)[iTrack];
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) return Failure;
// Otherwise keep going
return Success;
......@@ -53,7 +54,7 @@ StatusCode Alignment::run(Clipboard* clipboard){
// This method will move the detector in question, refit all of the tracks, and try to minimise the
// track chi2. If there were no clusters from this detector on any tracks then it would do nothing!
void MinimiseTrackChi2(Int_t &npar, Double_t *grad, Double_t &result, Double_t *par, Int_t flag) {
// Pick up new alignment conditions
globalParameters->detector[detectorToAlign]->displacementX(par[detNum*6 + 0]);
globalParameters->detector[detectorToAlign]->displacementY(par[detNum*6 + 1]);
......@@ -61,13 +62,13 @@ void MinimiseTrackChi2(Int_t &npar, Double_t *grad, Double_t &result, Double_t *
globalParameters->detector[detectorToAlign]->rotationX(par[detNum*6 + 3]);
globalParameters->detector[detectorToAlign]->rotationY(par[detNum*6 + 4]);
globalParameters->detector[detectorToAlign]->rotationZ(par[detNum*6 + 5]);
// Apply new alignment conditions
globalParameters->detector[detectorToAlign]->update();
// The chi2 value to be returned
result = 0.;
// Loop over all tracks
for(int iTrack=0; iTrack<globalTracks.size(); iTrack++){
// Get the track
......@@ -87,18 +88,18 @@ void MinimiseTrackChi2(Int_t &npar, Double_t *grad, Double_t &result, Double_t *
// Refit the track
track->fit();
// Add the new chi2
result += track->chi2();
}
}
// METHOD 1
// This method will move the detector in question and try to minimise the (unbiased) residuals. It uses
// the associated cluster container on the track (no refitting of the track)
void MinimiseResiduals(Int_t &npar, Double_t *grad, Double_t &result, Double_t *par, Int_t flag) {
// Pick up new alignment conditions
globalParameters->detector[globalParameters->detectorToAlign]->displacementX(par[0]);
globalParameters->detector[globalParameters->detectorToAlign]->displacementY(par[1]);
......@@ -106,13 +107,13 @@ void MinimiseResiduals(Int_t &npar, Double_t *grad, Double_t &result, Double_t *
globalParameters->detector[globalParameters->detectorToAlign]->rotationX(par[3]);
globalParameters->detector[globalParameters->detectorToAlign]->rotationY(par[4]);
globalParameters->detector[globalParameters->detectorToAlign]->rotationZ(par[5]);
// Apply new alignment conditions
globalParameters->detector[globalParameters->detectorToAlign]->update();
// The chi2 value to be returned
result = 0.;
// Loop over all tracks
for(int iTrack=0; iTrack<globalTracks.size(); iTrack++){
// Get the track
......@@ -146,17 +147,17 @@ void MinimiseResiduals(Int_t &npar, Double_t *grad, Double_t &result, Double_t *
// ==================================================================
void Alignment::finalise(){
// If not enough tracks were produced, do nothing
//if(m_alignmenttracks.size() < m_numberOfTracksForAlignment) return;
// Make the fitting object
TVirtualFitter* residualFitter = TVirtualFitter::Fitter(0,50);
// Tell it what to minimise
if(parameters->alignmentMethod == 0)residualFitter->SetFCN(MinimiseTrackChi2);
if(parameters->alignmentMethod == 1)residualFitter->SetFCN(MinimiseResiduals);
// Set the global parameters
globalTracks = m_alignmenttracks;
globalParameters = parameters;
......@@ -165,18 +166,18 @@ void Alignment::finalise(){
Double_t arglist[10];
arglist[0] = 3;
residualFitter->ExecuteCommand("SET PRINT",arglist,1);
// Set some fitter parameters
arglist[0] = 1000; // number of function calls
arglist[1] = 0.001; // tolerance
// This has been inserted in a temporary way. If the alignment method is 1 then it will align the single detector and then
// return. This should be made into separate functions.
if(parameters->alignmentMethod == 1){
// Get the name of the detector to align
string detectorID = parameters->detectorToAlign;
// Add the parameters to the fitter (z displacement not allowed to move!)
residualFitter->SetParameter(0,(detectorID+"_displacementX").c_str(),parameters->detector[detectorID]->displacementX(),0.01,-50,50);
residualFitter->SetParameter(1,(detectorID+"_displacementY").c_str(),parameters->detector[detectorID]->displacementY(),0.01,-50,50);
......@@ -184,12 +185,12 @@ void Alignment::finalise(){
residualFitter->SetParameter(3,(detectorID+"_rotationX").c_str(),parameters->detector[detectorID]->rotationX(),0.001,-6.30,6.30);
residualFitter->SetParameter(4,(detectorID+"_rotationY").c_str(),parameters->detector[detectorID]->rotationY(),0.001,-6.30,6.30);
residualFitter->SetParameter(5,(detectorID+"_rotationZ").c_str(),parameters->detector[detectorID]->rotationZ(),0.001,-6.30,6.30);
for(int iteration=0;iteration<nIterations;iteration++){
// Fit this plane (minimising global track chi2)
residualFitter->ExecuteCommand("MIGRAD",arglist,2);
// Set the alignment parameters of this plane to be the optimised values from the alignment
parameters->detector[detectorID]->displacementX(residualFitter->GetParameter(0));
parameters->detector[detectorID]->displacementY(residualFitter->GetParameter(1));
......@@ -199,17 +200,17 @@ void Alignment::finalise(){
parameters->detector[detectorID]->rotationZ(residualFitter->GetParameter(5));
}
// Write the output alignment file
parameters->writeConditions();
return;
}
// Loop over all planes. For each plane, set the plane alignment parameters which will be varied, and
// then minimise the track chi2 (sum of biased residuals). This means that tracks are refitted with
// each minimisation step.
int det = 0;
for(int iteration=0;iteration<nIterations;iteration++){
......@@ -229,10 +230,10 @@ void Alignment::finalise(){
residualFitter->SetParameter(det*6+3,(detectorID+"_rotationX").c_str(),parameters->detector[detectorID]->rotationX(),0.001,-6.30,6.30);
residualFitter->SetParameter(det*6+4,(detectorID+"_rotationY").c_str(),parameters->detector[detectorID]->rotationY(),0.001,-6.30,6.30);
residualFitter->SetParameter(det*6+5,(detectorID+"_rotationZ").c_str(),parameters->detector[detectorID]->rotationZ(),0.001,-6.30,6.30);
// Fit this plane (minimising global track chi2)
residualFitter->ExecuteCommand("MIGRAD",arglist,2);
// Now that this device is fitted, set parameter errors to 0 so that they are not fitted again
residualFitter->SetParameter(det*6+0,(detectorID+"_displacementX").c_str(),residualFitter->GetParameter(det*6+0),0,-50,50);
residualFitter->SetParameter(det*6+1,(detectorID+"_displacementY").c_str(),residualFitter->GetParameter(det*6+1),0,-50,50);
......@@ -240,7 +241,7 @@ void Alignment::finalise(){
residualFitter->SetParameter(det*6+3,(detectorID+"_rotationX").c_str(),residualFitter->GetParameter(det*6+3),0,-6.30,6.30);
residualFitter->SetParameter(det*6+4,(detectorID+"_rotationY").c_str(),residualFitter->GetParameter(det*6+4),0,-6.30,6.30);
residualFitter->SetParameter(det*6+5,(detectorID+"_rotationZ").c_str(),residualFitter->GetParameter(det*6+5),0,-6.30,6.30);
// Set the alignment parameters of this plane to be the optimised values from the alignment
parameters->detector[detectorID]->displacementX(residualFitter->GetParameter(det*6+0));
parameters->detector[detectorID]->displacementY(residualFitter->GetParameter(det*6+1));
......@@ -251,7 +252,7 @@ void Alignment::finalise(){
parameters->detector[detectorID]->update();
det++;
}
}
det = 0;
......@@ -270,14 +271,11 @@ void Alignment::finalise(){
double rotationY = residualFitter->GetParameter(det*6+4);
double rotationZ = residualFitter->GetParameter(det*6+5);
tcout<<" Detector "<<detectorID<<" new alignment parameters: T("<<displacementX<<","<<displacementY<<","<<displacementZ<<") R("<<rotationX<<","<<rotationY<<","<<rotationZ<<")"<<endl;
LOG(INFO) <<" Detector "<<detectorID<<" new alignment parameters: T("<<displacementX<<","<<displacementY<<","<<displacementZ<<") R("<<rotationX<<","<<rotationY<<","<<rotationZ<<")";
det++;
}
// Write the output alignment file
parameters->writeConditions();
}
src/algorithms/Alignment.h
View file @ 36f8f5d2
......@@ -10,23 +10,24 @@
#include "objects/Cluster.h"
#include "objects/Track.h"
class Alignment : public Algorithm {
namespace corryvreckan {
class Alignment : public Algorithm {
public:
// Constructors and destructors
Alignment(bool);
~Alignment(){}
public:
// Constructors and destructors
Alignment(Configuration config, Clipboard* clipboard);
~Alignment(){}
// Functions
void initialise(Parameters*);
StatusCode run(Clipboard*);
void finalise();
// Functions
void initialise(Parameters*);
StatusCode run(Clipboard*);
void finalise();
// Member variables
Tracks m_alignmenttracks;
int nIterations;
int m_numberOfTracksForAlignment;
};
// Member variables
Tracks m_alignmenttracks;
int nIterations;
int m_numberOfTracksForAlignment;
};
}
#endif // ALIGNMENT_H
src/algorithms/BasicTracking.C
View file @ 36f8f5d2
......@@ -2,15 +2,14 @@
#include "objects/KDTree.h"
#include "TCanvas.h"
BasicTracking::BasicTracking(bool debugging)
: Algorithm("BasicTracking"){
debug = false;
using namespace corryvreckan;
BasicTracking::BasicTracking(Configuration config, Clipboard* clipboard)
: Algorithm(std::move(config), clipboard){
// Default values for cuts
timingCut = 200./1000000000.; // 200 ns
spatialCut = 0.2; // 200 um
minHitsOnTrack = 6;
}
......@@ -43,7 +42,7 @@ void BasicTracking::initialise(Parameters* par){
StatusCode BasicTracking::run(Clipboard* clipboard){
if(debug) tcout<<"Start of event"<<endl;
LOG(DEBUG) <<"Start of event";
// Container for all clusters, and detectors in tracking
map<string,KDTree*> trees;
vector<string> detectors;
......@@ -63,10 +62,10 @@ StatusCode BasicTracking::run(Clipboard* clipboard){
// Get the clusters
Clusters* tempClusters = (Clusters*)clipboard->get(detectorID,"clusters");
if(tempClusters == NULL){
if(debug) tcout<<"Detector "<<detectorID<<" does not have any clusters on the clipboard"<<endl;
LOG(DEBUG) <<"Detector "<<detectorID<<" does not have any clusters on the clipboard";
}else{
// Store them
if(debug) tcout<<"Picked up "<<tempClusters->size()<<" clusters from "<<detectorID<<endl;
LOG(DEBUG) <<"Picked up "<<tempClusters->size()<<" clusters from "<<detectorID;
if(firstDetector){referenceClusters = tempClusters; seedPlane = det;}
firstDetector = false;
if(tempClusters->size() == 0) continue;
......@@ -86,7 +85,7 @@ StatusCode BasicTracking::run(Clipboard* clipboard){
for(int iSeedCluster=0;iSeedCluster<nSeedClusters;iSeedCluster++){
// Make a new track
if(debug) tcout<<"Looking at seed cluster "<<iSeedCluster<<endl;
LOG(DEBUG) <<"Looking at seed cluster "<<iSeedCluster;
Track* track = new Track();
// Get the cluster
Cluster* cluster = (*referenceClusters)[iSeedCluster];
......@@ -136,12 +135,12 @@ StatusCode BasicTracking::run(Clipboard* clipboard){
if(parameters->excludedFromTracking.count(detectors[det]) != 0) continue;
// Get all neighbours within 200 ns
if(debug) tcout<<"Searching for neighbouring cluster on "<<detectors[det]<<endl;
if(debug) tcout<<"- cluster time is "<<cluster->timestamp()<<endl;
LOG(DEBUG)<<"Searching for neighbouring cluster on "<<detectors[det];
LOG(DEBUG) <<"- cluster time is "<<cluster->timestamp();
Cluster* closestCluster = NULL; double closestClusterDistance = spatialCut;
Clusters neighbours = trees[detectors[det]]->getAllClustersInTimeWindow(cluster,timingCut);
if(debug) tcout<<"- found "<<neighbours.size()<<" neighbours"<<endl;
LOG(DEBUG) <<"- found "<<neighbours.size()<<" neighbours";
// Now look for the spatially closest cluster on the next plane
double interceptX, interceptY;
......@@ -172,7 +171,7 @@ StatusCode BasicTracking::run(Clipboard* clipboard){
if(closestCluster == NULL) continue;
// Add the cluster to the track
if(debug) tcout<<"- added cluster to track"<<endl;
LOG(DEBUG) <<"- added cluster to track";
track->addCluster(closestCluster);
}//*/
......@@ -221,7 +220,7 @@ StatusCode BasicTracking::run(Clipboard* clipboard){
if(trees.count(detectorID) != 0) delete trees[detectorID];
}
if(debug) tcout<<"End of event"<<endl;
LOG(DEBUG) <<"End of event";
return Success;
}
......
src/algorithms/BasicTracking.h
View file @ 36f8f5d2
......@@ -10,39 +10,40 @@
#include "objects/Cluster.h"
#include "objects/Track.h"
class BasicTracking : public Algorithm {
public:
// Constructors and destructors
BasicTracking(bool);
~BasicTracking(){}
// Functions
void initialise(Parameters*);
StatusCode run(Clipboard*);
void finalise();
// Cluster* getNearestCluster(Cluster*, map<Cluster*, bool>, Clusters*);
Cluster* getNearestCluster(long long int, Clusters);
// Member variables
// Histograms
TH1F* trackChi2;
TH1F* clustersPerTrack;
TH1F* trackChi2ndof;
TH1F* tracksPerEvent;
TH1F* trackAngleX;
TH1F* trackAngleY;
map<string,TH1F*> residualsX;
map<string,TH1F*> residualsY;
// Cuts for tracking
double timingCut;
double spatialCut;
int minHitsOnTrack;
double nTracksTotal;
};
namespace corryvreckan {
class BasicTracking : public Algorithm {
public:
// Constructors and destructors
BasicTracking(Configuration config, Clipboard* clipboard);
~BasicTracking(){}
// Functions
void initialise(Parameters*);
StatusCode run(Clipboard*);
void finalise();
// Cluster* getNearestCluster(Cluster*, map<Cluster*, bool>, Clusters*);
Cluster* getNearestCluster(long long int, Clusters);
// Member variables
// Histograms
TH1F* trackChi2;
TH1F* clustersPerTrack;
TH1F* trackChi2ndof;
TH1F* tracksPerEvent;
TH1F* trackAngleX;
TH1F* trackAngleY;
map<string,TH1F*> residualsX;
map<string,TH1F*> residualsY;
// Cuts for tracking
double timingCut;
double spatialCut;
int minHitsOnTrack;
double nTracksTotal;
};
}
#endif // BASICTRACKING_H
src/algorithms/CLICpixEventLoader.C
View file @ 36f8f5d2
#include "CLICpixEventLoader.h"
CLICpixEventLoader::CLICpixEventLoader(bool debugging)
: Algorithm("CLICpixEventLoader"){
debug = debugging;
using namespace corryvreckan;
CLICpixEventLoader::CLICpixEventLoader(Configuration config, Clipboard* clipboard)
: Algorithm(std::move(config), clipboard){
m_filename = "";
}
void CLICpixEventLoader::initialise(Parameters* par){
parameters = par;
// File structure is RunX/CLICpix/RunX.dat
// Take input directory from global parameters
string inputDirectory = parameters->inputDirectory + "/CLICpix";
// Open the root directory
DIR* directory = opendir(inputDirectory.c_str());
if (directory == NULL){tcout<<"Directory "<<inputDirectory<<" does not exist"<<endl; return;}
if (directory == NULL){LOG(ERROR)<<"Directory "<<inputDirectory<<" does not exist"; return;}
dirent* entry; dirent* file;
// Read the entries in the folder
while (entry = readdir(directory)){
// Check for the data file
......@@ -29,13 +30,13 @@ void CLICpixEventLoader::initialise(Parameters* par){
m_filename = filename;
}
}
// If no data was loaded, give a warning
if(m_filename.length() == 0) tcout<<"Warning: no data file was found for CLICpix in "<<inputDirectory<<endl;
if(m_filename.length() == 0) LOG(WARNING) <<"No data file was found for CLICpix in "<<inputDirectory;
// Open the data file for later
m_file.open(m_filename.c_str());
// Make histograms for debugging
hHitMap = new TH2F("hitMap","hitMap",64,0,64,64,0,64);
hPixelToT = new TH1F("pixelToT","pixelToT",20,0,20);
......@@ -45,33 +46,33 @@ void CLICpixEventLoader::initialise(Parameters* par){
StatusCode CLICpixEventLoader::run(Clipboard* clipboard){
// tcout<<"Running"<<endl;
// LOG(TRACE) <<"Running";
// Assume that the CLICpix is the DUT (if running this algorithm
string detectorID = parameters->DUT;
// If have reached the end of file, close it and exit program running
if(m_file.eof()){
m_file.close();
return Failure;
}
// Otherwise load a new frame
// Pixel container, shutter information
Pixels* pixels = new Pixels();
long double shutterStartTime, shutterStopTime;
string data;
int npixels=0;
// Read file and load data
while(getline(m_file,data)){
// tcout<<"Data: "<<data<<endl;
// LOG(TRACE) <<"Data: "<<data;
// If line is empty then we have finished this event, stop looping
if(data.length() < 5) break;
// Check if this is a header/shutter/power info
if(data.find("PWR_RISE") != string::npos || data.find("PWR_FALL") != string::npos) continue;
if(data.find("SHT_RISE") != string::npos){
......@@ -80,7 +81,7 @@ StatusCode CLICpixEventLoader::run(Clipboard* clipboard){
istringstream header(data);
header >> name >> timeInt;
shutterStartTime = (double)timeInt/(40000000.);
// tcout<<"Shutter rise time: "<<shutterStartTime<<endl;
LOG(TRACE) <<"Shutter rise time: "<<shutterStartTime;
continue;
}
if(data.find("SHT_FALL") != string::npos){
......@@ -89,19 +90,19 @@ StatusCode CLICpixEventLoader::run(Clipboard* clipboard){
istringstream header(data);
header >> name >> timeInt;
shutterStopTime = (double)timeInt/(40000000.);
// tcout<<"Shutter fall time: "<<shutterStopTime<<endl;
LOG(TRACE) <<"Shutter fall time: "<<shutterStopTime;
continue;
}
// Otherwise load data
int row, col, counter, tot(0);
long int time;
// tcout<<"Pixel data: "<<data<<endl;
LOG(TRACE) <<"Pixel data: "<<data;
istringstream pixelData(data);
pixelData >> col >> row >> counter >> tot;
tot++;
row=63-row;
// tcout<<"New pixe: "<<col<<","<<row<<" with tot "<<tot<<endl;
LOG(TRACE) <<"New pixel: "<<col<<","<<row<<" with tot "<<tot;
// If this pixel is masked, do not save it
if(parameters->detector[detectorID]->masked(col,row)) continue;
......@@ -110,27 +111,27 @@ StatusCode CLICpixEventLoader::run(Clipboard* clipboard){
npixels++;
hHitMap->Fill(col,row);
hPixelToT->Fill(tot);
}
// Now set the event time so that the Timepix3 data is loaded correctly
parameters->currentTime = shutterStartTime;
parameters->eventLength = (shutterStopTime-shutterStartTime);
// tcout<<"Loaded "<<npixels<<" pixels"<<endl;
LOG(TRACE) <<"Loaded "<<npixels<<" pixels";
// Put the data on the clipboard
if(pixels->size() > 0) clipboard->put(detectorID,"pixels",(TestBeamObjects*)pixels);
// Fill histograms
hPixelsPerFrame->Fill(npixels);
hShutterLength->Fill(shutterStopTime-shutterStartTime);
// Return value telling analysis to keep running
return Success;
}
void CLICpixEventLoader::finalise(){
if(debug) tcout<<"Analysed "<<m_eventNumber<<" events"<<endl;
LOG(DEBUG) <<"Analysed "<<m_eventNumber<<" events";
}
src/algorithms/CLICpixEventLoader.h
View file @ 36f8f5d2
......@@ -14,28 +14,29 @@
#include <stdio.h>
#include <dirent.h>
class CLICpixEventLoader : public Algorithm {
public:
// Constructors and destructors
CLICpixEventLoader(bool);
~CLICpixEventLoader(){}
// Functions