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"
#include "TStyle.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
}

double rhoZ( double r, double z, double R=6.38, double a=0.535, double w = 0.0 ){
    return rho( sqrt( r*r + z*z ), R, a, w);
}
double rhoZIntegral( double z ){
    return rho( sqrt( _r*_r + z*z ), _R, _a, _w );
}
double f1RhoZIntegral( double *x, double *p ) {
    _r = x[0];
    _R = p[0];
    _a = p[1];
    _w = p[2];

    return igZ.Integral( 0, 100 );
    // return rhoZIntegral( 0 );

}

double FormFactorRhoZIntegral(double y ){
    double q = sqrt(_q2);
    double p[] = {_R, _a, _w};
    double x[] = {y};

    double rv =  y * f1RhoZIntegral( x, p ) * sin( q * y );
    return rv;
}

double FormFactorRhoZ( double *x, double *p ){
    // _R = p[0];
    // _a = p[1];
    // _w = p[2];
    double q2 = x[0];
    _q2 = q2 / ( hbarc * hbarc );
    double q = sqrt(_q2);


    double a0 = 4 * pi / (Z * q);
    
    double v = igFZ.Integral( 0, 100 );
    
    // self normalizing
    if ( FFZatZeroRhoZ < 0 ){
        FFZatZeroRhoZ = 1.0;
        double xatzero[] = {1e-9};
        FFZatZeroRhoZ = FormFactorRhoZ( xatzero, p );
    }

    // cout << "FormFactorRhoZ( " << q << ") = " << a0 * v / FFatZero << endl;
    // cout << "v=" << a0 * v / FFZatZeroRhoZ << endl;
    return a0 * v / FFZatZeroRhoZ;
}

// 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];
    double A = p[0];
    _R = p[1];
    _a = p[2];
    _w = p[3];
    
    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];
    double p2 = p[2];
    
    return A * exp( t/Q02 ) * pow(t, 1.5) * (1.0 / (1.0 + p2*t));
    // return (A/Q02)/pow(1.0+t/Q02,2.) * (log(1.0 + t/Q02)) * pow( t*t, p2 ); ;
}

double f1Total( double *x, double *p ){
    double t = x[0];
    
    // coherent parameters
    double ACoh = p[0];
    _R = p[1];
    _a = p[2];
    _w = p[3];
    
    // incoherent parameters
    double AInco = p[4];
    double Q02 = p[5];
    double p2 = p[6];

    double cp[] = { ACoh, _R, _a, _w };
    double pp[] = { AInco, Q02, p2 };
    return f1Coherent( x, cp ) + f1Incoherent( x, pp );
}

int main(){

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

    gStyle->SetOptFit(1111);

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

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


    TF1 * tf1Total = new TF1("tf1Total", &f1Total, 1e-9, 1, 7);
    
    tf1Total->SetNpx(2000);
    tf1Total->SetLineColor( kRed );
    tf1Total->Draw("");
    tf1Total->SetParNames("ACoh", "R", "a", "w", "AInco", "Q0^2", "p");
    tf1Total->SetParameters( 1e4, 4.38, 0.53, 1, 1, 1, 1, 1 );
    tf1Total->SetParLimits( 0, 0, 1e9 );
    // tf1Total->SetParLimits( 1, 0, 1 );
    // tf1Total->SetParLimits( 2, 0, 1 );
    // tf1Total->SetParLimits( 3, 0, 1e9 );
    // tf1Total->SetParLimits( 4, 0, 1e9 );
    // tf1Total->SetParLimits( 5, 0, 1e4 );
    // tf1Total->FixParameter( 1, 4.38 );
    tf1Total->FixParameter( 2, 0.53 );
    // tf1Total->FixParameter( 3, 0.0 );

    tf1Total->FixParameter( 4, 0 );

    TFile *fIn = new TFile( "1D_cophipt2_hist.root" );
    TH1 * ht = (TH1*)fIn->Get("ht");
    ht->GetXaxis()->SetRangeUser(0, 0.2);
    ht->GetYaxis()->SetRangeUser(1e-1, 1e5);
    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 );
    tf1Total->FixParameter( 1, tf1Total->GetParameter(1) );
    tf1Total->FixParameter( 2, tf1Total->GetParameter(2) );
    tf1Total->FixParameter( 3, tf1Total->GetParameter(3) );
    // tf1Total->FixParameter( 4, 0 );
    // tf1Total->FixParameter( 2, tf1Total->GetParameter(1) );
    tf1Total->ReleaseParameter(4);
    ht->Fit( tf1Total, "R", "", 0, 0.5 );
    gPad->Print("wsfit2.pdf");

    // ht->Fit( tf1Incoherent, "R", "", .02, 0.9 );

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

    gStyle->SetOptFit(11111);

    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;
}