#ifndef __CINT__

#include "TFile.h"
#include "TProfile.h"
#include "TRandom.h"
#include "TSystem.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TDirectory.h"
#include <TVirtualFitter.h>
#include <TH1F.h>
#include <TStyle.h>
#include <TLatex.h>
#include <TText.h>
#include <TF1.h>
#include <TF2.h>
#include <TROOT.h>

#endif

#include "TCanvas.h"
#include <TFumili.h>
#include <TGraph.h>

//#define FUMILI 

void FCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
double fitfun(Double_t *xx, Double_t *par);

  double ch2;

const int NSPEC=14;
const int NFITPAR=10;
//const int NFITPAR=6;
TGraph *cutg[NSPEC];

float Z[NSPEC]={1,1,1,2,2,2,2,3,3,3,4, 5, 6, 7};
float A[NSPEC]={1,2,3,3,4,6,8,6,7,8,7,11,13,15};
float AZ2[NSPEC];


void fit_rl_fcn_inv_parl()
{
gROOT->Reset();
gStyle->SetOptTitle(0);
gStyle->SetOptStat(0);
gStyle->SetOptFit(0);

#ifdef FUMILI
TVirtualFitter::SetDefaultFitter("Fumili");
#else
TVirtualFitter::SetDefaultFitter("Minuit");
#endif
TVirtualFitter *fitter=TVirtualFitter::Fitter(0,NFITPAR);


for(int i=0;i<NSPEC;i++){
  AZ2[i]=A[i]*Z[i]*Z[i];
}

//                       F        S
   cutg[0] = new TGraph(10);
   cutg[0]->SetPoint(0,52.1849,22.5474);
   cutg[0]->SetPoint(1,92.0483,62.2827);
   cutg[0]->SetPoint(2,155.467,124.196);
   cutg[0]->SetPoint(3,218.887,198.122);
   cutg[0]->SetPoint(4,296.801,289.605);
   cutg[0]->SetPoint(5,400.084,408.811);
   cutg[0]->SetPoint(6,479.811,498.446);
   cutg[0]->SetPoint(7,550.478,583.461);
   cutg[0]->SetPoint(8,628.393,674.945);
   cutg[0]->SetPoint(9,700.872,759.959);

   cutg[1] = new TGraph(11);
   cutg[1]->SetPoint(0,104.474,59.3203);
   cutg[1]->SetPoint(1,152.046,101.809);
   cutg[1]->SetPoint(2,241.699,190.054);
   cutg[1]->SetPoint(3,338.672,294.641);
   cutg[1]->SetPoint(4,428.326,392.691);
   cutg[1]->SetPoint(5,527.128,503.814);
   cutg[1]->SetPoint(6,635.079,628.011);
   cutg[1]->SetPoint(7,730.221,739.135);
   cutg[1]->SetPoint(8,810.727,835.551);
   cutg[1]->SetPoint(9,902.21,941.772);
   cutg[1]->SetPoint(10,962.59,1016.94);
   
   cutg[2] = new TGraph(12);
   cutg[2]->SetPoint(0,238.487,159.311);
   cutg[2]->SetPoint(1,314.886,229.274);
   cutg[2]->SetPoint(2,393.562,310.982);
   cutg[2]->SetPoint(3,483.216,410.667);
   cutg[2]->SetPoint(4,563.721,500.546);
   cutg[2]->SetPoint(5,668.013,614.938);
   cutg[2]->SetPoint(6,748.518,709.72 );
   cutg[2]->SetPoint(7,823.535,794.696);
   cutg[2]->SetPoint(8,927.826,918.893);
   cutg[2]->SetPoint(9,1022.97,1026.75);
   cutg[2]->SetPoint(10,1099.81,1118.26);
   cutg[2]->SetPoint(11,1158.36,1186.9);  

   cutg[3] = new TGraph(15);
   cutg[3]->SetPoint(0,473.603,218.215);
   cutg[3]->SetPoint(1,667.432,357.819);
   cutg[3]->SetPoint(2,819.611,476.956);
   cutg[3]->SetPoint(3,1051.89,673.57 );
   cutg[3]->SetPoint(4,1279.35,876.032);
   cutg[3]->SetPoint(5,1482.79,1064.61);
   cutg[3]->SetPoint(6,1660.6 ,1242.22);
   cutg[3]->SetPoint(7,1841.62,1420.56);
   cutg[3]->SetPoint(8,2014.62,1603.28);
   cutg[3]->SetPoint(9,2195.64,1786.01);
   cutg[3]->SetPoint(10,2322.18,1924.88);
   cutg[3]->SetPoint(11,2447.13,2061.56);
   cutg[3]->SetPoint(12,2586.5 ,2218.71);
   cutg[3]->SetPoint(13,2735.47,2380.24);
   cutg[3]->SetPoint(14,2900.47,2565.16);
  
   cutg[4] = new TGraph(18);
   cutg[4]->SetPoint(0,186.562,47.7042);   
   cutg[4]->SetPoint(1,345.937,121.474);   
   cutg[4]->SetPoint(2,528.75 ,224.751);    
   cutg[4]->SetPoint(3,716.25 ,345.241);    
   cutg[4]->SetPoint(4,955.312,522.287);   
   cutg[4]->SetPoint(5,1170.94,699.334);   
   cutg[4]->SetPoint(6,1353.75,854.25 );    
   cutg[4]->SetPoint(7,1541.25,1016.54);   
   cutg[4]->SetPoint(8,1663.12,1132.11);   
   cutg[4]->SetPoint(9,1822.5 ,1277.19);    
   cutg[4]->SetPoint(10,2019.38,1468.99);  
   cutg[4]->SetPoint(11,2197.5 ,1643.58);   
   cutg[4]->SetPoint(12,2370.94,1818.17);  
   cutg[4]->SetPoint(13,2549.06,2009.97);  
   cutg[4]->SetPoint(14,2717.81,2187.02);  
   cutg[4]->SetPoint(15,2872.5 ,2349.31);   
   cutg[4]->SetPoint(16,3041.25,2528.81);  
   cutg[4]->SetPoint(17,3228.75,2737.83); 
  
   cutg[5] = new TGraph(20);
   cutg[5]->SetPoint(0, 356.36 ,100.172);
   cutg[5]->SetPoint(1, 455.717,151.617);
   cutg[5]->SetPoint(2, 654.43 ,260.444);
   cutg[5]->SetPoint(3, 788.346,339.591);
   cutg[5]->SetPoint(4, 1004.34,472.162);
   cutg[5]->SetPoint(5, 1190.09,610.669);
   cutg[5]->SetPoint(6, 1384.49,743.239);
   cutg[5]->SetPoint(7, 1557.28,879.768);
   cutg[5]->SetPoint(8, 1768.95,1051.91);
   cutg[5]->SetPoint(9, 1976.31,1224.06);
   cutg[5]->SetPoint(10,2183.66,1408.07);
   cutg[5]->SetPoint(11,2382.37,1584.17);
   cutg[5]->SetPoint(12,2602.68,1784.02);
   cutg[5]->SetPoint(13,2784.12,1970.01);
   cutg[5]->SetPoint(14,2969.87,2140.18);
   cutg[5]->SetPoint(15,3155.63,2332.11);
   cutg[5]->SetPoint(16,3371.62,2537.89);
   cutg[5]->SetPoint(17,3583.29,2763.46);
   cutg[5]->SetPoint(18,3855.44,3050.37);
   cutg[5]->SetPoint(19,4110.31,3321.45);

   cutg[6] = new TGraph(11);
   cutg[6]->SetPoint(0 ,1099.12,483.048);
   cutg[6]->SetPoint(1 ,1235.61,570.318);
   cutg[6]->SetPoint(2 ,1392.78,674.211);
   cutg[6]->SetPoint(3 ,1607.85,821.74 );
   cutg[6]->SetPoint(4 ,1822.92,977.58 );
   cutg[6]->SetPoint(5 ,2112.44,1204.07);
   cutg[6]->SetPoint(6 ,2373.01,1407.7 );
   cutg[6]->SetPoint(7 ,2608.77,1613.41);
   cutg[6]->SetPoint(8 ,3113.36,2047.68);
   cutg[6]->SetPoint(9 ,3514.56,2434.16);
   cutg[6]->SetPoint(10,3944.71,2849.74);

   cutg[7] = new TGraph(15);
   cutg[7]->SetPoint(0,1420.96,587.416);
   cutg[7]->SetPoint(1,1801.47,808.954);
   cutg[7]->SetPoint(2,2159.93,1030.49);
   cutg[7]->SetPoint(3,2573.53,1301.76);
   cutg[7]->SetPoint(4,2970.59,1566.25);
   cutg[7]->SetPoint(5,3362.13,1853.35);
   cutg[7]->SetPoint(6,3748.16,2144.96);
   cutg[7]->SetPoint(7,4134.19,2438.84);
   cutg[7]->SetPoint(8,4586.4 ,2805.06);
   cutg[7]->SetPoint(9,4966.91,3126.06);
   cutg[7]->SetPoint(10,5363.97,3465.15);
   cutg[7]->SetPoint(11,5788.6 ,3831.36);
   cutg[7]->SetPoint(12,6086.4 ,4095.85);
   cutg[7]->SetPoint(13,6345.59,4337.74);
   cutg[7]->SetPoint(14,6643.38,4611.27); 
  
   cutg[8] = new TGraph(22);
   cutg[8]->SetPoint(0,459.727,121.082);  
   cutg[8]->SetPoint(1,728.503,229.846);  
   cutg[8]->SetPoint(2,1044.02,375.302);  
   cutg[8]->SetPoint(3,1373.57,533.862);  
   cutg[8]->SetPoint(4,1616.63,667.524);  
   cutg[8]->SetPoint(5,2013.95,882.431);  
   cutg[8]->SetPoint(6,2376.22,1109.13);  
   cutg[8]->SetPoint(7,2726.8 ,1330.59);   
   cutg[8]->SetPoint(8,3016.61,1528.46);  
   cutg[8]->SetPoint(9,3322.78,1743.37);  
   cutg[8]->SetPoint(10,3722.44,2026.42); 
   cutg[8]->SetPoint(11,4122.1 ,2326.5 );   
   cutg[8]->SetPoint(12,4458.65,2579.41); 
   cutg[8]->SetPoint(13,4823.25,2862.46); 
   cutg[8]->SetPoint(14,5155.13,3120.61); 
   cutg[8]->SetPoint(15,5423.91,3348.63); 
   cutg[8]->SetPoint(16,5741.76,3619.88); 
   cutg[8]->SetPoint(17,6092.34,3916.03); 
   cutg[8]->SetPoint(18,6438.25,4214.81); 
   cutg[8]->SetPoint(19,6749.09,4489.99); 
   cutg[8]->SetPoint(20,7066.95,4776.97); 
   cutg[8]->SetPoint(21,7380.13,5041.68); 

   cutg[9] = new TGraph(16);
   cutg[9]->SetPoint(0 ,1936.87,794.358);
   cutg[9]->SetPoint(1 ,2196.3 ,939.203);
   cutg[9]->SetPoint(2 ,2476.76,1106.77);
   cutg[9]->SetPoint(3 ,2764.23,1291.37);
   cutg[9]->SetPoint(4 ,3086.76,1498.7 );
   cutg[9]->SetPoint(5 ,3430.33,1728.75);
   cutg[9]->SetPoint(6 ,3872.06,2032.64);
   cutg[9]->SetPoint(7 ,4271.72,2310.97);
   cutg[9]->SetPoint(8 ,4685.4 ,2626.22);
   cutg[9]->SetPoint(9 ,5078.05,2918.75);
   cutg[9]->SetPoint(10,5470.69,3222.64);
   cutg[9]->SetPoint(11,5870.35,3537.89);
   cutg[9]->SetPoint(12,6277.02,3873.02);
   cutg[9]->SetPoint(13,6725.76,4253.59);
   cutg[9]->SetPoint(14,7111.4 ,4602.93);
   cutg[9]->SetPoint(15,7419.91,4872.74);
  
   cutg[10] = new TGraph(23);
   cutg[10]->SetPoint(0,1128   ,363.952);
   cutg[10]->SetPoint(1,1331.19,447.723);
   cutg[10]->SetPoint(2,1629.43,579.179);
   cutg[10]->SetPoint(3,1966.99,744.143);
   cutg[10]->SetPoint(4,2320.94,915.551);
   cutg[10]->SetPoint(5,2691.28,1114.02);
   cutg[10]->SetPoint(6,3150.1 ,1360.18);
   cutg[10]->SetPoint(7,3477.83,1553.5 );
   cutg[10]->SetPoint(8,3936.66,1812.54);
   cutg[10]->SetPoint(9,4270.94,2014.88);
   cutg[10]->SetPoint(10,4641.28,2239.13);
   cutg[10]->SetPoint(11,5047.66,2489.16);
   cutg[10]->SetPoint(12,5418   ,2732.74);
   cutg[10]->SetPoint(13,5814.55,2989.2 );
   cutg[10]->SetPoint(14,6158.67,3217.32);
   cutg[10]->SetPoint(15,6522.45,3471.21);
   cutg[10]->SetPoint(16,6945.23,3768.92);
   cutg[10]->SetPoint(17,7430.27,4111.73);
   cutg[10]->SetPoint(18,7961.19,4506.1 );
   cutg[10]->SetPoint(19,8446.24,4861.8 );
   cutg[10]->SetPoint(20,8957.5 ,5235.55);
   cutg[10]->SetPoint(21,9458.92,5632.5 );
   cutg[10]->SetPoint(22,9917.75,5992.07);


   cutg[11] = new TGraph(18);
   cutg[11]->SetPoint(0 ,710.084,143.646 );
   cutg[11]->SetPoint(1 ,1235.29,318.168 );
   cutg[11]->SetPoint(2 ,2117.65,667.214 );
   cutg[11]->SetPoint(3 ,3084.03,1069.96 );
   cutg[11]->SetPoint(4 ,3861.34,1425.72 );
   cutg[11]->SetPoint(5 ,4554.62,1741.2  );
   cutg[11]->SetPoint(6 ,5415.97,2177.51 );
   cutg[11]->SetPoint(7 ,6193.28,2560.11 );
   cutg[11]->SetPoint(8 ,6886.55,2942.72 );
   cutg[11]->SetPoint(9 ,7558.82,3298.48 );
   cutg[11]->SetPoint(10,8084.03,3607.25);
   cutg[11]->SetPoint(11,8903.36,4083.83);
   cutg[11]->SetPoint(12,9785.71,4580.55);
   cutg[11]->SetPoint(13,10752.1,5191.38);
   cutg[11]->SetPoint(14,11760.5,5829.06);
   cutg[11]->SetPoint(15,12579.8,6366.05);
   cutg[11]->SetPoint(16,13252.1,6835.92);
   cutg[11]->SetPoint(17,14071.4,7433.32);

   cutg[12] = new TGraph(13);
   cutg[12]->SetPoint (0,1991.6 ,566.527);
   cutg[12]->SetPoint( 1,3252.1 ,1056.53);
   cutg[12]->SetPoint( 2,4260.5 ,1479.42);
   cutg[12]->SetPoint( 3,5268.91,1922.43);
   cutg[12]->SetPoint( 4,6424.37,2459.43);
   cutg[12]->SetPoint( 5,7558.82,3009.84);
   cutg[12]->SetPoint( 6,8567.23,3506.56);
   cutg[12]->SetPoint( 7,9743.7 ,4077.12);
   cutg[12]->SetPoint( 8,10668.1,4567.12);
   cutg[12]->SetPoint( 9,11655.5,5077.27);
   cutg[12]->SetPoint(10,12516.8,5533.71);
   cutg[12]->SetPoint(11,13546.2,6097.55);
   cutg[12]->SetPoint(12,14953.8,6882.91);

   cutg[13] = new TGraph(7);
   cutg[13]->SetPoint(0,5920.17,2130.52);
   cutg[13]->SetPoint(1,6844.54,2519.84);
   cutg[13]->SetPoint(2,7852.94,2976.28);
   cutg[13]->SetPoint(3,9386.55,3687.8 );
   cutg[13]->SetPoint(4,11445.4,4614.11);
   cutg[13]->SetPoint(5,13399.2,5587.41);
   cutg[13]->SetPoint(6,15416  ,6594.27);

// TF1 *fifu= new TF1("fifu",fitfun,0,5000,4);
TF1 *fifu= new TF1("fifu",fitfun,0,16000,NFITPAR+3);
fitter->SetUserFunc(fifu);
fitter->SetFitMethod("chisquare");
fitter->SetFCN(FCN);

TCanvas *c1 = new TCanvas("c1","c1",0,0,743,525);
c1->cd();

TH1F *null1 = new TH1F("null1"," ",200,0,16000);
null1->SetMinimum(0);
null1->SetMaximum(8000);
null1->SetDirectory(0);
null1->Draw();

  Double_t arglist[100];
  arglist[0] = 0;
  fitter->ExecuteCommand("SET PRINT", arglist, 1); // no print
//  arglist[0] = 0;
//  fitter->ExecuteCommand("SET NOWAR", arglist, 1); // no warnings

//    fitter->SetParameter(0,"a0",0.6,.01,0.,300);
//    fitter->SetParameter(1,"a1",0.7,.1,0.,300);
//    fitter->SetParameter(2,"a2",1.545,.01,0.001,3000);
//    fitter->SetParameter(3,"a3",0.5,.01,0,30);
//    fitter->SetParameter(4,"a4",2.2,.01,0.,150);
//    fitter->SetParameter(5,"a5",0.842,.01,0.,5);
//    fitter->SetParameter(6,"a6",1.889,.01,0.,30);
//    fitter->SetParameter(0,"a0",3.128761,.01,0.,3000);
//    fitter->SetParameter(1,"a1",4.415996,.01,0.,30000);
//    fitter->SetParameter(2,"a2",9.880118,.01,0.001,30000);
//    fitter->SetParameter(3,"a3",0.773359,.01,0,30);
//    fitter->SetParameter(4,"a4",1.261736,.01,0.,1500);
//    fitter->SetParameter(5,"a5",1.010950,.01,0.,5);
//    fitter->SetParameter(6,"a6",2.236419,.01,0.,30);
//    fitter->SetParameter(7,"a7",0.669835,.01,-10,300);
//    fitter->SetParameter(8,"a8",0.906181,.01,0.2,300);
// 
//    fitter->SetParameter(0,"a0",0.250920,.01,0.,3000);
//    fitter->SetParameter(1,"a1",0.040218,.01,0.,30000);
//    fitter->SetParameter(2,"a2",0.099969,.01,0.001,30000);
//    fitter->SetParameter(3,"a3",0.808692,.01,0,30);
//    fitter->SetParameter(4,"a4",1.392486,.01,0.,1500);
//    fitter->SetParameter(5,"a5",1.175370,.01,0.,5);
//    fitter->SetParameter(6,"a6",4.428919,.01,0.,30);
//    fitter->SetParameter(7,"a7",0,.01,-10,300);
//    //fitter->SetParameter(8,"a8",0.200000,.01,-10,300);
// 
//z1-7
//    fitter->SetParameter(0,"a0",0.864609,.01,0.,100);
//    fitter->SetParameter(1,"a1",1.260373,.01,0.,100);
//    fitter->SetParameter(2,"a2",4.637384,.01,0.001,300);
//    fitter->SetParameter(3,"a3",7.434025,.01,0.,300);
//    fitter->SetParameter(4,"a4",0.708746,.01,0.,2);
//    fitter->SetParameter(5,"a5",0.660143,.01,0.,2);
//    fitter->SetParameter(6,"a6",0.893421,.01,0.,2);
//    fitter->SetParameter(7,"a7",2.333764,.01,0.,5);
			       
    fitter->SetParameter(0,"a0",0.88931,.01,0.,10000);
    fitter->SetParameter(1,"a1",9.09410,.01,0.,10000);
    fitter->SetParameter(2,"a2",0.80758,.01,0.001,15000);
    fitter->SetParameter(3,"a3",34.2071,.01,0.,15000);
    fitter->SetParameter(4,"a4",0.79167,.01,0.,2);
    fitter->SetParameter(5,"a5",0.74398,.01,0.,2);
    fitter->SetParameter(6,"a6",0.95612,.01,0.,2);
    fitter->SetParameter(7,"a7",2.40034,.01,0.,5);
    fitter->SetParameter(8,"a8",12.6894,.01,0.,200);
    fitter->SetParameter(9,"a9",0.5,.01,0.,5);








   fitter->SetParameter(NFITPAR,"az2",1,.01,0.,30000);
   fitter->SetParameter(NFITPAR+1,"a",1,.01,0.,30000);
   fitter->SetParameter(NFITPAR+2,"z",1,.01,0.,30000);
   fitter->FixParameter(NFITPAR);
   fitter->FixParameter(NFITPAR+1);
   fitter->FixParameter(NFITPAR+2);

// fitter->FixParameter(0);
// //fitter->FixParameter(1);
// fitter->FixParameter(2);
// fitter->FixParameter(3);
// //fitter->FixParameter(4);
// //fitter->FixParameter(5);
// fitter->FixParameter(6);
// fitter->FixParameter(7);
fitter->FixParameter(9);
//fitter->FixParameter(2);
fitter->FixParameter(8);

  if(1){
    arglist[0] =  1000;
    arglist[1] =     3;
    fitter->ExecuteCommand("SEEK", arglist, 2);
  }

  if(1){
    arglist[0] = 10000;
    arglist[1] = 0.0001;
    //printf("SIMPLEX\n");
    fitter->ExecuteCommand("SIMPLEX", arglist, 2);
  }

  if(1){
    arglist[0] = 10000;
    arglist[1] = 0.001;
    fitter->ExecuteCommand("MIGRAD", arglist, 2);
  }
    
  Double_t eplus, eminus, eparab, globcc;
  Double_t value, verr, vlow, vhigh;
  Double_t ppar[12], epar[12];
  char name[80];
  for(int i=0;i<NFITPAR;i++){
    fitter->GetParameter(i, &name[0], value, verr, vlow, vhigh);
    fitter->GetErrors(i, eplus, eminus, eparab, globcc);
    ppar[i]=value;
    epar[i]=eparab;
    printf("%d %s %f %f %f %f %f\n",i,name,value,eplus, eminus, eparab, globcc);
  }  

printf("chi2/ndf=%f %d\n",ch2,NFITPAR);
    

TF1 *fq[NSPEC];
for(int i=0;i<NSPEC;i++){
//   fq[i] = new TF1(Form("fq%d",i),Form("[0]*sqrt(x*x+2*[1]*%f*x*(1+log(1+x/([1]*[2]*%f))))",AZ2[i],AZ2[i]),0,5000);
  double xx,yy;
  cutg[i]->GetPoint(cutg[i]->GetN()-1,xx,yy);
  fq[i] = new TF1(Form("fq%d",i),fitfun,0,xx,NFITPAR+3);
  //printf("%d %2d %f\n",int(Z[i]),int(A[i]),xx);
  for(int j=0;j<NFITPAR;j++) fq[i]->SetParameter(j,ppar[j]);
  fq[i]->SetParameter(NFITPAR,AZ2[i]);
  fq[i]->SetParameter(NFITPAR+1,A[i]);
  fq[i]->SetParameter(NFITPAR+2,Z[i]);
  cutg[i]->Draw("lp");
  fq[i]->SetLineWidth(1);
  fq[i]->Draw("same");
}

// for(int i=0;i<NSPEC;i++){
//   int npt= cutg[i]->GetN();
//   for(int ipt=0;ipt<npt;ipt++){
//     Double_t xx,yy;
//     cutg[i]->GetPoint(ipt,xx,yy);
//     printf("{%f,%f,%f},\n",xx,double(AZ2[i]),yy);
//   }
// }

// for(int i=0;i<NSPEC;i++){
//   int npt= cutg[i]->GetN();
//   for(int ipt=0;ipt<npt;ipt++){
//     Double_t xx,yy;
//     cutg[i]->GetPoint(ipt,xx,yy);
//     printf("{%f,%f,%f,%f},\n",xx,double(A[i]),double(Z[i]),yy);
//   }
// }


}
///////////////////////////////////////////////////////////////////////////////////////////////////////////

