]>
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), | |
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), | |
83b33650 | 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), | |
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 | 178 | AliHBTCorrectQInvCorrelFctn::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 | |
195 | AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn() | |
196 | { | |
e09fa876 | 197 | //dtor |
83b33650 | 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 | { | |
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 | ||
245 | void 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 | ||
255 | void 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 | |
285 | void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r) | |
286 | { | |
e09fa876 | 287 | //Sets Initial Values |
83b33650 | 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 | |
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 | ||
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 | { | |
e09fa876 | 425 | //Returns Fitted radius |
83b33650 | 426 | if (fFittedR <= 0.0) Fit(); |
427 | return fFittedR; | |
428 | } | |
429 | /******************************************************************/ | |
430 | ||
431 | Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda() | |
432 | { | |
e09fa876 | 433 | //Returns Fitted Intercept paramter |
83b33650 | 434 | if (fFittedR <= 0.0) Fit(); |
435 | return fFittedLambda; | |
436 | } | |
437 | /******************************************************************/ | |
438 | ||
439 | void 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 |