void DataPointStatNewEPEQDCSProbColorWithMRF02()
{
 const Int_t nof = 250;					// number of functions to be generated
 TCanvas *can = new TCanvas("c1", "c1", 1000,1000);	// setting the canvas
 //gStyle->SetOptStat("nerf");
 //gStyle->SetOptStat();
 gStyle->SetOptStat(0);
 can->Draw();

 TH1F *avf = new TH1F("avf", "Average of Functions", 100, 0.0025, 1.0);
 avf->Draw();

 TF1 *fpf[nof], *qbs[nof], *qbc[nof], *qes[nof], *qec[nof], *fpa[nof], *dbs[nof], *abc[nof], *aes[nof], *aec[nof];
 TF1 *fpb[nof], *bbs[nof], *bbc[nof], *bes[nof], *bec[nof], *fpc[nof], *cbs[nof], *cbc[nof], *ces[nof], *cec[nof];
 TF1 *fpe[nof], *ebs[nof], *ebc[nof], *ees[nof], *eec[nof], *fpg[nof], *gbs[nof], *gbc[nof], *ges[nof], *gec[nof];

 TString fna;
 TRandom3 r;
 Double_t rdm = r.Rndm();

//=== generating the functions ===//
//=======================Theta=27.00
 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  //choice: transfer momentum square; electron(finite)proton
  fpf[i] = new TF1("fna","(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[3]*[3])/(2.0*[3]*[3]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((4.0*[3]*[3])/(4.0*[3]*[3]+x))*(((0.71)/(0.71+x))*((0.71)/(0.71+x))*((0.71)/(0.71+x))*((0.71)/(0.71+x))*([7]*[7]*[8]*[8]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+[6]*[6]*[7]*[7]*(2.0*[3]*[3]*x*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(8.0*[3]*[3]*[3]*[3]))))", 0.0025, 1.0);
  fpf[i]->SetParameter(0,1.0);
  fpf[i]->SetParameter(1,1.0);
  fpf[i]->SetParameter(2,1.0/137.0);
  fpf[i]->SetParameter(3,gRandom->Gaus(0.938272081,0.000000006));
  fpf[i]->SetParameter(4,27.0);
  fpf[i]->SetParameter(5,TMath::Pi());
  fpf[i]->SetParameter(6,2.792847351);
  fpf[i]->SetParameter(7,gRandom->Gaus(0.979,0.022)+gRandom->Gaus(0.0,0.010));
  fpf[i]->SetParameter(8,gRandom->Gaus(0.914,0.024)+gRandom->Gaus(0.0,0.010));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 2/3*(e-u) and 1/3*(e-d) bare mass with modified recoil factor
  qbs[i] = new TF1("fna","(2.0/3.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(1.0/3.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  qbs[i]->SetParameter(0,1.0);
  qbs[i]->SetParameter(1,2.0/3.0);
  qbs[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   qbs[i]->SetParameter(3,gRandom->Gaus(0.0017,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        qbs[i]->SetParameter(3,gRandom->Gaus(0.0033,0.0008));
                                       }
                                       else {
                                             qbs[i]->SetParameter(3,(gRandom->Gaus(0.0017,0.0008)+gRandom->Gaus(0.0033,0.0008))/2.0);
					    }
  qbs[i]->SetParameter(4,27.0);
  qbs[i]->SetParameter(5,TMath::Pi());
  qbs[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   qbs[i]->SetParameter(7,gRandom->Gaus(0.0041,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        qbs[i]->SetParameter(7,gRandom->Gaus(0.0058,0.0008));
                                       }
                                       else {
                                             qbs[i]->SetParameter(7,(gRandom->Gaus(0.0041,0.0008)+gRandom->Gaus(0.0058,0.0008))/2.0);
                                            }
  qbs[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000023));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 72/81*(e-u) and 9/81*(e-d) bare mass with modified recoil factor
  qbc[i] = new TF1("fna","(72.0/81.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(9.0/81.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  qbc[i]->SetParameter(0,1.0);
  qbc[i]->SetParameter(1,2.0/3.0);
  qbc[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   qbc[i]->SetParameter(3,gRandom->Gaus(0.0017,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        qbc[i]->SetParameter(3,gRandom->Gaus(0.0033,0.0008));
                                       }
                                       else {
                                             qbc[i]->SetParameter(3,(gRandom->Gaus(0.0017,0.0008)+gRandom->Gaus(0.0033,0.0008))/2.0);
                                            }
  qbc[i]->SetParameter(4,27.0);
  qbc[i]->SetParameter(5,TMath::Pi());
  qbc[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   qbc[i]->SetParameter(7,gRandom->Gaus(0.0041,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        qbc[i]->SetParameter(7,gRandom->Gaus(0.0058,0.0008));
                                       }
                                       else {
                                             qbc[i]->SetParameter(7,(gRandom->Gaus(0.0041,0.0008)+gRandom->Gaus(0.0058,0.0008))/2.0);
                                            }
  qbc[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000023));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 2/3*(e-u) and 1/3*(e-d) effective mass with modified recoil factor
  qes[i] = new TF1("fna","(2.0/3.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(1.0/3.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  qes[i]->SetParameter(0,1.0);
  qes[i]->SetParameter(1,2.0/3.0);
  qes[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   qes[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                  }
                  else if (rdm > 0.75) {
                                        qes[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                       }
                                       else {
                                             qes[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                            }
  qes[i]->SetParameter(4,27.0);
  qes[i]->SetParameter(5,TMath::Pi());
  qes[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   qes[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                  }
                  else if (rdm > 0.75) {
                                        qes[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                       }
					else {
                                             qes[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                            }
  qes[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000006));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 72/81*(e-u) and 9/81*(e-d) effective mass with modified recoil factor
  qec[i] = new TF1("fna","(72.0/81.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(9.0/81.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  qec[i]->SetParameter(0,1.0);
  qec[i]->SetParameter(1,2.0/3.0);
  qec[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   qec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                  }
                  else if (rdm > 0.75) {
                                        qec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                       }
                                       else {
                                             qec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                            }
  qec[i]->SetParameter(4,27.0);
  qec[i]->SetParameter(5,TMath::Pi());
  qec[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   qec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                  }
                  else if (rdm > 0.75) {
                                        qec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                       }
                                       else {
                                             qec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                            }
  qec[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000006));
 }

//=== getting the average ===//

 Int_t nop=400;
 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 av1 = new TGraph(nop);
 TGraph *av1 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av1->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av1->GetPoint(j, xpt, ypt);
   av1->SetPoint(j, xpt, ypt + fpf[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av2 = new TGraph(nop);
 TGraph *av2 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av2->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av2->GetPoint(j, xpt, ypt);
   av2->SetPoint(j, xpt, ypt + qbs[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av3 = new TGraph(nop);
 TGraph *av3 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av3->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av3->GetPoint(j, xpt, ypt);
   av3->SetPoint(j, xpt, ypt + qbc[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av4 = new TGraph(nop);
 TGraph *av4 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av4->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av4->GetPoint(j, xpt, ypt);
   av4->SetPoint(j, xpt, ypt + qes[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av5 = new TGraph(nop);
 TGraph *av5 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av5->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av5->GetPoint(j, xpt, ypt);
   av5->SetPoint(j, xpt, ypt + qec[i]->Eval(xpt)/(Double_t)nof);
  }
 }

//=== generating the functions ===//
//========================Theta=27.0
 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  //choice: transfer momentum square; electron(finite)proton
  fpa[i] = new TF1("fna","(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[3]*[3])/(2.0*[3]*[3]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((4.0*[3]*[3])/(4.0*[3]*[3]+x))*(((0.71)/(0.71+x))*((0.71)/(0.71+x))*((0.71)/(0.71+x))*((0.71)/(0.71+x))*([7]*[7]*[8]*[8]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+[6]*[6]*[7]*[7]*(2.0*[3]*[3]*x*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(8.0*[3]*[3]*[3]*[3]))))", 0.0025, 1.0);
  fpa[i]->SetParameter(0,1.0);
  fpa[i]->SetParameter(1,1.0);
  fpa[i]->SetParameter(2,1.0/137.0);
  fpa[i]->SetParameter(3,gRandom->Gaus(0.938272081,0.000000006));
  fpa[i]->SetParameter(4,27.0);
  fpa[i]->SetParameter(5,TMath::Pi());
  fpa[i]->SetParameter(6,2.792847351);
  fpa[i]->SetParameter(7,gRandom->Gaus(0.974,0.012)+gRandom->Gaus(0.0,0.006));
  fpa[i]->SetParameter(8,gRandom->Gaus(0.959,0.017)+gRandom->Gaus(0.0,0.006));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 2/3*(e-u) and 1/3*(e-d) bare mass with modified recoil factor
  dbs[i] = new TF1("fna","(2.0/3.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(1.0/3.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  dbs[i]->SetParameter(0,1.0);
  dbs[i]->SetParameter(1,2.0/3.0);
  dbs[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   dbs[i]->SetParameter(3,gRandom->Gaus(0.0017,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        dbs[i]->SetParameter(3,gRandom->Gaus(0.0033,0.0008));
                                       }
                                       else {
                                             dbs[i]->SetParameter(3,(gRandom->Gaus(0.0017,0.0008)+gRandom->Gaus(0.0033,0.0008))/2.0);
					    }
  dbs[i]->SetParameter(4,27.0);
  dbs[i]->SetParameter(5,TMath::Pi());
  dbs[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   dbs[i]->SetParameter(7,gRandom->Gaus(0.0041,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        dbs[i]->SetParameter(7,gRandom->Gaus(0.0058,0.0008));
                                       }
                                       else {
                                             dbs[i]->SetParameter(7,(gRandom->Gaus(0.0041,0.0008)+gRandom->Gaus(0.0058,0.0008))/2.0);
                                            }
  dbs[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000023));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 72/81*(e-u) and 9/81*(e-d) bare mass with modified recoil factor
  abc[i] = new TF1("fna","(72.0/81.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(9.0/81.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  abc[i]->SetParameter(0,1.0);
  abc[i]->SetParameter(1,2.0/3.0);
  abc[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   abc[i]->SetParameter(3,gRandom->Gaus(0.0017,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        abc[i]->SetParameter(3,gRandom->Gaus(0.0033,0.0008));
                                       }
                                       else {
                                             abc[i]->SetParameter(3,(gRandom->Gaus(0.0017,0.0008)+gRandom->Gaus(0.0033,0.0008))/2.0);
                                            }
  abc[i]->SetParameter(4,27.0);
  abc[i]->SetParameter(5,TMath::Pi());
  abc[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   abc[i]->SetParameter(7,gRandom->Gaus(0.0041,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        abc[i]->SetParameter(7,gRandom->Gaus(0.0058,0.0008));
                                       }
                                       else {
                                             abc[i]->SetParameter(7,(gRandom->Gaus(0.0041,0.0008)+gRandom->Gaus(0.0058,0.0008))/2.0);
                                            }
  abc[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000023));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 2/3*(e-u) and 1/3*(e-d) effective mass with modified recoil factor
  aes[i] = new TF1("fna","(2.0/3.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(1.0/3.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  aes[i]->SetParameter(0,1.0);
  aes[i]->SetParameter(1,2.0/3.0);
  aes[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   aes[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                  }
                  else if (rdm > 0.75) {
                                        aes[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                       }
                                       else {
                                             aes[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                            }
  aes[i]->SetParameter(4,27.0);
  aes[i]->SetParameter(5,TMath::Pi());
  aes[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   aes[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                  }
                  else if (rdm > 0.75) {
                                        aes[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                       }
					else {
                                             aes[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                            }
  aes[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000006));
 }
 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 72/81*(e-u) and 9/81*(e-d) effective mass with modified recoil factor
  aec[i] = new TF1("fna","(72.0/81.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(9.0/81.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  aec[i]->SetParameter(0,1.0);
  aec[i]->SetParameter(1,2.0/3.0);
  aec[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   aec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                  }
                  else if (rdm > 0.75) {
                                        aec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                       }
                                       else {
                                             aec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                            }
  aec[i]->SetParameter(4,27.0);
  aec[i]->SetParameter(5,TMath::Pi());
  aec[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   aec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                  }
                  else if (rdm > 0.75) {
                                        aec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                       }
                                       else {
                                             aec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                            }
  aec[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000006));
 }

//=== getting the average ===//

 //auto av6 = new TGraph(nop);
 TGraph *av6 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av6->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av6->GetPoint(j, xpt, ypt);
   av6->SetPoint(j, xpt, ypt + fpa[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av7 = new TGraph(nop);
 TGraph *av7 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av7->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av7->GetPoint(j, xpt, ypt);
   av7->SetPoint(j, xpt, ypt + dbs[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av8 = new TGraph(nop);
 TGraph *av8 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av8->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av8->GetPoint(j, xpt, ypt);
   av8->SetPoint(j, xpt, ypt + abc[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av9 = new TGraph(nop);
 TGraph *av9 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av9->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av9->GetPoint(j, xpt, ypt);
   av9->SetPoint(j, xpt, ypt + aes[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av10 = new TGraph(nop);
 TGraph *av10 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av10->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av10->GetPoint(j, xpt, ypt);
   av10->SetPoint(j, xpt, ypt + aec[i]->Eval(xpt)/(Double_t)nof);
  }
 }

//=== generating the functions ===//
//========================Theta=35.0
 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  //choice: transfer momentum square; electron(finite)proton
  fpb[i] = new TF1("fna","(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[3]*[3])/(2.0*[3]*[3]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((4.0*[3]*[3])/(4.0*[3]*[3]+x))*(((0.71)/(0.71+x))*((0.71)/(0.71+x))*((0.71)/(0.71+x))*((0.71)/(0.71+x))*([7]*[7]*[8]*[8]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+[6]*[6]*[7]*[7]*(2.0*[3]*[3]*x*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(8.0*[3]*[3]*[3]*[3]))))", 0.0025, 1.0);
  fpb[i]->SetParameter(0,1.0);
  fpb[i]->SetParameter(1,1.0);
  fpb[i]->SetParameter(2,1.0/137.0);
  fpb[i]->SetParameter(3,gRandom->Gaus(0.938272081,0.000000006));
  fpb[i]->SetParameter(4,35.0);
  fpb[i]->SetParameter(5,TMath::Pi());
  fpb[i]->SetParameter(6,2.792847351);
  fpb[i]->SetParameter(7,gRandom->Gaus(0.994,0.008)+gRandom->Gaus(0.0,0.006));
  fpb[i]->SetParameter(8,gRandom->Gaus(0.959,0.016)+gRandom->Gaus(0.0,0.006));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 2/3*(e-u) and 1/3*(e-d) bare mass with modified recoil factor
  bbs[i] = new TF1("fna","(2.0/3.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(1.0/3.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  bbs[i]->SetParameter(0,1.0);
  bbs[i]->SetParameter(1,2.0/3.0);
  bbs[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   bbs[i]->SetParameter(3,gRandom->Gaus(0.0017,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        bbs[i]->SetParameter(3,gRandom->Gaus(0.0033,0.0008));
                                       }
                                       else {
                                             bbs[i]->SetParameter(3,(gRandom->Gaus(0.0017,0.0008)+gRandom->Gaus(0.0033,0.0008))/2.0);
					    }
  bbs[i]->SetParameter(4,35.0);
  bbs[i]->SetParameter(5,TMath::Pi());
  bbs[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   bbs[i]->SetParameter(7,gRandom->Gaus(0.0041,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        bbs[i]->SetParameter(7,gRandom->Gaus(0.0058,0.0008));
                                       }
                                       else {
                                             bbs[i]->SetParameter(7,(gRandom->Gaus(0.0041,0.0008)+gRandom->Gaus(0.0058,0.0008))/2.0);
                                            }
  bbs[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000023));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 72/81*(e-u) and 9/81*(e-d) bare mass with modified recoil factor
  bbc[i] = new TF1("fna","(72.0/81.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(9.0/81.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  bbc[i]->SetParameter(0,1.0);
  bbc[i]->SetParameter(1,2.0/3.0);
  bbc[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   bbc[i]->SetParameter(3,gRandom->Gaus(0.0017,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        bbc[i]->SetParameter(3,gRandom->Gaus(0.0033,0.0008));
                                       }
                                       else {
                                             bbc[i]->SetParameter(3,(gRandom->Gaus(0.0017,0.0008)+gRandom->Gaus(0.0033,0.0008))/2.0);
                                            }
  bbc[i]->SetParameter(4,35.0);
  bbc[i]->SetParameter(5,TMath::Pi());
  bbc[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   bbc[i]->SetParameter(7,gRandom->Gaus(0.0041,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        bbc[i]->SetParameter(7,gRandom->Gaus(0.0058,0.0008));
                                       }
                                       else {
                                             bbc[i]->SetParameter(7,(gRandom->Gaus(0.0041,0.0008)+gRandom->Gaus(0.0058,0.0008))/2.0);
                                            }
  bbc[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000023));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 2/3*(e-u) and 1/3*(e-d) effective mass with modified recoil factor
  bes[i] = new TF1("fna","(2.0/3.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(1.0/3.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  bes[i]->SetParameter(0,1.0);
  bes[i]->SetParameter(1,2.0/3.0);
  bes[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   bes[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                  }
                  else if (rdm > 0.75) {
                                        bes[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                       }
                                       else {
                                             bes[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                            }
  bes[i]->SetParameter(4,35.0);
  bes[i]->SetParameter(5,TMath::Pi());
  bes[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   bes[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                  }
                  else if (rdm > 0.75) {
                                        bes[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                       }
					else {
                                             bes[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                            }
  bes[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000006));
 }
 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 72/81*(e-u) and 9/81*(e-d) effective mass with modified recoil factor
  bec[i] = new TF1("fna","(72.0/81.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(9.0/81.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  bec[i]->SetParameter(0,1.0);
  bec[i]->SetParameter(1,2.0/3.0);
  bec[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   bec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                  }
                  else if (rdm > 0.75) {
                                        bec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                       }
                                       else {
                                             bec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                            }
  bec[i]->SetParameter(4,35.0);
  bec[i]->SetParameter(5,TMath::Pi());
  bec[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   bec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                  }
                  else if (rdm > 0.75) {
                                        bec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                       }
                                       else {
                                             bec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                            }
  bec[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000006));
 }

//=== getting the average ===//

 //auto av11 = new TGraph(nop);
 TGraph *av11 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av11->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av11->GetPoint(j, xpt, ypt);
   av11->SetPoint(j, xpt, ypt + fpb[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av12 = new TGraph(nop);
 TGraph *av12 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av12->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av12->GetPoint(j, xpt, ypt);
   av12->SetPoint(j, xpt, ypt + bbs[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av13 = new TGraph(nop);
 TGraph *av13 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av13->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av13->GetPoint(j, xpt, ypt);
   av13->SetPoint(j, xpt, ypt + bbc[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av14 = new TGraph(nop);
 TGraph *av14 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av14->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av14->GetPoint(j, xpt, ypt);
   av14->SetPoint(j, xpt, ypt + bes[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av15 = new TGraph(nop);
 TGraph *av15 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av15->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av15->GetPoint(j, xpt, ypt);
   av15->SetPoint(j, xpt, ypt + bec[i]->Eval(xpt)/(Double_t)nof);
  }
 }

//=== generating the functions ===//
//========================Theta=39.0
 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  //choice: transfer momentum square; electron(finite)proton
  fpc[i] = new TF1("fna","(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[3]*[3])/(2.0*[3]*[3]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((4.0*[3]*[3])/(4.0*[3]*[3]+x))*(((0.71)/(0.71+x))*((0.71)/(0.71+x))*((0.71)/(0.71+x))*((0.71)/(0.71+x))*([7]*[7]*[8]*[8]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+[6]*[6]*[7]*[7]*(2.0*[3]*[3]*x*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(8.0*[3]*[3]*[3]*[3]))))", 0.0025, 1.0);
  fpc[i]->SetParameter(0,1.0);
  fpc[i]->SetParameter(1,1.0);
  fpc[i]->SetParameter(2,1.0/137.0);
  fpc[i]->SetParameter(3,gRandom->Gaus(0.938272081,0.000000006));
  fpc[i]->SetParameter(4,39.0);
  fpc[i]->SetParameter(5,TMath::Pi());
  fpc[i]->SetParameter(6,2.792847351);
  fpc[i]->SetParameter(7,gRandom->Gaus(0.996,0.006)+gRandom->Gaus(0.0,0.006));
  fpc[i]->SetParameter(8,gRandom->Gaus(0.989,0.015)+gRandom->Gaus(0.0,0.006));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 2/3*(e-u) and 1/3*(e-d) bare mass with modified recoil factor
  cbs[i] = new TF1("fna","(2.0/3.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(1.0/3.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  cbs[i]->SetParameter(0,1.0);
  cbs[i]->SetParameter(1,2.0/3.0);
  cbs[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   cbs[i]->SetParameter(3,gRandom->Gaus(0.0017,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        cbs[i]->SetParameter(3,gRandom->Gaus(0.0033,0.0008));
                                       }
                                       else {
                                             cbs[i]->SetParameter(3,(gRandom->Gaus(0.0017,0.0008)+gRandom->Gaus(0.0033,0.0008))/2.0);
					    }
  cbs[i]->SetParameter(4,39.0);
  cbs[i]->SetParameter(5,TMath::Pi());
  cbs[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   cbs[i]->SetParameter(7,gRandom->Gaus(0.0041,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        cbs[i]->SetParameter(7,gRandom->Gaus(0.0058,0.0008));
                                       }
                                       else {
                                             cbs[i]->SetParameter(7,(gRandom->Gaus(0.0041,0.0008)+gRandom->Gaus(0.0058,0.0008))/2.0);
                                            }
  cbs[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000023));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 72/81*(e-u) and 9/81*(e-d) bare mass with modified recoil factor
  cbc[i] = new TF1("fna","(72.0/81.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(9.0/81.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  cbc[i]->SetParameter(0,1.0);
  cbc[i]->SetParameter(1,2.0/3.0);
  cbc[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   cbc[i]->SetParameter(3,gRandom->Gaus(0.0017,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        cbc[i]->SetParameter(3,gRandom->Gaus(0.0033,0.0008));
                                       }
                                       else {
                                             cbc[i]->SetParameter(3,(gRandom->Gaus(0.0017,0.0008)+gRandom->Gaus(0.0033,0.0008))/2.0);
                                            }
  cbc[i]->SetParameter(4,39.0);
  cbc[i]->SetParameter(5,TMath::Pi());
  cbc[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   cbc[i]->SetParameter(7,gRandom->Gaus(0.0041,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        cbc[i]->SetParameter(7,gRandom->Gaus(0.0058,0.0008));
                                       }
                                       else {
                                             cbc[i]->SetParameter(7,(gRandom->Gaus(0.0041,0.0008)+gRandom->Gaus(0.0058,0.0008))/2.0);
                                            }
  cbc[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000023));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 2/3*(e-u) and 1/3*(e-d) effective mass with modified recoil factor
  ces[i] = new TF1("fna","(2.0/3.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(1.0/3.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  ces[i]->SetParameter(0,1.0);
  ces[i]->SetParameter(1,2.0/3.0);
  ces[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   ces[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                  }
                  else if (rdm > 0.75) {
                                        ces[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                       }
                                       else {
                                             ces[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                            }
  ces[i]->SetParameter(4,39.0);
  ces[i]->SetParameter(5,TMath::Pi());
  ces[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   ces[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                  }
                  else if (rdm > 0.75) {
                                        ces[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                       }
					else {
                                             ces[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                            }
  ces[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000006));
 }
 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 72/81*(e-u) and 9/81*(e-d) effective mass with modified recoil factor
  cec[i] = new TF1("fna","(72.0/81.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(9.0/81.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  cec[i]->SetParameter(0,1.0);
  cec[i]->SetParameter(1,2.0/3.0);
  cec[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   cec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                  }
                  else if (rdm > 0.75) {
                                        cec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                       }
                                       else {
                                             cec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                            }
  cec[i]->SetParameter(4,39.0);
  cec[i]->SetParameter(5,TMath::Pi());
  cec[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   cec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                  }
                  else if (rdm > 0.75) {
                                        cec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                       }
                                       else {
                                             cec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                            }
  cec[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000006));
 }

//=== getting the average ===//

 //auto av16 = new TGraph(nop);
 TGraph *av16 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av16->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av16->GetPoint(j, xpt, ypt);
   av16->SetPoint(j, xpt, ypt + fpc[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av17 = new TGraph(nop);
 TGraph *av17 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av17->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av17->GetPoint(j, xpt, ypt);
   av17->SetPoint(j, xpt, ypt + cbs[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av18 = new TGraph(nop);
 TGraph *av18 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av18->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av18->GetPoint(j, xpt, ypt);
   av18->SetPoint(j, xpt, ypt + cbc[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av19 = new TGraph(nop);
 TGraph *av19 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av19->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av19->GetPoint(j, xpt, ypt);
   av19->SetPoint(j, xpt, ypt + ces[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av20 = new TGraph(nop);
 TGraph *av20 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av20->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av20->GetPoint(j, xpt, ypt);
   av20->SetPoint(j, xpt, ypt + cec[i]->Eval(xpt)/(Double_t)nof);
  }
 }

//=== generating the functions ===//
//========================Theta=43.0

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  //choice: transfer momentum square; electron(finite)proton
  fpe[i] = new TF1("fna","(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[3]*[3])/(2.0*[3]*[3]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((4.0*[3]*[3])/(4.0*[3]*[3]+x))*(((0.71)/(0.71+x))*((0.71)/(0.71+x))*((0.71)/(0.71+x))*((0.71)/(0.71+x))*([7]*[7]*[8]*[8]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+[6]*[6]*[7]*[7]*(2.0*[3]*[3]*x*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(8.0*[3]*[3]*[3]*[3]))))", 0.0025, 1.0);
  fpe[i]->SetParameter(0,1.0);
  fpe[i]->SetParameter(1,1.0);
  fpe[i]->SetParameter(2,1.0/137.0);
  fpe[i]->SetParameter(3,gRandom->Gaus(0.938272081,0.000000006));
  fpe[i]->SetParameter(4,43.0);
  fpe[i]->SetParameter(5,TMath::Pi());
  fpe[i]->SetParameter(6,2.792847351);
  fpe[i]->SetParameter(7,gRandom->Gaus(1.010,0.005)+gRandom->Gaus(0.0,0.005));
  fpe[i]->SetParameter(8,gRandom->Gaus(0.973,0.014)+gRandom->Gaus(0.0,0.005));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 2/3*(e-u) and 1/3*(e-d) bare mass with modified recoil factor
  ebs[i] = new TF1("fna","(2.0/3.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(1.0/3.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  ebs[i]->SetParameter(0,1.0);
  ebs[i]->SetParameter(1,2.0/3.0);
  ebs[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   ebs[i]->SetParameter(3,gRandom->Gaus(0.0017,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        ebs[i]->SetParameter(3,gRandom->Gaus(0.0033,0.0008));
                                       }
                                       else {
                                             ebs[i]->SetParameter(3,(gRandom->Gaus(0.0017,0.0008)+gRandom->Gaus(0.0033,0.0008))/2.0);
					    }
  ebs[i]->SetParameter(4,43.0);
  ebs[i]->SetParameter(5,TMath::Pi());
  ebs[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   ebs[i]->SetParameter(7,gRandom->Gaus(0.0041,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        ebs[i]->SetParameter(7,gRandom->Gaus(0.0058,0.0008));
                                       }
                                       else {
                                             ebs[i]->SetParameter(7,(gRandom->Gaus(0.0041,0.0008)+gRandom->Gaus(0.0058,0.0008))/2.0);
                                            }
  ebs[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000023));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 72/81*(e-u) and 9/81*(e-d) bare mass with modified recoil factor
  ebc[i] = new TF1("fna","(72.0/81.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(9.0/81.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  ebc[i]->SetParameter(0,1.0);
  ebc[i]->SetParameter(1,2.0/3.0);
  ebc[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   ebc[i]->SetParameter(3,gRandom->Gaus(0.0017,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        ebc[i]->SetParameter(3,gRandom->Gaus(0.0033,0.0008));
                                       }
                                       else {
                                             ebc[i]->SetParameter(3,(gRandom->Gaus(0.0017,0.0008)+gRandom->Gaus(0.0033,0.0008))/2.0);
                                            }
  ebc[i]->SetParameter(4,43.0);
  ebc[i]->SetParameter(5,TMath::Pi());
  ebc[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   ebc[i]->SetParameter(7,gRandom->Gaus(0.0041,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        ebc[i]->SetParameter(7,gRandom->Gaus(0.0058,0.0008));
                                       }
                                       else {
                                             ebc[i]->SetParameter(7,(gRandom->Gaus(0.0041,0.0008)+gRandom->Gaus(0.0058,0.0008))/2.0);
                                            }
  ebc[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000023));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 2/3*(e-u) and 1/3*(e-d) effective mass with modified recoil factor
  ees[i] = new TF1("fna","(2.0/3.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(1.0/3.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  ees[i]->SetParameter(0,1.0);
  ees[i]->SetParameter(1,2.0/3.0);
  ees[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   ees[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                  }
                  else if (rdm > 0.75) {
                                        ees[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                       }
                                       else {
                                             ees[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                            }
  ees[i]->SetParameter(4,43.0);
  ees[i]->SetParameter(5,TMath::Pi());
  ees[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   ees[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                  }
                  else if (rdm > 0.75) {
                                        ees[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                       }
					else {
                                             ees[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                            }
  ees[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000006));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 72/81*(e-u) and 9/81*(e-d) effective mass with modified recoil factor
  eec[i] = new TF1("fna","(72.0/81.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(9.0/81.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  eec[i]->SetParameter(0,1.0);
  eec[i]->SetParameter(1,2.0/3.0);
  eec[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   eec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                  }
                  else if (rdm > 0.75) {
                                        eec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                       }
                                       else {
                                             eec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                            }
  eec[i]->SetParameter(4,43.0);
  eec[i]->SetParameter(5,TMath::Pi());
  eec[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   eec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                  }
                  else if (rdm > 0.75) {
                                        eec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                       }
                                       else {
                                             eec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                            }
  eec[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000006));
 }

//=== getting the average ===//

 //auto av21 = new TGraph(nop);
 TGraph *av21 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av21->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av21->GetPoint(j, xpt, ypt);
   av21->SetPoint(j, xpt, ypt + fpe[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av22 = new TGraph(nop);
 TGraph *av22 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av22->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av22->GetPoint(j, xpt, ypt);
   av22->SetPoint(j, xpt, ypt + ebs[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av23 = new TGraph(nop);
 TGraph *av23 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av23->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av23->GetPoint(j, xpt, ypt);
   av23->SetPoint(j, xpt, ypt + ebc[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av24 = new TGraph(nop);
 TGraph *av24 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av24->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av24->GetPoint(j, xpt, ypt);
   av24->SetPoint(j, xpt, ypt + ees[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av25 = new TGraph(nop);
 TGraph *av25 = new TGraph();

 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av25->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av25->GetPoint(j, xpt, ypt);
   av25->SetPoint(j, xpt, ypt + eec[i]->Eval(xpt)/(Double_t)nof);
  }
 }

//=== generating the functions ===//
//========================Theta=33.0

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  //choice: transfer momentum square; electron(finite)proton
  fpg[i] = new TF1("fna","(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[3]*[3])/(2.0*[3]*[3]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((4.0*[3]*[3])/(4.0*[3]*[3]+x))*(((0.71)/(0.71+x))*((0.71)/(0.71+x))*((0.71)/(0.71+x))*((0.71)/(0.71+x))*([7]*[7]*[8]*[8]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+[6]*[6]*[7]*[7]*(2.0*[3]*[3]*x*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(8.0*[3]*[3]*[3]*[3]))))", 0.0025, 1.0);
  fpg[i]->SetParameter(0,1.0);
  fpg[i]->SetParameter(1,1.0);
  fpg[i]->SetParameter(2,1.0/137.0);
  fpg[i]->SetParameter(3,gRandom->Gaus(0.938272081,0.000000006));
  fpg[i]->SetParameter(4,33.0);
  fpg[i]->SetParameter(5,TMath::Pi());
  fpg[i]->SetParameter(6,2.792847351);
  fpg[i]->SetParameter(7,gRandom->Gaus(1.017,0.004)+gRandom->Gaus(0.0,0.007));
  fpg[i]->SetParameter(8,gRandom->Gaus(0.994,0.013)+gRandom->Gaus(0.0,0.007));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 2/3*(e-u) and 1/3*(e-d) bare mass with modified recoil factor
  gbs[i] = new TF1("fna","(2.0/3.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(1.0/3.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  gbs[i]->SetParameter(0,1.0);
  gbs[i]->SetParameter(1,2.0/3.0);
  gbs[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   gbs[i]->SetParameter(3,gRandom->Gaus(0.0017,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        gbs[i]->SetParameter(3,gRandom->Gaus(0.0033,0.0008));
                                       }
                                       else {
                                             gbs[i]->SetParameter(3,(gRandom->Gaus(0.0017,0.0008)+gRandom->Gaus(0.0033,0.0008))/2.0);
					    }
  gbs[i]->SetParameter(4,33.0);
  gbs[i]->SetParameter(5,TMath::Pi());
  gbs[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   gbs[i]->SetParameter(7,gRandom->Gaus(0.0041,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        gbs[i]->SetParameter(7,gRandom->Gaus(0.0058,0.0008));
                                       }
                                       else {
                                             gbs[i]->SetParameter(7,(gRandom->Gaus(0.0041,0.0008)+gRandom->Gaus(0.0058,0.0008))/2.0);
                                            }
  gbs[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000023));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 72/81*(e-u) and 9/81*(e-d) bare mass with modified recoil factor
  gbc[i] = new TF1("fna","(72.0/81.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(9.0/81.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  gbc[i]->SetParameter(0,1.0);
  gbc[i]->SetParameter(1,2.0/3.0);
  gbc[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   gbc[i]->SetParameter(3,gRandom->Gaus(0.0017,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        gbc[i]->SetParameter(3,gRandom->Gaus(0.0033,0.0008));
                                       }
                                       else {
                                             gbc[i]->SetParameter(3,(gRandom->Gaus(0.0017,0.0008)+gRandom->Gaus(0.0033,0.0008))/2.0);
                                            }
  gbc[i]->SetParameter(4,33.0);
  gbc[i]->SetParameter(5,TMath::Pi());
  gbc[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   gbc[i]->SetParameter(7,gRandom->Gaus(0.0041,0.0008));
                  }
                  else if (rdm > 0.75) {
                                        gbc[i]->SetParameter(7,gRandom->Gaus(0.0058,0.0008));
                                       }
                                       else {
                                             gbc[i]->SetParameter(7,(gRandom->Gaus(0.0041,0.0008)+gRandom->Gaus(0.0058,0.0008))/2.0);
                                            }
  gbc[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000023));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 2/3*(e-u) and 1/3*(e-d) effective mass with modified recoil factor
  ges[i] = new TF1("fna","(2.0/3.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(1.0/3.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  ges[i]->SetParameter(0,1.0);
  ges[i]->SetParameter(1,2.0/3.0);
  ges[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   ges[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                  }
                  else if (rdm > 0.75) {
                                        ges[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                       }
                                       else {
                                             ges[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                            }
  ges[i]->SetParameter(4,33.0);
  ges[i]->SetParameter(5,TMath::Pi());
  ges[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   ges[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                  }
                  else if (rdm > 0.75) {
                                        ges[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                       }
					else {
                                             ges[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                            }
  ges[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000006));
 }

 for (Int_t i=0; i<nof; i++)
 {
  fna = "f"; fna += i;
  // 72/81*(e-u) and 9/81*(e-d) effective mass with modified recoil factor
  gec[i] = new TF1("fna","(72.0/81.0)*(((2.0*[0]*[0]*[1]*[1]*[2]*[2]*[3]*[3])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[3]*[3]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[3]*[3]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[3]*[3])))+(9.0/81.0)*(((2.0*[0]*[0]*[6]*[6]*[2]*[2]*[7]*[7])/(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+2.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+(x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[7]*[7]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)))))*((2.0*[8]*[8])/(2.0*[8]*[8]+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+sqrt(x*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)*sin([4]*[5]/360.0)+4.0*[8]*[8]*x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))))*((2.0*[7]*[7]*cos([4]*[5]/360.0)*cos([4]*[5]/360.0)+x*sin([4]*[5]/360.0)*sin([4]*[5]/360.0))/(2.0*[7]*[7])))", 0.0025, 1.0);
  gec[i]->SetParameter(0,1.0);
  gec[i]->SetParameter(1,2.0/3.0);
  gec[i]->SetParameter(2,1.0/137.0);
  if (rdm < 0.25) {
                   gec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                  }
                  else if (rdm > 0.75) {
                                        gec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                       }
                                       else {
                                             gec[i]->SetParameter(3,gRandom->Gaus(0.336,0.010));
                                            }
  gec[i]->SetParameter(4,33.0);
  gec[i]->SetParameter(5,TMath::Pi());
  gec[i]->SetParameter(6,-1.0/3.0);
  if (rdm < 0.25) {
                   gec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                  }
                  else if (rdm > 0.75) {
                                        gec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                       }
                                       else {
                                             gec[i]->SetParameter(7,gRandom->Gaus(0.340,0.010));
                                            }
  gec[i]->SetParameter(8,gRandom->Gaus(0.938272081,0.000000006));
 }

//=== getting the average ===//

 //auto av26 = new TGraph(nop);
 TGraph *av26 = new TGraph();

 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av26->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av26->GetPoint(j, xpt, ypt);
   av26->SetPoint(j, xpt, ypt + fpg[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av27 = new TGraph(nop);
 TGraph *av27 = new TGraph();

 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av27->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av27->GetPoint(j, xpt, ypt);
   av27->SetPoint(j, xpt, ypt + gbs[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av28 = new TGraph(nop);
 TGraph *av28 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av28->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av28->GetPoint(j, xpt, ypt);
   av28->SetPoint(j, xpt, ypt + gbc[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av29 = new TGraph(nop);
 TGraph *av29 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av29->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av29->GetPoint(j, xpt, ypt);
   av29->SetPoint(j, xpt, ypt + ges[i]->Eval(xpt)/(Double_t)nof);
  }
 }

 //auto av30 = new TGraph(nop);
 TGraph *av30 = new TGraph();
 for (Int_t j=0; j < nop; j++)
 {
  xpt = del*j + 0.0025*del;
  av30->SetPoint(j, xpt, 0.0);
 }
 for (Int_t i=0; i<nof; i++)
 {
  for (Int_t j=0; j < nop; j++)
  {
   av30->GetPoint(j, xpt, ypt);
   av30->SetPoint(j, xpt, ypt + gec[i]->Eval(xpt)/(Double_t)nof);
  }
 }

//======================

 av1->SetMarkerStyle(7);
 av1->SetMarkerColor(kBlack);
 av1->SetMarkerSize(1);

 av2->SetMarkerStyle(7);
 av2->SetMarkerColor(kRed);
 av2->SetMarkerSize(1);

 av3->SetMarkerStyle(7);
 av3->SetMarkerColor(kGreen);
 av3->SetMarkerSize(1);

 av4->SetMarkerStyle(7);
 av4->SetMarkerColor(kBlue);
 av4->SetMarkerSize(1);

 av5->SetMarkerStyle(7);
 av5->SetMarkerColor(kCyan);
 av5->SetMarkerSize(1);

 av6->SetMarkerStyle(7);
 av6->SetMarkerColor(kBlack);
 av6->SetMarkerSize(1);

 av7->SetMarkerStyle(7);
 av7->SetMarkerColor(kRed);
 av7->SetMarkerSize(1);

 av8->SetMarkerStyle(7);
 av8->SetMarkerColor(kGreen);
 av8->SetMarkerSize(1);

 av9->SetMarkerStyle(7);
 av9->SetMarkerColor(kBlue);
 av9->SetMarkerSize(1);

 av10->SetMarkerStyle(7);
 av10->SetMarkerColor(kCyan);
 av10->SetMarkerSize(1);

 av11->SetMarkerStyle(7);
 av11->SetMarkerColor(kBlack);
 av11->SetMarkerSize(1);

 av12->SetMarkerStyle(7);
 av12->SetMarkerColor(kRed);
 av12->SetMarkerSize(1);

 av13->SetMarkerStyle(7);
 av13->SetMarkerColor(kGreen);
 av13->SetMarkerSize(1);

 av14->SetMarkerStyle(7);
 av14->SetMarkerColor(kBlue);
 av14->SetMarkerSize(1);

 av15->SetMarkerStyle(7);
 av15->SetMarkerColor(kCyan);
 av15->SetMarkerSize(1);

 av16->SetMarkerStyle(7);
 av16->SetMarkerColor(kBlack);
 av16->SetMarkerSize(1);

 av17->SetMarkerStyle(7);
 av17->SetMarkerColor(kRed);
 av17->SetMarkerSize(1);

 av18->SetMarkerStyle(7);
 av18->SetMarkerColor(kGreen);
 av18->SetMarkerSize(1);

 av19->SetMarkerStyle(7);
 av19->SetMarkerColor(kBlue);
 av19->SetMarkerSize(1);

 av20->SetMarkerStyle(7);
 av20->SetMarkerColor(kCyan);
 av20->SetMarkerSize(1);

 av21->SetMarkerStyle(7);
 av21->SetMarkerColor(kBlack);
 av21->SetMarkerSize(1);

 av22->SetMarkerStyle(7);
 av22->SetMarkerColor(kRed);
 av22->SetMarkerSize(1);

 av23->SetMarkerStyle(7);
 av23->SetMarkerColor(kGreen);
 av23->SetMarkerSize(1);

 av24->SetMarkerStyle(7);
 av24->SetMarkerColor(kBlue);
 av24->SetMarkerSize(1);

 av25->SetMarkerStyle(7);
 av25->SetMarkerColor(kCyan);
 av25->SetMarkerSize(1);

 av26->SetMarkerStyle(7);
 av26->SetMarkerColor(kBlack);
 av26->SetMarkerSize(1);

 av27->SetMarkerStyle(7);
 av27->SetMarkerColor(kRed);
 av27->SetMarkerSize(1);

 av28->SetMarkerStyle(7);
 av28->SetMarkerColor(kGreen);
 av28->SetMarkerSize(1);

 av29->SetMarkerStyle(7);
 av29->SetMarkerColor(kBlue);
 av29->SetMarkerSize(1);

 av30->SetMarkerStyle(7);
 av30->SetMarkerColor(kCyan);
 av30->SetMarkerSize(1);

//=== Plotting the points on canvas ===//

TGraph *av0000 = new TGraph();
Int_t j = 0;

//================== dcs of ep 

 for (Int_t i=0; i<nof; i++)
 {
   av1->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399006) < 0.00001)
   {
    av0000->GetPoint(i, xpt, ypt);
    av0000->SetPoint(j, xpt, ypt + fpf[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av6->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.498756) < 0.00001)
   {
    av0000->GetPoint(i, xpt, ypt);
    av0000->SetPoint(j, xpt, ypt + fpa[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av11->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.601) < 0.001)
   {
    av0000->GetPoint(i, xpt, ypt);
    av0000->SetPoint(j, xpt, ypt + fpb[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

TGraph *av0001 = new TGraph();
//====================Q^2=0.50
double q21 = 0.498758;
double ept = 0.00001;
double eps = 0.00001;

 for (Int_t i=0; i<nof; i++)
 {
   av21->GetPoint(i, xpt, ypt);
   if ( abs(xpt - q21) < ept )
   {
    av0001->GetPoint(i, xpt, ypt);
    av0001->SetPoint(j, xpt, ypt + fpe[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av26->GetPoint(i, xpt, ypt);
   if ( abs(xpt - q21) < ept )
   {
    av0001->GetPoint(i, xpt, ypt);
    av0001->SetPoint(j, xpt, ypt + fpg[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

TGraph *av0002 = new TGraph();
//===============================Q^2=0.60
double q22 = 0.59851;
double epu = 0.00001;

 for (Int_t i=0; i<nof; i++)
 {
   av21->GetPoint(i, xpt, ypt);
   if ( abs(xpt - q22) < epu )
   {
    av0001->GetPoint(i, xpt, ypt);
    av0001->SetPoint(j, xpt, ypt + fpe[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av26->GetPoint(i, xpt, ypt);
   if ( abs(xpt - q22) < epu )
   {
    av0001->GetPoint(i, xpt, ypt);
    av0001->SetPoint(j, xpt, ypt + fpg[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

//==========================
TGraph *av001 = new TGraph();

 for (Int_t i=0; i<nof; i++)
 {
   av2->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av001->GetPoint(i, xpt, ypt);
    av001->SetPoint(j, xpt, ypt + qbs[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av7->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av001->GetPoint(i, xpt, ypt);
    av001->SetPoint(j, xpt, ypt + dbs[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av12->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av001->GetPoint(i, xpt, ypt);
    av001->SetPoint(j, xpt, ypt + bbs[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av17->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av001->GetPoint(i, xpt, ypt);
    av001->SetPoint(j, xpt, ypt + cbs[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

//==========================
TGraph *av002 = new TGraph();

 for (Int_t i=0; i<nof; i++)
 {
   av3->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av002->GetPoint(i, xpt, ypt);
    av002->SetPoint(j, xpt, ypt + qbc[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av8->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av002->GetPoint(i, xpt, ypt);
    av002->SetPoint(j, xpt, ypt + abc[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av13->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av002->GetPoint(i, xpt, ypt);
    av002->SetPoint(j, xpt, ypt + bbc[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av18->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av002->GetPoint(i, xpt, ypt);
    av002->SetPoint(j, xpt, ypt + cbc[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }
//==========================
TGraph *av003 = new TGraph();

 for (Int_t i=0; i<nof; i++)
 {
   av4->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av003->GetPoint(i, xpt, ypt);
    av003->SetPoint(j, xpt, ypt + qes[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av9->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av003->GetPoint(i, xpt, ypt);
    av003->SetPoint(j, xpt, ypt + aes[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av14->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av003->GetPoint(i, xpt, ypt);
    av003->SetPoint(j, xpt, ypt + bes[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av19->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av003->GetPoint(i, xpt, ypt);
    av003->SetPoint(j, xpt, ypt + ces[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

//==========================
TGraph *av004 = new TGraph();

 for (Int_t i=0; i<nof; i++)
 {
   av5->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av004->GetPoint(i, xpt, ypt);
    av004->SetPoint(j, xpt, ypt + qec[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av10->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av004->GetPoint(i, xpt, ypt);
    av004->SetPoint(j, xpt, ypt + aec[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av15->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av004->GetPoint(i, xpt, ypt);
    av004->SetPoint(j, xpt, ypt + bec[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

 for (Int_t i=0; i<nof; i++)
 {
   av20->GetPoint(i, xpt, ypt);
   if ( abs(xpt - 0.399012) < eps)
   {
    av004->GetPoint(i, xpt, ypt);
    av004->SetPoint(j, xpt, ypt + cec[i]->Eval(xpt)/(Double_t)nof);
    j++;
   }
 }

//=== drawing the canvas ===//

 gPad->SetLogy(1);

 av0000->SetMarkerStyle(22);
 av0000->SetMarkerColor(3); //green
 av0000->SetMarkerSize(1);

 av0001->SetMarkerStyle(22);
 av0001->SetMarkerColor(1); //black
 av0001->SetMarkerSize(1);

 av001->SetMarkerStyle(22);
 av001->SetMarkerColor(3); //green
 av001->SetMarkerSize(1);

 av002->SetMarkerStyle(22);
 av002->SetMarkerColor(2); //red
 av002->SetMarkerSize(1);

 av003->SetMarkerStyle(22);
 av003->SetMarkerColor(25); //brown
 av003->SetMarkerSize(1);

 av004->SetMarkerStyle(22);
 av004->SetMarkerColor(38); //blue-violet
 av004->SetMarkerSize(1);

//===Unified multigraph===//

 TMultiGraph *mgp = new TMultiGraph();
 mgp->SetTitle("dcs of ep, eq_s, eq_c (with mrf); Transfer Momentum in GeV^2; Differential Cross Section");
 mgp->SetMinimum(1.e-5);
 mgp->SetMaximum(1.e0);
 mgp->Add(av0000,"PL");
 //mgp->Add(av0001,"P");
 //mgp->Add(av001,"P");
 //mgp->Add(av002,"P");
 //mgp->Add(av003,"P");
 //mgp->Add(av004,"P");
 //mgp->Add(av1,"PL");
 //mgp->Add(av2,"PL");
 //mgp->Add(av3,"PL");
 //mgp->Add(av4,"PL");
 //mgp->Add(av5,"PL");
 //mgp->Add(av6,"PL");
 //mgp->Add(av7,"PL");
 //mgp->Add(av8,"PL");
 //mgp->Add(av9,"PL");
 //mgp->Add(av10,"PL");
 //mgp->Add(av11,"PL");
 //mgp->Add(av12,"PL");
 //mgp->Add(av13,"PL");
 //mgp->Add(av14,"PL");
 //mgp->Add(av15,"PL");
 //mgp->Add(av16,"PL");
 //mgp->Add(av17,"PL");
 //mgp->Add(av18,"PL");
 //mgp->Add(av19,"PL");
 //mgp->Add(av20,"PL");
 //mgp->Add(av21,"PL");
 //mgp->Add(av22,"PL");
 //mgp->Add(av23,"PL");
 //mgp->Add(av24,"PL");
 //mgp->Add(av25,"PL");
 //mgp->Add(av26,"PL");
 //mgp->Add(av27,"PL");
 //mgp->Add(av28,"PL");
 //mgp->Add(av29,"PL");
 //mgp->Add(av30,"PL");

 mgp->Draw("a");

 avf->SetMaximum(1.e0);
 can->Modified();
 can->Update();

 av1->Print();
 av6->Print();
 av11->Print();
 av16->Print();
 av21->Print();
 av26->Print();
}