[u/mrichter/AliRoot.git] / HBTAN / AliHBTCorrectQInvCorrelFctn.cxx

#include "AliHBTCorrectQInvCorrelFctn.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://aliweb.cern.ch/people/skowron/analyzer    //
//                                                   //
///////////////////////////////////////////////////////

// dPt/Pt
// root [19] rms->Fit(f,"","",0.1,2)
//  FCN=7.0017 FROM MIGRAD    STATUS=CONVERGED     126 CALLS         127 TOTAL
//                      EDM=2.28804e-15    STRATEGY= 1      ERROR MATRIX ACCURATE
//    EXT PARAMETER                                   STEP         FIRST
//  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE
//    1  p0           5.78220e-03   3.14576e-05   4.97822e-09  -1.90059e-05
//    2  p1           3.98063e-05   1.61877e-06   1.04380e-10   1.91454e-04
//    3  p2          -2.78008e+00   2.13581e-02   1.66031e-06   3.16574e-06
//    4  p3           5.07594e-04   4.79619e-05   1.29793e-08  -2.29242e-05


// Phi
// root [17] rms->Fit(f,"","",0.15,2.5)
// Warning in <TH1D::Fit>: Abnormal termination of minimization.
// FCN=33.4898 FROM MIGRAD    STATUS=FAILED         91 CALLS          92 TOTAL
//                     EDM=1.19154e-15    STRATEGY= 1      ERR MATRIX APPROXIMATE
//    EXT PARAMETER                APPROXIMATE        STEP         FIRST
//   NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE
//    1  p0           5.87693e-04   5.04254e-06   2.49187e-09   5.84546e-04
//    2  p1           2.16488e-06   3.68880e-07   6.41507e-11  -7.36564e-02
//    3  p2          -3.10218e+00   1.01695e-01   2.00177e-05   7.54285e-07
//    4  p3          -1.79892e-05   5.44067e-06   2.15870e-09   4.11441e-04


// Theta
// root [14] rms->Fit(f,"","",0.1,3)
//  FCN=8.9981 FROM MIGRAD    STATUS=CONVERGED      79 CALLS          80 TOTAL
//                      EDM=3.03049e-17    STRATEGY= 1      ERR MATRIX NOT POS-DEF
//  EXT PARAMETER                APPROXIMATE        STEP         FIRST
//  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE
//     1  p0          -1.68773e-02   2.67644e-05   8.04770e-09   4.48079e-05
//     2  p1           1.78440e-02   2.65467e-05   8.50867e-09   6.43012e-05
//     3  p2          -5.26559e-02   5.06308e-04   2.28595e-07  -1.63963e-05
//     4  p3           2.00940e-04   1.14440e-05   3.98737e-09   1.78198e-05


#include <TH1.h>
#include <TH3.h>
#include <TF1.h>
#include <TRandom.h>
#include <AliAODParticle.h>


ClassImp(AliHBTCorrectedCorrelFctn)

AliHBTCorrectedCorrelFctn::AliHBTCorrectedCorrelFctn():
  fDPtOverPtA(5.78220e-03),
  fDPtOverPtB(3.98063e-05),
  fDPtOverPtAlpha(-2.78008),
  fDPtOverPtC(5.07594e-04),
  fThetaA(5.87693e-04),
  fThetaB(2.16488e-06),
  fThetaAlpha(-3.10218e+00),
  fThetaC(-1.79892e-05),
  fPhiA(-1.68773e-02),
  fPhiB(1.78440e-02),
  fPhiAlpha(-5.26559e-02),
  fPhiC(2.00940e-04)
{
  //ctor
}

/******************************************************************/
void AliHBTCorrectedCorrelFctn::Smear(AliHBTPair* pair,AliHBTPair& smeared)
{
//Smears pair
  Smear(pair->Particle1(),smeared.Particle1());
  Smear(pair->Particle2(),smeared.Particle2());
  smeared.Changed();
}
/******************************************************************/

