[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    //
//                                                   //
///////////////////////////////////////////////////////

//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>
#include <AliAODParticle.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)
{
//ctor
}	            
/******************************************************************/
AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const AliHBTCorrectQInvCorrelFctn& in):
  AliHBTOnePairFctn1D(in),
  fMeasCorrelFctn(0x0),
  fMeasNumer(0x0),
  fMeasDenom(0x0),
  fSmearedNumer(0x0),
  fSmearedDenom(0x0),
  fDPtOverPtRMS(0),
  fThetaA(0),
  fThetaB(0),
  fThetaAlpha(0),
  fPhiA(0),
  fPhiB(0),
  fPhiAlpha(0),
  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::Smear(AliHBTPair* pair,AliHBTPair& smeared)
{
//Smears pair
  Smear(pair->Particle1(),smeared.Particle1());
  Smear(pair->Particle2(),smeared.Particle2());
  smeared.Changed();
}
/******************************************************************/

void AliHBTCorrectQInvCorrelFctn::Smear(AliVAODParticle* part, AliVAODParticle* smeared)
{
 //Smears momenta
  Double_t sin2theta = TMath::Sin(part->Theta());
  sin2theta = sin2theta*sin2theta;
  Double_t pt = part->Pt();

  double dPtDivPt = 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+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));
}
/******************************************************************/

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
	18	//Parameters fit from pi+ pi+ resolution analysis for pair with qinv<50MeV
	19	// chi2/NFD = 97/157
	20	// FCN=98.0971 FROM MIGRAD STATUS=CONVERGED 332 CALLS 333 TOTAL
	21	// EDM=2.13364e-12 STRATEGY= 1 ERROR MATRIX UNCERTAINTY 2.4 per cent
	22	// EXT PARAMETER STEP FIRST
	23	// NO. NAME VALUE ERROR SIZE DERIVATIVE
	24	// 1 fThetaA 2.72546e-03 1.43905e-04 -0.00000e+00 5.54375e-02
	25	// 2 fThetaB 1.87116e-04 5.11862e-05 0.00000e+00 3.66500e-01
	26	// 3 fThetaAlpha -2.36868e+00 1.83230e-01 0.00000e+00 -7.01301e-05
	27
	28	// FCN=120.603 FROM MIGRAD STATUS=CONVERGED 117 CALLS 118 TOTAL
	29	// EDM=2.5153e-12 STRATEGY= 1 ERROR MATRIX UNCERTAINTY 2.4 per cent
	30	// EXT PARAMETER STEP FIRST
	31	// NO. NAME VALUE ERROR SIZE DERIVATIVE
	32	// 1 fPhiA 1.93913e-03 1.31059e-04 -0.00000e+00 -5.87280e-02
	33	// 2 fPhiA 2.48687e-04 5.41251e-05 0.00000e+00 -1.87361e-01
	34	// 3 fPhiAlpha -2.22649e+00 1.44503e-01 0.00000e+00 8.97538e-06
	35
	36	#include <TH1.h>
	37	#include <TH3.h>
	38	#include <TF1.h>
	39	#include <TRandom.h>
78d7c6d3	40	#include <AliAODParticle.h>
	41
	42
83b33650	43	ClassImp(AliHBTCorrectQInvCorrelFctn)
	44
	45	AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name,const char* title):
	46	AliHBTOnePairFctn1D(name,title),
	47	fMeasCorrelFctn(0x0),
	48	fMeasNumer(0x0),
	49	fMeasDenom(0x0),
	50	fSmearedNumer(0x0),
	51	fSmearedDenom(0x0),
	52	fDPtOverPtRMS(0.004),
	53	fThetaA(2.72e-03),
	54	fThetaB(1.87e-04),
	55	fThetaAlpha(-2.4),
	56	fPhiA(1.94e-03),
	57	fPhiB(2.5e-04),
	58	fPhiAlpha(-2.2),
	59	fR2(0.0),
	60	fLambda(0.0),
	61	fRConvergenceTreshold(0.3),
	62	fLambdaConvergenceTreshold(0.05)
	63	{
	64
	65	}
	66	/******************************************************************/
	67
	68	AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(TH1D* measqinv,const char* name,const char* title):
	69	AliHBTOnePairFctn1D(name,title,
	70	measqinv->GetNbinsX(),
	71	measqinv->GetXaxis()->GetXmax(),
	72	measqinv->GetXaxis()->GetXmin()),
	73	fMeasCorrelFctn(measqinv),
	74	fMeasNumer(0x0),
	75	fMeasDenom(0x0),
	76	fSmearedNumer(0x0),
	77	fSmearedDenom(0x0),
	78	fDPtOverPtRMS(0.004),
	79	fThetaA(2.72e-03),
	80	fThetaB(1.87e-04),
	81	fThetaAlpha(-2.4),
	82	fPhiA(1.94e-03),
	83	fPhiB(2.5e-04),
	84	fPhiAlpha(-2.2),
	85	fR2(0.0),
	86	fLambda(0.0),
	87	fRConvergenceTreshold(0.3),
	88	fLambdaConvergenceTreshold(0.05)
	89	{
	90
	91	}
	92	/******************************************************************/
	93
	94	AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name, const char* title,
	95	Int_t nbins, Float_t maxXval, Float_t minXval):
	96	AliHBTOnePairFctn1D(name,title,nbins,maxXval,minXval),
	97	fMeasCorrelFctn(0x0),
	98	fMeasNumer(0x0),
	99	fMeasDenom(0x0),
	100	fSmearedNumer(0x0),
	101	fSmearedDenom(0x0),
	102	fDPtOverPtRMS(0.004),
	103	fThetaA(2.72e-03),
	104	fThetaB(1.87e-04),
	105	fThetaAlpha(-2.4),
	106	fPhiA(1.94e-03),
