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     if (!esd) {\r
 354         AliError("No ESD Event");\r
 355         return;\r
 356     }\r
 357     \r
 358     if (fRunNo<=0) {\r
 359       if (SetupRun(esd)<0)\r
 360          AliFatal("Centrality File not available for this run");\r
 361     }\r
 362 \r
 363     esdCent = esd->GetCentrality();\r
 364 \r
 365     // ***** V0 info    \r
 366     AliESDVZERO* esdV0 = esd->GetVZEROData();\r
 367     multV0A=esdV0->GetMTotV0A();\r
 368     multV0C=esdV0->GetMTotV0C();\r
 369 \r
 370     float v0CorrR;\r
 371     v0Corr = (Short_t)AliESDUtils::GetCorrV0(esd,v0CorrR);\r
 372     v0CorrResc = (Short_t)v0CorrR;\r
 373 \r
 374     // ***** Vertex Info\r
 375     const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();\r
 376     zvtx        = vtxESD->GetZ(); \r
 377 \r
 378     // ***** CB info (tracklets, clusters, chips)\r
 379     //nTracks    = event->GetNumberOfTracks();     \r
 380     nTracks    = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;\r
 381 \r
 382     const AliMultiplicity *mult = esd->GetMultiplicity();\r
 383 \r
 384     nTracklets = mult->GetNumberOfTracklets();\r
 385 \r
 386     for(Int_t ilay=0; ilay<6; ilay++){\r
 387       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);\r
 388     }\r
 389 \r
 390     for(Int_t ilay=0; ilay<2; ilay++){\r
 391       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);\r
 392     }\r
 393 \r
 394     spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);    \r
 395 \r
 396     // ***** FMD info\r
 397     AliESDFMD *fmd = esd->GetFMDData();\r
 398     Float_t totalMultA = 0;\r
 399     Float_t totalMultC = 0;\r
 400     const Float_t fFMDLowCut = 0.4;\r
 401     \r
 402     for(UShort_t det=1;det<=3;det++) {\r
 403       Int_t nRings = (det==1 ? 1 : 2);\r
 404       for (UShort_t ir = 0; ir < nRings; ir++) {          \r
 405         Char_t   ring = (ir == 0 ? 'I' : 'O');\r
 406         UShort_t nsec = (ir == 0 ? 20  : 40);\r
 407         UShort_t nstr = (ir == 0 ? 512 : 256);\r
 408         for(UShort_t sec =0; sec < nsec;  sec++)  {\r
 409           for(UShort_t strip = 0; strip < nstr; strip++) {\r
 410             \r
 411             Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);\r
 412             if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;\r
 413             \r
 414             Float_t nParticles=0;\r
 415             \r
 416             if(FMDmult > fFMDLowCut) {\r
 417               nParticles = 1.;\r
 418             }\r
 419             \r
 420             if (det<3) totalMultA = totalMultA + nParticles;\r
 421             else totalMultC = totalMultC + nParticles;\r
 422             \r
 423           }\r
 424         }\r
 425       }\r
 426     }\r
 427     multFMDA = totalMultA;\r
 428     multFMDC = totalMultC;\r
 429     \r
 430     // ***** ZDC info\r
 431     AliESDZDC *esdZDC = esd->GetESDZDC();\r
 432     zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;\r
 433     zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;\r
 434     znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;\r
 435     zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;\r
 436     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;\r
 437     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;\r
 438     \r
 439   }   \r
 440   else if(fAnalysisInput.CompareTo("AOD")==0){\r
 441     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());\r
 442     // to be implemented\r
 443     printf("  AOD analysis not yet implemented!!!\n\n");\r
 444     return;\r
 445   }\r
 446 \r
 447   // ***** Centrality Selection\r
 448   if(fHtempV0M)  fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));\r
 449   ///  else     printf("  Centrality by V0 not available!!!\n\n");\r
 450   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));\r
 451   //  else     printf("  Centrality by FMD not available!!!\n\n");\r
 452   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));\r
 453   //  else     printf("  Centrality by TRK not available!!!\n\n");\r
 454   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));\r
 455   //  else     printf("  Centrality by TKL not available!!!\n\n");\r
 456   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));\r
 457   //  else     printf("  Centrality by CL0 not available!!!\n\n");\r
 458   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));\r
 459   ///  else     printf("  Centrality by CL1 not available!!!\n\n");\r
 460   \r
 461   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));\r
 462   //  else     printf("  Centrality by V0 vs FMD not available!!!\n\n");\r
 463   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));\r
 464   //  else     printf("  Centrality by V0 vs TKL not available!!!\n\n");\r
 465   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));\r
 466   //  else     printf("  Centrality by ZEM vs ZDC not available!!!\n\n");\r
 467   if (esdCent) {\r
 468       esdCent->SetCentralityV0M(fCentV0M);\r
 469       esdCent->SetCentralityFMD(fCentFMD);\r
 470       esdCent->SetCentralityTRK(fCentTRK);\r
 471       esdCent->SetCentralityTKL(fCentTKL);\r
 472       esdCent->SetCentralityCL0(fCentCL0);\r
 473       esdCent->SetCentralityCL1(fCentCL1);\r
 474       esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);\r
 475       esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);\r
 476       esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);\r
 477   }\r
 478   \r
 479   fHOutCentV0M->Fill(fCentV0M);\r
 480   fHOutCentFMD->Fill(fCentFMD);\r
 481   fHOutCentTRK->Fill(fCentTRK);\r
 482   fHOutCentTKL->Fill(fCentTKL);\r
 483   fHOutCentCL0->Fill(fCentCL0);\r
 484   fHOutCentCL1->Fill(fCentCL1);\r
 485   fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);\r
 486   fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);\r
 487   fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);\r
 488   fHOutMultV0M->Fill(v0Corr);\r
 489   fHOutMultV0R->Fill(multV0A+multV0C);\r
 490   fHOutMultFMD->Fill((multFMDA+multFMDC));\r
 491   fHOutMultTRK->Fill(nTracks);\r
 492   fHOutMultTKL->Fill(nTracklets);\r
 493   fHOutMultCL0->Fill(nClusters[0]);\r
 494   fHOutMultCL1->Fill(spdCorr);\r
 495   fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));\r
 496   fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));\r
 497   fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);\r
 498   fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);\r
 499   fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);\r
 500 \r
 501   PostData(1, fOutputList); \r
 502 }\r
 503 \r
 504 //________________________________________________________________________\r
 505 void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename) \r
 506 {\r
 507   //  Read centrality histograms\r
 508   TDirectory *owd = gDirectory;\r
 509   // Check if the file is present\r
 510   TString path = gSystem->ExpandPathName(fCentfilename.Data());\r
 511   if (gSystem->AccessPathName(path)) {\r
 512      AliError(Form("File %s does not exist", path.Data()));\r
 513      return;\r
 514   }\r
 515   fFile  = TFile::Open(fCentfilename);\r
 516   owd->cd();\r
 517   fHtempV0M  = (TH1F*) (fFile->Get("hmultV0_percentile"));\r
 518   fHtempFMD  = (TH1F*) (fFile->Get("hmultFMD_percentile"));\r
 519   fHtempTRK  = (TH1F*) (fFile->Get("hNtracks_percentile"));\r
 520   fHtempTKL  = (TH1F*) (fFile->Get("hNtracklets_percentile"));\r
 521   fHtempCL0  = (TH1F*) (fFile->Get("hNclusters0_percentile"));\r
 522   fHtempCL1  = (TH1F*) (fFile->Get("hNclusters1_percentile"));\r
 523   owd->cd();\r
 524 }  \r
 525 \r
 526 //________________________________________________________________________\r
 527 void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2) \r
 528 {\r
 529   //  Read centrality histograms\r
 530   TDirectory *owd = gDirectory;\r
 531   TString path = gSystem->ExpandPathName(fCentfilename2.Data());\r
 532   if (gSystem->AccessPathName(path)) {\r
 533      AliError(Form("File %s does not exist", path.Data()));\r
 534      return;\r
 535   }   \r
 536   fFile2  = TFile::Open(fCentfilename2);\r
 537   owd->cd();\r
 538   fHtempV0MvsFMD =  (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));\r
 539   fHtempTKLvsV0M  = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));\r
 540   fHtempZEMvsZDC  = (TH1F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));\r
 541   owd->cd();\r
 542 }\r
 543 \r
 544 //________________________________________________________________________\r
 545 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)\r
 546 {\r
 547   // Terminate analysis\r
 548   if (fFile && fFile->IsOpen())\r
 549     fFile->Close();  \r
 550   if (fFile2 && fFile2->IsOpen())\r
 551     fFile2->Close();  \r
 552 }\r
 553 //________________________________________________________________________\r
 554 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* esd)\r
 555 {\r
 556   // Setup files for run\r
 557 \r
 558   if (!esd)\r
 559     return -1;\r
 560 \r
 561   // check if something to be done\r
 562   if (fCurrentRun == esd->GetRunNumber())\r
 563     return 0;\r
 564   else\r
 565     fCurrentRun = esd->GetRunNumber();\r
 566   \r
 567   AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));\r
 568 \r
 569   // CHANGE HERE FOR RUN RANGES\r
 570   if ( fCurrentRun <= 137165 ) fRunNo = 137161;\r
 571   else fRunNo = 137366;\r
 572   // CHANGE HERE FOR RUN RANGES\r
 573   \r
 574   TString fileName(Form("%s/COMMON/CENTRALITY/data/AliCentralityBy1D_%d.root", AliAnalysisManager::GetOADBPath(), fRunNo));\r
 575   TString fileName2(Form("%s/COMMON/CENTRALITY/data/AliCentralityByFunction_%d.root", AliAnalysisManager::GetOADBPath(), fRunNo));\r
 576   \r
 577   AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data()));\r
 578   ReadCentralityHistos(fileName.Data());\r
 579   ReadCentralityHistos2(fileName2.Data());\r
 580   if (!fFile && !fFile2) {\r
 581      AliFatal(Form("Run %d not known to centrality selection!", fCurrentRun));       \r
 582      return -1;\r
 583   }   \r
 584   return 0;\r
 585 }\r