void AliHBTCorrectedCorrelFctn::Smear(AliVAODParticle* part, AliVAODParticle* smeared)
{
 //Smears momenta
  Double_t sin2theta = TMath::Sin(part->Theta());
  sin2theta = sin2theta*sin2theta;
  Double_t pt = part->Pt();
  
  Double_t sigmapt = fDPtOverPtA + fDPtOverPtB*TMath::Power(pt,fDPtOverPtAlpha) + fDPtOverPtC*pt;
  Double_t dPtDivPt = gRandom->Gaus(0.0,sigmapt);
  Double_t dphi = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fPhiAlpha) + fPhiC*pt);
  Double_t dtheta = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fThetaAlpha) +fThetaC*pt);
  
  Double_t smearedPx = part->Px()*(1.0+dPtDivPt) - part->Py()*dphi;
//  fourmom.setX(px*(1.0+dPtDivPt) - py*dphi);
  Double_t smearedPy = part->Py()*(1.0+dPtDivPt) - part->Px()*dphi;
//  fourmom.setY(py*(1.0+dPtDivPt) + px*dphi);
  Double_t smearedPz = part->Pz()*(1.0+dPtDivPt) - pt*dtheta/sin2theta;
//  fourmom.setZ(pz*(1.0+dPtDivPt) - pT*dtheta/sin2theta);
  
  Double_t mass2 = part->Mass()*part->Mass();
  Double_t e = mass2 + smearedPx*smearedPx + 
                       smearedPy*smearedPy + 
                       smearedPz*smearedPz;
	   
  smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(e));
}

/****************************************************************/
/****************************************************************/
/****************************************************************/


ClassImp(AliHBTCorrectQInvCorrelFctn)

AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name,const char* title):
  AliHBTOnePairFctn1D(name,title),
  fMeasCorrelFctn(0x0),
  fMeasNumer(0x0),
  fMeasDenom(0x0),
  fSmearedNumer(0x0),
  fSmearedDenom(0x0),
  fR2(0.0),
  fLambda(0.0),
  fFittedR(0),
  fFittedLambda(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),
  fR2(0.0),
  fLambda(0.0),
  fFittedR(0),
  fFittedLambda(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),
  AliHBTCorrectedCorrelFctn(),
  fMeasCorrelFctn(0x0),
  fMeasNumer(0x0),
  fMeasDenom(0x0),
  fSmearedNumer(0x0),
  fSmearedDenom(0x0),
  fR2(0.0),
  fLambda(0.0),
  fFittedR(0),
  fFittedLambda(0),
  fRConvergenceTreshold(0.3),
  fLambdaConvergenceTreshold(0.05)
{
//ctor
}	            
/******************************************************************/
AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const AliHBTCorrectQInvCorrelFctn& in):
  AliHBTOnePairFctn1D(in),
  AliHBTCorrectedCorrelFctn(),
  fMeasCorrelFctn(0x0),
  fMeasNumer(0x0),
  fMeasDenom(0x0),
  fSmearedNumer(0x0),
  fSmearedDenom(0x0),
  fR2(0.0),
  fLambda(0.0),
  fFittedR(0),
  fFittedLambda(0),
  fRConvergenceTreshold(0),
  fLambdaConvergenceTreshold(0)
{
//cpy ;ctor
 in.Copy(*this);
}

AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
{
//dtor
  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()
{
//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)
{
 //Processes particles that originates from the same event
  if (fMeasNumer == 0x0) return;
  pair = CheckPair(pair);
  if( pair == 0x0) return;
  fMeasNumer->Fill(pair->GetQInv());
}
/******************************************************************/

void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair)
{
//Process different events 
  static AliAODParticle 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::SetInitialValues(Double_t lambda, Double_t r)
{
 //Sets Initial Values
  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()
{
 //Returns Fitted radius
  if (fFittedR <= 0.0) Fit();
  return fFittedR;
}
/******************************************************************/

Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
{
 //Returns Fitted Intercept paramter
  if (fFittedR <= 0.0) Fit();
  return fFittedLambda;
}
/******************************************************************/

void AliHBTCorrectQInvCorrelFctn::WriteAll()
{
//Writes function and all additional information
  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);
   }
}
Commit	Line	Data
e09fa876	1	#include "AliHBTCorrectQInvCorrelFctn.h"
83b33650	2	//____________________
	3	///////////////////////////////////////////////////////
	4	// //
	5	// AliHBTCorrectQInvCorrelFctn //
	6	// //
	7	// Class for calculating Q Invariant correlation //
	8	// taking to the account resolution of the //
	9	// detector and coulomb effects. //
	10	// Implemented on basis of STAR procedure //
	11	// implemented and described by Michael A. Lisa //
	12	// //
	13	// Piotr.Skowronski@cern.ch //
c7ffd78f	14	// http://aliweb.cern.ch/people/skowron/analyzer //
83b33650	15	// //
	16	///////////////////////////////////////////////////////
	17
62e1b4fe	18	// dPt/Pt
	19	// root [19] rms->Fit(f,"","",0.1,2)
	20	// FCN=7.0017 FROM MIGRAD STATUS=CONVERGED 126 CALLS 127 TOTAL
	21	// EDM=2.28804e-15 STRATEGY= 1 ERROR MATRIX ACCURATE
	22	// EXT PARAMETER STEP FIRST
	23	// NO. NAME VALUE ERROR SIZE DERIVATIVE
	24	// 1 p0 5.78220e-03 3.14576e-05 4.97822e-09 -1.90059e-05
	25	// 2 p1 3.98063e-05 1.61877e-06 1.04380e-10 1.91454e-04
	26	// 3 p2 -2.78008e+00 2.13581e-02 1.66031e-06 3.16574e-06
	27	// 4 p3 5.07594e-04 4.79619e-05 1.29793e-08 -2.29242e-05
	28
	29
	30	// Phi
	31	// root [17] rms->Fit(f,"","",0.15,2.5)
	32	// Warning in <TH1D::Fit>: Abnormal termination of minimization.
	33	// FCN=33.4898 FROM MIGRAD STATUS=FAILED 91 CALLS 92 TOTAL
	34	// EDM=1.19154e-15 STRATEGY= 1 ERR MATRIX APPROXIMATE
	35	// EXT PARAMETER APPROXIMATE STEP FIRST
	36	// NO. NAME VALUE ERROR SIZE DERIVATIVE
	37	// 1 p0 5.87693e-04 5.04254e-06 2.49187e-09 5.84546e-04
	38	// 2 p1 2.16488e-06 3.68880e-07 6.41507e-11 -7.36564e-02
	39	// 3 p2 -3.10218e+00 1.01695e-01 2.00177e-05 7.54285e-07
	40	// 4 p3 -1.79892e-05 5.44067e-06 2.15870e-09 4.11441e-04
	41
	42
	43
	44	// Theta
	45	// root [14] rms->Fit(f,"","",0.1,3)
	46	// FCN=8.9981 FROM MIGRAD STATUS=CONVERGED 79 CALLS 80 TOTAL
	47	// EDM=3.03049e-17 STRATEGY= 1 ERR MATRIX NOT POS-DEF
	48	// EXT PARAMETER APPROXIMATE STEP FIRST
	49	// NO. NAME VALUE ERROR SIZE DERIVATIVE
	50	// 1 p0 -1.68773e-02 2.67644e-05 8.04770e-09 4.48079e-05
	51	// 2 p1 1.78440e-02 2.65467e-05 8.50867e-09 6.43012e-05
	52	// 3 p2 -5.26559e-02 5.06308e-04 2.28595e-07 -1.63963e-05
	53	// 4 p3 2.00940e-04 1.14440e-05 3.98737e-09 1.78198e-05
83b33650	54
83b33650	55
	56	#include <TH1.h>
	57	#include <TH3.h>
	58	#include <TF1.h>
	59	#include <TRandom.h>
78d7c6d3	60	#include <AliAODParticle.h>
	61
	62
62e1b4fe	63	ClassImp(AliHBTCorrectedCorrelFctn)
	64
	65	AliHBTCorrectedCorrelFctn::AliHBTCorrectedCorrelFctn():
	66	fDPtOverPtA(5.78220e-03),
	67	fDPtOverPtB(3.98063e-05),
	68	fDPtOverPtAlpha(-2.78008),
	69	fDPtOverPtC(5.07594e-04),
	70	fThetaA(5.87693e-04),
	71	fThetaB(2.16488e-06),
	72	fThetaAlpha(-3.10218e+00),
	73	fThetaC(-1.79892e-05),
	74	fPhiA(-1.68773e-02),
	75	fPhiB(1.78440e-02),
	76	fPhiAlpha(-5.26559e-02),
	77	fPhiC(2.00940e-04)
	78	{
	79	//ctor
	80	}
	81
	82	/******************************************************************/
	83	void AliHBTCorrectedCorrelFctn::Smear(AliHBTPair* pair,AliHBTPair& smeared)
	84	{
	85	//Smears pair
	86	Smear(pair->Particle1(),smeared.Particle1());
	87	Smear(pair->Particle2(),smeared.Particle2());
	88	smeared.Changed();
	89	}
	90	/******************************************************************/
	91
	92	void AliHBTCorrectedCorrelFctn::Smear(AliVAODParticle* part, AliVAODParticle* smeared)
	93	{
	94	//Smears momenta
	95	Double_t sin2theta = TMath::Sin(part->Theta());
	96	sin2theta = sin2theta*sin2theta;
	97	Double_t pt = part->Pt();
	98
	99	Double_t sigmapt = fDPtOverPtA + fDPtOverPtBTMath::Power(pt,fDPtOverPtAlpha) + fDPtOverPtCpt;
	100	Double_t dPtDivPt = gRandom->Gaus(0.0,sigmapt);
	101	Double_t dphi = gRandom->Gaus(0.0,fPhiA+fPhiBTMath::Power(pt,fPhiAlpha) + fPhiCpt);
	102	Double_t dtheta = gRandom->Gaus(0.0,fPhiA+fPhiBTMath::Power(pt,fThetaAlpha) +fThetaCpt);
	103
	104	Double_t smearedPx = part->Px()(1.0+dPtDivPt) - part->Py()dphi;
	105	// fourmom.setX(px(1.0+dPtDivPt) - pydphi);
	106	Double_t smearedPy = part->Py()(1.0+dPtDivPt) - part->Px()dphi;
	107	// fourmom.setY(py(1.0+dPtDivPt) + pxdphi);
	108	Double_t smearedPz = part->Pz()(1.0+dPtDivPt) - ptdtheta/sin2theta;
	109	// fourmom.setZ(pz(1.0+dPtDivPt) - pTdtheta/sin2theta);
	110
	111	Double_t mass2 = part->Mass()*part->Mass();
	112	Double_t e = mass2 + smearedPx*smearedPx +
	113	smearedPy*smearedPy +
	114	smearedPz*smearedPz;
	115
	116	smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(e));
	117	}
	118
	119	/****************************************************************/
	120	/****************************************************************/
	121	/****************************************************************/
	122
	123
