EvalQA.C

Go to the documentation of this file.

#include <Rtypes.h>
#include <TCanvas.h>
#include <TChain.h>
#include <TClonesArray.h>
#include <TF1.h>
#include <TH1.h>
#include <TH2.h>
#include <TMath.h>
#include <TPad.h>
#include <TString.h>
#include <TVirtualPad.h>
R__ADD_INCLUDE_PATH($VMCWORKDIR)
#include "macro/run/bmnloadlibs.C"
 
using namespace std;
using namespace TMath;
 
const Int_t dim = 2;
// Redefine histo limits (supposing exclusion - less(idx0) or more(idx1)) due to cuts obtained (useCutHistoLimits should
// be kTRUE)
 
// Ramin background function
Double_t background_new(Double_t* x, Double_t* par)
{
    return par[0] * TMath::Sqrt(x[0] - par[2]) * TMath::Exp(-par[1] * (x[0] - par[2]));
}
// signal
Double_t sig_gaus(Double_t* x, Double_t* par)
{
    //                return par[0] * Exp(-0.5 * pow((x[0] - par[1]) / par[2], 2))/(par[2] * sqrt(2 * M_PI)) ;
    return par[0] * Exp(-0.5 * pow((x[0] - par[1]) / par[2], 2)) /*/(par[2] * sqrt(2 * M_PI))*/;
}
 
