Classes for Corrections with Lisa algorithm added

diff --git a/HBTAN/AliHBTCorrectCorrelFctn.cxx b/HBTAN/AliHBTCorrectCorrelFctn.cxx
new file mode 100644 (file)
index 0000000..c50174c
--- /dev/null
+++ b/HBTAN/AliHBTCorrectCorrelFctn.cxx
@@ -0,0 +1,441 @@
+#include "AliHBTCorrectCorrelFctn.h"
+//____________________
+///////////////////////////////////////////////////////
+//                                                   //
+// AliHBTCorrectQInvCorrelFctn                       //
+//                                                   //
+// Class for calculating Q Invariant correlation     //
+// taking to the account resolution of the           //
+// detector and coulomb effects.                     //
+// Implemented on basis of STAR procedure            //
+// implemented and described by Michael A. Lisa      //
+//                                                   //
+// Piotr.Skowronski@cern.ch                          //
+// http://alisoft.cern.ch/people/skowron/analyzer    //
+//                                                   //
+///////////////////////////////////////////////////////
+
+//Parameters fit from pi+ pi+ resolution analysis for pair with qinv<50MeV
+// chi2/NFD = 97/157
+// FCN=98.0971 FROM MIGRAD    STATUS=CONVERGED     332 CALLS         333 TOTAL
+//   EDM=2.13364e-12    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   2.4 per cent
+//  EXT PARAMETER                                   STEP         FIRST
+//  NO.   NAME         VALUE          ERROR         SIZE      DERIVATIVE
+//   1   fThetaA     2.72546e-03   1.43905e-04  -0.00000e+00   5.54375e-02
+//   2   fThetaB     1.87116e-04   5.11862e-05   0.00000e+00   3.66500e-01
+//   3  fThetaAlpha -2.36868e+00   1.83230e-01   0.00000e+00  -7.01301e-05
+
+// FCN=120.603 FROM MIGRAD    STATUS=CONVERGED     117 CALLS         118 TOTAL
+//  EDM=2.5153e-12    STRATEGY= 1  ERROR MATRIX UNCERTAINTY   2.4 per cent
+//  EXT PARAMETER                                   STEP         FIRST
+//  NO.   NAME         VALUE          ERROR         SIZE      DERIVATIVE
+//   1   fPhiA       1.93913e-03   1.31059e-04  -0.00000e+00  -5.87280e-02
+//   2   fPhiA       2.48687e-04   5.41251e-05   0.00000e+00  -1.87361e-01
+//   3  fPhiAlpha   -2.22649e+00   1.44503e-01   0.00000e+00   8.97538e-06
+
+#include <TH1.h>
+#include <TH3.h>
+#include <TF1.h>
+#include <TRandom.h>
+ClassImp(AliHBTCorrectQInvCorrelFctn)
+
+AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name,const char* title):
+  AliHBTOnePairFctn1D(name,title),
+  fMeasCorrelFctn(0x0),
+  fMeasNumer(0x0),
+  fMeasDenom(0x0),
+  fSmearedNumer(0x0),
+  fSmearedDenom(0x0),
+  fDPtOverPtRMS(0.004),
+  fThetaA(2.72e-03),
+  fThetaB(1.87e-04),
+  fThetaAlpha(-2.4),
+  fPhiA(1.94e-03),
+  fPhiB(2.5e-04),
+  fPhiAlpha(-2.2),
+  fR2(0.0),
+  fLambda(0.0),
+  fRConvergenceTreshold(0.3),
+  fLambdaConvergenceTreshold(0.05)
+{
+
+}
+/******************************************************************/
+
+AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(TH1D* measqinv,const char* name,const char* title):
+  AliHBTOnePairFctn1D(name,title,
+                      measqinv->GetNbinsX(),
+         measqinv->GetXaxis()->GetXmax(),
+         measqinv->GetXaxis()->GetXmin()),
+  fMeasCorrelFctn(measqinv),
+  fMeasNumer(0x0),
+  fMeasDenom(0x0),
+  fSmearedNumer(0x0),
+  fSmearedDenom(0x0),
+  fDPtOverPtRMS(0.004),
+  fThetaA(2.72e-03),
+  fThetaB(1.87e-04),
+  fThetaAlpha(-2.4),
+  fPhiA(1.94e-03),
+  fPhiB(2.5e-04),
+  fPhiAlpha(-2.2),
+  fR2(0.0),
+  fLambda(0.0),
+  fRConvergenceTreshold(0.3),
+  fLambdaConvergenceTreshold(0.05)
+{
+  
+}
+/******************************************************************/
+
+AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name, const char* title,
+                   Int_t nbins, Float_t maxXval, Float_t minXval):
+  AliHBTOnePairFctn1D(name,title,nbins,maxXval,minXval),
+  fMeasCorrelFctn(0x0),
+  fMeasNumer(0x0),
+  fMeasDenom(0x0),
+  fSmearedNumer(0x0),
+  fSmearedDenom(0x0),
+  fDPtOverPtRMS(0.004),
+  fThetaA(2.72e-03),
+  fThetaB(1.87e-04),
+  fThetaAlpha(-2.4),
+  fPhiA(1.94e-03),
+  fPhiB(2.5e-04),
+  fPhiAlpha(-2.2),
+  fR2(0.0),
+  fLambda(0.0),
+  fRConvergenceTreshold(0.3),
+  fLambdaConvergenceTreshold(0.05)
+{
+}                  
+/******************************************************************/
+
+AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
+{
+  delete fMeasCorrelFctn;
+  delete fSmearedNumer;
+  delete fSmearedDenom;
+  delete fMeasNumer;
+  delete fMeasDenom;
+}
+/******************************************************************/
+
+void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min)
+{
+  AliHBTFunction1D::BuildHistos(nbins,max,min);
+  
+  TString numstr = fName + " Smeared Numerator";  //title and name of the numerator histogram
+  TString denstr = fName + " Smeared Denominator";//title and name of the denominator histogram
+  
+  fSmearedNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
+  fSmearedDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
+  fSmearedNumer->Sumw2();
+  fSmearedDenom->Sumw2();
+  
+  if (fMeasCorrelFctn == 0x0)
+   { 
+     numstr = fName + " Measured Numerator";  //title and name of the numerator histogram
+     denstr = fName + " Measured Denominator";//title and name of the denominator histogram
+    
+     fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
+     fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
+     fMeasNumer->Sumw2();
+     fMeasDenom->Sumw2();
+   }
+}
+/******************************************************************/
+
+void AliHBTCorrectQInvCorrelFctn::Init()
+{
+  AliHBTOnePairFctn1D::Init();
+  Info("Init","");
+  fSmearedNumer->Reset();
+  fSmearedDenom->Reset();
+  if (fMeasNumer) fMeasNumer->Reset();
+  if (fMeasDenom) fMeasDenom->Reset();
+  fFittedR = -1.0;
+  fFittedLambda = 0.0;
+}
+/******************************************************************/
+
+void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair)
+{
+  if (fMeasNumer == 0x0) return;
+  pair = CheckPair(pair);
+  if( pair == 0x0) return;
+  fMeasNumer->Fill(pair->GetQInv());
+}
+/******************************************************************/
+
+void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair)
+{
+//Process different events 
+  static AliHBTParticle part1, part2;
+  static AliHBTPair smearedpair(&part1,&part2);
+  
+  pair = CheckPair(pair);
+  if( pair == 0x0) return;
+  Double_t cc = GetCoulombCorrection(pair);
+  Double_t qinv = pair->GetQInv();
+
+  //measured histogram -> if we are interested 
+  //only if fMeasCorrelFctn is not specified by user
+  if (fMeasDenom) fMeasDenom->Fill(qinv,cc);
+
+  Smear(pair,smearedpair);
+  Double_t modelqinv = GetModelValue(qinv);  
+  //Ideal histogram
+  fNumerator->Fill(qinv,modelqinv*cc);
+  fDenominator->Fill(qinv,cc);
+  
+  //Smeared histogram
+  Double_t smearedqinv = smearedpair.GetQInv();
+  fSmearedNumer->Fill(smearedqinv,modelqinv);
+  Double_t smearedcc = GetCoulombCorrection(&smearedpair);
+  fSmearedDenom->Fill(smearedqinv,smearedcc);
+  
+}
+/******************************************************************/
+void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTPair* pair,AliHBTPair& smeared)
+{
+//Smears pair
+  Smear(pair->Particle1(),smeared.Particle1());
+  Smear(pair->Particle2(),smeared.Particle2());
+  smeared.Changed();
+}
+/******************************************************************/
+
+void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTParticle* part, AliHBTParticle* smeared)
+{
+  Double_t sin2theta = TMath::Sin(part->Theta());
+  sin2theta = sin2theta*sin2theta;
+  Double_t pt = part->Pt();
+
+  double DpT_div_pT = gRandom->Gaus(0.0,fDPtOverPtRMS);
+  double Dphi = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fPhiAlpha));
+  double Dtheta = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fThetaAlpha));
+  
+  Double_t smearedPx = part->Px()*(1.0+DpT_div_pT) - part->Py()*Dphi;
+//  fourmom.setX(px*(1.0+DpT_div_pT) - py*Dphi);
+  Double_t smearedPy = part->Py()*(1.0+DpT_div_pT) - part->Px()*Dphi;
+//  fourmom.setY(py*(1.0+DpT_div_pT) + px*Dphi);
+  Double_t smearedPz = part->Pz()*(1.0+DpT_div_pT) - pt*Dtheta/sin2theta;
+//  fourmom.setZ(pz*(1.0+DpT_div_pT) - pT*Dtheta/sin2theta);
+  
+  Double_t mass2 = part->GetMass()*part->GetMass();
+  Double_t E = mass2 + smearedPx*smearedPx + 
+                       smearedPy*smearedPy + 
+                       smearedPz*smearedPz;
+          
+  smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(E));
+}
+/******************************************************************/
+
+void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
+{
+  fLambda = lambda;
+  fR2 = r*r;
+}
+/******************************************************************/
+void AliHBTCorrectQInvCorrelFctn::MakeMeasCF()
+{
+  //makes measured correlation function
+  delete fMeasCorrelFctn;
+  fMeasCorrelFctn = 0x0;
+  
+  if (fMeasNumer&&fMeasDenom)
+   {
+     Double_t measscale = Scale(fMeasNumer,fMeasDenom);
+     if (measscale == 0.0)
+      {
+        Error("GetResult","GetRatio for measured CF returned 0.0");
+        return;
+      }
+     TString str = fName + "measured ratio";
+     fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data());
+     fMeasCorrelFctn->SetTitle(str.Data());
+     fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale);
+   }
+}
+
+TH1* AliHBTCorrectQInvCorrelFctn::GetResult()
+{
+  //In case we don't have yet Measured Correlation Function
+  //Try to get it
+  //result is
+  //        N[meas]   N[ideal]/D[ideal]
+  //C(Q) =  ------- * -----------------
+  //        D[meas]   N[smear]/D[smear]
+ 
+  TString str;
+  if (fMeasCorrelFctn == 0x0) MakeMeasCF();
+  
+  if (fMeasCorrelFctn == 0x0)
+   {
+     Error("GetResult",  
+           "Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null");
+     return 0x0;
+   }
+  
+  TH1D* ideal = (TH1D*)GetRatio(Scale());
+  if (ideal == 0x0)
+   {
+     Error("GetResult","Ratio of ideal histograms is null");
+     return 0x0;
+   }
+  str = fName + " smeared ratio";
+  TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
+  smearedCF->SetTitle(str.Data());
+  Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
+  smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
+  
+  str = fName + " product meas ideal CF";
+  TH1D* measideal = (TH1D*)ideal->Clone(str.Data());
+  measideal->Multiply(ideal,fMeasCorrelFctn);
+  
+  str = fName + " Corrected Result";
+  TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data());
+  result->SetTitle(str.Data());
+  
+  Double_t resultscale = Scale(measideal,smearedCF);
+  result->Divide(measideal,smearedCF,resultscale);
+  return result;
+  
+}
+/******************************************************************/
+
+void AliHBTCorrectQInvCorrelFctn::Fit()
+{
+//fits resuting histogram with function 1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))
+//where [0] is lambda
+//      [1] is radius
+//   0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI
+
+  Info("Fit","Before fFittedLambda = %f",fFittedLambda);
+  Info("Fit","Before fFittedR = %f",fFittedR);
+  TH1D* result = (TH1D*)GetResult();
+  if (result == 0x0)
+   {
+     Error("Fit","Can not get result");
+     return;
+   }
+  TF1* fitfctn = new TF1("fitfctn","1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))");
+  fitfctn->SetParameter(0,1.0);
+  fitfctn->SetParameter(1,6.0);
+  Float_t max = result->GetXaxis()->GetXmax();
+  Info("Fit","Max is %f",max);
+  result->Fit(fitfctn,"0","",0.008,max);
+  fFittedLambda = fitfctn->GetParameter(0);
+  fFittedR = fitfctn->GetParameter(1);
+  Info("Fit","After fFittedLambda = %f",fFittedLambda);
+  Info("Fit","After fFittedR = %f",fFittedR);
+  delete fitfctn;
+  delete result;
+}
+/******************************************************************/
+
+Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged()
+{
+  //check if fitting was performed
+  if (fFittedR <= 0.0)
+   {
+     Fit();//if not do fit first
+     if (fFittedR <= 0.0)
+      {
+        Error("IsConverged","Fitting failed");
+        return kFALSE;
+      }
+   }
+
+  Double_t guessedR = TMath::Sqrt(fR2);
+  
+  Info("IsConverged","Fitted   lambda            : %8f Fitted  Radius             : %8f",fFittedLambda,fFittedR);
+  Info("IsConverged","Guessed  lambda            : %8f Guessed Radius             : %8f",fLambda,guessedR);
+  Info("IsConverged","Demanded lambda convergence: %8f Demanded Radius convergence: %8f",
+       fLambdaConvergenceTreshold,fRConvergenceTreshold);
+  
+  if ( (TMath::Abs(fLambda-fFittedLambda)<fLambdaConvergenceTreshold) && 
+       (TMath::Abs(fFittedR-guessedR)<fRConvergenceTreshold) )
+    {
+      Info("IsConverged","Cnvergence reached");
+      return kTRUE;
+    }
+  else
+   {
+      Info("IsConverged","Cnvergence NOT reached");
+      return kFALSE;
+   }
+}
+/******************************************************************/
+
+Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius()
+{
+  if (fFittedR <= 0.0) Fit();
+  return fFittedR;
+}
+/******************************************************************/
+
+Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
+{
+  if (fFittedR <= 0.0) Fit();
+  return fFittedLambda;
+}
+/******************************************************************/
+
+void AliHBTCorrectQInvCorrelFctn::WriteAll()
+{
+  Write();
+  if (fMeasCorrelFctn) fMeasCorrelFctn->Write();
+  if (fMeasNumer ) fMeasNumer->Write();
+  if (fMeasDenom ) fMeasDenom->Write();
+  if (fSmearedNumer) fSmearedNumer->Write();
+  if (fSmearedDenom) fSmearedDenom->Write();
+  if (fSmearedNumer && fSmearedDenom)
+   {
+    TString str = fName + " smeared ratio";
+    TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
+    smearedCF->SetTitle(str.Data());
+    Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
+    smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
+   }
+}
+
+/******************************************************************/
+/******************************************************************/
+/******************************************************************/
+
+//____________________
+///////////////////////////////////////////////////////
+//                                                   //
+// AliHBTCorrectQ3DCorrelFctn                        //
+//                                                   //
+// Class for calculating Q Invariant correlation     //
+// taking to the account resolution of the           //
+// detector and coulomb effects.                     //
+//                                                   //
+///////////////////////////////////////////////////////
+
+
+AliHBTCorrectQ3DCorrelFctn::AliHBTCorrectQ3DCorrelFctn(const char* name, const char* title):
+ AliHBTOnePairFctn3D(name,title),
+ fMeasCorrelFctn(0x0),
+ fSmearedNumer(0x0),
+ fSmearedDenom(0x0),
+ fMeasNumer(0x0),
+ fMeasDenom(0x0)
+{
+
+}
+/******************************************************************/
+
+AliHBTCorrectQ3DCorrelFctn::~AliHBTCorrectQ3DCorrelFctn()
+{
+ delete fMeasCorrelFctn;
+ delete fSmearedNumer;
+ delete fSmearedDenom;
+ delete fMeasNumer;
+ delete fMeasDenom;
+}
+
+/******************************************************************/

