[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>
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 AliHBTParticle part1, part2;
  static AliHBTPair smearedpair(&part1,&part2);
  
  pair = CheckPair(pair);
  if( pair == 0x0) return;
  Double_t cc = GetCoulombCorrection(pair);
  Double_t qinv = pair->GetQInv();

  //measured histogram -> if we are interested 
  //only if fMeasCorrelFctn is not specified by user
  if (fMeasDenom) fMeasDenom->Fill(qinv,cc);

  Smear(pair,smearedpair);
  Double_t modelqinv = GetModelValue(qinv);  
  //Ideal histogram
  fNumerator->Fill(qinv,modelqinv*cc);
  fDenominator->Fill(qinv,cc);
  
  //Smeared histogram
  Double_t smearedqinv = smearedpair.GetQInv();
  fSmearedNumer->Fill(smearedqinv,modelqinv);
  Double_t smearedcc = GetCoulombCorrection(&smearedpair);
  fSmearedDenom->Fill(smearedqinv,smearedcc);
  
}
/******************************************************************/
void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTPair* pair,AliHBTPair& smeared)
{
//Smears pair
  Smear(pair->Particle1(),smeared.Particle1());
  Smear(pair->Particle2(),smeared.Particle2());
  smeared.Changed();
}
/******************************************************************/

void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTParticle* part, AliHBTParticle* smeared)
{
 //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->GetMass()*part->GetMass();
  Double_t e = mass2 + smearedPx*smearedPx + 
                       smearedPy*smearedPy + 
                       smearedPz*smearedPz;
	   
  smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(e));
}
/******************************************************************/

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