#include "TCanvas.h"
#include "TPad.h"
#include "TStyle.h"
#include "TH1.h"
#include "TF1.h"
#include "TString.h"
#include "TRandom3.h"
#include "TGraph.h"
#include "TMultiGraph.h"
#include "TMath.h"

TRandom3 r;
//const rdm = r.Rndm();
Int_t nop=250;				// number of points on the curve
Double_t xpt, ypt, del = 0.9975/nop;	// del is the space between two consecutive points from x-min to x-max (the domain of the function)

//auto avf = new TH1F("avf", "Average of Functions", 100, 0.0025, 1.0);
//avf->Draw();
//avf->SetMaximum(1.e0);

Double_t dcsep(Double_t x,Double_t *p)
{
 Double_t fcn = (((2.0*p[0]*p[0]*p[1]*p[1]*p[2]*p[2]*p[3]*p[3])/(x*x*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0)+2.0*p[3]*p[3]*x*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0)+(x*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0)*sqrt(x*x*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0)+4.0*p[3]*p[3]*x*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0)))))*((2.0*p[3]*p[3])/(2.0*p[3]*p[3]+x*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0)+sqrt(x*x*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0)+4.0*p[3]*p[3]*x*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0))))*((4.0*p[3]*p[3])/(4.0*p[3]*p[3]+x))*(((0.71)/(0.71+x))*((0.71)/(0.71+x))*((0.71)/(0.71+x))*((0.71)/(0.71+x))*(8.0*p[3]*p[3]*p[3]*p[3]*cos(p[4]*p[5]/360.0)*cos(p[4]*p[5]/360.0)+(2.0*p[3]*p[3]*x*cos(p[4]*p[5]/360.0)*cos(p[4]*p[5]/360.0)+4.0*p[3]*p[3]*x*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0)+x*x*sin(p[4]*p[5]/360.0)*sin(p[4]*p[5]/360.0))*((p[6]*p[6])/(p[7]*p[7])))/(8.0*p[3]*p[3]*p[3]*p[3])));

 return fcn;
}

void DataPointsCode()
{
 const Int_t nof = 250;			// number of functions to be generated

 TF1 *f001[nof];
 TString fna;

 //=== generating the average of functions for ep scattering ===//

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;

  TF1 *f001 = new TF1("fna", dcsep, 0.0, 0.0025, 1.0);
  f001[i]->SetParameter(0,1.0);
  f001[i]->SetParameter(1,1.0);
  f001[i]->SetParameter(2,1.0/137.0);
  f001[i]->SetParameter(3,gRandom->Gaus(0.938272081,0.000000006));
  f001[i]->SetParameter(4,60.6018);
  f001[i]->SetParameter(5,TMath::Pi());
  f001[i]->SetParameter(6,2.792847351);
  f001[i]->SetParameter(7,gRandom->Gaus(0.963,0.010));
 }
 //f001->Draw();

 //=== getting the average of the above function===//

 TGraph *av001 = new TGraph();

 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av001->SetPoint(j, xpt, 0.0);
 }

 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av001->GetPoint(j, xpt, ypt);
   av001->SetPoint(j, xpt, ypt + f001[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 av001->Print();
}