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