Default constructor added to ITS separation cut

diff --git a/HBTAN/AliHBTFits.cxx b/HBTAN/AliHBTFits.cxx
index 31f727571a248f4d995be259ff0cc7327731b59b..9cd1e30ecb2609a9cbe0bf3812a90af16c1ddfe0 100644 (file)
--- a/HBTAN/AliHBTFits.cxx
+++ b/HBTAN/AliHBTFits.cxx
@@ -12,6 +12,7 @@
 // Piotr.Skowronski@cern.ch
 //
 ///////////////////////////////////////////////////////////////////////////////////
+#include "AliHBTAnalysis.h"
 
 #include <TROOT.h>
 #include <TH1.h>
@@ -28,6 +29,24 @@ ClassImp(AliHBTFits)
 
 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()
@@ -84,8 +103,84 @@ void AliHBTFits::FitQOutQSideQLongCylSurf(const TString& hname,Option_t* fopt, F
  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;
 
 }
 
@@ -93,9 +188,38 @@ void AliHBTFits::FitQOutQSideQLongCylSurf(const TString& hname,Option_t* fopt, F
 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];
@@ -103,69 +227,176 @@ Double_t AliHBTFits::QOutQSideQLongCylSurf(Double_t *x, Double_t *par)
   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;
 }
 /******************************************************************************************/
@@ -232,7 +463,7 @@ Double_t AliHBTFits::QOutQSideCylSurf(Double_t *x, Double_t *par)
   
   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;
    }
@@ -261,8 +492,10 @@ Double_t AliHBTFits::QOutQSideCylSurf(Double_t *x, Double_t *par)
     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");
@@ -302,9 +535,15 @@ void AliHBTFits::FitQOutCylSurf(const TString& hname, Option_t* fopt, Float_t ma
  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");
@@ -323,8 +562,8 @@ Double_t AliHBTFits::QOutCylSurf(Double_t *x, Double_t *par)
   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);
   
@@ -347,8 +586,10 @@ Double_t AliHBTFits::QOutCylSurf(Double_t *x, Double_t *par)
     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");
@@ -389,6 +630,11 @@ void AliHBTFits::FitQSideCylSurf(const TString& hname, Option_t* fopt, Float_t m
  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");
@@ -408,7 +654,8 @@ Double_t AliHBTFits::QSideCylSurf(Double_t *x, Double_t *par)
 
 //  ::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);
@@ -425,19 +672,188 @@ Double_t AliHBTFits::QSideCylSurf(Double_t *x, Double_t *par)
 
     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;
 
+}

diff --git a/HBTAN/AliHBTFits.h b/HBTAN/AliHBTFits.h
index 479c63ed811e0345092713628e2d7f53334bae67..c9a15760542af4b8df4813d8bd684a7943bf1973 100644 (file)
--- a/HBTAN/AliHBTFits.h
+++ b/HBTAN/AliHBTFits.h
@@ -25,6 +25,7 @@ class AliHBTFits: public TObject
     virtual ~AliHBTFits();
     static void FitQOutCylSurf (const TString& hname, Option_t* fopt = "R", Float_t max = 0.0);
     static void FitQSideCylSurf(const TString& hname, Option_t* fopt = "R", Float_t max = 0.0);
+    static void FitQLongCylSurf(const TString& hname, Option_t* fopt = "R", Float_t max = 0.0);
 
     static void FitQOutQSideQLongCylSurf(const TString& hname,Option_t* fopt = "R",
                  Float_t xmax = 0,Float_t ymax = 0, Float_t zmax = 0);
@@ -34,6 +35,7 @@ class AliHBTFits: public TObject
     
     static Double_t QOutCylSurf(Double_t *x, Double_t *par);
     static Double_t QSideCylSurf(Double_t *x, Double_t *par);
+    static Double_t QLongCylSurf(Double_t *x, Double_t *par);
     
     static Double_t QOutQSideQLongCylSurf(Double_t *x, Double_t *par);
     static Double_t QOutQSideCylSurf(Double_t *x, Double_t *par);
@@ -41,11 +43,32 @@ class AliHBTFits: public TObject
     
 //    static double Qromo(double (*fun)(double), double a, double b,double(*choose)())
 //    static double Midexp(double (*fun)(double), double aa, double bb, int n)
