1 #include "AliHBTCorrectQInvCorrelFctn.h"
3 ///////////////////////////////////////////////////////
5 // AliHBTCorrectQInvCorrelFctn //
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 //
13 // Piotr.Skowronski@cern.ch //
14 // http://aliweb.cern.ch/people/skowron/analyzer //
16 ///////////////////////////////////////////////////////
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
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
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
60 #include <AliAODParticle.h>
63 ClassImp(AliHBTCorrectedCorrelFctn)
65 AliHBTCorrectedCorrelFctn::AliHBTCorrectedCorrelFctn():
66 fDPtOverPtA(5.78220e-03),
67 fDPtOverPtB(3.98063e-05),
68 fDPtOverPtAlpha(-2.78008),
69 fDPtOverPtC(5.07594e-04),
72 fThetaAlpha(-3.10218e+00),
73 fThetaC(-1.79892e-05),
76 fPhiAlpha(-5.26559e-02),
82 /******************************************************************/
83 void AliHBTCorrectedCorrelFctn::Smear(AliHBTPair* pair,AliHBTPair& smeared)
86 Smear(pair->Particle1(),smeared.Particle1());
87 Smear(pair->Particle2(),smeared.Particle2());
90 /******************************************************************/
92 void AliHBTCorrectedCorrelFctn::Smear(AliVAODParticle* part, AliVAODParticle* smeared)
95 Double_t sin2theta = TMath::Sin(part->Theta());
96 sin2theta = sin2theta*sin2theta;
97 Double_t pt = part->Pt();
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);
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);
111 Double_t mass2 = part->Mass()*part->Mass();
112 Double_t e = mass2 + smearedPx*smearedPx +
113 smearedPy*smearedPy +
116 smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(e));
119 /****************************************************************/
120 /****************************************************************/
121 /****************************************************************/
124 ClassImp(AliHBTCorrectQInvCorrelFctn)
126 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name,const char* title):
127 AliHBTOnePairFctn1D(name,title),
128 fMeasCorrelFctn(0x0),
135 fRConvergenceTreshold(0.3),
136 fLambdaConvergenceTreshold(0.05)
140 /******************************************************************/
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),
154 fRConvergenceTreshold(0.3),
155 fLambdaConvergenceTreshold(0.05)
159 /******************************************************************/
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),
164 AliHBTCorrectedCorrelFctn(),
165 fMeasCorrelFctn(0x0),
172 fRConvergenceTreshold(0.3),
173 fLambdaConvergenceTreshold(0.05)
177 /******************************************************************/
178 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const AliHBTCorrectQInvCorrelFctn& in):
179 AliHBTOnePairFctn1D(in),
180 AliHBTCorrectedCorrelFctn(),
181 fMeasCorrelFctn(0x0),
188 fRConvergenceTreshold(0),
189 fLambdaConvergenceTreshold(0)
195 AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
198 delete fMeasCorrelFctn;
199 delete fSmearedNumer;
200 delete fSmearedDenom;
204 /******************************************************************/
206 void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min)
208 AliHBTFunction1D::BuildHistos(nbins,max,min);
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
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();
218 if (fMeasCorrelFctn == 0x0)
220 numstr = fName + " Measured Numerator"; //title and name of the numerator histogram
221 denstr = fName + " Measured Denominator";//title and name of the denominator histogram
223 fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
224 fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
229 /******************************************************************/
231 void AliHBTCorrectQInvCorrelFctn::Init()
234 AliHBTOnePairFctn1D::Init();
236 fSmearedNumer->Reset();
237 fSmearedDenom->Reset();
238 if (fMeasNumer) fMeasNumer->Reset();
239 if (fMeasDenom) fMeasDenom->Reset();
243 /******************************************************************/
245 void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair)
247 //Processes particles that originates from the same event
248 if (fMeasNumer == 0x0) return;
249 pair = CheckPair(pair);
250 if( pair == 0x0) return;
251 fMeasNumer->Fill(pair->GetQInv());
253 /******************************************************************/
255 void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair)
257 //Process different events
258 static AliAODParticle part1, part2;
259 static AliHBTPair smearedpair(&part1,&part2);
261 pair = CheckPair(pair);
262 if( pair == 0x0) return;
263 Double_t cc = GetCoulombCorrection(pair);
264 Double_t qinv = pair->GetQInv();
266 //measured histogram -> if we are interested
267 //only if fMeasCorrelFctn is not specified by user
268 if (fMeasDenom) fMeasDenom->Fill(qinv,cc);
270 Smear(pair,smearedpair);
271 Double_t modelqinv = GetModelValue(qinv);
273 fNumerator->Fill(qinv,modelqinv*cc);
274 fDenominator->Fill(qinv,cc);
277 Double_t smearedqinv = smearedpair.GetQInv();
278 fSmearedNumer->Fill(smearedqinv,modelqinv);
279 Double_t smearedcc = GetCoulombCorrection(&smearedpair);
280 fSmearedDenom->Fill(smearedqinv,smearedcc);
283 /******************************************************************/
285 void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
287 //Sets Initial Values
291 /******************************************************************/
292 void AliHBTCorrectQInvCorrelFctn::MakeMeasCF()
294 //makes measured correlation function
295 delete fMeasCorrelFctn;
296 fMeasCorrelFctn = 0x0;
298 if (fMeasNumer&&fMeasDenom)
300 Double_t measscale = Scale(fMeasNumer,fMeasDenom);
301 if (measscale == 0.0)
303 Error("GetResult","GetRatio for measured CF returned 0.0");
306 TString str = fName + "measured ratio";
307 fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data());
308 fMeasCorrelFctn->SetTitle(str.Data());
309 fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale);
313 TH1* AliHBTCorrectQInvCorrelFctn::GetResult()
315 //In case we don't have yet Measured Correlation Function
318 // N[meas] N[ideal]/D[ideal]
319 //C(Q) = ------- * -----------------
320 // D[meas] N[smear]/D[smear]
323 if (fMeasCorrelFctn == 0x0) MakeMeasCF();
325 if (fMeasCorrelFctn == 0x0)
328 "Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null");
332 TH1D* ideal = (TH1D*)GetRatio(Scale());
335 Error("GetResult","Ratio of ideal histograms is null");
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);
344 str = fName + " product meas ideal CF";
345 TH1D* measideal = (TH1D*)ideal->Clone(str.Data());
346 measideal->Multiply(ideal,fMeasCorrelFctn);
348 str = fName + " Corrected Result";
349 TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data());
350 result->SetTitle(str.Data());
352 Double_t resultscale = Scale(measideal,smearedCF);
353 result->Divide(measideal,smearedCF,resultscale);
357 /******************************************************************/
359 void AliHBTCorrectQInvCorrelFctn::Fit()
361 //fits resuting histogram with function 1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))
362 //where [0] is lambda
364 // 0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI
366 Info("Fit","Before fFittedLambda = %f",fFittedLambda);
367 Info("Fit","Before fFittedR = %f",fFittedR);
368 TH1D* result = (TH1D*)GetResult();
371 Error("Fit","Can not get result");
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);
387 /******************************************************************/
389 Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged()
391 //check if fitting was performed
394 Fit();//if not do fit first
397 Error("IsConverged","Fitting failed");
402 Double_t guessedR = TMath::Sqrt(fR2);
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);
409 if ( (TMath::Abs(fLambda-fFittedLambda)<fLambdaConvergenceTreshold) &&
410 (TMath::Abs(fFittedR-guessedR)<fRConvergenceTreshold) )
412 Info("IsConverged","Cnvergence reached");
417 Info("IsConverged","Cnvergence NOT reached");
421 /******************************************************************/
423 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius()
425 //Returns Fitted radius
426 if (fFittedR <= 0.0) Fit();
429 /******************************************************************/
431 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
433 //Returns Fitted Intercept paramter
434 if (fFittedR <= 0.0) Fit();
435 return fFittedLambda;
437 /******************************************************************/
439 void AliHBTCorrectQInvCorrelFctn::WriteAll()
441 //Writes function and all additional information
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)
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);