]>
Commit | Line | Data |
---|---|---|
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> | |
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 | /******************************************************************/ |