-    
+   static Double_t XXX(Double_t *x, Double_t *par);
+   static Bool_t fLongOnly;
+   
   protected:
   private:
     static TF1*     fgF1; //functions
     static TF1*     fgF2; //functions
+    
+    static Double_t *fXbins;
+    static Double_t *fYbins;
+    static Double_t *fZbins;
+    
+    static Double_t *fExpT;
+    static Double_t *fLongT;
+    static Double_t *fIntegrT;
+    
+    static Int_t fNX;
+    static Int_t fNY;
+    static Int_t fNZ;
+    static Int_t fNXY;
+    static Int_t fNT;
+    
+    static Int_t fBinX;
+    static Int_t fBinY;
+    static Int_t fBinZ;
+    
     ClassDef(AliHBTFits,1)
 };
 #endif

diff --git a/HBTAN/AliHBTPairCut.cxx b/HBTAN/AliHBTPairCut.cxx
index 1acdc8e098233ab8578e30416bfb281710af5b5d..aa7ea4af381bf11f136e3b5d78e2b8e6eafa9976 100644 (file)
--- a/HBTAN/AliHBTPairCut.cxx
+++ b/HBTAN/AliHBTPairCut.cxx
@@ -27,6 +27,10 @@ AliHBTPairCut::AliHBTPairCut():
   fSecondPartCut= new AliHBTEmptyParticleCut(); //empty cuts
     
   fCuts = new AliHbtBasePairCut*[fgkMaxCuts];
+  for (Int_t i = 0;i<fNCuts;i++)
+   {
+     fCuts[i] = 0x0;
+   }
 }
 /**********************************************************/
 
@@ -262,6 +266,7 @@ void AliHBTPairCut::SetITSSeparation(Int_t layer, Double_t drphi, Double_t dz)
       }
    }
   fCuts[fNCuts++] = new AliHBTITSSeparationCut(layer,drphi,dz);
+//  Info("SetITSSeparation","Added %d at address %#x",fNCuts-1,fCuts[fNCuts-1]);
 }
 /**********************************************************/
 

