]>
Commit | Line | Data |
---|---|---|
bde49c8a | 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 | ||
410ff8cc | 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 | ||
bde49c8a | 25 | |
26 | #include "TH1D.h" | |
27 | #include "TFile.h" | |
28 | #include "TList.h" | |
410ff8cc | 29 | #include "TCanvas.h" |
30 | #include "TLegend.h" | |
bde49c8a | 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() | |
410ff8cc | 56 | ,fInputList(0) |
bde49c8a | 57 | ,fResultsList(0x0) |
58 | ||
59 | { | |
60 | // Constructor | |
61 | ||
011e988c | 62 | Printf("Running %s!\n", name); |
63 | ||
410ff8cc | 64 | fInputList = new TList(); |
bde49c8a | 65 | fResultsList = new TList(); |
011e988c | 66 | |
bde49c8a | 67 | } |
011e988c | 68 | |
bde49c8a | 69 | AliAnalysisCentralExtrapolate::~AliAnalysisCentralExtrapolate() { |
011e988c | 70 | |
bde49c8a | 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; | |
410ff8cc | 94 | |
95 | ||
bde49c8a | 96 | TFile *file1 = new TFile("$ALICE_ROOT/PWG2/data/AliAnalysisCentralEfficiency.root", "read"); |
97 | ||
98 | AliCFEffGrid *eff = 0; | |
99 | ||
410ff8cc | 100 | if(fPartType.Contains("kPi")){ |
101 | data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_Pi")); | |
bde49c8a | 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")){ | |
410ff8cc | 109 | data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_K")); |
bde49c8a | 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")){ | |
410ff8cc | 117 | data = dynamic_cast<AliCFContainer*>(fInputList->FindObject("TaskCentral_CFCont_P")); |
bde49c8a | 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"); | |
8a993eac | 128 | return; |
bde49c8a | 129 | } |
130 | ||
131 | if(!eff){ | |
132 | printf("No Eff Grid found! \n"); | |
8a993eac | 133 | return; |
bde49c8a | 134 | } |
410ff8cc | 135 | |
136 | // TCanvas *ccorrdata = new TCanvas(); | |
bde49c8a | 137 | // ccorrdata->Divide(1,2); |
138 | // ccorrdata->cd(1); | |
410ff8cc | 139 | // ccorrdata->cd(1)->SetLogy(); |
bde49c8a | 140 | |
141 | ||
67a4103d | 142 | AliCFDataGrid *corrdata = new AliCFDataGrid("corrdata","corrected data",*data, stepRec); |
bde49c8a | 143 | |
144 | //correct selection step "reconstructed" | |
67a4103d | 145 | // corrdata->SetMeasured(stepRec); //set data to be corrected |
bde49c8a | 146 | corrdata->ApplyEffCorrection(*eff);//apply the correction for efficiency |
147 | ||
410ff8cc | 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}}"); | |
bde49c8a | 153 | // hPtMC->Draw("p e1"); |
410ff8cc | 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"); | |
bde49c8a | 159 | // hPtESDI->Draw("p e1 same"); |
160 | ||
7faab71b | 161 | TH1D *hPtESDCorr = (TH1D*)corrdata->Project(0); //corrected data |
bde49c8a | 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(); | |
410ff8cc | 176 | // |
177 | // TH1D *efficiency = eff->Project(0); //the efficiency vs pt (ipt = 0) | |
178 | // efficiency->GetXaxis()->SetTitle("p_{T}(GeV/c)"); | |
bde49c8a | 179 | // efficiency->Draw(); |
180 | ||
181 | ||
0fb1f0cf | 182 | TH1D *hNoEvt = dynamic_cast<TH1D*>(fInputList->FindObject("TaskCentral_NoEvt")); |
bde49c8a | 183 | if(!hNoEvt){ |
184 | printf("Unable to get the number of events! \n"); | |
8a993eac | 185 | return; |
bde49c8a | 186 | } |
187 | ||
410ff8cc | 188 | |
bde49c8a | 189 | Int_t noEvt = (Int_t)(hNoEvt->GetEntries()); |
190 | ||
191 | printf("\n** No of processed events = %i **\n",noEvt); | |
8a993eac | 192 | |
bde49c8a | 193 | Double_t scale = 1.0/noEvt; |
8a993eac | 194 | |
bde49c8a | 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")); | |
8a993eac | 277 | if(!ptSpectra){ |
410ff8cc | 278 | printf("TsallisFit: Can't get the normalized spectrum\n"); |
8a993eac | 279 | return; |
280 | } | |
281 | ||
282 | if(!ptSpectra->GetEntries()){ | |
410ff8cc | 283 | printf("TsallisFit: The fit data is empty!\n"); |
8a993eac | 284 | return; |
285 | } | |
bde49c8a | 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 | ||
410ff8cc | 293 | // ptFit->SetParName(0,"T"); |
294 | // ptFit->SetParName(1,"betaS"); | |
295 | // ptFit->SetParName(2,"q"); | |
bde49c8a | 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 | |
0fb1f0cf | 324 | TH1D *hPtExtTsallis = new TH1D("PtExtTsallis","Pt Corr Norm Ext",25,0.0,5.0); |
bde49c8a | 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 | ||
410ff8cc | 365 | |
bde49c8a | 366 | Double_t rho = TMath::ATanH(betaR); |
367 | ||
410ff8cc | 368 | |
011e988c | 369 | mt = sqrt(mass*mass + pt*pt);// !!! par[0]= pt !!!!! |
410ff8cc | 370 | |
bde49c8a | 371 | Double_t argI0 = (pt*TMath::SinH(rho))/T; //T = par[1] |
410ff8cc | 372 | if(argI0>700.0) return 0.0; // !! floating point exception protection |
bde49c8a | 373 | |
374 | Double_t argK1 = (mt*TMath::CosH(rho))/T; | |
375 | ||
410ff8cc | 376 | |
bde49c8a | 377 | Double_t i0 = TMath::BesselI0(argI0); |
378 | ||
410ff8cc | 379 | |
bde49c8a | 380 | Double_t k1 = TMath::BesselK1(argK1); |
381 | ||
bde49c8a | 382 | |
410ff8cc | 383 | Double_t f = r*mt*i0*k1; |
384 | ||
bde49c8a | 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); | |
8a993eac | 421 | |
bde49c8a | 422 | timer.Start(); |
423 | ||
424 | ||
425 | TH1D *ptSpectra = dynamic_cast<TH1D*>(fResultsList->FindObject("hPtESDCorrNorm")); | |
8a993eac | 426 | if(!ptSpectra){ |
427 | printf("BoltzmannFit: Can't get the normalized spectrum\n"); | |
428 | return; | |
429 | } | |
430 | ||
410ff8cc | 431 | printf("pt spectra get entries: %f\n",ptSpectra->GetEntries()); |
432 | ||
8a993eac | 433 | if(!ptSpectra->GetEntries()){ |
434 | printf("BoltzmannFit: The fit data is empty!\n"); | |
435 | return; | |
436 | } | |
bde49c8a | 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 | ||
410ff8cc | 446 | // ptFit->SetParName(0,"T"); |
447 | // ptFit->SetParName(1,"betaS"); | |
bde49c8a | 448 | |
449 | ptFit->SetParameters(0.1,0.5); //GeV | |
410ff8cc | 450 | ptFit->SetParLimits(0,0.001,0.5);//GeV |
bde49c8a | 451 | ptFit->SetParLimits(1,0.0,1.0); |
410ff8cc | 452 | |
bde49c8a | 453 | ptSpectra->Fit("ptFit","Q R B ME I N"); |
454 | ||
410ff8cc | 455 | timer.Stop(); |
bde49c8a | 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 | ||
0fb1f0cf | 471 | TH1D *hPtExtBoltzmann = new TH1D("PtExtBoltzmann","Pt Corr Norm Ext",25,0.0,5.0); |
bde49c8a | 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 |