107	fPhiB(2.5e-04),
108	fPhiAlpha(-2.2),
109	fR2(0.0),
110	fLambda(0.0),
111	fRConvergenceTreshold(0.3),
112	fLambdaConvergenceTreshold(0.05)
113	{
e09fa876	114	//ctor
83b33650	115	}
83b33650	116	/******************************************************************/
e09fa876	117	AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const AliHBTCorrectQInvCorrelFctn& in):
	118	AliHBTOnePairFctn1D(in),
	119	fMeasCorrelFctn(0x0),
	120	fMeasNumer(0x0),
	121	fMeasDenom(0x0),
	122	fSmearedNumer(0x0),
	123	fSmearedDenom(0x0),
	124	fDPtOverPtRMS(0),
	125	fThetaA(0),
	126	fThetaB(0),
	127	fThetaAlpha(0),
	128	fPhiA(0),
	129	fPhiB(0),
	130	fPhiAlpha(0),
	131	fR2(0.0),
	132	fLambda(0.0),
	133	fRConvergenceTreshold(0),
	134	fLambdaConvergenceTreshold(0)
	135	{
	136	//cpy ;ctor
	137	in.Copy(*this);
	138	}
83b33650	139
	140	AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
	141	{
e09fa876	142	//dtor
83b33650	143	delete fMeasCorrelFctn;
	144	delete fSmearedNumer;
	145	delete fSmearedDenom;
	146	delete fMeasNumer;
	147	delete fMeasDenom;
	148	}
	149	/******************************************************************/
	150
	151	void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min)
	152	{
	153	AliHBTFunction1D::BuildHistos(nbins,max,min);
	154
	155	TString numstr = fName + " Smeared Numerator"; //title and name of the numerator histogram
	156	TString denstr = fName + " Smeared Denominator";//title and name of the denominator histogram
	157
	158	fSmearedNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
	159	fSmearedDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
	160	fSmearedNumer->Sumw2();
	161	fSmearedDenom->Sumw2();
	162
	163	if (fMeasCorrelFctn == 0x0)
	164	{
	165	numstr = fName + " Measured Numerator"; //title and name of the numerator histogram
	166	denstr = fName + " Measured Denominator";//title and name of the denominator histogram
	167
	168	fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
	169	fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
	170	fMeasNumer->Sumw2();
	171	fMeasDenom->Sumw2();
	172	}
	173	}
	174	/******************************************************************/
	175
	176	void AliHBTCorrectQInvCorrelFctn::Init()
	177	{
e09fa876	178	//Init
83b33650	179	AliHBTOnePairFctn1D::Init();
	180	Info("Init","");
	181	fSmearedNumer->Reset();
	182	fSmearedDenom->Reset();
	183	if (fMeasNumer) fMeasNumer->Reset();
	184	if (fMeasDenom) fMeasDenom->Reset();
	185	fFittedR = -1.0;
	186	fFittedLambda = 0.0;
	187	}
	188	/******************************************************************/
	189
	190	void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair)
	191	{
e09fa876	192	//Processes particles that originates from the same event
83b33650	193	if (fMeasNumer == 0x0) return;
	194	pair = CheckPair(pair);
	195	if( pair == 0x0) return;
	196	fMeasNumer->Fill(pair->GetQInv());
	197	}
	198	/******************************************************************/
	199
	200	void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair)
	201	{
	202	//Process different events
78d7c6d3	203	static AliAODParticle part1, part2;
83b33650	204	static AliHBTPair smearedpair(&part1,&part2);
	205
	206	pair = CheckPair(pair);
	207	if( pair == 0x0) return;
	208	Double_t cc = GetCoulombCorrection(pair);
	209	Double_t qinv = pair->GetQInv();
	210
	211	//measured histogram -> if we are interested
	212	//only if fMeasCorrelFctn is not specified by user
	213	if (fMeasDenom) fMeasDenom->Fill(qinv,cc);
	214
	215	Smear(pair,smearedpair);
	216	Double_t modelqinv = GetModelValue(qinv);
	217	//Ideal histogram
	218	fNumerator->Fill(qinv,modelqinv*cc);
	219	fDenominator->Fill(qinv,cc);
	220
	221	//Smeared histogram
	222	Double_t smearedqinv = smearedpair.GetQInv();
	223	fSmearedNumer->Fill(smearedqinv,modelqinv);
	224	Double_t smearedcc = GetCoulombCorrection(&smearedpair);
	225	fSmearedDenom->Fill(smearedqinv,smearedcc);
	226
	227	}
	228	/******************************************************************/
	229	void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTPair* pair,AliHBTPair& smeared)
	230	{
	231	//Smears pair
	232	Smear(pair->Particle1(),smeared.Particle1());
	233	Smear(pair->Particle2(),smeared.Particle2());
	234	smeared.Changed();
	235	}
	236	/******************************************************************/
	237
