return norm * pt * integral;
}
+Double_t
+BGBlastWave_Func_OneOverPt(const Double_t *x, const Double_t *p)
+{
+ /* 1/pt dN/dpt */
+
+ Double_t pt = x[0];
+ Double_t mass = p[0];
+ Double_t mt = TMath::Sqrt(pt * pt + mass * mass);
+ Double_t beta_max = p[1];
+ Double_t temp = p[2];
+ Double_t n = p[3];
+ Double_t norm = p[4];
+
+ if (!fBGBlastWave_Integrand)
+ fBGBlastWave_Integrand = new TF1("fBGBlastWave_Integrand", BGBlastWave_Integrand, 0., 1., 5);
+ fBGBlastWave_Integrand->SetParameters(mt, pt, beta_max, temp, n);
+ Double_t integral = fBGBlastWave_Integrand->Integral(0., 1.);
+
+ return norm * integral;
+}
+
+
TF1 *
BGBlastWave(const Char_t *name, Double_t mass, Double_t beta_max = 0.9, Double_t temp = 0.1, Double_t n = 1., Double_t norm = 1.e6)
{
fBGBlastWave->FixParameter(0, mass);
fBGBlastWave->SetParLimits(1, 0.01, 0.99);
fBGBlastWave->SetParLimits(2, 0.01, 1.);
- fBGBlastWave->SetParLimits(3, 0.1, 10.);
+ fBGBlastWave->SetParLimits(3, 0.01, 10.);
+ return fBGBlastWave;
+}
+
+TF1 * BGBlastWave_OneOverPT(const Char_t *name, Double_t mass, Double_t beta_max = 0.9, Double_t temp = 0.1, Double_t n = 1., Double_t norm = 1.e6)
+{
+
+ TF1 *fBGBlastWave = new TF1(name, BGBlastWave_Func_OneOverPt, 0., 10., 5);
+ fBGBlastWave->SetParameters(mass, beta_max, temp, n, norm);
+ fBGBlastWave->SetParNames("mass", "beta_max", "T", "n", "norm");
+ fBGBlastWave->FixParameter(0, mass);
+ fBGBlastWave->SetParLimits(1, 0.01, 0.99);
+ fBGBlastWave->SetParLimits(2, 0.01, 1.);
+ fBGBlastWave->SetParLimits(3, 0.01, 10.);
return fBGBlastWave;
}
TGraphErrors *gBW[1000];
TObjArray *
-BGBlastWave_GlobalFit(TObjArray *data, Double_t *mass, Double_t profile = 1., Bool_t fixProfile = kFALSE)
+BGBlastWave_GlobalFit(TObjArray *data, Double_t *mass, Double_t profile = .9, Bool_t fixProfile = kFALSE)
{
/* get data */
minuit->mnexcm("SET ERR", arglist, 1, ierflg);
for (Int_t idata = 0; idata < nBW; idata++)
minuit->mnparm(idata, Form("norm%d", idata), 1.e6, 1., 0., 0., ierflg);
- minuit->mnparm(nBW + 0, "<beta>", 0.5, 0.1, 0., 1., ierflg);
- minuit->mnparm(nBW + 1, "T", 0.2, 0.1, 0., 1., ierflg);
+ // minuit->mnparm(nBW + 0, "<beta>", 0.65, 0.01, 0., 1., ierflg);
+ // minuit->mnparm(nBW + 1, "T", 0.1, 0.01, 0., 1., ierflg);
+ minuit->mnparm(nBW + 0, "<beta>", 0.55, 0.01, 0., 1., ierflg);
+ minuit->mnparm(nBW + 1, "T", 0.13, 0.01, 0., 1., ierflg);
minuit->mnparm(nBW + 2, "n", profile, 0.1, 0., 10., ierflg);
if (fixProfile) minuit->FixParameter(nBW + 2);
/* 1-sigma contour */
minuit->SetErrorDef(1);
TGraph *gCont1 = NULL;
- gCont1 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
+ // gCont1 = (TGraph *) minuit->Contour(50, nBW + 0, nBW + 1);
if (gCont1) gCont1->SetName("gCont1");
/* 2-sigma contour */
/*****************************************************************/
-GetYieldAndMean(TH1 *h, TF1 *f, Double_t &yield, Double_t &yielderr, Double_t &mean, Double_t &meanerr, Double_t min, Double_t max, Double_t *partyield, Double_t *partyielderr)
+IntegratedProduction(TH1 *h, Double_t mass, Option_t *opt = "")
+{
+
+ Double_t yield, yielderr, yielderrcorr, mean, meanerr, meanerrcorr, partyield[3], partyielderr[3], partyielderrcorr[3];
+ TF1 *f = BGBlastWave_SingleFit(h, mass, opt);
+ GetYieldAndMean(h, f, yield, yielderr, yielderrcorr, mean, meanerr, meanerrcorr, 0., 10., partyield, partyielderr, partyielderrcorr);
+
+ // Double_t fint = f->Integral(0.,10.);
+ // Double_t finte = f->IntegralError(0.,10.);
+ // Double_t fmean = f->Mean(0., 10.);
+
+ printf("dN/dy = %f +- %f (%f)\n", yield, yielderr, yielderrcorr);
+ printf("<pt> = %f +- %f (%f)\n", mean, meanerr, meanerrcorr);
+ printf("dN/dy (data) = %f +- %f (%f)\n", partyield[0], partyielderr[0], partyielderrcorr[0]);
+ printf("dN/dy (low) = %f +- %f (%f)\n", partyield[1], partyielderr[1], partyielderrcorr[1]);
+ printf("dN/dy (high) = %f +- %f (%f)\n", partyield[2], partyielderr[2], partyielderrcorr[2]);
+ // printf("----\n");
+ // printf("dN/dy (func) = %f +- %f\n", fint, finte);
+ // printf("<pT> (func) = %f +- %f\n", fmean, 0.);
+
+ // TH1 *hr = (TH1 *)h->Clone("hr");
+ // hr->Divide(f);
+ // new TCanvas;
+ // hr->Draw();
+
+ // TProfile *p = new TProfile("p", "", 100, 0., 10.);
+ // gROOT->LoadMacro("HistoUtils.C");
+ // HistoUtils_Function2Profile(f, p);
+ // p->Draw();
+}
+
+GetYieldAndMean(TH1 *h, TF1 *f, Double_t &yield, Double_t &yielderr, Double_t &yielderrcorr, Double_t &mean, Double_t &meanerr, Double_t &meanerrcorr, Double_t min, Double_t max, Double_t *partyield, Double_t *partyielderr, Double_t *partyielderrcorr)
{
/* find lowest edge in histo */
}
/* integrate the data */
- Double_t cont, err, width, cent, integral_data = 0., integralerr_data = 0., meanintegral_data = 0., meanintegralerr_data = 0.;
+ Double_t cont, err, width, cent, integral_data = 0., integralerr_data = 0., integralerrcorr_data = 0., meanintegral_data = 0., meanintegralerr_data = 0., meanintegralerrcorr_data = 0.;
for (Int_t ibin = binlo; ibin < binhi; ibin++) {
cent = h->GetBinCenter(ibin);
width = h->GetBinWidth(ibin);
/* all right, use data */
integral_data += cont * width;
integralerr_data += err * err * width * width;
+ integralerrcorr_data += err * width;
meanintegral_data += cont * width * cent;
meanintegralerr_data += err * err * width * width * cent * cent;
+ meanintegralerrcorr_data += err * width * cent;
}
else {
/* missing data-point, complain and use function */
yielderr = TMath::Sqrt(integralerr_data * integralerr_data +
integralerr_lo * integralerr_lo +
integralerr_hi * integralerr_hi);
+ yielderrcorr = TMath::Sqrt(integralerrcorr_data * integralerrcorr_data +
+ integralerr_lo * integralerr_lo +
+ integralerr_hi * integralerr_hi);
/* compute integrated mean */
mean = (meanintegral_data + meanintegral_lo + meanintegral_hi) / yield;
meanerr = TMath::Sqrt(meanintegralerr_data * meanintegralerr_data +
meanintegralerr_lo * meanintegralerr_lo +
meanintegralerr_hi * meanintegralerr_hi) / yield;
+ meanerrcorr = TMath::Sqrt(meanintegralerrcorr_data * meanintegralerrcorr_data +
+ meanintegralerr_lo * meanintegralerr_lo +
+ meanintegralerr_hi * meanintegralerr_hi) / yield;
/* set partial yields */
partyield[0] = integral_data;
partyielderr[0] = integralerr_data;
+ partyielderrcorr[0] = integralerrcorr_data;
partyield[1] = integral_lo;
partyielderr[1] = integralerr_lo;
+ partyielderrcorr[1] = integralerr_lo;
partyield[2] = integral_hi;
partyielderr[2] = integralerr_hi;
+ partyielderrcorr[2] = integralerr_hi;
}
}
TH1 *
-Convert_dNdy_1over2pipt_dNdeta(TH1 *hin, Double_t mass, Double_t y = 0.5)
+Convert_dNdy_1over2pipt_dNdeta(TH1 *hin, Double_t mass, Double_t eta = 0.8)
+{
+
+ TH1 *hout = hin->Clone("hout");
+ hout->Reset();
+ Double_t pt, mt, conv, val, vale;
+ for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+ pt = hin->GetBinCenter(ibin + 1);
+ conv = eta2y(pt, mass, eta) / eta;
+ val = hin->GetBinContent(ibin + 1);
+ vale = hin->GetBinError(ibin + 1);
+ val /= (2. * TMath::Pi() * pt);
+ vale /= (2. * TMath::Pi() * pt);
+ val *= conv;
+ vale *= conv;
+ hout->SetBinContent(ibin + 1, val);
+ hout->SetBinError(ibin + 1, vale);
+ }
+ return hout;
+}
+
+TH1 *
+Convert_dNdy_1over2pipt_dNdy(TH1 *hin)
{
TH1 *hout = hin->Clone("hout");
Double_t pt, mt, conv, val, vale;
for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
pt = hin->GetBinCenter(ibin + 1);
- conv = y / y2eta(pt, mass, y);
val = hin->GetBinContent(ibin + 1);
vale = hin->GetBinError(ibin + 1);
val /= (2. * TMath::Pi() * pt);
vale /= (2. * TMath::Pi() * pt);
+ hout->SetBinContent(ibin + 1, val);
+ hout->SetBinError(ibin + 1, vale);
+ }
+ return hout;
+}
+
+TH1 *
+Convert_dNdy_dNdeta(TH1 *hin, Double_t mass, Double_t eta = 0.8)
+{
+
+ TH1 *hout = hin->Clone("hout");
+ hout->Reset();
+ Double_t pt, mt, conv, val, vale;
+ for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
+ pt = hin->GetBinCenter(ibin + 1);
+ conv = eta2y(pt, mass, eta) / eta;
+ val = hin->GetBinContent(ibin + 1);
+ vale = hin->GetBinError(ibin + 1);
val *= conv;
vale *= conv;
hout->SetBinContent(ibin + 1, val);
return hout;
}
+TGraph *
+Convert_dNdy_dNdeta(TGraph *hin, Double_t mass, Double_t eta = 0.8)
+{
+
+ TGraph *hout = hin->Clone("hout");
+ // hout->Reset();
+ Double_t pt, mt, conv, val, valelo, valehi;
+ for (Int_t ibin = 0; ibin < hin->GetN(); ibin++) {
+ pt = hin->GetX()[ibin];
+ conv = eta2y(pt, mass, eta) / eta;
+ val = hin->GetY()[ibin];
+ valelo = hin->GetEYlow()[ibin];
+ valehi = hin->GetEYhigh()[ibin];
+ val *= conv;
+ valelo *= conv;
+ valehi *= conv;
+ hout->GetX()[ibin] = pt;
+ hout->GetY()[ibin] = val;
+ hout->GetEYlow()[ibin] = valelo;
+ hout->GetEYhigh()[ibin] = valehi;
+ }
+ return hout;
+}
+
TH1 *
SummedId_1over2pipt_dNdeta(const Char_t *filename, Int_t icent)
{
return hsum;
}
+TH1 *
+SummedId_dNdeta(const Char_t *filename, Int_t icent)
+{
+
+ const Char_t *chargeName[2] = {
+ "plus", "minus"
+ };
+
+ TFile *filein = TFile::Open(filename);
+ TH1 *hy[AliPID::kSPECIES][2];
+ TH1 *heta[AliPID::kSPECIES][2];
+ for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
+ for (Int_t icharge = 0; icharge < 2; icharge++) {
+ hy[ipart][icharge] = (TH1 *)filein->Get(Form("cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
+ if (!hy[ipart][icharge]) {
+ printf("cannot find cent%d_%s_%s\n", icent, AliPID::ParticleName(ipart), chargeName[icharge]);
+ return NULL;
+ }
+ heta[ipart][icharge] = Convert_dNdy_dNdeta(hy[ipart][icharge], AliPID::ParticleMass(ipart));
+ }
+
+ /* sum */
+ TH1D *hsum = heta[2][0]->Clone("hsum");
+ hsum->Reset();
+ for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++)
+ for (Int_t icharge = 0; icharge < 2; icharge++)
+ hsum->Add(heta[ipart][icharge]);
+
+ return hsum;
+}
+
/*****************************************************************/
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff