Incrementing class versions
[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://aliweb.cern.ch/people/skowron/analyzer    //
15 //                                                   //
16 ///////////////////////////////////////////////////////
17
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
54
55
56 #include <TH1.h>
57 #include <TH3.h>
58 #include <TF1.h>
59 #include <TRandom.h>
60 #include <AliAODParticle.h>
61
62
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 + fDPtOverPtB*TMath::Power(pt,fDPtOverPtAlpha) + fDPtOverPtC*pt;
100   Double_t dPtDivPt = gRandom->Gaus(0.0,sigmapt);
101   Double_t dphi = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fPhiAlpha) + fPhiC*pt);
102   Double_t dtheta = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fThetaAlpha) +fThetaC*pt);
103   
104   Double_t smearedPx = part->Px()*(1.0+dPtDivPt) - part->Py()*dphi;
105 //  fourmom.setX(px*(1.0+dPtDivPt) - py*dphi);
106   Double_t smearedPy = part->Py()*(1.0+dPtDivPt) - part->Px()*dphi;
107 //  fourmom.setY(py*(1.0+dPtDivPt) + px*dphi);
108   Double_t smearedPz = part->Pz()*(1.0+dPtDivPt) - pt*dtheta/sin2theta;
109 //  fourmom.setZ(pz*(1.0+dPtDivPt) - pT*dtheta/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
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),
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),
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),
164   AliHBTCorrectedCorrelFctn(),
165   fMeasCorrelFctn(0x0),
166   fMeasNumer(0x0),
167   fMeasDenom(0x0),
168   fSmearedNumer(0x0),
169   fSmearedDenom(0x0),
170   fR2(0.0),
171   fLambda(0.0),
172   fRConvergenceTreshold(0.3),
173   fLambdaConvergenceTreshold(0.05)
174 {
175 //ctor
176 }                   
177 /******************************************************************/
178 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const AliHBTCorrectQInvCorrelFctn& in):
179   AliHBTOnePairFctn1D(in),
180   AliHBTCorrectedCorrelFctn(),
181   fMeasCorrelFctn(0x0),
182   fMeasNumer(0x0),
183   fMeasDenom(0x0),
184   fSmearedNumer(0x0),
185   fSmearedDenom(0x0),
186   fR2(0.0),
187   fLambda(0.0),
188   fRConvergenceTreshold(0),
189   fLambdaConvergenceTreshold(0)
190 {
191 //cpy ;ctor
192  in.Copy(*this);
193 }
194
195 AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
196 {
197 //dtor
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 {
233 //Init
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 {
247  //Processes particles that originates from the same event
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 
258   static AliAODParticle part1, part2;
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 /******************************************************************/
284
285 void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
286 {
287  //Sets Initial Values
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
316   //Try to get it 
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]*x*x/(-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]*x*x/(-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 {
425  //Returns Fitted radius
426   if (fFittedR <= 0.0) Fit();
427   return fFittedR;
428 }
429 /******************************************************************/
430
431 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
432 {
433  //Returns Fitted Intercept paramter
434   if (fFittedR <= 0.0) Fit();
435   return fFittedLambda;
436 }
437 /******************************************************************/
438
439 void AliHBTCorrectQInvCorrelFctn::WriteAll()
440 {
441 //Writes function and all additional information
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