q vector from loop on stack
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / MultEvShape / AliAnalysisCentralExtrapolate.cxx
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