Updated to the new way AOD stores particles
[u/mrichter/AliRoot.git] / HBTAN / AliHBTCorrectQInvCorrelFctn.cxx
1 #include "AliHBTCorrectQInvCorrelFctn.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 #include <AliAODParticle.h>
41
42
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 {
114 //ctor
115 }                   
116 /******************************************************************/
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 }
139
140 AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
141 {
142 //dtor
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 {
178 //Init
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 {
192  //Processes particles that originates from the same event
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 
203   static AliAODParticle part1, part2;
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
238 void AliHBTCorrectQInvCorrelFctn::Smear(AliVAODParticle* part, AliVAODParticle* smeared)
239 {
240  //Smears momenta
241   Double_t sin2theta = TMath::Sin(part->Theta());
242   sin2theta = sin2theta*sin2theta;
243   Double_t pt = part->Pt();
244
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));
248   
249   Double_t smearedPx = part->Px()*(1.0+dPtDivPt) - part->Py()*dphi;
250 //  fourmom.setX(px*(1.0+dPtDivPt) - py*dphi);
251   Double_t smearedPy = part->Py()*(1.0+dPtDivPt) - part->Px()*dphi;
252 //  fourmom.setY(py*(1.0+dPtDivPt) + px*dphi);
253   Double_t smearedPz = part->Pz()*(1.0+dPtDivPt) - pt*dtheta/sin2theta;
254 //  fourmom.setZ(pz*(1.0+dPtDivPt) - pT*dtheta/sin2theta);
255   
256   Double_t mass2 = part->Mass()*part->Mass();
257   Double_t e = mass2 + smearedPx*smearedPx + 
258                        smearedPy*smearedPy + 
259                        smearedPz*smearedPz;
260            
261   smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(e));
262 }
263 /******************************************************************/
264
265 void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
266 {
267  //Sets Initial Values
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
296   //Try to get it 
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]*x*x/(-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]*x*x/(-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 {
405  //Returns Fitted radius
406   if (fFittedR <= 0.0) Fit();
407   return fFittedR;
408 }
409 /******************************************************************/
410
411 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
412 {
413  //Returns Fitted Intercept paramter
414   if (fFittedR <= 0.0) Fit();
415   return fFittedLambda;
416 }
417 /******************************************************************/
418
419 void AliHBTCorrectQInvCorrelFctn::WriteAll()
420 {
421 //Writes function and all additional information
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