filtrator.C

Go to the documentation of this file.

#include "TCanvas.h"
#include "TF1.h"
#include "TH1.h"
#include "TLegend.h"
#include "TMath.h"
 
#include <TFile.h>
#include <TKey.h>
#include <TStyle.h>
 
using namespace std;
using namespace TMath;
 
// signal region
//  Float_t leftSig = 2.980;
//  Float_t rightSig = 3.002;
 
const Float_t massTheor = 2.992;
 
// Float_t range = 0.012;
// Float_t leftSig = massTheor - range;
// Float_t rightSig = massTheor + range;
// full region
const Float_t minBound = 2.96;
const Float_t maxBound = 3.02;
 
Double_t backgroundChe(Double_t* xx, Double_t* p)
{
    Double_t fA = minBound;
    Double_t fB = maxBound;
 
    vector<Double_t> fT;
    fT.resize(6);   // polynomial
 
    Double_t x = (2.0 * xx[0] - fA - fB) / (fB - fA);
    int order = fT.size();
    if (order == 1)
        return p[0];
    if (order == 2)
        return p[0] + x * p[1];
    // build the polynomials
    fT[0] = 1;
    fT[1] = x;
    for (Int_t i = 1; i < order; ++i) {
        fT[i + 1] = 2 * x * fT[i] - fT[i - 1];
    }
    Double_t sum = p[0] * fT[0];
    for (int i = 1; i <= order; ++i) {
        sum += p[i] * fT[i];
    }
    return sum;
}
 
Double_t backgroundChe3(Double_t* xx, Double_t* p)
{
 
    // convert from (A,B) to (-1,1), Chebyshev defined on (-1, 1)
    Double_t x = (2.0 * xx[0] - minBound - maxBound) / (maxBound - minBound);
 
    Double_t T0 = 1.0;
    Double_t T1 = 1.0 * x;
    Double_t T2 = 2.0 * x * x - 1.0;
    Double_t T3 = 4.0 * x * x * x - 3.0 * x;
    Double_t result = p[0] * T0 + p[1] * T1 + p[2] * T2 + p[3] * T3;
 
    return result;
}
 
Double_t backgroundChe4(Double_t* xx, Double_t* p)
{
 
    // convert from (A,B) to (-1,1), Chebyshev defined on (-1, 1)
    Double_t x = (2.0 * xx[0] - minBound - maxBound) / (maxBound - minBound);
 
    Double_t T0 = 1.0;
    Double_t T1 = 1.0 * x;
    Double_t T2 = 2.0 * x * x - 1.0;
    Double_t T3 = 4.0 * x * x * x - 3.0 * x;
    Double_t T4 = 8.0 * x * x * x * x - 8.0 * x * x + 1.0;
    Double_t result = p[0] * T0 + p[1] * T1 + p[2] * T2 + p[3] * T3 + p[4] * T4;
 
    return result;
}
 
Double_t backgroundChe6(Double_t* xx, Double_t* p)
{
 
    // convert from (A,B) to (-1,1), Chebyshev defined on (-1, 1)
    Double_t x = (2.0 * xx[0] - minBound - maxBound) / (maxBound - minBound);
 
    Double_t T0 = 1.0;
    Double_t T1 = 1.0 * x;
    Double_t T2 = 2.0 * x * x - 1.0;
    Double_t T3 = 4.0 * x * x * x - 3.0 * x;
    Double_t T4 = 8.0 * x * x * x * x - 8.0 * x * x + 1.0;
    Double_t T5 = 16.0 * x * x * x * x * x - 20.0 * x * x * x + 5.0 * x;
    Double_t T6 = 32.0 * x * x * x * x * x * x - 48.0 * x * x * x * x + 18.0 * x * x - 1.0;
    Double_t result = p[0] * T0 + p[1] * T1 + p[2] * T2 + p[3] * T3 + p[4] * T4 + p[5] * T5 + p[6] * T6;
 
    return result;
}
 
