L1AlgoPulls.h

Go to the documentation of this file.

#ifndef L1AlgoPulls_h
#define L1AlgoPulls_h
 
// #define BUILD_HISTO_FOR_EACH_STANTION
 
#ifdef BUILD_HISTO_FOR_EACH_STANTION
const int NStations = 8;
#else
const int NStations = 0;
#endif   // BUILD_HISTO_FOR_EACH_STANTION
 
#include "CbmL1.h"
#include "L1Algo/CbmL1Def.h"
#include "L1Algo/L1Algo.h"
#include "L1Algo/L1StsHit.h"
#include "L1Algo/L1TrackPar.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TStyle.h"
 
#include <iostream>
#include <vector>
 
using std::cout;
using std::endl;
using std::vector;
 
struct TL1TrackParameters
{
    double x, y, tx, ty, qp;
 
    static const int NParameters = 5;
 
    TL1TrackParameters() {};
    TL1TrackParameters(L1TrackPar& T, int i)
        : x(T.x[i])
        , y(T.y[i])
        , tx(T.tx[i])
        , ty(T.ty[i])
        , qp(T.qp[i]) {};
    TL1TrackParameters(CbmL1MCPoint& T)
        : x(T.x)
        , y(T.y)
        , tx(T.px / T.pz)
        , ty(T.py / T.pz)
        , qp(T.q / T.p) {};
 
    double operator[](int i)
    {
        switch (i) {
            case 0:
                return x;
            case 1:
                return y;
            case 2:
                return tx;
            case 3:
                return ty;
            case 4:
                return qp;
        };
    }
 
    TL1TrackParameters operator-(TL1TrackParameters& b)
    {
        TL1TrackParameters c;
        c.x = x - b.x;
        c.y = y - b.y;
        c.tx = tx - b.tx;
        c.ty = ty - b.ty;
        c.qp = qp - b.qp;
        return c;
    }
 
    TL1TrackParameters operator/(TL1TrackParameters& b)
    {
        TL1TrackParameters c;
        c.x = x / b.x;
        c.y = y / b.y;
        c.tx = tx / b.tx;
        c.ty = ty / b.ty;
        c.qp = qp / b.qp;
        return c;
    }
 
    void Print() { cout << x << " " << y << " " << tx << " " << ty << " " << qp << endl; }
};
 
const TString L1TrackParametersNames[TL1TrackParameters::NParameters] = {"x", "y", "tx", "ty", "qp"};
 
class L1AlgoPulls
{
  public:
    L1AlgoPulls()
        : fNAllPulls(0)
    {
        fGPulls.clear();
        for (int i = 0; i < NStations; i++)
            fStaPulls[i].clear();
    };
    void Init();
 
    //     void AddVec(L1TrackPar& T, THitI ih);
    void AddOne(L1TrackPar& T, int i, THitI ih);
    void Print();   // fast method to see pulls :)
    void Build(bool draw = 1);
 
    int fNAllPulls;   // number of AddOne calls
  private:
    void makeUpHisto(TH1* hist, TString title, float& sigma);
 
    CbmL1* fL1;
    vector<TL1TrackParameters> fGPulls;                // pulls for all stations
    vector<TL1TrackParameters> fStaPulls[NStations];   // pulls for each station
    TH1F* histoPull[(1 + NStations) * TL1TrackParameters::NParameters];
 
    vector<TL1TrackParameters> fGRes;   // residuals for all stations
    TH1F* histoRes[(1 + NStations) * TL1TrackParameters::NParameters];
 
    static const float TailCut = 5000.;   //
    static const float csCut = 5.;        // chi-square cut
    static const int textFont = 22;       // TNewRoman
    TStyle* histoStyle;
};
 
// ===================================================================================
 
