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