Removing dependences on AliDAQ class (in raw)
[u/mrichter/AliRoot.git] / HBTAN / AliHBTCorrectQInvCorrelFctn.cxx
CommitLineData
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 //
c7ffd78f 14// http://aliweb.cern.ch/people/skowron/analyzer //
83b33650 15// //
16///////////////////////////////////////////////////////
17
62e1b4fe 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
83b33650 54
83b33650 55
56#include <TH1.h>
57#include <TH3.h>
58#include <TF1.h>
59#include <TRandom.h>
78d7c6d3 60#include <AliAODParticle.h>
61
62
62e1b4fe 63ClassImp(AliHBTCorrectedCorrelFctn)
64
65AliHBTCorrectedCorrelFctn::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/******************************************************************/
83void 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
92void 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
83b33650 124ClassImp(AliHBTCorrectQInvCorrelFctn)
125
126AliHBTCorrectQInvCorrelFctn::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),
83b33650 133 fR2(0.0),
134 fLambda(0.0),
135 fRConvergenceTreshold(0.3),
136 fLambdaConvergenceTreshold(0.05)
137{
138
139}
140/******************************************************************/
141
142AliHBTCorrectQInvCorrelFctn::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),
83b33650 152 fR2(0.0),
153 fLambda(0.0),
154 fRConvergenceTreshold(0.3),
155 fLambdaConvergenceTreshold(0.05)
156{
157
158}
159/******************************************************************/
160
161AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name, const char* title,
162 Int_t nbins, Float_t maxXval, Float_t minXval):
163 AliHBTOnePairFctn1D(name,title,nbins,maxXval,minXval),
62e1b4fe 164 AliHBTCorrectedCorrelFctn(),
83b33650 165 fMeasCorrelFctn(0x0),
166 fMeasNumer(0x0),
167 fMeasDenom(0x0),
168 fSmearedNumer(0x0),
169 fSmearedDenom(0x0),
83b33650 170 fR2(0.0),
171 fLambda(0.0),
172 fRConvergenceTreshold(0.3),
173 fLambdaConvergenceTreshold(0.05)
174{
e09fa876 175//ctor
83b33650 176}
177/******************************************************************/
e09fa876 178AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const AliHBTCorrectQInvCorrelFctn& in):
179 AliHBTOnePairFctn1D(in),
62e1b4fe 180 AliHBTCorrectedCorrelFctn(),
e09fa876 181 fMeasCorrelFctn(0x0),
182 fMeasNumer(0x0),
183 fMeasDenom(0x0),
184 fSmearedNumer(0x0),
185 fSmearedDenom(0x0),
e09fa876 186 fR2(0.0),
187 fLambda(0.0),
188 fRConvergenceTreshold(0),
189 fLambdaConvergenceTreshold(0)
190{
191//cpy ;ctor
192 in.Copy(*this);
193}
83b33650 194
195AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
196{
e09fa876 197//dtor
83b33650 198 delete fMeasCorrelFctn;
199 delete fSmearedNumer;
200 delete fSmearedDenom;
201 delete fMeasNumer;
202 delete fMeasDenom;
203}
204/******************************************************************/
205
206void 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
231void AliHBTCorrectQInvCorrelFctn::Init()
232{
e09fa876 233//Init
83b33650 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
245void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair)
246{
e09fa876 247 //Processes particles that originates from the same event
83b33650 248 if (fMeasNumer == 0x0) return;
249 pair = CheckPair(pair);
250 if( pair == 0x0) return;
251 fMeasNumer->Fill(pair->GetQInv());
252}
253/******************************************************************/
254
255void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair)
256{
257//Process different events
78d7c6d3 258 static AliAODParticle part1, part2;
83b33650 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/******************************************************************/
83b33650 284
285void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
286{
e09fa876 287 //Sets Initial Values
83b33650 288 fLambda = lambda;
289 fR2 = r*r;
290}
291/******************************************************************/
292void 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
313TH1* AliHBTCorrectQInvCorrelFctn::GetResult()
314{
315 //In case we don't have yet Measured Correlation Function
e09fa876 316 //Try to get it
83b33650 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
359void 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
389Bool_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
423Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius()
424{
e09fa876 425 //Returns Fitted radius
83b33650 426 if (fFittedR <= 0.0) Fit();
427 return fFittedR;
428}
429/******************************************************************/
430
431Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
432{
e09fa876 433 //Returns Fitted Intercept paramter
83b33650 434 if (fFittedR <= 0.0) Fit();
435 return fFittedLambda;
436}
437/******************************************************************/
438
439void AliHBTCorrectQInvCorrelFctn::WriteAll()
440{
e09fa876 441//Writes function and all additional information
83b33650 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