PWG2/SPECTRA/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);
 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 = 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