resonance_continuum_fit_nophase.C

#include <TMinuit.h>
#include <TVirtualFitter.h>
#include <TGraph.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TStyle.h>
#include <TF1.h>
#include <TF1Convolution.h>
#include <TLegend.h>
#include <iostream>
#include <cmath>
#include <TROOT.h>
#include <TMath.h>
#include <TColor.h>
#include <TH1F.h>


TString pName = "fit_nophase";
double Gamma( double Mkk, double Mphi, double Gamma0 ){
    double Mkk2 = Mkk*Mkk;
    double Mphi2 = Mphi*Mphi;
    double mk = 495.2;
    double mk2 = mk*mk;
    double t0 = Gamma0 * ( Mphi / Mkk );
    double t1 = pow( (Mkk2 - 4*mk2)/(Mphi2-4*mk2), (3.0/2.0) );
    return t0 * t1;
}

// Relativistic Breit-Wigner (amplitude level) for the phi
// Real or Imaginary part of the amplitude
double Amplitude_BW( double A, double Mkk, double Mphi, double Gamma0, bool real = true ){
    double Mkk2 = Mkk*Mkk;
    double Mphi2 = Mphi*Mphi;
    double G = Gamma( Mkk, Mphi, Gamma0 );
    double G2 = G*G;

    double pre = A * sqrt(Mkk * Mphi * G);
    
    double t0 = Mkk2 - Mphi2;
    double t1 = pow( Mkk2 - Mphi2, 2 ) + Mphi2*G2;

    if (real){
        return pre * t0 / t1;
    } else {
        return pre * (-Mphi*G) / t1;
    }
}

double Amplitude_BW( double *x, double *p ){
    double Mkk = x[0]; // x-variable
    double Mphi = p[0];  // mass of the phi
    double Gamma0 = p[1]; // width of the phi
    double A = p[2]; // alpha
    double Br = p[3]; // beta (Real)
    double Bi = p[4]; // beta (Imaginary)
    double RI = p[5]; // real or imaginary
    return Amplitude_BW( A, Mkk, Mphi, Gamma0, (int)(RI) );
}


double ResidualBackground( double *x, double *p ){
    
    double Mkk = x[0]; // x-variable
    double Bg0 = p[0]; // background
    double Bg1 = p[1]; // background
    double Bg2 = p[2]; // background
    // double Bg3 = p[3]; // background
    // double Bg4 = p[4]; // background


    return Bg0 / ( 1 + exp( -Bg1 * (Mkk - Bg2) ) );

    // return fabs(Bg0) + fabs(Bg1 * Mkk) + fabs(Bg2 * Mkk * Mkk);
    // return fabs(Bg0 + (Bg1 * Mkk) + Bg2 * Mkk * Mkk);
}

double xs( double *x, double *p ){
    double Mkk = x[0]; // x-variable
    // fit parameters
    double Mphi = p[0];  // mass of the phi
    double Gamma0 = p[1]; // width of the phi
    double alpha = p[2]; // alpha
    double betaR = p[3]; // beta (Real)
    double betaI = p[4]; // beta (Imaginary)

    double ReA = Amplitude_BW( alpha, Mkk, Mphi, Gamma0, true );
    double ImA = Amplitude_BW( alpha, Mkk, Mphi, Gamma0, false );
    
    // Let the continuum be complex B = betaR + i*betaI
    double ReB = betaR;
    double ImB = betaI;

    // compute each (real) term in the matrix element squared
    // (ReA + ImA) * (ReA + ImA) = ReA*ReA + ImA*ImA + 2*ReA*ImA
    // (ReB + ImB) * (ReB + ImB) = ReB*ReB + ImB*ImB + 2*ReB*ImB
    // (ReA + ImA) * (ReB + ImB) = ReA*ReB + ImA*ImB + ReA*ImB + ReB*ImA
    double t0 = ReA*ReA + ImA*ImA;
    double t1 = ReA*ReB + ImA*ImB;
    double t2 = 2*(ReB*ReB + ImB*ImB);

    return t0 + t1 + t2 + ResidualBackground( x, &(p[5]) );
}


