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

//_____________________________________________________________________________
TGraph*  mrst_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();

  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->Mrst2003c(xx,Q,ISet,uval,dval,usea,dsea,str,chm,bot,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] = 0;
    part[6][1] = 0;

    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 mrst2003_test(Int_t IPart) {
  // IPart=1: U
  // IPart=2: D
  // IPart=6: G

  TGraph* g_mrst;

  TH2F* h_mrst2003 = new TH2F("h_mrst2003","MRST2003, Q^2 = 10 GeV", 
			      10000,0,1.,1,0,0.4);

  h_mrst2003->Draw();

  g_mrst = mrst_x_53_pdf(1,IPart,3.16);
  g_mrst->Draw();
}

//_____________________________________________________________________________
void mrst2003_test1() {
  // IPart=1: U
  // IPart=2: D
  // IPart=6: G

  TGraph *g, *u, *d;

  TH2F* h = new TH2F("h_mrst2003","MRST2003, Q^2 = 10 GeV", 
			      10000,0,1.,1,0,0.4);

  //  h->GetXaxis()->SetTitle("X^{1/3}");

  h->Draw();

  g = mrst_x_53_pdf(1,0,3.16);
  g->Draw();

  u = mrst_x_53_pdf(1,1,3.16);
  u->Draw();

  d = mrst_x_53_pdf(1,2,3.16);
  d->Draw();

}