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

//_____________________________________________________________________________
TGraph*  pdflib_x_53_pdf(Int_t ISet, Int_t IPart, Double_t Q)
{
  // distribution for x^(5/3)*G(x,Q), values of IPart:
  // 0: gluon, 1-6:quarks, -1:-6: antiquarks
  // Structm takes Q, not Q2

  // part[0][0] : g , [0][1]: not used
  // part[1][0] : d , [1][1]: dbar
  // part[2][0] : u , [2][1]: ubar
  // part[3][0] : s , [3][1]: sbar
  // part[4][0] : c , [4][1]: cbar
  // part[5][0] : b , [5][1]: bbar
  // part[6][0] : t , [6][1]: tbar

  int ind1, ind2;

  const int np = 10000;
  double xx, x[np], y[np];
  double uval, dval, usea, dsea, str, chm, bot, top, glu;
  double ngroup, nset;

  double part[7][2];

  TG3Pythia6* py = TG3Pythia6::Instance();

  ngroup = ISet / 1000;
  nset   = ISet % 1000;

  printf (" ngroup, nset: %3i %3i\n",ngroup, nset);

  py->Pdfset("NPTYPE",1);
  py->Pdfset("NGROUP",ngroup);
  py->Pdfset("NSET"  ,nset);

  if (IPart >= 0) {
				// quark / gluon
    ind1 = IPart;
    ind2 = 0;
  }
  else {
				// antiquark
    ind1 = -IPart;
    ind2 = 1;
  }

  for (int i=0; i<np; i++) {
    xx  = (i+0.5)/np;
    py->Structm(xx,Q,uval,dval,usea,dsea,str,chm,bot,top,glu);

    part[0][0] = glu;
    part[0][1] = 0;
    part[1][0] = uval+usea;
    part[1][1] = usea;
    part[2][0] = dval+dsea;
    part[2][1] = dsea;
    part[3][0] = str;
    part[3][1] = str;
    part[4][0] = chm;
    part[4][1] = chm;
    part[5][0] = bot;
    part[5][1] = bot;
    part[6][0] = top;
    part[6][1] = top;

    double p = part[ind1][ind2];

    x[i] = xx; // pow(xx,1./3.);
    y[i] = pow(xx,5./3.-1)*p;

    //    printf("x,y = %10.5f  %10.5f \n",x[i],y[i]);
  }

  TGraph* gr = new TGraph(np,x,y);

  return gr;
}

//_____________________________________________________________________________
void pdflib_test1() {

  TH2F* h = new TH2F("h_4046","PDFLIB: CTEQ5L(4046)",
		     10000,0,1.,1,0,0.4);

  float x[10001];

//    for (int i=0; i<10001; i++) {
//      float xx = float(i)/10000.;
//      x[i] = pow(xx,1./3.);
//    }

//    TAxis* xax = h->GetXaxis();

//    xax->SetTitle("X^{1/3}");
//    //  xax->Set(10000,x);

  h->Draw();

  TGraph* glu_4046 = pdflib_x_53_pdf(4046,0,3.16);

  xax = glu_4046->GetXaxis();

//    xax->SetTitle("X^{1/3}");
//    xax->Set(10000,x);


  glu_4046->Draw();

//    TGraph* d_4046   = pdflib_x_53_pdf(4046,1,3.16);
//    d_4046->Draw();

  TGraph* u_4046   = pdflib_x_53_pdf(4046,2,3.16);
  u_4046->Draw();

  TGraph* d_4046   = pdflib_x_53_pdf(4046,1,3.16);
  d_4046->Draw();
}