// Gausian Peak function
Double_t gausPeak(Double_t* x, Double_t* par)
{
    Double_t arg = (x[0] - par[1]) / par[2];
    Double_t fitval = par[0] * TMath::Exp(-0.5 * arg * arg);
    return fitval;
}
 
// Sum of background and peak function
Double_t fitFunction(Double_t* x, Double_t* par)
{
    return gausPeak(x, &par[0]) + backgroundChe4(x, &par[3]);
}
 
void DrawNicePlot(TH1F* h, Int_t S, Float_t B)
{
    gStyle->SetOptStat(kFALSE);
    auto c = new TCanvas(h->GetName(), h->GetName(), 800, 500);
    h->SetTitle("");
    h->GetXaxis()->SetTitle("M_{^{3}He + #pi^{-}}, GeV/c^{2}");
    h->GetXaxis()->SetTitleSize(0.045);
    h->GetXaxis()->SetTitleOffset(0.9);
    h->GetXaxis()->SetTitleFont(62);
    h->GetXaxis()->SetLabelFont(62);
 
    h->GetYaxis()->SetTitle(
        Form("Entries / (%d MeV/c^{2})", (int)(1000 * (3.05 - 2.95) / h->GetXaxis()->GetNbins() + 0.5)));
    h->GetYaxis()->SetTitleSize(0.045);
    h->GetYaxis()->SetTitleOffset(1.1);
    h->GetYaxis()->SetTitleFont(62);
    h->GetYaxis()->SetLabelFont(62);
 
    h->SetMarkerStyle(20);
    h->SetMarkerSize(1.0);
    h->SetMarkerColor(kBlack);
    h->SetLineColor(kBlack);
    h->SetMinimum(0);
 
    TF1* fitS = h->GetFunction("fitFcn");
    if (!fitS)
        return;
    fitS->SetNpx(500);
    fitS->SetLineColor(kBlack);
    TF1* fitB = h->GetFunction("fitBack");
    if (!fitB)
        return;
    fitB->SetLineColor(kBlack);
    fitB->SetLineStyle(7);
    auto leg = new TLegend(0.2, 0.2, 0.6, 0.5);
    // auto leg = new TLegend(0.5, 0.5, 0.9, 0.9);
 
    Float_t T = S + B;
 
    Float_t dS = TMath::Sqrt(T + B);
    Float_t dSB = TMath::Sqrt((T * B + 1) / B / B / B);
    Float_t dSIGN = 0.5 * TMath::Sqrt(1 + 12 * T * B / (T + B) / (T + B));
 
    leg->SetHeader("Parameters:", "L");
    leg->AddEntry((TObject*)0, Form("S = %d #pm %d", S, Int_t(dS)), "");
    leg->AddEntry((TObject*)0, Form("S / B = %4.3f #pm %4.3f", S / B, dSB), "");
    leg->AddEntry((TObject*)0, Form("S / #sqrt{S + 2B} = %4.1f #pm %4.1f", S / dS, dSIGN), "");
    leg->AddEntry((TObject*)0, Form("#mu = %4.4f #pm %4.4f", fitS->GetParameter(1), fitS->GetParError(1)), "");
    leg->AddEntry((TObject*)0,
                  Form("#sigma = %4.4f #pm %4.4f", TMath::Abs(fitS->GetParameter(2)), fitS->GetParError(2)), "");
    h->Draw("pe1");
    leg->Draw();
 
    // TH1F* hSignal = new TH1F("hSignal", "hSignal", h->GetNbinsX(), h->GetXaxis()->GetXmin(),
    // h->GetXaxis()->GetXmax()); hSignal->SetMarkerStyle(29); hSignal->SetMarkerSize(1.5);
    // hSignal->SetMarkerColor(kBlack);
    // hSignal->SetLineColor(kBlack);
    // for (Int_t i = 0; i < h->GetNbinsX(); ++i) {
    //     hSignal->SetBinContent(i, h->GetBinContent(i) - fitB->Eval(h->GetBinCenter(i)));
    // }
    // hSignal->Draw("pe1 same");
    // c->SaveAs(Form("%s.png", h->GetName()));
}
 
