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