[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://alisoft.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),
  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),
  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),
  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),
  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 //
	14	// http://alisoft.cern.ch/people/skowron/analyzer //
	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),
	134	fLambda(0.0),
	135	fRConvergenceTreshold(0.3),
	136	fLambdaConvergenceTreshold(0.05)
	137	{
	138
	139	}
	140	/******************************************************************/
	141
	142	AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(TH1D* measqinv,const char* name,const char* title):
	143	AliHBTOnePairFctn1D(name,title,
	144	measqinv->GetNbinsX(),
	145	measqinv->GetXaxis()->GetXmax(),
	146	measqinv->GetXaxis()->GetXmin()),
	147	fMeasCorrelFctn(measqinv),
	148	fMeasNumer(0x0),
	149	fMeasDenom(0x0),
	150	fSmearedNumer(0x0),
	151	fSmearedDenom(0x0),
83b33650	152	fR2(0.0),
	153	fLambda(0.0),
	154	fRConvergenceTreshold(0.3),
	155	fLambdaConvergenceTreshold(0.05)
	156	{
	157
	158	}
	159	/******************************************************************/
	160
	161	AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name, const char* title,
	162	Int_t nbins, Float_t maxXval, Float_t minXval):
	163	AliHBTOnePairFctn1D(name,title,nbins,maxXval,minXval),
62e1b4fe	164	AliHBTCorrectedCorrelFctn(),
83b33650	165	fMeasCorrelFctn(0x0),
	166	fMeasNumer(0x0),
	167	fMeasDenom(0x0),
	168	fSmearedNumer(0x0),
	169	fSmearedDenom(0x0),
83b33650	170	fR2(0.0),
	171	fLambda(0.0),
	172	fRConvergenceTreshold(0.3),
	173	fLambdaConvergenceTreshold(0.05)
	174	{
e09fa876	175	//ctor
83b33650	176	}
83b33650	177	/******************************************************************/
e09fa876	178	AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const AliHBTCorrectQInvCorrelFctn& in):
e09fa876	179	AliHBTOnePairFctn1D(in),
62e1b4fe	180	AliHBTCorrectedCorrelFctn(),
e09fa876	181	fMeasCorrelFctn(0x0),
	182	fMeasNumer(0x0),
	183	fMeasDenom(0x0),
	184	fSmearedNumer(0x0),
	185	fSmearedDenom(0x0),
e09fa876	186	fR2(0.0),
	187	fLambda(0.0),
	188	fRConvergenceTreshold(0),
	189	fLambdaConvergenceTreshold(0)
	190	{
	191	//cpy ;ctor
	192	in.Copy(*this);
	193	}
