Commit e86292e5 authored by Daniel Hynds's avatar Daniel Hynds
Browse files

First commit of new testbeam software

parent 2a481b60
BINDIR = bin/
SRCDIR = src/
INCLUDE_PATH = -I core \
-I algorithms \
-I objs
# Compiler
CC = g++
# Compiler flags
ROOTCFLAGS := $(shell root-config --cflags)
ROOTLIBS := $(shell root-config --glibs)
CFLAGS = -fPIC -w -g -W ${ROOTCFLAGS}
LFLAGS = -O3 ${ROOTLIBS} -g -lGenVector -lMinuit
# Automatically decide what to compile
CORE = $(notdir $(wildcard core/*.C))
OBJS = $(CORE:.C=.o)
OBJS := $(addprefix ${BINDIR}, ${OBJS})
EXE = ${BINDIR}tbAnalysis
ALGORITHMS = $(notdir $(wildcard algorithms/*.C))
ALGOBJS = $(ALGORITHMS:.C=.ao)
ALGOBJS := $(addprefix ${BINDIR}, ${ALGOBJS})
DATAOBJS = $(wildcard ${PWD}/objs/*.h)
# Compile core, user algorithms and make an executable for each user algorithm
all : ${OBJS} ${ALGOBJS} ${EXE}
@echo "Done"
# Compile core algorithms (everything .C)
${BINDIR}%.o : core/%.C core/%.h
@echo "Compiling $(notdir $<)"
@$(CC) $(CFLAGS) ${INCLUDE_PATH} -c $< -o $@
${BINDIR}%.ao : algorithms/%.C algorithms/%.h
@echo "Compiling $(notdir $<)"
@$(CC) $(CFLAGS) ${INCLUDE_PATH} -c $< -o $@
${BINDIR}Steering.o : core/Steering.C
@echo "Compiling steering file"
@$(CC) $(CFLAGS) ${INCLUDE_PATH} -c core/Steering.C -o ${BINDIR}Steering.o
${BINDIR}EventDict.o : ${DATAOBJS}
@echo "Making event dictionary"
@rm -f core/EventDict.C core/EventDict.h
@rootcint core/EventDict.C -c ${DATAOBJS}
@$(CC) $(CFLAGS) ${INCLUDE_PATH} -c core/EventDict.C -o ${BINDIR}EventDict.o
${EXE} : ${OBJS} ${ALGOBJS} ${BINDIR}EventDict.o
@echo "Making executable tbAnalysis"
@$(CC) -o ${BINDIR}tbAnalysis $(OBJS) ${ALGOBJS} $(LFLAGS)
# Remove all executables and object files
clean:
@rm -f ${BINDIR}*
@echo "Cleaned"
#include "Alignment.h"
#include <TVirtualFitter.h>
Alignment::Alignment(bool debugging)
: Algorithm("Alignment"){
debug = debugging;
}
// Global container declarations
Timepix3Tracks globalTracks;
string detectorToAlign;
Parameters* globalParameters;
int detNum;
void Alignment::initialise(Parameters* par){
parameters = par;
}
int Alignment::run(Clipboard* clipboard){
// Get the tracks
Timepix3Tracks* tracks = (Timepix3Tracks*)clipboard->get("Timepix3","tracks");
if(tracks == NULL){
return 1;
}
// Make a local copy and store it
for(int iTrack=0; iTrack<tracks->size(); iTrack++){
Timepix3Track* track = (*tracks)[iTrack];
Timepix3Track* alignmentTrack = new Timepix3Track(track);
m_alignmenttracks.push_back(alignmentTrack);
}
return 1;
}
void Alignment::finalise(){
// Make the fitting object
TVirtualFitter* residualFitter = TVirtualFitter::Fitter(0,50);
// Tell it what to minimise
residualFitter->SetFCN(SumDistance2);
// Set the global parameters
globalTracks = m_alignmenttracks;
globalParameters = parameters;
// Set the printout arguments of the fitter
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
// 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 ndet = 0; ndet<parameters->nDetectors; ndet++){
// Check if they are a Timepix3
string detectorID = parameters->detectors[ndet];
if(parameters->detector[detectorID]->type() != "Timepix3") continue;
// Do not align the reference plane
if(detectorID == parameters->reference) continue;
if(detectorID == "W0013_F09") continue;
if(detectorID == "W0013_G02") continue;
// Say that this is the detector we align
detectorToAlign = detectorID;
detNum=det;
// Add the parameters to the fitter (z displacement not allowed to move!)
residualFitter->SetParameter(det*6+0,(detectorID+"_displacementX").c_str(),parameters->detector[detectorID]->displacementX(),0.01,-50,50);
residualFitter->SetParameter(det*6+1,(detectorID+"_displacementY").c_str(),parameters->detector[detectorID]->displacementY(),0.01,-50,50);
residualFitter->SetParameter(det*6+2,(detectorID+"_displacementZ").c_str(),parameters->detector[detectorID]->displacementZ(),0,-10,500);
residualFitter->SetParameter(det*6+3,(detectorID+"_rotationX").c_str(),parameters->detector[detectorID]->rotationX(),0.01,-6.30,6.30);
residualFitter->SetParameter(det*6+4,(detectorID+"_rotationY").c_str(),parameters->detector[detectorID]->rotationY(),0.01,-6.30,6.30);
residualFitter->SetParameter(det*6+5,(detectorID+"_rotationZ").c_str(),parameters->detector[detectorID]->rotationZ(),0.01,-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);
residualFitter->SetParameter(det*6+2,(detectorID+"_displacementZ").c_str(),residualFitter->GetParameter(det*6+2),0,-10,500);
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));
parameters->detector[detectorID]->displacementZ(residualFitter->GetParameter(det*6+2));
parameters->detector[detectorID]->rotationX(residualFitter->GetParameter(det*6+3));
parameters->detector[detectorID]->rotationY(residualFitter->GetParameter(det*6+4));
parameters->detector[detectorID]->rotationZ(residualFitter->GetParameter(det*6+5));
parameters->detector[detectorID]->update();
det++;
}
det = 0;
// Now list the new alignment parameters
for(int ndet = 0; ndet<parameters->nDetectors; ndet++){
// Check if they are a Timepix3
string detectorID = parameters->detectors[ndet];
if(parameters->detector[detectorID]->type() != "Timepix3") continue;
// Do not align the reference plane
if(detectorID == parameters->reference) continue;
if(detectorID == "W0013_F09") continue;
if(detectorID == "W0013_G02") continue;
// Get the alignment parameters
double displacementX = residualFitter->GetParameter(det*6+0);
double displacementY = residualFitter->GetParameter(det*6+1);
double displacementZ = residualFitter->GetParameter(det*6+2);
double rotationX = residualFitter->GetParameter(det*6+3);
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;
det++;
}
// Write the output alignment file
parameters->writeConditions();
}
void SumDistance2(Int_t &npar, Double_t *grad, Double_t &result, Double_t *par, Int_t flag) {
// cout<<"Parameter size "<<par.size()<<endl;
// Pick up new alignment conditions
globalParameters->detector[detectorToAlign]->displacementX(par[detNum*6 + 0]);
globalParameters->detector[detectorToAlign]->displacementY(par[detNum*6 + 1]);
globalParameters->detector[detectorToAlign]->displacementZ(par[detNum*6 + 2]);
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();
// cout<<"displacement x "<<par[detNum*6 + 0]<<endl;
// cout<<"displacement y "<<par[detNum*6 + 1]<<endl;
// cout<<"displacement z "<<par[detNum*6 + 2]<<endl;
// cout<<"rotation x "<<par[detNum*6 + 3]<<endl;
// cout<<"rotation y "<<par[detNum*6 + 4]<<endl;
// cout<<"rotation z "<<par[detNum*6 + 5]<<endl;
// cout<<"Updated alignmnet parameters for detector "<<detectorToAlign<<endl;
// The chi2 value to be returned
result = 0.;
// Loop over all tracks
for(int iTrack=0; iTrack<globalTracks.size(); iTrack++){
// Get the track
Timepix3Track* track = globalTracks[iTrack];
// Get all clusters on the track
Timepix3Clusters trackClusters = track->clusters();
// Find the cluster that needs to have its position recalculated
for(int iTrackCluster=0; iTrackCluster<trackClusters.size(); iTrackCluster++){
Timepix3Cluster* trackCluster = trackClusters[iTrackCluster];
string detectorID = trackCluster->detectorID();
// Recalculate the global position from the local
PositionVector3D<Cartesian3D<double> > positionLocal(trackCluster->localX(),trackCluster->localY(),trackCluster->localZ());
PositionVector3D<Cartesian3D<double> > positionGlobal = *(globalParameters->detector[detectorID]->m_localToGlobal) * positionLocal;
trackCluster->setClusterCentre(positionGlobal.X(), positionGlobal.Y(),positionGlobal.Z());
}
// Refit the track
track->fit();
// Add the new chi2
result += track->chi2();
}
// cout<<"Chi2: "<<result<<endl;
}
#ifndef ALIGNMENT_H
#define ALIGNMENT_H 1
#include "Algorithm.h"
#include "Timepix3Cluster.h"
#include "Timepix3Track.h"
// Fitting function for minuit
void SumDistance2(Int_t &, Double_t *, Double_t &, Double_t *, Int_t );
class Alignment : public Algorithm {
public:
// Constructors and destructors
Alignment(bool);
~Alignment(){}
// Functions
void initialise(Parameters*);
int run(Clipboard*);
void finalise();
// Member variables
Timepix3Tracks m_alignmenttracks;
int nIterations;
};
#endif // ALIGNMENT_H
#include "BasicTracking.h"
BasicTracking::BasicTracking(bool debugging)
: Algorithm("BasicTracking"){
debug = debugging;
}
void BasicTracking::initialise(Parameters* par){
parameters = par;
// Set up histograms
trackChi2 = new TH1F("trackChi2","trackChi2",150,0,150);
// Loop over all Timepix3
for(int det = 0; det<parameters->nDetectors; det++){
// Check if they are a Timepix3
string detectorID = parameters->detectors[det];
if(parameters->detector[detectorID]->type() != "Timepix3") continue;
string name = "residualsX_"+detectorID;
residualsX[detectorID] = new TH1F(name.c_str(),name.c_str(),400,-0.2,0.2);
name = "residualsY_"+detectorID;
residualsY[detectorID] = new TH1F(name.c_str(),name.c_str(),400,-0.2,0.2);
}
}
int BasicTracking::run(Clipboard* clipboard){
// Container for all clusters
map<string,Timepix3Clusters*> clusters;
vector<string> detectors;
// Track container
Timepix3Tracks* tracks = new Timepix3Tracks();
// Loop over all Timepix3 and get clusters
for(int det = 0; det<parameters->nDetectors; det++){
// Check if they are a Timepix3
string detectorID = parameters->detectors[det];
if(parameters->detector[detectorID]->type() != "Timepix3") continue;
// Get the clusters
Timepix3Clusters* tempClusters = (Timepix3Clusters*)clipboard->get(detectorID,"clusters");
if(tempClusters == NULL){
tcout<<"Detector "<<detectorID<<" does not have any clusters on the clipboard"<<endl;
}
// Store them
clusters[detectorID] = tempClusters;
detectors.push_back(detectorID);
}
// Use the first plane as a seeding plane. For something quick, look a cluster in < 100 ns in the next plane, and continue
string reference = parameters->reference;
map<Timepix3Cluster*, bool> used;
// Loop over all clusters
for(int iSeedCluster=0;iSeedCluster<clusters[reference]->size();iSeedCluster++){
// Make a new track
Timepix3Track* track = new Timepix3Track();
// Get the cluster
Timepix3Cluster* cluster = (*clusters[reference])[iSeedCluster];
// Add the cluster to the track
track->addCluster(cluster);
used[cluster] = true;
// Get the cluster time
long long int timestamp = cluster->timestamp();
// Loop over each subsequent plane and look for a cluster within 100 ns
for(int det=0; det<detectors.size(); det++){
if(detectors[det] == reference) continue;
Timepix3Cluster* newCluster = getNearestCluster(timestamp, (*clusters[detectors[det]]) );
if( ((newCluster->timestamp() - timestamp) / (4096.*40000000.)) > (10./1000000000.) ) continue;
// Check if spatially more than 200 um
if( abs(cluster->globalX() - newCluster->globalX()) > 0.3 || abs(cluster->globalY() - newCluster->globalY()) > 0.3 ) continue;
track->addCluster(newCluster);
}
// Now should have a track with one cluster from each plane
if(track->nClusters() < 5) continue;
// Fit the track and save it
track->fit();
tracks->push_back(track);
// Fill histograms
trackChi2->Fill(track->chi2());
// Make residuals
Timepix3Clusters trackClusters = track->clusters();
for(int iTrackCluster=0; iTrackCluster<trackClusters.size(); iTrackCluster++){
Timepix3Cluster* trackCluster = trackClusters[iTrackCluster];
string detectorID = trackCluster->detectorID();
ROOT::Math::XYZPoint intercept = track->intercept(trackCluster->globalZ());
residualsX[detectorID]->Fill(intercept.X() - trackCluster->globalX());
residualsY[detectorID]->Fill(intercept.Y() - trackCluster->globalY());
}
}
tcout<<"Made "<<tracks->size()<<" tracks"<<endl;
clipboard->put("Timepix3","tracks",(TestBeamObjects*)tracks);
return 1;
}
Timepix3Cluster* BasicTracking::getNearestCluster(long long int timestamp, Timepix3Clusters clusters){
Timepix3Cluster* bestCluster = NULL;
// Loop over all clusters and return the one with the closest timestamp
for(int iCluster=0;iCluster<clusters.size();iCluster++){
Timepix3Cluster* cluster = clusters[iCluster];
if(bestCluster == NULL) bestCluster = cluster;
if(abs(cluster->timestamp() - timestamp) < abs(bestCluster->timestamp()-timestamp)) bestCluster = cluster;
}
return bestCluster;
}
void BasicTracking::finalise(){
}
#ifndef BASICTRACKING_H
#define BASICTRACKING_H 1
#include "Algorithm.h"
#include <iostream>
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "Timepix3Pixel.h"
#include "Timepix3Cluster.h"
#include "Timepix3Track.h"
class BasicTracking : public Algorithm {
public:
// Constructors and destructors
BasicTracking(bool);
~BasicTracking(){}
// Functions
void initialise(Parameters*);
int run(Clipboard*);
void finalise();
Timepix3Cluster* getNearestCluster(long long int, Timepix3Clusters);
// Member variables
TH1F* trackChi2;
map<string,TH1F*> residualsX;
map<string,TH1F*> residualsY;
};
#endif // BASICTRACKING_H
#include "EventDisplay.h"
#include "TApplication.h"
#include "TPolyLine3D.h"
EventDisplay::EventDisplay(bool debugging)
: Algorithm("EventDisplay"){
debug = debugging;
}
void EventDisplay::initialise(Parameters* par){
parameters = par;
// Set up histograms
eventMap = new TH3F("eventMap","eventMap",100,-5,5,100,-5,5,105,-10,200);
}
int EventDisplay::run(Clipboard* clipboard){
// Get the tracks
Timepix3Tracks* tracks = (Timepix3Tracks*)clipboard->get("Timepix3","tracks");
if(tracks == NULL){
return 1;
}
TApplication* rootapp = new TApplication("example",0, 0);
TCanvas* canv = new TCanvas();
eventMap->DrawCopy("");
// Loop over all tracks. Add the cluster points to the histogram and draw a line for the tracks
for(int iTrack=0; iTrack<tracks->size(); iTrack++){
// Get the track
Timepix3Track* track = (*tracks)[iTrack];
// Get all clusters on the track
Timepix3Clusters trackClusters = track->clusters();
// Fill the event map
for(int iTrackCluster=0; iTrackCluster<trackClusters.size(); iTrackCluster++){
Timepix3Cluster* trackCluster = trackClusters[iTrackCluster];
eventMap->Fill(trackCluster->globalX(), trackCluster->globalY(), trackCluster->globalZ());
}
ROOT::Math::XYZPoint linestart = track->intercept(0);
ROOT::Math::XYZPoint linestop = track->intercept(200);
TPolyLine3D *line = new TPolyLine3D(2);
line->SetPoint(0,linestart.X(),linestart.Y(),linestart.Z());
line->SetPoint(1,linestop.X(),linestop.Y(),linestop.Z());
line->SetLineColor(2);
line->Draw("same");
}
eventMap->DrawCopy("same,box");
canv->Update();
rootapp->Run();
return 1;
}
void EventDisplay::finalise(){
}
#ifndef EVENTDISPLAY_H
#define EVENTDISPLAY_H 1
#include "Algorithm.h"
#include <iostream>
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TCanvas.h"
#include "Timepix3Pixel.h"
#include "Timepix3Cluster.h"
#include "Timepix3Track.h"
class EventDisplay : public Algorithm {
public:
// Constructors and destructors
EventDisplay(bool);
~EventDisplay(){}
// Functions
void initialise(Parameters*);
int run(Clipboard*);
void finalise();
TH3F* eventMap;
};
#endif // EVENTDISPLAY_H
#include "EventLoader.h"
#include "Pixel.h"
#include <dirent.h>
#include <sstream>
#include <string>
EventLoader::EventLoader(bool debugging)
: Algorithm("EventLoader"){
debug = debugging;
m_inputDirectory = "testDir";
m_nFiles = 0;
m_fileNumber = 0;
}
void EventLoader::initialise(Parameters* par){
parameters = par;
// Open the input directory and store all filenames
DIR* directory = opendir(m_inputDirectory.c_str());
dirent* file;
while (file = readdir(directory)){
if (string(file->d_name).find(".txt") != string::npos){
m_files.push_back(m_inputDirectory + "/" + (string)file->d_name);
m_nFiles++;
}
}
}
int EventLoader::run(Clipboard* clipboard){
cout<<"\rLooping over files: "<<m_fileNumber<<"/"<<m_nFiles<<flush;
// If there are no more files to open, finish the event loop
if(!m_currentFile.is_open() && (m_fileNumber == m_nFiles)){
cout<<endl;
return 0;
}
// If no files are open, look at the next input file
if(!m_currentFile.is_open()){
m_currentFile.open(m_files[m_fileNumber]);
m_fileNumber++;
}