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://alisoft.cern.ch/people/skowron/analyzer //
16 ///////////////////////////////////////////////////////
18 //Parameters fit from pi+ pi+ resolution analysis for pair with qinv<50MeV
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
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
40 #include <AliAODParticle.h>
43 ClassImp(AliHBTCorrectQInvCorrelFctn)
45 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name,const char* title):
46 AliHBTOnePairFctn1D(name,title),
61 fRConvergenceTreshold(0.3),
62 fLambdaConvergenceTreshold(0.05)
66 /******************************************************************/
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),
87 fRConvergenceTreshold(0.3),
88 fLambdaConvergenceTreshold(0.05)
92 /******************************************************************/
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),
102 fDPtOverPtRMS(0.004),
111 fRConvergenceTreshold(0.3),
112 fLambdaConvergenceTreshold(0.05)
116 /******************************************************************/
117 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const AliHBTCorrectQInvCorrelFctn& in):
118 AliHBTOnePairFctn1D(in),
119 fMeasCorrelFctn(0x0),
133 fRConvergenceTreshold(0),
134 fLambdaConvergenceTreshold(0)
140 AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
143 delete fMeasCorrelFctn;
144 delete fSmearedNumer;
145 delete fSmearedDenom;
149 /******************************************************************/
151 void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min)
153 AliHBTFunction1D::BuildHistos(nbins,max,min);
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
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();
163 if (fMeasCorrelFctn == 0x0)
165 numstr = fName + " Measured Numerator"; //title and name of the numerator histogram
166 denstr = fName + " Measured Denominator";//title and name of the denominator histogram
168 fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
169 fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
174 /******************************************************************/
176 void AliHBTCorrectQInvCorrelFctn::Init()
179 AliHBTOnePairFctn1D::Init();
181 fSmearedNumer->Reset();
182 fSmearedDenom->Reset();
183 if (fMeasNumer) fMeasNumer->Reset();
184 if (fMeasDenom) fMeasDenom->Reset();
188 /******************************************************************/
190 void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair)
192 //Processes particles that originates from the same event
193 if (fMeasNumer == 0x0) return;
194 pair = CheckPair(pair);
195 if( pair == 0x0) return;
196 fMeasNumer->Fill(pair->GetQInv());
198 /******************************************************************/
200 void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair)
202 //Process different events
203 static AliAODParticle part1, part2;
204 static AliHBTPair smearedpair(&part1,&part2);
206 pair = CheckPair(pair);
207 if( pair == 0x0) return;
208 Double_t cc = GetCoulombCorrection(pair);
209 Double_t qinv = pair->GetQInv();
211 //measured histogram -> if we are interested
212 //only if fMeasCorrelFctn is not specified by user
213 if (fMeasDenom) fMeasDenom->Fill(qinv,cc);
215 Smear(pair,smearedpair);
216 Double_t modelqinv = GetModelValue(qinv);
218 fNumerator->Fill(qinv,modelqinv*cc);
219 fDenominator->Fill(qinv,cc);
222 Double_t smearedqinv = smearedpair.GetQInv();
223 fSmearedNumer->Fill(smearedqinv,modelqinv);
224 Double_t smearedcc = GetCoulombCorrection(&smearedpair);
225 fSmearedDenom->Fill(smearedqinv,smearedcc);
228 /******************************************************************/
229 void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTPair* pair,AliHBTPair& smeared)
232 Smear(pair->Particle1(),smeared.Particle1());
233 Smear(pair->Particle2(),smeared.Particle2());
236 /******************************************************************/
238 void AliHBTCorrectQInvCorrelFctn::Smear(AliVAODParticle* part, AliVAODParticle* smeared)
241 Double_t sin2theta = TMath::Sin(part->Theta());
242 sin2theta = sin2theta*sin2theta;
243 Double_t pt = part->Pt();
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));
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);
256 Double_t mass2 = part->Mass()*part->Mass();
257 Double_t e = mass2 + smearedPx*smearedPx +
258 smearedPy*smearedPy +
261 smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(e));
263 /******************************************************************/
265 void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
267 //Sets Initial Values
271 /******************************************************************/
272 void AliHBTCorrectQInvCorrelFctn::MakeMeasCF()
274 //makes measured correlation function
275 delete fMeasCorrelFctn;
276 fMeasCorrelFctn = 0x0;
278 if (fMeasNumer&&fMeasDenom)
280 Double_t measscale = Scale(fMeasNumer,fMeasDenom);
281 if (measscale == 0.0)
283 Error("GetResult","GetRatio for measured CF returned 0.0");
286 TString str = fName + "measured ratio";
287 fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data());
288 fMeasCorrelFctn->SetTitle(str.Data());
289 fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale);
293 TH1* AliHBTCorrectQInvCorrelFctn::GetResult()
295 //In case we don't have yet Measured Correlation Function
298 // N[meas] N[ideal]/D[ideal]
299 //C(Q) = ------- * -----------------
300 // D[meas] N[smear]/D[smear]
303 if (fMeasCorrelFctn == 0x0) MakeMeasCF();
305 if (fMeasCorrelFctn == 0x0)
308 "Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null");
312 TH1D* ideal = (TH1D*)GetRatio(Scale());
315 Error("GetResult","Ratio of ideal histograms is null");
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);
324 str = fName + " product meas ideal CF";
325 TH1D* measideal = (TH1D*)ideal->Clone(str.Data());
326 measideal->Multiply(ideal,fMeasCorrelFctn);
328 str = fName + " Corrected Result";
329 TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data());
330 result->SetTitle(str.Data());
332 Double_t resultscale = Scale(measideal,smearedCF);
333 result->Divide(measideal,smearedCF,resultscale);
337 /******************************************************************/
339 void AliHBTCorrectQInvCorrelFctn::Fit()
341 //fits resuting histogram with function 1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))
342 //where [0] is lambda
344 // 0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI
346 Info("Fit","Before fFittedLambda = %f",fFittedLambda);
347 Info("Fit","Before fFittedR = %f",fFittedR);
348 TH1D* result = (TH1D*)GetResult();
351 Error("Fit","Can not get result");
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);
367 /******************************************************************/
369 Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged()
371 //check if fitting was performed
374 Fit();//if not do fit first
377 Error("IsConverged","Fitting failed");
382 Double_t guessedR = TMath::Sqrt(fR2);
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);
389 if ( (TMath::Abs(fLambda-fFittedLambda)<fLambdaConvergenceTreshold) &&
390 (TMath::Abs(fFittedR-guessedR)<fRConvergenceTreshold) )
392 Info("IsConverged","Cnvergence reached");
397 Info("IsConverged","Cnvergence NOT reached");
401 /******************************************************************/
403 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius()
405 //Returns Fitted radius
406 if (fFittedR <= 0.0) Fit();
409 /******************************************************************/
411 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
413 //Returns Fitted Intercept paramter
414 if (fFittedR <= 0.0) Fit();
415 return fFittedLambda;
417 /******************************************************************/
419 void AliHBTCorrectQInvCorrelFctn::WriteAll()
421 //Writes function and all additional information
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)
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);