double xs_( double *x, double *p ){
    double Mkk = x[0]; // x-variable
    // fit parameters
    double Mphi = p[0];  // mass of the phi
    double Gamma0 = p[1]; // width of the phi
    double alpha = p[2]; // alpha
    double betaR = p[3]; // beta (Real)
    double betaI = p[4]; // beta (Imaginary)
    double term0 = p[5]; // background

    double ReA = Amplitude_BW( alpha, Mkk, Mphi, Gamma0, true );
    double ImA = Amplitude_BW( alpha, Mkk, Mphi, Gamma0, false );
    
    // Let the continuum be complex B = betaR + i*betaI
    double ReB = betaR;
    double ImB = betaI;

    // compute each (real) term in the matrix element squared
    // (ReA + ImA) * (ReA + ImA) = ReA*ReA + ImA*ImA + 2*ReA*ImA
    // (ReB + ImB) * (ReB + ImB) = ReB*ReB + ImB*ImB + 2*ReB*ImB
    // (ReA + ImA) * (ReB + ImB) = ReA*ReB + ImA*ImB + ReA*ImB + ReB*ImA
    double t0 = ReA*ReA + ImA*ImA;
    double t1 = ReA*ReB + ImA*ImB;
    double t2 = 2*(ReB*ReB + ImB*ImB);

    if (term0< -1)
        return t0 ;
    else if ( fabs(term0) < 1 )
        return t1;
    else
        return t2;
    
}

TGraph* contour_for(int nSigma = 1, int color = kBlue) {
    
    gMinuit->SetErrorDef(nSigma * nSigma); // Set the error definition for the contour

    TGraph *contour_BrBi = (TGraph*)gMinuit->Contour(100, 3, 4); 
    TGraph *contour_ABr  = (TGraph*)gMinuit->Contour(100, 2, 3); 

    if (!contour_BrBi || !contour_ABr) {
      printf("Error: Contour could not be generated.\n");
      return nullptr;
    }

    int nPoints = contour_BrBi->GetN();
    double *Br_vals = contour_BrBi->GetX();
    double *Bi_vals = contour_BrBi->GetY();
    double *A_vals = contour_ABr->GetX();
    double BoA_vals[10000]; // Array for B over A values
    double delta_vals[10000]; // Array for delta values

    for ( int i = 0; i < nPoints; i++ ){
        double Br = Br_vals[i];
        double Bi = Bi_vals[i];
        double A  = A_vals[i];

        // Calculate B over A and delta
        BoA_vals[i]   = sqrt( pow( Br, 2 ) + pow( Bi, 2 ) ) / A;
        delta_vals[i] = atan2( Bi, Br );
        // cout << "Br = " << Br << ", Bi = " << Bi << ", A = " << A << endl;
    }

    TGraph *contour_ALICE = new TGraph(nPoints, BoA_vals, delta_vals);
    contour_ALICE->SetTitle(TString::Format("%d-sigma Contour;|B/A|;#delta", nSigma));
    contour_ALICE->SetLineColor(color);
    contour_ALICE->SetLineColor(color);
    contour_ALICE->SetLineWidth(2);
    contour_ALICE->SetMarkerStyle(20);
    contour_ALICE->SetMarkerSize(0.001);
    return contour_ALICE;
}