78d7c6d3	238	void AliHBTCorrectQInvCorrelFctn::Smear(AliVAODParticle* part, AliVAODParticle* smeared)
83b33650	239	{
e09fa876	240	//Smears momenta
83b33650	241	Double_t sin2theta = TMath::Sin(part->Theta());
	242	sin2theta = sin2theta*sin2theta;
	243	Double_t pt = part->Pt();
	244
e09fa876	245	double dPtDivPt = gRandom->Gaus(0.0,fDPtOverPtRMS);
	246	double dphi = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fPhiAlpha));
	247	double dtheta = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fThetaAlpha));
83b33650	248
e09fa876	249	Double_t smearedPx = part->Px()(1.0+dPtDivPt) - part->Py()dphi;
	250	// fourmom.setX(px(1.0+dPtDivPt) - pydphi);
	251	Double_t smearedPy = part->Py()(1.0+dPtDivPt) - part->Px()dphi;
	252	// fourmom.setY(py(1.0+dPtDivPt) + pxdphi);
	253	Double_t smearedPz = part->Pz()(1.0+dPtDivPt) - ptdtheta/sin2theta;
	254	// fourmom.setZ(pz(1.0+dPtDivPt) - pTdtheta/sin2theta);
83b33650	255
78d7c6d3	256	Double_t mass2 = part->Mass()*part->Mass();
e09fa876	257	Double_t e = mass2 + smearedPx*smearedPx +
83b33650	258	smearedPy*smearedPy +
	259	smearedPz*smearedPz;
	260
