]>
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 // | |
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> | |
78d7c6d3 | 40 | #include <AliAODParticle.h> |
41 | ||
42 | ||
83b33650 | 43 | ClassImp(AliHBTCorrectQInvCorrelFctn) |
44 | ||
45 | AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name,const char* title): | |
46 | AliHBTOnePairFctn1D(name,title), | |
47 | fMeasCorrelFctn(0x0), | |
48 | fMeasNumer(0x0), | |
49 | fMeasDenom(0x0), | |
50 | fSmearedNumer(0x0), | |
51 | fSmearedDenom(0x0), | |
52 | fDPtOverPtRMS(0.004), | |
53 | fThetaA(2.72e-03), | |
54 | fThetaB(1.87e-04), | |
55 | fThetaAlpha(-2.4), | |
56 | fPhiA(1.94e-03), | |
57 | fPhiB(2.5e-04), | |
58 | fPhiAlpha(-2.2), | |
59 | fR2(0.0), | |
60 | fLambda(0.0), | |
61 | fRConvergenceTreshold(0.3), | |
62 | fLambdaConvergenceTreshold(0.05) | |
63 | { | |
64 | ||
65 | } | |
66 | /******************************************************************/ | |
67 | ||
68 | AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(TH1D* measqinv,const char* name,const char* title): | |
69 | AliHBTOnePairFctn1D(name,title, | |
70 | measqinv->GetNbinsX(), | |
71 | measqinv->GetXaxis()->GetXmax(), | |
72 | measqinv->GetXaxis()->GetXmin()), | |
73 | fMeasCorrelFctn(measqinv), | |
74 | fMeasNumer(0x0), | |
75 | fMeasDenom(0x0), | |
76 | fSmearedNumer(0x0), | |
77 | fSmearedDenom(0x0), | |
78 | fDPtOverPtRMS(0.004), | |
79 | fThetaA(2.72e-03), | |
80 | fThetaB(1.87e-04), | |
81 | fThetaAlpha(-2.4), | |
82 | fPhiA(1.94e-03), | |
83 | fPhiB(2.5e-04), | |
84 | fPhiAlpha(-2.2), | |
85 | fR2(0.0), | |
86 | fLambda(0.0), | |
87 | fRConvergenceTreshold(0.3), | |
88 | fLambdaConvergenceTreshold(0.05) | |
89 | { | |
90 | ||
91 | } | |
92 | /******************************************************************/ | |
93 | ||
94 | AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name, const char* title, | |
95 | Int_t nbins, Float_t maxXval, Float_t minXval): | |
96 | AliHBTOnePairFctn1D(name,title,nbins,maxXval,minXval), | |
97 | fMeasCorrelFctn(0x0), | |
98 | fMeasNumer(0x0), | |
99 | fMeasDenom(0x0), | |
100 | fSmearedNumer(0x0), | |
101 | fSmearedDenom(0x0), | |
102 | fDPtOverPtRMS(0.004), | |
103 | fThetaA(2.72e-03), | |
104 | fThetaB(1.87e-04), | |
105 | fThetaAlpha(-2.4), | |
106 | fPhiA(1.94e-03), | |
107 | fPhiB(2.5e-04), | |
108 | fPhiAlpha(-2.2), | |
109 | fR2(0.0), | |
110 | fLambda(0.0), | |
111 | fRConvergenceTreshold(0.3), | |
112 | fLambdaConvergenceTreshold(0.05) | |
113 | { | |
e09fa876 | 114 | //ctor |
83b33650 | 115 | } |
116 | /******************************************************************/ | |
e09fa876 | 117 | AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const AliHBTCorrectQInvCorrelFctn& in): |
118 | AliHBTOnePairFctn1D(in), | |
119 | fMeasCorrelFctn(0x0), | |
120 | fMeasNumer(0x0), | |
121 | fMeasDenom(0x0), | |
122 | fSmearedNumer(0x0), | |
123 | fSmearedDenom(0x0), | |
124 | fDPtOverPtRMS(0), | |
125 | fThetaA(0), | |
126 | fThetaB(0), | |
127 | fThetaAlpha(0), | |
128 | fPhiA(0), | |
129 | fPhiB(0), | |
130 | fPhiAlpha(0), | |
131 | fR2(0.0), | |
132 | fLambda(0.0), | |
133 | fRConvergenceTreshold(0), | |
134 | fLambdaConvergenceTreshold(0) | |
135 | { | |
136 | //cpy ;ctor | |
137 | in.Copy(*this); | |
138 | } | |
83b33650 | 139 | |
140 | AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn() | |
141 | { | |
e09fa876 | 142 | //dtor |
83b33650 | 143 | delete fMeasCorrelFctn; |
144 | delete fSmearedNumer; | |
145 | delete fSmearedDenom; | |
146 | delete fMeasNumer; | |
147 | delete fMeasDenom; | |
148 | } | |
149 | /******************************************************************/ | |
150 | ||
151 | void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min) | |
152 | { | |
153 | AliHBTFunction1D::BuildHistos(nbins,max,min); | |
154 | ||
155 | TString numstr = fName + " Smeared Numerator"; //title and name of the numerator histogram | |
156 | TString denstr = fName + " Smeared Denominator";//title and name of the denominator histogram | |
157 | ||
158 | fSmearedNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max); | |
159 | fSmearedDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max); | |
160 | fSmearedNumer->Sumw2(); | |
161 | fSmearedDenom->Sumw2(); | |
162 | ||
163 | if (fMeasCorrelFctn == 0x0) | |
164 | { | |
165 | numstr = fName + " Measured Numerator"; //title and name of the numerator histogram | |
166 | denstr = fName + " Measured Denominator";//title and name of the denominator histogram | |
167 | ||
168 | fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max); | |
169 | fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max); | |
170 | fMeasNumer->Sumw2(); | |
171 | fMeasDenom->Sumw2(); | |
172 | } | |
173 | } | |
174 | /******************************************************************/ | |
175 | ||
176 | void AliHBTCorrectQInvCorrelFctn::Init() | |
177 | { | |
e09fa876 | 178 | //Init |
83b33650 | 179 | AliHBTOnePairFctn1D::Init(); |
180 | Info("Init",""); | |
181 | fSmearedNumer->Reset(); | |
182 | fSmearedDenom->Reset(); | |
183 | if (fMeasNumer) fMeasNumer->Reset(); | |
184 | if (fMeasDenom) fMeasDenom->Reset(); | |
185 | fFittedR = -1.0; | |
186 | fFittedLambda = 0.0; | |
187 | } | |
188 | /******************************************************************/ | |
189 | ||
190 | void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair) | |
191 | { | |
e09fa876 | 192 | //Processes particles that originates from the same event |
83b33650 | 193 | if (fMeasNumer == 0x0) return; |
194 | pair = CheckPair(pair); | |
195 | if( pair == 0x0) return; | |
196 | fMeasNumer->Fill(pair->GetQInv()); | |
197 | } | |
198 | /******************************************************************/ | |
199 | ||
200 | void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair) | |
201 | { | |
202 | //Process different events | |
78d7c6d3 | 203 | static AliAODParticle part1, part2; |
83b33650 | 204 | static AliHBTPair smearedpair(&part1,&part2); |
205 | ||
206 | pair = CheckPair(pair); | |
207 | if( pair == 0x0) return; | |
208 | Double_t cc = GetCoulombCorrection(pair); | |
209 | Double_t qinv = pair->GetQInv(); | |
210 | ||
211 | //measured histogram -> if we are interested | |
212 | //only if fMeasCorrelFctn is not specified by user | |
213 | if (fMeasDenom) fMeasDenom->Fill(qinv,cc); | |
214 | ||
215 | Smear(pair,smearedpair); | |
216 | Double_t modelqinv = GetModelValue(qinv); | |
217 | //Ideal histogram | |
218 | fNumerator->Fill(qinv,modelqinv*cc); | |
219 | fDenominator->Fill(qinv,cc); | |
220 | ||
221 | //Smeared histogram | |
222 | Double_t smearedqinv = smearedpair.GetQInv(); | |
223 | fSmearedNumer->Fill(smearedqinv,modelqinv); | |
224 | Double_t smearedcc = GetCoulombCorrection(&smearedpair); | |
225 | fSmearedDenom->Fill(smearedqinv,smearedcc); | |
226 | ||
227 | } | |
228 | /******************************************************************/ | |
229 | void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTPair* pair,AliHBTPair& smeared) | |
230 | { | |
231 | //Smears pair | |
232 | Smear(pair->Particle1(),smeared.Particle1()); | |
233 | Smear(pair->Particle2(),smeared.Particle2()); | |
234 | smeared.Changed(); | |
235 | } | |
236 | /******************************************************************/ | |
237 | ||
78d7c6d3 | 238 | void AliHBTCorrectQInvCorrelFctn::Smear(AliVAODParticle* part, AliVAODParticle* smeared) |
83b33650 | 239 | { |
e09fa876 | 240 | //Smears momenta |
83b33650 | 241 | Double_t sin2theta = TMath::Sin(part->Theta()); |
242 | sin2theta = sin2theta*sin2theta; | |
243 | Double_t pt = part->Pt(); | |
244 | ||
e09fa876 | 245 | double dPtDivPt = gRandom->Gaus(0.0,fDPtOverPtRMS); |
246 | double dphi = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fPhiAlpha)); | |
247 | double dtheta = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fThetaAlpha)); | |
83b33650 | 248 | |
e09fa876 | 249 | Double_t smearedPx = part->Px()*(1.0+dPtDivPt) - part->Py()*dphi; |
250 | // fourmom.setX(px*(1.0+dPtDivPt) - py*dphi); | |
251 | Double_t smearedPy = part->Py()*(1.0+dPtDivPt) - part->Px()*dphi; | |
252 | // fourmom.setY(py*(1.0+dPtDivPt) + px*dphi); | |
253 | Double_t smearedPz = part->Pz()*(1.0+dPtDivPt) - pt*dtheta/sin2theta; | |
254 | // fourmom.setZ(pz*(1.0+dPtDivPt) - pT*dtheta/sin2theta); | |
83b33650 | 255 | |
78d7c6d3 | 256 | Double_t mass2 = part->Mass()*part->Mass(); |
e09fa876 | 257 | Double_t e = mass2 + smearedPx*smearedPx + |
83b33650 | 258 | smearedPy*smearedPy + |
259 | smearedPz*smearedPz; | |
260 | ||
e09fa876 | 261 | smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(e)); |
83b33650 | 262 | } |
263 | /******************************************************************/ | |
264 | ||
265 | void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r) | |
266 | { | |
e09fa876 | 267 | //Sets Initial Values |
83b33650 | 268 | fLambda = lambda; |
269 | fR2 = r*r; | |
270 | } | |
271 | /******************************************************************/ | |
272 | void AliHBTCorrectQInvCorrelFctn::MakeMeasCF() | |
273 | { | |
274 | //makes measured correlation function | |
275 | delete fMeasCorrelFctn; | |
276 | fMeasCorrelFctn = 0x0; | |
277 | ||
278 | if (fMeasNumer&&fMeasDenom) | |
279 | { | |
280 | Double_t measscale = Scale(fMeasNumer,fMeasDenom); | |
281 | if (measscale == 0.0) | |
282 | { | |
283 | Error("GetResult","GetRatio for measured CF returned 0.0"); | |
284 | return; | |
285 | } | |
286 | TString str = fName + "measured ratio"; | |
287 | fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data()); | |
288 | fMeasCorrelFctn->SetTitle(str.Data()); | |
289 | fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale); | |
290 | } | |
291 | } | |
292 | ||
293 | TH1* AliHBTCorrectQInvCorrelFctn::GetResult() | |
294 | { | |
295 | //In case we don't have yet Measured Correlation Function | |
e09fa876 | 296 | //Try to get it |
83b33650 | 297 | //result is |
298 | // N[meas] N[ideal]/D[ideal] | |
299 | //C(Q) = ------- * ----------------- | |
300 | // D[meas] N[smear]/D[smear] | |
301 | ||
302 | TString str; | |
303 | if (fMeasCorrelFctn == 0x0) MakeMeasCF(); | |
304 | ||
305 | if (fMeasCorrelFctn == 0x0) | |
306 | { | |
307 | Error("GetResult", | |
308 | "Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null"); | |
309 | return 0x0; | |
310 | } | |
311 | ||
312 | TH1D* ideal = (TH1D*)GetRatio(Scale()); | |
313 | if (ideal == 0x0) | |
314 | { | |
315 | Error("GetResult","Ratio of ideal histograms is null"); | |
316 | return 0x0; | |
317 | } | |
318 | str = fName + " smeared ratio"; | |
319 | TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data()); | |
320 | smearedCF->SetTitle(str.Data()); | |
321 | Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom); | |
322 | smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale); | |
323 | ||
324 | str = fName + " product meas ideal CF"; | |
325 | TH1D* measideal = (TH1D*)ideal->Clone(str.Data()); | |
326 | measideal->Multiply(ideal,fMeasCorrelFctn); | |
327 | ||
328 | str = fName + " Corrected Result"; | |
329 | TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data()); | |
330 | result->SetTitle(str.Data()); | |
331 | ||
332 | Double_t resultscale = Scale(measideal,smearedCF); | |
333 | result->Divide(measideal,smearedCF,resultscale); | |
334 | return result; | |
335 | ||
336 | } | |
337 | /******************************************************************/ | |
338 | ||
339 | void AliHBTCorrectQInvCorrelFctn::Fit() | |
340 | { | |
341 | //fits resuting histogram with function 1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329)) | |
342 | //where [0] is lambda | |
343 | // [1] is radius | |
344 | // 0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI | |
345 | ||
346 | Info("Fit","Before fFittedLambda = %f",fFittedLambda); | |
347 | Info("Fit","Before fFittedR = %f",fFittedR); | |
348 | TH1D* result = (TH1D*)GetResult(); | |
349 | if (result == 0x0) | |
350 | { | |
351 | Error("Fit","Can not get result"); | |
352 | return; | |
353 | } | |
354 | TF1* fitfctn = new TF1("fitfctn","1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))"); | |
355 | fitfctn->SetParameter(0,1.0); | |
356 | fitfctn->SetParameter(1,6.0); | |
357 | Float_t max = result->GetXaxis()->GetXmax(); | |
358 | Info("Fit","Max is %f",max); | |
359 | result->Fit(fitfctn,"0","",0.008,max); | |
360 | fFittedLambda = fitfctn->GetParameter(0); | |
361 | fFittedR = fitfctn->GetParameter(1); | |
362 | Info("Fit","After fFittedLambda = %f",fFittedLambda); | |
363 | Info("Fit","After fFittedR = %f",fFittedR); | |
364 | delete fitfctn; | |
365 | delete result; | |
366 | } | |
367 | /******************************************************************/ | |
368 | ||
369 | Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged() | |
370 | { | |
371 | //check if fitting was performed | |
372 | if (fFittedR <= 0.0) | |
373 | { | |
374 | Fit();//if not do fit first | |
375 | if (fFittedR <= 0.0) | |
376 | { | |
377 | Error("IsConverged","Fitting failed"); | |
378 | return kFALSE; | |
379 | } | |
380 | } | |
381 | ||
382 | Double_t guessedR = TMath::Sqrt(fR2); | |
383 | ||
384 | Info("IsConverged","Fitted lambda : %8f Fitted Radius : %8f",fFittedLambda,fFittedR); | |
385 | Info("IsConverged","Guessed lambda : %8f Guessed Radius : %8f",fLambda,guessedR); | |
386 | Info("IsConverged","Demanded lambda convergence: %8f Demanded Radius convergence: %8f", | |
387 | fLambdaConvergenceTreshold,fRConvergenceTreshold); | |
388 | ||
389 | if ( (TMath::Abs(fLambda-fFittedLambda)<fLambdaConvergenceTreshold) && | |
390 | (TMath::Abs(fFittedR-guessedR)<fRConvergenceTreshold) ) | |
391 | { | |
392 | Info("IsConverged","Cnvergence reached"); | |
393 | return kTRUE; | |
394 | } | |
395 | else | |
396 | { | |
397 | Info("IsConverged","Cnvergence NOT reached"); | |
398 | return kFALSE; | |
399 | } | |
400 | } | |
401 | /******************************************************************/ | |
402 | ||
403 | Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius() | |
404 | { | |
e09fa876 | 405 | //Returns Fitted radius |
83b33650 | 406 | if (fFittedR <= 0.0) Fit(); |
407 | return fFittedR; | |
408 | } | |
409 | /******************************************************************/ | |
410 | ||
411 | Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda() | |
412 | { | |
e09fa876 | 413 | //Returns Fitted Intercept paramter |
83b33650 | 414 | if (fFittedR <= 0.0) Fit(); |
415 | return fFittedLambda; | |
416 | } | |
417 | /******************************************************************/ | |
418 | ||
419 | void AliHBTCorrectQInvCorrelFctn::WriteAll() | |
420 | { | |
e09fa876 | 421 | //Writes function and all additional information |
83b33650 | 422 | Write(); |
423 | if (fMeasCorrelFctn) fMeasCorrelFctn->Write(); | |
424 | if (fMeasNumer ) fMeasNumer->Write(); | |
425 | if (fMeasDenom ) fMeasDenom->Write(); | |
426 | if (fSmearedNumer) fSmearedNumer->Write(); | |
427 | if (fSmearedDenom) fSmearedDenom->Write(); | |
428 | if (fSmearedNumer && fSmearedDenom) | |
429 | { | |
430 | TString str = fName + " smeared ratio"; | |
431 | TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data()); | |
432 | smearedCF->SetTitle(str.Data()); | |
433 | Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom); | |
434 | smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale); | |
435 | } | |
436 | } | |
437 |