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