ClusterFitter.h

Go to the documentation of this file.

#pragma once
 
#include<TGraph.h>
#include<TF1.h>
#include<TFitResult.h>
 
#include<vector>
 
#include "TRDConstants.h"
#include "Line.h"
 
using std::cerr;
using std::endl;
 
 
namespace ClusterFitter
{
    //__________________________________________________________________________________________________________________________
    std::vector<Line::FitLine2D> GetClusterProjectionLinearFits(const std::vector<double>& TrkClusters_x,    
                                                                const std::vector<double>& TrkClusters_y,    
                                                                const std::vector<double>& TrkClusters_z,    
                                                                const int& trkID,
                                                                const int& g3PID,
                                                                const bool& debug = false,
                                                                const double threshold = 1000           
                                                                )
    {
        if (debug == true)
        {
            printf("\nClusterFitter:");
            printf("\nX\n"); VectorUtilities::printVector(TrkClusters_x);
            printf("\nY\n"); VectorUtilities::printVector(TrkClusters_y);
            printf("\nZ\n"); VectorUtilities::printVector(TrkClusters_z);
        }
        
        
        if (TrkClusters_x.size() != TrkClusters_y.size() || TrkClusters_x.size() != TrkClusters_z.size())
        {
            cerr<<"\n\nERROR! Argument cluster vector lengths are not equal.\n";
            cerr<<"SizeX = "<<TrkClusters_x.size()<<", SizeY = "<<TrkClusters_y.size()<<", SizeZ = "<<TrkClusters_z.size()<<endl;
            cerr<<"Returning empty vector\n\n";
            return {};
        }
        
        if (debug == true)
            cerr<<"\nClusterFitter: number of clusters to be fitted for this tracks = "<<TrkClusters_x.size()<<endl;
        if (TrkClusters_x.size() < 2)
            return {}; 
        
        
        TGraph * graph_clusters_XZ = new TGraph(TrkClusters_x.size());
        TGraph * graph_clusters_YZ = new TGraph(TrkClusters_x.size());
        TGraph * graph_clusters_XY = new TGraph(TrkClusters_x.size());
        for(unsigned int iclust = 0; iclust < TrkClusters_x.size(); iclust++)
        {
            graph_clusters_XZ->SetPoint((int)iclust, TrkClusters_x[iclust], TrkClusters_z[iclust]);
            graph_clusters_YZ->SetPoint((int)iclust, TrkClusters_y[iclust], TrkClusters_z[iclust]);
            graph_clusters_XY->SetPoint((int)iclust, TrkClusters_x[iclust], TrkClusters_y[iclust]);
        }
        TF1 * lin_fit_XZ = new TF1("track_XZ", "[0]*x+[1]", TRD_XYMIN, TRD_XYMAX);  
        TF1 * lin_fit_YZ = new TF1("track_YZ", "[0]*x+[1]", TRD_XYMIN, TRD_XYMAX);  
        TF1 * lin_fit_XY = new TF1("track_XY", "[0]*x+[1]", TRD_XYMIN, TRD_XYMAX);  
 
        graph_clusters_XZ->Fit(lin_fit_XZ, "SQ"); //"V" for verbose fitting
        graph_clusters_YZ->Fit(lin_fit_YZ, "SQ");
        graph_clusters_XY->Fit(lin_fit_XY, "SQ");
        TFitResultPtr r1 = graph_clusters_XZ->Fit(lin_fit_XZ, "SQ");
        TFitResultPtr r2 = graph_clusters_YZ->Fit(lin_fit_YZ, "SQ");
        TFitResultPtr r3 = graph_clusters_XY->Fit(lin_fit_XY, "SQ");
        
        double fit_xz_chi2ndf = r1->Chi2()/r1->Ndf();
        double fit_yz_chi2ndf = r2->Chi2()/r2->Ndf();
        double fit_xy_chi2ndf = r3->Chi2()/r3->Ndf();
        
        delete graph_clusters_XZ;
        delete graph_clusters_YZ;
        delete graph_clusters_XY;
        delete lin_fit_XZ;
        delete lin_fit_YZ;
        delete lin_fit_XY;
        
        double mXZ = r1->Parameter(0);
        double mPerErrXZ = 100*fabs(r1->ParError(0)/mXZ);
        double sXZ = r1->Parameter(1);
        double sPerErrXZ = 100*fabs(r1->ParError(1)/sXZ);
        
        double mYZ = r2->Parameter(0);
        double mPerErrYZ = 100*fabs(r2->ParError(0)/mYZ);
        double sYZ = r2->Parameter(1);
        double sPerErrYZ = 100*fabs(r2->ParError(1)/sYZ);
        
        double mXY = r3->Parameter(0);
        double mPerErrXY = 100*fabs(r3->ParError(0)/mXY);
        double sXY = r3->Parameter(1);
        double sPerErrXY = 100*fabs(r3->ParError(1)/sXY);
            
        if (debug == true)
        {            
            printf("Fit XZ: M = %14.7f, %%errM = %14.7f, C = %14.7f, %%errC = %14.7f, Chi-sq/ndf = %14.7f\n", mXZ, mPerErrXZ, sXZ, sPerErrXZ, fit_xz_chi2ndf);
            printf("Fit YZ: M = %14.7f, %%errM = %14.7f, C = %14.7f, %%errC = %14.7f, Chi-sq/ndf = %14.7f\n", mYZ, mPerErrYZ, sYZ, sPerErrYZ, fit_yz_chi2ndf);
            printf("Fit XY: M = %14.7f, %%errM = %14.7f, C = %14.7f, %%errC = %14.7f, Chi-sq/ndf = %14.7f\n", mXY, mPerErrXY, sXY, sPerErrXY, fit_xy_chi2ndf);
        }
        
        if (fit_xz_chi2ndf > threshold || fit_yz_chi2ndf > threshold || fit_xy_chi2ndf > threshold)
        {
            if (debug == true)
                cerr<<"\n\nINFO: Fit is poorer than provided threshold. Returning empty vector.\n\n";
            return {}; 
        }
        
        Line::FitLine2D lineXZ (mXZ, sXZ, r1->ParError(0), r1->ParError(1), fit_xz_chi2ndf, trkID, g3PID);
        Line::FitLine2D lineYZ (mYZ, sYZ, r2->ParError(0), r2->ParError(1), fit_yz_chi2ndf, trkID, g3PID);
        Line::FitLine2D lineXY (mXY, sXY, r3->ParError(0), r3->ParError(1), fit_xy_chi2ndf, trkID, g3PID);
        return {lineXZ, lineYZ, lineXY};
    }
}