void fitator(TH1F* h, Double_t& Tout, Double_t& Bout, Double_t& sigmaOut)
{
 
    if (h->GetEntries() == 0)
        return;
 
    Int_t leftBin = 9;     // h->FindBin(leftSig);
    Int_t rightBin = 13;   // h->FindBin(rightSig);
 
    vector<Float_t> values;
    vector<Float_t> errors;
    for (Int_t i = leftBin; i < rightBin; i++) {
        values.push_back(h->GetBinContent(i));
        errors.push_back(h->GetBinError(i));
        h->SetBinContent(i, 0);
        h->SetBinError(i, 1e6);
    }
    TF1* fitBack = new TF1("fitBack", backgroundChe4, minBound, maxBound, 5);
    // TF1* fitBack = new TF1("fitBack", backgroundChe4, minBound, maxBound, 5);
    // TF1* fitBack = new TF1("fitBack", backgroundChe3, minBound, maxBound, 4);
    h->Fit(fitBack, "QIER", "", minBound, maxBound);
 
    Double_t B = 0;
    Double_t T = 0;
 
    for (Int_t i = leftBin; i < rightBin; i++) {
        h->SetBinContent(i, values[i - leftBin]);
        h->SetBinError(i, errors[i - leftBin]);
        T += values[i - leftBin];
        B += fitBack->Eval(h->GetBinCenter(i));
    }
 
    // if (T - B < 1000) return;
    // cout << T << "   " << B << endl;
 
    TF1* fitFcn = new TF1("fitFcn", fitFunction, minBound, maxBound, 8);
    fitFcn->SetLineColor(kBlue);
 
    // Устанавливаем параметры гаусса
    fitFcn->SetParameter(0, (T - B) / 2);
    fitFcn->SetParameter(1, massTheor);
    fitFcn->SetParameter(2, 0.004);   // Ожидаемая ширина
    // Ограничения для стабильности
    fitFcn->SetParLimits(0, 0.1 * (T - B), 1.5 * (T - B));   // Амплитуда не может быть отрицательной
    fitFcn->SetParLimits(1, massTheor - 0.010, massTheor + 0.010);
    fitFcn->SetParLimits(2, 0.003, 0.010);
    fitFcn->FixParameter(3, fitBack->GetParameter(0));
    fitFcn->FixParameter(4, fitBack->GetParameter(1));
    fitFcn->FixParameter(5, fitBack->GetParameter(2));
    fitFcn->FixParameter(6, fitBack->GetParameter(3));
    fitFcn->FixParameter(7, fitBack->GetParameter(4));
    // fitFcn->FixParameter(7, fitBack->GetParameter(4));
    // fitFcn->FixParameter(8, fitBack->GetParameter(5));
    // fitFcn->FixParameter(9, fitBack->GetParameter(6));
 
    h->Fit(fitFcn, "QIER+", "e1p", minBound, maxBound);
    // cout << h->GetName() << endl;
    Tout = T;
    Bout = B;
    sigmaOut = fitFcn->GetParameter(2);
}
 
