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