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