//рис.1b верхний средний 6<pT(trig)<8, 3<pT(assoc)<4

#include "TAxis.h"
#include <iostream>
#include <fstream>

#include "TF1.h"
#include "TMath.h"
#include "TMinuit.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TFrame.h"

#include "TGraphErrors.h"


using namespace std;


Double_t radius;

Double_t Sx(Double_t *x, Double_t *par)// my integrand function
{

    Double_t e= cos(par[3]);
    Double_t z= e * par[0];
    Double_t S = ( pow(  par[2] - par[0]* par[1] , 2. ) * (par[3]*par[3]) + 2.
                * (par[1] - par[0] * par[2] ) * ( par[2] - par[0] * par[1])
                * e * par[3] + pow( par[1] - par[0]* par[2],2.) )
                / ( (1. - pow(par[0],2))*(1. + pow(par[3],2) - 2.* z * par[3]) );
    Double_t G = 2 * ( 1 - pow(par[0],2))
                * TMath::Exp( - (  pow(par[1],2) - 2. * par[0] * par[1] * par[3] + pow(par[2],2))
                / (1 - pow(par[0],2))  ) /TMath::Pi();
    Double_t podint = G * x[0] * ( 1. + S ) * TMath::Exp(S) / ( 1 + pow(x[0],2) - 2 * z * x[0]);
 return podint;

}


struct MyIntegFunc
{

    MyIntegFunc(TF1 *f):
        fFunc(f) {}
        double operator() (double *x, double *p)
        {

            double radius = 100;

        /* std::copy(p,p+3,param); -- commented out because otherwise  I get  occurrence of the Error:
          Function copy(p,p+3,param) is not defined in current scope  C:\root\obshaya1.C(44) */

        // the similar way:
            param[0]= p[0];
            param[1]= p[1];
            param[2]= p[2];
            param[3]= *x;


            return fFunc->Integral(0.,radius,param);

        }
        TF1 *fFunc;
        double param[4];

};



void obshaya1()
{

    // the data set
    const Int_t n = 35;
    Double_t ey[n];
    //Double_t ex[n];
    Double_t radius;
    Double_t y[n]  = {  0.169565217391304,	0.166086956521739,	0.157391304347826,	0.166086956521739,	0.169565217391304,
                        0.157391304347826,	0.169565217391304,	0.18,	            0.204347826086957,	0.26,
                        0.200869565217391,	0.171304347826087,	0.174782608695652,	0.159130434782609,	0.164347826086957,
                        0.159130434782609,	0.171304347826087,	0.166086956521739,	0.17304347826087,	0.160869565217391,
                        0.171304347826087,	0.166086956521739,	0.167826086956522,	0.166086956521739,	0.174782608695652,
                        0.181739130434783,	0.185217391304348,	0.193913043478261,	0.18,	            0.176521739130435,
                        0.169565217391304,	0.174782608695652,	0.169565217391304,	0.171304347826087,	0.167826086956522};



    Double_t x[n] ={ -1.51325301204819,	-1.34301204819277,	-1.17277108433735,	-1.02144578313253,	-0.870120481927711,
                        -0.718795180722892,	-0.567469879518072,	-0.397228915662651,	-0.226987951807229,	    0,
                        0.226987951807229,	0.378313253012048,	0.529638554216867,	0.699879518072289,	0.870120481927711,
                        1.02144578313253,	1.17277108433735,	1.32409638554217,	1.47542168674699,	1.62674698795181,
                        1.79698795180723,	1.96722891566265,	2.11855421686747,	2.3455421686747,	2.66710843373494,
                        2.91301204819277,	3.06433734939759,	3.21566265060241,	3.36698795180723,	3.51831325301205,
                        3.76421686746988,	4.01012048192771,	4.23710843373494,	4.46409638554217,	4.59650602409639};


    Double_t ex[n] = { 0.0756626506024101,	0.0945783132530099,	0.0756626506024101,	0.0756626506024095,	0.07566265060241,
                            0.0756626506024095,	0.0756626506024095,	0.094578313253012,	0.07566265060241,	0.151325301204819,
                            0.07566265060241,	0.0756626506024095,	0.0756626506024095,	0.094578313253012,	0.07566265060241,
                            0.0756626506024095,	0.0756626506024101,	0.07566265060241,	0.0756626506024101,	0.07566265060241,
                            0.09457831325301,	0.0756626506024101,	0.0756626506024098,	0.15132530120482,	0.17024096385542,
                            0.0756626506024101,	0.0756626506024101,	0.0756626506024098,	0.0756626506024101,	0.0756626506024101,
                            0.17024096385542,	0.0756626506024103,	0.15132530120482,	0.0756626506024101,	0.0567469879518052};

    for(int i=0; i < n;  i++)
    {

        ey[i]=0.05*y[i];
    }


    TCanvas *c1 = new TCanvas();
	c1->SetBorderMode(0);
	c1->SetFillColor(10);
	c1->SetGrid();
	c1->GetFrame()->SetBorderSize(10);




    TGraphErrors *histo = new TGraphErrors(n,x,y,ex,ey);
    histo->SetMarkerStyle(3);
    histo->SetMarkerSize(1.0);

    histo->GetYaxis()->SetRangeUser(0.,0.33);
    histo->GetXaxis()->SetRangeUser(-1.513253012,4.653253012);
    histo->GetYaxis()->SetTitle("#frac{1}{N_{trigger}} #frac{dN}{d(#Delta#phi)}");
    histo->GetXaxis()->SetTitle("#Delta#phi");
    histo->GetYaxis()->SetTitleOffset(1.4);
    histo->GetXaxis()->SetTitleOffset(0.55);
    histo->GetYaxis()->SetNoExponent(kTRUE);
    histo->GetYaxis()->SetTitleSize(0.03);
    histo->GetYaxis()->SetLabelSize(0.03);
    histo->GetXaxis()->SetLabelSize(0.03);

    histo->Draw("AP");

    gStyle->SetOptFit(kTRUE);

    TF1 *funcx = new TF1("funcx", Sx, -1.513253012,4.653253012,4);

    MyIntegFunc *integ = new MyIntegFunc(funcx);
    TF1 *ifuncx = new TF1("ifuncx", integ, -1.513253012,4.653253012, 3, "MyIntegFunc");

    ifuncx->SetParameter(2, 0.75);
    ifuncx->SetParameter(1, 0.118);
    ifuncx->SetParameter(0, 0.398);

    ifuncx->SetParLimits(2, -0.999, 0.999);
    ifuncx->SetParLimits(1, -0.999, 0.999);
    ifuncx->SetParLimits(0, 0.001, 0.999);

    ifuncx->SetLineColor(kBlue);

    histo->Fit("ifuncx","EVMRL","elpv");
    histo->Draw("same");


}