//const int nR1 = 42;
const int nR1 = 35;


Double_t cteq6(Double_t x) {
  Double_t p[] = 
  {1.22024e+13, -1.74694e+01, 1.82338e+01, -1.22027e+01, 3.86174e+00, -4.74714e-01, 5.10307e+00};

  Double_t val;
  val = p[0]*TMath::Power(10,(p[1]*TMath::log10(x)+
			      p[2]*TMath::Power(TMath::log10(x),2)+p[3]*TMath::Power(TMath::log10(x),3)+
			      p[4]*TMath::Power(TMath::log10(x),4)+p[5]*TMath::Power(TMath::log10(x),5)))*
    TMath::Power((1-x/900),p[6]);
  return val;
}


void plotPar () {

   Double_t CTEQ6_Et[  42]={  
     20.0,  25.0,  31.6,  37.6,  43.7,  49.6,  55.4,  61.2,  66.9,  72.5,  78.1,  83.7,  89.2,  94.7, 100.2, 105.7, 111.1, 116.5, 122.0, 127.4, 132.8, 138.2, 146.0, 156.8, 167.6, 178.4, 189.2, 200.0, 210.8, 226.1, 247.9, 269.7, 291.5, 313.2, 334.8, 363.8, 415.4, 490.7, 552.9, 600.0, 650.0, 700.0};

   Double_t CTEQ6_EtErr[  42]={  
     20.0,  25.0,  31.6,  37.6,  43.7,  49.6,  55.4,  61.2,  66.9,  72.5,  78.1,  83.7,  89.2,  94.7, 100.2, 105.7, 111.1, 116.5, 122.0, 127.4, 132.8, 138.2, 146.0, 156.8, 167.6, 178.4, 189.2, 200.0, 210.8, 226.1, 247.9, 269.7, 291.5, 313.2, 334.8, 363.8, 415.4, 490.7, 552.9, 600.0, 650.0, 700.0};

   Double_t CTEQ6_sigma[  42]={  
     3.989e+03, 1.286e+03, 3.823e+02, 1.555e+02, 7.019e+01, 3.516e+01, 1.900e+01, 1.087e+01, 6.582e+00, 4.161e+00, 2.711e+00, 1.810e+00, 1.244e+00, 8.721e-01, 6.207e-01, 4.486e-01, 3.308e-01, 2.470e-01, 1.854e-01, 1.413e-01, 1.087e-01, 8.422e-02, 5.901e-02, 3.689e-02, 2.361e-02, 1.545e-02, 1.028e-02, 6.943e-03, 4.747e-03, 2.820e-03, 1.384e-03, 7.001e-04, 3.625e-04, 1.916e-04, 1.028e-04, 4.505e-05, 1.051e-05, 1.229e-06, 1.946e-07, 4.524e-08, 8.965e-09, 1.640e-09};
   

   Double_t xCorEta1[] = {
     43.350, 49.399, 55.305, 61.101, 66.817, 72.471, 78.078, 83.646, 89.182, 94.691, 100.177, 105.644, 111.095, 116.531, 121.956, 127.371, 132.779, 138.180, 145.953, 156.751, 167.542, 178.335, 189.135, 199.945, 210.769, 226.078, 247.849, 269.645, 291.429, 313.145, 334.723, 363.625, 415.051, 490.243, 552.487};

   Double_t eylCorEta1[] = {     
     0.14559, 0.10206, 0.0746, 0.056725, 0.044515, 0.034989, 0.0059305, 0.0048382,  0.0040188, 0.0033665, 0.0028459, 0.0014336, 0.0012259, 0.0010613, 0.00091571, 0.00080592, 0.0006999, 0.00061667, 0.00012953, 0.00010255, 8.2076e-05, 6.65e-05, 5.3819e-05, 4.4363e-05, 3.6862e-05, 2.0497e-05, 1.4403e-05, 1.05e-05, 7.6994e-06, 5.9408e-06, 4.4442e-06, 2.1201e-06, 8.1898e-07, 3.2548e-07, 1.2951e-07};



   Double_t xR1[nR1];
   Double_t yR1[nR1];
   Double_t yRh1[nR1];
   Double_t yRl1[nR1];
   Double_t exlR1[nR1];
   Double_t exhR1[nR1];
   Double_t eylR1[nR1];
   Double_t eyhR1[nR1];

   Double_t yR1Err[nR1];

   for (int i=0; i<nR1; i++) {

     Double_t *x,*y,xx,z;
     x = &xx;

     //     xx        = CTEQ6_Et[i];
     //     xR1[i]    = CTEQ6_Et[i];
     //     yR1[i]    = CTEQ6_sigma[i];

     xx        = xCorEta1[i];
     xR1[i]    = xCorEta1[i];
     //     yR1[i]    = cteq61(x, y);
     yR1[i]    = cteq6(xx);

     exlR1[i]  = 0.;
     exhR1[i]  = 0.;
     eylR1[i]  = eylCorEta1[i];
     eyhR1[i]  = eylCorEta1[i];

     yR1Err[i] = 0;

     // Underlying event
     //     yRh1[i] = yR1[i] + sysUEDiffP(x, y)*yR1[i];
     // el/ph
     //     yRh1[i] = yR1[i] + syselphDiffP(x, y)*yR1[i];
     // high pt
     //     yRh1[i] = yR1[i] + syshiptDiffP(x, y)*yR1[i];
     // low pt
     //     yRh1[i] = yR1[i] + sysloptDiffP(x, y)*yR1[i];
     // fragmentation
     //     yRh1[i] = yR1[i] + sysFragDiffP(x, y)*yR1[i];
     // resolution
     //     yRh1[i] = yR1[i] + sysResDiffP(x, y)*yR1[i];
     // energy scale
     //     yRh1[i] = yR1[i] + sysEScaleDiffP(x, y)*yR1[i];
     yRh1[i] = yR1[i] + sysTotDiffP(x, y)*yR1[i];



     // Underlying event
     //     yRl1[i] = yR1[i] + sysUEDiffN(x, y)*yR1[i];
     // el/ph
     //     yRl1[i] = yR1[i] + syselphDiffN(x, y)*yR1[i];
     // high pt
     //     yRl1[i] = yR1[i] + syshiptDiffN(x, y)*yR1[i];
     // low pt
     //     yRl1[i] = yR1[i] + sysloptDiffN(x, y)*yR1[i];
     // fragmentation
     //     yRl1[i] = yR1[i] + sysFragDiffN(x, y)*yR1[i];
     // resolution
     //     yRl1[i] = yR1[i] + sysResDiffN(x, y)*yR1[i];
     // energy scale
     //     yRl1[i] = yR1[i] + sysEScaleDiffN(x, y)*yR1[i];
     yRl1[i] = yR1[i] + sysTotDiffN(x, y)*yR1[i];

     cout << xR1[i] << " " << yR1[i] << " " 
	  << eylR1[i] << endl;

     cout << " "
	  << sysUEDiffP(x, y)     << " " << sysUEDiffN(x, y)   << " "
	  << syselphDiffP(x, y)   << " " << syselphDiffN(x, y) << " "
	  << syshiptDiffP(x, y)   << " " << syshiptDiffN(x, y) << " "
	  << sysloptDiffP(x, y)   << " " << sysloptDiffN(x, y) << " "
	  << sysFragDiffP(x, y)   << " " << sysFragDiffN(x, y) << " "
	  << sysResDiffP(x, y)    << " " << sysResDiffN(x, y)  << " "
	  << sysEScaleDiffP(x, y) << " " << sysEScaleDiffN(x, y) 
	  << endl;  

   }

   TGraphAsymmErrors *grR1  = new TGraphAsymmErrors(nR1, xR1, yR1, exlR1, 
						    exhR1, eylR1, eyhR1);

   //   TGraph *grR1      = new TGraph(nR1, xR1, yR1);
   TGraph *grRh1     = new TGraph(nR1, xR1, yRh1);
   TGraph *grRl1     = new TGraph(nR1, xR1, yRl1);

   //***************
   //***************
   can_Par1 = new TCanvas("par1");
   can_Par1->Divide(1,1);

   can_Par1->cd(1); 

   gPad->SetLogy();

   gPad->SetBottomMargin(0.18);
   gPad->SetTopMargin(0.05);
   gPad->SetLeftMargin(0.12);
   gPad->SetRightMargin(0.05);

   grR1->SetLineColor(kRed);
   grR1->SetLineWidth(3);
   grR1->SetMarkerColor(kRed);
   grR1->SetMarkerStyle(20);

   grR1->SetMaximum(1000);
   grR1->SetMinimum(.00000001);
   grR1->Draw("AL");
   grRh1->Draw("LS");
   grRl1->Draw("LS");
  
}