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