// double fitfun(Double_t *xx, Double_t *par){
// 
// double x = xx[0];
// double az2 = par[3];
// double fu=1000;
// double arglog = 1.+x/(par[1]*par[2]*az2);
// if(arglog<=0) return fu;
// double argsqrt = x*x+2.*par[1]*az2*x*(1.+log(arglog));
// if(argsqrt<0) return fu;
// fu=par[0]*sqrt(argsqrt);
// return fu;
// }

///////////////////////////////////////////////////////////////////////////////////////////////////////////

// double fitfun(Double_t *xx, Double_t *par){
// 
// double az2 = par[NFITPAR];
// double a   = par[NFITPAR+1];
// double z2  = pow(par[NFITPAR+2],par[4]);
// double S   = pow(xx[0]/a,par[3]);
// double F=1000;
// double arglog = 1.+S/(par[1]*par[2]*z2);
// if(arglog<=0) return F;
// double argsqrt = S*S+2.*par[1]*z2*S*(1.+log(arglog));
// if(argsqrt<0) return F;
// F=par[0]*pow(a,par[3])*sqrt(argsqrt);
// return F;
// }

///////////////////////////////////////////////////////////////////////////////////////////////////////////

// double fitfun(Double_t *xx, Double_t *par){
// 
// double az2 = par[NFITPAR];
// double a   = par[NFITPAR+1];
// double z2  = pow(par[NFITPAR+2],2);
// double S   = xx[0];
// double F=1000;
// double arglog = 1.+S/(par[1]*par[2]*a*z2);
// if(arglog<=0) return F;
// double argsqrt = pow(S,par[3])+par[1]*z2*a*S*(par[4]+log(arglog));
// if(argsqrt<0) return F;
// F=par[0]*sqrt(argsqrt);
// return F;
// }

