]>
Commit | Line | Data |
---|---|---|
1 | /************************************************************************* | |
2 | * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | // ------------------------------------------------------- | |
17 | // pt spectra extrapolation in the 0. - 0.2 region using | |
18 | // Boltzmann-Gibbs Blast Wave model or Tsallis Blast | |
19 | // Wave model for azimuthal isotropic expansion in | |
20 | // highly central collisions analysis | |
21 | // author: Cristian Andrei | |
22 | // acristian@niham.nipne.ro | |
23 | // ---------------------------------------------------------- | |
24 | ||
25 | ||
26 | #include "TH1D.h" | |
27 | #include "TFile.h" | |
28 | #include "TList.h" | |
29 | #include "TCanvas.h" | |
30 | #include "TLegend.h" | |
31 | // #include "TString.h" | |
32 | // #include "TPaveStats.h" | |
33 | #include "TMath.h" | |
34 | #include "TStopwatch.h" | |
35 | #include "TStyle.h" | |
36 | #include "TF1.h" | |
37 | ||
38 | #include "TVirtualFitter.h" | |
39 | #include "Math/WrappedFunction.h" | |
40 | #include "Math/Integrator.h" | |
41 | #include "Math/IntegratorMultiDim.h" | |
42 | ||
43 | #include "AliCFContainer.h" | |
44 | #include "AliCFDataGrid.h" | |
45 | #include "AliCFEffGrid.h" | |
46 | ||
47 | #include "AliAnalysisCentralExtrapolate.h" | |
48 | ||
49 | ||
50 | ClassImp(AliAnalysisCentralExtrapolate) | |
51 | ||
52 | //________________________________________________________________________ | |
53 | AliAnalysisCentralExtrapolate::AliAnalysisCentralExtrapolate(const char *name) | |
54 | : TObject() | |
55 | ,fPartType() | |
56 | ,fInputList(0) | |
57 | ,fResultsList(0x0) | |
58 | ||
59 | { | |
60 | // Constructor | |
61 | ||
62 | Printf("Running %s!\n", name); | |
63 | ||
64 | fInputList = new TList(); | |
65 | fResultsList = new TList(); | |
66 | ||
67 | } | |
68 | ||
69 | AliAnalysisCentralExtrapolate::~AliAnalysisCentralExtrapolate() { | |
70 | ||
71 | if(fInputList) delete fInputList; | |
72 | if(fResultsList) delete fResultsList; | |
73 | ||
74 | } | |
75 | ||
76 | //************************************************************************************ | |
77 | //____________________________________________________________________________________ | |
78 | // using global variables to avoid the limitations of ROOT::Math::WrappedMultiFunction<> | |
79 | static Double_t mass = 0.; | |
80 | static Double_t pt = 0.; | |
81 | static Double_t T = 0.; | |
82 | static Double_t betaS = 0.; | |
83 | static Double_t q = 0.; | |
84 | static Double_t mt = 0.; | |
85 | ||
86 | void AliAnalysisCentralExtrapolate::ApplyEff(){ | |
87 | // applies the efficiency map to the pt spectra | |
88 | ||
89 | TH1::SetDefaultSumw2(); | |
90 | // Int_t stepGen = 0; | |
91 | Int_t stepRec = 1; | |
92 | ||
93 | AliCFContainer *data = 0; | |
94 | ||
95 | ||
96 | TFile *file1 = new TFile("$ALICE_ROOT/PWG2/data/AliAnalysisCentralEfficiency.root", "read"); | |
97 | ||
98 | AliCFEffGrid *eff = 0; | |
99 | ||
100 | if(fPartType.Contains("kPi")){ | |
101 | data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_Pi")); | |
102 | eff = dynamic_cast<AliCFEffGrid*>(file1->Get("eff_Pi")); | |
103 | mass= 0.13957;//(GeV - pions) | |
104 | // ccorrdata =new TCanvas("Pions ccorrdata","Pions - corrected data",0,0,600,800); | |
105 | printf("\n\n**************************************\n"); | |
106 | printf("\tRunning for pions!\n"); | |
107 | } | |
108 | else if(fPartType.Contains("kK")){ | |
109 | data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_K")); | |
110 | eff = dynamic_cast<AliCFEffGrid*>(file1->Get("eff_K")); | |
111 | mass= 0.49368;//(GeV - kaons) | |
112 | // ccorrdata =new TCanvas("Kaons ccorrdata","Kaons - corrected data",0,0,600,800); | |
113 | printf("\n\n**************************************\n"); | |
114 | printf("\tRunning for kaons!\n"); | |
115 | } | |
116 | else if(fPartType.Contains("kProton")){ | |
117 | data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_P")); | |
118 | eff = dynamic_cast<AliCFEffGrid*>(file1->Get("eff_P")); | |
119 | mass = 0.93827;//(GeV - protons) | |
120 | // ccorrdata =new TCanvas("Protons ccorrdata","Protons - corrected data",0,0,600,800); | |
121 | printf("\n\n**************************************\n"); | |
122 | printf("\tRunning for protons!\n"); | |
123 | } | |
124 | else printf("Unsupported particle type!\n"); | |
125 | ||
126 | if(!data){ | |
127 | printf("Unable to get CFContainer! \n"); | |
128 | return; | |
129 | } | |
130 | ||
131 | if(!eff){ | |
132 | printf("No Eff Grid found! \n"); | |
133 | return; | |
134 | } | |
135 | ||
136 | // TCanvas *ccorrdata = new TCanvas(); | |
137 | // ccorrdata->Divide(1,2); | |
138 | // ccorrdata->cd(1); | |
139 | // ccorrdata->cd(1)->SetLogy(); | |
140 | ||
141 | ||
142 | AliCFDataGrid *corrdata = new AliCFDataGrid("corrdata","corrected data",*data, stepRec); | |
143 | ||
144 | //correct selection step "reconstructed" | |
145 | // corrdata->SetMeasured(stepRec); //set data to be corrected | |
146 | corrdata->ApplyEffCorrection(*eff);//apply the correction for efficiency | |
147 | ||
148 | // TH1D *hPtMC = data->ShowProjection(0,0); //MC distribution ShowProjection(ivar, istep) | |
149 | // hPtMC->SetMarkerStyle(20); | |
150 | // hPtMC->SetMarkerColor(kGreen-3); | |
151 | // hPtMC->GetXaxis()->SetTitle("p_{T}(GeV/c)"); | |
152 | // hPtMC->GetYaxis()->SetTitle("#frac{dN}{p_{T}dp_{T}}"); | |
153 | // hPtMC->Draw("p e1"); | |
154 | // | |
155 | // TH1D *hPtESDI = corrdata->GetData()->Project(0); //uncorrected ESD Project(ivar) | |
156 | // hPtESDI->SetMarkerStyle(25); | |
157 | // hPtESDI->SetMarkerColor(kBlue); | |
158 | // hPtESDI->GetXaxis()->SetTitle("Pt"); | |
159 | // hPtESDI->Draw("p e1 same"); | |
160 | ||
161 | TH1D *hPtESDCorr = (TH1D*)corrdata->Project(0); //corrected data | |
162 | hPtESDCorr->SetMarkerStyle(26); | |
163 | hPtESDCorr->SetMarkerColor(kRed); //ESD corrected | |
164 | hPtESDCorr->SetName("hPtESDCorr"); | |
165 | // hPtESDCorr->Draw("p e1 same"); | |
166 | ||
167 | // TLegend *leg2 = new TLegend(0.40,0.65,0.87,0.8); | |
168 | // leg2->AddEntry(hPtMC," generated","p"); | |
169 | // leg2->AddEntry(hPtESDI," reconstructed (not corrected)","p"); | |
170 | // leg2->AddEntry(hPtESDCorr," reconstructed & corrected","p"); | |
171 | // leg2->Draw(); | |
172 | ||
173 | ||
174 | // ccorrdata->cd(2); | |
175 | // ccorrdata->cd(2)->SetLogy(); | |
176 | // | |
177 | // TH1D *efficiency = eff->Project(0); //the efficiency vs pt (ipt = 0) | |
178 | // efficiency->GetXaxis()->SetTitle("p_{T}(GeV/c)"); | |
179 | // efficiency->Draw(); | |
180 | ||
181 | ||
182 | TH1D *hNoEvt = dynamic_cast<TH1D*>(fInputList->FindObject("TaskCentral_NoEvt")); | |
183 | if(!hNoEvt){ | |
184 | printf("Unable to get the number of events! \n"); | |
185 | return; | |
186 | } | |
187 | ||
188 | ||
189 | Int_t noEvt = (Int_t)(hNoEvt->GetEntries()); | |
190 | ||
191 | printf("\n** No of processed events = %i **\n",noEvt); | |
192 | ||
193 | Double_t scale = 1.0/noEvt; | |
194 | ||
195 | TH1D *hPtESDCorrNorm = (TH1D*)hPtESDCorr->Clone("hPtESDCorrNorm"); | |
196 | hPtESDCorrNorm->Scale(scale); | |
197 | ||
198 | fResultsList->SetName(fPartType); | |
199 | fResultsList->Add(hPtESDCorr); | |
200 | fResultsList->Add(hPtESDCorrNorm); | |
201 | ||
202 | } | |
203 | ||
204 | ||
205 | //********************************************************************************* | |
206 | //_________________________________________________________________________________ | |
207 | // fit Tsallis | |
208 | ||
209 | ||
210 | Double_t Integ1(const Double_t *x){ | |
211 | //the function under the integral - TBW | |
212 | ||
213 | const Double_t rMax = 11.0;//(fm) | |
214 | ||
215 | Double_t betaR = betaS*(x[0]/rMax); | |
216 | ||
217 | ||
218 | Double_t rho = TMath::ATanH(betaR); | |
219 | ||
220 | ||
221 | Double_t temp = 1.0+((q-1.0)*(mt*TMath::CosH(x[2])*TMath::CosH(rho) - pt*TMath::SinH(rho)*TMath::Cos(x[1]))/T); | |
222 | ||
223 | Double_t power = -1.0/(q-1.0); | |
224 | ||
225 | Double_t temp1 = pow(temp, power); | |
226 | ||
227 | Double_t f = TMath::CosH(x[1])*x[0]*temp1; | |
228 | ||
229 | return f; | |
230 | } | |
231 | ||
232 | ||
233 | // TF1 requires the function to have the ( )( Double_t *, Double_t *) signature | |
234 | Double_t Tsallis(Double_t* const x, Double_t* const par){ | |
235 | //computes the triple integral in the TBW formula | |
236 | ||
237 | pt = x[0]; | |
238 | ||
239 | T = par[0]; | |
240 | betaS = par[1]; | |
241 | q = par[2]; | |
242 | ||
243 | mt = sqrt(mass*mass+ pt*pt); | |
244 | ||
245 | ||
246 | ROOT::Math::WrappedMultiFunction<> f1(Integ1,3); | |
247 | ||
248 | ROOT::Math::IntegratorMultiDim ig(f1, ROOT::Math::IntegrationMultiDim::kPLAIN); | |
249 | ||
250 | ||
251 | Double_t xmin[3] = {0.0, -TMath::Pi(), -0.5}; //rmin, phi_min, ymin | |
252 | Double_t xmax[3] = {11.0, TMath::Pi(), 0.5}; //rmax, phi_max, ymax | |
253 | ||
254 | Double_t rez = mt*ig.Integral(xmin,xmax); | |
255 | ||
256 | return rez; | |
257 | } | |
258 | ||
259 | ||
260 | void AliAnalysisCentralExtrapolate::TsallisFit(){ | |
261 | // fits and extrpolates the pt spectrum using TBW model | |
262 | ||
263 | printf("---------------------------------------------\n"); | |
264 | printf("***** Tsallis Blast Wave Fit *****\n"); | |
265 | ||
266 | printf("AliAnalysisCentralExtrapolate::Tsallis mass = %g\n", mass); | |
267 | ||
268 | ||
269 | TStopwatch timer; | |
270 | ||
271 | TH1::SetDefaultSumw2(); | |
272 | TH1::AddDirectory(0); | |
273 | ||
274 | timer.Start(); | |
275 | ||
276 | TH1D *ptSpectra = dynamic_cast<TH1D*>(fResultsList->FindObject("hPtESDCorrNorm")); | |
277 | if(!ptSpectra){ | |
278 | printf("TsallisFit: Can't get the normalized spectrum\n"); | |
279 | return; | |
280 | } | |
281 | ||
282 | if(!ptSpectra->GetEntries()){ | |
283 | printf("TsallisFit: The fit data is empty!\n"); | |
284 | return; | |
285 | } | |
286 | ||
287 | TVirtualFitter::SetDefaultFitter("Minuit2");//options Minuit, Minuit2, Fumili, Fumili2 | |
288 | ||
289 | TF1 *ptFit = new TF1("ptFit",Tsallis,0.2,4.0,3); | |
290 | ||
291 | gStyle->SetOptFit(1112); | |
292 | ||
293 | // ptFit->SetParName(0,"T"); | |
294 | // ptFit->SetParName(1,"betaS"); | |
295 | // ptFit->SetParName(2,"q"); | |
296 | ||
297 | ptFit->SetParameters(0.1,0.5,1.05); //GeV | |
298 | ptFit->SetParLimits(0,0.0,0.5);//GeV | |
299 | ptFit->SetParLimits(1,0.0,1.0); | |
300 | ptFit->SetParLimits(2,1.00001,1.5); | |
301 | ||
302 | ptSpectra->Fit("ptFit","Q R B ME N"); | |
303 | ||
304 | timer.Stop(); | |
305 | timer.Print(); | |
306 | ||
307 | // ----------- print the fit result ---------- | |
308 | Double_t chi2 = ptFit->GetChisquare(); | |
309 | Double_t ndf = ptFit->GetNDF(); | |
310 | Double_t p1 = ptFit->GetParameter(0); | |
311 | Double_t param1Err = ptFit->GetParError(0); | |
312 | Double_t p2 = ptFit->GetParameter(1); | |
313 | Double_t param2Err = ptFit->GetParError(1); | |
314 | Double_t p3 = ptFit->GetParameter(2); | |
315 | Double_t param3Err = ptFit->GetParError(2); | |
316 | ||
317 | printf("chi2/ndf = %f/%f\n", chi2,ndf); | |
318 | printf("p1 = %f\tparam1Err = %f\n", p1, param1Err); | |
319 | printf("p2 = %f\tparam2Err = %f\n", p2, param2Err); | |
320 | printf("p3 = %f\tparam3Err = %f\n", p3, param3Err); | |
321 | ||
322 | ||
323 | // create a new, "extended" histogram | |
324 | TH1D *hPtExtTsallis = new TH1D("PtExtTsallis","Pt Corr Norm Ext",25,0.0,5.0); | |
325 | ||
326 | Double_t bin, binerr, test; | |
327 | ||
328 | for(Int_t i=0; i<ptSpectra->GetNbinsX()+1;i++){ | |
329 | ||
330 | bin = ptSpectra->GetBinContent(i); | |
331 | binerr = ptSpectra->GetBinError(i); | |
332 | test = ptSpectra->GetBinLowEdge(i); | |
333 | ||
334 | hPtExtTsallis->SetBinContent(i, bin); | |
335 | hPtExtTsallis->SetBinError(i, binerr); | |
336 | } | |
337 | ||
338 | Double_t eval; | |
339 | eval = ptFit->Eval(0.1); | |
340 | printf("extrapolated pt value = %g \n", eval); | |
341 | printf("****************************************************\n\n"); | |
342 | ||
343 | hPtExtTsallis->SetBinContent(1, eval); | |
344 | ||
345 | hPtExtTsallis->SetMarkerStyle(24); | |
346 | hPtExtTsallis->SetMarkerColor(kRed); //ESD corrected | |
347 | hPtExtTsallis->GetXaxis()->SetTitle("Pt"); | |
348 | ||
349 | ||
350 | fResultsList->Add(hPtExtTsallis); | |
351 | ||
352 | } | |
353 | ||
354 | //********************************************************************************* | |
355 | //_________________________________________________________________________________ | |
356 | // fit Boltzmann | |
357 | ||
358 | Double_t func(Double_t r){ | |
359 | //the function under the integral - BGBW | |
360 | ||
361 | const Double_t rMax = 11.0;//(fm) | |
362 | ||
363 | Double_t betaR = betaS*(r/rMax); | |
364 | ||
365 | ||
366 | Double_t rho = TMath::ATanH(betaR); | |
367 | ||
368 | ||
369 | mt = sqrt(mass*mass + pt*pt);// !!! par[0]= pt !!!!! | |
370 | ||
371 | Double_t argI0 = (pt*TMath::SinH(rho))/T; //T = par[1] | |
372 | if(argI0>700.0) return 0.0; // !! floating point exception protection | |
373 | ||
374 | Double_t argK1 = (mt*TMath::CosH(rho))/T; | |
375 | ||
376 | ||
377 | Double_t i0 = TMath::BesselI0(argI0); | |
378 | ||
379 | ||
380 | Double_t k1 = TMath::BesselK1(argK1); | |
381 | ||
382 | ||
383 | Double_t f = r*mt*i0*k1; | |
384 | ||
385 | return f; | |
386 | } | |
387 | ||
388 | ||
389 | Double_t Boltzmann(Double_t* const x, Double_t* const par){ | |
390 | //computes the integral in the BGBW formula | |
391 | ||
392 | pt = x[0]; | |
393 | T = par[0]; | |
394 | ||
395 | betaS = par[1]; | |
396 | ||
397 | ||
398 | ROOT::Math::WrappedFunction<> f1(func); | |
399 | ||
400 | ROOT::Math::Integrator ig(f1, ROOT::Math::IntegrationOneDim::kGAUSS,1.E-12,1.E-12,10000); | |
401 | ||
402 | Double_t rez = ig.Integral(0.0,11.0); | |
403 | ||
404 | return rez; | |
405 | } | |
406 | ||
407 | ||
408 | void AliAnalysisCentralExtrapolate::BoltzmannFit(){ | |
409 | //fits and extrapoates the pt spectrum using the BGBW model | |
410 | ||
411 | printf("---------------------------------------------------\n"); | |
412 | printf("*** Boltzmann-Gibbs Blast Wave Fit ***\n"); | |
413 | ||
414 | printf("AliAnalysisCentralExtrapolate::Boltzmann mass = %g\n", mass); | |
415 | ||
416 | ||
417 | TStopwatch timer; | |
418 | ||
419 | TH1::SetDefaultSumw2(); | |
420 | TH1::AddDirectory(0); | |
421 | ||
422 | timer.Start(); | |
423 | ||
424 | ||
425 | TH1D *ptSpectra = dynamic_cast<TH1D*>(fResultsList->FindObject("hPtESDCorrNorm")); | |
426 | if(!ptSpectra){ | |
427 | printf("BoltzmannFit: Can't get the normalized spectrum\n"); | |
428 | return; | |
429 | } | |
430 | ||
431 | printf("pt spectra get entries: %f\n",ptSpectra->GetEntries()); | |
432 | ||
433 | if(!ptSpectra->GetEntries()){ | |
434 | printf("BoltzmannFit: The fit data is empty!\n"); | |
435 | return; | |
436 | } | |
437 | ||
438 | gStyle->SetOptStat("neRM"); | |
439 | ||
440 | TVirtualFitter::SetDefaultFitter("Minuit2"); | |
441 | ||
442 | TF1 *ptFit = new TF1("ptFit",Boltzmann,0.2,4.0,2); | |
443 | ||
444 | gStyle->SetOptFit(1112); | |
445 | ||
446 | // ptFit->SetParName(0,"T"); | |
447 | // ptFit->SetParName(1,"betaS"); | |
448 | ||
449 | ptFit->SetParameters(0.1,0.5); //GeV | |
450 | ptFit->SetParLimits(0,0.001,0.5);//GeV | |
451 | ptFit->SetParLimits(1,0.0,1.0); | |
452 | ||
453 | ptSpectra->Fit("ptFit","Q R B ME I N"); | |
454 | ||
455 | timer.Stop(); | |
456 | timer.Print(); | |
457 | ||
458 | // ----------- print the fit results ---------- | |
459 | Double_t chi2 = ptFit->GetChisquare(); | |
460 | Double_t ndf = ptFit->GetNDF(); | |
461 | Double_t p1 = ptFit->GetParameter(0); | |
462 | Double_t param1Err = ptFit->GetParError(0); | |
463 | Double_t p2 = ptFit->GetParameter(1); | |
464 | Double_t param2Err = ptFit->GetParError(1); | |
465 | ||
466 | printf("chi2/ndf = %f/%f = %f\n", chi2,ndf,chi2/ndf); | |
467 | printf("p1 = %f (GeV)\tparam1Err = %f\n", p1, param1Err); | |
468 | printf("p2 = %f\tparam2Err = %f\n", p2, param2Err); | |
469 | ||
470 | ||
471 | TH1D *hPtExtBoltzmann = new TH1D("PtExtBoltzmann","Pt Corr Norm Ext",25,0.0,5.0); | |
472 | ||
473 | Double_t bin, binerr, test; | |
474 | ||
475 | for(Int_t i=0; i<ptSpectra->GetNbinsX()+1;i++){ | |
476 | ||
477 | bin = ptSpectra->GetBinContent(i); | |
478 | binerr = ptSpectra->GetBinError(i); | |
479 | test = ptSpectra->GetBinLowEdge(i); | |
480 | ||
481 | hPtExtBoltzmann->SetBinContent(i, bin); | |
482 | hPtExtBoltzmann->SetBinError(i, binerr); | |
483 | } | |
484 | ||
485 | Double_t eval; | |
486 | eval = ptFit->Eval(0.1); | |
487 | printf("extrapolated pt value = %g \n", eval); | |
488 | printf("********************************************\n\n"); | |
489 | ||
490 | hPtExtBoltzmann->SetBinContent(1, eval); | |
491 | ||
492 | hPtExtBoltzmann->SetMarkerStyle(24); | |
493 | hPtExtBoltzmann->SetMarkerColor(kRed); //ESD corrected | |
494 | hPtExtBoltzmann->GetXaxis()->SetTitle("Pt"); | |
495 | ||
496 | ||
497 | fResultsList->Add(hPtExtBoltzmann); | |
498 | ||
499 | } | |
500 | ||
501 |