diff --git a/HBTAN/AliHBTCorrectCorrelFctn.h b/HBTAN/AliHBTCorrectCorrelFctn.h
new file mode 100644 (file)
index 0000000..01f0eb7
--- /dev/null
+++ b/HBTAN/AliHBTCorrectCorrelFctn.h
@@ -0,0 +1,132 @@
+#ifndef ALIHBTCORRECTCORRELFCTN_H
+#define ALIHBTCORRECTCORRELFCTN_H
+//____________________
+///////////////////////////////////////////////////////
+//                                                   //
+// AliHBTCorrectQInvCorrelFctn                       //
+//                                                   //
+// Class for calculating Q Invariant correlation     //
+// taking to the account resolution of the           //
+// detector and coulomb effects.                     //
+//                                                   //
+///////////////////////////////////////////////////////
+
+#include "AliHBTFunction.h"
+
+  class AliHBTCorrectQInvCorrelFctn: public AliHBTOnePairFctn1D
+{
+  public:
+    AliHBTCorrectQInvCorrelFctn(const char* name = "qinvcorrectedCF", 
+                                const char* title= "Corrected Q_{inv} Correlation Fonction");
+
+    AliHBTCorrectQInvCorrelFctn(const char* name, const char* title,
+                   Int_t nbins, Float_t maxXval, Float_t minXval);
+
+    AliHBTCorrectQInvCorrelFctn(TH1D* measqinv, 
+                                const char* name = "qinvcorrectedCF", 
+                                const char* title= "Corrected Q_{inv} Correlation Fonction");
+    virtual ~AliHBTCorrectQInvCorrelFctn();
+    void SetInitialValues(Double_t lambda, Double_t r);
+    void Init();
+    void ProcessSameEventParticles(AliHBTPair* pair);//process particles from same event (real pair)
+    void ProcessDiffEventParticles(AliHBTPair* pair);//process particles coming from different events (mixed pairs)
+    void SetMeasuredHistogram(TH1D* meas){fMeasCorrelFctn = meas;}
+    TH1* GetResult();//returns the result histogram
+    Double_t GetRadius()const{ return TMath::Sqrt(fR2);}//returns assumed radius
+    Double_t GetLambda()const{ return fLambda;}//retutrns assumed intercept parameter
+    void     SetRadiusConvergenceTreshold(Double_t ct){fRConvergenceTreshold=ct;}//if fitted and assumed R us different less then that number con
+    void     SetLambdaConvergenceTreshold(Double_t ct){fLambdaConvergenceTreshold=ct;}
+    Bool_t   IsConverged();
+    void     Fit();
+    Double_t GetFittedRadius();
+    Double_t GetFittedLambda();
+    void     WriteAll();
+    void     SetMeasNum(TH1D* measnum){fMeasNumer = measnum;}
+    void     SetMeasDen(TH1D* h){fMeasDenom = h;}
+    void     MakeMeasCF();
+  protected:
+    virtual void BuildHistos(Int_t nbins, Float_t max, Float_t min);
+    Double_t GetCoulombCorrection(AliHBTPair* pair){return 1.0;}
+    Double_t GetValue(AliHBTPair * pair){return pair->GetQInv();}
+    void Smear(AliHBTPair* pair,AliHBTPair& smeared);
+    void Smear(AliHBTParticle* part, AliHBTParticle* smeared);
+    Double_t GetModelValue(Double_t qinv);
+
+    //Our ideal numerator 
+    TH1D* fMeasCorrelFctn;
+    TH1D* fMeasNumer;
+    TH1D* fMeasDenom;  
+    
+    TH1D* fSmearedNumer; //! Numerator of smeard q
+    TH1D* fSmearedDenom; //! Denominator of smeard q
+    
+    //Parameters of Pt RMS
+    //linear dependence dPt/Pt from Pt itself 
+    Float_t fDPtOverPtRMS; //RMS of dPt/Pt 
+    
+    //We assume that RMS of Theta and Phisangle depends on Pt Like A+B*(Pt)^Alpha
+    //Idea copied from Star HBT Maker (Fabrice Retiere)
+    //Parameters comes from Monte Carlo Resolution Analysis
+
+    Float_t fThetaA; //"A" parameter of theta RMS dependence
+    Float_t fThetaB; //"B" parameter of theta RMS dependence
+    Float_t fThetaAlpha; //"Alpha" parameter (power) of theta RMS dependence
+
+    Float_t fPhiA;//"A" parameter of phi RMS dependence
+    Float_t fPhiB;//"B" parameter of phi RMS dependence
+    Float_t fPhiAlpha;//"Alpha" parameter (power) of phi RMS dependence
+    
+    Double_t fR2;//square of radius
+    Double_t fLambda;//Interception parameter
+
+    Double_t fFittedR;//fitted radius
+    Double_t fFittedLambda;//fitted Interception parameter
+        
+    Float_t  fRConvergenceTreshold;
+    Float_t  fLambdaConvergenceTreshold;
+  private:
+  public:
+    ClassDef(AliHBTCorrectQInvCorrelFctn,1)
+};
+
+inline
+Double_t AliHBTCorrectQInvCorrelFctn::GetModelValue(Double_t qinv)
+{
+  //factor 0.038936366329 conected with units change GeV<->SI
+  return 1.0 + fLambda*TMath::Exp(-fR2*qinv*qinv/0.038936366329);
+}
+
+//____________________
+///////////////////////////////////////////////////////
+//                                                   //
+// AliHBTCorrectQ3DCorrelFctn                        //
+//                                                   //
+// Class for calculating Q Invariant correlation     //
+// taking to the account resolution of the           //
+// detector and coulomb effects.                     //
+//                                                   //
+///////////////////////////////////////////////////////
+
+class AliHBTCorrectQ3DCorrelFctn: public AliHBTOnePairFctn3D
+{
+  public:
+   AliHBTCorrectQ3DCorrelFctn(const char* name = "qinvcorrectedCF", 
+                               const char* title= "Corrected Q_{inv} Correlation Fonction");
+   virtual ~AliHBTCorrectQ3DCorrelFctn();
+   
+  protected:
+    TH3D* fMeasCorrelFctn;
+    
+    TH3D* fSmearedNumer; //! Numerator of smeard q
+    TH3D* fSmearedDenom; //! Denominator of smeard q
+    TH3D* fMeasNumer;  //! Numerator of ideal q calculated on basis of model equation
+    TH3D* fMeasDenom;  //! Denominator of ideal q calculated on basis of model equation
+
+    
+  private:
+  
+  public:
+    ClassDef(AliHBTCorrectQ3DCorrelFctn,1)
+};
+
+#endif