// ///////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
// double fitfun(Double_t *xx, Double_t *par){
// 
// double a0  = par[0];
// double a1  = par[1];
// double a2  = par[2];
// double a3  = par[3];
// double a4  = par[4];
// double a5  = par[5];
// double a6  = par[6];
// double a7  = par[7];
// double a   = pow(par[NFITPAR+1],a6);
// double z2  = pow(par[NFITPAR+2],a7);
// double F   = pow(xx[0],a4);
// //double az2 = pow(par[NFITPAR],par[5]);
// // double arg = (F-par[7]*pow(F,0.64)-par[1]*a*z2*log(1.+F/(par[2]*a*z2)));
// //double arg = (F-par[7]*pow(F,par[8])-par[1]*a*z2*log(1.+F/(par[2]*a*z2)));
// //double arg = (F-par[7]*a-par[1]*a*z2*log(1.+F/(par[2]*a*z2)));
// 
// double arg = (F-a1*a-a2*a*z2*log(1.+F/(a3*a*z2)));
// double S=0;
// if(arg>0) S= a0*pow(arg,a5);
// else      S=-a0*pow(-arg,a5);
// // printf(">>  %f %f %f %g\n",par[NFITPAR+1],par[NFITPAR+2],F ,
// // S
// // );
// 
// return S;
// }
// 
// ///////////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////

