1 #include "AliHBTCorrectCorrelFctn.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 ClassImp(AliHBTCorrectQInvCorrelFctn)
42 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name,const char* title):
43 AliHBTOnePairFctn1D(name,title),
58 fRConvergenceTreshold(0.3),
59 fLambdaConvergenceTreshold(0.05)
63 /******************************************************************/
65 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(TH1D* measqinv,const char* name,const char* title):
66 AliHBTOnePairFctn1D(name,title,
67 measqinv->GetNbinsX(),
68 measqinv->GetXaxis()->GetXmax(),
69 measqinv->GetXaxis()->GetXmin()),
70 fMeasCorrelFctn(measqinv),
84 fRConvergenceTreshold(0.3),
85 fLambdaConvergenceTreshold(0.05)
89 /******************************************************************/
91 AliHBTCorrectQInvCorrelFctn::AliHBTCorrectQInvCorrelFctn(const char* name, const char* title,
92 Int_t nbins, Float_t maxXval, Float_t minXval):
93 AliHBTOnePairFctn1D(name,title,nbins,maxXval,minXval),
108 fRConvergenceTreshold(0.3),
109 fLambdaConvergenceTreshold(0.05)
112 /******************************************************************/
114 AliHBTCorrectQInvCorrelFctn::~AliHBTCorrectQInvCorrelFctn()
116 delete fMeasCorrelFctn;
117 delete fSmearedNumer;
118 delete fSmearedDenom;
122 /******************************************************************/
124 void AliHBTCorrectQInvCorrelFctn::BuildHistos(Int_t nbins, Float_t max, Float_t min)
126 AliHBTFunction1D::BuildHistos(nbins,max,min);
128 TString numstr = fName + " Smeared Numerator"; //title and name of the numerator histogram
129 TString denstr = fName + " Smeared Denominator";//title and name of the denominator histogram
131 fSmearedNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
132 fSmearedDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
133 fSmearedNumer->Sumw2();
134 fSmearedDenom->Sumw2();
136 if (fMeasCorrelFctn == 0x0)
138 numstr = fName + " Measured Numerator"; //title and name of the numerator histogram
139 denstr = fName + " Measured Denominator";//title and name of the denominator histogram
141 fMeasNumer = new TH1D(numstr.Data(),numstr.Data(),nbins,min,max);
142 fMeasDenom = new TH1D(denstr.Data(),denstr.Data(),nbins,min,max);
147 /******************************************************************/
149 void AliHBTCorrectQInvCorrelFctn::Init()
151 AliHBTOnePairFctn1D::Init();
153 fSmearedNumer->Reset();
154 fSmearedDenom->Reset();
155 if (fMeasNumer) fMeasNumer->Reset();
156 if (fMeasDenom) fMeasDenom->Reset();
160 /******************************************************************/
162 void AliHBTCorrectQInvCorrelFctn::ProcessSameEventParticles(AliHBTPair* pair)
164 if (fMeasNumer == 0x0) return;
165 pair = CheckPair(pair);
166 if( pair == 0x0) return;
167 fMeasNumer->Fill(pair->GetQInv());
169 /******************************************************************/
171 void AliHBTCorrectQInvCorrelFctn::ProcessDiffEventParticles(AliHBTPair* pair)
173 //Process different events
174 static AliHBTParticle part1, part2;
175 static AliHBTPair smearedpair(&part1,&part2);
177 pair = CheckPair(pair);
178 if( pair == 0x0) return;
179 Double_t cc = GetCoulombCorrection(pair);
180 Double_t qinv = pair->GetQInv();
182 //measured histogram -> if we are interested
183 //only if fMeasCorrelFctn is not specified by user
184 if (fMeasDenom) fMeasDenom->Fill(qinv,cc);
186 Smear(pair,smearedpair);
187 Double_t modelqinv = GetModelValue(qinv);
189 fNumerator->Fill(qinv,modelqinv*cc);
190 fDenominator->Fill(qinv,cc);
193 Double_t smearedqinv = smearedpair.GetQInv();
194 fSmearedNumer->Fill(smearedqinv,modelqinv);
195 Double_t smearedcc = GetCoulombCorrection(&smearedpair);
196 fSmearedDenom->Fill(smearedqinv,smearedcc);
199 /******************************************************************/
200 void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTPair* pair,AliHBTPair& smeared)
203 Smear(pair->Particle1(),smeared.Particle1());
204 Smear(pair->Particle2(),smeared.Particle2());
207 /******************************************************************/
209 void AliHBTCorrectQInvCorrelFctn::Smear(AliHBTParticle* part, AliHBTParticle* smeared)
211 Double_t sin2theta = TMath::Sin(part->Theta());
212 sin2theta = sin2theta*sin2theta;
213 Double_t pt = part->Pt();
215 double DpT_div_pT = gRandom->Gaus(0.0,fDPtOverPtRMS);
216 double Dphi = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fPhiAlpha));
217 double Dtheta = gRandom->Gaus(0.0,fPhiA+fPhiB*TMath::Power(pt,fThetaAlpha));
219 Double_t smearedPx = part->Px()*(1.0+DpT_div_pT) - part->Py()*Dphi;
220 // fourmom.setX(px*(1.0+DpT_div_pT) - py*Dphi);
221 Double_t smearedPy = part->Py()*(1.0+DpT_div_pT) - part->Px()*Dphi;
222 // fourmom.setY(py*(1.0+DpT_div_pT) + px*Dphi);
223 Double_t smearedPz = part->Pz()*(1.0+DpT_div_pT) - pt*Dtheta/sin2theta;
224 // fourmom.setZ(pz*(1.0+DpT_div_pT) - pT*Dtheta/sin2theta);
226 Double_t mass2 = part->GetMass()*part->GetMass();
227 Double_t E = mass2 + smearedPx*smearedPx +
228 smearedPy*smearedPy +
231 smeared->SetMomentum(smearedPx,smearedPy,smearedPz,TMath::Sqrt(E));
233 /******************************************************************/
235 void AliHBTCorrectQInvCorrelFctn::SetInitialValues(Double_t lambda, Double_t r)
240 /******************************************************************/
241 void AliHBTCorrectQInvCorrelFctn::MakeMeasCF()
243 //makes measured correlation function
244 delete fMeasCorrelFctn;
245 fMeasCorrelFctn = 0x0;
247 if (fMeasNumer&&fMeasDenom)
249 Double_t measscale = Scale(fMeasNumer,fMeasDenom);
250 if (measscale == 0.0)
252 Error("GetResult","GetRatio for measured CF returned 0.0");
255 TString str = fName + "measured ratio";
256 fMeasCorrelFctn = (TH1D*)fMeasNumer->Clone(str.Data());
257 fMeasCorrelFctn->SetTitle(str.Data());
258 fMeasCorrelFctn->Divide(fMeasNumer,fMeasDenom,measscale);
262 TH1* AliHBTCorrectQInvCorrelFctn::GetResult()
264 //In case we don't have yet Measured Correlation Function
267 // N[meas] N[ideal]/D[ideal]
268 //C(Q) = ------- * -----------------
269 // D[meas] N[smear]/D[smear]
272 if (fMeasCorrelFctn == 0x0) MakeMeasCF();
274 if (fMeasCorrelFctn == 0x0)
277 "Measured Correlation Function is not defined and measured numeraor and/or denominator are/is null");
281 TH1D* ideal = (TH1D*)GetRatio(Scale());
284 Error("GetResult","Ratio of ideal histograms is null");
287 str = fName + " smeared ratio";
288 TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
289 smearedCF->SetTitle(str.Data());
290 Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
291 smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
293 str = fName + " product meas ideal CF";
294 TH1D* measideal = (TH1D*)ideal->Clone(str.Data());
295 measideal->Multiply(ideal,fMeasCorrelFctn);
297 str = fName + " Corrected Result";
298 TH1D* result = (TH1D*)fSmearedNumer->Clone(str.Data());
299 result->SetTitle(str.Data());
301 Double_t resultscale = Scale(measideal,smearedCF);
302 result->Divide(measideal,smearedCF,resultscale);
306 /******************************************************************/
308 void AliHBTCorrectQInvCorrelFctn::Fit()
310 //fits resuting histogram with function 1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))
311 //where [0] is lambda
313 // 0.038936366329 - constant needed for units transformation eV(c=1,etc.) -> SI
315 Info("Fit","Before fFittedLambda = %f",fFittedLambda);
316 Info("Fit","Before fFittedR = %f",fFittedR);
317 TH1D* result = (TH1D*)GetResult();
320 Error("Fit","Can not get result");
323 TF1* fitfctn = new TF1("fitfctn","1.0 + [0]*exp([1]*[1]*x*x/(-0.038936366329))");
324 fitfctn->SetParameter(0,1.0);
325 fitfctn->SetParameter(1,6.0);
326 Float_t max = result->GetXaxis()->GetXmax();
327 Info("Fit","Max is %f",max);
328 result->Fit(fitfctn,"0","",0.008,max);
329 fFittedLambda = fitfctn->GetParameter(0);
330 fFittedR = fitfctn->GetParameter(1);
331 Info("Fit","After fFittedLambda = %f",fFittedLambda);
332 Info("Fit","After fFittedR = %f",fFittedR);
336 /******************************************************************/
338 Bool_t AliHBTCorrectQInvCorrelFctn::IsConverged()
340 //check if fitting was performed
343 Fit();//if not do fit first
346 Error("IsConverged","Fitting failed");
351 Double_t guessedR = TMath::Sqrt(fR2);
353 Info("IsConverged","Fitted lambda : %8f Fitted Radius : %8f",fFittedLambda,fFittedR);
354 Info("IsConverged","Guessed lambda : %8f Guessed Radius : %8f",fLambda,guessedR);
355 Info("IsConverged","Demanded lambda convergence: %8f Demanded Radius convergence: %8f",
356 fLambdaConvergenceTreshold,fRConvergenceTreshold);
358 if ( (TMath::Abs(fLambda-fFittedLambda)<fLambdaConvergenceTreshold) &&
359 (TMath::Abs(fFittedR-guessedR)<fRConvergenceTreshold) )
361 Info("IsConverged","Cnvergence reached");
366 Info("IsConverged","Cnvergence NOT reached");
370 /******************************************************************/
372 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedRadius()
374 if (fFittedR <= 0.0) Fit();
377 /******************************************************************/
379 Double_t AliHBTCorrectQInvCorrelFctn::GetFittedLambda()
381 if (fFittedR <= 0.0) Fit();
382 return fFittedLambda;
384 /******************************************************************/
386 void AliHBTCorrectQInvCorrelFctn::WriteAll()
389 if (fMeasCorrelFctn) fMeasCorrelFctn->Write();
390 if (fMeasNumer ) fMeasNumer->Write();
391 if (fMeasDenom ) fMeasDenom->Write();
392 if (fSmearedNumer) fSmearedNumer->Write();
393 if (fSmearedDenom) fSmearedDenom->Write();
394 if (fSmearedNumer && fSmearedDenom)
396 TString str = fName + " smeared ratio";
397 TH1D* smearedCF = (TH1D*)fSmearedNumer->Clone(str.Data());
398 smearedCF->SetTitle(str.Data());
399 Double_t smearedscale = Scale(fSmearedNumer,fSmearedDenom);
400 smearedCF->Divide(fSmearedNumer,fSmearedDenom,smearedscale);
404 /******************************************************************/
405 /******************************************************************/
406 /******************************************************************/
408 //____________________
409 ///////////////////////////////////////////////////////
411 // AliHBTCorrectQ3DCorrelFctn //
413 // Class for calculating Q Invariant correlation //
414 // taking to the account resolution of the //
415 // detector and coulomb effects. //
417 ///////////////////////////////////////////////////////
420 AliHBTCorrectQ3DCorrelFctn::AliHBTCorrectQ3DCorrelFctn(const char* name, const char* title):
421 AliHBTOnePairFctn3D(name,title),
422 fMeasCorrelFctn(0x0),
430 /******************************************************************/
432 AliHBTCorrectQ3DCorrelFctn::~AliHBTCorrectQ3DCorrelFctn()
434 delete fMeasCorrelFctn;
435 delete fSmearedNumer;
436 delete fSmearedDenom;
441 /******************************************************************/