void EvalQA(TString fileName = "tmp.root",
            Bool_t useCuts = true,
            Bool_t drawFit = true,
            TString outNameInvMass = "im.pdf",
            TString outNameFitRes = "fr.root",
            Double_t thr = 150.0,
            int32_t pdg = 3122)   // Lambda0 3122, K_S0 310
{
    gStyle->SetOptStat(1);
    // -----   Timer   ---------------------------------------------------------
    TStopwatch timer;
    timer.Start();
 
    TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdg);
    Double_t mass = particle->Mass();
    Double_t width = 0.004;                     // particle->Width();
    Double_t rL_Resonance = mass - 3 * width;   // 1.1;
    Double_t rR_Resonance = mass + 3 * width;   // 1.13;
    //    particle->DecayChannel()
    printf("m %f  w %f  l %f  r %f\n", mass, width, rL_Resonance, rR_Resonance);
 
    //    TChain* out = new TChain("bmndata");
    //    out->Add(fileName.Data());
    //    cout << "#recorded entries = " << out->GetEntries() << endl;
    //    if (out->GetEntries() == 0) {
    //        cout << "Something is wrong! Please, check file you passed as input!" << endl;
    //        exit(-1);
    //    }
 
    Double_t rL = 1.08;
    Double_t rR = 1.2;
    //    invMassSpectrum->SetTitle("Invariant mass: #Lambda^{0} #rightarrow p + #pi^{-}");
    unique_ptr<TFile> file{TFile::Open(fileName, "READ")};
    TTree* tree = dynamic_cast<TTree*>(file->Get("event"));
    // TH1F* invMassSpectrum = dynamic_cast<TH1F*>(file->Get("h"));
    TH1F* invMassSpectrum = new TH1F("invMassSpectrum", "invMassSpectrum", 75, rL, rR);
 
    vector<TString> cutStrVec;
    //        cutStrVec.push_back(Form("%.1f < Path < %.1f", Path[0], Path[1]));
    //        cutStrVec.push_back(Form("%.1f < Dca12 < %.1f", Dca12[0], Dca12[1]));
    tree->Draw(
        Form("l0.massh>>%s", invMassSpectrum->GetName()),
        useCuts
            ? "l0.chi2s[][0]>12.5&&l0.chi2s[][1]>10.5&&l0.c2pv<9&&l0.path>2&&l0.tmvas[][0]>-0.03&&l0.tmvas[][1]>-0.03"
            : "");
    if (useCuts)
        printf("Cuts applied: %3.1f/%lld\n", invMassSpectrum->GetEntries(), tree->GetEntries());
 
    Double_t par[8] = {};
    TF1* bg = new TF1("BG", "pol4", rL, rR);
    //    TF1* bg = new TF1("BG", background_new, rL, rR, 5);
    Double_t aver = invMassSpectrum->Integral() / invMassSpectrum->GetNbinsX();   //(rR - rL);
    printf("average %f\n", aver);
    //    par[0] = aver;
    bg->SetParameters(&par[0]);
    invMassSpectrum->Fit(bg, "RS", "", rL, rR);
    bg->GetParameters(&par[0]);
 
    TF1* sig =
        new TF1("SIG", sig_gaus,
                //            [](Double_t *x, Double_t * p) {
                //                return p[0] * Exp(-0.5 * pow((x[0] - p[1]) / p[2], 2))/*/(p[2] * sqrt(2 * M_PI))*/ ;
                //            },
                rL_Resonance, rR_Resonance, 3);
 
    Double_t peak_val = invMassSpectrum->GetBinContent(invMassSpectrum->FindBin((rR_Resonance + rL_Resonance) / 2.0));
    printf("peak_val %f\n", peak_val);
 
    Double_t temp[3] = {peak_val, mass, width};
    sig->SetParameters(temp);
    invMassSpectrum->Fit(sig, "R", "", rL_Resonance, rR_Resonance);
    sig->GetParameters(&par[5]);
 
    //    TF1 *f = new TF1("f", "pol4(0)+gaus(5)", rL, rR);
    TF1* f = new TF1(
        "f",
        [&](Double_t* x, Double_t* p) {
            //        return bg->EvalPar(x, p) + sig->EvalPar(x, p+5);
            return bg->EvalPar(x, p) + sig_gaus(x, &p[5]);
            //        return background_new(x, p) + sig_gaus(x, &p[5]);
        },
        rL, rR, 8);
    f->SetNpx(500);
 
    //    par[0] = aver;
    //    par[5] = peak_val - aver;
    for (int i = 0; i < 8; i++)
        printf("p[%d] = %f\n", i, par[i]);
 
    f->SetParameters(par);
    //    f->SetParLimits(0, 0.0, 1e12);
    //    f->SetParLimits(5, 0.0, 1e12);
    f->SetParLimits(7, 0.0, 1);
    f->SetLineColor(kMagenta + 1);
    f->SetLineWidth(3);
    TString fitOpt = "RS";
    TFitResultPtr fitRes =
        invMassSpectrum->Fit(f, (drawFit) ? fitOpt.Data() : TString(fitOpt + "0").Data(), "", rL, rR);
    Double_t mean = fitRes->Parameter(6);
    Double_t sigma = fitRes->Parameter(7);
    Double_t T = invMassSpectrum->Integral(invMassSpectrum->FindBin(mean - 3 * sigma),
                                           invMassSpectrum->FindBin(mean + 3 * sigma));
    Double_t B = bg->Integral(mean - 3 * sigma, mean + 3 * sigma) / invMassSpectrum->GetBinWidth(1);
    Double_t S = sig->Integral(mean - 3 * sigma, mean + 3 * sigma) / invMassSpectrum->GetBinWidth(1);
    Double_t S_to_B = T / B - 1;
    Double_t Signif = (T - B) / T;
    printf("Tot %f  Bg %f sig %f\n", T, B, S);
    f->GetParameters(par);
    //    delete sig;
    //    sig = new TF1("SIG", sig_gaus, rL_Resonance, rR_Resonance, 3);
    sig->SetParameters(&par[5]);
    bg->SetParameters(&par[0]);
    T = invMassSpectrum->Integral(invMassSpectrum->FindBin(mean - 3 * sigma),
                                  invMassSpectrum->FindBin(mean + 3 * sigma));
    B = bg->Integral(mean - 3 * sigma, mean + 3 * sigma) / invMassSpectrum->GetBinWidth(1);
    S = sig->Integral(mean - 3 * sigma, mean + 3 * sigma) / invMassSpectrum->GetBinWidth(1);
    S_to_B = T / B - 1;
    Signif = (T - B) / T;
    printf("Tot %f  Bg %f sig %f after refit\n", T, B, S);
    invMassSpectrum->SetMarkerStyle(20);
    invMassSpectrum->SetMarkerColor(kSpring - 6);
    invMassSpectrum->SetLineColor(kSpring - 6);
    invMassSpectrum->SetLineWidth(1);
 
    TFile out(outNameFitRes, "RECREATE");
 
    unique_ptr<TCanvas> canvLambda = make_unique<TCanvas>(Form("im_thr%1.1f", thr), "Inv mass", 1600, 900);
    canvLambda->Divide(1, 1);
    TVirtualPad* pad = canvLambda->cd();
    pad->SetLeftMargin(0.2);
    canvLambda->cd(1);
 
    invMassSpectrum->Draw("PE1X0");
    // invMassSpectrum->GetYaxis()->SetRangeUser(0., 40.);
    invMassSpectrum->GetXaxis()->SetTitle("M_{(p + #pi^{-})}, GeV/c^{2}");
    invMassSpectrum->GetYaxis()->SetTitle("Entries / 2 MeV/c^{2}");
    invMassSpectrum->GetXaxis()->SetLabelSize(0.035);
    invMassSpectrum->GetYaxis()->SetLabelSize(0.035);
    invMassSpectrum->GetXaxis()->SetTitleSize(0.05);
    invMassSpectrum->GetYaxis()->SetTitleSize(0.05);
    invMassSpectrum->GetXaxis()->SetTitleOffset(0.85);
    invMassSpectrum->GetYaxis()->SetTitleOffset(1.5);
    invMassSpectrum->GetYaxis()->CenterTitle();
    invMassSpectrum->SetMinimum(0);
 
    // Draw the legend
    Double_t legendX = 0.55;
    Double_t legendY = 0.85;
    TLegend* legend = new TLegend(legendX, legendY - (4 + cutStrVec.size()) * 0.05, legendX + 0.05, legendY);
    legend->SetTextFont(72);
    legend->SetTextSize(0.03);
    legend->AddEntry("", Form("Mass = %.4f", mean), 0);
    legend->AddEntry("", Form("  #sigma    = %.4f", sigma), 0);
    legend->AddEntry("", Form("S/B = %.4f", S_to_B), 0);
    legend->AddEntry("", Form("S/(S + B) = %.4f", Signif), 0);
    legend->AddEntry("", Form("Sig integral = %.4f", S), 0);
    if (cutStrVec.size() > 0)
        legend->AddEntry("", "", 0);
    for (TString& cutS : cutStrVec) {
        legend->AddEntry("", cutS.Data(), 0);
    }
    legend->SetLineColor(0);
    //    if (useCuts){
    legend->Draw();
    sig->SetLineColor(kBlue);
    sig->Draw("SAME");
    //    delete sig;
    //    delete bg;
    //    }
    // -----   Finish   --------------------------------------------------------
    canvLambda->SaveAs(outNameInvMass.Data());
    canvLambda->Write();
 
    TTree* tree_res = new TTree("FitRes", "FitRes");
    tree_res->Branch("Mean", &mean);
    tree_res->Branch("S_to_B", &S_to_B);
    tree_res->Branch("Signif", &Signif);
    tree_res->Branch("Signal", &S);
    tree_res->Branch("Thr", &thr);
    tree_res->Branch("Tot", &T);
    tree_res->Branch("Bg", &B);
    tree_res->Fill();
    out.Write();
    out.Close();
    printf("Results written to %s\n", outNameFitRes.Data());
 
    timer.Stop();
    Double_t rtime = timer.RealTime();
    Double_t ctime = timer.CpuTime();
    cout << endl << endl;
    cout << "Macro finished successfully." << endl;
    cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
    cout << endl;
}