83b33650	124	ClassImp(AliHBTCorrectQInvCorrelFctn)
	125
	126	AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name,const char* title):
	127	AliHBTOnePairFctn1D(name,title),
	128	fMeasCorrelFctn(0x0),
	129	fMeasNumer(0x0),
	130	fMeasDenom(0x0),
	131	fSmearedNumer(0x0),
	132	fSmearedDenom(0x0),
83b33650	133	fR2(0.0),
83b33650	134	fLambda(0.0),
4b1c9620	135	fFittedR(0),
4b1c9620	136	fFittedLambda(0),
83b33650	137	fRConvergenceTreshold(0.3),
	138	fLambdaConvergenceTreshold(0.05)
	139	{
	140
	141	}
	142	/******************************************************************/
	143
	144	AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(TH1D* measqinv,const char* name,const char* title):
	145	AliHBTOnePairFctn1D(name,title,
	146	measqinv->GetNbinsX(),
	147	measqinv->GetXaxis()->GetXmax(),
	148	measqinv->GetXaxis()->GetXmin()),
	149	fMeasCorrelFctn(measqinv),
	150	fMeasNumer(0x0),
	151	fMeasDenom(0x0),
	152	fSmearedNumer(0x0),
	153	fSmearedDenom(0x0),
83b33650	154	fR2(0.0),
83b33650	155	fLambda(0.0),
4b1c9620	156	fFittedR(0),
4b1c9620	157	fFittedLambda(0),
83b33650	158	fRConvergenceTreshold(0.3),
	159	fLambdaConvergenceTreshold(0.05)
	160	{
	161
	162	}
	163	/******************************************************************/
	164
	165	AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name, const char* title,
	166	Int_t nbins, Float_t maxXval, Float_t minXval):
	167	AliHBTOnePairFctn1D(name,title,nbins,maxXval,minXval),