void L1AlgoPulls::Init()
{
    fL1 = CbmL1::Instance();
 
    static bool first_call = 1;
    if (first_call) {
        //     TDirectory *curdir = gDirectory;
        //     fL1->histodir->cd();
        //     gDirectory->mkdir("L1AlgoPulls");
        //     gDirectory->cd("L1AlgoPulls");
 
        // add global pulls
        for (int i = 0; i < TL1TrackParameters::NParameters; i++) {
            TString name = "pull_";
            name += L1TrackParametersNames[i];
            histoPull[i] = new TH1F(name, name, 50, -10, 10);
        }
 
#ifdef BUILD_HISTO_FOR_EACH_STANTION
        // add station pulls
        for (int i = TL1TrackParameters::NParameters; i < (NStations + 1) * TL1TrackParameters::NParameters; i++) {
            int ista = i / TL1TrackParameters::NParameters - 1;
            TString name = "pull_sta";
            name += ista;
            name += "_";
            name += L1TrackParametersNames[i % TL1TrackParameters::NParameters];
            histoPull[i] = new TH1F(name, name, 50, -10, 10);
        }
#endif   // BUILD_HISTO_FOR_EACH_STANTION
         // add global residuals
        for (int i = 0; i < TL1TrackParameters::NParameters; i++) {
            TString name = "residual_";
            name += L1TrackParametersNames[i];
            float size = 10;
            switch (i) {
                case 0:
                    size = .01;
                    break;
                case 1:
                    size = .01;
                    break;
                case 2:
                    size = 0.003;
                    break;
                case 3:
                    size = 0.003;
                    break;
                case 4:
                    size = 0.1;
                    break;
            };
            histoRes[i] = new TH1F(name, name, 50, -size, size);
        }
 
        // define style
        histoStyle = new TStyle("histoStyle", "Plain Style(no colors/fill areas)");
 
        histoStyle->SetTextFont(textFont);
        histoStyle->SetPadColor(0);
        histoStyle->SetCanvasColor(0);
        histoStyle->SetTitleColor(0);
        histoStyle->SetStatColor(0);
 
        histoStyle->SetOptTitle(0);   // without main up title
        histoStyle->SetOptStat(1000001010);
        //   The parameter mode can be = ksiourmen  (default = 000001111)
        //   k = 1;  kurtosis printed
        //   k = 2;  kurtosis and kurtosis error printed
        //   s = 1;  skewness printed
        //   s = 2;  skewness and skewness error printed
        //   i = 1;  integral of bins printed
        //   o = 1;  number of overflows printed
        //   u = 1;  number of underflows printed
        //   r = 1;  rms printed
        //   r = 2;  rms and rms error printed
        //   m = 1;  mean value printed
        //   m = 2;  mean and mean error values printed
        //   e = 1;  number of entries printed
        //   n = 1;  name of histogram is printed
 
        histoStyle->SetOptFit(10001);
        //   The parameter mode can be = pcev  (default = 0111)
        //   p = 1;  print Probability
        //   c = 1;  print Chisquare/Number of degress of freedom
        //   e = 1;  print errors (if e=1, v must be 1)
        //   v = 1;  print name/values of parameters
 
        histoStyle->SetStatW(0.175);
        histoStyle->SetStatH(0.02);
        histoStyle->SetStatX(0.95);
        histoStyle->SetStatY(0.97);
        histoStyle->SetStatFontSize(0.05);
 
        histoStyle->SetStatFont(textFont);
 
        //     gDirectory->cd("..");
        //     curdir->cd();
        first_call = 0;
    }
};
 
// inline void L1AlgoPulls::AddVec(L1TrackPar& T_, THitI ih)
// {
//   for (int i = 0; i < fvecLen; i++)
//     AddOne(T_,i,ih);
// }
 
inline void L1AlgoPulls::AddOne(L1TrackPar& T_, int i, THitI ih)
{
    fNAllPulls++;
    TL1TrackParameters T(T_, i);
 
    if (T_.chi2[i] > csCut * T_.NDF[i])
        return;
    // get error
    TL1TrackParameters err;
    if (!(finite(T_.C00[i]) && T_.C00[i] > 0))
        return;
    if (!(finite(T_.C11[i]) && T_.C11[i] > 0))
        return;
    if (!(finite(T_.C22[i]) && T_.C22[i] > 0))
        return;
    if (!(finite(T_.C33[i]) && T_.C33[i] > 0))
        return;
    if (!(finite(T_.C44[i]) && T_.C44[i] > 0))
        return;
    err.x = sqrt(T_.C00[i]);
    err.y = sqrt(T_.C11[i]);
    err.tx = sqrt(T_.C22[i]);
    err.ty = sqrt(T_.C33[i]);
    err.qp = sqrt(T_.C44[i]);
 
    // mc data
    int iMCP = fL1->vHitMCRef[ih];
    if (iMCP < 0)
        return;
    CbmL1MCPoint& mcP = fL1->vMCPoints[iMCP];
    TL1TrackParameters mc(mcP);
 
    // fill residuals
    TL1TrackParameters res = (mc - T);
    fGRes.push_back(res);
 
    // fill pulls
    TL1TrackParameters P = res / err;
    fGPulls.push_back(P);
 
#ifdef BUILD_HISTO_FOR_EACH_STANTION
    int ista = mcP.iStation - 2;   // last hit
                                   //   int ista = fL1->algo->vSFlag[ fL1->algo->vStsHits[ih].f ]/4 - 2; // last hit
    fStaPulls[ista].push_back(P);
#endif   // BUILD_HISTO_FOR_EACH_STANTION
};
 
inline void L1AlgoPulls::Print()
{   // TODO: renew
    cout << "All pulls:     " << fNAllPulls << endl;
    cout << "Correct pulls: " << fGPulls.size() << endl;
    cout << "x     y    tx    ty    qp" << endl;
    for (int i = 0; i < fGPulls.size(); i++) {
        TL1TrackParameters& pull = fGPulls[i];
        pull.Print();
    }
};
 
