ANALYSIS/AliCentralitySelectionTask.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 //   Class AliCentralitySelectionTask
  18 //   Class to analyze determine centrality
  19 //   author: Alberica Toia
  20 //*****************************************************
  21
  22 #include <TTree.h>
  23 #include <TList.h>
  24 #include <TH1F.h>
  25 #include <TH2F.h>
  26 #include <TProfile.h>
  27 #include <TFile.h>
  28 #include <TObjString.h>
  29 #include <TString.h>
  30 #include <TCanvas.h>
  31 #include <TROOT.h>
  32 #include <TDirectory.h>
  33 #include <iostream>
  34
  35 #include "AliAnalysisManager.h"
  36 #include "AliVEvent.h"
  37 #include "AliESD.h"
  38 #include "AliESDEvent.h"
  39 #include "AliESDHeader.h"
  40 #include "AliESDInputHandler.h"
  41 #include "AliESDZDC.h"
  42 #include "AliESDFMD.h"
  43 #include "AliESDVZERO.h"
  44 #include "AliESDCentrality.h"
  45 #include "AliMultiplicity.h"
  46 #include "AliAODHandler.h"
  47 #include "AliAODEvent.h"
  48 #include "AliESDVertex.h"
  49 #include "AliAODVertex.h"
  50 #include "AliAODMCHeader.h"
  51 #include "AliMCEvent.h"
  52 #include "AliMCEventHandler.h"
  53 #include "AliMCParticle.h"
  54 #include "AliStack.h"
  55 #include "AliHeader.h"
  56 #include "AliAODMCParticle.h"
  57 #include "AliAnalysisTaskSE.h"
  58 #include "AliGenEventHeader.h"
  59 #include "AliGenHijingEventHeader.h"
  60 #include "AliPhysicsSelectionTask.h"
  61 #include "AliPhysicsSelection.h"
  62 #include "AliBackgroundSelection.h"
  63 #include "AliCentralitySelectionTask.h"
  64
  65 ClassImp(AliCentralitySelectionTask)
  66
  67
  68 //________________________________________________________________________
  69 AliCentralitySelectionTask::AliCentralitySelectionTask():
  70 AliAnalysisTaskSE(),
  71   fDebug(0),
  72   fAnalysisInput("ESD"),
  73   fIsMCInput(kFALSE),
  74   fFile(0),
  75   fFile2(0),
  76   fCentfilename(""),
  77   fCentfilename2(""),
  78   fFileList(new TList),
  79   fFileList2(new TList),
  80   fCurrentRun(-1),
  81   fRunNo(-1),
  82   fCentV0M(0),
  83   fCentFMD(0),
  84   fCentTRK(0),
  85   fCentTKL(0),
  86   fCentCL0(0),
  87   fCentCL1(0),
  88   fCentV0MvsFMD(0),
  89   fCentTKLvsV0M(0),
  90   fCentZEMvsZDC(0),
  91   fHtempV0M(0),
  92   fHtempFMD(0),
  93   fHtempTRK(0),
  94   fHtempTKL(0),
  95   fHtempCL0(0),
  96   fHtempCL1(0),
  97   fHtempV0MvsFMD(0),
  98   fHtempTKLvsV0M(0),
  99   fHtempZEMvsZDC(0)
 100 {
 101   // Default constructor
 102   fFileList->SetOwner();
 103   fFileList2->SetOwner();
 104
 105   AliInfo("Centrality Selection enabled.");
 106 }
 107
 108 //________________________________________________________________________
 109 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
 110   AliAnalysisTaskSE(name),
 111   fDebug(0),
 112   fAnalysisInput("ESD"),
 113   fIsMCInput(kFALSE),
 114   fFile(0),
 115   fFile2(0),
 116   fCentfilename(""),
 117   fCentfilename2(""),
 118   fFileList(new TList),
 119   fFileList2(new TList),
 120   fCurrentRun(-1),
 121   fRunNo(-1),
 122   fCentV0M(0),
 123   fCentFMD(0),
 124   fCentTRK(0),
 125   fCentTKL(0),
 126   fCentCL0(0),
 127   fCentCL1(0),
 128   fCentV0MvsFMD(0),
 129   fCentTKLvsV0M(0),
 130   fCentZEMvsZDC(0),
 131   fHtempV0M(0),
 132   fHtempFMD(0),
 133   fHtempTRK(0),
 134   fHtempTKL(0),
 135   fHtempCL0(0),
 136   fHtempCL1(0),
 137   fHtempV0MvsFMD(0),
 138   fHtempTKLvsV0M(0),
 139   fHtempZEMvsZDC(0)
 140 {
 141   // Default constructor
 142   fFileList->SetOwner();
 143   fFileList2->SetOwner();
 144
 145   AliInfo("Centrality Selection enabled.");
 146 }
 147
 148 //________________________________________________________________________
 149 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
 150 {
 151   // Assignment operator
 152   if (this!=&c) {
 153     AliAnalysisTaskSE::operator=(c);
 154   }
 155   return *this;
 156 }
 157
 158 //________________________________________________________________________
 159 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
 160   AliAnalysisTaskSE(ana),
 161   fDebug(ana.fDebug),
 162   fAnalysisInput(ana.fDebug),
 163   fIsMCInput(ana.fIsMCInput),
 164   fFile(ana.fFile),
 165   fFile2(ana.fFile2),
 166   fCentfilename(ana.fCentfilename),
 167   fCentfilename2(ana.fCentfilename2),
 168   fFileList(ana.fFileList),
 169   fFileList2(ana.fFileList2),
 170   fCurrentRun(ana.fCurrentRun),
 171   fRunNo(ana.fRunNo),
 172   fCentV0M(ana.fCentV0M),
 173   fCentFMD(ana.fCentFMD),
 174   fCentTRK(ana.fCentTRK),
 175   fCentTKL(ana.fCentTKL),
 176   fCentCL0(ana.fCentCL0),
 177   fCentCL1(ana.fCentCL1),
 178   fCentV0MvsFMD(ana.fCentV0MvsFMD),
 179   fCentTKLvsV0M(ana.fCentTKLvsV0M),
 180   fCentZEMvsZDC(ana.fCentZEMvsZDC),
 181   fHtempV0M(ana.fHtempV0M),
 182   fHtempFMD(ana.fHtempFMD),
 183   fHtempTRK(ana.fHtempTRK),
 184   fHtempTKL(ana.fHtempTKL),
 185   fHtempCL0(ana.fHtempCL0),
 186   fHtempCL1(ana.fHtempCL1),
 187   fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
 188   fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
 189   fHtempZEMvsZDC(ana.fHtempZEMvsZDC)
 190 {
 191   // Copy Constructor
 192 }
 193
 194 //________________________________________________________________________
 195 AliCentralitySelectionTask::~AliCentralitySelectionTask()
 196 {
 197   // Destructor
 198
 199   if (fFileList) {
 200     fFileList->Clear();
 201     delete fFileList;
 202   }
 203   fFileList = NULL;
 204
 205   if (fFileList2) {
 206     fFileList2->Clear();
 207     delete fFileList2;
 208   }
 209   fFileList2 = NULL;
 210 }
 211
 212 //________________________________________________________________________
 213 void AliCentralitySelectionTask::UserCreateOutputObjects()
 214 {
 215   // Create the output containers
 216   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
 217
 218   if (fFileList->GetEntries() < 1) {
 219     AliError("No Inputfiles Added");
 220   }
 221
 222   AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
 223 }
 224
 225 //________________________________________________________________________
 226 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
 227 {
 228   // Execute analysis for current event:
 229   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
 230
 231   Float_t  zncEnergy;          //  ZNC Energy
 232   Float_t  zpcEnergy;          //  ZPC Energy
 233   Float_t  znaEnergy;          //  ZNA Energy
 234   Float_t  zpaEnergy;          //  ZPA Energy
 235   Float_t  zem1Energy = 0.;         //  ZEM1 Energy
 236   Float_t  zem2Energy = 0.;         //  ZEM2 Energy
 237
 238   Int_t    nTracks = 0;             //  no. tracks
 239   Int_t    nTracklets = 0;          //  no. tracklets
 240   Int_t    nClusters[6];            //  no. clusters on 6 ITS layers
 241   Int_t    nChips[2];               //  no. chips on 2 SPD layers
 242   Float_t  spdCorr =0;              //  corrected spd2 multiplicity
 243
 244   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
 245   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
 246   Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A
 247   Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C
 248
 249   Short_t v0Corr = 0;               // corrected V0 multiplicity
 250   Short_t v0CorrResc = 0;           // corrected and rescaled V0 multiplicity
 251
 252   Float_t zvtx =0;                  // z-vertex SPD
 253
 254   AliESDCentrality *esdCent = 0;
 255
 256   if(fAnalysisInput.CompareTo("ESD")==0){
 257
 258     AliVEvent* event = InputEvent();
 259     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 260
 261     if (SetupRun(esd))
 262       return;
 263
 264     esdCent = esd->GetCentrality();
 265
 266     // ***** V0 info
 267     AliESDVZERO* esdV0 = esd->GetVZEROData();
 268     multV0A=esdV0->GetMTotV0A();
 269     multV0C=esdV0->GetMTotV0C();
 270
 271     float v0CorrR;
 272     v0Corr = (Short_t) (GetCorrV0(esd,v0CorrR));
 273     v0CorrResc = (Short_t)v0CorrR;
 274
 275     // ***** Vertex Info
 276     const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
 277     zvtx        = vtxESD->GetZ();
 278
 279     // ***** CB info (tracklets, clusters, chips)
 280     nTracks    = event->GetNumberOfTracks();
 281
 282     const AliMultiplicity *mult = esd->GetMultiplicity();
 283
 284     nTracklets = mult->GetNumberOfTracklets();
 285
 286     for(Int_t ilay=0; ilay<6; ilay++){
 287       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
 288     }
 289
 290     for(Int_t ilay=0; ilay<2; ilay++){
 291       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
 292     }
 293
 294     spdCorr = GetCorrSPD2(nClusters[1],zvtx);
 295
 296     // ***** FMD info
 297     AliESDFMD *fmd = esd->GetFMDData();
 298     Float_t totalMultA = 0;
 299     Float_t totalMultC = 0;
 300     const Float_t fFMDLowCut = 0.4;
 301
 302     for(UShort_t det=1;det<=3;det++) {
 303       Int_t nRings = (det==1 ? 1 : 2);
 304       for (UShort_t ir = 0; ir < nRings; ir++) {
 305         Char_t   ring = (ir == 0 ? 'I' : 'O');
 306         UShort_t nsec = (ir == 0 ? 20  : 40);
 307         UShort_t nstr = (ir == 0 ? 512 : 256);
 308         for(UShort_t sec =0; sec < nsec;  sec++)  {
 309           for(UShort_t strip = 0; strip < nstr; strip++) {
 310
 311             Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
 312             if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
 313
 314             Float_t nParticles=0;
 315
 316             if(FMDmult > fFMDLowCut) {
 317               nParticles = 1.;
 318             }
 319
 320             if (det<3) totalMultA = totalMultA + nParticles;
 321             else totalMultC = totalMultC + nParticles;
 322
 323           }
 324         }
 325       }
 326     }
 327     multFMDA = totalMultA;
 328     multFMDC = totalMultC;
 329
 330     // ***** ZDC info
 331     AliESDZDC *esdZDC = esd->GetESDZDC();
 332     zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
 333     zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
 334     znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
 335     zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
 336     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
 337     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
 338
 339   }
 340   else if(fAnalysisInput.CompareTo("AOD")==0){
 341     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
 342     // to be implemented
 343     printf("  AOD analysis not yet implemented!!!\n\n");
 344     return;
 345   }
 346
 347   // ***** Centrality Selection
 348   if(fHtempV0M)  fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
 349   ///  else     printf("  Centrality by V0 not available!!!\n\n");
 350   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
 351   //  else     printf("  Centrality by FMD not available!!!\n\n");
 352   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
 353   //  else     printf("  Centrality by TRK not available!!!\n\n");
 354   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
 355   //  else     printf("  Centrality by TKL not available!!!\n\n");
 356   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
 357   //  else     printf("  Centrality by CL0 not available!!!\n\n");
 358   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
 359   ///  else     printf("  Centrality by CL1 not available!!!\n\n");
 360
 361   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
 362   //  else     printf("  Centrality by V0 vs FMD not available!!!\n\n");
 363   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
 364   //  else     printf("  Centrality by V0 vs TKL not available!!!\n\n");
 365   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));
 366   //  else     printf("  Centrality by ZEM vs ZDC not available!!!\n\n");
 367
 368   esdCent->SetCentralityV0M(fCentV0M);
 369   esdCent->SetCentralityFMD(fCentFMD);
 370   esdCent->SetCentralityTRK(fCentTRK);
 371   esdCent->SetCentralityTKL(fCentTKL);
 372   esdCent->SetCentralityCL0(fCentCL0);
 373   esdCent->SetCentralityCL1(fCentCL1);
 374   esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
 375   esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
 376   esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
 377 }
 378
 379 //________________________________________________________________________
 380 void AliCentralitySelectionTask::ReadCentralityHistos()
 381 {
 382 //  Read centrality histograms
 383     TDirectory *owd = gDirectory;
 384     fFile  = TFile::Open(fCentfilename);
 385     owd->cd();
 386     fHtempV0M  = (TH1F*) (fFile->Get("hmultV0_percentile"));
 387     fHtempFMD  = (TH1F*) (fFile->Get("hmultFMD_percentile"));
 388     fHtempTRK  = (TH1F*) (fFile->Get("hNtracks_percentile"));
 389     fHtempTKL  = (TH1F*) (fFile->Get("hNtracklets_percentile"));
 390     fHtempCL0  = (TH1F*) (fFile->Get("hNclusters0_percentile"));
 391     fHtempCL1  = (TH1F*) (fFile->Get("hNclusters1_percentile"));
 392 }
 393
 394 //________________________________________________________________________
 395 void AliCentralitySelectionTask::ReadCentralityHistos2()
 396 {
 397 //  Read centrality histograms
 398     TDirectory *owd = gDirectory;
 399     fFile2  = TFile::Open(fCentfilename2);
 400     owd->cd();
 401     fHtempV0MvsFMD =  (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
 402     fHtempTKLvsV0M  = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
 403     fHtempZEMvsZDC  = (TH1F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
 404 }
 405
 406 //________________________________________________________________________
 407 void AliCentralitySelectionTask::SetPercentileFile(TString filename)
 408 {
 409 // Set the percentile file name
 410   fCentfilename = filename;
 411   ReadCentralityHistos();
 412 }
 413
 414 //________________________________________________________________________
 415 void AliCentralitySelectionTask::SetPercentileFile2(TString filename)
 416 {
 417 // Set the percentile file name
 418   fCentfilename2 = filename;
 419   ReadCentralityHistos2();
 420 }
 421
 422 //________________________________________________________________________
 423 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
 424 {
 425   // Terminate analysis
 426   if (fFile && fFile->IsOpen())
 427     fFile->Close();
 428   if (fFile2 && fFile2->IsOpen())
 429     fFile2->Close();
 430 }
 431 //________________________________________________________________________
 432 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* esd)
 433 {
 434   // Setup files for run
 435
 436   if (!esd)
 437     return -1;
 438
 439   // check if something to be done
 440   if (fCurrentRun == esd->GetRunNumber())
 441     return 0;
 442   else
 443     fCurrentRun = esd->GetRunNumber();
 444
 445   AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
 446
 447   fRunNo = fCurrentRun;
 448
 449   // CHANGE HERE FOR RUN RANGES
 450   if ( fRunNo == 137162 ) fRunNo = 137161;
 451   else if ( fRunNo == 137365 ) fRunNo = 137366;
 452   else if ( fRunNo > 137366 ) fRunNo = 137366;
 453   // CHANGE HERE FOR RUN RANGES
 454
 455   TString runName(Form("%d", fRunNo));
 456   TString fileName("");
 457   Bool_t isRunKnown = kFALSE;
 458
 459   // Check if run is in fileList
 460   // if not, take the last name in the list
 461   for ( Int_t idx=0 ; idx < fFileList->GetEntries(); ++idx ) {
 462
 463     TString str((dynamic_cast<TObjString*>(fFileList->At(idx)))->GetString());
 464     if (str.Contains(runName)) {
 465       fileName += str;
 466       isRunKnown = kTRUE;
 467       break;
 468     }
 469   }
 470
 471   if (!isRunKnown) {
 472     if (fFileList->Last()) {
 473       fileName += (dynamic_cast<TObjString*>(fFileList->Last()))->GetString();
 474       AliError(Form("Run %d not known to centrality selection!", fCurrentRun));
 475     }
 476   }
 477
 478   if (fileName.Contains(".root")) {
 479     AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data()));
 480     SetPercentileFile(fileName.Data());
 481     return 0;
 482   }
 483
 484   return -1;
 485 }
 486 //________________________________________________________________________
 487 Float_t AliCentralitySelectionTask::GetCorrV0(const AliESDEvent* esd, float &v0CorrResc)
 488 {
 489   // correct V0 non-linearity, prepare a version rescaled to SPD2 corr
 490   Double_t *par0;
 491   Double_t *par1;
 492   Double_t *par2;
 493
 494   Double_t par0_137161[64] = { 6.71e-02 , 6.86e-02 , 7.06e-02 , 6.32e-02 ,
 495                         5.91e-02 , 6.07e-02 , 5.78e-02 , 5.73e-02 , 5.91e-02 , 6.22e-02 ,
 496                         5.90e-02 , 6.11e-02 , 5.55e-02 , 5.29e-02 , 5.19e-02 , 5.56e-02 ,
 497                         6.25e-02 , 7.03e-02 , 5.64e-02 , 5.81e-02 , 4.57e-02 , 5.30e-02 ,
 498                         5.13e-02 , 6.43e-02 , 6.27e-02 , 6.48e-02 , 6.07e-02 , 1.01e-01 ,
 499                         6.68e-02 , 7.16e-02 , 6.36e-02 , 5.95e-02 , 2.52e-02 , 2.82e-02 ,
 500                         2.56e-02 , 2.86e-02 , 2.82e-02 , 2.10e-02 , 2.13e-02 , 2.32e-02 ,
 501                         2.75e-02 , 4.34e-02 , 3.78e-02 , 4.52e-02 , 4.11e-02 , 3.89e-02 ,
 502                         4.10e-02 , 3.73e-02 , 4.51e-02 , 5.07e-02 , 5.42e-02 , 4.74e-02 ,
 503                         4.33e-02 , 4.44e-02 , 4.64e-02 , 3.01e-02 , 6.38e-02 , 5.26e-02 ,
 504                         4.99e-02 , 5.26e-02 , 5.47e-02 , 3.84e-02 , 5.00e-02 , 5.20e-02 };
 505   Double_t par1_137161[64] = { -6.68e-05 , -7.78e-05 , -6.88e-05 , -5.92e-05 ,
 506                         -2.43e-05 , -3.54e-05 , -2.91e-05 , -1.99e-05 , -1.40e-05 , -4.01e-05 ,
 507                         -2.29e-05 , -3.68e-05 , -2.53e-05 , -2.44e-06 , -9.22e-06 , -1.51e-05 ,
 508                         -2.80e-05 , -2.34e-05 , -1.72e-05 , -1.81e-05 , -1.29e-05 , -2.65e-05 ,
 509                         -1.61e-05 , -2.86e-05 , -1.74e-05 , -4.23e-05 , -3.41e-05 , -1.05e-04 ,
 510                         -2.76e-05 , -4.71e-05 , -3.06e-05 , -2.32e-05 , -1.55e-06 , 2.15e-05 ,
 511                         1.40e-05 , 2.16e-05 , 1.21e-05 , 3.05e-06 , 1.67e-05 , -3.84e-06 ,
 512                         3.09e-06 , 1.50e-05 , 3.47e-06 , 4.87e-06 , -3.71e-07 , -1.75e-06 ,
 513                         -1.80e-06 , 9.99e-06 , -6.46e-06 , -4.91e-06 , 1.33e-05 , -2.52e-07 ,
 514                         -3.85e-06 , 4.94e-06 , -2.48e-07 , -1.20e-05 , 2.07e-06 , 6.12e-06 ,
 515                         -1.18e-06 , 4.54e-06 , -1.54e-05 , -1.25e-05 , 1.46e-06 , -6.67e-06 };
 516   Double_t par2_137161[64] = { 1.29e-08 , 1.51e-08 , 1.43e-08 , 1.11e-08 ,
 517                         5.04e-09 , 6.99e-09 , 5.58e-09 , 4.15e-09 , 4.00e-09 , 8.22e-09 ,
 518                         4.97e-09 , 7.66e-09 , 4.91e-09 , 1.10e-09 , 2.64e-09 , 3.64e-09 ,
 519                         5.76e-09 , 5.46e-09 , 3.38e-09 , 3.47e-09 , 2.43e-09 , 4.13e-09 ,
 520                         2.80e-09 , 5.80e-09 , 3.86e-09 , 7.46e-09 , 5.98e-09 , 2.58e-08 ,
 521                         5.50e-09 , 8.72e-09 , 5.23e-09 , 4.37e-09 , 2.33e-09 , -6.01e-10 ,
 522                         3.99e-11 , -2.02e-10 , 7.67e-10 , 2.03e-09 , 1.17e-10 , 2.56e-09 ,
 523                         1.16e-09 , -4.75e-10 , 1.28e-09 , 1.23e-09 , 1.62e-09 , 1.61e-09 ,
 524                         1.93e-09 , 2.97e-10 , 2.21e-09 , 2.16e-09 , 5.22e-10 , 1.03e-09 ,
 525                         1.56e-09 , 5.00e-10 , 1.01e-09 , 2.93e-09 , 1.05e-09 , 9.96e-11 ,
 526                         1.21e-09 , 7.45e-10 , 3.07e-09 , 2.31e-09 , 6.70e-10 , 1.89e-09 };
 527
 528   Double_t par0_137366[64] = { 7.12e-02 , 7.34e-02 , 7.39e-02 , 6.54e-02 , 6.11e-02 , 6.31e-02 , 6.15e-02 ,
 529                                6.00e-02 , 6.10e-02 , 6.49e-02 , 6.17e-02 , 6.33e-02 , 6.00e-02 , 5.48e-02 ,
 530                                5.44e-02 , 5.81e-02 , 6.49e-02 , 7.07e-02 , 5.91e-02 , 6.18e-02 , 4.82e-02 ,
 531                                5.67e-02 , 5.36e-02 , 6.60e-02 , 6.37e-02 , 6.78e-02 , 6.31e-02 , 1.04e-01 ,
 532                                6.91e-02 , 7.32e-02 , 6.61e-02 , 6.16e-02 , 2.64e-02 , 2.81e-02 , 2.64e-02 ,
 533                                2.85e-02 , 2.87e-02 , 2.18e-02 , 2.19e-02 , 2.43e-02 , 2.81e-02 , 4.37e-02 ,
 534                                3.90e-02 , 4.66e-02 , 4.24e-02 , 4.09e-02 , 4.21e-02 , 3.88e-02 , 4.83e-02 ,
 535                                5.23e-02 , 5.44e-02 , 4.85e-02 , 4.42e-02 , 4.58e-02 , 4.74e-02 , 3.14e-02 ,
 536                                6.31e-02 , 5.30e-02 , 5.01e-02 , 5.33e-02 , 5.70e-02 , 3.95e-02 , 4.98e-02 , 5.31e-02 };
 537   Double_t par1_137366[64] = { -6.99e-05 , -6.99e-05 , -6.94e-05 , -6.55e-05 , -3.55e-05 , -4.50e-05 ,
 538                                -3.10e-05 , -2.81e-05 , -2.29e-05 , -3.89e-05 , -2.53e-05 , -4.25e-05 ,
 539                                -1.87e-05 , -2.01e-05 , -1.53e-05 , -2.14e-05 , -2.86e-05 , -4.70e-05 ,
 540                                -2.23e-05 , -3.30e-05 ,-9.74e-06 , -2.62e-05 , -1.76e-05 , -2.38e-05 ,
 541                                -2.40e-05 , -3.43e-05 , -2.75e-05 , -6.86e-05 ,-2.35e-05 , -4.45e-05 ,
 542                                -2.51e-05 , -2.20e-05 , -1.25e-16 , -2.04e-17 , -2.06e-17 , -3.74e-19 ,
 543                                -1.18e-18 , -2.02e-15 , -3.78e-06 , -1.26e-06 , -2.71e-06 , -6.23e-17 ,
 544                                -7.39e-08 , -1.76e-16 , -8.98e-06 , -4.10e-18 , -1.34e-05 , -1.06e-16 ,
 545                                -3.34e-06 , -1.04e-05 , -5.28e-06 , -7.34e-06 , -1.05e-05 , -7.68e-06 ,
 546                                -1.78e-05 , -1.19e-05 , -1.78e-05 , -1.34e-06 , -9.23e-06 , -3.34e-06 ,
 547                                -8.02e-06 , -1.39e-05 , -1.38e-05 , -1.40e-05 };
 548   Double_t par2_137366[64] = { 1.41e-08 , 1.47e-08 , 1.48e-08 , 1.24e-08 , 6.82e-09 , 8.73e-09 , 6.26e-09 ,
 549                                5.53e-09 , 5.40e-09 , 7.93e-09 , 5.49e-09 , 8.77e-09 , 4.21e-09 , 3.93e-09 ,
 550                                3.60e-09 , 4.67e-09 , 5.59e-09 , 8.81e-09 , 3.89e-09 , 6.19e-09 , 1.97e-09 ,
 551                                4.38e-09 , 3.26e-09 , 5.00e-09 , 4.58e-09 , 6.39e-09 , 5.03e-09 , 1.30e-08 ,
 552                                4.95e-09 , 8.26e-09 , 4.57e-09 , 4.10e-09 , 2.35e-09 , 2.30e-09 , 2.15e-09 ,
 553                                2.27e-09 , 2.17e-09 , 2.27e-09 , 2.97e-09 , 2.25e-09 , 1.69e-09 , 1.44e-09 ,
 554                                1.66e-09 , 1.75e-09 , 2.88e-09 , 1.82e-09 , 3.64e-09 , 1.80e-09 , 1.71e-09 ,
 555                                2.66e-09 , 3.01e-09 , 1.95e-09 , 2.64e-09 , 2.42e-09 , 3.68e-09 , 2.66e-09 ,
 556                                3.92e-09 , 1.18e-09 , 2.26e-09 , 1.57e-09 , 2.02e-09 , 2.71e-09 , 2.99e-09 , 3.04e-09 };
 557
 558
 559   if (esd->GetRunNumber() <= 137165) {
 560     par0=par0_137161;
 561     par1=par1_137161;
 562     par2=par2_137161;
 563   }  else  {
 564     par0=par0_137366;
 565     par1=par1_137366;
 566     par2=par2_137366;
 567  }
 568   //
 569   Float_t multCorr = 0;
 570   Float_t multCorr2 = 0;
 571   Float_t multChCorr[64];
 572   AliESDVZERO* esdV0 = esd->GetVZEROData();
 573   for(Int_t i = 0; i < 64; ++i) {
 574     Double_t b = (esdV0->GetMultiplicity(i)*par1[i]-par0[i]);
 575     Double_t s = (b*b-4.*par2[i]*esdV0->GetMultiplicity(i)*esdV0->GetMultiplicity(i));
 576     Double_t n = (s<0) ? -b : (-b + TMath::Sqrt(s));
 577     multChCorr[i] = 2.*esdV0->GetMultiplicity(i)/n*par0[i];
 578     multCorr += multChCorr[i];
 579     multCorr2 += (multChCorr[i]/par0[i]/64.);
 580   }
 581   v0CorrResc =  multCorr2;
 582   return multCorr;
 583 }
 584
 585 //____________________________________________________________________
 586 Float_t AliCentralitySelectionTask::GetCorrSPD2(Float_t spd2raw,Float_t zv)
 587 {
 588   // renormalize N spd2 clusters at given Zv to acceptance at Zv=0
 589   const double pars[] = {8.10030e-01,-2.80364e-03,-7.19504e-04};
 590   zv -= pars[0];
 591   float corr = 1 + zv*(pars[1] + zv*pars[2]);
 592   return corr>0 ? spd2raw/corr : -1;
 593 }