double fitfun(Double_t *xx, Double_t *par){

double a0  = par[0];
double a1  = par[1];
double a2  = par[2];
double a3  = par[3];
double a4  = par[4];
double a5  = par[5];
double a6  = par[6];
double a7  = par[7];
double a8  = par[8];
double a9  = par[9];
double a   = pow(par[NFITPAR+1],a6);
double F   = pow(xx[0],a4);
double z2  = pow(par[NFITPAR+2],a7);
//double az2 = pow(par[NFITPAR],par[5]);
// double arg = (F-par[7]*pow(F,0.64)-par[1]*a*z2*log(1.+F/(par[2]*a*z2)));
//double arg = (F-par[7]*pow(F,par[8])-par[1]*a*z2*log(1.+F/(par[2]*a*z2)));
//double arg = (F-par[7]*a-par[1]*a*z2*log(1.+F/(par[2]*a*z2)));

// double arg = (F-a1*a-a2*a*z2*log(1.+F/(a3*a*z2)));
// double arg = (F-a1*a*z2*log(1.+F/(a8*a*z2))+0.5*a1*a*z2*log((F+a8*a*z2)/(a3*a+a8*a*z2)));
double arg = (F-a1*a*z2*log(1.+F/(a1*a*z2))+a2*a1*a*z2*log((F+a1*a*z2)/(a3*a+a1*a*z2)));

double S=0;
if(arg>0) S= a0*pow(arg,1./a5);
else      S=-a0*pow(-arg,1./a5);
// printf(">>  %f %f %f %g\n",par[NFITPAR+1],par[NFITPAR+2],F ,
// S
// );

return S;
}