62e1b4fe	168	AliHBTCorrectedCorrelFctn(),
83b33650	169	fMeasCorrelFctn(0x0),
	170	fMeasNumer(0x0),
	171	fMeasDenom(0x0),
	172	fSmearedNumer(0x0),
	173	fSmearedDenom(0x0),
83b33650	174	fR2(0.0),
83b33650	175	fLambda(0.0),
4b1c9620	176	fFittedR(0),
4b1c9620	177	fFittedLambda(0),
83b33650	178	fRConvergenceTreshold(0.3),
	179	fLambdaConvergenceTreshold(0.05)
	180	{
e09fa876	181	//ctor
83b33650	182	}
83b33650	183	/******************************************************************/
e09fa876	184	AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const AliHBTCorrectQInvCorrelFctn& in):
e09fa876	185	AliHBTOnePairFctn1D(in),
62e1b4fe	186	AliHBTCorrectedCorrelFctn(),
e09fa876	187	fMeasCorrelFctn(0x0),
	188	fMeasNumer(0x0),
	189	fMeasDenom(0x0),
	190	fSmearedNumer(0x0),
	191	fSmearedDenom(0x0),
e09fa876	192	fR2(0.0),
e09fa876	193	fLambda(0.0),
4b1c9620	194	fFittedR(0),
4b1c9620	195	fFittedLambda(0),
e09fa876	196	fRConvergenceTreshold(0),
	197	fLambdaConvergenceTreshold(0)
	198	{
	199	//cpy ;ctor
	200	in.Copy(*this);
	201	}
