Commit 05b43e73 authored by Simon Spannagel's avatar Simon Spannagel
Browse files

Merge branch 'master' of ssh://gitlab.cern.ch:7999/CLICdp/DAQs/tbAnalysis


Former-commit-id: 45071ad8005e5253659f8bed5697f1dc36bd0485
parents f0a832df 05a1146b
......@@ -60,6 +60,7 @@ ${BINDIR}EventDict.o : ${DATAOBJS}
@rootcint core/EventDict.C -c ${DATAOBJS}
@$(CC) $(CFLAGS) ${INCLUDE_PATH} -c core/EventDict.C -o ${BINDIR}EventDict.o
@rm -f core/EventDict.C core/EventDict.h
@mv core/EventDict_rdict.pcm bin/
${EXE} : ${BINDIR}EventDict.o ${OBJS} ${DATACLASSOBJS} ${ALGOBJS}
@echo "Making executable tbAnalysis"
......
......@@ -5,8 +5,8 @@
Alignment::Alignment(bool debugging)
: Algorithm("Alignment"){
debug = debugging;
m_numberOfTracksForAlignment = 400000;
nIterations = 5;
m_numberOfTracksForAlignment = 20000;
nIterations = 3;
}
// Global container declarations
......@@ -16,12 +16,11 @@ Parameters* globalParameters;
int detNum;
void Alignment::initialise(Parameters* par){
// Pick up the global parameters
parameters = par;
}
// During run, just pick up tracks and save them till the end
StatusCode Alignment::run(Clipboard* clipboard){
// Get the tracks
......@@ -50,7 +49,7 @@ StatusCode Alignment::run(Clipboard* clipboard){
// Minimisation functions for Minuit
// ========================================
// METHOD 0
// 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) {
......@@ -95,6 +94,7 @@ void MinimiseTrackChi2(Int_t &npar, Double_t *grad, Double_t &result, Double_t *
}
// 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) {
......@@ -148,7 +148,7 @@ 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;
//if(m_alignmenttracks.size() < m_numberOfTracksForAlignment) return;
// Make the fitting object
TVirtualFitter* residualFitter = TVirtualFitter::Fitter(0,50);
......
......@@ -10,9 +10,6 @@
#include "Cluster.h"
#include "Track.h"
// Fitting function for minuit
void SumDistance2(Int_t &, Double_t *, Double_t &, Double_t *, Int_t );
class Alignment : public Algorithm {
public:
......
......@@ -32,9 +32,9 @@ void BasicTracking::initialise(Parameters* par){
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);
residualsX[detectorID] = new TH1F(name.c_str(),name.c_str(),100,-0.02,0.02);
name = "residualsY_"+detectorID;
residualsY[detectorID] = new TH1F(name.c_str(),name.c_str(),400,-0.2,0.2);
residualsY[detectorID] = new TH1F(name.c_str(),name.c_str(),100,-0.02,0.02);
}
nTracksTotal = 0.;
......
#include "Clicpix2Correlator.h"
Clicpix2Correlator::Clicpix2Correlator(bool debugging)
: Algorithm("Clicpix2Correlator"){
debug = debugging;
}
template <typename T>
std::string makeString(T number) {
std::ostringstream ss;
ss << number;
return ss.str();
}
void Clicpix2Correlator::initialise(Parameters* par){
parameters = par;
// Get the DUT ID
dutID = parameters->DUT;
// Initialise histograms
hTrackDiffX["standard"] = new TH1F("hTrackDiffX_standard","hTrackDiffX_standard",4000,-20,20);
hTrackDiffY["standard"] = new TH1F("hTrackDiffY_standard","hTrackDiffY_standard",4000,-20,20);
hTrackDiffX["prevEvent"] = new TH1F("hTrackDiffX_prevEvent","hTrackDiffX_prevEvent",4000,-20,20);
hTrackDiffY["prevEvent"] = new TH1F("hTrackDiffY_prevEvent","hTrackDiffY_prevEvent",4000,-20,20);
// Rotatation histograms
angleStart = 0;
angleStep = 0.6;
angleStop = 2.*M_PI;
for(double angle=angleStart;angle<angleStop;angle+=angleStep){
string name = "rotated" + makeString(angle);
string histo = "hTrackDiffX_" + name;
hTrackDiffX[name] = new TH1F(histo.c_str(),histo.c_str(),4000,-20,20);
}
// Initialise member variables
m_eventNumber = 0;
}
StatusCode Clicpix2Correlator::run(Clipboard* clipboard){
// Get the clicpix clusters in this event
Clusters* clusters = (Clusters*)clipboard->get(dutID,"clusters");
if(clusters == NULL){
if(debug) tcout<<"No clusters for "<<dutID<<" on the clipboard"<<endl;
m_eventNumber++;
return Success;
}
// Get the tracks
Tracks* tracks = (Tracks*)clipboard->get("tracks");
if(tracks == NULL){
m_eventNumber++;
return Success;
}
// Make local copies of these objects
for(int iTrack=0; iTrack<tracks->size(); iTrack++){
Track* track = (*tracks)[iTrack];
Track* storageTrack = new Track(track);
m_eventTracks[m_eventNumber].push_back(storageTrack);
}
for(int iCluster=0; iCluster<clusters->size(); iCluster++){
Cluster* cluster = (*clusters)[iCluster];
Cluster* storageCluster = new Cluster(cluster);
m_eventClusters[m_eventNumber].push_back(storageCluster);
}
// Increment event counter
m_eventNumber++;
// Return value telling analysis to keep running
return Success;
}
void Clicpix2Correlator::finalise(){
if(debug) tcout<<"Analysed "<<m_eventNumber<<" events"<<endl;
// Now make all of the correlations. For each event loop over the clusters
// and tracks and make correlations between the two
// Will rotate the detector and look for correlations
for(double angle=angleStart;angle<angleStop;angle+=angleStep){
// Set the angle
parameters->detector[dutID]->rotationX(angle);
parameters->detector[dutID]->update();
int event;
for(event=0;event<m_eventNumber;event++){
// Get the clusters and tracks
Tracks tracks = m_eventTracks[m_eventNumber];
Clusters clusters = m_eventClusters[m_eventNumber];
// Loop over tracks and make correlations
for(int iTrack=0;iTrack<tracks.size();iTrack++){
// Get the track
Track* track = tracks[iTrack];
// Get the track intercept with the clicpix plane (global co-ordinates)
PositionVector3D< Cartesian3D<double> > trackIntercept = parameters->detector[dutID]->getIntercept(track);
// Loop over all clusters from this event
for(int iCluster=0;iCluster<clusters.size();iCluster++){
// Get the cluster
Cluster* cluster = clusters[iCluster];
// Get the distance between this cluster and the track intercept (global)
double xcorr = cluster->globalX()-trackIntercept.X();
double ycorr = cluster->globalY()-trackIntercept.Y();
// Fill histograms on correlations
string name = "rotated" + makeString(angle);
hTrackDiffX[name]->Fill(xcorr);
hTrackDiffY[name]->Fill(ycorr);
}
}
}
}
}
#ifndef Clicpix2Correlator_H
#define Clicpix2Correlator_H 1
#include "Algorithm.h"
#include <iostream>
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "Pixel.h"
#include "Cluster.h"
#include "Track.h"
#include <sstream>
class Clicpix2Correlator : public Algorithm {
public:
// Constructors and destructors
Clicpix2Correlator(bool);
~Clicpix2Correlator(){}
// Functions
void initialise(Parameters*);
StatusCode run(Clipboard*);
void finalise();
// Member variables
int m_eventNumber;
string dutID;
map<int, Clusters> m_eventClusters;
map<int, Tracks> m_eventTracks;
double angleStart, angleStop, angleStep;
// Histograms
map<string,TH1F*> hTrackDiffX;
map<string,TH1F*> hTrackDiffY;
};
#endif // Clicpix2Correlator_H
#include "Clicpix2EventLoader.h"
Clicpix2EventLoader::Clicpix2EventLoader(bool debugging)
: Algorithm("Clicpix2EventLoader"){
debug = debugging;
}
void Clicpix2EventLoader::initialise(Parameters* par){
parameters = par;
// File structure is RunX/CLICpix2/data.csv
// Take input directory from global parameters
string inputDirectory = parameters->inputDirectory + "/CLICpix2";
// Open the root directory
DIR* directory = opendir(inputDirectory.c_str());
if (directory == NULL){tcout<<"Directory "<<inputDirectory<<" does not exist"<<endl; return;}
dirent* entry; dirent* file;
// Read the entries in the folder
while (entry = readdir(directory)){
// Check for the data file
string filename = inputDirectory + "/" + entry->d_name;
if(filename.find(".csv") != string::npos){
m_filename = filename;
}
}
// If no data was loaded, give a warning
if(m_filename.length() == 0) tcout<<"Warning: no data file was found for CLICpix2 in "<<inputDirectory<<endl;
// Open the data file for later
m_file.open(m_filename.c_str());
// Make histograms for debugging
hHitMap = new TH2F("hitMap","hitMap",128,0,128,128,0,128);
hPixelToT = new TH1F("pixelToT","pixelToT",100,0,100);
hPixelsPerFrame = new TH1F("pixelsPerFrame","pixelsPerFrame",1000,0,1000);
// Initialise member variables
m_eventNumber = 0;
}
StatusCode Clicpix2EventLoader::run(Clipboard* clipboard){
// 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 long int shutterStartTimeInt, shutterStopTimeInt;
long double shutterStartTime, shutterStopTime;
string data; int npixels=0;
bool shutterOpen = false;
// Read file and load data
while(getline(m_file,data)){
// Check if this is a header
if(data.find("=====") != string::npos){
istringstream header(data);
char guff; int frameNumber;
header >> guff >> frameNumber >> guff;
continue;
}
// If we are at the end of a frame then the next character will be a "="
char c = m_file.peek();
if(strcmp(&c,"=") == 0) break;
// Otherwise load data
// If there is a colon, then this is a timestamp
if(data.find(":") != string::npos){
istringstream timestamp(data);
char colon; int value; long long int time;
timestamp >> value >> colon >> time;
if(value == 3){
shutterOpen = true;
shutterStartTimeInt = time;
}else if(value == 1 && shutterOpen){
shutterOpen = false;
shutterStopTimeInt = time;
}
continue;
}
// Otherwise pixel data
int row(0), col(0), hitFlag(0), tot(0), toa(0); char comma; long int time;
istringstream pixelData(data);
pixelData >> col >> comma >> row >> comma >> hitFlag >> comma >> tot >> comma >> toa;
tot++;
// If this pixel is masked, do not save it
if(parameters->detector[detectorID]->masked(col,row)) continue;
Pixel* pixel = new Pixel(detectorID,row,col,tot,0);
pixels->push_back(pixel);
npixels++;
hHitMap->Fill(col,row);
hPixelToT->Fill(tot);
}
// Now set the event time so that the Timepix3 data is loaded correctly
shutterStartTime = shutterStartTimeInt * 25./1000000000.;
shutterStopTime = shutterStopTimeInt * 25./1000000000.;
parameters->currentTime = shutterStartTime;
parameters->eventLength = (shutterStopTime-shutterStartTime);
// Put the data on the clipboard
if(pixels->size() > 0) clipboard->put(detectorID,"pixels",(TestBeamObjects*)pixels);
// Fill histograms
hPixelsPerFrame->Fill(npixels);
// Return value telling analysis to keep running
return Success;
}
void Clicpix2EventLoader::finalise(){
if(debug) tcout<<"Analysed "<<m_eventNumber<<" events"<<endl;
}
#ifndef Clicpix2EventLoader_H
#define Clicpix2EventLoader_H 1
#include "Algorithm.h"
#include <iostream>
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "Pixel.h"
#include "Cluster.h"
#include "Track.h"
#include <fstream>
#include <sstream>
#include <stdio.h>
#include <dirent.h>
#include <string.h>
class Clicpix2EventLoader : public Algorithm {
public:
// Constructors and destructors
Clicpix2EventLoader(bool);
~Clicpix2EventLoader(){}
// Functions
void initialise(Parameters*);
StatusCode run(Clipboard*);
void finalise();
// Histograms for several devices
map<string, TH2F*> plotPerDevice;
// Single histograms
TH1F* singlePlot;
// Member variables
int m_eventNumber;
string m_filename;
ifstream m_file;
TH2F* hHitMap;
TH1F* hPixelToT;
TH1F* hPixelsPerFrame;
};
#endif // Clicpix2EventLoader_H
......@@ -6,16 +6,25 @@
ClicpixAnalysis::ClicpixAnalysis(bool debugging)
: Algorithm("ClicpixAnalysis"){
debug = debugging;
m_associationCut = 0.1; // 100 um
m_proximityCut = 0.00001; // 125 um
m_associationCut = 0.05; // 100 um
m_proximityCut = 0.0005; // 125 um
timepix3Telescope = false;
}
template <typename T>
std::string convertToString(T number) {
std::ostringstream ss;
ss << number;
return ss.str();
}
void ClicpixAnalysis::initialise(Parameters* par){
parameters = par;
// Initialise member variables
m_eventNumber = 0; m_triggerNumber = 0; dutID = parameters->DUT; m_lostHits = 0.;
// Cluster/pixel histograms
hHitPixels = new TH2F("hHitPixels","hHitPixels",64,0,64,64,0,64);
hColumnHits = new TH1F("hColumnHits","hColumnHits",64,0,64);
......@@ -23,7 +32,7 @@ void ClicpixAnalysis::initialise(Parameters* par){
hClusterSizeAll = new TH1F("hClusterSizeAll","hClusterSizeAll",20,0,20);
hClusterTOTAll = new TH1F("hClusterTOTAll","hClusterTOTAll",50,0,50);
hClustersPerEvent = new TH1F("hClustersPerEvent","hClustersPerEvent",50,0,50);
hClustersPerEvent = new TH1F("hClustersPerEvent","hClustersPerEvent",500,0,500);
hClustersVersusEventNo = new TH1F("hClustersVersusEventNo","hClustersVersusEventNo",60000,0,60000);
hGlobalClusterPositions = new TH2F("hGlobalClusterPositions","hGlobalClusterPositions",200,-2.0,2.0,300,-1.,2);
......@@ -47,6 +56,9 @@ void ClicpixAnalysis::initialise(Parameters* par){
hGlobalResidualsYversusColWidth = new TH2F("hGlobalResidualsYversusColWidth","hGlobalResidualsYversusColWidth",30,0,30,600,-0.3,0.3);
hGlobalResidualsYversusRowWidth = new TH2F("hGlobalResidualsYversusRowWidth","hGlobalResidualsYversusRowWidth",30,0,30,600,-0.3,0.3);
hTrackInterceptRow = new TH1F("hTrackInterceptRow","hTrackInterceptRow",660,-1.,65.);
hTrackInterceptCol = new TH1F("hTrackInterceptCol","hTrackInterceptCol",660,-1.,65.);
hAbsoluteResidualMap = new TH2F("hAbsoluteResidualMap","hAbsoluteResidualMap",50,0,50,25,0,25);
hXresidualVersusYresidual = new TH2F("hXresidualVersusYresidual","hXresidualVersusYresidual",600,-0.3,0.3,600,-0.3,0.3);
......@@ -110,14 +122,25 @@ void ClicpixAnalysis::initialise(Parameters* par){
hChipEfficiencyMap = new TH2F("hChipEfficiencyMap","hChipEfficiencyMap",65,-0.5,64.5,65,-0.5,64.5);
hGlobalEfficiencyMap = new TH2F("hGlobalEfficiencyMap","hGlobalEfficiencyMap",200,-2.0,2.0,300,-1.,2);
hInterceptClusterSize1 = new TH2F("hInterceptClusterSize1","hInterceptClusterSize1",50,0,50,25,0,25);
hInterceptClusterSize2 = new TH2F("hInterceptClusterSize2","hInterceptClusterSize2",50,0,50,25,0,25);
hInterceptClusterSize3 = new TH2F("hInterceptClusterSize3","hInterceptClusterSize3",50,0,50,25,0,25);
hInterceptClusterSize4 = new TH2F("hInterceptClusterSize4","hInterceptClusterSize4",50,0,50,25,0,25);
hInterceptClusterSize1 = new TH2F("hInterceptClusterSize1","hInterceptClusterSize1",25,0,25,25,0,25);
hInterceptClusterSize2 = new TH2F("hInterceptClusterSize2","hInterceptClusterSize2",25,0,25,25,0,25);
hInterceptClusterSize3 = new TH2F("hInterceptClusterSize3","hInterceptClusterSize3",25,0,25,25,0,25);
hInterceptClusterSize4 = new TH2F("hInterceptClusterSize4","hInterceptClusterSize4",25,0,25,25,0,25);
m_nBinsX=32; m_nBinsY=32;
hMapClusterSizeAssociated = new TH2F("hMapClusterSizeAssociated","hMapClusterSizeAssociated",m_nBinsX,0,parameters->detector[dutID]->nPixelsX(),m_nBinsY,0,parameters->detector[dutID]->nPixelsY());
for(int x=0;x<m_nBinsX;x++){
for(int y=0;y<m_nBinsY;y++){
int id = x+y*m_nBinsX;
std::string name = "hMapClusterTOTAssociated1pix"+convertToString(id);
TH1F* hMapEntryClusterTOTAssociated1pix = new TH1F(name.c_str(),name.c_str(),50,0,50);
hMapClusterTOTAssociated1pix[id] = hMapEntryClusterTOTAssociated1pix;
}
}
// Initialise member variables
m_eventNumber = 0; m_triggerNumber = 0; dutID = parameters->DUT; m_lostHits = 0.;
}
StatusCode ClicpixAnalysis::run(Clipboard* clipboard){
......@@ -157,6 +180,9 @@ StatusCode ClicpixAnalysis::run(Clipboard* clipboard){
Track* track = (*itTrack);
if (!track) continue;
// Cut on the track chi2/ndof
if(track->chi2ndof() < 3.0) continue;
// Get the track intercept with the clicpix plane (global and local co-ordinates)
PositionVector3D< Cartesian3D<double> > trackIntercept = parameters->detector[dutID]->getIntercept(track);
PositionVector3D< Cartesian3D<double> > trackInterceptLocal = *(parameters->detector[dutID]->m_globalToLocal) * trackIntercept;
......@@ -259,10 +285,14 @@ StatusCode ClicpixAnalysis::run(Clipboard* clipboard){
hClusterSizeAssociated->Fill((*bestCluster)->size());
// hClusterWidthColAssociated->Fill((*bestCluster)->colWidth());
// hClusterWidthRowAssociated->Fill((*bestCluster)->rowWidth());
hMapClusterSizeAssociated->Fill(chipInterceptCol,chipInterceptRow,(*bestCluster)->tot());
if((*bestCluster)->size() == 1){
hClusterTOTAssociated1pix->Fill((*bestCluster)->tot());
hInterceptClusterSize1->Fill(pixelInterceptX,pixelInterceptY);
int id = floor(chipInterceptCol*m_nBinsX/parameters->detector[dutID]->nPixelsX())+floor(chipInterceptRow*m_nBinsY/parameters->detector[dutID]->nPixelsY())*m_nBinsX;
hMapClusterTOTAssociated1pix[id]->Fill((*bestCluster)->tot());
}
if((*bestCluster)->size() == 2){
hClusterTOTAssociated2pix->Fill((*bestCluster)->tot());
......@@ -302,10 +332,10 @@ StatusCode ClicpixAnalysis::run(Clipboard* clipboard){
}
}
// Increment counter for number of hits found this way
if(pixelMatch){
*/ if(pixelMatch){
m_lostHits++;
hTrackInterceptsChipLost->Fill(chipInterceptCol,chipInterceptRow);
}*/
}//*/
if(!pixelMatch) hTrackInterceptsChipUnassociated->Fill(chipInterceptCol,chipInterceptRow);
}
......@@ -404,9 +434,9 @@ void ClicpixAnalysis::fillClusterHistos(Clusters* clusters){
for (itc = clusters->begin(); itc != clusters->end(); ++itc) {
// Loop over pixels and check if there are pixels not known
Pixels pixels = (*itc)->pixels();
Pixels* pixels = (*itc)->pixels();
Pixels::iterator itp;
for (itp = pixels.begin(); itp != pixels.end(); itp++) {
for (itp = pixels->begin(); itp != pixels->end(); itp++) {
// Check if this clicpix frame is still the current
int pixelID = (*itp)->m_column + nCols*(*itp)->m_row;
if( m_hitPixels[pixelID] != (*itp)->m_adc){
......@@ -437,9 +467,9 @@ void ClicpixAnalysis::fillClusterHistos(Clusters* clusters){
void ClicpixAnalysis::fillResponseHistos(double trackInterceptX, double trackInterceptY, Cluster* cluster){
// Loop over pixels in the cluster and show their distance from the track intercept
Pixels pixels = cluster->pixels();
Pixels* pixels = cluster->pixels();
Pixels::iterator itp;
for (itp = pixels.begin(); itp != pixels.end(); itp++) {
for (itp = pixels->begin(); itp != pixels->end(); itp++) {
// Get the pixel
Pixel* pixel = (*itp);
......
......@@ -3,6 +3,7 @@
#include "Algorithm.h"
#include <iostream>
#include <sstream>
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
......@@ -49,6 +50,8 @@ public:
TH2F* hGlobalResidualsXversusRowWidth;
TH2F* hGlobalResidualsYversusColWidth;
TH2F* hGlobalResidualsYversusRowWidth;
TH1F* hTrackInterceptRow;
TH1F* hTrackInterceptCol;
TH2F* hAbsoluteResidualMap;
TH2F* hXresidualVersusYresidual;
TH1F* hAssociatedClustersPerEvent;<