void filtrator(TString inFileName = "/home/merz/data/h3l/h3l_plots_25090_A_withKplus/h3l_plots.root")
{
 
    TFile* f = TFile::Open(inFileName);
    TIter keyList(f->GetListOfKeys());
    TKey* key;
 
    Int_t nFiltered = 0;
 
    cout << "Number of total hist = " << f->GetListOfKeys()->GetSize() << endl;
 
    Float_t maxS = 0.0;
    Float_t maxS_B = 0.0;
    Float_t maxS_S = 0.0;
    Float_t maxSB = 0.0;
    Float_t maxSB_B = 0.0;
    Float_t maxSB_S = 0.0;
    Float_t maxSig = 0.0;
    Float_t maxSig_B = 0.0;
    Float_t maxSig_S = 0.0;
    TH1F* hMaxS = nullptr;
    TH1F* hMaxSB = nullptr;
    TH1F* hMaxSig = nullptr;
 
    TFile* fOut = new TFile("/home/merz/data/h3l/h3l_plots_25090_A_withKplus/h3l_plots_filtered.root", "recreate");
 
    vector<Float_t> dSignArray;
    vector<Float_t> dSArray;
    vector<Float_t> dBArray;
    vector<Float_t> signArray;
    vector<Float_t> sArray;
    vector<Float_t> bArray;
    while ((key = (TKey*)keyList())) {
        TH1F* h = (TH1F*)key->ReadObj();
        TString name = h->GetName();
 
        h->Rebin(2);
        h->Rebin(2);
 
        Double_t T = 0;
        Double_t B = 0;
        Double_t sigma = 0;
        fitator(h, T, B, sigma);
 
        if (T == 0 || B == 0 || sigma == 0)
            continue;
 
        Float_t S = T - B;
        Float_t SB = S / B;
        Float_t SIGN = S / TMath::Sqrt(T + B);
 
        h->SetTitle(Form("S:%4.2f_B:%4.2f_SB:%4.2f_SIGN:%4.2f", S, B, SB, SIGN));
 
        if (S > maxS) {
            maxS = S;
            maxS_B = B;
            maxS_S = S;
            hMaxS = h;
            hMaxS->SetTitle(Form("maxS___S:%4.2f_B:%4.2f_SB:%4.2f_SIGN:%4.2f", S, B, SB, SIGN));
        }
 
        if (SB > maxSB) {
            maxSB = SB;
            maxSB_B = B;
            maxSB_S = S;
            hMaxSB = h;
            hMaxSB->SetTitle(Form("maxSB___S:%4.2f_B:%4.2f_SB:%4.2f_SIGN:%4.2f", S, B, SB, SIGN));
        }
 
        if (SIGN > maxSig) {
            maxSig = SIGN;
            maxSig_B = B;
            maxSig_S = S;
            hMaxSig = h;
            hMaxSig->SetTitle(Form("maxSig___S:%4.2f_B:%4.2f_SB:%4.2f_SIGN:%4.2f", S, B, SB, SIGN));
        }
 
        // if (S < 2000)
        // continue;
        // if (SB < 0.035)
        //     continue;
        if (SIGN < 3.9)
            continue;
        // if (sigma > 0.0050)
        //     continue;
        Float_t dS = TMath::Sqrt(T + B);
        Float_t dB = TMath::Sqrt(B);
        Float_t dSIGN = 0.5 * TMath::Sqrt(1 + 12 * T * B / (T + B) / (T + B));
        // cout << h->GetName() << ":\t" << sigma << "\t" << S << "\t" << dS << "\t" << SIGN << "\t" << dSIGN << "\t" <<
        // B << "\t" << dB << endl;
        cout << h->GetName() << ":\t" << sigma << "\t" << S << "\t" << SIGN << "\t" << SB << endl;
        signArray.push_back(SIGN);
        sArray.push_back(S);
        bArray.push_back(B);
        dSignArray.push_back(dSIGN);
        dBArray.push_back(dB);
        dSArray.push_back(dS);
        h->Write();
        nFiltered++;
    }
 
    cout << endl;
    for (auto s : signArray) {
        printf("%4.1f, ", s);
    }
    cout << endl;
    for (auto s : dSignArray) {
        printf("%4.1f, ", s);
    }
    cout << endl;
    for (auto s : sArray) {
        printf("%4.0f, ", s);
    }
    cout << endl;
    for (auto s : dSArray) {
        printf("%4.0f, ", s);
    }
    cout << endl;
    for (auto s : bArray) {
        printf("%4.0f, ", s);
    }
    cout << endl;
    for (auto s : dBArray) {
        printf("%4.0f, ", s);
    }
    cout << endl;
    cout << endl;
 
    DrawNicePlot(hMaxS, maxS_S, maxS_B);
    DrawNicePlot(hMaxSB, maxSB_S, maxSB_B);
    DrawNicePlot(hMaxSig, maxSig_S, maxSig_B);
 
    cout << hMaxS->GetName() << "\t" << hMaxSB->GetName() << "\t" << hMaxSig->GetName() << endl;
 
    cout << "Number of filtered histo = " << nFiltered << endl;
 
    fOut->Close();
}