83b33650	202
	203	AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
	204	{
e09fa876	205	//dtor
83b33650	206	delete fMeasCorrelFctn;
	207	delete fSmearedNumer;
	208	delete fSmearedDenom;
	209	delete fMeasNumer;
	210	delete fMeasDenom;
	211	}
	212	/******************************************************************/
	213
	214	void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min)
	215	{
	216	AliHBTFunction1D::BuildHistos(nbins,max,min);
	217
	218	TString numstr = fName + " Smeared Numerator"; //title and name of the numerator histogram
	219	TString denstr = fName + " Smeared Denominator";//title and name of the denominator histogram
	220
	221	fSmearedNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
	222	fSmearedDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
	223	fSmearedNumer->Sumw2();
	224	fSmearedDenom->Sumw2();
	225
	226	if (fMeasCorrelFctn == 0x0)
	227	{
	228	numstr = fName + " Measured Numerator"; //title and name of the numerator histogram
	229	denstr = fName + " Measured Denominator";//title and name of the denominator histogram
	230
	231	fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
	232	fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
	233	fMeasNumer->Sumw2();
	234	fMeasDenom->Sumw2();
	235	}
	236	}
	237	/******************************************************************/
	238
	239	void AliHBTCorrectQInvCorrelFctn::Init()
	240	{
e09fa876	241	//Init
83b33650	242	AliHBTOnePairFctn1D::Init();
	243	Info("Init","");
	244	fSmearedNumer->Reset();
	245	fSmearedDenom->Reset();
	246	if (fMeasNumer) fMeasNumer->Reset();
	247	if (fMeasDenom) fMeasDenom->Reset();
	248	fFittedR = -1.0;
	249	fFittedLambda = 0.0;
	250	}
	251	/******************************************************************/
	252
	253	void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair)
	254	{
e09fa876	255	//Processes particles that originates from the same event
83b33650	256	if (fMeasNumer == 0x0) return;
	257	pair = CheckPair(pair);
	258	if( pair == 0x0) return;
	259	fMeasNumer->Fill(pair->GetQInv());
	260	}
	261	/******************************************************************/
	262
	263	void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair)
	264	{
	265	//Process different events
78d7c6d3	266	static AliAODParticle part1, part2;
83b33650	267	static AliHBTPair smearedpair(&part1,&part2);
	268
	269	pair = CheckPair(pair);
	270	if( pair == 0x0) return;
	271	Double_t cc = GetCoulombCorrection(pair);
	272	Double_t qinv = pair->GetQInv();
	273
	274	//measured histogram -> if we are interested
	275	//only if fMeasCorrelFctn is not specified by user
	276	if (fMeasDenom) fMeasDenom->Fill(qinv,cc);
	277
	278	Smear(pair,smearedpair);
	279	Double_t modelqinv = GetModelValue(qinv);
	280	//Ideal histogram
	281	fNumerator->Fill(qinv,modelqinv*cc);
	282	fDenominator->Fill(qinv,cc);
	283
	284	//Smeared histogram
	285	Double_t smearedqinv = smearedpair.GetQInv();
	286	fSmearedNumer->Fill(smearedqinv,modelqinv);
	287	Double_t smearedcc = GetCoulombCorrection(&smearedpair);
	288	fSmearedDenom->Fill(smearedqinv,smearedcc);
	289
	290	}
	291	/******************************************************************/
