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