inline void L1AlgoPulls::Build(bool draw)
{
    // --- fill histograms ---
    // global pulls
    for (int i = 0; i < fGPulls.size(); i++) {
        TL1TrackParameters& pull = fGPulls[i];
        for (int ih = 0; ih < TL1TrackParameters::NParameters; ih++) {
            if (TailCut > fabs(pull[ih]))
                histoPull[ih]->Fill(pull[ih]);
        }
    }
#ifdef BUILD_HISTO_FOR_EACH_STANTION
    // station pulls
    for (int iSta = 0; iSta < NStations; iSta++) {
        vector<TL1TrackParameters>& Pulls = fStaPulls[iSta];
        for (int i = 0; i < Pulls.size(); i++) {
            TL1TrackParameters& pull = Pulls[i];
            for (int ih = 0; ih < TL1TrackParameters::NParameters; ih++) {
                if (TailCut > fabs(pull[ih]))
                    histoPull[(iSta + 1) * TL1TrackParameters::NParameters + ih]->Fill(pull[ih]);
            }
        }
    }
#endif   // BUILD_HISTO_FOR_EACH_STANTION
 
    // global residuals
    for (int i = 0; i < fGRes.size(); i++) {
        TL1TrackParameters& res = fGRes[i];
        for (int ih = 0; ih < TL1TrackParameters::NParameters; ih++) {
            if (TailCut > fabs(res[ih]))
                histoRes[ih]->Fill(res[ih]);
        }
    }
 
    // --- draw histograms --- and save info
    float pulls[(NStations + 1) * TL1TrackParameters::NParameters][2],   // 0 - sigma, 1 - RMS
        residuals[(NStations + 1) * TL1TrackParameters::NParameters][2];
 
    system("mkdir L1_Pulls -p");
    chdir("L1_Pulls");
    TCanvas* c2 = new TCanvas("c2", "c2", 0, 0, 600, 400);
    c2->cd();
    histoStyle->cd();
    for (int ih = 0; ih < (NStations + 1) * TL1TrackParameters::NParameters; ih++) {
        makeUpHisto(histoPull[ih], histoPull[ih]->GetName(), pulls[ih][0]);
        pulls[ih][1] = histoPull[ih]->GetRMS();
        if (draw) {
            histoPull[ih]->Draw();
            TString name = histoPull[ih]->GetName();
            name += ".pdf";
            c2->SaveAs(name);
        }
    }
    for (int ih = 0; ih < (0 + 1) * TL1TrackParameters::NParameters; ih++) {
        makeUpHisto(histoRes[ih], histoRes[ih]->GetName(), residuals[ih][0]);
        residuals[ih][1] = histoRes[ih]->GetRMS();
        if (draw) {
            histoRes[ih]->Draw();
            TString name = histoRes[ih]->GetName();
            name += ".pdf";
            c2->SaveAs(name);
        }
    }
    chdir("..");
 
    // --- print information ---
    cout << "All entries:     " << fNAllPulls << endl;
    cout << "Correct entries: " << fGPulls.size() << endl;
    cout << "Pulls sigma & RMS: " << endl;
    for (int ih = 0; ih < (NStations + 1) * TL1TrackParameters::NParameters; ih++) {
        int ipar = ih % TL1TrackParameters::NParameters;
        int ista = ih / TL1TrackParameters::NParameters;
        if ((ista > 0) && (ipar == 0))
            cout << "Station " << ista - 1 << endl;
        cout << L1TrackParametersNames[ipar] << "\t" << pulls[ih][0] << "\t" << pulls[ih][1] << endl;
    }
    cout << "Residuals sigma & RMS: " << endl;
    for (int ih = 0; ih < (0 + 1) * TL1TrackParameters::NParameters; ih++) {
        int ipar = ih % TL1TrackParameters::NParameters;
        cout << L1TrackParametersNames[ipar] << "\t" << residuals[ih][0] << "\t" << residuals[ih][1] << endl;
    }
};
 
inline void L1AlgoPulls::makeUpHisto(TH1* hist, TString title, float& sigma)
{
    if (hist && (hist->GetEntries() != 0)) {
        TF1* fit = new TF1("fit", "gaus");
        fit->SetLineColor(2);
        fit->SetLineWidth(3);
        hist->Fit("fit", "", "", hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax());
        sigma = fit->GetParameter(2);
 
        hist->GetXaxis()->SetLabelFont(textFont);
        hist->GetXaxis()->SetTitleFont(textFont);
        hist->GetYaxis()->SetLabelFont(textFont);
        hist->GetYaxis()->SetTitleFont(textFont);
 
        hist->GetXaxis()->SetTitle(title);
        hist->GetXaxis()->SetTitleOffset(1);
        hist->GetYaxis()->SetTitle("Entries");
        hist->GetYaxis()->SetTitleOffset(1.05);   // good then entries per bit <= 9999
    } else {
        std::cout << " E: Read hists error! " << std::endl;
    }
}
 
#endif