///////////////////////////////////////////////////////////////////////////////////////////////////////////

void FCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag){

Double_t chi2=0.;

   Double_t cu,eu,fu,fsum;
   Double_t *zik=0;
   Double_t *pl0=0;
   Int_t i;

#ifdef FUMILI
TFumili *fitter = (TFumili*)TVirtualFitter::GetFitter();
#else
TVirtualFitter *fitter = TVirtualFitter::GetFitter();
#endif

//printf("(1) %f %f %f %f\n",par[0],par[1],par[2],par[3]);

TF1 *f1 = (TF1*)fitter->GetUserFunc();
Double_t x[3]={0.,0.,0.};
npar = f1->GetNpar();

#ifdef FUMILI
fitter->SetParNumber(npar);
#endif
if(iflag == 9) return;

#ifdef FUMILI
zik = fitter->GetZ();
pl0 = fitter->GetPL0();

Double_t *df=new Double_t[npar];
#endif
f1->InitArgs(x,par);

f = 0;

int np=0;
//for(int icut=0;icut<7;icut++){
//for(int icut=3;icut<7;icut++){
//for(int icut=3;icut<10;icut++){
//for(int icut=10;icut<13;icut++){
//for(int icut=0;icut<NSPEC;icut++){

//---------
//for(int icut=3;icut<NSPEC;icut++){
//  for(int icut=0;icut<7;icut++){
//---------
 
 for(int icut=0;icut<NSPEC;icut++){
//  for(int icut=0;icut<3;icut++){
//  for(int icut=3;icut<7;icut++){
//for(int icut=7;icut<10;icut++){
// for(int icut=10;icut<NSPEC;icut++){
// for(int icut=12;icut<=12;icut++){

 int npt= cutg[icut]->GetN();
  par[NFITPAR]   = AZ2[icut];
  par[NFITPAR+1] = A[icut];
  par[NFITPAR+2] = Z[icut];
  //printf(">> %d %f %f %f\n",icut,par[NFITPAR],par[NFITPAR+1],par[NFITPAR+2]);
  for(int ipt=0;ipt<npt;ipt++){
    Double_t xx,yy;
    cutg[icut]->GetPoint(ipt,xx,yy);
    x[0]=xx;
    cu=yy;
//    eu=(Double_t)hist->GetBinError(ix);
//    eu=1;
  eu=sqrt(fabs(cu));
//   if(cu>0 && eu>0){
  //printf(">> %d %f %f %f\n",ipt,xx,yy,AZ2[icut]);
  if(eu>0){
//  printf(">>> %d %f %f %f\n",ipt,xx,yy,par[0]*sqrt(xx*xx+2.*par[1]*AZ2[icut]*xx*(1.+log(1.+xx/(par[1]*par[2]*AZ2[icut])))));
//     fu=par[0]*sqrt(xx*xx+2.*par[1]*AZ2[icut]*xx*(1.+log(1.+xx/(par[1]*par[2]*AZ2[icut]))));
    fu=f1->EvalPar(x,par);
  //printf(">>> %d %f %f %f\n",ipt,xx,yy,fu);
//     if(fu==-1.){
//       f=1.e9;
//       return;
//     }
    np++;
    fsum = (fu-cu)/eu;
#ifdef FUMILI
    fitter->Derivatives(df,x);
  //printf(">>> %d %f %f %f\n",ipt,xx,yy,fu);
    Int_t N = 0;
    if(iflag!=1){
      for(i=0;i<npar;i++) 
  	if(pl0[i]>0){
  	  df[N] = df[i]/eu; 
  	  // left only non-fixed param derivatives / by Sigma
  	  gin[i] += df[N]*fsum;
  	  N++;
  	}
      Int_t L = 0;
      for (i=0;i<N;i++)
  	for (Int_t j=0;j<=i;j++) 
  	  zik[L++] += df[i]*df[j];
    }
#endif
//     f += .5*fsum*fsum;
    f += fsum*fsum;
 }
 }
}

#ifdef FUMILI
delete[] df;
#endif

chi2=f;

// take care of the fixed parameters...
//Int_t  ndof=np-1-(NFITPAR-1);
Int_t  ndof=np-1;
//printf("chi2/ndf=%f %d\n",chi2/ndof,NFITPAR);
ch2 = chi2/ndof;
//ch2 = chi2;
if(iflag==3){
  //printf("%c\n",7);
  //printf("%c\n",7);

//   double pr_chi2 = 1.-TMath::Gamma((Double_t)(0.5*ndof),(Double_t)(0.5*chi2));  
//   double ch2 = (ndof)? chi2/ndof : chi2;
//   pr_chi2 = 1.-TMath::Gamma((Double_t)(0.5*ndof),(Double_t)(0.5*chi2));  
//   ch2 = (ndof>0)? chi2/ndof : chi2;
  //printf("chi2=%f chi2/ndf=%f pr_chi2=%f\n",chi2,ch2,pr_chi2);

//   printf("%7.1f %5.3f %5.3f %7.1f %e\n",par[0],par[1],par[2],par[3],par[4]);
//   printf("%7.1f %5.3f %5.3f %7.1f %e\n",
//   fitter->fParamError[0],
//   fitter->fParamError[1],
//   fitter->fParamError[2],
//   fitter->fParamError[3],
//   fitter->fParamError[4]
//   );
  
// printf("%8.3f %8.3f  %7.3f %7.3f  %7.3f %7.3f  %8.3f %8.3f  %7.3f\n",
// par[0],fitter->fParamError[0],
// par[1],fitter->fParamError[1],
// par[2],fitter->fParamError[2],
// par[3],fitter->fParamError[3],ch2);
  
//  fitter->PrintResults(4,f);

//   for(int i=0;i<5;i++){
//     fitter->GetParameter(i, &name[0], value, verr, vlow, vhigh);
//     fitter->GetErrors(i, eplus, eminus, eparab, globcc);
//     ppar[i]=value;
//     epar[i]=eparab;
//     printf("%d %s %f %f %f %f %f\n",i,name,value,eplus, eminus, eparab, globcc);
//   } 
}
//printf("chi2/ndf=%f\n",chi2/ndof);

}