e09fa876	261	smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(e));
83b33650	262	}
	263	/******************************************************************/
	264
	265	void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
	266	{
e09fa876	267	//Sets Initial Values
83b33650	268	fLambda = lambda;
	269	fR2 = r*r;
	270	}
	271	/******************************************************************/
	272	void AliHBTCorrectQInvCorrelFctn::MakeMeasCF()
	273	{
	274	//makes measured correlation function
	275	delete fMeasCorrelFctn;
	276	fMeasCorrelFctn = 0x0;
	277
	278	if (fMeasNumer&&fMeasDenom)
	279	{
	280	Double_t measscale = Scale(fMeasNumer,fMeasDenom);
	281	if (measscale == 0.0)
	282	{
	283	Error("GetResult","GetRatio for measured CF returned 0.0");
	284	return;
	285	}
	286	TString str = fName + "measured ratio";
	287	fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data());
	288	fMeasCorrelFctn->SetTitle(str.Data());
	289	fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale);
	290	}
	291	}
	292
	293	TH1* AliHBTCorrectQInvCorrelFctn::GetResult()
	294	{
	295	//In case we don't have yet Measured Correlation Function
e09fa876	296	//Try to get it
83b33650	297	//result is
	298	// N[meas] N[ideal]/D[ideal]
	299	//C(Q) = ------- * -----------------
	300	// D[meas] N[smear]/D[smear]
	301
	302	TString str;
	303	if (fMeasCorrelFctn == 0x0) MakeMeasCF();
	304
	305	if (fMeasCorrelFctn == 0x0)
	306	{
	307	Error("GetResult",
	308	"Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null");
	309	return 0x0;
	310	}
	311
	312	TH1D* ideal = (TH1D*)GetRatio(Scale());
	313	if (ideal == 0x0)
	314	{
	315	Error("GetResult","Ratio of ideal histograms is null");
	316	return 0x0;
	317	}
	318	str = fName + " smeared ratio";
	319	TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
	320	smearedCF->SetTitle(str.Data());
	321	Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
	322	smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
	323
	324	str = fName + " product meas ideal CF";
	325	TH1D* measideal = (TH1D*)ideal->Clone(str.Data());
	326	measideal->Multiply(ideal,fMeasCorrelFctn);
	327
	328	str = fName + " Corrected Result";
	329	TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data());
	330	result->SetTitle(str.Data());
	331
	332	Double_t resultscale = Scale(measideal,smearedCF);
	333	result->Divide(measideal,smearedCF,resultscale);
	334	return result;
	335
	336	}
	337	/******************************************************************/
	338
	339	void AliHBTCorrectQInvCorrelFctn::Fit()
	340	{
	341	//fits resuting histogram with function 1.0 + [0]exp([1][1]xx/(-0.038936366329))
	342	//where [0] is lambda
	343	// [1] is radius
	344	// 0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI
	345
	346	Info("Fit","Before fFittedLambda = %f",fFittedLambda);
	347	Info("Fit","Before fFittedR = %f",fFittedR);
	348	TH1D* result = (TH1D*)GetResult();
	349	if (result == 0x0)
	350	{
	351	Error("Fit","Can not get result");
	352	return;
	353	}
	354	TF1* fitfctn = new TF1("fitfctn","1.0 + [0]exp([1][1]xx/(-0.038936366329))");
	355	fitfctn->SetParameter(0,1.0);
	356	fitfctn->SetParameter(1,6.0);
	357	Float_t max = result->GetXaxis()->GetXmax();
	358	Info("Fit","Max is %f",max);
	359	result->Fit(fitfctn,"0","",0.008,max);
	360	fFittedLambda = fitfctn->GetParameter(0);
361	fFittedR = fitfctn->GetParameter(1);
362	Info("Fit","After fFittedLambda = %f",fFittedLambda);
363	Info("Fit","After fFittedR = %f",fFittedR);
364	delete fitfctn;
365	delete result;
366	}
367	/******************************************************************/
368
369	Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged()
370	{
371	//check if fitting was performed
372	if (fFittedR <= 0.0)
373	{
374	Fit();//if not do fit first
375	if (fFittedR <= 0.0)
376	{
377	Error("IsConverged","Fitting failed");
378	return kFALSE;
379	}
380	}
381
382	Double_t guessedR = TMath::Sqrt(fR2);
383
384	Info("IsConverged","Fitted lambda : %8f Fitted Radius : %8f",fFittedLambda,fFittedR);
385	Info("IsConverged","Guessed lambda : %8f Guessed Radius : %8f",fLambda,guessedR);
386	Info("IsConverged","Demanded lambda convergence: %8f Demanded Radius convergence: %8f",
387	fLambdaConvergenceTreshold,fRConvergenceTreshold);
388
389	if ( (TMath::Abs(fLambda-fFittedLambda)<fLambdaConvergenceTreshold) &&
390	(TMath::Abs(fFittedR-guessedR)<fRConvergenceTreshold) )
391	{
392	Info("IsConverged","Cnvergence reached");
393	return kTRUE;
394	}
395	else
396	{
397	Info("IsConverged","Cnvergence NOT reached");
398	return kFALSE;
399	}
400	}
401	/******************************************************************/
402
403	Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius()
404	{
e09fa876	405	//Returns Fitted radius
83b33650	406	if (fFittedR <= 0.0) Fit();
	407	return fFittedR;
	408	}
	409	/******************************************************************/
	410
	411	Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
	412	{
e09fa876	413	//Returns Fitted Intercept paramter
83b33650	414	if (fFittedR <= 0.0) Fit();
	415	return fFittedLambda;
	416	}
	417	/******************************************************************/
	418
	419	void AliHBTCorrectQInvCorrelFctn::WriteAll()
	420	{
e09fa876	421	//Writes function and all additional information
83b33650	422	Write();
	423	if (fMeasCorrelFctn) fMeasCorrelFctn->Write();
	424	if (fMeasNumer ) fMeasNumer->Write();
	425	if (fMeasDenom ) fMeasDenom->Write();
	426	if (fSmearedNumer) fSmearedNumer->Write();
	427	if (fSmearedDenom) fSmearedDenom->Write();
	428	if (fSmearedNumer && fSmearedDenom)
	429	{
	430	TString str = fName + " smeared ratio";
	431	TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
	432	smearedCF->SetTitle(str.Data());
	433	Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
	434	smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
	435	}
	436	}
	437