//
// plot unpolarized PDF distributions (GRV98, CTEQ) with various assumptions
//
int Plot_GRV98(char * InFile1 = "GRV98_Q10.txt"){

  
  gROOT->Reset();
  
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat("Plain");
  //  gStyle->SetPadLeftMargin(0.15);
  //  gStyle->SetPadBottomMargin(0.2);
  //  gStyle->SetLineWidth(1.5);

  m1 = new TCanvas("m1","GRV98 NLO/MS PDF @Q^2 = 16GeV^2", 10, 20,600,500);


  FILE *F_PDF1   = fopen(InFile1,"r");   // unpolarized_PDF 



  //input variables

  const Int_t n = 800;
  Int_t ncol1 =0;
  
  Float_t  x[n], uv[n],Puv[n],dv[n],Pdv[n],us[n],Pus[n],ds[n],Pds[n],gx[n],dgx_std[n],dgx_max[n];
  Float_t        us[n],Pus[n],ds[n],Pds[n],ss[n],Pss[n];
  Float_t        xgx[n], xdgx_std[n],xdgx_max[n];
  Float_t        xuv[n], Pxuv[n],xdv[n],Pxdv[n], xus[n],Pxus[n],xds[n],Pxds[n], xss[n],Pxss[n]  ;
  Float_t  x[n], ddgx_std[n], ddgx_max[n];
  Float_t        duv[n], ddv[n], dus[n],dds[n], dss[n];
 

  
  // ------- POL PDF STD  -----
  cout << "GRV98 PDF, NLO/MS" << endl;
  for (int i=0; i<n;i++){
      
    ncol1 = fscanf(F_PDF1, "%f %f %f %f %f %f %f", &x[i], &uv[i],&dv[i],&gx[i],&us[i],&ds[i],&ss[i]);

    //valence quarks
    xuv[i] = x[i]*uv[i];
     
    xdv[i] = x[i]*dv[i];

    //sea quarks
    xus[i] = x[i]*us[i];

    xds[i] = x[i]*ds[i];

    ss[i] = us[i];
    xss[i] = x[i]*ss[i];


    //gluon
    xgx[i] = x[i]*gx[i];    
    xdgx_std[i] = x[i]*dgx_std[i];
    xdgx_max[i] = x[i]*dgx_max[i];




    if (ncol1<0) break;

    cout << "x = " << x[i] << "\t uv = " << uv[i] << "\t dv= " << dv[i] << "\t g(x)=" << gx[i] << endl; 
  } // read in all data points    


  // ----- make plots ----
  
  // --- Pol-PDF--- x*dF(x)---
  m1.cd();
  	
  m1.SetLogy();
  m1.SetLogx();

  //--- GRV98-PDF-STD
  PDF_STD = new TGraph(n, x, xgx);
  PDF_STD->SetTitle("x*f(x) vs x");
  PDF_STD->GetXaxis()->SetTitle("x");
  PDF_STD->GetYaxis()->SetTitle("x*f(x)");
  
  PDF_STD->GetXaxis()->CenterTitle();
  PDF_STD->GetYaxis()->CenterTitle();
  
  PDF_STD->GetXaxis()->SetLimits(0.001,0.8);
  PDF_STD->SetLineColor(1);
  PDF_STD->SetLineWidth(2);
  PDF_STD->SetLineStyle(2);
  
//  PDF_STD->SetMaximum(2.00);
//  PDF_STD->SetMinimum(-0.5);
  PDF_STD->Draw("AC");
  


  PDF_UV = new TGraph(n, x, xuv);
  PDF_UV->SetLineColor(2);
  PDF_UV->SetLineWidth(2);
  PDF_UV->SetLineStyle(2);

  PDF_DV = new TGraph(n, x, xdv);
  PDF_DV->SetLineColor(8);
  PDF_DV->SetLineWidth(2);
  PDF_DV->SetLineStyle(2);

  PDF_US = new TGraph(n, x, xus);
  PDF_US->SetLineColor(4);
  PDF_US->SetLineWidth(2);
  PDF_US->SetLineStyle(2);



  PDF_UV->Draw();
  PDF_DV->Draw();
  PDF_US->Draw();



  
  //draw text
  TLatex lx;
  
  lx.SetTextSize(0.035);

  lx.SetTextColor(1);
  lx.DrawLatex(0.002,3.7,"Q = 10.0GeV");

  lx.SetTextColor(1);
  lx.DrawLatex(0.002,1.5,"GRV98A: xg(x)");



// fit functions

// TF1 *pdf_xgx = new TF1("pdf_xgx", "[0]+[1]*log(x) + [2]*(log(x))**2 + [3]*(log(x))**3");
// TF1 *pdf_xgx = new TF1("pdf_xgx", "[0]+[1]*log(x) + [2]*(log(x))**2");

 TF1 *pdf_xgx = new TF1("pdf_xgx", "[0]+[1]*log10(x) + [2]*(log10(x))**2 + [3]*(log10(x))**3");

//PDF_STD->Fit("pol8");

PDF_STD->Fit("pdf_xgx");
}