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