#include <math.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <time.h>

double weight[22];
double sig = 0.065;

TRandom2 *r = new TRandom2();
r->SetSeed(1045789);

double L_Ev[20] = { 23.315,
    28.865,
    31.135,
    33.153,
    34.667,
    36.432,
    38.198,
    39.964,
    42.234,
    44.505,
    47.532,
    50.559,
    54.090,
    57.874,
    61.405,
    65.694,
    70.234,
    75.784,
    82.847,
    96.721 };

double PKam[20] = {	0.595,
    0.641,
    0.803,
    0.882,
    0.752,
    0.595,
    0.618,
    0.650,
    0.632,
    0.493,
    0.424,
    0.289,
    0.479,
    0.414,
    0.470,
    0.609,
    0.674,
    0.734,
    0.595,
    0.539 };

double L_Length[22] = {
    160.,179.,191.,88.,138.,214.,146.,349.,351.,295.,401.,431.,561.,754.,830.,783.,138.,141.,735.,709.,986.,712.};

double L_Flux[22] = {
    7.62,3.40,2.22,1.95,1.88,1.84,1.83,0.93,0.85,0.30,0.19,0.18,0.15,0.094,0.061,0.043,0.21,0.32,0.14,0.13,0.12,0.19};

void fit_test_f() {
    
    double sum_flux=0;
    for(int i=0; i<22; i++) {
        sum_flux += L_Flux[i];
    }
    
    for(int i=0;i<22;i++) {
        weight[i] = (L_Flux[i])/sum_flux;
    }
    
    double OscSmear(double *L_E, double *par) {
        double EiSmearSum = 0;
        for(int i=0; i<1100; i++) {
            
            double Esig = sig*(sqrt(180./(L_E[0])));
            double Eo = 180./L_E[0];
            
            double Ej = r->Gaus(Eo,Esig); ///*** THIS IS WHERE THE RANDOM NUMBER IS GENERATED ***
            double Lj = L_Length[i/50];
            
            double P_current = 1 - (par[0] * sin(1.27*par[1]*Lj/Ej)*sin(1.27*par[1]*Lj/Ej)); // *** NB: If replace Ej in this line with Eo of fixed value, the fit will converge
            EiSmearSum = EiSmearSum + P_current*weight[i/50];
            
        }
        return EiSmearSum/(50);
    }
    
    TGraph *Original2 = new TGraph(20,L_Ev,PKam);
    TF1 *osc = new TF1("Oscillation", OscSmear,20,100,2);
    
    osc->SetParameters(.88,.070);
    //osc->SetParLimits(0, 0, 1.);
    //osc->SetParLimits(1, 0.060, 0.099);
    
    TCanvas *c1 = new TCanvas("c1", "test", 600,600);
    
    Original2->SetMarkerStyle(33);
    Original2->Fit(osc,"R","EX0",20,100);
    Original2->Draw("AP");
    
    
}