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