diff --git a/HBTAN/HBTAnalysisLinkDef.h b/HBTAN/HBTAnalysisLinkDef.h
index 4c8c7c3fc95f0e0fc0ff2ac3f29c064733e34718..c0efe59b2d8c4e53c36f69e17388157e8914e9ef 100644 (file)
--- a/HBTAN/HBTAnalysisLinkDef.h
+++ b/HBTAN/HBTAnalysisLinkDef.h
@@ -85,6 +85,9 @@
 #pragma link C++ class AliHBTInvMassCorrelFctn+;
 #pragma link C++ class AliHBTCorrFitFctn+;
 
+#pragma link C++ class AliHBTCorrectQInvCorrelFctn+;
+#pragma link C++ class AliHBTCorrectQ3DCorrelFctn+;
+
 #pragma link C++ class AliHBTKtResolVsQInvFctn+;
 
 #pragma link C++ class AliHBTQInvResolVsQInvFctn+;

diff --git a/HBTAN/hbtcorrections.C b/HBTAN/hbtcorrections.C
new file mode 100644 (file)
index 0000000..bd9ebf8
--- /dev/null
+++ b/HBTAN/hbtcorrections.C
@@ -0,0 +1,81 @@
+void hbtcorrections(Int_t first = -1,Int_t last = -1)
+{  
+  
+  const char* basedir=".";
+  const char* serie="";
+  const char* field = "";
+
+  AliHBTReader* reader = new AliHBTReaderInternal("ESD.old.root");
+  TObjArray* dirs=0;
+  if ( (first >= 0) && (last>=0) && ( (last-first)>=0 ) )
+   {//read from many dirs dirs
+     char buff[50];
+     dirs = new TObjArray(last-first+1);
+     for (Int_t i = first; i<=last; i++)
+      {
+        sprintf(buff,"%s/%s/%s/%d",basedir,field,serie,i);
+        TObjString *odir= new TObjString(buff);
+        dirs->Add(odir);
+      }
+    }
+  reader->SetDirs(dirs);
+
+  AliHBTParticleCut* readerpartcut= new AliHBTParticleCut();
+  readerpartcut->SetPtRange(0.0,1.0);
+  readerpartcut->SetPID(kPiPlus);
+  reader->AddParticleCut(readerpartcut);//read this particle type with this cut
+  
+  AliHBTAnalysis* analysis = new AliHBTAnalysis();
+  analysis->SetReader(reader);
+  analysis->SetDisplayInfo(100000);
+  AliHBTPairCut *paircut = new AliHBTPairCut();
+  paircut->SetQInvRange(0.0,0.15);
+  analysis->SetGlobalPairCut(paircut);
+  
+
+
+  AliHBTCorrectQInvCorrelFctn* correctqinvCF = new AliHBTCorrectQInvCorrelFctn();
+  
+  analysis->AddTrackFunction(correctqinvCF);
+
+  TFile* outf = TFile::Open("hbtcorrected.root","recreate");
+  if (outf == 0x0)
+   {
+     cout<<"\n\nERROR: can not open file"<<endl;
+     return;
+   }
+  Int_t iter = 0;
+  TString dirname;
+  TString ietration("Iteration");
+  correctqinvCF->SetInitialValues(0.8,12.0);
+  correctqinvCF->SetRadiusConvergenceTreshold(0.005);
+  correctqinvCF->SetLambdaConvergenceTreshold(0.01);
+  while (1)
+   {
+    iter++;
+    dirname = ietration+iter;
+    TDirectory* dir = outf->mkdir(dirname,"results after "+dirname);
+    outf->Write();
+    if (dir) dir->cd();
+    else 
+     {
+       cout<<"\n\nERROR: can not make an directory in file named"<<dirname<<endl;
+       return;
+     }
+
+    analysis->Process("Tracks");
+
+/*****************************************/
+    dir->cd();
+    correctqinvCF->WriteAll();
+/*****************************************/
+    if (correctqinvCF->IsConverged()) break;
+    else
+     {
+       correctqinvCF->SetInitialValues(correctqinvCF->GetFittedLambda(),correctqinvCF->GetFittedRadius());
+     }
+   }
+  outf->cd();
+  analysis->Write();
+  delete outf;
+}

diff --git a/HBTAN/libHBTAN.pkg b/HBTAN/libHBTAN.pkg
index 219216a50d00e75f7c1c107108d65377f1ac0974..6dda49eb9809e04e3b79da59a554dc80d362d5c4 100644 (file)
--- a/HBTAN/libHBTAN.pkg
+++ b/HBTAN/libHBTAN.pkg
@@ -13,7 +13,7 @@ AliHBTMonDistributionFctns.cxx  AliHBTMonResolutionFctns.cxx \
 AliHBTLLWeights.cxx             AliHBTLLWeightFctn.cxx \
 AliHBTLLWeightsPID.cxx          AliHBTLLWeightTheorFctn.cxx \
 AliHBTPositionRandomizer.cxx AliHBTEventBuffer.cxx \
-AliHBTCorrFitFctn.cxx
+AliHBTCorrFitFctn.cxx AliHBTCorrectCorrelFctn.cxx
 
 FSRCS   = fsiini.F  fsiw.F  led_bldata.F  ltran12.F