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

//_____________________________________________________________________________
TGraph*  cteq6_x_53_pdf(Int_t ISet, Int_t IPart, Double_t Q) {
  //  distribution for x^(5/3)*G(x,Q)
  int const np = 10000;
  double xx, x[np], y[np];

  TCteq6* cteq6 = TCteq6::Instance();

  cteq6->Setctq6(ISet);

  for (int i=0; i<np; i++) {
    xx  = (i+0.5)/np;
    x[i] = xx; // pow(xx,1./3.);
    //    y[i] = pow(xx,5./3.)*cteq6->Ctq6pdf(IPart,xx,Q);
    y[i] = xx*cteq6->Ctq6pdf(IPart,xx,Q);

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

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

  return gr;
}

//_____________________________________________________________________________
void test1() {
  // this test reproduces Fig 2. from hep-ph/0303013 v1

//    TH2F* h = new TH2F("h","Fig.2 from hep-ph/0303013: x^{5/3}*G(x,Q), Q^{2} = 10 GeV",
//  		     10000,0,1.,1,0,0.4);

  TH2F* h = new TH2F("h","CTEQ6M PDF: x*G(x,Q), Q^{2} = 10 GeV, arbitrary normalization",
		     10000,0,1.,1,0,5.);

  h->GetXaxis()->SetTitle("X: fraction of the proton momentum");

  TGraph* gr_200_g  = cteq6_x_53_pdf(200, 0,3.16);
  TGraph* gr_200_u  = cteq6_x_53_pdf(200, 1,3.16);
  TGraph* gr_200_d  = cteq6_x_53_pdf(200, 2,3.16);
  TGraph* gr_200_ub = cteq6_x_53_pdf(200,-1,3.16);
  TGraph* gr_200_db = cteq6_x_53_pdf(200,-2,3.16);
  TGraph* gr_200_s  = cteq6_x_53_pdf(200, 3,3.16);

//    TGraph* gr_229 = cteq6_x_53_pdf(229,0,3.16);
//    TGraph* gr_230 = cteq6_x_53_pdf(230,0,3.16);

  h->SetStats(0);
  h->Draw();

  

  gr_200_g->Draw();

  gr_200_u->SetLineColor(8);
  gr_200_u->Draw();

  gr_200_d->SetLineColor(6);
  gr_200_d->Draw();
  //  gr_200_ub->Draw();

  gr_200_db->SetLineColor(4);
  gr_200_db->Draw();

  gr_200_s->SetLineColor(2);
  gr_200_s->Draw();

//    gr_229->Draw();
//    gr_230->Draw();
}