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 "AliMultiplicity.h"
  53 #include "AliAODHandler.h"
  54 #include "AliAODHeader.h"
  55 #include "AliAODEvent.h"
  56 #include "AliAODVertex.h"
  57 #include "AliAODVZERO.h"
  58 #include "AliAODTracklets.h"
  59 #include "AliAODMCHeader.h"
  60 #include "AliMCEventHandler.h"
  61 #include "AliMCEvent.h"
  62 #include "AliAODMCParticle.h"
  63 #include "AliMCParticle.h"
  64 #include "AliStack.h"
  65 #include "AliAnalysisTaskSE.h"
  66 #include "AliGenEventHeader.h"
  67 #include "AliGenHijingEventHeader.h"
  68 #include "AliPhysicsSelectionTask.h"
  69 #include "AliPhysicsSelection.h"
  70 #include "AliBackgroundSelection.h"
  71 #include "AliESDUtils.h"
  72
  73 ClassImp(AliCentralitySelectionTask)
  74
  75
  76 //________________________________________________________________________
  77 AliCentralitySelectionTask::AliCentralitySelectionTask():
  78 AliAnalysisTaskSE(),
  79   fAnalysisInput("ESD"),
  80   fIsMCInput(kFALSE),
  81   fPass(0),
  82   fFile(0),
  83   fFile2(0),
  84   fCurrentRun(-1),
  85   fRunNo(-1),
  86   fLowRunN(0),
  87   fHighRunN(0),
  88   fUseScaling(0),
  89   fUseCleaning(0),
  90   fTrackCuts(0),
  91   fZVCut(10),
  92   fOutliersCut(5),
  93   fQuality(999),
  94   fCentV0M(0),
  95   fCentFMD(0),
  96   fCentTRK(0),
  97   fCentTKL(0),
  98   fCentCL0(0),
  99   fCentCL1(0),
 100   fCentV0MvsFMD(0),
 101   fCentTKLvsV0M(0),
 102   fCentZEMvsZDC(0),
 103   fHtempV0M(0),
 104   fHtempFMD(0),
 105   fHtempTRK(0),
 106   fHtempTKL(0),
 107   fHtempCL0(0),
 108   fHtempCL1(0),
 109   fHtempV0MvsFMD(0),
 110   fHtempTKLvsV0M(0),
 111   fHtempZEMvsZDC(0),
 112   fOutputList(0),
 113   fHOutCentV0M     (0),
 114   fHOutCentFMD     (0),
 115   fHOutCentTRK     (0),
 116   fHOutCentTKL     (0),
 117   fHOutCentCL0     (0),
 118   fHOutCentCL1     (0),
 119   fHOutCentV0MvsFMD(0),
 120   fHOutCentTKLvsV0M(0),
 121   fHOutCentZEMvsZDC(0),
 122   fHOutCentV0MvsCentCL1(0),
 123   fHOutCentV0MvsCentTRK(0),
 124   fHOutCentTRKvsCentCL1(0),
 125   fHOutCentV0MvsCentZDC(0),
 126   fHOutMultV0M(0),
 127   fHOutMultV0R(0),
 128   fHOutMultFMD(0),
 129   fHOutMultTRK(0),
 130   fHOutMultTKL(0),
 131   fHOutMultCL0(0),
 132   fHOutMultCL1(0),
 133   fHOutMultV0MvsZDN(0),
 134   fHOutMultZEMvsZDN(0),
 135   fHOutMultV0MvsZDC(0),
 136   fHOutMultZEMvsZDC(0),
 137   fHOutMultV0MvsCL1(0),
 138   fHOutMultV0MvsTRK(0),
 139   fHOutMultTRKvsCL1(0),
 140   fHOutCentV0Mqual1(0),
 141   fHOutCentTRKqual1(0),
 142   fHOutCentCL1qual1(0),
 143   fHOutCentV0Mqual2(0),
 144   fHOutCentTRKqual2(0),
 145   fHOutCentCL1qual2(0),
 146   fHOutQuality(0),
 147   fHOutVertex(0)
 148 {
 149   // Default constructor
 150   AliInfo("Centrality Selection enabled.");
 151   fLowRunN =136851;
 152   fHighRunN=139517;
 153
 154   for (Int_t i=0; i < 2667; i++) {
 155     fV0MScaleFactor[i]=0.0;
 156     fSPDScaleFactor[i]=0.0;
 157     fTPCScaleFactor[i]=0.0;
 158     fV0MScaleFactorMC[i]=0.0;
 159   }
 160   fUseScaling=kTRUE;
 161   fUseCleaning=kTRUE;
 162   fBranchNames="ESD:AliESDRun.,AliESDHeader.,AliESDZDC.,AliESDFMD.,AliESDVZERO."
 163                ",SPDVertex.,TPCVertex.,PrimaryVertex.,AliMultiplicity.,Tracks ";
 164 }
 165
 166 //________________________________________________________________________
 167 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
 168   AliAnalysisTaskSE(name),
 169   fAnalysisInput("ESD"),
 170   fIsMCInput(kFALSE),
 171   fPass(0),
 172   fFile(0),
 173   fFile2(0),
 174   fCurrentRun(-1),
 175   fRunNo(-1),
 176   fLowRunN(0),
 177   fHighRunN(0),
 178   fUseScaling(0),
 179   fUseCleaning(0),
 180   fTrackCuts(0),
 181   fZVCut(10),
 182   fOutliersCut(5),
 183   fQuality(999),
 184   fCentV0M(0),
 185   fCentFMD(0),
 186   fCentTRK(0),
 187   fCentTKL(0),
 188   fCentCL0(0),
 189   fCentCL1(0),
 190   fCentV0MvsFMD(0),
 191   fCentTKLvsV0M(0),
 192   fCentZEMvsZDC(0),
 193   fHtempV0M(0),
 194   fHtempFMD(0),
 195   fHtempTRK(0),
 196   fHtempTKL(0),
 197   fHtempCL0(0),
 198   fHtempCL1(0),
 199   fHtempV0MvsFMD(0),
 200   fHtempTKLvsV0M(0),
 201   fHtempZEMvsZDC(0),
 202   fOutputList(0),
 203   fHOutCentV0M     (0),
 204   fHOutCentFMD     (0),
 205   fHOutCentTRK     (0),
 206   fHOutCentTKL     (0),
 207   fHOutCentCL0     (0),
 208   fHOutCentCL1     (0),
 209   fHOutCentV0MvsFMD(0),
 210   fHOutCentTKLvsV0M(0),
 211   fHOutCentZEMvsZDC(0),
 212   fHOutCentV0MvsCentCL1(0),
 213   fHOutCentV0MvsCentTRK(0),
 214   fHOutCentTRKvsCentCL1(0),
 215   fHOutCentV0MvsCentZDC(0),
 216   fHOutMultV0M(0),
 217   fHOutMultV0R(0),
 218   fHOutMultFMD(0),
 219   fHOutMultTRK(0),
 220   fHOutMultTKL(0),
 221   fHOutMultCL0(0),
 222   fHOutMultCL1(0),
 223   fHOutMultV0MvsZDN(0),
 224   fHOutMultZEMvsZDN(0),
 225   fHOutMultV0MvsZDC(0),
 226   fHOutMultZEMvsZDC(0),
 227   fHOutMultV0MvsCL1(0),
 228   fHOutMultV0MvsTRK(0),
 229   fHOutMultTRKvsCL1(0),
 230   fHOutCentV0Mqual1(0),
 231   fHOutCentTRKqual1(0),
 232   fHOutCentCL1qual1(0),
 233   fHOutCentV0Mqual2(0),
 234   fHOutCentTRKqual2(0),
 235   fHOutCentCL1qual2(0),
 236   fHOutQuality(0),
 237   fHOutVertex(0)
 238 {
 239     // Default constructor
 240     fLowRunN =136851;
 241     fHighRunN=139517;
 242
 243   AliInfo("Centrality Selection enabled.");
 244   DefineOutput(1, TList::Class());
 245   for (Int_t i=0; i<2667; i++) {
 246     fV0MScaleFactor[i]=0.0;
 247     fSPDScaleFactor[i]=0.0;
 248     fTPCScaleFactor[i]=0.0;
 249     fV0MScaleFactorMC[i]=0.0;
 250   }
 251   fUseScaling=kTRUE;
 252   fUseCleaning=kTRUE;
 253   fBranchNames="ESD:AliESDRun.,AliESDHeader.,AliESDZDC.,AliESDFMD.,AliESDVZERO."
 254                ",SPDVertex.,TPCVertex.,PrimaryVertex.,AliMultiplicity.,Tracks ";
 255 }
 256
 257 //________________________________________________________________________
 258 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
 259 {
 260   // Assignment operator
 261   if (this!=&c) {
 262     AliAnalysisTaskSE::operator=(c);
 263   }
 264   return *this;
 265 }
 266
 267 //________________________________________________________________________
 268 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
 269   AliAnalysisTaskSE(ana),
 270   fAnalysisInput(ana.fDebug),
 271   fIsMCInput(ana.fIsMCInput),
 272   fPass(ana.fPass),
 273   fFile(ana.fFile),
 274   fFile2(ana.fFile2),
 275   fCurrentRun(ana.fCurrentRun),
 276   fRunNo(ana.fRunNo),
 277   fLowRunN(ana.fLowRunN),
 278   fHighRunN(ana.fHighRunN),
 279   fUseScaling(ana.fUseScaling),
 280   fUseCleaning(ana.fUseCleaning),
 281   fTrackCuts(ana.fTrackCuts),
 282   fZVCut(ana.fZVCut),
 283   fOutliersCut(ana.fOutliersCut),
 284   fQuality(ana.fQuality),
 285   fCentV0M(ana.fCentV0M),
 286   fCentFMD(ana.fCentFMD),
 287   fCentTRK(ana.fCentTRK),
 288   fCentTKL(ana.fCentTKL),
 289   fCentCL0(ana.fCentCL0),
 290   fCentCL1(ana.fCentCL1),
 291   fCentV0MvsFMD(ana.fCentV0MvsFMD),
 292   fCentTKLvsV0M(ana.fCentTKLvsV0M),
 293   fCentZEMvsZDC(ana.fCentZEMvsZDC),
 294   fHtempV0M(ana.fHtempV0M),
 295   fHtempFMD(ana.fHtempFMD),
 296   fHtempTRK(ana.fHtempTRK),
 297   fHtempTKL(ana.fHtempTKL),
 298   fHtempCL0(ana.fHtempCL0),
 299   fHtempCL1(ana.fHtempCL1),
 300   fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
 301   fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
 302   fHtempZEMvsZDC(ana.fHtempZEMvsZDC),
 303   fOutputList(ana.fOutputList),
 304   fHOutCentV0M     (ana.fHOutCentV0M     ),
 305   fHOutCentFMD     (ana.fHOutCentFMD     ),
 306   fHOutCentTRK     (ana.fHOutCentTRK     ),
 307   fHOutCentTKL     (ana.fHOutCentTKL     ),
 308   fHOutCentCL0     (ana.fHOutCentCL0     ),
 309   fHOutCentCL1     (ana.fHOutCentCL1     ),
 310   fHOutCentV0MvsFMD(ana.fHOutCentV0MvsFMD),
 311   fHOutCentTKLvsV0M(ana.fHOutCentTKLvsV0M),
 312   fHOutCentZEMvsZDC(ana.fHOutCentZEMvsZDC),
 313   fHOutCentV0MvsCentCL1(ana.fHOutCentV0MvsCentCL1),
 314   fHOutCentV0MvsCentTRK(ana.fHOutCentV0MvsCentTRK),
 315   fHOutCentTRKvsCentCL1(ana.fHOutCentTRKvsCentCL1),
 316   fHOutCentV0MvsCentZDC(ana.fHOutCentV0MvsCentZDC),
 317   fHOutMultV0M(ana.fHOutMultV0M),
 318   fHOutMultV0R(ana.fHOutMultV0R),
 319   fHOutMultFMD(ana.fHOutMultFMD),
 320   fHOutMultTRK(ana.fHOutMultTRK),
 321   fHOutMultTKL(ana.fHOutMultTKL),
 322   fHOutMultCL0(ana.fHOutMultCL0),
 323   fHOutMultCL1(ana.fHOutMultCL1),
 324   fHOutMultV0MvsZDN(ana.fHOutMultV0MvsZDN),
 325   fHOutMultZEMvsZDN(ana.fHOutMultZEMvsZDN),
 326   fHOutMultV0MvsZDC(ana.fHOutMultV0MvsZDC),
 327   fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC),
 328   fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1),
 329   fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK),
 330   fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1),
 331   fHOutCentV0Mqual1(ana.fHOutCentV0Mqual1),
 332   fHOutCentTRKqual1(ana.fHOutCentTRKqual1),
 333   fHOutCentCL1qual1(ana.fHOutCentCL1qual1),
 334   fHOutCentV0Mqual2(ana.fHOutCentV0Mqual2),
 335   fHOutCentTRKqual2(ana.fHOutCentTRKqual2),
 336   fHOutCentCL1qual2(ana.fHOutCentCL1qual2),
 337   fHOutQuality(ana.fHOutQuality),
 338   fHOutVertex(ana.fHOutVertex)
 339 {
 340   // Copy Constructor
 341     for (Int_t i=0; i<2667; i++) {
 342         fV0MScaleFactor[i]=0.0;
 343         fSPDScaleFactor[i]=0.0;
 344         fTPCScaleFactor[i]=0.0;
 345         fV0MScaleFactorMC[i]=0.0;
 346     }
 347
 348 }
 349
 350 //________________________________________________________________________
 351 AliCentralitySelectionTask::~AliCentralitySelectionTask()
 352 {
 353   // Destructor
 354   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
 355   if (fTrackCuts) delete fTrackCuts;
 356 }
 357
 358 //________________________________________________________________________
 359 void AliCentralitySelectionTask::UserCreateOutputObjects()
 360 {
 361   // Create the output containers
 362   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
 363   AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
 364
 365   fOutputList = new TList();
 366   fOutputList->SetOwner();
 367   fHOutCentV0M     = new TH1F("fHOutCentV0M","fHOutCentV0M; Centrality V0",505,0,101);
 368   fHOutCentFMD     = new TH1F("fHOutCentFMD","fHOutCentFMD; Centrality FMD",505,0,101);
 369   fHOutCentTRK     = new TH1F("fHOutCentTRK","fHOutCentTRK; Centrality TPC",505,0,101);
 370   fHOutCentTKL     = new TH1F("fHOutCentTKL","fHOutCentTKL; Centrality tracklets",505,0,101);
 371   fHOutCentCL0     = new TH1F("fHOutCentCL0","fHOutCentCL0; Centrality SPD inner",505,0,101);
 372   fHOutCentCL1     = new TH1F("fHOutCentCL1","fHOutCentCL1; Centrality SPD outer",505,0,101);
 373   fHOutCentV0MvsFMD= new TH1F("fHOutCentV0MvsFMD","fHOutCentV0MvsFMD; Centrality V0 vs FMD",505,0,101);
 374   fHOutCentTKLvsV0M= new TH1F("fHOutCentTKLvsV0M","fHOutCentTKLvsV0M; Centrality tracklets vs V0",505,0,101);
 375   fHOutCentZEMvsZDC= new TH1F("fHOutCentZEMvsZDC","fHOutCentZEMvsZDC; Centrality ZEM vs ZDC",505,0,101);
 376   fHOutCentV0MvsCentCL1= new TH2F("fHOutCentV0MvsCentCL1","fHOutCentV0MvsCentCL1; Cent V0; Cent SPD",505,0,101,505,0,101);
 377   fHOutCentV0MvsCentTRK= new TH2F("fHOutCentV0MvsCentTRK","fHOutCentV0MvsCentTRK; Cent V0; Cent TPC",505,0,101,505,0,101);
 378   fHOutCentTRKvsCentCL1= new TH2F("fHOutCentTRKvsCentCL1","fHOutCentTRKvsCentCL1; Cent TPC; Cent SPD",505,0,101,505,0,101);
 379   fHOutCentV0MvsCentZDC= new TH2F("fHOutCentV0MvsCentZDC","fHOutCentV0MvsCentZDC; Cent V0; Cent ZDC",505,0,101,505,0,101);
 380
 381   fHOutMultV0M = new TH1F("fHOutMultV0M","fHOutMultV0M; Multiplicity V0",25000,0,25000);
 382   fHOutMultV0R = new TH1F("fHOutMultV0R","fHOutMultV0R; Multiplicity V0",30000,0,30000);
 383   fHOutMultFMD = new TH1F("fHOutMultFMD","fHOutMultFMD; Multiplicity FMD",24000,0,24000);
 384   fHOutMultTRK = new TH1F("fHOutMultTRK","fHOutMultTRK; Multiplicity TPC",4000,0,4000);
 385   fHOutMultTKL = new TH1F("fHOutMultTKL","fHOutMultTKL; Multiplicity tracklets",5000,0,5000);
 386   fHOutMultCL0 = new TH1F("fHOutMultCL0","fHOutMultCL0; Multiplicity SPD inner",7000,0,7000);
 387   fHOutMultCL1 = new TH1F("fHOutMultCL1","fHOutMultCL1; Multiplicity SPD outer",7000,0,7000);
 388   fHOutMultV0MvsZDN = new TH2F("fHOutMultV0MvsZDN","fHOutMultV0MvsZDN; Multiplicity V0; Energy ZDC-N",500,0,25000,500,0,180000);
 389   fHOutMultZEMvsZDN = new TH2F("fHOutMultZEMvsZDN","fHOutMultZEMvsZDN; Energy ZEM; Energy ZDC-N",500,0,2500,500,0,180000);
 390   fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,25000,500,0,200000);
 391   fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,200000);
 392   fHOutMultV0MvsCL1 = new TH2F("fHOutMultV0MvsCL1","fHOutMultV0MvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
 393   fHOutMultV0MvsTRK = new TH2F("fHOutMultV0MvsTRK","fHOutMultV0MvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
 394   fHOutMultTRKvsCL1 = new TH2F("fHOutMultTRKvsCL1","fHOutMultTRKvsCL1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
 395
 396   fHOutCentV0Mqual1 = new TH1F("fHOutCentV0M_qual1","fHOutCentV0M_qual1; Centrality V0",505,0,101);
 397   fHOutCentTRKqual1 = new TH1F("fHOutCentTRK_qual1","fHOutCentTRK_qual1; Centrality TPC",505,0,101);
 398   fHOutCentCL1qual1 = new TH1F("fHOutCentCL1_qual1","fHOutCentCL1_qual1; Centrality SPD outer",505,0,101);
 399
 400   fHOutCentV0Mqual2 = new TH1F("fHOutCentV0M_qual2","fHOutCentV0M_qual2; Centrality V0",505,0,101);
 401   fHOutCentTRKqual2 = new TH1F("fHOutCentTRK_qual2","fHOutCentTRK_qual2; Centrality TPC",505,0,101);
 402   fHOutCentCL1qual2 = new TH1F("fHOutCentCL1_qual2","fHOutCentCL1_qual2; Centrality SPD outer",505,0,101);
 403
 404   fHOutQuality = new TH1F("fHOutQuality", "fHOutQuality", 100,-0.5,99.5);
 405   fHOutVertex  = new TH1F("fHOutVertex", "fHOutVertex", 100,-20,20);
 406
 407   fOutputList->Add(  fHOutCentV0M     );
 408   fOutputList->Add(  fHOutCentFMD     );
 409   fOutputList->Add(  fHOutCentTRK     );
 410   fOutputList->Add(  fHOutCentTKL     );
 411   fOutputList->Add(  fHOutCentCL0     );
 412   fOutputList->Add(  fHOutCentCL1     );
 413   fOutputList->Add(  fHOutCentV0MvsFMD);
 414   fOutputList->Add(  fHOutCentTKLvsV0M);
 415   fOutputList->Add(  fHOutCentZEMvsZDC);
 416   fOutputList->Add(  fHOutCentV0MvsCentCL1);
 417   fOutputList->Add(  fHOutCentV0MvsCentTRK);
 418   fOutputList->Add(  fHOutCentTRKvsCentCL1);
 419   fOutputList->Add(  fHOutCentV0MvsCentZDC);
 420   fOutputList->Add(  fHOutMultV0M);
 421   fOutputList->Add(  fHOutMultV0R);
 422   fOutputList->Add(  fHOutMultFMD);
 423   fOutputList->Add(  fHOutMultTRK);
 424   fOutputList->Add(  fHOutMultTKL);
 425   fOutputList->Add(  fHOutMultCL0);
 426   fOutputList->Add(  fHOutMultCL1);
 427   fOutputList->Add(  fHOutMultV0MvsZDN);
 428   fOutputList->Add(  fHOutMultZEMvsZDN);
 429   fOutputList->Add(  fHOutMultV0MvsZDC);
 430   fOutputList->Add(  fHOutMultZEMvsZDC);
 431   fOutputList->Add(  fHOutMultV0MvsCL1);
 432   fOutputList->Add(  fHOutMultV0MvsTRK);
 433   fOutputList->Add(  fHOutMultTRKvsCL1);
 434   fOutputList->Add(  fHOutCentV0Mqual1 );
 435   fOutputList->Add(  fHOutCentTRKqual1 );
 436   fOutputList->Add(  fHOutCentCL1qual1 );
 437   fOutputList->Add(  fHOutCentV0Mqual2 );
 438   fOutputList->Add(  fHOutCentTRKqual2 );
 439   fOutputList->Add(  fHOutCentCL1qual2 );
 440   fOutputList->Add(  fHOutQuality );
 441   fOutputList->Add(  fHOutVertex );
 442
 443
 444   fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 445
 446   PostData(1, fOutputList);
 447
 448   MyInitScaleFactor();
 449   if (fIsMCInput) MyInitScaleFactorMC();
 450
 451   if (fPass==0) AliFatal("Which pass are you analyzing? You should set it via taskCentrality->SetPass(N)");
 452 }
 453
 454 //________________________________________________________________________
 455 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
 456 {
 457   // Execute analysis for current event:
 458   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
 459
 460   Float_t  zncEnergy = 0.;          //  ZNC Energy
 461   Float_t  zpcEnergy = 0.;          //  ZPC Energy
 462   Float_t  znaEnergy = 0.;          //  ZNA Energy
 463   Float_t  zpaEnergy = 0.;          //  ZPA Energy
 464   Float_t  zem1Energy = 0.;         //  ZEM1 Energy
 465   Float_t  zem2Energy = 0.;         //  ZEM2 Energy
 466   Bool_t   zdcEnergyCal = kFALSE;   // if zdc is calibrated (in pass2)
 467
 468   Int_t    nTracks = 0;             //  no. tracks
 469   Int_t    nTracklets = 0;          //  no. tracklets
 470   Int_t    nClusters[6] = {0};      //  no. clusters on 6 ITS layers
 471   Int_t    nChips[2];               //  no. chips on 2 SPD layers
 472   Float_t  spdCorr =0;              //  corrected spd2 multiplicity
 473
 474   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
 475   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
 476   Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A
 477   Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C
 478
 479   Short_t v0Corr = 0;               // corrected V0 multiplicity
 480   Short_t v0CorrResc = 0;           // corrected and rescaled V0 multiplicity
 481
 482   Float_t zvtx =0;                  // z-vertex SPD
 483   Int_t zvtxNcont =0;               // contributors to z-vertex SPD
 484
 485
 486   AliCentrality *esdCent = 0;
 487
 488   if(fAnalysisInput.CompareTo("ESD")==0){
 489
 490     AliVEvent* event = InputEvent();
 491     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 492     if (!esd) {
 493         AliError("No ESD Event");
 494         return;
 495     }
 496
 497     LoadBranches();
 498
 499     if (fRunNo<=0) {
 500         if (SetupRun(esd)<0) {
 501             AliError("Centrality File not available for this run");
 502             return;
 503         }
 504     }
 505
 506     esdCent = esd->GetCentrality();
 507
 508     // ***** V0 info
 509     AliESDVZERO* esdV0 = esd->GetVZEROData();
 510     multV0A=esdV0->GetMTotV0A();
 511     multV0C=esdV0->GetMTotV0C();
 512
 513     Float_t v0CorrR;
 514     v0Corr = (Short_t)AliESDUtils::GetCorrV0(esd,v0CorrR);
 515     v0CorrResc = (Short_t)v0CorrR;
 516
 517     // ***** Vertex Info
 518     const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
 519     zvtx        = vtxESD->GetZ();
 520     zvtxNcont   = vtxESD->GetNContributors();
 521
 522     // ***** CB info (tracklets, clusters, chips)
 523     //nTracks    = event->GetNumberOfTracks();
 524     nTracks    = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
 525
 526     const AliMultiplicity *mult = esd->GetMultiplicity();
 527
 528     nTracklets = mult->GetNumberOfTracklets();
 529
 530     for(Int_t ilay=0; ilay<6; ilay++){
 531       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
 532     }
 533
 534     for(Int_t ilay=0; ilay<2; ilay++){
 535       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
 536     }
 537
 538     spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);
 539
 540     // ***** FMD info
 541     AliESDFMD *fmd = esd->GetFMDData();
 542     Float_t totalMultA = 0;
 543     Float_t totalMultC = 0;
 544     const Float_t fFMDLowCut = 0.4;
 545
 546     for(UShort_t det=1;det<=3;det++) {
 547       Int_t nRings = (det==1 ? 1 : 2);
 548       for (UShort_t ir = 0; ir < nRings; ir++) {
 549         Char_t   ring = (ir == 0 ? 'I' : 'O');
 550         UShort_t nsec = (ir == 0 ? 20  : 40);
 551         UShort_t nstr = (ir == 0 ? 512 : 256);
 552         for(UShort_t sec =0; sec < nsec;  sec++)  {
 553           for(UShort_t strip = 0; strip < nstr; strip++) {
 554
 555             Float_t fmdMult = fmd->Multiplicity(det,ring,sec,strip);
 556             if(fmdMult == 0 || fmdMult == AliESDFMD::kInvalidMult) continue;
 557
 558             Float_t nParticles=0;
 559
 560             if(fmdMult > fFMDLowCut) {
 561               nParticles = 1.;
 562             }
 563
 564             if (det<3) totalMultA = totalMultA + nParticles;
 565             else totalMultC = totalMultC + nParticles;
 566
 567           }
 568         }
 569       }
 570     }
 571     multFMDA = totalMultA;
 572     multFMDC = totalMultC;
 573
 574     // ***** ZDC info
 575     AliESDZDC *esdZDC = esd->GetESDZDC();
 576     zdcEnergyCal = esdZDC->AliESDZDC::TestBit(AliESDZDC::kEnergyCalibratedSignal);
 577     if (zdcEnergyCal) {
 578       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
 579       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
 580       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
 581       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
 582     } else {
 583       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;
 584       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;
 585       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;
 586       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;
 587     }
 588     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;
 589     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;
 590
 591   }
 592
 593   else if(fAnalysisInput.CompareTo("AOD")==0){
 594     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
 595     // to be implemented
 596     printf("  AOD analysis not yet implemented!!!\n\n");
 597     return;
 598   }
 599
 600   // ***** Scaling
 601   // ***** Scaling for pass2
 602   if (fPass==2) {
 603     fUseScaling=kFALSE;
 604   }
 605   // ***** Scaling for MC
 606   if (fIsMCInput) {
 607     fUseScaling=kFALSE;
 608     fUseCleaning=kFALSE;
 609     Float_t tempScalefactorV0M = MyGetScaleFactorMC(fCurrentRun);
 610     v0Corr  = Short_t((multV0A+multV0C)  * tempScalefactorV0M);
 611   }
 612   // ***** Scaling for Data
 613   if (fUseScaling) {
 614     Float_t tempScalefactorV0M = MyGetScaleFactor(fCurrentRun,0);
 615     Float_t tempScalefactorSPD = MyGetScaleFactor(fCurrentRun,1);
 616     Float_t tempScalefactorTPC = MyGetScaleFactor(fCurrentRun,2);
 617     v0Corr  = Short_t(v0Corr / tempScalefactorV0M);
 618     spdCorr = spdCorr / tempScalefactorSPD;
 619     nTracks = Int_t(nTracks / tempScalefactorTPC);
 620   }
 621
 622   // ***** Centrality Selection
 623   if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
 624   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
 625   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
 626   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
 627   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
 628   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
 629
 630   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
 631   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
 632   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin(zem1Energy+zem2Energy,zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
 633
 634   // ***** Cleaning
 635   if (fUseCleaning) {
 636       fQuality=0;
 637       fZVCut=10;
 638       fOutliersCut=6;
 639
 640       // ***** vertex
 641       if (TMath::Abs(zvtx)>fZVCut || zvtxNcont<1) fQuality += 1;
 642
 643       // ***** outliers
 644       // **** V0 vs SPD
 645       if (IsOutlierV0MSPD(spdCorr, v0Corr, int(fCentV0M))) fQuality  += 2;
 646       // ***** V0 vs TPC
 647       if (IsOutlierV0MTPC(nTracks, v0Corr, int(fCentV0M))) fQuality  += 4;
 648       // ***** V0 vs ZDC
 649        if (IsOutlierV0MZDC((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
 650            (zdcEnergyCal==kFALSE) && !(fIsMCInput)) fQuality  += 8;
 651        if (IsOutlierV0MZDCECal((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
 652            ((zdcEnergyCal==kTRUE) || (fIsMCInput))) fQuality  += 8;
 653   } else {
 654       fQuality = 0;
 655   }
 656
 657
 658   if (esdCent) {
 659       esdCent->SetQuality(fQuality);
 660       esdCent->SetCentralityV0M(fCentV0M);
 661       esdCent->SetCentralityFMD(fCentFMD);
 662       esdCent->SetCentralityTRK(fCentTRK);
 663       esdCent->SetCentralityTKL(fCentTKL);
 664       esdCent->SetCentralityCL0(fCentCL0);
 665       esdCent->SetCentralityCL1(fCentCL1);
 666       esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
 667       esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
 668       esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
 669   }
 670
 671   fHOutQuality->Fill(fQuality);
 672   fHOutVertex->Fill(zvtx);
 673
 674   fHOutMultV0MvsZDN->Fill(v0Corr,(zncEnergy+znaEnergy));
 675   fHOutMultZEMvsZDN->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy));
 676   fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
 677   fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
 678   fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
 679   fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
 680
 681   if (fQuality==0) {
 682     fHOutCentV0M->Fill(fCentV0M);
 683     fHOutCentFMD->Fill(fCentFMD);
 684     fHOutCentTRK->Fill(fCentTRK);
 685     fHOutCentTKL->Fill(fCentTKL);
 686     fHOutCentCL0->Fill(fCentCL0);
 687     fHOutCentCL1->Fill(fCentCL1);
 688     fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);
 689     fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);
 690     fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);
 691     fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentCL1);
 692     fHOutCentV0MvsCentTRK->Fill(fCentV0M,fCentTRK);
 693     fHOutCentTRKvsCentCL1->Fill(fCentTRK,fCentCL1);
 694     fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentZEMvsZDC);
 695     fHOutMultV0M->Fill(v0Corr);
 696     fHOutMultV0R->Fill(multV0A+multV0C);
 697     fHOutMultFMD->Fill((multFMDA+multFMDC));
 698     fHOutMultTRK->Fill(nTracks);
 699     fHOutMultTKL->Fill(nTracklets);
 700     fHOutMultCL0->Fill(nClusters[0]);
 701     fHOutMultCL1->Fill(spdCorr);
 702     fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
 703   } else if (fQuality ==1) {
 704     fHOutCentV0Mqual1->Fill(fCentV0M);
 705     fHOutCentTRKqual1->Fill(fCentTRK);
 706     fHOutCentCL1qual1->Fill(fCentCL1);
 707   } else {
 708     fHOutCentV0Mqual2->Fill(fCentV0M);
 709     fHOutCentTRKqual2->Fill(fCentTRK);
 710     fHOutCentCL1qual2->Fill(fCentCL1);
 711   }
 712
 713   PostData(1, fOutputList);
 714 }
 715
 716 //________________________________________________________________________
 717 void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename)
 718 {
 719   //  Read centrality histograms
 720   TDirectory *owd = gDirectory;
 721   // Check if the file is present
 722   TString path = gSystem->ExpandPathName(fCentfilename.Data());
 723   if (gSystem->AccessPathName(path)) {
 724      AliError(Form("File %s does not exist", path.Data()));
 725      return;
 726   }
 727   fFile  = TFile::Open(fCentfilename);
 728   owd->cd();
 729   fHtempV0M  = (TH1F*) (fFile->Get("hmultV0_percentile"));
 730   fHtempFMD  = (TH1F*) (fFile->Get("hmultFMD_percentile"));
 731   fHtempTRK  = (TH1F*) (fFile->Get("hNtracks_percentile"));
 732   fHtempTKL  = (TH1F*) (fFile->Get("hNtracklets_percentile"));
 733   fHtempCL0  = (TH1F*) (fFile->Get("hNclusters0_percentile"));
 734   fHtempCL1  = (TH1F*) (fFile->Get("hNclusters1_percentile"));
 735
 736   if (!fHtempV0M) AliWarning(Form("Calibration for V0M does not exist in %s", path.Data()));
 737   if (!fHtempFMD) AliWarning(Form("Calibration for FMD does not exist in %s", path.Data()));
 738   if (!fHtempTRK) AliWarning(Form("Calibration for TRK does not exist in %s", path.Data()));
 739   if (!fHtempTKL) AliWarning(Form("Calibration for TKL does not exist in %s", path.Data()));
 740   if (!fHtempCL0) AliWarning(Form("Calibration for CL0 does not exist in %s", path.Data()));
 741   if (!fHtempCL1) AliWarning(Form("Calibration for CL1 does not exist in %s", path.Data()));
 742
 743   owd->cd();
 744 }
 745
 746 //________________________________________________________________________
 747 void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2)
 748 {
 749   //  Read centrality histograms
 750   TDirectory *owd = gDirectory;
 751   TString path = gSystem->ExpandPathName(fCentfilename2.Data());
 752   if (gSystem->AccessPathName(path)) {
 753      AliError(Form("File %s does not exist", path.Data()));
 754      return;
 755   }
 756   fFile2  = TFile::Open(fCentfilename2);
 757   owd->cd();
 758   fHtempV0MvsFMD =  (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
 759   fHtempTKLvsV0M  = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
 760   fHtempZEMvsZDC  = (TH2F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
 761
 762   if (!fHtempV0MvsFMD) AliWarning(Form("Calibration for V0MvsFMD does not exist in %s", path.Data()));
 763   if (!fHtempTKLvsV0M) AliWarning(Form("Calibration for TKLvsV0M does not exist in %s", path.Data()));
 764   if (!fHtempZEMvsZDC) AliWarning(Form("Calibration for ZEMvsZDC does not exist in %s", path.Data()));
 765
 766   owd->cd();
 767 }
 768
 769 //________________________________________________________________________
 770 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
 771 {
 772   // Terminate analysis
 773   if (fFile && fFile->IsOpen())
 774     fFile->Close();
 775   if (fFile2 && fFile2->IsOpen())
 776     fFile2->Close();
 777 }
 778 //________________________________________________________________________
 779 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* const esd)
 780 {
 781   // Setup files for run
 782
 783   if (!esd)
 784     return -1;
 785
 786   // check if something to be done
 787   if (fCurrentRun == esd->GetRunNumber())
 788     return 0;
 789   else
 790     fCurrentRun = esd->GetRunNumber();
 791
 792   AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
 793
 794   // CHANGE HERE FOR RUN RANGES
 795   switch (fPass) {
 796     case 1:
 797       if ( fCurrentRun <= 137165 ) fRunNo = 137161;
 798       else fRunNo = 137366;
 799       break;
 800     case 2:
 801       if      ( fCurrentRun >= 136851  && fCurrentRun <= 137165 ) fRunNo = 137161;
 802       else if ( fCurrentRun >= 138190  && fCurrentRun <= 138275 ) fRunNo = 137722;
 803       else if ( fCurrentRun >= 137230  && fCurrentRun <= 137848 ) fRunNo = 138200;
 804       else if ( fCurrentRun >= 138125  && fCurrentRun <= 138154 ) fRunNo = 138150;
 805       else fRunNo = 139172;
 806       break;
 807     default:
 808       AliError(Form("Run %d not known to centrality selection!", fCurrentRun));
 809       return -1;
 810   }
 811   TString fileName =(Form("%s/COMMON/CENTRALITY/data/pass%d/AliCentralityBy1D_%d.root", AliAnalysisManager::GetOADBPath(), fPass, fRunNo));
 812   TString fileName2=(Form("%s/COMMON/CENTRALITY/data/pass%d/AliCentralityByFunction_%d.root", AliAnalysisManager::GetOADBPath(), fPass, fRunNo));
 813
 814   AliInfo(Form("Centrality Selection for run %d and pass %d is initialized with %s", fCurrentRun, fPass, fileName.Data()));
 815   ReadCentralityHistos(fileName.Data());
 816   ReadCentralityHistos2(fileName2.Data());
 817   if (!fFile && !fFile2) {
 818      AliError(Form("Run %d not known to centrality selection!", fCurrentRun));
 819      return -1;
 820   }
 821   return 0;
 822 }
 823
 824 //________________________________________________________________________
 825 Bool_t AliCentralitySelectionTask::IsOutlierV0MSPD(Float_t spd, Float_t v0, Int_t cent) const
 826 {
 827 // Clean outliers
 828   Float_t val= -0.579579 +  0.273949 * v0;
 829   Float_t spdSigma[101]={216.473, 194.377, 190.432, 185.019, 179.787, 175.99, 170.695, 167.276, 162.229, 159.016, 155.592, 151.909, 148.967, 146.367, 142.752, 139.718, 136.838, 134.133, 131.374, 128.616, 126.195, 123.49, 120.76, 118.368, 115.671, 113.354, 110.986, 108.761, 105.816, 103.373, 101.525, 99.3941, 96.8891, 95.1021, 92.728, 90.6434, 88.6393, 86.4624, 84.19, 82.3937, 80.4803, 78.5633, 76.5256, 74.7059, 73.0146, 70.9321, 69.1671, 67.2923, 65.5859, 63.9027, 61.7826, 60.1855, 58.6887, 56.8523, 55.1695, 53.422, 51.9776, 50.3506, 48.5304, 47.1517, 45.4786, 43.9914, 42.4579, 41.0889, 39.6201, 38.3046, 36.8624, 35.3697, 33.9223, 32.4637, 30.859, 29.9707, 28.764, 27.5677, 26.3437, 24.8656, 23.5253, 22.9878, 21.5172, 20.4645, 19.7614, 19.0363, 18.0875, 17.3675, 16.7313, 16.4101, 15.8235, 15.4795, 14.897, 14.578, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 14.3506, 22.5965, 22.5965};
 830
 831   if ( TMath::Abs(spd-val) > fOutliersCut*spdSigma[cent] )
 832     return kTRUE;
 833   else
 834     return kFALSE;
 835 }
 836
 837 //________________________________________________________________________
 838 Bool_t AliCentralitySelectionTask::IsOutlierV0MTPC(Int_t tracks, Float_t v0, Int_t cent) const
 839 {
 840 // Clean outliers
 841   Float_t val = -1.03873 +  0.125691 * v0;
 842   Float_t tpcSigma[101]={105.1, 95.7594, 94.5777, 91.7369, 89.655, 87.9469, 85.2747, 83.8137, 81.5663, 79.9871, 78.3491, 76.664, 75.1352, 73.9838, 72.0122, 70.7627, 69.1832, 67.8707, 66.3974, 65.0222, 64.1065, 62.7248, 61.498, 60.11, 58.7405, 57.681, 56.2979, 55.3153, 53.983, 52.6813, 51.5321, 50.5531, 49.3309, 48.2269, 47.1665, 46.1924, 45.1143, 44.0394, 42.8571, 41.8471, 40.9594, 39.8298, 38.8128, 37.8845, 36.8788, 35.9896, 34.9477, 34.1293, 33.1267, 32.161, 31.1132, 30.3521, 29.5237, 28.5992, 27.6283, 26.8885, 26.0719, 25.1334, 24.2125, 23.5376, 22.5383, 21.8396, 21.0384, 20.2942, 19.4746, 18.5545, 18.26, 17.2997, 16.7109, 15.9324, 15.2755, 14.764, 14.2358, 13.6728, 13.2606, 12.846, 12.2636, 11.722, 11.1491, 10.5952, 10.235, 9.87852, 9.48122, 9.17571, 8.90995, 8.81675, 8.62595, 8.59371, 8.52936, 8.58058, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 8.5543, 14.2584, 14.2584};
 843
 844   if ( TMath::Abs(tracks-val) > fOutliersCut*tpcSigma[cent] )
 845     return kTRUE;
 846   else
 847     return kFALSE;
 848 }
 849
 850 //________________________________________________________________________
 851 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDC(Float_t zdc, Float_t v0) const
 852 {
 853 // Clean outliers
 854   Float_t val = 235000. - 9.5 * v0;
 855   if (zdc >  val)
 856     return kTRUE;
 857   else
 858     return kFALSE;
 859 }
 860
 861 //________________________________________________________________________
 862 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDCECal(Float_t /*zdc*/, Float_t /*v0*/) const
 863 {
 864 // Clean outliers
 865     return kFALSE;
 866 }
 867
 868 //________________________________________________________________________
 869 Float_t AliCentralitySelectionTask::MyGetScaleFactor(Int_t runnumber, Int_t flag) const
 870 {
 871 // Get scaling factor
 872   if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
 873     cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
 874     return 0.0;
 875   }
 876
 877   Float_t scalefactor=0.0;
 878   if (flag==0)
 879     scalefactor = fV0MScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
 880   else if (flag==1)
 881     scalefactor = fSPDScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
 882   else if (flag==2)
 883     scalefactor = fTPCScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
 884
 885   return scalefactor;
 886
 887 }
 888
 889 //________________________________________________________________________
 890 Float_t AliCentralitySelectionTask::MyGetScaleFactorMC(Int_t runnumber) const
 891 {
 892 // Get MC scaling factor
 893   if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
 894     cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
 895     return 0.0;
 896   }
 897
 898   Float_t scalefactor= fV0MScaleFactorMC[runnumber - fLowRunN]; // subtracting reference offset index
 899   return scalefactor;
 900
 901 }
 902
 903 //________________________________________________________________________
 904 void AliCentralitySelectionTask::MyInitScaleFactor ()
 905 {
 906 // Initialize the scaling factors
 907   for (int i=0; i<=(fHighRunN-fLowRunN); i++) fV0MScaleFactor[i] = 0.0;
 908   for (int i=0; i<=(fHighRunN-fLowRunN); i++) fSPDScaleFactor[i] = 0.0;
 909   for (int i=0; i<=(fHighRunN-fLowRunN); i++) fTPCScaleFactor[i] = 0.0;
 910
 911   // scale factors determined from <V0 charge> on a run-by-run basis
 912    fV0MScaleFactor[310] = 1.;
 913    fV0MScaleFactor[311] = 1.;
 914    fV0MScaleFactor[514] = 1.0046;
 915    fV0MScaleFactor[515] = 0.983535;
 916    fV0MScaleFactor[579] = 0.988185;
 917    fV0MScaleFactor[580] = 0.983351;
 918    fV0MScaleFactor[581] = 0.989013;
 919    fV0MScaleFactor[583] = 0.990056;
 920    fV0MScaleFactor[588] = 0.974438;
 921    fV0MScaleFactor[589] = 0.981572;
 922    fV0MScaleFactor[590] = 0.989316;
 923    fV0MScaleFactor[592] = 0.98345;
 924    fV0MScaleFactor[688] = 0.993647;
 925    fV0MScaleFactor[690] = 0.994758;
 926    fV0MScaleFactor[693] = 0.989569;
 927    fV0MScaleFactor[698] = 0.993119;
 928    fV0MScaleFactor[744] = 0.989583;
 929    fV0MScaleFactor[757] = 0.990377;
 930    fV0MScaleFactor[787] = 0.990176;
 931    fV0MScaleFactor[788] = 0.98723;
 932    fV0MScaleFactor[834] = 1.00403;
 933    fV0MScaleFactor[835] = 0.994376;
 934    fV0MScaleFactor[840] = 0.99759;
 935    fV0MScaleFactor[842] = 1.01031;
 936    fV0MScaleFactor[900] = 0.996216;
 937    fV0MScaleFactor[901] = 0.994205;
 938    fV0MScaleFactor[992] = 0.998479;
 939    fV0MScaleFactor[997] = 1.00078;
 940    fV0MScaleFactor[1299] = 1.00515;
 941    fV0MScaleFactor[1303] = 1.00094;
 942    fV0MScaleFactor[1339] = 0.986596;
 943    fV0MScaleFactor[1346] = 0.972226;
 944    fV0MScaleFactor[1349] = 0.960358;
 945    fV0MScaleFactor[1350] = 0.970023;
 946    fV0MScaleFactor[1374] = 1.00575;
 947    fV0MScaleFactor[1545] = 1.00471;
 948    fV0MScaleFactor[1587] = 1.00611;
 949    fV0MScaleFactor[1588] = 1.00976;
 950    fV0MScaleFactor[1618] = 1.00771;
 951    fV0MScaleFactor[1682] = 1.01622;
 952    fV0MScaleFactor[1727] = 1.01305;
 953    fV0MScaleFactor[1728] = 1.00388;
 954    fV0MScaleFactor[1731] = 1.00673;
 955    fV0MScaleFactor[1732] = 1.00916;
 956    fV0MScaleFactor[1770] = 1.0092;
 957    fV0MScaleFactor[1773] = 1.00728;
 958    fV0MScaleFactor[1786] = 1.01655;
 959    fV0MScaleFactor[1787] = 1.00672;
 960    fV0MScaleFactor[1801] = 0.983339;
 961    fV0MScaleFactor[1802] = 1.00754;
 962    fV0MScaleFactor[1811] = 1.00608;
 963    fV0MScaleFactor[1815] = 1.01227;
 964    fV0MScaleFactor[1879] = 0.99907;
 965    fV0MScaleFactor[1881] = 0.995696;
 966    fV0MScaleFactor[2186] = 1.00559;
 967    fV0MScaleFactor[2187] = 1.00631;
 968    fV0MScaleFactor[2254] = 1.01512;
 969    fV0MScaleFactor[2256] = 0.998727;
 970    fV0MScaleFactor[2321] = 1.00701;
 971    fSPDScaleFactor[310] = 1.00211;
 972    fSPDScaleFactor[311] = 1.00067;
 973    fSPDScaleFactor[514] = 1.02268;
 974    fSPDScaleFactor[515] = 0.994902;
 975    fSPDScaleFactor[579] = 1.00215;
 976    fSPDScaleFactor[580] = 0.993421;
 977    fSPDScaleFactor[581] = 1.00129;
 978    fSPDScaleFactor[583] = 1.00242;
 979    fSPDScaleFactor[588] = 0.984762;
 980    fSPDScaleFactor[589] = 0.994355;
 981    fSPDScaleFactor[590] = 1.00073;
 982    fSPDScaleFactor[592] = 0.995889;
 983    fSPDScaleFactor[688] = 0.994532;
 984    fSPDScaleFactor[690] = 0.998307;
 985    fSPDScaleFactor[693] = 0.994052;
 986    fSPDScaleFactor[698] = 0.993224;
 987    fSPDScaleFactor[744] = 0.993279;
 988    fSPDScaleFactor[757] = 0.992494;
 989    fSPDScaleFactor[787] = 0.992678;
 990    fSPDScaleFactor[788] = 0.996563;
 991    fSPDScaleFactor[834] = 1.01116;
 992    fSPDScaleFactor[835] = 0.993108;
 993    fSPDScaleFactor[840] = 0.997574;
 994    fSPDScaleFactor[842] = 1.01829;
 995    fSPDScaleFactor[900] = 0.999438;
 996    fSPDScaleFactor[901] = 0.995849;
 997    fSPDScaleFactor[992] = 0.999227;
 998    fSPDScaleFactor[997] = 1.00575;
 999    fSPDScaleFactor[1299] = 0.99877;
1000    fSPDScaleFactor[1303] = 0.999682;
1001    fSPDScaleFactor[1339] = 0.978198;
1002    fSPDScaleFactor[1346] = 0.964178;
1003    fSPDScaleFactor[1349] = 0.959439;
1004    fSPDScaleFactor[1350] = 0.956945;
1005    fSPDScaleFactor[1374] = 0.994434;
1006    fSPDScaleFactor[1545] = 1.0016;
1007    fSPDScaleFactor[1587] = 1.00153;
1008    fSPDScaleFactor[1588] = 1.00698;
1009    fSPDScaleFactor[1618] = 1.00554;
1010    fSPDScaleFactor[1682] = 1.0123;
1011    fSPDScaleFactor[1727] = 1.011;
1012    fSPDScaleFactor[1728] = 1.00143;
1013    fSPDScaleFactor[1731] = 1.00486;
1014    fSPDScaleFactor[1732] = 1.00646;
1015    fSPDScaleFactor[1770] = 1.00515;
1016    fSPDScaleFactor[1773] = 1.00485;
1017    fSPDScaleFactor[1786] = 1.01215;
1018    fSPDScaleFactor[1787] = 1.00419;
1019    fSPDScaleFactor[1801] = 0.983327;
1020    fSPDScaleFactor[1802] = 1.00529;
1021    fSPDScaleFactor[1811] = 1.00367;
1022    fSPDScaleFactor[1815] = 1.01045;
1023    fSPDScaleFactor[1879] = 0.996374;
1024    fSPDScaleFactor[1881] = 0.988827;
1025    fSPDScaleFactor[2186] = 1.00354;
1026    fSPDScaleFactor[2187] = 1.00397;
1027    fSPDScaleFactor[2254] = 1.01138;
1028    fSPDScaleFactor[2256] = 0.996641;
1029    fSPDScaleFactor[2321] = 1.00357;
1030    fTPCScaleFactor[310] = 1.00434;
1031    fTPCScaleFactor[311] = 1.0056;
1032    fTPCScaleFactor[514] = 1.02185;
1033    fTPCScaleFactor[515] = 1.0036;
1034    fTPCScaleFactor[579] = 1.00607;
1035    fTPCScaleFactor[580] = 1.00183;
1036    fTPCScaleFactor[581] = 1.00693;
1037    fTPCScaleFactor[583] = 1.00746;
1038    fTPCScaleFactor[588] = 0.990524;
1039    fTPCScaleFactor[589] = 0.998582;
1040    fTPCScaleFactor[590] = 1.00618;
1041    fTPCScaleFactor[592] = 1.00088;
1042    fTPCScaleFactor[688] = 1.00598;
1043    fTPCScaleFactor[690] = 1.00658;
1044    fTPCScaleFactor[693] = 1.00144;
1045    fTPCScaleFactor[698] = 1.00279;
1046    fTPCScaleFactor[744] = 1.00122;
1047    fTPCScaleFactor[757] = 1.002;
1048    fTPCScaleFactor[787] = 0.997818;
1049    fTPCScaleFactor[788] = 0.994583;
1050    fTPCScaleFactor[834] = 1.01508;
1051    fTPCScaleFactor[835] = 1.00218;
1052    fTPCScaleFactor[840] = 1.00569;
1053    fTPCScaleFactor[842] = 1.01789;
1054    fTPCScaleFactor[900] = 1.00739;
1055    fTPCScaleFactor[901] = 1.00462;
1056    fTPCScaleFactor[992] = 1.00502;
1057    fTPCScaleFactor[997] = 1.00943;
1058    fTPCScaleFactor[1299] = 1.00438;
1059    fTPCScaleFactor[1303] = 0.996701;
1060    fTPCScaleFactor[1339] = 0.978641;
1061    fTPCScaleFactor[1346] = 0.968906;
1062    fTPCScaleFactor[1349] = 0.954311;
1063    fTPCScaleFactor[1350] = 0.958764;
1064    fTPCScaleFactor[1374] = 0.997899;
1065    fTPCScaleFactor[1545] = 0.992;
1066    fTPCScaleFactor[1587] = 0.992635;
1067    fTPCScaleFactor[1588] = 1.00207;
1068    fTPCScaleFactor[1618] = 1.00126;
1069    fTPCScaleFactor[1682] = 1.00324;
1070    fTPCScaleFactor[1727] = 1.00042;
1071    fTPCScaleFactor[1728] = 0.978881;
1072    fTPCScaleFactor[1731] = 0.999818;
1073    fTPCScaleFactor[1732] = 1.00109;
1074    fTPCScaleFactor[1770] = 0.99935;
1075    fTPCScaleFactor[1773] = 0.998531;
1076    fTPCScaleFactor[1786] = 0.999125;
1077    fTPCScaleFactor[1787] = 0.998479;
1078    fTPCScaleFactor[1801] = 0.9775;
1079    fTPCScaleFactor[1802] = 0.999095;
1080    fTPCScaleFactor[1811] = 0.998197;
1081    fTPCScaleFactor[1815] = 1.00413;
1082    fTPCScaleFactor[1879] = 0.990916;
1083    fTPCScaleFactor[1881] = 0.987241;
1084    fTPCScaleFactor[2186] = 1.00048;
1085    fTPCScaleFactor[2187] = 1.00057;
1086    fTPCScaleFactor[2254] = 1.00588;
1087    fTPCScaleFactor[2256] = 0.991503;
1088    fTPCScaleFactor[2321] = 1.00057;
1089
1090
1091   // set all missing values to the value of the run before it ....
1092   for (int i=0; i<=(fHighRunN-fLowRunN); i++) {
1093     if (fV0MScaleFactor[i] == 0.0) {
1094       if (i==0) {
1095         fV0MScaleFactor[i] = 1.00;
1096       } else {
1097         // search for last run number with non-zero value
1098         for (int j=i-1;j>=0;j--) {
1099           if (fV0MScaleFactor[j] != 0.0) {
1100             fV0MScaleFactor[i] = fV0MScaleFactor[j];
1101             break;
1102           }
1103         }
1104       }
1105     }
1106   } // end loop over checking all run-numbers
1107
1108   for (int i=0; i<=(fHighRunN-fLowRunN); i++) {
1109     if (fSPDScaleFactor[i] == 0.0) {
1110       if (i==0) {
1111         fSPDScaleFactor[i] = 1.00;
1112       } else {
1113         for (int j=i-1;j>=0;j--) {
1114           if (fSPDScaleFactor[j] != 0.0) {
1115             fSPDScaleFactor[i] = fSPDScaleFactor[j];
1116             break;
1117           }
1118         }
1119       }
1120     }
1121   }
1122
1123   for (int i=0; i<=(fHighRunN-fLowRunN); i++) {
1124     if (fTPCScaleFactor[i] == 0.0) {
1125       if (i==0) {
1126         fTPCScaleFactor[i] = 1.00;
1127       } else {
1128         for (int j=i-1;j>=0;j--) {
1129           if (fTPCScaleFactor[j] != 0.0) {
1130             fTPCScaleFactor[i] = fTPCScaleFactor[j];
1131             break;
1132           }
1133         }
1134       }
1135     }
1136   }
1137
1138
1139   //    for (int i=0; i<=(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << fV0MScaleFactor[i] << " for Run " << i+fLowRunN << endl;
1140
1141   return;
1142
1143 }
1144
1145
1146 //________________________________________________________________________
1147 void AliCentralitySelectionTask::MyInitScaleFactorMC()
1148 {
1149 // Initialize the MC scaling factors
1150   for (int i=0; i<(fHighRunN-fLowRunN); i++) fV0MScaleFactorMC[i] = 0.0;
1151   // scale factors determined from <V0 charge> on a run-by-run basis
1152   switch (fPass) {
1153   case 1:
1154     fV0MScaleFactorMC[0] = 0.75108;
1155     break;
1156   case 2:
1157     fV0MScaleFactorMC[0] = 0.8;
1158     break;
1159   default:
1160     AliError(Form("Unknown reconstruction pass (%d)! Setting MC scale in V0 to 1", fPass));
1161     fV0MScaleFactorMC[0] = 1.0;
1162     break;
1163   }
1164
1165   // set all missing values to the value of the run before it ....
1166   for (int i=0; i<=(fHighRunN-fLowRunN); i++) {
1167     if (fV0MScaleFactorMC[i] == 0.0) {
1168       if (i==0) {
1169         fV0MScaleFactorMC[i] = 1.00;
1170       } else {
1171         // search for last run number with non-zero value
1172         for (int j=i-1;j>=0;j--) {
1173           if (fV0MScaleFactorMC[j] != 0.0) {
1174             fV0MScaleFactorMC[i] = fV0MScaleFactorMC[j];
1175             break;
1176           }
1177         }
1178       }
1179     }
1180   } // end loop over checking all run-numbers
1181
1182
1183   //    for (int i=0; i<=(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << fV0MScaleFactorMC[i] << " for Run " << i+fLowRunN << endl;
1184
1185   return;
1186
1187 }
1188