83b33650	292
	293	void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
	294	{
e09fa876	295	//Sets Initial Values
83b33650	296	fLambda = lambda;
	297	fR2 = r*r;
	298	}
	299	/******************************************************************/
	300	void AliHBTCorrectQInvCorrelFctn::MakeMeasCF()
	301	{
	302	//makes measured correlation function
	303	delete fMeasCorrelFctn;
	304	fMeasCorrelFctn = 0x0;
	305
	306	if (fMeasNumer&&fMeasDenom)
	307	{
	308	Double_t measscale = Scale(fMeasNumer,fMeasDenom);
	309	if (measscale == 0.0)
	310	{
	311	Error("GetResult","GetRatio for measured CF returned 0.0");
	312	return;
	313	}
	314	TString str = fName + "measured ratio";
	315	fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data());
	316	fMeasCorrelFctn->SetTitle(str.Data());
	317	fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale);
	318	}
	319	}
	320
	321	TH1* AliHBTCorrectQInvCorrelFctn::GetResult()
	322	{
	323	//In case we don't have yet Measured Correlation Function
e09fa876	324	//Try to get it
83b33650	325	//result is
	326	// N[meas] N[ideal]/D[ideal]
	327	//C(Q) = ------- * -----------------
	328	// D[meas] N[smear]/D[smear]
	329
	330	TString str;
	331	if (fMeasCorrelFctn == 0x0) MakeMeasCF();
	332
	333	if (fMeasCorrelFctn == 0x0)
	334	{
	335	Error("GetResult",
	336	"Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null");
	337	return 0x0;
	338	}
	339
	340	TH1D* ideal = (TH1D*)GetRatio(Scale());
	341	if (ideal == 0x0)
	342	{
	343	Error("GetResult","Ratio of ideal histograms is null");
	344	return 0x0;
	345	}
	346	str = fName + " smeared ratio";
	347	TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
	348	smearedCF->SetTitle(str.Data());
	349	Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
	350	smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
	351
	352	str = fName + " product meas ideal CF";
	353	TH1D* measideal = (TH1D*)ideal->Clone(str.Data());
	354	measideal->Multiply(ideal,fMeasCorrelFctn);
	355
	356	str = fName + " Corrected Result";
	357	TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data());
	358	result->SetTitle(str.Data());
	359
	360	Double_t resultscale = Scale(measideal,smearedCF);
	361	result->Divide(measideal,smearedCF,resultscale);
	362	return result;
	363
	364	}
	365	/******************************************************************/
	366
	367	void AliHBTCorrectQInvCorrelFctn::Fit()
	368	{
	369	//fits resuting histogram with function 1.0 + [0]exp([1][1]xx/(-0.038936366329))
	370	//where [0] is lambda
	371	// [1] is radius
	372	// 0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI
	373
	374	Info("Fit","Before fFittedLambda = %f",fFittedLambda);
	375	Info("Fit","Before fFittedR = %f",fFittedR);
	376	TH1D* result = (TH1D*)GetResult();
	377	if (result == 0x0)
	378	{
	379	Error("Fit","Can not get result");
	380	return;
	381	}
	382	TF1* fitfctn = new TF1("fitfctn","1.0 + [0]exp([1][1]xx/(-0.038936366329))");
	383	fitfctn->SetParameter(0,1.0);
	384	fitfctn->SetParameter(1,6.0);
	385	Float_t max = result->GetXaxis()->GetXmax();
	386	Info("Fit","Max is %f",max);
	387	result->Fit(fitfctn,"0","",0.008,max);
	388	fFittedLambda = fitfctn->GetParameter(0);
389	fFittedR = fitfctn->GetParameter(1);
390	Info("Fit","After fFittedLambda = %f",fFittedLambda);
391	Info("Fit","After fFittedR = %f",fFittedR);
392	delete fitfctn;
393	delete result;
394	}
395	/******************************************************************/
396
397	Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged()
398	{
399	//check if fitting was performed
400	if (fFittedR <= 0.0)
401	{
402	Fit();//if not do fit first
403	if (fFittedR <= 0.0)
404	{
405	Error("IsConverged","Fitting failed");
406	return kFALSE;
407	}
408	}
409
410	Double_t guessedR = TMath::Sqrt(fR2);
411
412	Info("IsConverged","Fitted lambda : %8f Fitted Radius : %8f",fFittedLambda,fFittedR);
413	Info("IsConverged","Guessed lambda : %8f Guessed Radius : %8f",fLambda,guessedR);
414	Info("IsConverged","Demanded lambda convergence: %8f Demanded Radius convergence: %8f",
415	fLambdaConvergenceTreshold,fRConvergenceTreshold);
416
417	if ( (TMath::Abs(fLambda-fFittedLambda)<fLambdaConvergenceTreshold) &&
418	(TMath::Abs(fFittedR-guessedR)<fRConvergenceTreshold) )
419	{
420	Info("IsConverged","Cnvergence reached");
421	return kTRUE;
422	}
423	else
424	{
425	Info("IsConverged","Cnvergence NOT reached");
426	return kFALSE;
427	}
428	}
429	/******************************************************************/
430
431	Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius()
432	{
e09fa876	433	//Returns Fitted radius
83b33650	434	if (fFittedR <= 0.0) Fit();
	435	return fFittedR;
	436	}
	437	/******************************************************************/
	438
	439	Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
	440	{
e09fa876	441	//Returns Fitted Intercept paramter
83b33650	442	if (fFittedR <= 0.0) Fit();
	443	return fFittedLambda;
	444	}
	445	/******************************************************************/
	446
	447	void AliHBTCorrectQInvCorrelFctn::WriteAll()
	448	{
e09fa876	449	//Writes function and all additional information
83b33650	450	Write();
	451	if (fMeasCorrelFctn) fMeasCorrelFctn->Write();
	452	if (fMeasNumer ) fMeasNumer->Write();
	453	if (fMeasDenom ) fMeasDenom->Write();
	454	if (fSmearedNumer) fSmearedNumer->Write();
	455	if (fSmearedDenom) fSmearedDenom->Write();
	456	if (fSmearedNumer && fSmearedDenom)
	457	{
	458	TString str = fName + " smeared ratio";
	459	TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
	460	smearedCF->SetTitle(str.Data());
	461	Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
	462	smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
	463	}
	464	}
	465