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