diff --git a/HBTAN/AliHBTPairCut.h b/HBTAN/AliHBTPairCut.h
index 48a64ac6b38dd6c5ab86486b5b1ebbbfd2d7324c..2279bc6ae0d4a8eb23564468df92f12897a4e863 100644 (file)
--- a/HBTAN/AliHBTPairCut.h
+++ b/HBTAN/AliHBTPairCut.h
@@ -305,7 +305,7 @@ class AliHBTITSSeparationCut: public AliHbtBasePairCut
 {
 //Anti merging cut for the first layer of pixels
  public:
-  AliHBTITSSeparationCut(Int_t layer, Double_t deltarphi = 0.01, Double_t deltaz = 0.08):
+  AliHBTITSSeparationCut(Int_t layer = 0, Double_t deltarphi = 0.01, Double_t deltaz = 0.08):
     AliHbtBasePairCut(deltarphi,deltaz,kHbtPairCutPropPixelSepar),fLayer(layer){}
   virtual ~AliHBTITSSeparationCut(){}
   Bool_t   Pass(AliHBTPair* pair) const;

diff --git a/HBTAN/AliHBTTrackPoints.cxx b/HBTAN/AliHBTTrackPoints.cxx
index b7795cdbfa68166cbe884d660a1611b68394325d..6eae6222e4d75ca90bd7fb6b4472cb101b26eee9 100644 (file)
--- a/HBTAN/AliHBTTrackPoints.cxx
+++ b/HBTAN/AliHBTTrackPoints.cxx
@@ -45,9 +45,8 @@ AliHBTTrackPoints::AliHBTTrackPoints(AliHBTTrackPoints::ETypes type, AliESDtrack
   switch (type)
    {
      case kITS:
-       //Up to now we keep only point in pixels
        //Used only in non-id analysis
-       fN = 1;
+       fN = 6;
        fX = new Float_t[fN];
        fY = new Float_t[fN];
        fZ = new Float_t[fN];
@@ -278,13 +277,16 @@ void AliHBTTrackPoints::MakeITSPoints(AliESDtrack* track)
 // z=R*Pz/Pt
  AliITStrackV2 itstrack(*track,kTRUE);
  Double_t x,y,z;
- itstrack.GetGlobalXYZat(4.0,x,y,z);
- fX[0] = x;
- fY[0] = y;
- fZ[0] = z;
- 
-//  Info("MakeITSPoints","X %f Y %f Z %f R asked %f R obtained %f",
-//             fX[0],fY[0],fZ[0],4.0,TMath::Hypot(fX[0],fY[0]));
+ static const Double_t r[6] = {4.0, 7.0, 14.9, 23.8, 39.1, 43.6};
+ for (Int_t i = 0; i < 6; i++)
+  {
+    itstrack.GetGlobalXYZat(r[i],x,y,z);
+    fX[i] = x;
+    fY[i] = y;
+    fZ[i] = z;
+//    Info("MakeITSPoints","X %f Y %f Z %f R asked %f R obtained %f",
+//             fX[i],fY[i],fZ[i],r[i],TMath::Hypot(fX[i],fY[i]));
+  }   
  
 }
 

diff --git a/HBTAN/AliHBTWeightFctn.cxx b/HBTAN/AliHBTWeightFctn.cxx
index ea10f41a0b7e3d00f39fe9035084b3c4ab69890c..57e6e55a7fe3806873c3c4761cdd3f3b7fc1ebb6 100644 (file)
--- a/HBTAN/AliHBTWeightFctn.cxx
+++ b/HBTAN/AliHBTWeightFctn.cxx
@@ -496,14 +496,24 @@ void AliHBTWeightQOutSQideQLongFctn::ProcessSameEventParticles(AliHBTPair* track
     Double_t side = TMath::Abs(trackpair->GetQSideCMSLC());
     Double_t lon = TMath::Abs(trackpair->GetQLongCMSLC());
     
-    if (out < 0.005)
+    if (out < 0.001)
      if (side < 0.005)
        if (lon < 0.005)
         {
-          trackpair->Particle1()->Print();
-          trackpair->Particle2()->Print();
+          Info("","================================================");
+          Info("","BumBumBum");
           Info("","Delta Theta %f, Delta Phi %f",
               trackpair->GetDeltaTheta(),trackpair->GetDeltaPhi());
+          Info("\n","Track1");
+          trackpair->Particle1()->Print();
+          Info("\n","Track2");
+          trackpair->Particle2()->Print();
+          Info("\n","Particle1");
+          partpair->Particle1()->Print();
+          Info("\n","Particle2");
+          partpair->Particle2()->Print();
+          fflush(0);
+
         }
     
     fNumerator->Fill(out,side,lon,weight);//here we fill in q's corresponding to track pair 

diff --git a/HBTAN/AliHBTWeightTheorFctn.cxx b/HBTAN/AliHBTWeightTheorFctn.cxx
index 69338190bb625cac46128bb2e1fc4f601899f8b7..b838434231867248471fa2c3c47a0c19f4ca8c8a 100644 (file)
--- a/HBTAN/AliHBTWeightTheorFctn.cxx
+++ b/HBTAN/AliHBTWeightTheorFctn.cxx
@@ -185,6 +185,21 @@ void AliHBTWeightTheorOSLFctn::ProcessSameEventParticles(AliHBTPair* partpair)
   Double_t out = TMath::Abs(partpair->GetQOutCMSLC());
   Double_t side = TMath::Abs(partpair->GetQSideCMSLC());
   Double_t lon = TMath::Abs(partpair->GetQLongCMSLC());
+
+/*  
+  if (out < 0.01)
+    if (side < 0.01)
+      if (lon < 0.01)
+       {
+         Info("TheorOSL Num","================================================================");
+         Info("TheorOSL Num","o:%f, s:%f, l:%f, w%f",out,side,lon,weight);
+         Info("TheorOSL Num","First");
+         partpair->Particle1()->Print();
+         Info("TheorOSL Num","Second");
+         partpair->Particle2()->Print();
+         fflush(0);
+       }
+*/ 
   fNumerator->Fill(out,side,lon,weight);
 }
 /*************************************************************/