83b33650	194
	195	AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
	196	{
e09fa876	197	//dtor
83b33650	198	delete fMeasCorrelFctn;
	199	delete fSmearedNumer;
	200	delete fSmearedDenom;
	201	delete fMeasNumer;
	202	delete fMeasDenom;
	203	}
	204	/******************************************************************/
	205
	206	void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min)
	207	{
	208	AliHBTFunction1D::BuildHistos(nbins,max,min);
	209
	210	TString numstr = fName + " Smeared Numerator"; //title and name of the numerator histogram
	211	TString denstr = fName + " Smeared Denominator";//title and name of the denominator histogram
	212
	213	fSmearedNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
	214	fSmearedDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
	215	fSmearedNumer->Sumw2();
	216	fSmearedDenom->Sumw2();
	217
	218	if (fMeasCorrelFctn == 0x0)
	219	{
	220	numstr = fName + " Measured Numerator"; //title and name of the numerator histogram
	221	denstr = fName + " Measured Denominator";//title and name of the denominator histogram
	222
	223	fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
	224	fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
	225	fMeasNumer->Sumw2();
	226	fMeasDenom->Sumw2();
	227	}
	228	}
	229	/******************************************************************/
	230
	231	void AliHBTCorrectQInvCorrelFctn::Init()
	232	{
e09fa876	233	//Init
83b33650	234	AliHBTOnePairFctn1D::Init();
	235	Info("Init","");
	236	fSmearedNumer->Reset();
	237	fSmearedDenom->Reset();
	238	if (fMeasNumer) fMeasNumer->Reset();
	239	if (fMeasDenom) fMeasDenom->Reset();
	240	fFittedR = -1.0;
	241	fFittedLambda = 0.0;
	242	}
	243	/******************************************************************/
	244
	245	void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair)
	246	{
e09fa876	247	//Processes particles that originates from the same event
83b33650	248	if (fMeasNumer == 0x0) return;
	249	pair = CheckPair(pair);
	250	if( pair == 0x0) return;
	251	fMeasNumer->Fill(pair->GetQInv());
	252	}
	253	/******************************************************************/
	254
	255	void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair)
	256	{
	257	//Process different events
78d7c6d3	258	static AliAODParticle part1, part2;
83b33650	259	static AliHBTPair smearedpair(&part1,&part2);
	260
	261	pair = CheckPair(pair);
	262	if( pair == 0x0) return;
	263	Double_t cc = GetCoulombCorrection(pair);
	264	Double_t qinv = pair->GetQInv();
	265
	266	//measured histogram -> if we are interested
	267	//only if fMeasCorrelFctn is not specified by user
	268	if (fMeasDenom) fMeasDenom->Fill(qinv,cc);
	269
	270	Smear(pair,smearedpair);
	271	Double_t modelqinv = GetModelValue(qinv);
	272	//Ideal histogram
	273	fNumerator->Fill(qinv,modelqinv*cc);
	274	fDenominator->Fill(qinv,cc);
	275
	276	//Smeared histogram
	277	Double_t smearedqinv = smearedpair.GetQInv();
	278	fSmearedNumer->Fill(smearedqinv,modelqinv);
	279	Double_t smearedcc = GetCoulombCorrection(&smearedpair);
	280	fSmearedDenom->Fill(smearedqinv,smearedcc);
	281
	282	}
	283	/******************************************************************/
