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 <TF1.h>
  29 #include <TProfile.h>
  30 #include <TFile.h>
  31 #include <TObjString.h>
  32 #include <TString.h>
  33 #include <TCanvas.h>
  34 #include <TROOT.h>
  35 #include <TDirectory.h>
  36 #include <TSystem.h>
  37 #include <iostream>
  38
  39 #include "AliAnalysisManager.h"
  40 #include "AliHeader.h"
  41 #include "AliVEvent.h"
  42 #include "AliESD.h"
  43 #include "AliESDEvent.h"
  44 #include "AliESDHeader.h"
  45 #include "AliESDInputHandler.h"
  46 #include "AliESDZDC.h"
  47 #include "AliESDFMD.h"
  48 #include "AliESDTZERO.h"
  49 #include "AliESDVZERO.h"
  50 #include "AliESDtrackCuts.h"
  51 #include "AliESDVertex.h"
  52 #include "AliCentrality.h"
  53 #include "AliOADBCentrality.h"
  54 #include "AliOADBContainer.h"
  55 #include "AliMultiplicity.h"
  56 #include "AliAODHandler.h"
  57 #include "AliAODHeader.h"
  58 #include "AliAODEvent.h"
  59 #include "AliAODVertex.h"
  60 #include "AliAODVZERO.h"
  61 #include "AliAODTracklets.h"
  62 #include "AliAODMCHeader.h"
  63 #include "AliMCEventHandler.h"
  64 #include "AliMCEvent.h"
  65 #include "AliAODMCParticle.h"
  66 #include "AliMCParticle.h"
  67 #include "AliStack.h"
  68 #include "AliAnalysisTaskSE.h"
  69 #include "AliGenEventHeader.h"
  70 #include "AliGenHijingEventHeader.h"
  71 #include "AliPhysicsSelectionTask.h"
  72 #include "AliPhysicsSelection.h"
  73 #include "AliBackgroundSelection.h"
  74 #include "AliESDUtils.h"
  75
  76 ClassImp(AliCentralitySelectionTask)
  77
  78
  79 //________________________________________________________________________
  80 AliCentralitySelectionTask::AliCentralitySelectionTask():
  81 AliAnalysisTaskSE(),
  82   fAnalysisInput("ESD"),
  83   fIsMCInput(kFALSE),
  84   fCurrentRun(-1),
  85   fUseScaling(0),
  86   fUseCleaning(0),
  87   fFillHistos(0),
  88   fV0MScaleFactor(0),
  89   fSPDScaleFactor(0),
  90   fTPCScaleFactor(0),
  91   fV0MScaleFactorMC(0),
  92   fV0MSPDOutlierPar0(0),
  93   fV0MSPDOutlierPar1(0),
  94   fV0MTPCOutlierPar0(0),
  95   fV0MTPCOutlierPar1(0),
  96   fV0MSPDSigmaOutlierPar0(0),
  97   fV0MSPDSigmaOutlierPar1(0),
  98   fV0MSPDSigmaOutlierPar2(0),
  99   fV0MTPCSigmaOutlierPar0(0),
 100   fV0MTPCSigmaOutlierPar1(0),
 101   fV0MTPCSigmaOutlierPar2(0),
 102   fV0MZDCOutlierPar0(0),
 103   fV0MZDCOutlierPar1(0),
 104   fV0MZDCEcalOutlierPar0(0),
 105   fV0MZDCEcalOutlierPar1(0),
 106   fTrackCuts(0),
 107   fZVCut(10),
 108   fOutliersCut(5),
 109   fQuality(999),
 110   fIsSelected(0),
 111   fMSL(0),
 112   fMSH(0),
 113   fMUL(0),
 114   fMLL(0),
 115   fEJE(0),
 116   fEGA(0),
 117   fPHS(0),
 118   fCVHN(0),
 119   fCVLN(0),
 120   fCVHNbit(0),
 121   fCVLNbit(0),
 122   fCCENT(0),
 123   fCSEMI(0),
 124   fCCENTbit(0),
 125   fCSEMIbit(0),
 126   fCentV0M(0),
 127   fCentFMD(0),
 128   fCentTRK(0),
 129   fCentTKL(0),
 130   fCentCL0(0),
 131   fCentCL1(0),
 132   fCentV0MvsFMD(0),
 133   fCentTKLvsV0M(0),
 134   fCentZEMvsZDC(0),
 135   fHtempV0M(0),
 136   fHtempFMD(0),
 137   fHtempTRK(0),
 138   fHtempTKL(0),
 139   fHtempCL0(0),
 140   fHtempCL1(0),
 141   fHtempV0MvsFMD(0),
 142   fHtempTKLvsV0M(0),
 143   fHtempZEMvsZDC(0),
 144   fOutputList(0),
 145   fHOutCentV0M     (0),
 146   fHOutCentV0MCVHN(0),
 147   fHOutCentV0MCVLN(0),
 148   fHOutCentV0MCVHNinMB(0),
 149   fHOutCentV0MCVLNinMB(0),
 150   fHOutCentV0MCCENT(0),
 151   fHOutCentV0MCSEMI(0),
 152   fHOutCentV0MCCENTinMB(0),
 153   fHOutCentV0MCSEMIinMB(0),
 154   fHOutCentV0MMSL(0),
 155   fHOutCentV0MMSH(0),
 156   fHOutCentV0MMUL(0),
 157   fHOutCentV0MMLL(0),
 158   fHOutCentV0MEJE(0),
 159   fHOutCentV0MEGA(0),
 160   fHOutCentV0MPHS(0),
 161   fHOutCentV0MMSLinMB(0),
 162   fHOutCentV0MMSHinMB(0),
 163   fHOutCentV0MMULinMB(0),
 164   fHOutCentV0MMLLinMB(0),
 165   fHOutCentV0MEJEinMB(0),
 166   fHOutCentV0MEGAinMB(0),
 167   fHOutCentV0MPHSinMB(0),
 168   fHOutCentFMD     (0),
 169   fHOutCentTRK     (0),
 170   fHOutCentTKL     (0),
 171   fHOutCentCL0     (0),
 172   fHOutCentCL1     (0),
 173   fHOutCentV0MvsFMD(0),
 174   fHOutCentTKLvsV0M(0),
 175   fHOutCentZEMvsZDC(0),
 176   fHOutCentV0MvsCentCL1(0),
 177   fHOutCentV0MvsCentTRK(0),
 178   fHOutCentTRKvsCentCL1(0),
 179   fHOutCentV0MvsCentZDC(0),
 180   fHOutMultV0M(0),
 181   fHOutMultV0O(0),
 182   fHOutMultFMD(0),
 183   fHOutMultTRK(0),
 184   fHOutMultTKL(0),
 185   fHOutMultCL0(0),
 186   fHOutMultCL1(0),
 187   fHOutMultV0MvsZDN(0),
 188   fHOutMultZEMvsZDN(0),
 189   fHOutMultV0MvsZDC(0),
 190   fHOutMultZEMvsZDC(0),
 191   fHOutMultZEMvsZDCw(0),
 192   fHOutMultV0MvsCL1(0),
 193   fHOutMultV0MvsTRK(0),
 194   fHOutMultTRKvsCL1(0),
 195   fHOutMultV0MvsV0O(0),
 196   fHOutMultV0OvsCL1(0),
 197   fHOutMultV0OvsTRK(0),
 198   fHOutCentV0Mqual1(0),
 199   fHOutCentTRKqual1(0),
 200   fHOutCentCL1qual1(0),
 201   fHOutMultV0MvsCL1qual1(0),
 202   fHOutMultV0MvsTRKqual1(0),
 203   fHOutMultTRKvsCL1qual1(0),
 204   fHOutCentV0Mqual2(0),
 205   fHOutCentTRKqual2(0),
 206   fHOutCentCL1qual2(0),
 207   fHOutMultV0MvsCL1qual2(0),
 208   fHOutMultV0MvsTRKqual2(0),
 209   fHOutMultTRKvsCL1qual2(0),
 210   fHOutQuality(0),
 211   fHOutVertex(0),
 212   fHOutVertexT0(0)
 213 {
 214   // Default constructor
 215   AliInfo("Centrality Selection enabled.");
 216
 217   fUseScaling=kTRUE;
 218   fUseCleaning=kTRUE;
 219   fFillHistos=kFALSE;
 220   fBranchNames="ESD:AliESDRun.,AliESDHeader.,AliESDZDC.,AliESDFMD.,AliESDVZERO.,AliESDTZERO."
 221     ",SPDVertex.,TPCVertex.,PrimaryVertex.,AliMultiplicity.,Tracks ";
 222 }
 223
 224 //________________________________________________________________________
 225 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
 226   AliAnalysisTaskSE(name),
 227   fAnalysisInput("ESD"),
 228   fIsMCInput(kFALSE),
 229   fCurrentRun(-1),
 230   fUseScaling(0),
 231   fUseCleaning(0),
 232   fFillHistos(0),
 233   fV0MScaleFactor(0),
 234   fSPDScaleFactor(0),
 235   fTPCScaleFactor(0),
 236   fV0MScaleFactorMC(0),
 237   fV0MSPDOutlierPar0(0),
 238   fV0MSPDOutlierPar1(0),
 239   fV0MTPCOutlierPar0(0),
 240   fV0MTPCOutlierPar1(0),
 241   fV0MSPDSigmaOutlierPar0(0),
 242   fV0MSPDSigmaOutlierPar1(0),
 243   fV0MSPDSigmaOutlierPar2(0),
 244   fV0MTPCSigmaOutlierPar0(0),
 245   fV0MTPCSigmaOutlierPar1(0),
 246   fV0MTPCSigmaOutlierPar2(0),
 247   fV0MZDCOutlierPar0(0),
 248   fV0MZDCOutlierPar1(0),
 249   fV0MZDCEcalOutlierPar0(0),
 250   fV0MZDCEcalOutlierPar1(0),
 251   fTrackCuts(0),
 252   fZVCut(10),
 253   fOutliersCut(5),
 254   fQuality(999),
 255   fIsSelected(0),
 256   fMSL(0),
 257   fMSH(0),
 258   fMUL(0),
 259   fMLL(0),
 260   fEJE(0),
 261   fEGA(0),
 262   fPHS(0),
 263   fCVHN(0),
 264   fCVLN(0),
 265   fCVHNbit(0),
 266   fCVLNbit(0),
 267   fCCENT(0),
 268   fCSEMI(0),
 269   fCCENTbit(0),
 270   fCSEMIbit(0),
 271   fCentV0M(0),
 272   fCentFMD(0),
 273   fCentTRK(0),
 274   fCentTKL(0),
 275   fCentCL0(0),
 276   fCentCL1(0),
 277   fCentV0MvsFMD(0),
 278   fCentTKLvsV0M(0),
 279   fCentZEMvsZDC(0),
 280   fHtempV0M(0),
 281   fHtempFMD(0),
 282   fHtempTRK(0),
 283   fHtempTKL(0),
 284   fHtempCL0(0),
 285   fHtempCL1(0),
 286   fHtempV0MvsFMD(0),
 287   fHtempTKLvsV0M(0),
 288   fHtempZEMvsZDC(0),
 289   fOutputList(0),
 290   fHOutCentV0M     (0),
 291   fHOutCentV0MCVHN(0),
 292   fHOutCentV0MCVLN(0),
 293   fHOutCentV0MCVHNinMB(0),
 294   fHOutCentV0MCVLNinMB(0),
 295   fHOutCentV0MCCENT(0),
 296   fHOutCentV0MCSEMI(0),
 297   fHOutCentV0MCCENTinMB(0),
 298   fHOutCentV0MCSEMIinMB(0),
 299   fHOutCentV0MMSL(0),
 300   fHOutCentV0MMSH(0),
 301   fHOutCentV0MMUL(0),
 302   fHOutCentV0MMLL(0),
 303   fHOutCentV0MEJE(0),
 304   fHOutCentV0MEGA(0),
 305   fHOutCentV0MPHS(0),
 306   fHOutCentV0MMSLinMB(0),
 307   fHOutCentV0MMSHinMB(0),
 308   fHOutCentV0MMULinMB(0),
 309   fHOutCentV0MMLLinMB(0),
 310   fHOutCentV0MEJEinMB(0),
 311   fHOutCentV0MEGAinMB(0),
 312   fHOutCentV0MPHSinMB(0),
 313   fHOutCentFMD     (0),
 314   fHOutCentTRK     (0),
 315   fHOutCentTKL     (0),
 316   fHOutCentCL0     (0),
 317   fHOutCentCL1     (0),
 318   fHOutCentV0MvsFMD(0),
 319   fHOutCentTKLvsV0M(0),
 320   fHOutCentZEMvsZDC(0),
 321   fHOutCentV0MvsCentCL1(0),
 322   fHOutCentV0MvsCentTRK(0),
 323   fHOutCentTRKvsCentCL1(0),
 324   fHOutCentV0MvsCentZDC(0),
 325   fHOutMultV0M(0),
 326   fHOutMultV0O(0),
 327   fHOutMultFMD(0),
 328   fHOutMultTRK(0),
 329   fHOutMultTKL(0),
 330   fHOutMultCL0(0),
 331   fHOutMultCL1(0),
 332   fHOutMultV0MvsZDN(0),
 333   fHOutMultZEMvsZDN(0),
 334   fHOutMultV0MvsZDC(0),
 335   fHOutMultZEMvsZDC(0),
 336   fHOutMultZEMvsZDCw(0),
 337   fHOutMultV0MvsCL1(0),
 338   fHOutMultV0MvsTRK(0),
 339   fHOutMultTRKvsCL1(0),
 340   fHOutMultV0MvsV0O(0),
 341   fHOutMultV0OvsCL1(0),
 342   fHOutMultV0OvsTRK(0),
 343   fHOutCentV0Mqual1(0),
 344   fHOutCentTRKqual1(0),
 345   fHOutCentCL1qual1(0),
 346   fHOutMultV0MvsCL1qual1(0),
 347   fHOutMultV0MvsTRKqual1(0),
 348   fHOutMultTRKvsCL1qual1(0),
 349   fHOutCentV0Mqual2(0),
 350   fHOutCentTRKqual2(0),
 351   fHOutCentCL1qual2(0),
 352   fHOutMultV0MvsCL1qual2(0),
 353   fHOutMultV0MvsTRKqual2(0),
 354   fHOutMultTRKvsCL1qual2(0),
 355   fHOutQuality(0),
 356   fHOutVertex(0),
 357   fHOutVertexT0(0)
 358 {
 359   // Default constructor
 360   AliInfo("Centrality Selection enabled.");
 361   //DefineOutput(1, TList::Class());
 362   fUseScaling=kTRUE;
 363   fUseCleaning=kTRUE;
 364   fFillHistos=kFALSE;
 365   fBranchNames="ESD:AliESDRun.,AliESDHeader.,AliESDZDC.,AliESDFMD.,AliESDVZERO.,AliESDTZERO."
 366     ",SPDVertex.,TPCVertex.,PrimaryVertex.,AliMultiplicity.,Tracks ";
 367 }
 368
 369 //________________________________________________________________________
 370 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
 371 {
 372   // Assignment operator
 373   if (this!=&c) {
 374     AliAnalysisTaskSE::operator=(c);
 375   }
 376   return *this;
 377 }
 378
 379 //________________________________________________________________________
 380 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
 381   AliAnalysisTaskSE(ana),
 382   fAnalysisInput(ana.fAnalysisInput),
 383   fIsMCInput(ana.fIsMCInput),
 384   fCurrentRun(ana.fCurrentRun),
 385   fUseScaling(ana.fUseScaling),
 386   fUseCleaning(ana.fUseCleaning),
 387   fFillHistos(ana.fFillHistos),
 388   fV0MScaleFactor(ana.fV0MScaleFactor),
 389   fSPDScaleFactor(ana.fSPDScaleFactor),
 390   fTPCScaleFactor(ana.fTPCScaleFactor),
 391   fV0MScaleFactorMC(ana.fV0MScaleFactorMC),
 392   fV0MSPDOutlierPar0(ana.fV0MSPDOutlierPar0),
 393   fV0MSPDOutlierPar1(ana.fV0MSPDOutlierPar1),
 394   fV0MTPCOutlierPar0(ana.fV0MTPCOutlierPar0),
 395   fV0MTPCOutlierPar1(ana.fV0MTPCOutlierPar1),
 396   fV0MSPDSigmaOutlierPar0(ana.fV0MSPDSigmaOutlierPar0),
 397   fV0MSPDSigmaOutlierPar1(ana.fV0MSPDSigmaOutlierPar1),
 398   fV0MSPDSigmaOutlierPar2(ana.fV0MSPDSigmaOutlierPar2),
 399   fV0MTPCSigmaOutlierPar0(ana.fV0MTPCSigmaOutlierPar0),
 400   fV0MTPCSigmaOutlierPar1(ana.fV0MTPCSigmaOutlierPar1),
 401   fV0MTPCSigmaOutlierPar2(ana.fV0MTPCSigmaOutlierPar2),
 402   fV0MZDCOutlierPar0(ana.fV0MZDCOutlierPar0),
 403   fV0MZDCOutlierPar1(ana.fV0MZDCOutlierPar1),
 404   fV0MZDCEcalOutlierPar0(ana.fV0MZDCEcalOutlierPar0),
 405   fV0MZDCEcalOutlierPar1(ana.fV0MZDCEcalOutlierPar1),
 406   fTrackCuts(ana.fTrackCuts),
 407   fZVCut(ana.fZVCut),
 408   fOutliersCut(ana.fOutliersCut),
 409   fQuality(ana.fQuality),
 410   fIsSelected(ana.fIsSelected),
 411   fMSL(ana.fMSL),
 412   fMSH(ana.fMSH),
 413   fMUL(ana.fMUL),
 414   fMLL(ana.fMLL),
 415   fEJE(ana.fEJE),
 416   fEGA(ana.fEGA),
 417   fPHS(ana.fPHS),
 418   fCVHN(ana.fCVHN),
 419   fCVLN(ana.fCVLN),
 420   fCVHNbit(ana.fCVHNbit),
 421   fCVLNbit(ana.fCVLNbit),
 422   fCCENT(ana.fCCENT),
 423   fCSEMI(ana.fCSEMI),
 424   fCCENTbit(ana.fCCENTbit),
 425   fCSEMIbit(ana.fCSEMIbit),
 426   fCentV0M(ana.fCentV0M),
 427   fCentFMD(ana.fCentFMD),
 428   fCentTRK(ana.fCentTRK),
 429   fCentTKL(ana.fCentTKL),
 430   fCentCL0(ana.fCentCL0),
 431   fCentCL1(ana.fCentCL1),
 432   fCentV0MvsFMD(ana.fCentV0MvsFMD),
 433   fCentTKLvsV0M(ana.fCentTKLvsV0M),
 434   fCentZEMvsZDC(ana.fCentZEMvsZDC),
 435   fHtempV0M(ana.fHtempV0M),
 436   fHtempFMD(ana.fHtempFMD),
 437   fHtempTRK(ana.fHtempTRK),
 438   fHtempTKL(ana.fHtempTKL),
 439   fHtempCL0(ana.fHtempCL0),
 440   fHtempCL1(ana.fHtempCL1),
 441   fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
 442   fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
 443   fHtempZEMvsZDC(ana.fHtempZEMvsZDC),
 444   fOutputList(ana.fOutputList),
 445   fHOutCentV0M     (ana.fHOutCentV0M     ),
 446   fHOutCentV0MCVHN(ana.fHOutCentV0MCVHN),
 447   fHOutCentV0MCVLN(ana.fHOutCentV0MCVLN),
 448   fHOutCentV0MCVHNinMB(ana.fHOutCentV0MCVHNinMB),
 449   fHOutCentV0MCVLNinMB(ana.fHOutCentV0MCVLNinMB),
 450   fHOutCentV0MCCENT(ana.fHOutCentV0MCCENT),
 451   fHOutCentV0MCSEMI(ana.fHOutCentV0MCSEMI),
 452   fHOutCentV0MCCENTinMB(ana.fHOutCentV0MCCENTinMB),
 453   fHOutCentV0MCSEMIinMB(ana.fHOutCentV0MCSEMIinMB),
 454   fHOutCentV0MMSL(ana.fHOutCentV0MMSL),
 455   fHOutCentV0MMSH(ana.fHOutCentV0MMSH),
 456   fHOutCentV0MMUL(ana.fHOutCentV0MMUL),
 457   fHOutCentV0MMLL(ana.fHOutCentV0MMLL),
 458   fHOutCentV0MEJE(ana.fHOutCentV0MEJE),
 459   fHOutCentV0MEGA(ana.fHOutCentV0MEGA),
 460   fHOutCentV0MPHS(ana.fHOutCentV0MPHS),
 461   fHOutCentV0MMSLinMB(ana.fHOutCentV0MMSLinMB),
 462   fHOutCentV0MMSHinMB(ana.fHOutCentV0MMSHinMB),
 463   fHOutCentV0MMULinMB(ana.fHOutCentV0MMULinMB),
 464   fHOutCentV0MMLLinMB(ana.fHOutCentV0MMLLinMB),
 465   fHOutCentV0MEJEinMB(ana.fHOutCentV0MEJEinMB),
 466   fHOutCentV0MEGAinMB(ana.fHOutCentV0MEGAinMB),
 467   fHOutCentV0MPHSinMB(ana.fHOutCentV0MPHSinMB),
 468   fHOutCentFMD     (ana.fHOutCentFMD     ),
 469   fHOutCentTRK     (ana.fHOutCentTRK     ),
 470   fHOutCentTKL     (ana.fHOutCentTKL     ),
 471   fHOutCentCL0     (ana.fHOutCentCL0     ),
 472   fHOutCentCL1     (ana.fHOutCentCL1     ),
 473   fHOutCentV0MvsFMD(ana.fHOutCentV0MvsFMD),
 474   fHOutCentTKLvsV0M(ana.fHOutCentTKLvsV0M),
 475   fHOutCentZEMvsZDC(ana.fHOutCentZEMvsZDC),
 476   fHOutCentV0MvsCentCL1(ana.fHOutCentV0MvsCentCL1),
 477   fHOutCentV0MvsCentTRK(ana.fHOutCentV0MvsCentTRK),
 478   fHOutCentTRKvsCentCL1(ana.fHOutCentTRKvsCentCL1),
 479   fHOutCentV0MvsCentZDC(ana.fHOutCentV0MvsCentZDC),
 480   fHOutMultV0M(ana.fHOutMultV0M),
 481   fHOutMultV0O(ana.fHOutMultV0O),
 482   fHOutMultFMD(ana.fHOutMultFMD),
 483   fHOutMultTRK(ana.fHOutMultTRK),
 484   fHOutMultTKL(ana.fHOutMultTKL),
 485   fHOutMultCL0(ana.fHOutMultCL0),
 486   fHOutMultCL1(ana.fHOutMultCL1),
 487   fHOutMultV0MvsZDN(ana.fHOutMultV0MvsZDN),
 488   fHOutMultZEMvsZDN(ana.fHOutMultZEMvsZDN),
 489   fHOutMultV0MvsZDC(ana.fHOutMultV0MvsZDC),
 490   fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC),
 491   fHOutMultZEMvsZDCw(ana.fHOutMultZEMvsZDCw),
 492   fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1),
 493   fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK),
 494   fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1),
 495   fHOutMultV0MvsV0O(ana.fHOutMultV0MvsV0O),
 496   fHOutMultV0OvsCL1(ana.fHOutMultV0OvsCL1),
 497   fHOutMultV0OvsTRK(ana.fHOutMultV0OvsTRK),
 498   fHOutCentV0Mqual1(ana.fHOutCentV0Mqual1),
 499   fHOutCentTRKqual1(ana.fHOutCentTRKqual1),
 500   fHOutCentCL1qual1(ana.fHOutCentCL1qual1),
 501   fHOutMultV0MvsCL1qual1(ana.fHOutMultV0MvsCL1qual1),
 502   fHOutMultV0MvsTRKqual1(ana.fHOutMultV0MvsTRKqual1),
 503   fHOutMultTRKvsCL1qual1(ana.fHOutMultTRKvsCL1qual1),
 504   fHOutCentV0Mqual2(ana.fHOutCentV0Mqual2),
 505   fHOutCentTRKqual2(ana.fHOutCentTRKqual2),
 506   fHOutCentCL1qual2(ana.fHOutCentCL1qual2),
 507   fHOutMultV0MvsCL1qual2(ana.fHOutMultV0MvsCL1qual2),
 508   fHOutMultV0MvsTRKqual2(ana.fHOutMultV0MvsTRKqual2),
 509   fHOutMultTRKvsCL1qual2(ana.fHOutMultTRKvsCL1qual2),
 510   fHOutQuality(ana.fHOutQuality),
 511   fHOutVertex(ana.fHOutVertex),
 512   fHOutVertexT0(ana.fHOutVertexT0)
 513 {
 514   // Copy Constructor
 515
 516 }
 517
 518 //________________________________________________________________________
 519 AliCentralitySelectionTask::~AliCentralitySelectionTask()
 520 {
 521   // Destructor
 522   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
 523   if (fTrackCuts) delete fTrackCuts;
 524 }
 525
 526 //________________________________________________________________________
 527 void AliCentralitySelectionTask::UserCreateOutputObjects()
 528 {
 529   // Create the output containers
 530   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
 531   AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
 532
 533   if (fFillHistos) {
 534     fOutputList = new TList();
 535     fOutputList->SetOwner();
 536     fHOutCentV0M     = new TH1F("fHOutCentV0M","fHOutCentV0M; Centrality V0",505,0,101);
 537     fHOutCentV0MCVHN= new TH1F("fHOutCentV0M_CVHN","fHOutCentV0M_CVHN; Centrality V0",505,0,101);
 538     fHOutCentV0MCVLN= new TH1F("fHOutCentV0M_CVLN","fHOutCentV0M_CVLN; Centrality V0",505,0,101);
 539     fHOutCentV0MCVHNinMB= new TH1F("fHOutCentV0M_CVHNinMB","fHOutCentV0M_CVHN; Centrality V0",505,0,101);
 540     fHOutCentV0MCVLNinMB= new TH1F("fHOutCentV0M_CVLNinMB","fHOutCentV0M_CVLN; Centrality V0",505,0,101);
 541     fHOutCentV0MCCENT= new TH1F("fHOutCentV0M_CCENT","fHOutCentV0M_CCENT; Centrality V0",505,0,101);
 542     fHOutCentV0MCSEMI= new TH1F("fHOutCentV0M_CSEMI","fHOutCentV0M_CSEMI; Centrality V0",505,0,101);
 543     fHOutCentV0MCCENTinMB= new TH1F("fHOutCentV0M_CCENTinMB","fHOutCentV0M_CCENT; Centrality V0",505,0,101);
 544     fHOutCentV0MCSEMIinMB= new TH1F("fHOutCentV0M_CSEMIinMB","fHOutCentV0M_CSEMI; Centrality V0",505,0,101);
 545     fHOutCentV0MMSL= new TH1F("fHOutCentV0M_MSL","fHOutCentV0M_MSL; Centrality V0",505,0,101);
 546     fHOutCentV0MMSH= new TH1F("fHOutCentV0M_MSH","fHOutCentV0M_MSH; Centrality V0",505,0,101);
 547     fHOutCentV0MMUL= new TH1F("fHOutCentV0M_MUL","fHOutCentV0M_MUL; Centrality V0",505,0,101);
 548     fHOutCentV0MMLL= new TH1F("fHOutCentV0M_MLL","fHOutCentV0M_MLL; Centrality V0",505,0,101);
 549     fHOutCentV0MEJE= new TH1F("fHOutCentV0M_EJE","fHOutCentV0M_EJE; Centrality V0",505,0,101);
 550     fHOutCentV0MEGA= new TH1F("fHOutCentV0M_EGA","fHOutCentV0M_EGA; Centrality V0",505,0,101);
 551     fHOutCentV0MPHS= new TH1F("fHOutCentV0M_PHS","fHOutCentV0M_PHS; Centrality V0",505,0,101);
 552     fHOutCentV0MMSLinMB= new TH1F("fHOutCentV0M_MSLinMB","fHOutCentV0M_MSLinMB; Centrality V0",505,0,101);
 553     fHOutCentV0MMSHinMB= new TH1F("fHOutCentV0M_MSHinMB","fHOutCentV0M_MSHinMB; Centrality V0",505,0,101);
 554     fHOutCentV0MMULinMB= new TH1F("fHOutCentV0M_MULinMB","fHOutCentV0M_MULinMB; Centrality V0",505,0,101);
 555     fHOutCentV0MMLLinMB= new TH1F("fHOutCentV0M_MLLinMB","fHOutCentV0M_MLLinMB; Centrality V0",505,0,101);
 556     fHOutCentV0MEJEinMB= new TH1F("fHOutCentV0M_EJEinMB","fHOutCentV0M_EJEinMB; Centrality V0",505,0,101);
 557     fHOutCentV0MEGAinMB= new TH1F("fHOutCentV0M_EGAinMB","fHOutCentV0M_EGAinMB; Centrality V0",505,0,101);
 558     fHOutCentV0MPHSinMB= new TH1F("fHOutCentV0M_PHSinMB","fHOutCentV0M_PHSinMB; Centrality V0",505,0,101);
 559     fHOutCentFMD     = new TH1F("fHOutCentFMD","fHOutCentFMD; Centrality FMD",505,0,101);
 560     fHOutCentTRK     = new TH1F("fHOutCentTRK","fHOutCentTRK; Centrality TPC",505,0,101);
 561     fHOutCentTKL     = new TH1F("fHOutCentTKL","fHOutCentTKL; Centrality tracklets",505,0,101);
 562     fHOutCentCL0     = new TH1F("fHOutCentCL0","fHOutCentCL0; Centrality SPD inner",505,0,101);
 563     fHOutCentCL1     = new TH1F("fHOutCentCL1","fHOutCentCL1; Centrality SPD outer",505,0,101);
 564     fHOutCentV0MvsFMD= new TH1F("fHOutCentV0MvsFMD","fHOutCentV0MvsFMD; Centrality V0 vs FMD",505,0,101);
 565     fHOutCentTKLvsV0M= new TH1F("fHOutCentTKLvsV0M","fHOutCentTKLvsV0M; Centrality tracklets vs V0",505,0,101);
 566     fHOutCentZEMvsZDC= new TH1F("fHOutCentZEMvsZDC","fHOutCentZEMvsZDC; Centrality ZEM vs ZDC",505,0,101);
 567     fHOutCentV0MvsCentCL1= new TH2F("fHOutCentV0MvsCentCL1","fHOutCentV0MvsCentCL1; Cent V0; Cent SPD",505,0,101,505,0,101);
 568     fHOutCentV0MvsCentTRK= new TH2F("fHOutCentV0MvsCentTRK","fHOutCentV0MvsCentTRK; Cent V0; Cent TPC",505,0,101,505,0,101);
 569     fHOutCentTRKvsCentCL1= new TH2F("fHOutCentTRKvsCentCL1","fHOutCentTRKvsCentCL1; Cent TPC; Cent SPD",505,0,101,505,0,101);
 570     fHOutCentV0MvsCentZDC= new TH2F("fHOutCentV0MvsCentZDC","fHOutCentV0MvsCentZDC; Cent V0; Cent ZDC",505,0,101,505,0,101);
 571     fHOutMultV0M = new TH1F("fHOutMultV0M","fHOutMultV0M; Multiplicity V0",25000,0,30000);
 572     fHOutMultV0O = new TH1F("fHOutMultV0O","fHOutMultV0O; Multiplicity V0",40000,0,40000);
 573     fHOutMultFMD = new TH1F("fHOutMultFMD","fHOutMultFMD; Multiplicity FMD",24000,0,24000);
 574     fHOutMultTRK = new TH1F("fHOutMultTRK","fHOutMultTRK; Multiplicity TPC",4000,0,4000);
 575     fHOutMultTKL = new TH1F("fHOutMultTKL","fHOutMultTKL; Multiplicity tracklets",5000,0,5000);
 576     fHOutMultCL0 = new TH1F("fHOutMultCL0","fHOutMultCL0; Multiplicity SPD inner",7000,0,7000);
 577     fHOutMultCL1 = new TH1F("fHOutMultCL1","fHOutMultCL1; Multiplicity SPD outer",7000,0,7000);
 578     fHOutMultV0MvsZDN = new TH2F("fHOutMultV0MvsZDN","fHOutMultV0MvsZDN; Multiplicity V0; Energy ZDC-N",500,0,30000,500,0,180000);
 579     fHOutMultZEMvsZDN = new TH2F("fHOutMultZEMvsZDN","fHOutMultZEMvsZDN; Energy ZEM; Energy ZDC-N",500,0,2500,500,0,180000);
 580     fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,30000,500,0,200000);
 581     fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,200000);
 582     fHOutMultZEMvsZDCw = new TH2F("fHOutMultZEMvsZDCw","fHOutMultZEMvsZDCw; Energy ZEM; Energy ZDC (weigthed with V0 percentile)",500,0,2500,500,0,200000);
 583     fHOutMultV0MvsCL1 = new TH2F("fHOutMultV0MvsCL1","fHOutMultV0MvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,30000,700,0,7000);
 584     fHOutMultV0MvsTRK = new TH2F("fHOutMultV0MvsTRK","fHOutMultV0MvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,30000,400,0,4000);
 585     fHOutMultTRKvsCL1 = new TH2F("fHOutMultTRKvsCL1","fHOutMultTRKvsCL1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
 586     fHOutMultV0MvsV0O = new TH2F("fHOutMultV0MvsV0O","fHOutMultV0MvsV0O; Multiplicity V0; Multiplicity V0 Online",500,0,30000,500,0,40000);
 587     fHOutMultV0OvsCL1 = new TH2F("fHOutMultV0OvsCL1","fHOutMultV0OvsCL1; Multiplicity V0; Multiplicity SPD outer",500,0,40000,700,0,7000);
 588     fHOutMultV0OvsTRK = new TH2F("fHOutMultV0OvsTRK","fHOutMultV0OvsTRK; Multiplicity V0; Multiplicity TPC",500,0,40000,400,0,4000);
 589     fHOutMultV0MvsV0O = new TH2F("fHOutMultV0MvsV0O","fHOutMultV0MvsV0O; Multiplicity V0; Multiplicity V0 Online",500,0,30000,500,0,30000);
 590     fHOutMultV0OvsCL1 = new TH2F("fHOutMultV0OvsCL1","fHOutMultV0OvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,30000,700,0,7000);
 591     fHOutMultV0OvsTRK = new TH2F("fHOutMultV0OvsTRK","fHOutMultV0OvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,30000,400,0,4000);
 592
 593     fHOutCentV0Mqual1 = new TH1F("fHOutCentV0M_qual1","fHOutCentV0M_qual1; Centrality V0",505,0,101);
 594     fHOutCentTRKqual1 = new TH1F("fHOutCentTRK_qual1","fHOutCentTRK_qual1; Centrality TPC",505,0,101);
 595     fHOutCentCL1qual1 = new TH1F("fHOutCentCL1_qual1","fHOutCentCL1_qual1; Centrality SPD outer",505,0,101);
 596     fHOutMultV0MvsCL1qual1 = new TH2F("fHOutMultV0MvsCL1_qual1","fHOutMultV0MvsCL1_qual1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
 597     fHOutMultV0MvsTRKqual1 = new TH2F("fHOutMultV0MvsTRK_qual1","fHOutMultV0MvsTRK_qual1; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
 598     fHOutMultTRKvsCL1qual1 = new TH2F("fHOutMultTRKvsCL1_qual1","fHOutMultTRKvsCL1_qual1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
 599
 600     fHOutCentV0Mqual2 = new TH1F("fHOutCentV0M_qual2","fHOutCentV0M_qual2; Centrality V0",505,0,101);
 601     fHOutCentTRKqual2 = new TH1F("fHOutCentTRK_qual2","fHOutCentTRK_qual2; Centrality TPC",505,0,101);
 602     fHOutCentCL1qual2 = new TH1F("fHOutCentCL1_qual2","fHOutCentCL1_qual2; Centrality SPD outer",505,0,101);
 603     fHOutMultV0MvsCL1qual2 = new TH2F("fHOutMultV0MvsCL1_qual2","fHOutMultV0MvsCL1_qual2; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
 604     fHOutMultV0MvsTRKqual2 = new TH2F("fHOutMultV0MvsTRK_qual2","fHOutMultV0MvsTRK_qual2; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
 605     fHOutMultTRKvsCL1qual2 = new TH2F("fHOutMultTRKvsCL1_qual2","fHOutMultTRKvsCL1_qual2; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
 606
 607     fHOutQuality = new TH1F("fHOutQuality", "fHOutQuality", 100,-0.5,99.5);
 608     fHOutVertex  = new TH1F("fHOutVertex", "fHOutVertex", 100,-20,20);
 609     fHOutVertexT0  = new TH1F("fHOutVertexT0", "fHOutVertexT0", 100,-20,20);
 610
 611     fOutputList->Add(  fHOutCentV0M     );
 612     fOutputList->Add(  fHOutCentV0MCVHN);
 613     fOutputList->Add(  fHOutCentV0MCVLN);
 614     fOutputList->Add(  fHOutCentV0MCVHNinMB);
 615     fOutputList->Add(  fHOutCentV0MCVLNinMB);
 616     fOutputList->Add(  fHOutCentV0MCCENT);
 617     fOutputList->Add(  fHOutCentV0MCSEMI);
 618     fOutputList->Add(  fHOutCentV0MCCENTinMB);
 619     fOutputList->Add(  fHOutCentV0MCSEMIinMB);
 620     fOutputList->Add(  fHOutCentV0MMSL    );
 621     fOutputList->Add(  fHOutCentV0MMSH    );
 622     fOutputList->Add(  fHOutCentV0MMUL    );
 623     fOutputList->Add(  fHOutCentV0MMLL    );
 624     fOutputList->Add(  fHOutCentV0MEJE    );
 625     fOutputList->Add(  fHOutCentV0MEGA    );
 626     fOutputList->Add(  fHOutCentV0MPHS   );
 627     fOutputList->Add(  fHOutCentV0MMSLinMB);
 628     fOutputList->Add(  fHOutCentV0MMSHinMB);
 629     fOutputList->Add(  fHOutCentV0MMULinMB);
 630     fOutputList->Add(  fHOutCentV0MMLLinMB);
 631     fOutputList->Add(  fHOutCentV0MEJEinMB);
 632     fOutputList->Add(  fHOutCentV0MEGAinMB);
 633     fOutputList->Add(  fHOutCentV0MPHSinMB);
 634     fOutputList->Add(  fHOutCentFMD     );
 635     fOutputList->Add(  fHOutCentTRK     );
 636     fOutputList->Add(  fHOutCentTKL     );
 637     fOutputList->Add(  fHOutCentCL0     );
 638     fOutputList->Add(  fHOutCentCL1     );
 639     fOutputList->Add(  fHOutCentV0MvsFMD);
 640     fOutputList->Add(  fHOutCentTKLvsV0M);
 641     fOutputList->Add(  fHOutCentZEMvsZDC);
 642     fOutputList->Add(  fHOutCentV0MvsCentCL1);
 643     fOutputList->Add(  fHOutCentV0MvsCentTRK);
 644     fOutputList->Add(  fHOutCentTRKvsCentCL1);
 645     fOutputList->Add(  fHOutCentV0MvsCentZDC);
 646     fOutputList->Add(  fHOutMultV0M);
 647     fOutputList->Add(  fHOutMultV0O);
 648     fOutputList->Add(  fHOutMultFMD);
 649     fOutputList->Add(  fHOutMultTRK);
 650     fOutputList->Add(  fHOutMultTKL);
 651     fOutputList->Add(  fHOutMultCL0);
 652     fOutputList->Add(  fHOutMultCL1);
 653     fOutputList->Add(  fHOutMultV0MvsZDN);
 654     fOutputList->Add(  fHOutMultZEMvsZDN);
 655     fOutputList->Add(  fHOutMultV0MvsZDC);
 656     fOutputList->Add(  fHOutMultZEMvsZDC);
 657     fOutputList->Add(  fHOutMultZEMvsZDCw);
 658     fOutputList->Add(  fHOutMultV0MvsCL1);
 659     fOutputList->Add(  fHOutMultV0MvsTRK);
 660     fOutputList->Add(  fHOutMultTRKvsCL1);
 661     fOutputList->Add(  fHOutMultV0MvsV0O);
 662     fOutputList->Add(  fHOutMultV0OvsCL1);
 663     fOutputList->Add(  fHOutMultV0OvsTRK);
 664     fOutputList->Add(  fHOutCentV0Mqual1 );
 665     fOutputList->Add(  fHOutCentTRKqual1 );
 666     fOutputList->Add(  fHOutCentCL1qual1 );
 667     fOutputList->Add(  fHOutMultV0MvsCL1qual1);
 668     fOutputList->Add(  fHOutMultV0MvsTRKqual1);
 669     fOutputList->Add(  fHOutMultTRKvsCL1qual1);
 670     fOutputList->Add(  fHOutCentV0Mqual2 );
 671     fOutputList->Add(  fHOutCentTRKqual2 );
 672     fOutputList->Add(  fHOutCentCL1qual2 );
 673     fOutputList->Add(  fHOutMultV0MvsCL1qual2);
 674     fOutputList->Add(  fHOutMultV0MvsTRKqual2);
 675     fOutputList->Add(  fHOutMultTRKvsCL1qual2);
 676     fOutputList->Add(  fHOutQuality );
 677     fOutputList->Add(  fHOutVertex );
 678     fOutputList->Add(  fHOutVertexT0 );
 679
 680     PostData(1, fOutputList);
 681   }
 682
 683   fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 684 }
 685
 686 //________________________________________________________________________
 687 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
 688 {
 689   // Execute analysis for current event:
 690   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
 691
 692   Float_t  zncEnergy = 0.;          //  ZNC Energy
 693   Float_t  zpcEnergy = 0.;          //  ZPC Energy
 694   Float_t  znaEnergy = 0.;          //  ZNA Energy
 695   Float_t  zpaEnergy = 0.;          //  ZPA Energy
 696   Float_t  zem1Energy = 0.;         //  ZEM1 Energy
 697   Float_t  zem2Energy = 0.;         //  ZEM2 Energy
 698   Bool_t   zdcEnergyCal = kFALSE;   // if zdc is calibrated (in pass2)
 699
 700   Int_t    nTracks = 0;             //  no. tracks
 701   Int_t    nTracklets = 0;          //  no. tracklets
 702   Int_t    nClusters[6] = {0};      //  no. clusters on 6 ITS layers
 703   Int_t    nChips[2];               //  no. chips on 2 SPD layers
 704   Float_t  spdCorr =0;              //  corrected spd2 multiplicity
 705
 706   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
 707   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
 708   Short_t  multV0AOnline  = 0;      //  multiplicity from V0 reco side A
 709   Short_t  multV0COnline  = 0;      //  multiplicity from V0 reco side C
 710   Float_t  v0Corr = 0;               // corrected V0 multiplicity (used for MC)
 711
 712   Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A
 713   Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C
 714
 715   Float_t zvtx =0;                  // z-vertex SPD
 716   Int_t zvtxNcont =0;               // contributors to z-vertex SPD
 717
 718   Float_t zvtxT0 =0;                // z-vertex T0
 719
 720
 721   AliCentrality *esdCent = 0;
 722
 723   if(fAnalysisInput.CompareTo("ESD")==0){
 724
 725     AliVEvent* event = InputEvent();
 726     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 727     if (!esd) {
 728       AliError("No ESD Event");
 729       return;
 730     }
 731
 732     LoadBranches();
 733
 734     if (SetupRun(esd)<0) {
 735       AliError("Centrality File not available for this run");
 736       return;
 737     }
 738
 739     esdCent = esd->GetCentrality();
 740
 741     // ***** V0 info
 742     AliESDVZERO* esdV0 = esd->GetVZEROData();
 743     multV0A=esdV0->GetMTotV0A();
 744     multV0C=esdV0->GetMTotV0C();
 745     v0Corr = multV0A+multV0C;
 746
 747     multV0AOnline=esdV0->GetTriggerChargeA();
 748     multV0COnline=esdV0->GetTriggerChargeC();
 749
 750     // ***** T0 info
 751     const AliESDTZERO* esdT0 = esd->GetESDTZERO();
 752     if (!esdT0)
 753       {
 754         AliError("AliESDTZERO not available");
 755       }
 756     Int_t trig=esdT0->GetT0Trig();
 757     Bool_t kT0BB = kFALSE;
 758     if(trig&1) kT0BB=kTRUE;
 759     zvtxT0=esdT0->GetT0zVertex();
 760
 761
 762     // ***** Trigger info
 763     fIsSelected = ((esdV0->GetV0ADecision()==1) && (esdV0->GetV0CDecision()==1));
 764     TString trigStr(esd->GetFiredTriggerClasses());
 765
 766     fCVHN=kFALSE; fCVLN=kFALSE; fCCENT=kFALSE; fCSEMI=kFALSE;
 767     fMSL=kFALSE;  fMSH=kFALSE;  fMUL=kFALSE;   fMLL=kFALSE;
 768     fEJE=kFALSE;  fEGA=kFALSE;  fPHS=kFALSE;
 769
 770     if ( (trigStr.Contains("-B-")) &&  (trigStr.Contains("CVHN")) && (fIsSelected))
 771       fCVHN=kTRUE;
 772     if ( (trigStr.Contains("-B-")) &&  (trigStr.Contains("CVLN")) && (fIsSelected))
 773       fCVLN=kTRUE;
 774     if ( (trigStr.Contains("-B-")) &&  (trigStr.Contains("CCENT")) && (fIsSelected))
 775       fCCENT=kTRUE;
 776     if ( (trigStr.Contains("-B-")) &&  (trigStr.Contains("CSEMI")) && (fIsSelected))
 777       fCSEMI=kTRUE;
 778
 779     if ( (trigStr.Contains("-B-")) &&  (trigStr.Contains("CPBI1MSL")) && (fIsSelected))
 780       fMSL=kTRUE;
 781     if ( (trigStr.Contains("-B-")) &&  (trigStr.Contains("CPBI1MSH")) && (fIsSelected))
 782       fMSH=kTRUE;
 783     if ( (trigStr.Contains("-B-")) &&  (trigStr.Contains("CPBI1MUL")) && (fIsSelected))
 784       fMUL=kTRUE;
 785     if ( (trigStr.Contains("-B-")) &&  (trigStr.Contains("CPBI1MLL")) && (fIsSelected))
 786       fMLL=kTRUE;
 787     if ( (trigStr.Contains("-B-")) &&  (trigStr.Contains("CPBI1EJE")) && (fIsSelected))
 788       fEJE=kTRUE;
 789     if ( (trigStr.Contains("-B-")) &&  (trigStr.Contains("CPBI1EGA")) && (fIsSelected))
 790       fEGA=kTRUE;
 791     if ( (trigStr.Contains("-B-")) &&  (trigStr.Contains("CPBI1PHS")) && (fIsSelected))
 792       fPHS=kTRUE;
 793
 794     fCVHNbit=kFALSE;    fCVLNbit=kFALSE;       fCCENTbit=kFALSE;    fCSEMIbit=kFALSE;
 795     if (esdV0->GetTriggerBits() & (1<<8))
 796       fCVHNbit=kTRUE;
 797     if (esdV0->GetTriggerBits() & (1<<6))
 798       fCVLNbit=kTRUE;
 799
 800     if (kT0BB && fCVHNbit)
 801       fCCENTbit=kTRUE;
 802     if (kT0BB && fCVLNbit)
 803       fCSEMIbit=kTRUE;
 804
 805
 806     // ***** Vertex Info
 807     const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
 808     zvtx        = vtxESD->GetZ();
 809     zvtxNcont   = vtxESD->GetNContributors();
 810
 811     // ***** CB info (tracklets, clusters, chips)
 812     //nTracks    = event->GetNumberOfTracks();
 813     nTracks    = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
 814
 815     const AliMultiplicity *mult = esd->GetMultiplicity();
 816
 817     nTracklets = mult->GetNumberOfTracklets();
 818
 819     for(Int_t ilay=0; ilay<6; ilay++){
 820       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
 821     }
 822
 823     for(Int_t ilay=0; ilay<2; ilay++){
 824       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
 825     }
 826
 827     spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);
 828
 829     // ***** FMD info
 830     AliESDFMD *fmd = esd->GetFMDData();
 831     Float_t totalMultA = 0;
 832     Float_t totalMultC = 0;
 833     const Float_t fFMDLowCut = 0.4;
 834
 835     for(UShort_t det=1;det<=3;det++) {
 836       Int_t nRings = (det==1 ? 1 : 2);
 837       for (UShort_t ir = 0; ir < nRings; ir++) {
 838         Char_t   ring = (ir == 0 ? 'I' : 'O');
 839         UShort_t nsec = (ir == 0 ? 20  : 40);
 840         UShort_t nstr = (ir == 0 ? 512 : 256);
 841         for(UShort_t sec =0; sec < nsec;  sec++)  {
 842           for(UShort_t strip = 0; strip < nstr; strip++) {
 843
 844             Float_t fmdMult = fmd->Multiplicity(det,ring,sec,strip);
 845             if(fmdMult == 0 || fmdMult == AliESDFMD::kInvalidMult) continue;
 846
 847             Float_t nParticles=0;
 848
 849             if(fmdMult > fFMDLowCut) {
 850               nParticles = 1.;
 851             }
 852
 853             if (det<3) totalMultA = totalMultA + nParticles;
 854             else totalMultC = totalMultC + nParticles;
 855
 856           }
 857         }
 858       }
 859     }
 860     multFMDA = totalMultA;
 861     multFMDC = totalMultC;
 862
 863     // ***** ZDC info
 864     AliESDZDC *esdZDC = esd->GetESDZDC();
 865     zdcEnergyCal = esdZDC->AliESDZDC::TestBit(AliESDZDC::kEnergyCalibratedSignal);
 866     if (zdcEnergyCal) {
 867       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
 868       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
 869       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
 870       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
 871     } else {
 872       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;
 873       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;
 874       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;
 875       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;
 876     }
 877     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;
 878     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;
 879
 880   }
 881
 882   else if(fAnalysisInput.CompareTo("AOD")==0){
 883     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
 884     // to be implemented
 885     printf("  AOD analysis not yet implemented!!!\n\n");
 886     return;
 887   }
 888
 889   // ***** Scaling for MC
 890   if (fIsMCInput) {
 891     fUseScaling=kFALSE;
 892     v0Corr  = Short_t((multV0A+multV0C)  * fV0MScaleFactorMC);
 893   }
 894   // ***** Scaling for Data
 895   if (fUseScaling) {
 896     v0Corr  = Short_t(v0Corr /   fV0MScaleFactor);
 897     spdCorr = spdCorr / fSPDScaleFactor;
 898     nTracks = Int_t(nTracks /   fTPCScaleFactor);
 899   }
 900
 901   // ***** Centrality Selection
 902   if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
 903   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
 904   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
 905   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
 906   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
 907   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
 908
 909   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
 910   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
 911   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin(zem1Energy+zem2Energy,zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
 912
 913
 914   // ***** Cleaning
 915   if (fUseCleaning) {
 916     fQuality=0;
 917
 918     // ***** vertex
 919     if (TMath::Abs(zvtx)>fZVCut || zvtxNcont<1) fQuality += 1;
 920
 921     // ***** outliers, skip in case of MC input
 922     if (!fIsMCInput) {
 923       // **** V0 vs SPD
 924       if (IsOutlierV0MSPD(spdCorr, v0Corr, int(fCentV0M))) fQuality  += 2;
 925       // ***** V0 vs TPC
 926       if (IsOutlierV0MTPC(nTracks, v0Corr, int(fCentV0M))) fQuality  += 4;
 927       // ***** V0 vs ZDC
 928       if (IsOutlierV0MZDC((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
 929           (zdcEnergyCal==kFALSE) ) fQuality  += 8;
 930       if (IsOutlierV0MZDCECal((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
 931           (zdcEnergyCal==kTRUE) ) fQuality  += 8;
 932     }
 933   } else {
 934     fQuality = 0;
 935   }
 936
 937
 938   if (esdCent) {
 939     esdCent->SetQuality(fQuality);
 940     esdCent->SetCentralityV0M(fCentV0M);
 941     esdCent->SetCentralityFMD(fCentFMD);
 942     esdCent->SetCentralityTRK(fCentTRK);
 943     esdCent->SetCentralityTKL(fCentTKL);
 944     esdCent->SetCentralityCL0(fCentCL0);
 945     esdCent->SetCentralityCL1(fCentCL1);
 946     esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
 947     esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
 948     esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
 949   }
 950
 951   // filling QA histograms
 952   if (fFillHistos) {
 953     if (fCVHN)   fHOutCentV0MCVHN->Fill(fCentV0M);
 954     if (fCVLN)   fHOutCentV0MCVLN->Fill(fCentV0M);
 955     if (fCCENT)  fHOutCentV0MCCENT->Fill(fCentV0M);
 956     if (fCSEMI)  fHOutCentV0MCSEMI->Fill(fCentV0M);
 957     if (fMSL) fHOutCentV0MMSL->Fill(fCentV0M);
 958     if (fMSH) fHOutCentV0MMSH->Fill(fCentV0M);
 959     if (fMUL) fHOutCentV0MMUL->Fill(fCentV0M);
 960     if (fMLL) fHOutCentV0MMLL->Fill(fCentV0M);
 961     if (fEJE) fHOutCentV0MEJE->Fill(fCentV0M);
 962     if (fEGA) fHOutCentV0MEGA->Fill(fCentV0M);
 963     if (fPHS) fHOutCentV0MPHS->Fill(fCentV0M);
 964
 965     if (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB) { // fill the QA histograms only for MB events!
 966       fHOutQuality->Fill(fQuality);
 967       fHOutVertex->Fill(zvtx);
 968       fHOutVertexT0->Fill(zvtxT0);
 969
 970       if (fQuality==0) {
 971         fHOutCentV0M->Fill(fCentV0M);
 972
 973         if (fCVHNbit)  fHOutCentV0MCVHNinMB->Fill(fCentV0M);
 974         if (fCVLNbit)  fHOutCentV0MCVLNinMB->Fill(fCentV0M);
 975         if (fCCENTbit) fHOutCentV0MCCENTinMB->Fill(fCentV0M);
 976         if (fCSEMIbit) fHOutCentV0MCSEMIinMB->Fill(fCentV0M);
 977         if (fMSL) fHOutCentV0MMSLinMB->Fill(fCentV0M);
 978         if (fMSH) fHOutCentV0MMSHinMB->Fill(fCentV0M);
 979         if (fMUL) fHOutCentV0MMULinMB->Fill(fCentV0M);
 980         if (fMLL) fHOutCentV0MMLLinMB->Fill(fCentV0M);
 981         if (fEJE) fHOutCentV0MEJEinMB->Fill(fCentV0M);
 982         if (fEGA) fHOutCentV0MEGAinMB->Fill(fCentV0M);
 983         if (fPHS) fHOutCentV0MPHSinMB->Fill(fCentV0M);
 984
 985
 986         fHOutCentFMD->Fill(fCentFMD);
 987         fHOutCentTRK->Fill(fCentTRK);
 988         fHOutCentTKL->Fill(fCentTKL);
 989         fHOutCentCL0->Fill(fCentCL0);
 990         fHOutCentCL1->Fill(fCentCL1);
 991         fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);
 992         fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);
 993         fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);
 994         fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentCL1);
 995         fHOutCentV0MvsCentTRK->Fill(fCentV0M,fCentTRK);
 996         fHOutCentTRKvsCentCL1->Fill(fCentTRK,fCentCL1);
 997         fHOutCentV0MvsCentZDC->Fill(fCentV0M,fCentZEMvsZDC);
 998         fHOutMultV0M->Fill(multV0A+multV0C);
 999         fHOutMultV0O->Fill(multV0AOnline+multV0COnline);
1000         fHOutMultFMD->Fill(multFMDA+multFMDC);
1001         fHOutMultTRK->Fill(nTracks);
1002         fHOutMultTKL->Fill(nTracklets);
1003         fHOutMultCL0->Fill(nClusters[0]);
1004         fHOutMultCL1->Fill(spdCorr);
1005         fHOutMultV0MvsZDN->Fill(v0Corr,(zncEnergy+znaEnergy));
1006         fHOutMultZEMvsZDN->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy));
1007         fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
1008         fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
1009         fHOutMultZEMvsZDCw->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy),fCentV0M);
1010         fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
1011         fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
1012         fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
1013         fHOutMultV0MvsV0O->Fill(v0Corr,(multV0AOnline+multV0COnline));
1014         fHOutMultV0OvsCL1->Fill((multV0AOnline+multV0COnline),spdCorr);
1015         fHOutMultV0OvsTRK->Fill((multV0AOnline+multV0COnline),nTracks);
1016       } else if (fQuality%2 == 0) {
1017         fHOutCentV0Mqual1->Fill(fCentV0M);
1018         fHOutCentTRKqual1->Fill(fCentTRK);
1019         fHOutCentCL1qual1->Fill(fCentCL1);
1020         fHOutMultV0MvsCL1qual1->Fill(v0Corr,spdCorr);
1021         fHOutMultV0MvsTRKqual1->Fill(v0Corr,nTracks);
1022         fHOutMultTRKvsCL1qual1->Fill(nTracks,spdCorr);
1023       } else {
1024         fHOutCentV0Mqual2->Fill(fCentV0M);
1025         fHOutCentTRKqual2->Fill(fCentTRK);
1026         fHOutCentCL1qual2->Fill(fCentCL1);
1027         fHOutMultV0MvsCL1qual2->Fill(v0Corr,spdCorr);
1028         fHOutMultV0MvsTRKqual2->Fill(v0Corr,nTracks);
1029         fHOutMultTRKvsCL1qual2->Fill(nTracks,spdCorr);
1030       }
1031     }
1032     PostData(1, fOutputList);
1033   }
1034 }
1035 //________________________________________________________________________
1036 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
1037 {
1038   // Terminate analysis
1039 }
1040 //________________________________________________________________________
1041 Int_t AliCentralitySelectionTask::SetupRun(const AliESDEvent* const esd)
1042 {
1043   // Setup files for run
1044
1045   if (!esd)
1046     return -1;
1047
1048   // check if something to be done
1049   if (fCurrentRun == esd->GetRunNumber())
1050     return 0;
1051   else
1052     fCurrentRun = esd->GetRunNumber();
1053
1054   TString fileName =(Form("%s/COMMON/CENTRALITY/data/centrality.root", AliAnalysisManager::GetOADBPath()));
1055   AliInfo(Form("Setup Centrality Selection for run %d with file %s\n",fCurrentRun,fileName.Data()));
1056
1057   AliOADBContainer *con = new AliOADBContainer("OADB");
1058   con->InitFromFile(fileName,"Centrality");
1059
1060   AliOADBCentrality*  centOADB = 0;
1061   centOADB = (AliOADBCentrality*)(con->GetObject(fCurrentRun));
1062   if (!centOADB) {
1063     AliWarning(Form("Centrality OADB does not exist for run %d, using Default \n",fCurrentRun ));
1064     centOADB  = (AliOADBCentrality*)(con->GetDefaultObject("oadbDefault"));
1065   }
1066
1067   // modes
1068   fUseScaling   = centOADB->UseScaling();
1069   fUseCleaning  = centOADB->UseCleaning();
1070
1071   // cuts
1072   fZVCut        = centOADB->ZVCut();
1073   fOutliersCut  = centOADB->OutliersCut();
1074
1075   // centrality histos
1076   fHtempV0M       = centOADB->V0hist();
1077   fHtempTRK       = centOADB->TPChist();
1078   fHtempCL1       = centOADB->SPDhist();
1079   fHtempZEMvsZDC  = centOADB->ZEMvsZDChist();
1080
1081   TString path = gSystem->ExpandPathName(fileName.Data());
1082   if (!fHtempV0M) AliWarning(Form("Calibration for V0M does not exist in %s", path.Data()));
1083   if (!fHtempTRK) AliWarning(Form("Calibration for TRK does not exist in %s", path.Data()));
1084   if (!fHtempCL1) AliWarning(Form("Calibration for CL1 does not exist in %s", path.Data()));
1085   if (!fHtempZEMvsZDC) AliWarning(Form("Calibration for ZEMvsZDC does not exist in %s", path.Data()));
1086
1087   // scale factors
1088   fV0MScaleFactor    = centOADB->V0MScaleFactor();
1089   fSPDScaleFactor    = centOADB->SPDScaleFactor();
1090   fTPCScaleFactor    = centOADB->TPCScaleFactor();
1091   fV0MScaleFactorMC  = centOADB->V0MScaleFactor();
1092
1093   // outliers parameters
1094   fV0MSPDOutlierPar0 = centOADB->V0MSPDOutlierPar0();
1095   fV0MSPDOutlierPar1 = centOADB->V0MSPDOutlierPar1();
1096   fV0MTPCOutlierPar0 = centOADB->V0MTPCOutlierPar0();
1097   fV0MTPCOutlierPar1 = centOADB->V0MTPCOutlierPar1();
1098
1099   fV0MSPDSigmaOutlierPar0 = centOADB->V0MSPDSigmaOutlierPar0();
1100   fV0MSPDSigmaOutlierPar1 = centOADB->V0MSPDSigmaOutlierPar1();
1101   fV0MSPDSigmaOutlierPar2 = centOADB->V0MSPDSigmaOutlierPar2();
1102   fV0MTPCSigmaOutlierPar0 = centOADB->V0MTPCSigmaOutlierPar0();
1103   fV0MTPCSigmaOutlierPar1 = centOADB->V0MTPCSigmaOutlierPar1();
1104   fV0MTPCSigmaOutlierPar2 = centOADB->V0MTPCSigmaOutlierPar2();
1105
1106   fV0MZDCOutlierPar0 =      centOADB->V0MZDCOutlierPar0();
1107   fV0MZDCOutlierPar1 =      centOADB->V0MZDCOutlierPar1();
1108   fV0MZDCEcalOutlierPar0 =  centOADB->V0MZDCEcalOutlierPar0();
1109   fV0MZDCEcalOutlierPar1 =  centOADB->V0MZDCEcalOutlierPar1();
1110
1111
1112
1113   return 0;
1114 }
1115
1116
1117
1118 //________________________________________________________________________
1119 Bool_t AliCentralitySelectionTask::IsOutlierV0MSPD(Float_t spd, Float_t v0, Int_t cent) const
1120 {
1121   // Clean outliers
1122   Float_t val = fV0MSPDOutlierPar0 +  fV0MSPDOutlierPar1 * v0;
1123   Float_t spdSigma = fV0MSPDSigmaOutlierPar0 + fV0MSPDSigmaOutlierPar1*cent + fV0MSPDSigmaOutlierPar2*cent*cent;
1124   if ( TMath::Abs(spd-val) > fOutliersCut*spdSigma )
1125     return kTRUE;
1126   else
1127     return kFALSE;
1128 }
1129
1130 //________________________________________________________________________
1131 Bool_t AliCentralitySelectionTask::IsOutlierV0MTPC(Int_t tracks, Float_t v0, Int_t cent) const
1132 {
1133   // Clean outliers
1134   Float_t val = fV0MTPCOutlierPar0 +  fV0MTPCOutlierPar1 * v0;
1135   Float_t tpcSigma = fV0MTPCSigmaOutlierPar0 + fV0MTPCSigmaOutlierPar1*cent + fV0MTPCSigmaOutlierPar2*cent*cent;
1136   if ( TMath::Abs(tracks-val) > fOutliersCut*tpcSigma )
1137     return kTRUE;
1138   else
1139     return kFALSE;
1140 }
1141
1142 //________________________________________________________________________
1143 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDC(Float_t zdc, Float_t v0) const
1144 {
1145   // Clean outliers
1146   Float_t val = fV0MZDCOutlierPar0 + fV0MZDCOutlierPar1 * v0;
1147   if (zdc >  val)
1148     return kTRUE;
1149   else
1150   return kFALSE;
1151 }
1152
1153 //________________________________________________________________________
1154 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDCECal(Float_t zdc, Float_t v0) const
1155 {
1156   // Clean outliers
1157   Float_t val = fV0MZDCEcalOutlierPar0 + fV0MZDCEcalOutlierPar1 * v0;
1158   if (zdc >  val)
1159     return kTRUE;
1160   else
1161     return kFALSE;
1162 }
1163
1164