ANALYSIS/AliCentralitySelectionTask.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 \r
  16 //*****************************************************\r
  17 //   Class AliCentralitySelectionTask\r
  18 //   Class to analyze determine centrality            \r
  19 //   author: Alberica Toia\r
  20 //*****************************************************\r
  21 \r
  22 #include "AliCentralitySelectionTask.h"\r
  23 \r
  24 #include <TTree.h>\r
  25 #include <TList.h>\r
  26 #include <TH1F.h>\r
  27 #include <TH2F.h>\r
  28 #include <TProfile.h>\r
  29 #include <TFile.h>\r
  30 #include <TObjString.h>\r
  31 #include <TString.h>\r
  32 #include <TCanvas.h>\r
  33 #include <TROOT.h>\r
  34 #include <TDirectory.h>\r
  35 #include <TSystem.h>\r
  36 #include <iostream>\r
  37 \r
  38 #include "AliAnalysisManager.h"\r
  39 #include "AliVEvent.h"\r
  40 #include "AliESD.h"\r
  41 #include "AliESDEvent.h"\r
  42 #include "AliESDHeader.h"\r
  43 #include "AliESDInputHandler.h"\r
  44 #include "AliESDZDC.h"\r
  45 #include "AliESDFMD.h"\r
  46 #include "AliESDVZERO.h"\r
  47 #include "AliCentrality.h"\r
  48 #include "AliESDtrackCuts.h"\r
  49 #include "AliMultiplicity.h"\r
  50 #include "AliAODHandler.h"\r
  51 #include "AliAODEvent.h"\r
  52 #include "AliESDVertex.h"\r
  53 #include "AliAODVertex.h"\r
  54 #include "AliAODMCHeader.h"\r
  55 #include "AliMCEvent.h"\r
  56 #include "AliMCEventHandler.h"\r
  57 #include "AliMCParticle.h"\r
  58 #include "AliStack.h"\r
  59 #include "AliHeader.h"\r
  60 #include "AliAODMCParticle.h"\r
  61 #include "AliAnalysisTaskSE.h"\r
  62 #include "AliGenEventHeader.h"\r
  63 #include "AliGenHijingEventHeader.h"\r
  64 #include "AliPhysicsSelectionTask.h"\r
  65 #include "AliPhysicsSelection.h"\r
  66 #include "AliBackgroundSelection.h"\r
  67 #include "AliESDUtils.h"\r
  68 \r
  69 ClassImp(AliCentralitySelectionTask)\r
  70 \r
  71 \r
  72 //________________________________________________________________________\r
  73 AliCentralitySelectionTask::AliCentralitySelectionTask():\r
  74 AliAnalysisTaskSE(),\r
  75   fDebug(0),\r
  76   fAnalysisInput("ESD"),\r
  77   fIsMCInput(kFALSE),\r
  78   fFile(0),\r
  79   fFile2(0),\r
  80   fCurrentRun(-1),\r
  81   fRunNo(-1),\r
  82   fTrackCuts(0),\r
  83   fCentV0M(0),\r
  84   fCentFMD(0),\r
  85   fCentTRK(0),\r
  86   fCentTKL(0),\r
  87   fCentCL0(0),\r
  88   fCentCL1(0),\r
  89   fCentV0MvsFMD(0),\r
  90   fCentTKLvsV0M(0),\r
  91   fCentZEMvsZDC(0),\r
  92   fHtempV0M(0),\r
  93   fHtempFMD(0),\r
  94   fHtempTRK(0),\r
  95   fHtempTKL(0),\r
  96   fHtempCL0(0),\r
  97   fHtempCL1(0),\r
  98   fHtempV0MvsFMD(0),\r
  99   fHtempTKLvsV0M(0),\r
 100   fHtempZEMvsZDC(0),\r
 101   fOutputList(0),\r
 102   fHOutCentV0M     (0),\r
 103   fHOutCentFMD     (0),\r
 104   fHOutCentTRK     (0),\r
 105   fHOutCentTKL     (0),\r
 106   fHOutCentCL0     (0),\r
 107   fHOutCentCL1     (0),\r
 108   fHOutCentV0MvsFMD(0),\r
 109   fHOutCentTKLvsV0M(0),\r
 110   fHOutCentZEMvsZDC(0),\r
 111   fHOutMultV0M(0),\r
 112   fHOutMultV0R(0),\r
 113   fHOutMultFMD(0),\r
 114   fHOutMultTRK(0),\r
 115   fHOutMultTKL(0),\r
 116   fHOutMultCL0(0),\r
 117   fHOutMultCL1(0),\r
 118   fHOutMultV0MvsZDC(0),\r
 119   fHOutMultZEMvsZDC(0),\r
 120   fHOutMultV0MvsCL1(0),\r
 121   fHOutMultV0MvsTRK(0),\r
 122   fHOutMultTRKvsCL1(0)\r
 123 {   \r
 124   // Default constructor\r
 125   AliInfo("Centrality Selection enabled.");\r
 126 }   \r
 127 \r
 128 //________________________________________________________________________\r
 129 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):\r
 130   AliAnalysisTaskSE(name),\r
 131   fDebug(0),\r
 132   fAnalysisInput("ESD"),\r
 133   fIsMCInput(kFALSE),\r
 134   fFile(0),\r
 135   fFile2(0),\r
 136   fCurrentRun(-1),\r
 137   fRunNo(-1),\r
 138   fTrackCuts(0),\r
 139   fCentV0M(0),\r
 140   fCentFMD(0),\r
 141   fCentTRK(0),\r
 142   fCentTKL(0),\r
 143   fCentCL0(0),\r
 144   fCentCL1(0),\r
 145   fCentV0MvsFMD(0),\r
 146   fCentTKLvsV0M(0),\r
 147   fCentZEMvsZDC(0),\r
 148   fHtempV0M(0),\r
 149   fHtempFMD(0),\r
 150   fHtempTRK(0),\r
 151   fHtempTKL(0),\r
 152   fHtempCL0(0),\r
 153   fHtempCL1(0),\r
 154   fHtempV0MvsFMD(0),\r
 155   fHtempTKLvsV0M(0),\r
 156   fHtempZEMvsZDC(0),\r
 157   fOutputList(0),\r
 158   fHOutCentV0M     (0),\r
 159   fHOutCentFMD     (0),\r
 160   fHOutCentTRK     (0),\r
 161   fHOutCentTKL     (0),\r
 162   fHOutCentCL0     (0),\r
 163   fHOutCentCL1     (0),\r
 164   fHOutCentV0MvsFMD(0),\r
 165   fHOutCentTKLvsV0M(0),\r
 166   fHOutCentZEMvsZDC(0),\r
 167   fHOutMultV0M(0),\r
 168   fHOutMultV0R(0),\r
 169   fHOutMultFMD(0),\r
 170   fHOutMultTRK(0),\r
 171   fHOutMultTKL(0),\r
 172   fHOutMultCL0(0),\r
 173   fHOutMultCL1(0),\r
 174   fHOutMultV0MvsZDC(0),\r
 175   fHOutMultZEMvsZDC(0),\r
 176   fHOutMultV0MvsCL1(0),\r
 177   fHOutMultV0MvsTRK(0),\r
 178   fHOutMultTRKvsCL1(0)\r
 179 {\r
 180   // Default constructor\r
 181   AliInfo("Centrality Selection enabled.");\r
 182   DefineOutput(1, TList::Class());\r
 183 }\r
 184 \r
 185 //________________________________________________________________________\r
 186 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)\r
 187 {\r
 188   // Assignment operator\r
 189   if (this!=&c) {\r
 190     AliAnalysisTaskSE::operator=(c);\r
 191   }\r
 192   return *this;\r
 193 }\r
 194 \r
 195 //________________________________________________________________________\r
 196 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):\r
 197   AliAnalysisTaskSE(ana),\r
 198   fDebug(ana.fDebug),     \r
 199   fAnalysisInput(ana.fDebug),\r
 200   fIsMCInput(ana.fIsMCInput),\r
 201   fFile(ana.fFile),\r
 202   fFile2(ana.fFile2),\r
 203   fCurrentRun(ana.fCurrentRun),\r
 204   fRunNo(ana.fRunNo),\r
 205   fTrackCuts(ana.fTrackCuts),\r
 206   fCentV0M(ana.fCentV0M),\r
 207   fCentFMD(ana.fCentFMD),\r
 208   fCentTRK(ana.fCentTRK),\r
 209   fCentTKL(ana.fCentTKL),\r
 210   fCentCL0(ana.fCentCL0),\r
 211   fCentCL1(ana.fCentCL1),\r
 212   fCentV0MvsFMD(ana.fCentV0MvsFMD),\r
 213   fCentTKLvsV0M(ana.fCentTKLvsV0M),\r
 214   fCentZEMvsZDC(ana.fCentZEMvsZDC),\r
 215   fHtempV0M(ana.fHtempV0M),\r
 216   fHtempFMD(ana.fHtempFMD),\r
 217   fHtempTRK(ana.fHtempTRK),\r
 218   fHtempTKL(ana.fHtempTKL),\r
 219   fHtempCL0(ana.fHtempCL0),\r
 220   fHtempCL1(ana.fHtempCL1),\r
 221   fHtempV0MvsFMD(ana.fHtempV0MvsFMD),\r
 222   fHtempTKLvsV0M(ana.fHtempTKLvsV0M),\r
 223   fHtempZEMvsZDC(ana.fHtempZEMvsZDC),\r
 224   fOutputList(ana.fOutputList),\r
 225   fHOutCentV0M     (ana.fHOutCentV0M     ),\r
 226   fHOutCentFMD     (ana.fHOutCentFMD     ),\r
 227   fHOutCentTRK     (ana.fHOutCentTRK     ),\r
 228   fHOutCentTKL     (ana.fHOutCentTKL     ),\r
 229   fHOutCentCL0     (ana.fHOutCentCL0     ),\r
 230   fHOutCentCL1     (ana.fHOutCentCL1     ),\r
 231   fHOutCentV0MvsFMD(ana.fHOutCentV0MvsFMD),\r
 232   fHOutCentTKLvsV0M(ana.fHOutCentTKLvsV0M),\r
 233   fHOutCentZEMvsZDC(ana.fHOutCentZEMvsZDC),\r
 234   fHOutMultV0M(ana.fHOutMultV0M),\r
 235   fHOutMultV0R(ana.fHOutMultV0R),\r
 236   fHOutMultFMD(ana.fHOutMultFMD),\r
 237   fHOutMultTRK(ana.fHOutMultTRK),\r
 238   fHOutMultTKL(ana.fHOutMultTKL),\r
 239   fHOutMultCL0(ana.fHOutMultCL0),\r
 240   fHOutMultCL1(ana.fHOutMultCL1),\r
 241   fHOutMultV0MvsZDC(ana.fHOutMultV0MvsZDC),\r
 242   fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC),\r
 243   fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1),\r
 244   fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK),\r
 245   fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1)\r
 246 {\r
 247   // Copy Constructor   \r
 248 }\r
 249  \r
 250 //________________________________________________________________________\r
 251 AliCentralitySelectionTask::~AliCentralitySelectionTask()\r
 252 {\r
 253   // Destructor  \r
 254   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;\r
 255   if (fTrackCuts) delete fTrackCuts;\r
 256 }  \r
 257 \r
 258 //________________________________________________________________________\r
 259 void AliCentralitySelectionTask::UserCreateOutputObjects()\r
 260 {  \r
 261   // Create the output containers\r
 262   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");\r
 263   AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);\r
 264 \r
 265   fOutputList = new TList();\r
 266   fOutputList->SetOwner();\r
 267   fHOutCentV0M     = new TH1F("fHOutCentV0M","fHOutCentV0M; Centrality V0",101,-0.5,100.5);\r
 268   fHOutCentFMD     = new TH1F("fHOutCentFMD","fHOutCentFMD; Centrality FMD",101,-0.5,100.5);\r
 269   fHOutCentTRK     = new TH1F("fHOutCentTRK","fHOutCentTRK; Centrality TPC",101,-0.5,100.5);\r
 270   fHOutCentTKL     = new TH1F("fHOutCentTKL","fHOutCentTKL; Centrality tracklets",101,-0.5,100.5);\r
 271   fHOutCentCL0     = new TH1F("fHOutCentCL0","fHOutCentCL0; Centrality SPD inner",101,-0.5,100.5);\r
 272   fHOutCentCL1     = new TH1F("fHOutCentCL1","fHOutCentCL1; Centrality SPD outer",101,-0.5,100.5);\r
 273   fHOutCentV0MvsFMD= new TH1F("fHOutCentV0MvsFMD","fHOutCentV0MvsFMD; Centrality V0 vs FMD",101,-0.5,100.5);\r
 274   fHOutCentTKLvsV0M= new TH1F("fHOutCentTKLvsV0M","fHOutCentTKLvsV0M; Centrality tracklets vs V0",101,-0.5,100.5);\r
 275   fHOutCentZEMvsZDC= new TH1F("fHOutCentZEMvsZDC","fHOutCentZEMvsZDC; Centrality ZEM vs ZDC",101,-0.5,100.5);\r
 276 \r
 277   fHOutMultV0M = new TH1F("fHOutMultV0M","fHOutMultV0M; Multiplicity V0",25000,0,25000);\r
 278   fHOutMultV0R = new TH1F("fHOutMultV0R","fHOutMultV0R; Multiplicity V0",25000,0,25000);\r
 279   fHOutMultFMD = new TH1F("fHOutMultFMD","fHOutMultFMD; Multiplicity FMD",24000,0,24000);\r
 280   fHOutMultTRK = new TH1F("fHOutMultTRK","fHOutMultTRK; Multiplicity TPC",4000,0,4000);\r
 281   fHOutMultTKL = new TH1F("fHOutMultTKL","fHOutMultTKL; Multiplicity tracklets",5000,0,5000);\r
 282   fHOutMultCL0 = new TH1F("fHOutMultCL0","fHOutMultCL0; Multiplicity SPD inner",7000,0,7000);\r
 283   fHOutMultCL1 = new TH1F("fHOutMultCL1","fHOutMultCL1; Multiplicity SPD outer",7000,0,7000);\r
 284   fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,25000,500,0,6000);\r
 285   fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,6000);\r
 286   fHOutMultV0MvsCL1 = new TH2F("fHOutMultV0MvsCL1","fHOutMultV0MvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);\r
 287   fHOutMultV0MvsTRK = new TH2F("fHOutMultV0MvsTRK","fHOutMultV0MvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);\r
 288   fHOutMultTRKvsCL1 = new TH2F("fHOutMultTRKvsCL1","fHOutMultTRKvsCL1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);\r
 289 \r
 290   fOutputList->Add(  fHOutCentV0M     );\r
 291   fOutputList->Add(  fHOutCentFMD     );\r
 292   fOutputList->Add(  fHOutCentTRK     );\r
 293   fOutputList->Add(  fHOutCentTKL     );\r
 294   fOutputList->Add(  fHOutCentCL0     );\r
 295   fOutputList->Add(  fHOutCentCL1     );\r
 296   fOutputList->Add(  fHOutCentV0MvsFMD);\r
 297   fOutputList->Add(  fHOutCentTKLvsV0M);\r
 298   fOutputList->Add(  fHOutCentZEMvsZDC);\r
 299   fOutputList->Add(  fHOutMultV0M); \r
 300   fOutputList->Add(  fHOutMultV0R); \r
 301   fOutputList->Add(  fHOutMultFMD); \r
 302   fOutputList->Add(  fHOutMultTRK); \r
 303   fOutputList->Add(  fHOutMultTKL); \r
 304   fOutputList->Add(  fHOutMultCL0); \r
 305   fOutputList->Add(  fHOutMultCL1); \r
 306   fOutputList->Add(  fHOutMultV0MvsZDC);\r
 307   fOutputList->Add(  fHOutMultZEMvsZDC);\r
 308   fOutputList->Add(  fHOutMultV0MvsCL1);\r
 309   fOutputList->Add(  fHOutMultV0MvsTRK);\r
 310   fOutputList->Add(  fHOutMultTRKvsCL1);\r
 311 \r
 312 \r
 313   fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();\r
 314 \r
 315   PostData(1, fOutputList); \r
 316 }\r
 317 \r
 318 //________________________________________________________________________\r
 319 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)\r
 320 { \r
 321   // Execute analysis for current event:\r
 322   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");\r
 323   \r
 324   Float_t  zncEnergy = 0.;          //  ZNC Energy\r
 325   Float_t  zpcEnergy = 0.;          //  ZPC Energy\r
 326   Float_t  znaEnergy = 0.;          //  ZNA Energy\r
 327   Float_t  zpaEnergy = 0.;          //  ZPA Energy\r
 328   Float_t  zem1Energy = 0.;         //  ZEM1 Energy\r
 329   Float_t  zem2Energy = 0.;         //  ZEM2 Energy\r
 330   \r
 331   Int_t    nTracks = 0;             //  no. tracks\r
 332   Int_t    nTracklets = 0;          //  no. tracklets\r
 333   Int_t    nClusters[6] = {0};      //  no. clusters on 6 ITS layers\r
 334   Int_t    nChips[2];               //  no. chips on 2 SPD layers\r
 335   Float_t  spdCorr =0;              //  corrected spd2 multiplicity\r
 336 \r
 337   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A\r
 338   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C\r
 339   Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A\r
 340   Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C\r
 341 \r
 342   Short_t v0Corr = 0;               // corrected V0 multiplicity\r
 343   Short_t v0CorrResc = 0;           // corrected and rescaled V0 multiplicity\r
 344 \r
 345   Float_t zvtx =0;                  // z-vertex SPD\r
 346  \r
 347   AliCentrality *esdCent = 0;\r
 348 \r
 349   if(fAnalysisInput.CompareTo("ESD")==0){\r
 350 \r
 351     AliVEvent* event = InputEvent();\r
 352     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);\r
 353   \r
 354     if (fRunNo<=0) {\r
 355       if (SetupRun(esd)<0)\r
 356          AliFatal("Centrality File not available for this run");\r
 357     }\r
 358 \r
 359     esdCent = esd->GetCentrality();\r
 360 \r
 361     // ***** V0 info    \r
 362     AliESDVZERO* esdV0 = esd->GetVZEROData();\r
 363     multV0A=esdV0->GetMTotV0A();\r
 364     multV0C=esdV0->GetMTotV0C();\r
 365 \r
 366     float v0CorrR;\r
 367     v0Corr = (Short_t)AliESDUtils::GetCorrV0(esd,v0CorrR);\r
 368     v0CorrResc = (Short_t)v0CorrR;\r
 369 \r
 370     // ***** Vertex Info\r
 371     const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();\r
 372     zvtx        = vtxESD->GetZ(); \r
 373 \r
 374     // ***** CB info (tracklets, clusters, chips)\r
 375     //nTracks    = event->GetNumberOfTracks();     \r
 376     nTracks    = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;\r
 377 \r
 378     const AliMultiplicity *mult = esd->GetMultiplicity();\r
 379 \r
 380     nTracklets = mult->GetNumberOfTracklets();\r
 381 \r
 382     for(Int_t ilay=0; ilay<6; ilay++){\r
 383       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);\r
 384     }\r
 385 \r
 386     for(Int_t ilay=0; ilay<2; ilay++){\r
 387       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);\r
 388     }\r
 389 \r
 390     spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);    \r
 391 \r
 392     // ***** FMD info\r
 393     AliESDFMD *fmd = esd->GetFMDData();\r
 394     Float_t totalMultA = 0;\r
 395     Float_t totalMultC = 0;\r
 396     const Float_t fFMDLowCut = 0.4;\r
 397     \r
 398     for(UShort_t det=1;det<=3;det++) {\r
 399       Int_t nRings = (det==1 ? 1 : 2);\r
 400       for (UShort_t ir = 0; ir < nRings; ir++) {          \r
 401         Char_t   ring = (ir == 0 ? 'I' : 'O');\r
 402         UShort_t nsec = (ir == 0 ? 20  : 40);\r
 403         UShort_t nstr = (ir == 0 ? 512 : 256);\r
 404         for(UShort_t sec =0; sec < nsec;  sec++)  {\r
 405           for(UShort_t strip = 0; strip < nstr; strip++) {\r
 406             \r
 407             Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);\r
 408             if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;\r
 409             \r
 410             Float_t nParticles=0;\r
 411             \r
 412             if(FMDmult > fFMDLowCut) {\r
 413               nParticles = 1.;\r
 414             }\r
 415             \r
 416             if (det<3) totalMultA = totalMultA + nParticles;\r
 417             else totalMultC = totalMultC + nParticles;\r
 418             \r
 419           }\r
 420         }\r
 421       }\r
 422     }\r
 423     multFMDA = totalMultA;\r
 424     multFMDC = totalMultC;\r
 425     \r
 426     // ***** ZDC info\r
 427     AliESDZDC *esdZDC = esd->GetESDZDC();\r
 428     zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;\r
 429     zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;\r
 430     znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;\r
 431     zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;\r
 432     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;\r
 433     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;\r
 434     \r
 435   }   \r
 436   else if(fAnalysisInput.CompareTo("AOD")==0){\r
 437     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());\r
 438     // to be implemented\r
 439     printf("  AOD analysis not yet implemented!!!\n\n");\r
 440     return;\r
 441   }\r
 442 \r
 443   // ***** Centrality Selection\r
 444   if(fHtempV0M)  fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));\r
 445   ///  else     printf("  Centrality by V0 not available!!!\n\n");\r
 446   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));\r
 447   //  else     printf("  Centrality by FMD not available!!!\n\n");\r
 448   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));\r
 449   //  else     printf("  Centrality by TRK not available!!!\n\n");\r
 450   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));\r
 451   //  else     printf("  Centrality by TKL not available!!!\n\n");\r
 452   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));\r
 453   //  else     printf("  Centrality by CL0 not available!!!\n\n");\r
 454   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));\r
 455   ///  else     printf("  Centrality by CL1 not available!!!\n\n");\r
 456   \r
 457   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));\r
 458   //  else     printf("  Centrality by V0 vs FMD not available!!!\n\n");\r
 459   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));\r
 460   //  else     printf("  Centrality by V0 vs TKL not available!!!\n\n");\r
 461   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));\r
 462   //  else     printf("  Centrality by ZEM vs ZDC not available!!!\n\n");\r
 463 \r
 464   esdCent->SetCentralityV0M(fCentV0M);\r
 465   esdCent->SetCentralityFMD(fCentFMD);\r
 466   esdCent->SetCentralityTRK(fCentTRK);\r
 467   esdCent->SetCentralityTKL(fCentTKL);\r
 468   esdCent->SetCentralityCL0(fCentCL0);\r
 469   esdCent->SetCentralityCL1(fCentCL1);\r
 470   esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);\r
 471   esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);\r
 472   esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);\r
 473 \r
 474   fHOutCentV0M->Fill(fCentV0M);\r
 475   fHOutCentFMD->Fill(fCentFMD);\r
 476   fHOutCentTRK->Fill(fCentTRK);\r
 477   fHOutCentTKL->Fill(fCentTKL);\r
 478   fHOutCentCL0->Fill(fCentCL0);\r
 479   fHOutCentCL1->Fill(fCentCL1);\r
 480   fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);\r
 481   fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);\r
 482   fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);\r
 483   fHOutMultV0M->Fill(v0Corr);\r
 484   fHOutMultV0R->Fill(multV0A+multV0C);\r
 485   fHOutMultFMD->Fill((multFMDA+multFMDC));\r
 486   fHOutMultTRK->Fill(nTracks);\r
 487   fHOutMultTKL->Fill(nTracklets);\r
 488   fHOutMultCL0->Fill(nClusters[0]);\r
 489   fHOutMultCL1->Fill(spdCorr);\r
 490   fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));\r
 491   fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));\r
 492   fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);\r
 493   fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);\r
 494   fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);\r
 495 \r
 496   PostData(1, fOutputList); \r
 497 }\r
 498 \r
 499 //________________________________________________________________________\r
 500 void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename) \r
 501 {\r
 502   //  Read centrality histograms\r
 503   TDirectory *owd = gDirectory;\r
 504   // Check if the file is present\r
 505   TString path = gSystem->ExpandPathName(fCentfilename.Data());\r
 506   if (gSystem->AccessPathName(path)) {\r
 507      AliError(Form("File %s does not exist", path.Data()));\r
 508      return;\r
 509   }\r
 510   fFile  = TFile::Open(fCentfilename);\r
 511   owd->cd();\r
 512   fHtempV0M  = (TH1F*) (fFile->Get("hmultV0_percentile"));\r
 513   fHtempFMD  = (TH1F*) (fFile->Get("hmultFMD_percentile"));\r
 514   fHtempTRK  = (TH1F*) (fFile->Get("hNtracks_percentile"));\r
 515   fHtempTKL  = (TH1F*) (fFile->Get("hNtracklets_percentile"));\r
 516   fHtempCL0  = (TH1F*) (fFile->Get("hNclusters0_percentile"));\r
 517   fHtempCL1  = (TH1F*) (fFile->Get("hNclusters1_percentile"));\r
 518   owd->cd();\r
 519 }  \r
 520 \r
 521 //________________________________________________________________________\r
 522 void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2) \r
 523 {\r
 524   //  Read centrality histograms\r
 525   TDirectory *owd = gDirectory;\r
 526   TString path = gSystem->ExpandPathName(fCentfilename2.Data());\r
 527   if (gSystem->AccessPathName(path)) {\r
 528      AliError(Form("File %s does not exist", path.Data()));\r
 529      return;\r
 530   }   \r
 531   fFile2  = TFile::Open(fCentfilename2);\r
 532   owd->cd();\r
 533   fHtempV0MvsFMD =  (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));\r
 534   fHtempTKLvsV0M  = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));\r
 535   fHtempZEMvsZDC  = (TH1F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));\r
 536   owd->cd();\r
 537 }\r
 538 \r
 539 //________________________________________________________________________\r
 540 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)\r
 541 {\r
 542   // Terminate analysis\r
 543   if (fFile && fFile->IsOpen())\r
 544     fFile->Close();  \r
 545   if (fFile2 && fFile2->IsOpen())\r
 546     fFile2->Close();  \r
 547 }\r
 548 //________________________________________________________________________\r
 549 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* esd)\r
 550 {\r
 551   // Setup files for run\r
 552 \r
 553   if (!esd)\r
 554     return -1;\r
 555 \r
 556   // check if something to be done\r
 557   if (fCurrentRun == esd->GetRunNumber())\r
 558     return 0;\r
 559   else\r
 560     fCurrentRun = esd->GetRunNumber();\r
 561   \r
 562   AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));\r
 563 \r
 564   // CHANGE HERE FOR RUN RANGES\r
 565   if ( fCurrentRun <= 137165 ) fRunNo = 137161;\r
 566   else fRunNo = 137366;\r
 567   // CHANGE HERE FOR RUN RANGES\r
 568   \r
 569   TString fileName(Form("%s/COMMON/CENTRALITY/data/AliCentralityBy1D_%d.root", AliAnalysisManager::GetOADBPath(), fRunNo));\r
 570   TString fileName2(Form("%s/COMMON/CENTRALITY/data/AliCentralityByFunction_%d.root", AliAnalysisManager::GetOADBPath(), fRunNo));\r
 571   \r
 572   AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data()));\r
 573   ReadCentralityHistos(fileName.Data());\r
 574   ReadCentralityHistos2(fileName2.Data());\r
 575   if (!fFile && !fFile2) {\r
 576      AliFatal(Form("Run %d not known to centrality selection!", fCurrentRun));       \r
 577      return -1;\r
 578   }   \r
 579   return 0;\r
 580 }\r