q vector from loop on stack
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / MultEvShape / AliAnalysisCentralExtrapolate.cxx
CommitLineData
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
50ClassImp(AliAnalysisCentralExtrapolate)
51
52//________________________________________________________________________
53AliAnalysisCentralExtrapolate::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 69AliAnalysisCentralExtrapolate::~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<>
79static Double_t mass = 0.;
80static Double_t pt = 0.;
81static Double_t T = 0.;
82static Double_t betaS = 0.;
83static Double_t q = 0.;
84static Double_t mt = 0.;
85
86void 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
210Double_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
234Double_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
260void 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
358Double_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
389Double_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
408void 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