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