Classes for Corrections with Lisa algorithm added
[u/mrichter/AliRoot.git] / HBTAN / AliHBTCorrectCorrelFctn.cxx
1 #include "AliHBTCorrectCorrelFctn.h"
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 {
111 }                   
112 /******************************************************************/
113
114 AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
115 {
116   delete fMeasCorrelFctn;
117   delete fSmearedNumer;
118   delete fSmearedDenom;
119   delete fMeasNumer;
120   delete fMeasDenom;
121 }
122 /******************************************************************/
123
124 void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min)
125 {
126   AliHBTFunction1D::BuildHistos(nbins,max,min);
127   
128   TString numstr = fName + " Smeared Numerator";  //title and name of the numerator histogram
129   TString denstr = fName + " Smeared Denominator";//title and name of the denominator histogram
130   
131   fSmearedNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
132   fSmearedDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
133   fSmearedNumer->Sumw2();
134   fSmearedDenom->Sumw2();
135   
136   if (fMeasCorrelFctn == 0x0)
137    { 
138      numstr = fName + " Measured Numerator";  //title and name of the numerator histogram
139      denstr = fName + " Measured Denominator";//title and name of the denominator histogram
140     
141      fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
142      fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
143      fMeasNumer->Sumw2();
144      fMeasDenom->Sumw2();
145    }
146 }
147 /******************************************************************/
148
149 void AliHBTCorrectQInvCorrelFctn::Init()
150 {
151   AliHBTOnePairFctn1D::Init();
152   Info("Init","");
153   fSmearedNumer->Reset();
154   fSmearedDenom->Reset();
155   if (fMeasNumer) fMeasNumer->Reset();
156   if (fMeasDenom) fMeasDenom->Reset();
157   fFittedR = -1.0;
158   fFittedLambda = 0.0;
159 }
160 /******************************************************************/
161
162 void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair)
163 {
164   if (fMeasNumer == 0x0) return;
165   pair = CheckPair(pair);
166   if( pair == 0x0) return;
167   fMeasNumer->Fill(pair->GetQInv());
168 }
169 /******************************************************************/
170
171 void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair)
172 {
173 //Process different events 
174   static AliHBTParticle part1, part2;
175   static AliHBTPair smearedpair(&part1,&part2);
176   
177   pair = CheckPair(pair);
178   if( pair == 0x0) return;
179   Double_t cc = GetCoulombCorrection(pair);
180   Double_t qinv = pair->GetQInv();
181
182   //measured histogram -> if we are interested 
183   //only if fMeasCorrelFctn is not specified by user
184   if (fMeasDenom) fMeasDenom->Fill(qinv,cc);
185
186   Smear(pair,smearedpair);
187   Double_t modelqinv = GetModelValue(qinv);  
188   //Ideal histogram
189   fNumerator->Fill(qinv,modelqinv*cc);
190   fDenominator->Fill(qinv,cc);
191   
192   //Smeared histogram
193   Double_t smearedqinv = smearedpair.GetQInv();
194   fSmearedNumer->Fill(smearedqinv,modelqinv);
195   Double_t smearedcc = GetCoulombCorrection(&smearedpair);
196   fSmearedDenom->Fill(smearedqinv,smearedcc);
197   
198 }
199 /******************************************************************/
200 void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTPair* pair,AliHBTPair& smeared)
201 {
202 //Smears pair
203   Smear(pair->Particle1(),smeared.Particle1());
204   Smear(pair->Particle2(),smeared.Particle2());
205   smeared.Changed();
206 }
207 /******************************************************************/
208
209 void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTParticle* part, AliHBTParticle* smeared)
210 {
211   Double_t sin2theta = TMath::Sin(part->Theta());
212   sin2theta = sin2theta*sin2theta;
213   Double_t pt = part->Pt();
214
215   double DpT_div_pT = gRandom->Gaus(0.0,fDPtOverPtRMS);
216   double Dphi = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fPhiAlpha));
217   double Dtheta = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fThetaAlpha));
218   
219   Double_t smearedPx = part->Px()*(1.0+DpT_div_pT) - part->Py()*Dphi;
220 //  fourmom.setX(px*(1.0+DpT_div_pT) - py*Dphi);
221   Double_t smearedPy = part->Py()*(1.0+DpT_div_pT) - part->Px()*Dphi;
222 //  fourmom.setY(py*(1.0+DpT_div_pT) + px*Dphi);
223   Double_t smearedPz = part->Pz()*(1.0+DpT_div_pT) - pt*Dtheta/sin2theta;
224 //  fourmom.setZ(pz*(1.0+DpT_div_pT) - pT*Dtheta/sin2theta);
225   
226   Double_t mass2 = part->GetMass()*part->GetMass();
227   Double_t E = mass2 + smearedPx*smearedPx + 
228                        smearedPy*smearedPy + 
229                        smearedPz*smearedPz;
230            
231   smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(E));
232 }
233 /******************************************************************/
234
235 void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
236 {
237   fLambda = lambda;
238   fR2 = r*r;
239 }
240 /******************************************************************/
241 void AliHBTCorrectQInvCorrelFctn::MakeMeasCF()
242 {
243   //makes measured correlation function
244   delete fMeasCorrelFctn;
245   fMeasCorrelFctn = 0x0;
246   
247   if (fMeasNumer&&fMeasDenom)
248    {
249      Double_t measscale = Scale(fMeasNumer,fMeasDenom);
250      if (measscale == 0.0)
251       {
252         Error("GetResult","GetRatio for measured CF returned 0.0");
253         return;
254       }
255      TString str = fName + "measured ratio";
256      fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data());
257      fMeasCorrelFctn->SetTitle(str.Data());
258      fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale);
259    }
260 }
261
262 TH1* AliHBTCorrectQInvCorrelFctn::GetResult()
263 {
264   //In case we don't have yet Measured Correlation Function
265   //Try to get it
266   //result is
267   //        N[meas]   N[ideal]/D[ideal]
268   //C(Q) =  ------- * -----------------
269   //        D[meas]   N[smear]/D[smear]
270  
271   TString str;
272   if (fMeasCorrelFctn == 0x0) MakeMeasCF();
273   
274   if (fMeasCorrelFctn == 0x0)
275    {
276      Error("GetResult",  
277            "Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null");
278      return 0x0;
279    }
280   
281   TH1D* ideal = (TH1D*)GetRatio(Scale());
282   if (ideal == 0x0)
283    {
284      Error("GetResult","Ratio of ideal histograms is null");
285      return 0x0;
286    }
287   str = fName + " smeared ratio";
288   TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
289   smearedCF->SetTitle(str.Data());
290   Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
291   smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
292   
293   str = fName + " product meas ideal CF";
294   TH1D* measideal = (TH1D*)ideal->Clone(str.Data());
295   measideal->Multiply(ideal,fMeasCorrelFctn);
296   
297   str = fName + " Corrected Result";
298   TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data());
299   result->SetTitle(str.Data());
300   
301   Double_t resultscale = Scale(measideal,smearedCF);
302   result->Divide(measideal,smearedCF,resultscale);
303   return result;
304   
305 }
306 /******************************************************************/
307
308 void AliHBTCorrectQInvCorrelFctn::Fit()
309 {
310 //fits resuting histogram with function 1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))
311 //where [0] is lambda
312 //      [1] is radius
313 //   0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI
314
315   Info("Fit","Before fFittedLambda = %f",fFittedLambda);
316   Info("Fit","Before fFittedR = %f",fFittedR);
317   TH1D* result = (TH1D*)GetResult();
318   if (result == 0x0)
319    {
320      Error("Fit","Can not get result");
321      return;
322    }
323   TF1* fitfctn = new TF1("fitfctn","1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))");
324   fitfctn->SetParameter(0,1.0);
325   fitfctn->SetParameter(1,6.0);
326   Float_t max = result->GetXaxis()->GetXmax();
327   Info("Fit","Max is %f",max);
328   result->Fit(fitfctn,"0","",0.008,max);
329   fFittedLambda = fitfctn->GetParameter(0);
330   fFittedR = fitfctn->GetParameter(1);
331   Info("Fit","After fFittedLambda = %f",fFittedLambda);
332   Info("Fit","After fFittedR = %f",fFittedR);
333   delete fitfctn;
334   delete result;
335 }
336 /******************************************************************/
337
338 Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged()
339 {
340   //check if fitting was performed
341   if (fFittedR <= 0.0)
342    {
343      Fit();//if not do fit first
344      if (fFittedR <= 0.0)
345       {
346         Error("IsConverged","Fitting failed");
347         return kFALSE;
348       }
349    }
350
351   Double_t guessedR = TMath::Sqrt(fR2);
352   
353   Info("IsConverged","Fitted   lambda            : %8f Fitted  Radius             : %8f",fFittedLambda,fFittedR);
354   Info("IsConverged","Guessed  lambda            : %8f Guessed Radius             : %8f",fLambda,guessedR);
355   Info("IsConverged","Demanded lambda convergence: %8f Demanded Radius convergence: %8f",
356        fLambdaConvergenceTreshold,fRConvergenceTreshold);
357   
358   if ( (TMath::Abs(fLambda-fFittedLambda)<fLambdaConvergenceTreshold) && 
359        (TMath::Abs(fFittedR-guessedR)<fRConvergenceTreshold) )
360     {
361       Info("IsConverged","Cnvergence reached");
362       return kTRUE;
363     }
364   else
365    {
366       Info("IsConverged","Cnvergence NOT reached");
367       return kFALSE;
368    }
369 }
370 /******************************************************************/
371
372 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius()
373 {
374   if (fFittedR <= 0.0) Fit();
375   return fFittedR;
376 }
377 /******************************************************************/
378
379 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
380 {
381   if (fFittedR <= 0.0) Fit();
382   return fFittedLambda;
383 }
384 /******************************************************************/
385
386 void AliHBTCorrectQInvCorrelFctn::WriteAll()
387 {
388   Write();
389   if (fMeasCorrelFctn) fMeasCorrelFctn->Write();
390   if (fMeasNumer ) fMeasNumer->Write();
391   if (fMeasDenom ) fMeasDenom->Write();
392   if (fSmearedNumer) fSmearedNumer->Write();
393   if (fSmearedDenom) fSmearedDenom->Write();
394   if (fSmearedNumer && fSmearedDenom)
395    {
396     TString str = fName + " smeared ratio";
397     TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
398     smearedCF->SetTitle(str.Data());
399     Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
400     smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
401    }
402 }
403
404 /******************************************************************/
405 /******************************************************************/
406 /******************************************************************/
407
408 //____________________
409 ///////////////////////////////////////////////////////
410 //                                                   //
411 // AliHBTCorrectQ3DCorrelFctn                        //
412 //                                                   //
413 // Class for calculating Q Invariant correlation     //
414 // taking to the account resolution of the           //
415 // detector and coulomb effects.                     //
416 //                                                   //
417 ///////////////////////////////////////////////////////
418
419
420 AliHBTCorrectQ3DCorrelFctn::AliHBTCorrectQ3DCorrelFctn(const char* name, const char* title):
421  AliHBTOnePairFctn3D(name,title),
422  fMeasCorrelFctn(0x0),
423  fSmearedNumer(0x0),
424  fSmearedDenom(0x0),
425  fMeasNumer(0x0),
426  fMeasDenom(0x0)
427 {
428
429 }
430 /******************************************************************/
431
432 AliHBTCorrectQ3DCorrelFctn::~AliHBTCorrectQ3DCorrelFctn()
433 {
434  delete fMeasCorrelFctn;
435  delete fSmearedNumer;
436  delete fSmearedDenom;
437  delete fMeasNumer;
438  delete fMeasDenom;
439 }
440
441 /******************************************************************/