t-distribution fit, simple WS

Makefile

CXX=g++
RM=rm -f
CPPFLAGS=-g $(shell root-config --cflags)
LDFLAGS=-g $(shell root-config --ldflags)
LDLIBS=$(shell root-config --libs --glibs)

SRCS=wsfit.cxx
OBJS=$(subst .cxx,.o,$(SRCS))

all: wsfit

wsfit: $(OBJS)
	$(CXX) $(LDFLAGS) -o wsfit $(OBJS) $(LDLIBS)

%.o: %.cxx
	$(CXX) $(CPPFLAGS) -c $<

clean:
	$(RM) $(OBJS)

distclean: clean
	$(RM) wsfit

wsfit.cxx

#include "Math/Integrator.h"
#include "Math/Functor.h"

#include "TF1.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TFile.h"

ROOT::Math::Integrator ig(ROOT::Math::IntegrationOneDim::kADAPTIVE,0,0.001,1000);

double _R = 6.38, _a = 0.535, _w = 0.0, _A = 197;
double FFatZero = -1;

double Z = 79.;
double pi = 3.1415926;
double hbarc = 0.197327053;
double _q2 = 0;
double _r = 0;
double _t = 0;


double rho( double r, double R=6.38, double a=0.535, double w = 0.0 ) {
    // R is the radius
    // a is the skin depth
    // w is the shape parameter

    double b = exp( (r - R) / a );
    double c = 1 + (fabs(w)*w * r*r / R*R);
    double rho0 = (3 * _A) / ( 4 * pi * pow(R, 3) );
    // rho0=1.0;
    return (rho0 * c) / ( 1 + b );
    // this 0.936 is to normalize the FF to be = 1 at q^2 = 0
}

// perform 1D integral over r at a given q
double FormFactorIntegral(double x ){
    double r = x;
    double q = sqrt(_q2);

    double rv = r * rho(r, _R, _a, _w) * sin( q * r );
    return rv;
}

double FormFactor( double *x, double *par ){
    double q2 = x[0];
    _q2 = q2 / ( hbarc * hbarc );
    double q = sqrt(_q2);


    double a0 = 4 * pi / (Z * q);
    
    double v = ig.Integral( 0, 100 );

    // self normalizing
    if ( FFatZero < 0 ){
        FFatZero = 1.0;
        double xatzero[] = {1e-9};
        FFatZero = FormFactor( xatzero, par );
    }

    return a0 * v / FFatZero;
}

double f1Coherent( double *x, double *p ){
    double xp = x[0];
    _R = p[0];
    _a = p[1];
    _w = 0;
    double A = p[2];
    double parp = 0;

    double rrr= A * pow(FormFactor( &xp, &parp ), 2);
    // cout << "[" << x[0] << "] = " << rrr << endl;
    return rrr;
}

double f1Incoherent( double *x, double *p ){
    double t = x[0];
    double A = p[0];
    double Q02 = p[1];
    return (A/Q02)/pow(1.+t/Q02,2.);
}

double f1Total( double *x, double *p ){
    double t = x[0];
    _R = p[0];
    _a = p[1];
    double B = p[2];
    double A = p[3];
    double Q02 = p[4];
    double FF2 = pow(FormFactor( &t, &p[0] ), 2);
    // 
    return B * FF2 + (A/Q02)/pow(1.+t/Q02,2.) ;
}

int main(){

    ROOT::Math::Functor1D FFInt(&FormFactorIntegral);
    ig.SetRelTolerance(0.001);
    ig.SetAbsTolerance(0.0001);
    ig.SetFunction(FFInt);


    TF1 * tf1Coherent = new TF1("tf1Coherent", &f1Coherent, 1e-9, 1, 3);
    tf1Coherent->SetParameters( 6.38, 0.53, 1 );
    tf1Coherent->SetNpx(2000);
    tf1Coherent->SetLineColor( kRed );
    tf1Coherent->Draw("");

    TF1 * tf1Incoherent = new TF1("tf1Incoherent", &f1Incoherent, 1e-9, 1, 2);
    tf1Incoherent->SetParameters( 1, 1 );
    tf1Incoherent->SetNpx(2000);
    tf1Incoherent->SetLineColor( kBlue );


    TF1 * tf1Total = new TF1("tf1Total", &f1Total, 1e-9, 1, 5);
    tf1Total->SetParameters( 4.38, 0.53, 1, 1, 1 );
    tf1Total->SetNpx(2000);
    tf1Total->SetLineColor( kBlack );
    tf1Total->Draw("");
    tf1Total->SetParNames("R", "a", "w", "A", "Q0^2");
    tf1Total->SetParLimits( 0, 0, 10 );
    // tf1Total->SetParLimits( 1, 0, 1 );
    // tf1Total->SetParLimits( 2, 0, 1 );
    tf1Total->SetParLimits( 3, 0, 1e6 );
    tf1Total->SetParLimits( 4, 0, 1e6 );

    TFile *fIn = new TFile( "fpol_coho_wt.root" );
    TH1 * ht = (TH1*)fIn->Get("ht");
    ht->GetXaxis()->SetRangeUser(0, 0.05);
    ht->SetTitle( "Coherent Phi candidates; t (GeV/c)^{2};counts" );
    ht->Draw("pe");
    ht->Draw("same hist");;

    gPad->SetLogy();

    ht->Fit( tf1Total, "R", "", 0, 0.01 );
    ht->Fit( tf1Total, "R", "", 0, 0.05 );

    tf1Coherent->SetParameters( tf1Total->GetParameter(0), tf1Total->GetParameter(1), tf1Total->GetParameter(2) );
    tf1Incoherent->SetParameters( tf1Total->GetParameter(3), tf1Total->GetParameter(4) );

    tf1Coherent->Draw("SAME");
    tf1Incoherent->Draw("SAME");
    gPad->Print("wsfit.pdf");


    TFile *f = new TFile("wsfit.root", "RECREATE");
    tf1Coherent->Write();
    tf1Incoherent->Write();
    tf1Total->Write();
    f->Close();
    
    return 0;
}