So I want to fit the data with:
Dn/Dy = f ( y - y_beam) + f (y_beam - y)
where f is the Guass function.
So say, I want to fit the data of 40GeV energy I would have to fit with
Dn/Dy = f (y_40 - y_beam_40) +f (y_beam_40 - y_40)
In the following code I get the data and it also has a guass fit function. However, I don’t really know what the fit function (see function yFun1 towards the end of the code) is doing. I think its just fitting the y values. Can anyone tell me how to change the fit function so that I can get the fit that I want to. I have copied the macro and also attached it this post.
I would be really grateful if someone could help me out. Thanks !!
Note in the function yFun1 - the formula is for the Guassian fit.
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.
;
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.;
fPr80->SetParLimits(1,0.4,0.;
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.0ww))
+ TMath::Exp(-(y+s)(y+s)/(2.0ww)));
return f;
}
pr_pb_yspectra.C (16.7 KB)