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   // Execute analysis for current event:
 488   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
 489
 490   Float_t  zncEnergy = 0.;          //  ZNC Energy
 491   Float_t  zpcEnergy = 0.;          //  ZPC Energy
 492   Float_t  znaEnergy = 0.;          //  ZNA Energy
 493   Float_t  zpaEnergy = 0.;          //  ZPA Energy
 494   Float_t  zem1Energy = 0.;         //  ZEM1 Energy
 495   Float_t  zem2Energy = 0.;         //  ZEM2 Energy
 496   Bool_t   zdcEnergyCal = kFALSE;   // if zdc is calibrated (in pass2)
 497
 498   Int_t    nTracks = 0;             //  no. tracks
 499   Int_t    nTracklets = 0;          //  no. tracklets
 500   Int_t    nClusters[6] = {0};      //  no. clusters on 6 ITS layers
 501   Int_t    nChips[2];               //  no. chips on 2 SPD layers
 502   Float_t  spdCorr =0;              //  corrected spd2 multiplicity
 503
 504   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
 505   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
 506   Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A
 507   Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C
 508
 509   Short_t v0Corr = 0;               // corrected V0 multiplicity
 510   Short_t v0CorrResc = 0;           // corrected and rescaled V0 multiplicity
 511
 512   Float_t zvtx =0;                  // z-vertex SPD
 513   Int_t zvtxNcont =0;               // contributors to z-vertex SPD
 514
 515
 516   AliCentrality *esdCent = 0;
 517
 518   if(fAnalysisInput.CompareTo("ESD")==0){
 519
 520     AliVEvent* event = InputEvent();
 521     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 522     if (!esd) {
 523         AliError("No ESD Event");
 524         return;
 525     }
 526
 527     LoadBranches();
 528
 529     if (fRunNo<=0) {
 530         if (SetupRun(esd)<0) {
 531             AliError("Centrality File not available for this run");
 532             return;
 533         }
 534     }
 535
 536     esdCent = esd->GetCentrality();
 537
 538     // ***** V0 info
 539     AliESDVZERO* esdV0 = esd->GetVZEROData();
 540     multV0A=esdV0->GetMTotV0A();
 541     multV0C=esdV0->GetMTotV0C();
 542
 543     Float_t v0CorrR;
 544     v0Corr = (Short_t)AliESDUtils::GetCorrV0(esd,v0CorrR);
 545     v0CorrResc = (Short_t)v0CorrR;
 546
 547     // ***** Vertex Info
 548     const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
 549     zvtx        = vtxESD->GetZ();
 550     zvtxNcont   = vtxESD->GetNContributors();
 551
 552     // ***** CB info (tracklets, clusters, chips)
 553     //nTracks    = event->GetNumberOfTracks();
 554     nTracks    = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
 555
 556     const AliMultiplicity *mult = esd->GetMultiplicity();
 557
 558     nTracklets = mult->GetNumberOfTracklets();
 559
 560     for(Int_t ilay=0; ilay<6; ilay++){
 561       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
 562     }
 563
 564     for(Int_t ilay=0; ilay<2; ilay++){
 565       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
 566     }
 567
 568     spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);
 569
 570     // ***** FMD info
 571     AliESDFMD *fmd = esd->GetFMDData();
 572     Float_t totalMultA = 0;
 573     Float_t totalMultC = 0;
 574     const Float_t fFMDLowCut = 0.4;
 575
 576     for(UShort_t det=1;det<=3;det++) {
 577       Int_t nRings = (det==1 ? 1 : 2);
 578       for (UShort_t ir = 0; ir < nRings; ir++) {
 579         Char_t   ring = (ir == 0 ? 'I' : 'O');
 580         UShort_t nsec = (ir == 0 ? 20  : 40);
 581         UShort_t nstr = (ir == 0 ? 512 : 256);
 582         for(UShort_t sec =0; sec < nsec;  sec++)  {
 583           for(UShort_t strip = 0; strip < nstr; strip++) {
 584
 585             Float_t fmdMult = fmd->Multiplicity(det,ring,sec,strip);
 586             if(fmdMult == 0 || fmdMult == AliESDFMD::kInvalidMult) continue;
 587
 588             Float_t nParticles=0;
 589
 590             if(fmdMult > fFMDLowCut) {
 591               nParticles = 1.;
 592             }
 593
 594             if (det<3) totalMultA = totalMultA + nParticles;
 595             else totalMultC = totalMultC + nParticles;
 596
 597           }
 598         }
 599       }
 600     }
 601     multFMDA = totalMultA;
 602     multFMDC = totalMultC;
 603
 604     // ***** ZDC info
 605     AliESDZDC *esdZDC = esd->GetESDZDC();
 606     zdcEnergyCal = esdZDC->AliESDZDC::TestBit(AliESDZDC::kEnergyCalibratedSignal);
 607     if (zdcEnergyCal) {
 608       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
 609       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
 610       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
 611       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
 612     } else {
 613       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;
 614       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;
 615       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;
 616       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;
 617     }
 618     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;
 619     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;
 620
 621   }
 622
 623   else if(fAnalysisInput.CompareTo("AOD")==0){
 624     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
 625     // to be implemented
 626     printf("  AOD analysis not yet implemented!!!\n\n");
 627     return;
 628   }
 629
 630   // ***** Scaling
 631   // ***** Scaling for pass2
 632   if (fPass==2) {
 633     fUseScaling=kFALSE;
 634   }
 635   // ***** Scaling for MC
 636   if (fIsMCInput) {
 637     fUseScaling=kFALSE;
 638     fUseCleaning=kFALSE;
 639     Float_t tempScalefactorV0M = MyGetScaleFactorMC(fCurrentRun);
 640     v0Corr  = Short_t((multV0A+multV0C)  * tempScalefactorV0M);
 641   }
 642   // ***** Scaling for Data
 643   if (fUseScaling) {
 644     Float_t tempScalefactorV0M = MyGetScaleFactor(fCurrentRun,0);
 645     Float_t tempScalefactorSPD = MyGetScaleFactor(fCurrentRun,1);
 646     Float_t tempScalefactorTPC = MyGetScaleFactor(fCurrentRun,2);
 647     v0Corr  = Short_t(v0Corr / tempScalefactorV0M);
 648     spdCorr = spdCorr / tempScalefactorSPD;
 649     nTracks = Int_t(nTracks / tempScalefactorTPC);
 650   }
 651
 652   // ***** Centrality Selection
 653   if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
 654   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
 655   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
 656   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
 657   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
 658   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
 659
 660   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
 661   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
 662   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin(zem1Energy+zem2Energy,zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
 663
 664   // ***** Cleaning
 665   if (fUseCleaning) {
 666       fQuality=0;
 667       fZVCut=10;
 668       fOutliersCut=6;
 669
 670       // ***** vertex
 671       if (TMath::Abs(zvtx)>fZVCut || zvtxNcont<1) fQuality += 1;
 672
 673       // ***** outliers
 674       // **** V0 vs SPD
 675       if (IsOutlierV0MSPD(spdCorr, v0Corr, int(fCentV0M))) fQuality  += 2;
 676       // ***** V0 vs TPC
 677       if (IsOutlierV0MTPC(nTracks, v0Corr, int(fCentV0M))) fQuality  += 4;
 678       // ***** V0 vs ZDC
 679        if (IsOutlierV0MZDC((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
 680            (zdcEnergyCal==kFALSE) && !(fIsMCInput)) fQuality  += 8;
 681        if (IsOutlierV0MZDCECal((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
 682            ((zdcEnergyCal==kTRUE) || (fIsMCInput))) fQuality  += 8;
 683   } else {
 684       fQuality = 0;
 685   }
 686
 687
 688   if (esdCent) {
 689       esdCent->SetQuality(fQuality);
 690       esdCent->SetCentralityV0M(fCentV0M);
 691       esdCent->SetCentralityFMD(fCentFMD);
 692       esdCent->SetCentralityTRK(fCentTRK);
 693       esdCent->SetCentralityTKL(fCentTKL);
 694       esdCent->SetCentralityCL0(fCentCL0);
 695       esdCent->SetCentralityCL1(fCentCL1);
 696       esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
 697       esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
 698       esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
 699   }
 700
 701   fHOutQuality->Fill(fQuality);
 702   fHOutVertex->Fill(zvtx);
 703
 704   if (fQuality==0) {
 705     fHOutCentV0M->Fill(fCentV0M);
 706     fHOutCentFMD->Fill(fCentFMD);
 707     fHOutCentTRK->Fill(fCentTRK);
 708     fHOutCentTKL->Fill(fCentTKL);
 709     fHOutCentCL0->Fill(fCentCL0);
 710     fHOutCentCL1->Fill(fCentCL1);
 711     fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);
 712     fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);
 713     fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);
 714     fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentCL1);
 715     fHOutCentV0MvsCentTRK->Fill(fCentV0M,fCentTRK);
 716     fHOutCentTRKvsCentCL1->Fill(fCentTRK,fCentCL1);
 717     fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentZEMvsZDC);
 718     fHOutMultV0M->Fill(v0Corr);
 719     fHOutMultV0R->Fill(multV0A+multV0C);
 720     fHOutMultFMD->Fill((multFMDA+multFMDC));
 721     fHOutMultTRK->Fill(nTracks);
 722     fHOutMultTKL->Fill(nTracklets);
 723     fHOutMultCL0->Fill(nClusters[0]);
 724     fHOutMultCL1->Fill(spdCorr);
 725     fHOutMultV0MvsZDN->Fill(v0Corr,(zncEnergy+znaEnergy));
 726     fHOutMultZEMvsZDN->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy));
 727     fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
 728     fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
 729     fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
 730     fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
 731     fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
 732   } else if (fQuality%2 == 0) {
 733     fHOutCentV0Mqual1->Fill(fCentV0M);
 734     fHOutCentTRKqual1->Fill(fCentTRK);
 735     fHOutCentCL1qual1->Fill(fCentCL1);
 736     fHOutMultV0MvsCL1qual1->Fill(v0Corr,spdCorr);
 737     fHOutMultV0MvsTRKqual1->Fill(v0Corr,nTracks);
 738     fHOutMultTRKvsCL1qual1->Fill(nTracks,spdCorr);
 739   } else {
 740     fHOutCentV0Mqual2->Fill(fCentV0M);
 741     fHOutCentTRKqual2->Fill(fCentTRK);
 742     fHOutCentCL1qual2->Fill(fCentCL1);
 743     fHOutMultV0MvsCL1qual2->Fill(v0Corr,spdCorr);
 744     fHOutMultV0MvsTRKqual2->Fill(v0Corr,nTracks);
 745     fHOutMultTRKvsCL1qual2->Fill(nTracks,spdCorr);
 746   }
 747
 748   PostData(1, fOutputList);
 749 }
 750
 751 //________________________________________________________________________
 752 void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename)
 753 {
 754   //  Read centrality histograms
 755   TDirectory *owd = gDirectory;
 756   // Check if the file is present
 757   TString path = gSystem->ExpandPathName(fCentfilename.Data());
 758   if (gSystem->AccessPathName(path)) {
 759      AliError(Form("File %s does not exist", path.Data()));
 760      return;
 761   }
 762   fFile  = TFile::Open(fCentfilename);
 763   owd->cd();
 764   fHtempV0M  = (TH1F*) (fFile->Get("hmultV0_percentile"));
 765   fHtempFMD  = (TH1F*) (fFile->Get("hmultFMD_percentile"));
 766   fHtempTRK  = (TH1F*) (fFile->Get("hNtracks_percentile"));
 767   fHtempTKL  = (TH1F*) (fFile->Get("hNtracklets_percentile"));
 768   fHtempCL0  = (TH1F*) (fFile->Get("hNclusters0_percentile"));
 769   fHtempCL1  = (TH1F*) (fFile->Get("hNclusters1_percentile"));
 770
 771   if (!fHtempV0M) AliWarning(Form("Calibration for V0M does not exist in %s", path.Data()));
 772   if (!fHtempFMD) AliWarning(Form("Calibration for FMD does not exist in %s", path.Data()));
 773   if (!fHtempTRK) AliWarning(Form("Calibration for TRK does not exist in %s", path.Data()));
 774   if (!fHtempTKL) AliWarning(Form("Calibration for TKL does not exist in %s", path.Data()));
 775   if (!fHtempCL0) AliWarning(Form("Calibration for CL0 does not exist in %s", path.Data()));
 776   if (!fHtempCL1) AliWarning(Form("Calibration for CL1 does not exist in %s", path.Data()));
 777
 778   owd->cd();
 779 }
 780
 781 //________________________________________________________________________
 782 void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2)
 783 {
 784   //  Read centrality histograms
 785   TDirectory *owd = gDirectory;
 786   TString path = gSystem->ExpandPathName(fCentfilename2.Data());
 787   if (gSystem->AccessPathName(path)) {
 788      AliError(Form("File %s does not exist", path.Data()));
 789      return;
 790   }
 791   fFile2  = TFile::Open(fCentfilename2);
 792   owd->cd();
 793   fHtempV0MvsFMD =  (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
 794   fHtempTKLvsV0M  = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
 795   fHtempZEMvsZDC  = (TH2F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
 796
 797   if (!fHtempV0MvsFMD) AliWarning(Form("Calibration for V0MvsFMD does not exist in %s", path.Data()));
 798   if (!fHtempTKLvsV0M) AliWarning(Form("Calibration for TKLvsV0M does not exist in %s", path.Data()));
 799   if (!fHtempZEMvsZDC) AliWarning(Form("Calibration for ZEMvsZDC does not exist in %s", path.Data()));
 800
 801   owd->cd();
 802 }
 803
 804 //________________________________________________________________________
 805 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
 806 {
 807   // Terminate analysis
 808   if (fFile && fFile->IsOpen())
 809     fFile->Close();
 810   if (fFile2 && fFile2->IsOpen())
 811     fFile2->Close();
 812 }
 813 //________________________________________________________________________
 814 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* const esd)
 815 {
 816   // Setup files for run
 817
 818   if (!esd)
 819     return -1;
 820
 821   // check if something to be done
 822   if (fCurrentRun == esd->GetRunNumber())
 823     return 0;
 824   else
 825     fCurrentRun = esd->GetRunNumber();
 826
 827   AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
 828
 829   // CHANGE HERE FOR RUN RANGES
 830   switch (fPass) {
 831     case 1:
 832       if ( fCurrentRun <= 137165 ) fRunNo = 137161;
 833       else fRunNo = 137366;
 834       break;
 835     case 2:
 836       if      ( fCurrentRun >= 136851  && fCurrentRun <= 137165 ) fRunNo = 137161;
 837       else if ( fCurrentRun >= 137230  && fCurrentRun <= 137531 ) fRunNo = 137366;
 838       else if ( fCurrentRun >= 137539  && fCurrentRun <= 137848 ) fRunNo = 137722;
 839       else if ( fCurrentRun >= 138125  && fCurrentRun <= 138154 ) fRunNo = 138150;
 840       else if ( fCurrentRun >= 138190  && fCurrentRun <= 138275 ) fRunNo = 138200;
 841       else fRunNo = 139172;
 842       break;
 843     default:
 844       AliError(Form("Run %d not known to centrality selection!", fCurrentRun));
 845       return -1;
 846   }
 847   TString fileName =(Form("%s/COMMON/CENTRALITY/data/pass%d/AliCentralityBy1D_%d.root", AliAnalysisManager::GetOADBPath(), fPass, fRunNo));
 848   TString fileName2=(Form("%s/COMMON/CENTRALITY/data/pass%d/AliCentralityByFunction_%d.root", AliAnalysisManager::GetOADBPath(), fPass, fRunNo));
 849
 850   AliInfo(Form("Centrality Selection for run %d and pass %d is initialized with %s", fCurrentRun, fPass, fileName.Data()));
 851   ReadCentralityHistos(fileName.Data());
 852   ReadCentralityHistos2(fileName2.Data());
 853   if (!fFile && !fFile2) {
 854      AliError(Form("Run %d not known to centrality selection!", fCurrentRun));
 855      return -1;
 856   }
 857   return 0;
 858 }
 859
 860 //________________________________________________________________________
 861 Bool_t AliCentralitySelectionTask::IsOutlierV0MSPD(Float_t spd, Float_t v0, Int_t cent) const
 862 {
 863 // Clean outliers
 864   Float_t val= -0.579579 +  0.273949 * v0;
 865   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};
 866
 867   if ( TMath::Abs(spd-val) > fOutliersCut*spdSigma[cent] )
 868     return kTRUE;
 869   else
 870     return kFALSE;
 871 }
 872
 873 //________________________________________________________________________
 874 Bool_t AliCentralitySelectionTask::IsOutlierV0MTPC(Int_t tracks, Float_t v0, Int_t cent) const
 875 {
 876 // Clean outliers
 877   Float_t val = -1.03873 +  0.125691 * v0;
 878   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};
 879
 880   if ( TMath::Abs(tracks-val) > fOutliersCut*tpcSigma[cent] )
 881     return kTRUE;
 882   else
 883     return kFALSE;
 884 }
 885
 886 //________________________________________________________________________
 887 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDC(Float_t zdc, Float_t v0) const
 888 {
 889 // Clean outliers
 890   Float_t val = 235000. - 9.5 * v0;
 891   if (zdc >  val)
 892     return kTRUE;
 893   else
 894     return kFALSE;
 895 }
 896
 897 //________________________________________________________________________
 898 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDCECal(Float_t /*zdc*/, Float_t /*v0*/) const
 899 {
 900 // Clean outliers
 901     return kFALSE;
 902 }
 903
 904 //________________________________________________________________________
 905 Float_t AliCentralitySelectionTask::MyGetScaleFactor(Int_t runnumber, Int_t flag) const
 906 {
 907 // Get scaling factor
 908   if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
 909     cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
 910     return 0.0;
 911   }
 912
 913   Float_t scalefactor=0.0;
 914   if (flag==0)
 915     scalefactor = fV0MScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
 916   else if (flag==1)
 917     scalefactor = fSPDScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
 918   else if (flag==2)
 919     scalefactor = fTPCScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
 920
 921   return scalefactor;
 922
 923 }
 924
 925 //________________________________________________________________________
 926 Float_t AliCentralitySelectionTask::MyGetScaleFactorMC(Int_t runnumber) const
 927 {
 928 // Get MC scaling factor
 929   if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
 930     cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
 931     return 0.0;
 932   }
 933
 934   Float_t scalefactor= fV0MScaleFactorMC[runnumber - fLowRunN]; // subtracting reference offset index
 935   return scalefactor;
 936
 937 }
 938
 939 //________________________________________________________________________
 940 void AliCentralitySelectionTask::MyInitScaleFactor ()
 941 {
 942 // Initialize the scaling factors
 943   for (int i=0; i<=(fHighRunN-fLowRunN); i++) fV0MScaleFactor[i] = 0.0;
 944   for (int i=0; i<=(fHighRunN-fLowRunN); i++) fSPDScaleFactor[i] = 0.0;
 945   for (int i=0; i<=(fHighRunN-fLowRunN); i++) fTPCScaleFactor[i] = 0.0;
 946
 947   // scale factors determined from <V0 charge> on a run-by-run basis
 948    fV0MScaleFactor[310] = 1.;
 949    fV0MScaleFactor[311] = 1.;
 950    fV0MScaleFactor[514] = 1.0046;
 951    fV0MScaleFactor[515] = 0.983535;
 952    fV0MScaleFactor[579] = 0.988185;
 953    fV0MScaleFactor[580] = 0.983351;
 954    fV0MScaleFactor[581] = 0.989013;
 955    fV0MScaleFactor[583] = 0.990056;
 956    fV0MScaleFactor[588] = 0.974438;
 957    fV0MScaleFactor[589] = 0.981572;
 958    fV0MScaleFactor[590] = 0.989316;
 959    fV0MScaleFactor[592] = 0.98345;
 960    fV0MScaleFactor[688] = 0.993647;
 961    fV0MScaleFactor[690] = 0.994758;
 962    fV0MScaleFactor[693] = 0.989569;
 963    fV0MScaleFactor[698] = 0.993119;
 964    fV0MScaleFactor[744] = 0.989583;
 965    fV0MScaleFactor[757] = 0.990377;
 966    fV0MScaleFactor[787] = 0.990176;
 967    fV0MScaleFactor[788] = 0.98723;
 968    fV0MScaleFactor[834] = 1.00403;
 969    fV0MScaleFactor[835] = 0.994376;
 970    fV0MScaleFactor[840] = 0.99759;
 971    fV0MScaleFactor[842] = 1.01031;
 972    fV0MScaleFactor[900] = 0.996216;
 973    fV0MScaleFactor[901] = 0.994205;
 974    fV0MScaleFactor[992] = 0.998479;
 975    fV0MScaleFactor[997] = 1.00078;
 976    fV0MScaleFactor[1299] = 1.00515;
 977    fV0MScaleFactor[1303] = 1.00094;
 978    fV0MScaleFactor[1339] = 0.986596;
 979    fV0MScaleFactor[1346] = 0.972226;
 980    fV0MScaleFactor[1349] = 0.960358;
 981    fV0MScaleFactor[1350] = 0.970023;
 982    fV0MScaleFactor[1374] = 1.00575;
 983    fV0MScaleFactor[1545] = 1.00471;
 984    fV0MScaleFactor[1587] = 1.00611;
 985    fV0MScaleFactor[1588] = 1.00976;
 986    fV0MScaleFactor[1618] = 1.00771;
 987    fV0MScaleFactor[1682] = 1.01622;
 988    fV0MScaleFactor[1727] = 1.01305;
 989    fV0MScaleFactor[1728] = 1.00388;
 990    fV0MScaleFactor[1731] = 1.00673;
 991    fV0MScaleFactor[1732] = 1.00916;
 992    fV0MScaleFactor[1770] = 1.0092;
 993    fV0MScaleFactor[1773] = 1.00728;
 994    fV0MScaleFactor[1786] = 1.01655;
 995    fV0MScaleFactor[1787] = 1.00672;
 996    fV0MScaleFactor[1801] = 0.983339;
 997    fV0MScaleFactor[1802] = 1.00754;
 998    fV0MScaleFactor[1811] = 1.00608;
 999    fV0MScaleFactor[1815] = 1.01227;
1000    fV0MScaleFactor[1879] = 0.99907;
1001    fV0MScaleFactor[1881] = 0.995696;
1002    fV0MScaleFactor[2186] = 1.00559;
1003    fV0MScaleFactor[2187] = 1.00631;
1004    fV0MScaleFactor[2254] = 1.01512;
1005    fV0MScaleFactor[2256] = 0.998727;
1006    fV0MScaleFactor[2321] = 1.00701;
1007    fSPDScaleFactor[310] = 1.00211;
1008    fSPDScaleFactor[311] = 1.00067;
1009    fSPDScaleFactor[514] = 1.02268;
1010    fSPDScaleFactor[515] = 0.994902;
1011    fSPDScaleFactor[579] = 1.00215;
1012    fSPDScaleFactor[580] = 0.993421;
1013    fSPDScaleFactor[581] = 1.00129;
1014    fSPDScaleFactor[583] = 1.00242;
1015    fSPDScaleFactor[588] = 0.984762;
1016    fSPDScaleFactor[589] = 0.994355;
1017    fSPDScaleFactor[590] = 1.00073;
1018    fSPDScaleFactor[592] = 0.995889;
1019    fSPDScaleFactor[688] = 0.994532;
1020    fSPDScaleFactor[690] = 0.998307;
1021    fSPDScaleFactor[693] = 0.994052;
1022    fSPDScaleFactor[698] = 0.993224;
1023    fSPDScaleFactor[744] = 0.993279;
1024    fSPDScaleFactor[757] = 0.992494;
1025    fSPDScaleFactor[787] = 0.992678;
1026    fSPDScaleFactor[788] = 0.996563;
1027    fSPDScaleFactor[834] = 1.01116;
1028    fSPDScaleFactor[835] = 0.993108;
1029    fSPDScaleFactor[840] = 0.997574;
1030    fSPDScaleFactor[842] = 1.01829;
1031    fSPDScaleFactor[900] = 0.999438;
1032    fSPDScaleFactor[901] = 0.995849;
1033    fSPDScaleFactor[992] = 0.999227;
1034    fSPDScaleFactor[997] = 1.00575;
1035    fSPDScaleFactor[1299] = 0.99877;
1036    fSPDScaleFactor[1303] = 0.999682;
1037    fSPDScaleFactor[1339] = 0.978198;
1038    fSPDScaleFactor[1346] = 0.964178;
1039    fSPDScaleFactor[1349] = 0.959439;
1040    fSPDScaleFactor[1350] = 0.956945;
1041    fSPDScaleFactor[1374] = 0.994434;
1042    fSPDScaleFactor[1545] = 1.0016;
1043    fSPDScaleFactor[1587] = 1.00153;
1044    fSPDScaleFactor[1588] = 1.00698;
1045    fSPDScaleFactor[1618] = 1.00554;
1046    fSPDScaleFactor[1682] = 1.0123;
1047    fSPDScaleFactor[1727] = 1.011;
1048    fSPDScaleFactor[1728] = 1.00143;
1049    fSPDScaleFactor[1731] = 1.00486;
1050    fSPDScaleFactor[1732] = 1.00646;
1051    fSPDScaleFactor[1770] = 1.00515;
1052    fSPDScaleFactor[1773] = 1.00485;
1053    fSPDScaleFactor[1786] = 1.01215;
1054    fSPDScaleFactor[1787] = 1.00419;
1055    fSPDScaleFactor[1801] = 0.983327;
1056    fSPDScaleFactor[1802] = 1.00529;
1057    fSPDScaleFactor[1811] = 1.00367;
1058    fSPDScaleFactor[1815] = 1.01045;
1059    fSPDScaleFactor[1879] = 0.996374;
1060    fSPDScaleFactor[1881] = 0.988827;
1061    fSPDScaleFactor[2186] = 1.00354;
1062    fSPDScaleFactor[2187] = 1.00397;
1063    fSPDScaleFactor[2254] = 1.01138;
1064    fSPDScaleFactor[2256] = 0.996641;
1065    fSPDScaleFactor[2321] = 1.00357;
1066    fTPCScaleFactor[310] = 1.00434;
1067    fTPCScaleFactor[311] = 1.0056;
1068    fTPCScaleFactor[514] = 1.02185;
1069    fTPCScaleFactor[515] = 1.0036;
1070    fTPCScaleFactor[579] = 1.00607;
1071    fTPCScaleFactor[580] = 1.00183;
1072    fTPCScaleFactor[581] = 1.00693;
1073    fTPCScaleFactor[583] = 1.00746;
1074    fTPCScaleFactor[588] = 0.990524;
1075    fTPCScaleFactor[589] = 0.998582;
1076    fTPCScaleFactor[590] = 1.00618;
1077    fTPCScaleFactor[592] = 1.00088;
1078    fTPCScaleFactor[688] = 1.00598;
1079    fTPCScaleFactor[690] = 1.00658;
1080    fTPCScaleFactor[693] = 1.00144;
1081    fTPCScaleFactor[698] = 1.00279;
1082    fTPCScaleFactor[744] = 1.00122;
1083    fTPCScaleFactor[757] = 1.002;
1084    fTPCScaleFactor[787] = 0.997818;
1085    fTPCScaleFactor[788] = 0.994583;
1086    fTPCScaleFactor[834] = 1.01508;
1087    fTPCScaleFactor[835] = 1.00218;
1088    fTPCScaleFactor[840] = 1.00569;
1089    fTPCScaleFactor[842] = 1.01789;
1090    fTPCScaleFactor[900] = 1.00739;
1091    fTPCScaleFactor[901] = 1.00462;
1092    fTPCScaleFactor[992] = 1.00502;
1093    fTPCScaleFactor[997] = 1.00943;
1094    fTPCScaleFactor[1299] = 1.00438;
1095    fTPCScaleFactor[1303] = 0.996701;
1096    fTPCScaleFactor[1339] = 0.978641;
1097    fTPCScaleFactor[1346] = 0.968906;
1098    fTPCScaleFactor[1349] = 0.954311;
1099    fTPCScaleFactor[1350] = 0.958764;
1100    fTPCScaleFactor[1374] = 0.997899;
1101    fTPCScaleFactor[1545] = 0.992;
1102    fTPCScaleFactor[1587] = 0.992635;
1103    fTPCScaleFactor[1588] = 1.00207;
1104    fTPCScaleFactor[1618] = 1.00126;
1105    fTPCScaleFactor[1682] = 1.00324;
1106    fTPCScaleFactor[1727] = 1.00042;
1107    fTPCScaleFactor[1728] = 0.978881;
1108    fTPCScaleFactor[1731] = 0.999818;
1109    fTPCScaleFactor[1732] = 1.00109;
1110    fTPCScaleFactor[1770] = 0.99935;
1111    fTPCScaleFactor[1773] = 0.998531;
1112    fTPCScaleFactor[1786] = 0.999125;
1113    fTPCScaleFactor[1787] = 0.998479;
1114    fTPCScaleFactor[1801] = 0.9775;
1115    fTPCScaleFactor[1802] = 0.999095;
1116    fTPCScaleFactor[1811] = 0.998197;
1117    fTPCScaleFactor[1815] = 1.00413;
1118    fTPCScaleFactor[1879] = 0.990916;
1119    fTPCScaleFactor[1881] = 0.987241;
1120    fTPCScaleFactor[2186] = 1.00048;
1121    fTPCScaleFactor[2187] = 1.00057;
1122    fTPCScaleFactor[2254] = 1.00588;
1123    fTPCScaleFactor[2256] = 0.991503;
1124    fTPCScaleFactor[2321] = 1.00057;
1125
1126
1127   // set all missing values to the value of the run before it ....
1128   for (int i=0; i<=(fHighRunN-fLowRunN); i++) {
1129     if (fV0MScaleFactor[i] == 0.0) {
1130       if (i==0) {
1131         fV0MScaleFactor[i] = 1.00;
1132       } else {
1133         // search for last run number with non-zero value
1134         for (int j=i-1;j>=0;j--) {
1135           if (fV0MScaleFactor[j] != 0.0) {
1136             fV0MScaleFactor[i] = fV0MScaleFactor[j];
1137             break;
1138           }
1139         }
1140       }
1141     }
1142   } // end loop over checking all run-numbers
1143
1144   for (int i=0; i<=(fHighRunN-fLowRunN); i++) {
1145     if (fSPDScaleFactor[i] == 0.0) {
1146       if (i==0) {
1147         fSPDScaleFactor[i] = 1.00;
1148       } else {
1149         for (int j=i-1;j>=0;j--) {
1150           if (fSPDScaleFactor[j] != 0.0) {
1151             fSPDScaleFactor[i] = fSPDScaleFactor[j];
1152             break;
1153           }
1154         }
1155       }
1156     }
1157   }
1158
1159   for (int i=0; i<=(fHighRunN-fLowRunN); i++) {
1160     if (fTPCScaleFactor[i] == 0.0) {
1161       if (i==0) {
1162         fTPCScaleFactor[i] = 1.00;
1163       } else {
1164         for (int j=i-1;j>=0;j--) {
1165           if (fTPCScaleFactor[j] != 0.0) {
1166             fTPCScaleFactor[i] = fTPCScaleFactor[j];
1167             break;
1168           }
1169         }
1170       }
1171     }
1172   }
1173
1174
1175   //    for (int i=0; i<=(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << fV0MScaleFactor[i] << " for Run " << i+fLowRunN << endl;
1176
1177   return;
1178
1179 }
1180
1181
1182 //________________________________________________________________________
1183 void AliCentralitySelectionTask::MyInitScaleFactorMC()
1184 {
1185 // Initialize the MC scaling factors
1186   for (int i=0; i<(fHighRunN-fLowRunN); i++) fV0MScaleFactorMC[i] = 0.0;
1187   // scale factors determined from <V0 charge> on a run-by-run basis
1188   switch (fPass) {
1189   case 1:
1190     fV0MScaleFactorMC[0] = 0.75108;
1191     break;
1192   case 2:
1193     fV0MScaleFactorMC[0] = 0.8;
1194     break;
1195   default:
1196     AliError(Form("Unknown reconstruction pass (%d)! Setting MC scale in V0 to 1", fPass));
1197     fV0MScaleFactorMC[0] = 1.0;
1198     break;
1199   }
1200
1201   // set all missing values to the value of the run before it ....
1202   for (int i=0; i<=(fHighRunN-fLowRunN); i++) {
1203     if (fV0MScaleFactorMC[i] == 0.0) {
1204       if (i==0) {
1205         fV0MScaleFactorMC[i] = 1.00;
1206       } else {
1207         // search for last run number with non-zero value
1208         for (int j=i-1;j>=0;j--) {
1209           if (fV0MScaleFactorMC[j] != 0.0) {
1210             fV0MScaleFactorMC[i] = fV0MScaleFactorMC[j];
1211             break;
1212           }
1213         }
1214       }
1215     }
1216   } // end loop over checking all run-numbers
1217
1218
1219   //    for (int i=0; i<=(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << fV0MScaleFactorMC[i] << " for Run " << i+fLowRunN << endl;
1220
1221   return;
1222
1223 }
1224