void pr_pb_yspectra()
{
  //
  // (Anti-)proton rapidity distributions
  //

  gROOT->Reset();

  gStyle->SetOptDate(0);
  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  gStyle->SetPadBorderSize(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetPadLeftMargin(0.18);
  gStyle->SetPadRightMargin(0.08);
  gStyle->SetPadTopMargin(0.05);
  gStyle->SetPadBottomMargin(0.12);
  gStyle->SetLabelSize(0.045,"X");
  gStyle->SetLabelSize(0.045,"Y");
  gStyle->SetTitleOffset(1.1,"X");
  gStyle->SetTitleOffset(1.6,"Y");
  gStyle->SetTitleSize(0.05,"X");
  gStyle->SetTitleSize(0.05,"Y");

  gStyle->SetPadTickX(0);
  gStyle->SetPadTickY(0);

  Double_t marker_size = 2.5;

  Double_t kYbeam20    = 1.88;
  Double_t kYbeam30    = 2.08;
  Double_t kYbeam40    = 2.22;
  Double_t kYbeam80    = 2.56;

  Double_t scalePr20   = 3.0;
  Double_t scalePr30   = 2.0;
  Double_t scalePr40   = 1.5;
  Double_t scalePr80   = 1.0;

  Double_t scalePb40   = 1.0;
  Double_t scalePb80   = 1.0;

  //____________________________________________________________________________
  //

  //
  // 20AGeV
  //

  // Protons
  const Int_t nPr20  = 7;
  Double_t yPr20[nPr20]     = {  1.88,     2.3,       2.5,       2.7
                              ,  2.9,      3.1,       3.3      };
  Double_t eyPr20[nPr20]    = {  0.0,      0.0,       0.0,       0.0
                              ,  0.0,      0.0,       0.0      };
  Double_t dndyPr20[nPr20]  = { 46.1,     47.217,    44.2534,   41.519
                              , 37.7272,  32.257,    26.0018   };
  Double_t edndyPr20[nPr20] = {  2.1,      0.724,     0.628652,  0.550376
                              ,  0.496394, 0.444072,  0.378252 };
  Double_t ryPr20[nPr20];
  Double_t rdndyPr20[nPr20];
  Double_t redndyPr20[nPr20];
  for (Int_t i = 0; i < nPr20; i++) {
    yPr20[i]      = yPr20[i] - kYbeam20;
    dndyPr20[i]  *= scalePr20;
    edndyPr20[i] *= scalePr20;
  }
  for (Int_t i = 0; i < nPr20; i++) {
    ryPr20[i]     = 0.0 - yPr20[nPr20-1-i];
    rdndyPr20[i]  =  dndyPr20[nPr20-1-i];
    redndyPr20[i] = edndyPr20[nPr20-1-i];
  }

  TGraphErrors *grPr20 = new TGraphErrors(nPr20,ryPr20, rdndyPr20
                                               ,eyPr20,redndyPr20);
  TGraphErrors *gPr20  = new TGraphErrors(nPr20, yPr20,  dndyPr20
                                               ,eyPr20, edndyPr20);

  TF1 *fPr20 = new TF1("fPr20",yFun1,-2.5,2.5,3);
  fPr20->SetParameter(0,1.0);
  fPr20->SetParameter(1,0.5);
  fPr20->SetParameter(2,0.5);

  gPr20->Fit("fPr20","0");
  fPr20->SetRange(-2.5,2.5);

  Double_t yldPr20 = fPr20->Integral(-kYbeam20,kYbeam20) / scalePr20;
  Double_t extPr20 = scalePr20 * yldPr20 / fPr20->Integral(0.0,yPr20[nPr20-1]);
  Double_t errPr20 = 0.0;
  for (Int_t i = 0; i < nPr20; i++) {
    errPr20 += edndyPr20[i]*edndyPr20[i];
  }
  errPr20 = TMath::Sqrt(errPr20) * extPr20 / scalePr20;

  //____________________________________________________________________________
  //

  //
  // 30AGeV
  //

  // Protons
  const Int_t nPr30  = 7;
  Double_t yPr30[nPr30]      = {  2.08,     2.5,       2.7,       2.9
                               ,  3.1,      3.3,       3.5      };
  Double_t eyPr30[nPr30]     = {  0.0,      0.0,       0.0,       0.0
                               ,  0.0,      0.0,       0.0      };
  Double_t dndyPr30[nPr30]   = { 42.1,     40.5007,   41.0757,   39.8223
                               , 37.1432,  33.9603,   29.3874   };
  Double_t edndyPr30[nPr30]  = {  2.0,      0.53978,   0.510842,  0.482998
                               ,  0.457904, 0.426083,  0.403031 };
  Double_t ryPr30[nPr30];
  Double_t rdndyPr30[nPr30];
  Double_t redndyPr30[nPr30];
  for (Int_t i = 0; i < nPr30; i++) {
    yPr30[i]      = yPr30[i] - kYbeam30;
    dndyPr30[i]  *= scalePr30;
    edndyPr30[i] *= scalePr30;
  }
  for (Int_t i = 0; i < nPr30; i++) {
    ryPr30[i]     = 0.0 - yPr30[nPr30-1-i];
    rdndyPr30[i]  =  dndyPr30[nPr30-1-i];
    redndyPr30[i] = edndyPr30[nPr30-1-i];
  }

  TGraphErrors *grPr30 = new TGraphErrors(nPr30,ryPr30, rdndyPr30
                                               ,eyPr30,redndyPr30);
  TGraphErrors *gPr30  = new TGraphErrors(nPr30, yPr30,  dndyPr30
                                               ,eyPr30, edndyPr30);

  TF1 *fPr30 = new TF1("fPr30",yFun1,-2.5,2.5,3);
  fPr30->SetParameter(0,1.0);
  fPr30->SetParameter(1,0.5);
  fPr30->SetParameter(2,0.5);

  gPr30->Fit("fPr30","0");
  fPr30->SetRange(-2.5,2.5);

  Double_t yldPr30 = fPr30->Integral(-kYbeam30,kYbeam30) / scalePr30;
  Double_t extPr30 = scalePr30 * yldPr30 / fPr30->Integral(0.0,yPr30[nPr30-1]);
  Double_t errPr30 = 0.0;
  for (Int_t i = 0; i < nPr30; i++) {
    errPr30 += edndyPr30[i]*edndyPr30[i];
  }
  errPr30 = TMath::Sqrt(errPr30) * extPr30 / scalePr30;

  //____________________________________________________________________________
  //

  //
  // 40AGeV
  //

  // Protons
  const Int_t nPr40  = 6;
  Double_t yPr40[nPr40]      = {  2.22,       2.7,        2.9,        3.1
                               ,  3.3,        3.5       };
  Double_t eyPr40[nPr40]     = {  0.0,        0.0,        0.0,        0.0
                               ,  0.0,        0.0       };
  Double_t dndyPr40[nPr40]   = { 41.3,       40.1196,    41.0409,    40.2612
                               , 37.7089,    34.3201    };
  Double_t edndyPr40[nPr40]  = {  1.1,        0.501328,   0.49989,    0.506497
                               ,  0.475625,   0.484725  };
  Double_t ryPr40[nPr40];
  Double_t rdndyPr40[nPr40];
  Double_t redndyPr40[nPr40];
  for (Int_t i = 0; i < nPr40; i++) {
    yPr40[i]      = yPr40[i] - kYbeam40;
    dndyPr40[i]  *= scalePr40;
    edndyPr40[i] *= scalePr40;
  }
  for (Int_t i = 0; i < nPr40; i++) {
    ryPr40[i]     = 0.0 - yPr40[nPr40-1-i];
    rdndyPr40[i]  =  dndyPr40[nPr40-1-i];
    redndyPr40[i] = edndyPr40[nPr40-1-i];
  }

  TGraphErrors *grPr40 = new TGraphErrors(nPr40,ryPr40, rdndyPr40
                                               ,eyPr40,redndyPr40);
  TGraphErrors *gPr40  = new TGraphErrors(nPr40, yPr40,  dndyPr40
                                               ,eyPr40, edndyPr40);

  TF1 *fPr40 = new TF1("fPr40",yFun1,-2.5,2.5,3);
  fPr40->SetParameter(0,1.0);
  fPr40->SetParameter(1,0.5);
  fPr40->SetParameter(2,0.2);

  gPr40->Fit("fPr40","0");
  fPr40->SetRange(-2.5,2.5);

  Double_t yldPr40 = fPr40->Integral(-kYbeam40,kYbeam40) / scalePr40;
  Double_t extPr40 = scalePr40 * yldPr40 / fPr40->Integral(0.0,yPr40[nPr40-1]);
  Double_t errPr40 = 0.0;
  for (Int_t i = 0; i < nPr40; i++) {
    errPr40 += edndyPr40[i]*edndyPr40[i];
  }
  errPr40 = TMath::Sqrt(errPr40) * extPr40 / scalePr40;

  // Anti-protons
  const Int_t nPb40  = 5;
  Double_t yPb40[nPb40]      = {  2.22,       2.7,        2.9,        3.1
                               ,  3.3       };
  Double_t eyPb40[nPb40]     = {  0.0,        0.0,        0.0,        0.0
                               ,  0.0       };
  Double_t dndyPb40[nPb40]   = {  0.32,       0.281844,   0.253366,   0.170047
                               ,  0.105019  };
  Double_t edndyPb40[nPb40]  = {  0.03,       0.0424844,  0.0423575,  0.024879
                               ,  0.0183888 };
  Double_t ryPb40[nPb40];
  Double_t rdndyPb40[nPb40];
  Double_t redndyPb40[nPb40];
  for (Int_t i = 0; i < nPb40; i++) {
    yPb40[i]      = yPb40[i] - kYbeam40;
    dndyPb40[i]  *= scalePb40;
    edndyPb40[i] *= scalePb40;
  }
  for (Int_t i = 0; i < nPb40; i++) {
    ryPb40[i]     = 0.0 - yPb40[nPb40-1-i];
    rdndyPb40[i]  =  dndyPb40[nPb40-1-i];
    redndyPb40[i] = edndyPb40[nPb40-1-i];
  }

  TGraphErrors *grPb40 = new TGraphErrors(nPb40,ryPb40, rdndyPb40
                                               ,eyPb40,redndyPb40);
  TGraphErrors *gPb40  = new TGraphErrors(nPb40, yPb40,  dndyPb40
                                               ,eyPb40, edndyPb40);

  TF1 *fPb40 = new TF1("fPb40",yFun1,-2.5,2.5,3);
  fPb40->SetParameter(0,1.0);
  fPb40->SetParameter(1,0.5);
  fPb40->SetParameter(2,0.8);

  gPb40->Fit("fPb40","0");
  fPb40->SetRange(-2.5,2.5);

  Double_t yldPb40 = fPb40->Integral(-kYbeam40,kYbeam40) / scalePb40;
  Double_t extPb40 = scalePb40 * yldPb40 / fPb40->Integral(0.0,yPb40[nPb40-1]);
  Double_t errPb40 = 0.0;
  for (Int_t i = 0; i < nPb40; i++) {
    errPb40 += edndyPb40[i]*edndyPb40[i];
  }
  errPb40 = TMath::Sqrt(errPb40) * extPb40 / scalePb40;

  //____________________________________________________________________________
  //

  //
  // 80AGeV
  //

  // Protons
  const Int_t nPr80 = 5;
  Double_t yPr80[nPr80]      = {  2.56,       2.9,        3.1,        3.3
                               ,  3.5       };
  Double_t eyPr80[nPr80]     = {  0.0,        0.0,        0.0,        0.0
                               ,  0.0       };
  Double_t dndyPr80[nPr80]   = { 30.1,       30.3989,    31.3234,    31.3663
                               , 32.6559    };
  Double_t edndyPr80[nPr80]  = {  1.0,        0.400526,   0.404057,   0.402713
                               ,  0.457937  };
  Double_t ryPr80[nPr80];
  Double_t rdndyPr80[nPr80];
  Double_t redndyPr80[nPr80];
  for (Int_t i = 0; i < nPr80; i++) {
    yPr80[i]      = yPr80[i] - kYbeam80;
    dndyPr80[i]  *= scalePr80;
    edndyPr80[i] *= scalePr80;
  }
  for (Int_t i = 0; i < nPr80; i++) {
    ryPr80[i]     = 0.0 - yPr80[nPr80-1-i];
    rdndyPr80[i]  =  dndyPr80[nPr80-1-i];
    redndyPr80[i] = edndyPr80[nPr80-1-i];
  }

  TGraphErrors *grPr80 = new TGraphErrors(nPr80,ryPr80, rdndyPr80
                                               ,eyPr80,redndyPr80);
  TGraphErrors *gPr80  = new TGraphErrors(nPr80, yPr80,  dndyPr80
                                               ,eyPr80, edndyPr80);

  TF1 *fPr80 = new TF1("fPr80",yFun1,-2.5,2.5,3);
  fPr80->SetParameter(0,1.0);
  fPr80->SetParameter(1,0.5);
  fPr80->SetParameter(2,0.8);
  fPr80->SetParLimits(1,0.4,0.8);
  fPr80->SetParLimits(2,0.6,1.2);

  gPr80->Fit("fPr80","0");
  fPr80->SetRange(-2.5,2.5);

  Double_t yldPr80 = fPr80->Integral(-kYbeam80,kYbeam80) / scalePr80;
  Double_t extPr80 = scalePr80 * yldPr80 / fPr80->Integral(0.0,yPr80[nPr80-1]);
  Double_t errPr80 = 0.0;
  for (Int_t i = 0; i < nPr80; i++) {
    errPr80 += edndyPr80[i]*edndyPr80[i];
  }
  errPr80 = TMath::Sqrt(errPr80) * extPr80 / scalePr80;

  // Anti-protons
  const Int_t nPb80 = 5;
  Double_t yPb80[nPb80]      = {  2.56,       2.9,        3.1,        3.3
                               ,  3.5       };
  Double_t eyPb80[nPb80]     = {  0.0,        0.0,        0.0,        0.0
                               ,  0.0       };
  Double_t dndyPb80[nPb80]   = {  0.87,       0.899204,   0.879943,   0.669104
                               ,  0.4626    };
  Double_t edndyPb80[nPb80]  = {  0.07,       0.041448,   0.0416276,  0.0427043
                               ,  0.0356698 };
  Double_t ryPb80[nPb80];
  Double_t rdndyPb80[nPb80];
  Double_t redndyPb80[nPb80];
  for (Int_t i = 0; i < nPb80; i++) {
    yPb80[i]      = yPb80[i] - kYbeam80;
    dndyPb80[i]  *= scalePb80;
    edndyPb80[i] *= scalePb80;
  }
  for (Int_t i = 0; i < nPb80; i++) {
    ryPb80[i]     = 0.0 - yPb80[nPb80-1-i];
    rdndyPb80[i]  =  dndyPb80[nPb80-1-i];
    redndyPb80[i] = edndyPb80[nPb80-1-i];
  }

  TGraphErrors *grPb80 = new TGraphErrors(nPb80,ryPb80, rdndyPb80
                                               ,eyPb80,redndyPb80);
  TGraphErrors *gPb80  = new TGraphErrors(nPb80, yPb80,  dndyPb80
                                               ,eyPb80, edndyPb80);

  TF1 *fPb80 = new TF1("fPb80",yFun1,-2.5,2.5,3);
  fPb80->SetParameter(0,1.0);
  fPb80->SetParameter(1,0.5);
  fPb80->SetParameter(2,0.5);
  fPb80->SetParLimits(1,0.48,1.00);
  fPb80->SetParLimits(2,0.30,0.60);

  gPb80->Fit("fPb80","0");
  fPb80->SetRange(-2.5,2.5);

  Double_t yldPb80 = fPb80->Integral(-kYbeam80,kYbeam80) / scalePb80;
  Double_t extPb80 = scalePb80 * yldPb80 / fPb80->Integral(0.0,yPb80[nPb80-1]);
  Double_t errPb80 = 0.0;
  for (Int_t i = 0; i < nPb80; i++) {
    errPb80 += edndyPb80[i]*edndyPb80[i];
  }
  errPb80 = TMath::Sqrt(errPb80) * extPb80 / scalePb80;

  //____________________________________________________________________________
  //

  TCanvas *c1 = new TCanvas("c1","", 50, 50,600,700);

  TH2F *hN1 = new TH2F("hN1","",100,-2.5,2.5,100,0.0,190.0);
  hN1->SetNdivisions(504,"X");
  hN1->SetNdivisions(505,"Y");
  hN1->SetXTitle("y");
  hN1->SetYTitle("dn/dy");
  hN1->Draw();

  fPr20->Draw("SAME");
  fPr20->SetLineColor(35);
  fPr30->Draw("SAME");
  fPr30->SetLineColor(35);
  fPr40->Draw("SAME");
  fPr40->SetLineColor(35);
  fPr80->Draw("SAME");
  fPr80->SetLineColor(35);

  grPr20->SetMarkerStyle(29);
  grPr20->SetMarkerColor(0);
  grPr20->SetMarkerSize(marker_size);
  grPr20->SetLineColor(2);
  grPr20->Draw("P");
  grPr20 = (TGraphErrors *) grPr20->Clone();
  grPr20->SetMarkerStyle(30);
  grPr20->SetMarkerColor(2);
  grPr20->Draw("P");
  gPr20->SetMarkerStyle(29);
  gPr20->SetMarkerColor(2);
  gPr20->SetMarkerSize(marker_size);
  gPr20->SetLineColor(2);
  gPr20->Draw("P");

  grPr30->SetMarkerStyle(22);
  grPr30->SetMarkerColor(0);
  grPr30->SetMarkerSize(marker_size);
  grPr30->SetLineColor(2);
  grPr30->Draw("P");
  grPr30 = (TGraphErrors *) grPr30->Clone();
  grPr30->SetMarkerStyle(26);
  grPr30->SetMarkerColor(2);
  grPr30->Draw("P");
  gPr30->SetMarkerStyle(22);
  gPr30->SetMarkerColor(2);
  gPr30->SetMarkerSize(marker_size);
  gPr30->SetLineColor(2);
  gPr30->Draw("P");

  grPr40->SetMarkerStyle(21);
  grPr40->SetMarkerColor(0);
  grPr40->SetMarkerSize(marker_size);
  grPr40->SetLineColor(2);
  grPr40->Draw("P");
  grPr40 = (TGraphErrors *) grPr40->Clone();
  grPr40->SetMarkerStyle(25);
  grPr40->SetMarkerColor(2);
  grPr40->Draw("P");
  gPr40->SetMarkerStyle(21);
  gPr40->SetMarkerColor(2);
  gPr40->SetMarkerSize(marker_size);
  gPr40->SetLineColor(2);
  gPr40->Draw("P");

  grPr80->SetMarkerStyle(20);
  grPr80->SetMarkerColor(0);
  grPr80->SetMarkerSize(marker_size);
  grPr80->SetLineColor(2);
  grPr80->Draw("P");
  grPr80 = (TGraphErrors *) grPr80->Clone();
  grPr80->SetMarkerStyle(24);
  grPr80->SetMarkerColor(2);
  grPr80->Draw("P");
  gPr80->SetMarkerStyle(20);
  gPr80->SetMarkerColor(2);
  gPr80->SetMarkerSize(marker_size);
  gPr80->SetLineColor(2);
  gPr80->Draw("P");

  TLegendEntry *entry;
  TLegend      *leg1 = new TLegend(0.64,0.78,0.85,0.94);
  leg1->SetFillColor(0);
  leg1->SetBorderSize(0);
  leg1->SetTextSize(0.03);
  entry = leg1->AddEntry(gPr20,"20A GeV (#times 3.0)","P");
  entry->SetTextColor(1);
  entry = leg1->AddEntry(gPr30,"30A GeV (#times 2.0)","P");
  entry->SetTextColor(1);
  entry = leg1->AddEntry(gPr40,"40A GeV (#times 1.5)","P");
  entry->SetTextColor(1);
  entry = leg1->AddEntry(gPr80,"80A GeV (#times 1.0)","P");
  entry->SetTextColor(1);
  leg1->Draw();

  printf("\n");
  printf(" ---> pr-yield (20A GeV) = %f +- %f\n",yldPr20,errPr20);
  printf(" ---> pr-yield (30A GeV) = %f +- %f\n",yldPr30,errPr30);
  printf(" ---> pr-yield (40A GeV) = %f +- %f\n",yldPr40,errPr40);
  printf(" ---> pr-yield (80A GeV) = %f +- %f\n",yldPr80,errPr80);
  printf("\n");

  //____________________________________________________________________________
  //

  TCanvas *c2 = new TCanvas("c2","",100,100,600,700);

  TH2F *hN2 = new TH2F("hN2","",100,-2.5,2.5,100,0.0,1.2);
  hN2->SetNdivisions(504,"X");
  hN2->SetNdivisions(505,"Y");
  hN2->SetXTitle("y");
  hN2->SetYTitle("dn/dy");
  hN2->Draw();

  fPb40->Draw("SAME");
  fPb40->SetLineColor(35);
  fPb80->Draw("SAME");
  fPb80->SetLineColor(35);

  grPb40->SetMarkerStyle(21);
  grPb40->SetMarkerColor(0);
  grPb40->SetMarkerSize(marker_size);
  grPb40->SetLineColor(2);
  grPb40->Draw("P");
  grPb40 = (TGraphErrors *) grPb40->Clone();
  grPb40->SetMarkerStyle(25);
  grPb40->SetMarkerColor(2);
  grPb40->Draw("P");
  gPb40->SetMarkerStyle(21);
  gPb40->SetMarkerColor(2);
  gPb40->SetMarkerSize(marker_size);
  gPb40->SetLineColor(2);
  gPb40->Draw("P");

  grPb80->SetMarkerStyle(20);
  grPb80->SetMarkerColor(0);
  grPb80->SetMarkerSize(marker_size);
  grPb80->SetLineColor(2);
  grPb80->Draw("P");
  grPb80 = (TGraphErrors *) grPb80->Clone();
  grPb80->SetMarkerStyle(24);
  grPb80->SetMarkerColor(2);
  grPb80->Draw("P");
  gPb80->SetMarkerStyle(20);
  gPb80->SetMarkerColor(2);
  gPb80->SetMarkerSize(marker_size);
  gPb80->SetLineColor(2);
  gPb80->Draw("P");

  TLegend      *leg2 = new TLegend(0.72,0.86,0.85,0.94);
  leg2->SetFillColor(0);
  leg2->SetBorderSize(0);
  leg2->SetTextSize(0.03);
  entry = leg2->AddEntry(gPr40," 40A GeV","P");
  entry->SetTextColor(1);
  entry = leg2->AddEntry(gPr80," 80A GeV","P");
  entry->SetTextColor(1);
  leg2->Draw();

  printf("\n");
  printf(" ---> pb-yield (40A GeV) = %f +- %f\n",yldPb40,errPb40);
  printf(" ---> pb-yield (80A GeV) = %f +- %f\n",yldPb80,errPb80);
  printf("\n");

}

//______________________________________________________________________________
Double_t yFun1(Double_t *x, Double_t *p)
{

  Double_t y = x[0];
  Double_t c = p[0];
  Double_t w = p[1];
  Double_t s = p[2];

  Double_t f = c * (TMath::Exp(-(y-s)*(y-s)/(2.0*w*w))
                  + TMath::Exp(-(y+s)*(y+s)/(2.0*w*w)));

  return f;

}