// Piotr.Skowronski@cern.ch
//
///////////////////////////////////////////////////////////////////////////////////
+#include "AliHBTAnalysis.h"
#include <TROOT.h>
#include <TH1.h>
TF1* AliHBTFits::fgF1 = 0x0;
TF1* AliHBTFits::fgF2 = 0x0;
+Double_t *AliHBTFits::fXbins = 0;
+Double_t *AliHBTFits::fYbins = 0;
+Double_t *AliHBTFits::fZbins = 0;
+
+Double_t *AliHBTFits::fExpT = 0;
+Double_t *AliHBTFits::fLongT = 0;
+Double_t *AliHBTFits::fIntegrT = 0;
+
+Int_t AliHBTFits::fNX = 0;
+Int_t AliHBTFits::fNY = 0;
+Int_t AliHBTFits::fNZ = 0;
+ Int_t AliHBTFits::fNXY = 0;
+Int_t AliHBTFits::fNT = 0;
+Int_t AliHBTFits::fBinX = 0;
+Int_t AliHBTFits::fBinY = 0;
+Int_t AliHBTFits::fBinZ = 0;
+
+Bool_t AliHBTFits::fLongOnly = 1;
/******************************************************************************************/
AliHBTFits::~AliHBTFits()
theory->SetParName(2,"J");
theory->SetParName(3,"L");
theory->SetParName(4,"\\Phi");
+
+
+/***********************************/
+ /***********************************/
+ /***********************************/
+
+ Double_t xmin = 0.0;
+ Double_t ymin = 0.0;
+ Double_t zmin = 0.0;
+
+ Int_t hxfirst = h->GetXaxis()->FindFixBin(xmin); if (hxfirst < 1) hxfirst = 1;
+ Int_t hxlast = h->GetXaxis()->FindFixBin(xmax); if (hxlast > h->GetXaxis()->GetLast()) hxlast = h->GetXaxis()->GetLast();
+ Int_t hyfirst = h->GetYaxis()->FindFixBin(ymin); if (hyfirst < 1) hyfirst = 1;
+ Int_t hylast = h->GetYaxis()->FindFixBin(ymax); if (hylast > h->GetYaxis()->GetLast()) hylast = h->GetYaxis()->GetLast();
+ Int_t hzfirst = h->GetZaxis()->FindFixBin(zmin); if (hzfirst < 1) hzfirst = 1;
+ Int_t hzlast = h->GetZaxis()->FindFixBin(zmax); if (hzlast > h->GetZaxis()->GetLast()) hzlast = h->GetZaxis()->GetLast();
+
+ fNX = hxlast-hxfirst+1;
+ fNY = hylast-hyfirst+1;
+ fNZ = hzlast-hzfirst+1;
+
+ fXbins = new Double_t[fNX];
+ fYbins = new Double_t[fNY];
+ fZbins = new Double_t[fNZ];
+ Int_t c = 0;
+ for (Int_t k=hzfirst; k<=hzlast; k++)
+ {
+ if (h->GetZaxis()->GetBinCenter(k) > zmax) break;
+ fZbins[c++] = TMath::Abs(h->GetZaxis()->GetBinCenter(k));
+ }
+ fNZ = c;
+ c = 0;
+ for (Int_t j=hyfirst; j<=hylast; j++)
+ {
+ if (h->GetYaxis()->GetBinCenter(j) > ymax) break;
+ fYbins[c++] = TMath::Abs(h->GetYaxis()->GetBinCenter(j));
+ }
+ fNY = c;
+ c = 0;
+ for (Int_t i=hxfirst; i<=hxlast; i++)
+ {
+
+ if (h->GetXaxis()->GetBinCenter(i) > xmax) break;
+ fXbins[c++] = TMath::Abs(h->GetXaxis()->GetBinCenter(i));
+ ::Info("","%d %d %d %f",i,c-1,fNX,fXbins[c-1]);
+ }
+ fNX = c;
+
+ fNXY=fNX*fNY;
+ fNT=fNXY*fNZ;
+
+ fExpT = new Double_t[fNT];
+ fLongT = new Double_t[fNZ];
+ fIntegrT = new Double_t[fNXY];
+
+ fBinX = 0;
+ fBinY = 0;
+ fBinZ = 0;
+ /***********************************/
+ /***********************************/
+/***********************************/
- h->Fit(theory,fopt,"E");
+ h->Fit(theory,fopt,"E");
+
+ delete [] fXbins;
+ delete [] fYbins;
+ delete [] fZbins;
+
+ delete [] fExpT;
+ delete [] fLongT;
+ delete [] fIntegrT;
+
+ fXbins = 0x0;
+ fYbins = 0x0;
+ fZbins = 0x0;
+ fExpT = 0x0;
+ fLongT = 0x0;
+ fIntegrT = 0x0;
}
Double_t AliHBTFits::QOutQSideQLongCylSurf(Double_t *x, Double_t *par)
{
//QOutQSideQLongCylSurf Function
+
+ if (fXbins == 0x0)
+ {
+ ::Info("QOutQSideQLongCylSurf","arrays deleted");
+ return 1.0;
+ }
+
+
Double_t qout = x[0];
Double_t qside = x[1];
Double_t qlong = x[2];
+
+// ::Info("QOutQSideQLongCylSurf","out=%f, side=%f, long=%f %d %d %d",qout,qside,qlong,fBinX,fBinY,fBinZ);
+
+ if (fXbins[fBinX] != x[0])
+ {
+ ::Fatal("QOutQSideQLongCylSurf","expected X %f got %f",fXbins[fBinX],qout);
+ return 1.0;
+ }
+
+
+ if (fXbins[fBinY] != x[1])
+ {
+ ::Fatal("QOutQSideQLongCylSurf","expected Y %f got %f",fYbins[fBinY],qside);
+ return 1.0;
+ }
+
+ if (fXbins[fBinZ] != x[2])
+ {
+ ::Fatal("QOutQSideQLongCylSurf","expected Z %f got %f",fZbins[fBinZ],qlong);
+ return 1.0;
+ }
Double_t aLambda = par[0];
Double_t aR = par[1];
Double_t aL = par[3];
Double_t aPhi = par[4];
- static Double_t oldlam = aLambda;
- static Double_t oldPhi = aPhi;
-
- if (oldPhi != aPhi)
- {
- ::Info("QOutQSideQLongCylSurf","out=%f, side=%f, long=%f\n lambda=%f, R=%f, J=%f, L=%f Phi=%f",
- qout,qside,qlong,aLambda,aR,aJ,aL,aPhi);
- oldPhi = aPhi;
- }
+ static Double_t oldLambda = -1;
+ static Double_t oldR = -1;
+ static Double_t oldJ = -1;
+ static Double_t oldL = -1;
+ static Double_t oldPhi = -1;
- if (oldlam != aLambda)
+
+
+ if ( (oldLambda != aLambda) || (oldR != aR) || (oldL != aL))
{
::Info("QOutQSideQLongCylSurf","out=%f, side=%f, long=%f\n lambda=%f, R=%f, J=%f, L=%f Phi=%f",
qout,qside,qlong,aLambda,aR,aJ,aL,aPhi);
- oldlam = aLambda;
+ oldLambda = aLambda;
}
+
+ Double_t result = aLambda;
- Double_t result=aLambda*TMath::Exp(-(qout*qout+qside*qside+qlong*qlong)*aJ*aJ/0.0389273);
+ Double_t exppart;
- if ( (qlong != 0.0) && (aL != 0.0) )
+ if (aJ == oldJ)
{
- Double_t qlL = qlong*aL/2.0;
- Double_t sinqlL = TMath::Sin(qlL);
- Double_t sin2qlL = sinqlL*sinqlL;
- Double_t longpart = sin2qlL/(qlL*qlL);
- result*= longpart;
+ exppart = fExpT[fBinZ*fNXY + fBinY*fNX + fBinX];
}
+ else
+ {
+ Int_t cbin = 0;
+ Double_t aJ2=aJ*aJ/0.0389273;
+ for (Int_t k=0; k<fNZ; k++)
+ {
+ Double_t qL2 = fZbins[k]*fZbins[k];
+
+ for (Int_t j=0; j<fNY; j++)
+ {
+ Double_t qS2 = fYbins[j]*fYbins[j];
+
+ for (Int_t i=0; i<fNX; i++)
+ {
+ fExpT[cbin++] = TMath::Exp(-(fXbins[i]*fXbins[i]+qS2+qL2)*aJ2);
+ }
+ }
+ }
+ exppart = fExpT[0];
+ oldJ = aJ;
+ }
+
+ result = aLambda*exppart;
+// ::Info("QOSL","lambda*exp=%f",result);
+
+ Double_t longpart = 1.0;
- if (aPhi > 0.0)
+ if (oldL == aL)
{
- fgF1->SetParameter(0,aR);
- fgF2->SetParameter(0,aR);
+ longpart = fLongT[fBinZ];
+ }
+ else
+ {
+ if ( aL <= 0.0 )
+ {
+ for (Int_t i=0; i<fNZ; i++) fLongT[i] = 1.0;
+ }
+ else
+ {
+ for (Int_t i=0; i<fNZ; i++)
+ {
+ Double_t qlL = fZbins[i]*aL/(2.0*0.197);
+ Double_t sinqlL = TMath::Sin(qlL);
+ Double_t sin2qlL = sinqlL*sinqlL;
+ fLongT[i] = sin2qlL/(qlL*qlL);
+ }
+ }
+ longpart = fLongT[0];
+ oldL = aL;
+ }
+
+// ::Info("QOSL","longpart=%f",longpart);
+
+ result *=longpart;
+ Double_t itgrl;
- fgF1->SetParameter(1,aPhi);
- fgF2->SetParameter(1,aPhi);
+ if ( (oldR != aR) || (oldPhi !=aPhi) )
+ {
+ if (aPhi <= 0.0)
+ {
+ for (Int_t i=0; i<fNXY; i++)
+ {
+ fIntegrT[i] = 1.0;
+ }
+ }
+ else
+ {
+ Int_t cbin = 0;
+ for (Int_t j=0; j<fNY; j++)
+ {
+ for (Int_t i=0; i<fNX; i++)
+ {
+
+ fgF1->SetParameter(0,aR);
+ fgF2->SetParameter(0,aR);
+
+ fgF1->SetParameter(1,aPhi);
+ fgF2->SetParameter(1,aPhi);
+
+ fgF1->SetParameter(2,fXbins[i]);
+ fgF2->SetParameter(2,fXbins[i]);
+
+ fgF1->SetParameter(3,fYbins[j]);
+ fgF2->SetParameter(3,fYbins[j]);
+
+ Double_t sp = fgF1->Integral(-TMath::Pi(),TMath::Pi(),(Double_t*)0x0,1.e-8);
+ Double_t cp = fgF2->Integral(-TMath::Pi(),TMath::Pi(),(Double_t*)0x0,1.e-8);
+
+ // ::Info("QSideCylSurf","sp=%f, cp=%f",cp,sp);
+
+ fIntegrT[cbin] = cp*cp + sp*sp;
+ // ::Info("QSideCylSurf","2: result=%f",result);
+
+ Double_t oneOverPhi = 1./aPhi;
+ Double_t tmp = TMath::BesselI0(oneOverPhi)-TMath::StruveL0(oneOverPhi);
+ tmp = tmp*TMath::TwoPi();
+ tmp = tmp*tmp;
+ if (tmp == 0.0)
+ {
+ ::Error("QOutCylSurf","Div by 0");
+ return 1.0;
+ }
+ fIntegrT[cbin] /= tmp;
+// ::Info("integr calc","%d (%d %d) %f",cbin,i,j,fIntegrT[cbin]);
+ cbin++;
+ }
+ }
+
+ }
+ itgrl = fIntegrT[0];
+ oldR = aR;
+ oldPhi = aPhi;
+ }
+ else
+ {
+ itgrl = fIntegrT[fBinY*fNX+fBinX];
+ // ::Info("integr ret","%d (%d %d) %f",fBinY*fNX+fBinX,fBinX,fBinY,fIntegrT[fBinY*fNX+fBinX]);
+ }
- fgF1->SetParameter(2,qout);
- fgF2->SetParameter(2,qout);
+// ::Info("QOSL","itgrl/bessel=%f",itgrl);
- fgF1->SetParameter(3,qside);
- fgF2->SetParameter(3,qside);
-
- Double_t sp = fgF1->Integral(-TMath::Pi(),TMath::Pi(),(Double_t*)0x0,1.e-8);
- Double_t cp = fgF2->Integral(-TMath::Pi(),TMath::Pi(),(Double_t*)0x0,1.e-8);
-
-// ::Info("QSideCylSurf","sp=%f, cp=%f",cp,sp);
-
- result *= cp*cp + sp*sp;
-// ::Info("QSideCylSurf","2: result=%f",result);
-
- Double_t tmp = TMath::BesselI0(1./aPhi);
- tmp = tmp*tmp*TMath::TwoPi()*TMath::TwoPi();
- if (tmp == 0.0)
- {
- ::Error("QOutCylSurf","Div by 0");
- return 1.0;
- }
- result=result/tmp;
- }
+ result *= itgrl;
- result+=1.;
+ fBinX++;
+
+ if (fBinX == fNX)
+ {
+ fBinX = 0;
+ fBinY++;
+ if ( fBinY == fNY )
+ {
+ fBinY = 0;
+ fBinZ++;
+ if ( fBinZ == fNZ )
+ {
+ fBinZ = 0;
+ }
+ }
+ }
+// ::Info("QOutQSideQLongCylSurf Return","out=%f, side=%f, long=%f %d %d %d",qout,qside,qlong,fBinX,fBinY,fBinZ);
+
+ result+=1.;
+
// ::Info("QSideCylSurf","Final: result=%f",result);
+// ::Info("QSideCylSurf","XXX result=%f",XXX(x,par));
+// AliHBTAnalysis::PressAnyKey();
return result;
}
/******************************************************************************************/
if (oldlam != aLambda)
{
- ::Info("QOutQSideQLongCylSurf","out=%f, side=%f, lambda=%f, R=%f, J=%f, Phi=%f",
+ ::Info("QOutQSideCylSurf","out=%f, side=%f, lambda=%f, R=%f, J=%f, Phi=%f",
qout,qside,aLambda,aR,aJ,aPhi);
oldlam = aLambda;
}
result *= cp*cp + sp*sp;
// ::Info("QSideCylSurf","2: result=%f",result);
- Double_t tmp = TMath::BesselI0(1./aPhi);
- tmp = tmp*tmp*TMath::TwoPi()*TMath::TwoPi();
+ Double_t oneOverPhi = 1./aPhi;
+ Double_t tmp = TMath::BesselI0(oneOverPhi)-TMath::StruveL0(oneOverPhi);
+ tmp = tmp*TMath::TwoPi();
+ tmp = tmp*tmp;
if (tmp == 0.0)
{
::Error("QOutCylSurf","Div by 0");
if (max <= 0.0) max = h->GetXaxis()->GetXmax();
TF1 *theory = new TF1("BorisTomasikCylSurfQOut",&AliHBTFits::QOutCylSurf,0.0,max,3);//par : aR, aPhi
theory->SetNpx(1000);
+
theory->SetParameter(0,6.0);
theory->SetParameter(1,1.0);
theory->SetParameter(2,0.75);
+
+ theory->FixParameter(0,5.92);
+ theory->FixParameter(1,0.864);
+ theory->FixParameter(2,0.773);
+
theory->SetParName(0,"R_{out}");
theory->SetParName(1,"\\Phi");
theory->SetParName(2,"\\lambda");
Double_t aLambda = par[2];
// ::Info("QOutCylSurf","q=%f, lambda=%f, Rout=%f, Phi=%f",qout,aLambda,aR,aPhi);
-
- Double_t result=aLambda*TMath::Exp(-qout*qout/0.0389273);
+ static const Double_t j2=1.063*1.063;
+ Double_t result=aLambda*TMath::Exp(-qout*qout*j2/0.0389273);
// ::Info("QOutCylSurf","1: result=%f",result);
result = result * (cp*cp + sp*sp);
// ::Info("QOutCylSurf","2: result=%f",result);
- Double_t tmp = TMath::BesselI0(1./aPhi);
- tmp = tmp*tmp*TMath::TwoPi()*TMath::TwoPi();
+ Double_t oneOverPhi = 1./aPhi;
+ Double_t tmp = TMath::BesselI0(oneOverPhi)-TMath::StruveL0(oneOverPhi);
+ tmp = tmp*TMath::TwoPi();
+ tmp = tmp*tmp;
if (tmp == 0.0)
{
::Error("QOutCylSurf","Div by 0");
theory->SetParameter(0,6.0);
theory->SetParameter(1,1.0);
theory->SetParameter(2,0.75);
+
+ theory->FixParameter(0,5.92);
+ theory->FixParameter(1,0.864);
+ theory->FixParameter(2,0.773);
+
theory->SetParName(0,"R_{out}");
theory->SetParName(1,"\\Phi");
theory->SetParName(2,"\\lambda");
// ::Info("QSideCylSurf","q=%f, lambda=%f, Rside=%f, Phi=%f",qside,aLambda,aR,aPhi);
- Double_t result=aLambda*TMath::Exp(-qside*qside/0.0389273);
+ static const Double_t j2=1.063*1.063;
+ Double_t result=aLambda*TMath::Exp(-qside*qside*j2/0.0389273);
if (aPhi > 0.0)
{
fgF1->SetParameter(0,aR);
result *= cp*cp + sp*sp;
- Double_t tmp = TMath::BesselI0(1./aPhi);
- tmp = tmp*tmp*TMath::TwoPi()*TMath::TwoPi();
+ Double_t oneOverPhi = 1./aPhi;
+ Double_t tmp = TMath::BesselI0(oneOverPhi)-TMath::StruveL0(oneOverPhi);
+ tmp = tmp*TMath::TwoPi();
+ tmp = tmp*tmp;
if (tmp == 0.0)
{
- ::Error("QSideCylSurf","Div by 0");
+ ::Error("QOutCylSurf","Div by 0");
return 1.0;
}
+
result/=tmp;
}
return result + 1.;
}
/******************************************************************************************/
+void AliHBTFits::FitQLongCylSurf(const TString& hname, Option_t* fopt, Float_t max)
+{
+ //Fits QSide According to Boris Tomasik formula
+ TH1D* h = dynamic_cast<TH1D*>(gDirectory->Get(hname));
+ if (h == 0x0)
+ {
+ ::Error("FitQLongCylSurf","Can not find histogram %s",hname.Data());
+ return;
+ }
+
+ delete fgF1;
+ delete fgF2;
+ fgF1 = new TF1("ff1","sin([0]*[2]*sin(x)/0.197)*exp(cos(x)/[1])",-TMath::Pi(),TMath::Pi());
+ fgF2 = new TF1("ff2","cos([0]*[2]*sin(x)/0.197)*exp(cos(x)/[1])",-TMath::Pi(),TMath::Pi());
+ fgF1->SetNpx(1000);
+ fgF2->SetNpx(1000);
+
+ delete gROOT->GetListOfFunctions()->FindObject("FitQLongCylSurf");
+ if (max <= 0.0) max = h->GetXaxis()->GetXmax();
+ TF1 *theory = new TF1("BorisTomasikCylSurfQLong",&AliHBTFits::QLongCylSurf,0.0,max,3);//par : aL, aPhi,aLambda
+ theory->SetNpx(1000);
+ theory->SetParameter(0,2.0);//L
+ theory->SetParameter(1,1.0);//Phi
+ theory->SetParameter(2,0.75);//Lambda
+
+ theory->FixParameter(0,1.81);//L
+ theory->FixParameter(1,0.864);//Phi
+ theory->FixParameter(2,0.773);//Lambda
+
+ theory->SetParName(0,"R_{long}");
+ theory->SetParName(1,"\\Phi");
+ theory->SetParName(2,"\\lambda");
+
+
+ h->Fit(theory,fopt,"E");
+
+}
+
+Double_t AliHBTFits::QLongCylSurf(Double_t *x, Double_t *par)
+{
+ //Qout
+ Double_t qlong = x[0];
+ Double_t aL = par[0];
+ Double_t aPhi = par[1];
+ Double_t aLambda = par[2];
+
+ static const Double_t j2=1.063*1.063;
+
+ fLongOnly = 1;
+ Double_t pars[5];
+ Double_t xs[5];
+ xs[0] = 0.0;
+ xs[1] = 0.0;
+ xs[2] = qlong;
+ pars[0] = aLambda;
+ pars[1] = 1.0;
+ pars[2] = 1.063;
+ pars[3] = aL;
+ pars[4] = aPhi;
+ fLongOnly = 0;
+ return XXX(&(xs[0]),&(pars[0]));
+
+ Double_t result=aLambda*TMath::Exp(-qlong*qlong*j2/0.0389273);
+
+
+
+// ::Info("QLongCylSurf","l*exp %f",result);
+ Double_t longpart;
+ if ( (qlong != 0.0) && (aL != 0.0) )
+ {
+ Double_t qlL = qlong*aL/(2.0*0.197);
+ Double_t sinqlL = TMath::Sin(qlL);
+ Double_t sin2qlL = sinqlL*sinqlL;
+ longpart = sin2qlL/(qlL*qlL);
+ result*= longpart;
+ }
+
+// ::Info("QLongCylSurf","longpart %f",result);
+
+ if (aPhi > 0.0)
+ {
+ Double_t oneOverPhi = 1./aPhi;
+
+ Double_t tmp1 = TMath::StruveL0(oneOverPhi)/TMath::BesselI0(oneOverPhi);
+ Double_t tmp = 1.0 - tmp1*tmp1;
+ result*=tmp;
+ ::Info("QLongCylSurf","oneOverPhi %f, tmp1 %f, tmp %f",oneOverPhi ,tmp1, tmp);
+ }
+// ::Info("QLongCylSurf","Phipart %f",result);
+
+ result+=1.;
+
+// ::Info("QLongCylSurf","Final: result=%f",result);
+ return result;
+
+}
+Double_t AliHBTFits::XXX(Double_t *x, Double_t *par)
+{
+
+ //QOutQSideQLongCylSurf Function
+ Double_t qout = x[0];
+ Double_t qside = x[1];
+ Double_t qlong = x[2];
+
+ Double_t aLambda = par[0];
+ Double_t aR = par[1];
+ Double_t aJ = par[2];
+ Double_t aL = par[3];
+ Double_t aPhi = par[4];
+
+ Double_t result=aLambda*TMath::Exp(-(qout*qout+qside*qside+qlong*qlong)*aJ*aJ/0.0389273);
+
+// ::Info("XXX","lambda*exp=%f",result);
+ Double_t longpart;
+ if ( (qlong != 0.0) && (aL != 0.0) )
+ {
+ Double_t qlL = qlong*aL/(2.0*0.197);
+ Double_t sinqlL = TMath::Sin(qlL);
+ Double_t sin2qlL = sinqlL*sinqlL;
+ longpart = sin2qlL/(qlL*qlL);
+ result*= longpart;
+ }
+
+// ::Info("XXX","longpart=%f",longpart);
+
+ if (aPhi > 0.0)
+ {
+ fgF1->SetParameter(0,aR);
+ fgF2->SetParameter(0,aR);
+
+ fgF1->SetParameter(1,aPhi);
+ fgF2->SetParameter(1,aPhi);
+
+ fgF1->SetParameter(2,qout);
+ fgF2->SetParameter(2,qout);
+
+ fgF1->SetParameter(3,qside);
+ fgF2->SetParameter(3,qside);
+
+ Double_t sp = fgF1->Integral(-TMath::Pi(),TMath::Pi(),(Double_t*)0x0,1.e-8);
+ Double_t cp = fgF2->Integral(-TMath::Pi(),TMath::Pi(),(Double_t*)0x0,1.e-8);
+
+ ::Info("QSideCylSurf","sp=%f, cp=%f",cp,sp);
+
+ result *= cp*cp + sp*sp;
+// ::Info("QSideCylSurf","2: result=%f",result);
+
+
+ Double_t oneOverPhi = 1./aPhi;
+ Double_t tmp = TMath::BesselI0(oneOverPhi)-TMath::StruveL0(oneOverPhi);
+ tmp = tmp*TMath::TwoPi();
+ tmp = tmp*tmp;
+ if (tmp == 0.0)
+ {
+ ::Error("QOutCylSurf","Div by 0");
+ return 1.0;
+ }
+// ::Info("XXX","itgrl/bessel=%f",(cp*cp + sp*sp)/tmp);
+
+ result=result/tmp;
+ }
+
+ result+=1.;
+// ::Info("QSideCylSurf","Final: result=%f",result);
+ return result;
+}