void fit_histo( TH1 * hm){

    TCanvas *c1 = new TCanvas("c1", "c1", 1600, 900);
    c1->SetTopMargin(0.05);
    c1->SetRightMargin(0.05);
    c1->SetLeftMargin(0.15);
    c1->SetBottomMargin(0.15);
    c1->cd();


    hm->GetXaxis()->SetRangeUser( 990, 1450 );
    hm->SetFillColorAlpha( kBlue, 0.35 );
    hm->SetLineColor( kBlue );
    hm->GetYaxis()->SetTitle("Events / 10 MeV");
    hm->GetXaxis()->SetTitle("M_{KK} [MeV]");
    hm->SetTitle("M_{KK} Spectrum");
    hm->GetYaxis()->SetRangeUser( 10, hm->GetMaximum()*1.5 );
    hm->Draw( "pe" ) ;
    hm->Draw( "same h" ) ;

    // draw the candidates and background as well
    // hc->SetFillColorAlpha( kRed, 0.35 );
    // hc->SetLineColor( kRed );
    // hc->Draw( "same pe" ) ;
    // hbg->SetFillColorAlpha( kGreen, 0.35 );
    // hbg->SetLineColor( kGreen );
    // hbg->Draw( "same pe" ) ;

    gStyle->SetFitFormat( "5.3g" );
    gStyle->SetOptFit( 1111 );
    // hm->Draw();

     int colors[6] = {kRed, kBlue, kGreen+2, kMagenta, kOrange+7, kCyan+2};
    
    TF1 * tf1kkKernel = new TF1("xs", xs, 1000, 1400, 8);
    tf1kkKernel->SetNpx(1000);
    tf1kkKernel->SetLineColor( colors[0] );
    tf1kkKernel->SetLineWidth( 5 );
    tf1kkKernel->SetParNames( "Mphi", "Gamma", "A", "Br", "Bi", "Bg0", "Bg1", "Bg2" );
    tf1kkKernel->SetParameters( 1020, 5.09, 63, 1, 0, 100, 1.0, 1000  ); // setting initial guess for each parameter
    hm->Fit(tf1kkKernel, "MER", "", 1005, 1100);

    TF1 *gaus = new TF1("gaus", "exp(-0.5*(x/[0])*(x/[0]))", -100, 100);
    gaus->SetParameter(0, 1.0); // Initial sigma, e.g., 1 GeV, adjust as needed
    TF1Convolution *conv = new TF1Convolution(tf1kkKernel, gaus, true);
    conv->SetNofPointsFFT(10);
    // TF1 *tf1kk = new TF1("cxs", conv, 1000, 1200, tf1kkKernel->GetNpar() + gaus->GetNpar());
    TF1 * tf1kk = new TF1("xs", xs, 1000, 1400, 8);
    tf1kk->SetParameters( tf1kkKernel->GetParameters() );
    tf1kk->SetParNames( "Mphi", "Gamma", "A", "Br", "Bi", "Bg0", "Bg1", "Bg2", "RES" );
    // tf1kk->FixParameter( 8, 1.0 ); // Fix the resolution parameter
    // tf1kk->FixParameter( 6, 0 );
    // tf1kk->FixParameter( 7, 0 );
    // hm->Fit(tf1kk, "R", "", 1000, 1040);
    // hm->Fit(tf1kk, "R", "", 1000, 1200);

    // tf1kk->Draw("");
    // c1->Print(TString::Format("%s.png", pName.Data()).Data());
    // return;

    // Draw the residual background separately
    TF1 * tf1bg = new TF1("bg", ResidualBackground, 1000, 1400, 3);
    tf1bg->SetNpx(1000);
    tf1bg->SetLineWidth( 5 );
    tf1bg->SetParameters( tf1kk->GetParameter(5), tf1kk->GetParameter(6), tf1kk->GetParameter(7) );
    tf1bg->SetLineColor( colors[1] );
    tf1bg->Draw("same");

    gPad->SetLogy();

    
    double Mphi = tf1kk->GetParameter(0);
    double Gamma = tf1kk->GetParameter(1);
    double A = tf1kk->GetParameter(2);
    double Br = tf1kk->GetParameter(3);
    double Bi = tf1kk->GetParameter(4);
    double Bg0 = tf1kk->GetParameter(5);
    double Bg1 = tf1kk->GetParameter(6);
    double Bg2 = tf1kk->GetParameter(7);

    // lets plot each part of the amplitude separately
    TF1 * tf1xs0 = new TF1("xs0", xs_, 1000, 1400, 6);
    tf1xs0->SetNpx(1000);
    tf1xs0->SetLineWidth( 5 );
    tf1xs0->SetParameters( Mphi, Gamma, A, Br, Bi, -10 ); // first term
    tf1xs0->SetLineColor( colors[2] );
    tf1xs0->Draw("same");

    TF1 * tf1xs1 = new TF1("xs1", xs_, 1000, 1400, 6);
    tf1xs1->SetNpx(1000);
    tf1xs1->SetLineWidth( 5 );
    tf1xs1->SetParameters( Mphi, Gamma, A, Br, Bi, 0 ); // second term
    tf1xs1->SetLineColor( colors[3] );
    tf1xs1->Draw("same");

    TF1 * tf1xs2 = new TF1("xs2", xs_, 1000, 1400, 6);
    tf1xs2->SetNpx(1000);
    tf1xs2->SetLineWidth( 5 );
    tf1xs2->SetParameters( Mphi, Gamma, A, Br, Bi, 10 ); // third term
    tf1xs2->SetLineColor( colors[4] );
    tf1xs2->Draw("same");

    TLegend *leg = new TLegend(0.45, 0.75, 0.6, 0.95);
    leg->AddEntry( tf1kk, "Total Amplitude", "l" );
    leg->AddEntry( tf1xs0, "First Term", "l" );
    leg->AddEntry( tf1xs1, "Second Term", "l" );
    leg->AddEntry( tf1xs2, "Third Term", "l" );
    leg->AddEntry( tf1bg, "Residual Background", "l" );
    leg->Draw("same");

    c1->Print(TString::Format("%s.png", pName.Data()).Data());


    // make a new canvas
    TCanvas *c2 = new TCanvas("c2", "c2", 1600, 900);
    c2->SetTopMargin(0.05);
    c2->SetRightMargin(0.05);
    c2->SetLeftMargin(0.15);
    c2->SetBottomMargin(0.15);
    c2->cd();

    // lets make a pull plot between the tf1kk and the data
    TH1F *hpull = (TH1F*)hm->Clone("hpull");
    hpull->SetTitle("Pull Plot");

    for (int i = 1; i <= hm->GetNbinsX(); i++){
        double x = hm->GetBinCenter(i);
        double y = hm->GetBinContent(i);
        double yerr = hm->GetBinError(i);
        double yfit = tf1kk->Eval(x);
        double pull = (y - yfit) / yerr;
        if (pull!=pull || y < 10) continue;
        cout << "pull[" << i << "]: " << pull << endl;
        hpull->SetBinContent(i, pull);
        hpull->SetBinError(i, 1.0);
    }
    hpull->SetMarkerStyle(20);
    hpull->SetMarkerSize(1.5);
    hpull->SetLineWidth(5);
    hpull->SetLineColor( kBlue );
    hpull->GetYaxis()->SetRangeUser( -5, 5 );
    hpull->Draw( "pe");

    c1->Print(TString::Format("%s_pull.png", pName.Data()).Data());

    return;

    // new canvas
    TCanvas *c3 = new TCanvas("c3", "c3", 1600, 900);
    c3->SetTopMargin(0.05);
    c3->SetRightMargin(0.05);
    c3->SetLeftMargin(0.15);
    c3->SetBottomMargin(0.15);
    c3->cd();

    // TGraph *contour = (TGraph*)fitResult->GetConfidenceEllipse(0, 1, 50, 1); // Parameters 0, 1, 50 points, 1-sigma
    // contour->SetTitle("1-sigma Contour;Parameter 0;Parameter 1");
    // contour->Draw("AL");
    // TGraph *gr12 = (TGraph*)gMinuit->Contour(40,3,4);
    // gr12->Draw("alp");
    // c3->Draw();

    // c3->Print("fit_nophase_contour.png");

    double BoA = sqrt( 
        pow( tf1kk->GetParameter(3), 2 ) + pow( tf1kk->GetParameter(4), 2 )
    ) / tf1kk->GetParameter(2);
    double BoAerr = BoA * sqrt( 
        pow( tf1kk->GetParError(3) / tf1kk->GetParameter(3), 2 ) +
        pow( tf1kk->GetParError(4) / tf1kk->GetParameter(4), 2 ) +
        pow( tf1kk->GetParError(2) / tf1kk->GetParameter(2), 2 )
    );
    cout << "BoA = " << BoA << " +/- " << BoAerr << endl;

    double delta = atan2( tf1kk->GetParameter(4), tf1kk->GetParameter(3) );
    double deltaerr = sqrt( 
        pow( tf1kk->GetParError(4) / tf1kk->GetParameter(4), 2 ) +
        pow( tf1kk->GetParError(3) / tf1kk->GetParameter(3), 2 )
    );
    cout << "delta = " << delta << " +/- " << deltaerr << endl;
    return;

    // now we compute the contours
    TGraph * contour_ALICE3 = contour_for( /*nSigma=*/3, kRed );
    contour_ALICE3->SetTitle("Phase-free fit;|B/A|;#delta");
    contour_ALICE3->Draw("ALP");

    TGraph * contour_ALICE1 = contour_for( /*nSigma=*/1, kBlue );
    contour_ALICE1->Draw("same LP");


    TGraph * best_fit = new TGraph(1);
    best_fit->SetPoint(0, BoA, delta);
    best_fit->SetMarkerStyle(20);
    best_fit->SetMarkerSize(2);
    best_fit->SetMarkerColor(kBlack);
    best_fit->Draw("same P");

    TLegend *leg2 = new TLegend(0.15, 0.15, 0.35, 0.35);
    leg2->AddEntry( contour_ALICE1, "1-sigma Contour", "l" );
    leg2->AddEntry( contour_ALICE3, "3-sigma Contour", "l" );
    leg2->AddEntry( best_fit, "Best Fit", "p" );
    leg2->Draw("same");

    

    c1->Print(TString::Format("%s_contour.png", pName.Data()).Data());
    
    
    
    TFile *fout = new TFile("fit_nophase.root", "RECREATE");
    contour_ALICE1->SetName("contour_nophase1");
    contour_ALICE1->Write();
    contour_ALICE3->SetName("contour_nophase3");
    contour_ALICE3->Write();
    best_fit->SetName("central_nophase");
    best_fit->Write();
}

void fit_nophase( TString n = "FILTERED.root" ){

    //be sure default is Minuit since we will use gMinuit
   TVirtualFitter::SetDefaultFitter("Minuit");

    // open a root TFile for reading
    TFile *f = new TFile(TString::Format("%s", n.Data()).Data());

    auto hm0 = (TH1*)f->Get("hY"); // candidates
    fit_histo(hm0);

    for (int i = 0; i < 10; i++){
        cout << "i = " << i << endl;
        auto hm = (TH1*)f->Get(TString::Format("hY_%d", i)); // candidates
        pName = TString::Format("fit_nophase_%d", i);
        fit_histo(hm);
    }
    

}