]>
Commit | Line | Data |
---|---|---|
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 | 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 | ||
83b33650 | 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), | |
83b33650 | 133 | fR2(0.0), |
134 | fLambda(0.0), | |
4b1c9620 | 135 | fFittedR(0), |
136 | fFittedLambda(0), | |
83b33650 | 137 | fRConvergenceTreshold(0.3), |
138 | fLambdaConvergenceTreshold(0.05) | |
139 | { | |
140 | ||
141 | } | |
142 | /******************************************************************/ | |
143 | ||
144 | AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(TH1D* measqinv,const char* name,const char* title): | |
145 | AliHBTOnePairFctn1D(name,title, | |
146 | measqinv->GetNbinsX(), | |
147 | measqinv->GetXaxis()->GetXmax(), | |
148 | measqinv->GetXaxis()->GetXmin()), | |
149 | fMeasCorrelFctn(measqinv), | |
150 | fMeasNumer(0x0), | |
151 | fMeasDenom(0x0), | |
152 | fSmearedNumer(0x0), | |
153 | fSmearedDenom(0x0), | |
83b33650 | 154 | fR2(0.0), |
155 | fLambda(0.0), | |
4b1c9620 | 156 | fFittedR(0), |
157 | fFittedLambda(0), | |
83b33650 | 158 | fRConvergenceTreshold(0.3), |
159 | fLambdaConvergenceTreshold(0.05) | |
160 | { | |
161 | ||
162 | } | |
163 | /******************************************************************/ | |
164 | ||
165 | AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name, const char* title, | |
166 | Int_t nbins, Float_t maxXval, Float_t minXval): | |
167 | AliHBTOnePairFctn1D(name,title,nbins,maxXval,minXval), | |
62e1b4fe | 168 | AliHBTCorrectedCorrelFctn(), |
83b33650 | 169 | fMeasCorrelFctn(0x0), |
170 | fMeasNumer(0x0), | |
171 | fMeasDenom(0x0), | |
172 | fSmearedNumer(0x0), | |
173 | fSmearedDenom(0x0), | |
83b33650 | 174 | fR2(0.0), |
175 | fLambda(0.0), | |
4b1c9620 | 176 | fFittedR(0), |
177 | fFittedLambda(0), | |
83b33650 | 178 | fRConvergenceTreshold(0.3), |
179 | fLambdaConvergenceTreshold(0.05) | |
180 | { | |
e09fa876 | 181 | //ctor |
83b33650 | 182 | } |
183 | /******************************************************************/ | |
e09fa876 | 184 | AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const AliHBTCorrectQInvCorrelFctn& in): |
185 | AliHBTOnePairFctn1D(in), | |
62e1b4fe | 186 | AliHBTCorrectedCorrelFctn(), |
e09fa876 | 187 | fMeasCorrelFctn(0x0), |
188 | fMeasNumer(0x0), | |
189 | fMeasDenom(0x0), | |
190 | fSmearedNumer(0x0), | |
191 | fSmearedDenom(0x0), | |
e09fa876 | 192 | fR2(0.0), |
193 | fLambda(0.0), | |
4b1c9620 | 194 | fFittedR(0), |
195 | fFittedLambda(0), | |
e09fa876 | 196 | fRConvergenceTreshold(0), |
197 | fLambdaConvergenceTreshold(0) | |
198 | { | |
199 | //cpy ;ctor | |
200 | in.Copy(*this); | |
201 | } | |
83b33650 | 202 | |
203 | AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn() | |
204 | { | |
e09fa876 | 205 | //dtor |
83b33650 | 206 | delete fMeasCorrelFctn; |
207 | delete fSmearedNumer; | |
208 | delete fSmearedDenom; | |
209 | delete fMeasNumer; | |
210 | delete fMeasDenom; | |
211 | } | |
212 | /******************************************************************/ | |
213 | ||
214 | void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min) | |
215 | { | |
216 | AliHBTFunction1D::BuildHistos(nbins,max,min); | |
217 | ||
218 | TString numstr = fName + " Smeared Numerator"; //title and name of the numerator histogram | |
219 | TString denstr = fName + " Smeared Denominator";//title and name of the denominator histogram | |
220 | ||
221 | fSmearedNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max); | |
222 | fSmearedDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max); | |
223 | fSmearedNumer->Sumw2(); | |
224 | fSmearedDenom->Sumw2(); | |
225 | ||
226 | if (fMeasCorrelFctn == 0x0) | |
227 | { | |
228 | numstr = fName + " Measured Numerator"; //title and name of the numerator histogram | |
229 | denstr = fName + " Measured Denominator";//title and name of the denominator histogram | |
230 | ||
231 | fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max); | |
232 | fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max); | |
233 | fMeasNumer->Sumw2(); | |
234 | fMeasDenom->Sumw2(); | |
235 | } | |
236 | } | |
237 | /******************************************************************/ | |
238 | ||
239 | void AliHBTCorrectQInvCorrelFctn::Init() | |
240 | { | |
e09fa876 | 241 | //Init |
83b33650 | 242 | AliHBTOnePairFctn1D::Init(); |
243 | Info("Init",""); | |
244 | fSmearedNumer->Reset(); | |
245 | fSmearedDenom->Reset(); | |
246 | if (fMeasNumer) fMeasNumer->Reset(); | |
247 | if (fMeasDenom) fMeasDenom->Reset(); | |
248 | fFittedR = -1.0; | |
249 | fFittedLambda = 0.0; | |
250 | } | |
251 | /******************************************************************/ | |
252 | ||
253 | void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair) | |
254 | { | |
e09fa876 | 255 | //Processes particles that originates from the same event |
83b33650 | 256 | if (fMeasNumer == 0x0) return; |
257 | pair = CheckPair(pair); | |
258 | if( pair == 0x0) return; | |
259 | fMeasNumer->Fill(pair->GetQInv()); | |
260 | } | |
261 | /******************************************************************/ | |
262 | ||
263 | void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair) | |
264 | { | |
265 | //Process different events | |
78d7c6d3 | 266 | static AliAODParticle part1, part2; |
83b33650 | 267 | static AliHBTPair smearedpair(&part1,&part2); |
268 | ||
269 | pair = CheckPair(pair); | |
270 | if( pair == 0x0) return; | |
271 | Double_t cc = GetCoulombCorrection(pair); | |
272 | Double_t qinv = pair->GetQInv(); | |
273 | ||
274 | //measured histogram -> if we are interested | |
275 | //only if fMeasCorrelFctn is not specified by user | |
276 | if (fMeasDenom) fMeasDenom->Fill(qinv,cc); | |
277 | ||
278 | Smear(pair,smearedpair); | |
279 | Double_t modelqinv = GetModelValue(qinv); | |
280 | //Ideal histogram | |
281 | fNumerator->Fill(qinv,modelqinv*cc); | |
282 | fDenominator->Fill(qinv,cc); | |
283 | ||
284 | //Smeared histogram | |
285 | Double_t smearedqinv = smearedpair.GetQInv(); | |
286 | fSmearedNumer->Fill(smearedqinv,modelqinv); | |
287 | Double_t smearedcc = GetCoulombCorrection(&smearedpair); | |
288 | fSmearedDenom->Fill(smearedqinv,smearedcc); | |
289 | ||
290 | } | |
291 | /******************************************************************/ | |
83b33650 | 292 | |
293 | void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r) | |
294 | { | |
e09fa876 | 295 | //Sets Initial Values |
83b33650 | 296 | fLambda = lambda; |
297 | fR2 = r*r; | |
298 | } | |
299 | /******************************************************************/ | |
300 | void AliHBTCorrectQInvCorrelFctn::MakeMeasCF() | |
301 | { | |
302 | //makes measured correlation function | |
303 | delete fMeasCorrelFctn; | |
304 | fMeasCorrelFctn = 0x0; | |
305 | ||
306 | if (fMeasNumer&&fMeasDenom) | |
307 | { | |
308 | Double_t measscale = Scale(fMeasNumer,fMeasDenom); | |
309 | if (measscale == 0.0) | |
310 | { | |
311 | Error("GetResult","GetRatio for measured CF returned 0.0"); | |
312 | return; | |
313 | } | |
314 | TString str = fName + "measured ratio"; | |
315 | fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data()); | |
316 | fMeasCorrelFctn->SetTitle(str.Data()); | |
317 | fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale); | |
318 | } | |
319 | } | |
320 | ||
321 | TH1* AliHBTCorrectQInvCorrelFctn::GetResult() | |
322 | { | |
323 | //In case we don't have yet Measured Correlation Function | |
e09fa876 | 324 | //Try to get it |
83b33650 | 325 | //result is |
326 | // N[meas] N[ideal]/D[ideal] | |
327 | //C(Q) = ------- * ----------------- | |
328 | // D[meas] N[smear]/D[smear] | |
329 | ||
330 | TString str; | |
331 | if (fMeasCorrelFctn == 0x0) MakeMeasCF(); | |
332 | ||
333 | if (fMeasCorrelFctn == 0x0) | |
334 | { | |
335 | Error("GetResult", | |
336 | "Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null"); | |
337 | return 0x0; | |
338 | } | |
339 | ||
340 | TH1D* ideal = (TH1D*)GetRatio(Scale()); | |
341 | if (ideal == 0x0) | |
342 | { | |
343 | Error("GetResult","Ratio of ideal histograms is null"); | |
344 | return 0x0; | |
345 | } | |
346 | str = fName + " smeared ratio"; | |
347 | TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data()); | |
348 | smearedCF->SetTitle(str.Data()); | |
349 | Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom); | |
350 | smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale); | |
351 | ||
352 | str = fName + " product meas ideal CF"; | |
353 | TH1D* measideal = (TH1D*)ideal->Clone(str.Data()); | |
354 | measideal->Multiply(ideal,fMeasCorrelFctn); | |
355 | ||
356 | str = fName + " Corrected Result"; | |
357 | TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data()); | |
358 | result->SetTitle(str.Data()); | |
359 | ||
360 | Double_t resultscale = Scale(measideal,smearedCF); | |
361 | result->Divide(measideal,smearedCF,resultscale); | |
362 | return result; | |
363 | ||
364 | } | |
365 | /******************************************************************/ | |
366 | ||
367 | void AliHBTCorrectQInvCorrelFctn::Fit() | |
368 | { | |
369 | //fits resuting histogram with function 1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329)) | |
370 | //where [0] is lambda | |
371 | // [1] is radius | |
372 | // 0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI | |
373 | ||
374 | Info("Fit","Before fFittedLambda = %f",fFittedLambda); | |
375 | Info("Fit","Before fFittedR = %f",fFittedR); | |
376 | TH1D* result = (TH1D*)GetResult(); | |
377 | if (result == 0x0) | |
378 | { | |
379 | Error("Fit","Can not get result"); | |
380 | return; | |
381 | } | |
382 | TF1* fitfctn = new TF1("fitfctn","1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))"); | |
383 | fitfctn->SetParameter(0,1.0); | |
384 | fitfctn->SetParameter(1,6.0); | |
385 | Float_t max = result->GetXaxis()->GetXmax(); | |
386 | Info("Fit","Max is %f",max); | |
387 | result->Fit(fitfctn,"0","",0.008,max); | |
388 | fFittedLambda = fitfctn->GetParameter(0); | |
389 | fFittedR = fitfctn->GetParameter(1); | |
390 | Info("Fit","After fFittedLambda = %f",fFittedLambda); | |
391 | Info("Fit","After fFittedR = %f",fFittedR); | |
392 | delete fitfctn; | |
393 | delete result; | |
394 | } | |
395 | /******************************************************************/ | |
396 | ||
397 | Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged() | |
398 | { | |
399 | //check if fitting was performed | |
400 | if (fFittedR <= 0.0) | |
401 | { | |
402 | Fit();//if not do fit first | |
403 | if (fFittedR <= 0.0) | |
404 | { | |
405 | Error("IsConverged","Fitting failed"); | |
406 | return kFALSE; | |
407 | } | |
408 | } | |
409 | ||
410 | Double_t guessedR = TMath::Sqrt(fR2); | |
411 | ||
412 | Info("IsConverged","Fitted lambda : %8f Fitted Radius : %8f",fFittedLambda,fFittedR); | |
413 | Info("IsConverged","Guessed lambda : %8f Guessed Radius : %8f",fLambda,guessedR); | |
414 | Info("IsConverged","Demanded lambda convergence: %8f Demanded Radius convergence: %8f", | |
415 | fLambdaConvergenceTreshold,fRConvergenceTreshold); | |
416 | ||
417 | if ( (TMath::Abs(fLambda-fFittedLambda)<fLambdaConvergenceTreshold) && | |
418 | (TMath::Abs(fFittedR-guessedR)<fRConvergenceTreshold) ) | |
419 | { | |
420 | Info("IsConverged","Cnvergence reached"); | |
421 | return kTRUE; | |
422 | } | |
423 | else | |
424 | { | |
425 | Info("IsConverged","Cnvergence NOT reached"); | |
426 | return kFALSE; | |
427 | } | |
428 | } | |
429 | /******************************************************************/ | |
430 | ||
431 | Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius() | |
432 | { | |
e09fa876 | 433 | //Returns Fitted radius |
83b33650 | 434 | if (fFittedR <= 0.0) Fit(); |
435 | return fFittedR; | |
436 | } | |
437 | /******************************************************************/ | |
438 | ||
439 | Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda() | |
440 | { | |
e09fa876 | 441 | //Returns Fitted Intercept paramter |
83b33650 | 442 | if (fFittedR <= 0.0) Fit(); |
443 | return fFittedLambda; | |
444 | } | |
445 | /******************************************************************/ | |
446 | ||
447 | void AliHBTCorrectQInvCorrelFctn::WriteAll() | |
448 | { | |
e09fa876 | 449 | //Writes function and all additional information |
83b33650 | 450 | Write(); |
451 | if (fMeasCorrelFctn) fMeasCorrelFctn->Write(); | |
452 | if (fMeasNumer ) fMeasNumer->Write(); | |
453 | if (fMeasDenom ) fMeasDenom->Write(); | |
454 | if (fSmearedNumer) fSmearedNumer->Write(); | |
455 | if (fSmearedDenom) fSmearedDenom->Write(); | |
456 | if (fSmearedNumer && fSmearedDenom) | |
457 | { | |
458 | TString str = fName + " smeared ratio"; | |
459 | TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data()); | |
460 | smearedCF->SetTitle(str.Data()); | |
461 | Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom); | |
462 | smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale); | |
463 | } | |
464 | } | |
465 |