]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HBTAN/AliHBTCorrectCorrelFctn.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Using TMath instead of math.h
[u/mrichter/AliRoot.git] / HBTAN / AliHBTCorrectCorrelFctn.cxx
CommitLineData
83b33650 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>
40ClassImp(AliHBTCorrectQInvCorrelFctn)
41
42AliHBTCorrectQInvCorrelFctn::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
65AliHBTCorrectQInvCorrelFctn::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
91AliHBTCorrectQInvCorrelFctn::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
114AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
115{
116 delete fMeasCorrelFctn;
117 delete fSmearedNumer;
118 delete fSmearedDenom;
119 delete fMeasNumer;
120 delete fMeasDenom;
121}
122/******************************************************************/
123
124void 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
149void 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
162void 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
171void 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/******************************************************************/
200void 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
209void 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
235void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
236{
237 fLambda = lambda;
238 fR2 = r*r;
239}
240/******************************************************************/
241void 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
262TH1* 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
308void 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
338Bool_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
372Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius()
373{
374 if (fFittedR <= 0.0) Fit();
375 return fFittedR;
376}
377/******************************************************************/
378
379Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
380{
381 if (fFittedR <= 0.0) Fit();
382 return fFittedLambda;
383}
384/******************************************************************/
385
386void 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
420AliHBTCorrectQ3DCorrelFctn::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
432AliHBTCorrectQ3DCorrelFctn::~AliHBTCorrectQ3DCorrelFctn()
433{
434 delete fMeasCorrelFctn;
435 delete fSmearedNumer;
436 delete fSmearedDenom;
437 delete fMeasNumer;
438 delete fMeasDenom;
439}
440
441/******************************************************************/