83b33650	284
	285	void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
	286	{
e09fa876	287	//Sets Initial Values
83b33650	288	fLambda = lambda;
	289	fR2 = r*r;
	290	}
	291	/******************************************************************/
	292	void AliHBTCorrectQInvCorrelFctn::MakeMeasCF()
	293	{
	294	//makes measured correlation function
	295	delete fMeasCorrelFctn;
	296	fMeasCorrelFctn = 0x0;
	297
	298	if (fMeasNumer&&fMeasDenom)
	299	{
	300	Double_t measscale = Scale(fMeasNumer,fMeasDenom);
	301	if (measscale == 0.0)
	302	{
	303	Error("GetResult","GetRatio for measured CF returned 0.0");
	304	return;
	305	}
	306	TString str = fName + "measured ratio";
	307	fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data());
	308	fMeasCorrelFctn->SetTitle(str.Data());
	309	fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale);
	310	}
	311	}
	312
	313	TH1* AliHBTCorrectQInvCorrelFctn::GetResult()
	314	{
	315	//In case we don't have yet Measured Correlation Function
e09fa876	316	//Try to get it
83b33650	317	//result is
	318	// N[meas] N[ideal]/D[ideal]
	319	//C(Q) = ------- * -----------------
	320	// D[meas] N[smear]/D[smear]
	321
	322	TString str;
	323	if (fMeasCorrelFctn == 0x0) MakeMeasCF();
	324
	325	if (fMeasCorrelFctn == 0x0)
	326	{
	327	Error("GetResult",
	328	"Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null");
	329	return 0x0;
	330	}
	331
	332	TH1D* ideal = (TH1D*)GetRatio(Scale());
	333	if (ideal == 0x0)
	334	{
	335	Error("GetResult","Ratio of ideal histograms is null");
	336	return 0x0;
	337	}
	338	str = fName + " smeared ratio";
	339	TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
	340	smearedCF->SetTitle(str.Data());
	341	Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
	342	smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
	343
	344	str = fName + " product meas ideal CF";
	345	TH1D* measideal = (TH1D*)ideal->Clone(str.Data());
	346	measideal->Multiply(ideal,fMeasCorrelFctn);
	347
	348	str = fName + " Corrected Result";
	349	TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data());
	350	result->SetTitle(str.Data());
	351
	352	Double_t resultscale = Scale(measideal,smearedCF);
	353	result->Divide(measideal,smearedCF,resultscale);
	354	return result;
	355
	356	}
	357	/******************************************************************/
	358
	359	void AliHBTCorrectQInvCorrelFctn::Fit()
	360	{
	361	//fits resuting histogram with function 1.0 + [0]exp([1][1]xx/(-0.038936366329))
	362	//where [0] is lambda
	363	// [1] is radius
	364	// 0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI
	365
	366	Info("Fit","Before fFittedLambda = %f",fFittedLambda);
	367	Info("Fit","Before fFittedR = %f",fFittedR);
	368	TH1D* result = (TH1D*)GetResult();
	369	if (result == 0x0)
	370	{
	371	Error("Fit","Can not get result");
	372	return;
	373	}
	374	TF1* fitfctn = new TF1("fitfctn","1.0 + [0]exp([1][1]xx/(-0.038936366329))");
	375	fitfctn->SetParameter(0,1.0);
	376	fitfctn->SetParameter(1,6.0);
	377	Float_t max = result->GetXaxis()->GetXmax();
	378	Info("Fit","Max is %f",max);
	379	result->Fit(fitfctn,"0","",0.008,max);
	380	fFittedLambda = fitfctn->GetParameter(0);
381	fFittedR = fitfctn->GetParameter(1);
382	Info("Fit","After fFittedLambda = %f",fFittedLambda);
383	Info("Fit","After fFittedR = %f",fFittedR);
384	delete fitfctn;
385	delete result;
386	}
387	/******************************************************************/
388
389	Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged()
390	{
391	//check if fitting was performed
392	if (fFittedR <= 0.0)
393	{
394	Fit();//if not do fit first
395	if (fFittedR <= 0.0)
396	{
397	Error("IsConverged","Fitting failed");
398	return kFALSE;
399	}
400	}
401
402	Double_t guessedR = TMath::Sqrt(fR2);
403
404	Info("IsConverged","Fitted lambda : %8f Fitted Radius : %8f",fFittedLambda,fFittedR);
405	Info("IsConverged","Guessed lambda : %8f Guessed Radius : %8f",fLambda,guessedR);
406	Info("IsConverged","Demanded lambda convergence: %8f Demanded Radius convergence: %8f",
407	fLambdaConvergenceTreshold,fRConvergenceTreshold);
408
409	if ( (TMath::Abs(fLambda-fFittedLambda)<fLambdaConvergenceTreshold) &&
410	(TMath::Abs(fFittedR-guessedR)<fRConvergenceTreshold) )
411	{
412	Info("IsConverged","Cnvergence reached");
413	return kTRUE;
414	}
415	else
416	{
417	Info("IsConverged","Cnvergence NOT reached");
418	return kFALSE;
419	}
420	}
421	/******************************************************************/
422
423	Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius()
424	{
e09fa876	425	//Returns Fitted radius
83b33650	426	if (fFittedR <= 0.0) Fit();
	427	return fFittedR;
	428	}
	429	/******************************************************************/
	430
	431	Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
	432	{
e09fa876	433	//Returns Fitted Intercept paramter
83b33650	434	if (fFittedR <= 0.0) Fit();
	435	return fFittedLambda;
	436	}
	437	/******************************************************************/
	438
	439	void AliHBTCorrectQInvCorrelFctn::WriteAll()
	440	{
e09fa876	441	//Writes function and all additional information
83b33650	442	Write();
	443	if (fMeasCorrelFctn) fMeasCorrelFctn->Write();
	444	if (fMeasNumer ) fMeasNumer->Write();
	445	if (fMeasDenom ) fMeasDenom->Write();
	446	if (fSmearedNumer) fSmearedNumer->Write();
	447	if (fSmearedDenom) fSmearedDenom->Write();
	448	if (fSmearedNumer && fSmearedDenom)
	449	{
	450	TString str = fName + " smeared ratio";
	451	TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
	452	smearedCF->SetTitle(str.Data());
	453	Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
	454	smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
	455	}
	456	}
	457