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 "AliVEvent.h"
  41 #include "AliESD.h"
  42 #include "AliESDEvent.h"
  43 #include "AliESDHeader.h"
  44 #include "AliESDInputHandler.h"
  45 #include "AliESDZDC.h"
  46 #include "AliESDFMD.h"
  47 #include "AliESDVZERO.h"
  48 #include "AliCentrality.h"
  49 #include "AliESDtrackCuts.h"
  50 #include "AliMultiplicity.h"
  51 #include "AliAODHandler.h"
  52 #include "AliAODEvent.h"
  53 #include "AliESDVertex.h"
  54 #include "AliAODVertex.h"
  55 #include "AliAODMCHeader.h"
  56 #include "AliMCEvent.h"
  57 #include "AliMCEventHandler.h"
  58 #include "AliMCParticle.h"
  59 #include "AliStack.h"
  60 #include "AliHeader.h"
  61 #include "AliAODMCParticle.h"
  62 #include "AliAnalysisTaskSE.h"
  63 #include "AliGenEventHeader.h"
  64 #include "AliGenHijingEventHeader.h"
  65 #include "AliPhysicsSelectionTask.h"
  66 #include "AliPhysicsSelection.h"
  67 #include "AliBackgroundSelection.h"
  68 #include "AliESDUtils.h"
  69
  70 ClassImp(AliCentralitySelectionTask)
  71
  72
  73 //________________________________________________________________________
  74 AliCentralitySelectionTask::AliCentralitySelectionTask():
  75 AliAnalysisTaskSE(),
  76   fDebug(0),
  77   fAnalysisInput("ESD"),
  78   fIsMCInput(kFALSE),
  79   fFile(0),
  80   fFile2(0),
  81   fCurrentRun(-1),
  82   fRunNo(-1),
  83   fLowRunN(0),
  84   fHighRunN(0),
  85   fUseScaling(0),
  86   fUseCleaning(0),
  87   fTrackCuts(0),
  88   fZVCut(10),
  89   fOutliersCut(5),
  90   fQuality(0),
  91   fCentV0M(0),
  92   fCentFMD(0),
  93   fCentTRK(0),
  94   fCentTKL(0),
  95   fCentCL0(0),
  96   fCentCL1(0),
  97   fCentV0MvsFMD(0),
  98   fCentTKLvsV0M(0),
  99   fCentZEMvsZDC(0),
 100   fHtempV0M(0),
 101   fHtempFMD(0),
 102   fHtempTRK(0),
 103   fHtempTKL(0),
 104   fHtempCL0(0),
 105   fHtempCL1(0),
 106   fHtempV0MvsFMD(0),
 107   fHtempTKLvsV0M(0),
 108   fHtempZEMvsZDC(0),
 109   fOutputList(0),
 110   fHOutCentV0M     (0),
 111   fHOutCentFMD     (0),
 112   fHOutCentTRK     (0),
 113   fHOutCentTKL     (0),
 114   fHOutCentCL0     (0),
 115   fHOutCentCL1     (0),
 116   fHOutCentV0MvsFMD(0),
 117   fHOutCentTKLvsV0M(0),
 118   fHOutCentZEMvsZDC(0),
 119   fHOutCentV0MvsCentCL1(0),
 120   fHOutCentV0MvsCentTRK(0),
 121   fHOutCentTRKvsCentCL1(0),
 122   fHOutMultV0M(0),
 123   fHOutMultV0R(0),
 124   fHOutMultFMD(0),
 125   fHOutMultTRK(0),
 126   fHOutMultTKL(0),
 127   fHOutMultCL0(0),
 128   fHOutMultCL1(0),
 129   fHOutMultV0MvsZDN(0),
 130   fHOutMultZEMvsZDN(0),
 131   fHOutMultV0MvsZDC(0),
 132   fHOutMultZEMvsZDC(0),
 133   fHOutMultV0MvsCL1(0),
 134   fHOutMultV0MvsTRK(0),
 135   fHOutMultTRKvsCL1(0),
 136   fHOutCentV0M_qual1(0),
 137   fHOutCentTRK_qual1(0),
 138   fHOutCentCL1_qual1(0),
 139   fHOutCentV0M_qual2(0),
 140   fHOutCentTRK_qual2(0),
 141   fHOutCentCL1_qual2(0),
 142   fHOutQuality(0),
 143   fHOutVertex(0)
 144 {
 145   // Default constructor
 146   AliInfo("Centrality Selection enabled.");
 147   fLowRunN =136851;
 148   fHighRunN=139517;
 149
 150   for (Int_t i=0; i < 2667; i++) {
 151     V0MScaleFactor[i]=0.0;
 152     SPDScaleFactor[i]=0.0;
 153     TPCScaleFactor[i]=0.0;
 154     V0MScaleFactorMC[i]=0.0;
 155   }
 156   fUseScaling=kTRUE;
 157   fUseCleaning=kTRUE;
 158 }
 159
 160 //________________________________________________________________________
 161 AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
 162   AliAnalysisTaskSE(name),
 163   fDebug(0),
 164   fAnalysisInput("ESD"),
 165   fIsMCInput(kFALSE),
 166   fFile(0),
 167   fFile2(0),
 168   fCurrentRun(-1),
 169   fRunNo(-1),
 170   fLowRunN(0),
 171   fHighRunN(0),
 172   fUseScaling(0),
 173   fUseCleaning(0),
 174   fTrackCuts(0),
 175   fZVCut(10),
 176   fOutliersCut(5),
 177   fQuality(0),
 178   fCentV0M(0),
 179   fCentFMD(0),
 180   fCentTRK(0),
 181   fCentTKL(0),
 182   fCentCL0(0),
 183   fCentCL1(0),
 184   fCentV0MvsFMD(0),
 185   fCentTKLvsV0M(0),
 186   fCentZEMvsZDC(0),
 187   fHtempV0M(0),
 188   fHtempFMD(0),
 189   fHtempTRK(0),
 190   fHtempTKL(0),
 191   fHtempCL0(0),
 192   fHtempCL1(0),
 193   fHtempV0MvsFMD(0),
 194   fHtempTKLvsV0M(0),
 195   fHtempZEMvsZDC(0),
 196   fOutputList(0),
 197   fHOutCentV0M     (0),
 198   fHOutCentFMD     (0),
 199   fHOutCentTRK     (0),
 200   fHOutCentTKL     (0),
 201   fHOutCentCL0     (0),
 202   fHOutCentCL1     (0),
 203   fHOutCentV0MvsFMD(0),
 204   fHOutCentTKLvsV0M(0),
 205   fHOutCentZEMvsZDC(0),
 206   fHOutCentV0MvsCentCL1(0),
 207   fHOutCentV0MvsCentTRK(0),
 208   fHOutCentTRKvsCentCL1(0),
 209   fHOutMultV0M(0),
 210   fHOutMultV0R(0),
 211   fHOutMultFMD(0),
 212   fHOutMultTRK(0),
 213   fHOutMultTKL(0),
 214   fHOutMultCL0(0),
 215   fHOutMultCL1(0),
 216   fHOutMultV0MvsZDN(0),
 217   fHOutMultZEMvsZDN(0),
 218   fHOutMultV0MvsZDC(0),
 219   fHOutMultZEMvsZDC(0),
 220   fHOutMultV0MvsCL1(0),
 221   fHOutMultV0MvsTRK(0),
 222   fHOutMultTRKvsCL1(0),
 223   fHOutCentV0M_qual1(0),
 224   fHOutCentTRK_qual1(0),
 225   fHOutCentCL1_qual1(0),
 226   fHOutCentV0M_qual2(0),
 227   fHOutCentTRK_qual2(0),
 228   fHOutCentCL1_qual2(0),
 229   fHOutQuality(0),
 230   fHOutVertex(0)
 231 {
 232     // Default constructor
 233     fLowRunN =136851;
 234     fHighRunN=139517;
 235
 236   AliInfo("Centrality Selection enabled.");
 237   DefineOutput(1, TList::Class());
 238   for (Int_t i=0; i<2667; i++) {
 239     V0MScaleFactor[i]=0.0;
 240     SPDScaleFactor[i]=0.0;
 241     TPCScaleFactor[i]=0.0;
 242     V0MScaleFactorMC[i]=0.0;
 243   }
 244   fUseScaling=kTRUE;
 245   fUseCleaning=kTRUE;
 246 }
 247
 248 //________________________________________________________________________
 249 AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
 250 {
 251   // Assignment operator
 252   if (this!=&c) {
 253     AliAnalysisTaskSE::operator=(c);
 254   }
 255   return *this;
 256 }
 257
 258 //________________________________________________________________________
 259 AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
 260   AliAnalysisTaskSE(ana),
 261   fDebug(ana.fDebug),
 262   fAnalysisInput(ana.fDebug),
 263   fIsMCInput(ana.fIsMCInput),
 264   fFile(ana.fFile),
 265   fFile2(ana.fFile2),
 266   fCurrentRun(ana.fCurrentRun),
 267   fRunNo(ana.fRunNo),
 268   fLowRunN(ana.fLowRunN),
 269   fHighRunN(ana.fHighRunN),
 270   fUseScaling(ana.fUseScaling),
 271   fUseCleaning(ana.fUseCleaning),
 272   fTrackCuts(ana.fTrackCuts),
 273   fZVCut(ana.fZVCut),
 274   fOutliersCut(ana.fOutliersCut),
 275   fQuality(ana.fQuality),
 276   fCentV0M(ana.fCentV0M),
 277   fCentFMD(ana.fCentFMD),
 278   fCentTRK(ana.fCentTRK),
 279   fCentTKL(ana.fCentTKL),
 280   fCentCL0(ana.fCentCL0),
 281   fCentCL1(ana.fCentCL1),
 282   fCentV0MvsFMD(ana.fCentV0MvsFMD),
 283   fCentTKLvsV0M(ana.fCentTKLvsV0M),
 284   fCentZEMvsZDC(ana.fCentZEMvsZDC),
 285   fHtempV0M(ana.fHtempV0M),
 286   fHtempFMD(ana.fHtempFMD),
 287   fHtempTRK(ana.fHtempTRK),
 288   fHtempTKL(ana.fHtempTKL),
 289   fHtempCL0(ana.fHtempCL0),
 290   fHtempCL1(ana.fHtempCL1),
 291   fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
 292   fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
 293   fHtempZEMvsZDC(ana.fHtempZEMvsZDC),
 294   fOutputList(ana.fOutputList),
 295   fHOutCentV0M     (ana.fHOutCentV0M     ),
 296   fHOutCentFMD     (ana.fHOutCentFMD     ),
 297   fHOutCentTRK     (ana.fHOutCentTRK     ),
 298   fHOutCentTKL     (ana.fHOutCentTKL     ),
 299   fHOutCentCL0     (ana.fHOutCentCL0     ),
 300   fHOutCentCL1     (ana.fHOutCentCL1     ),
 301   fHOutCentV0MvsFMD(ana.fHOutCentV0MvsFMD),
 302   fHOutCentTKLvsV0M(ana.fHOutCentTKLvsV0M),
 303   fHOutCentZEMvsZDC(ana.fHOutCentZEMvsZDC),
 304   fHOutCentV0MvsCentCL1(ana.fHOutCentV0MvsCentCL1),
 305   fHOutCentV0MvsCentTRK(ana.fHOutCentV0MvsCentTRK),
 306   fHOutCentTRKvsCentCL1(ana.fHOutCentTRKvsCentCL1),
 307   fHOutMultV0M(ana.fHOutMultV0M),
 308   fHOutMultV0R(ana.fHOutMultV0R),
 309   fHOutMultFMD(ana.fHOutMultFMD),
 310   fHOutMultTRK(ana.fHOutMultTRK),
 311   fHOutMultTKL(ana.fHOutMultTKL),
 312   fHOutMultCL0(ana.fHOutMultCL0),
 313   fHOutMultCL1(ana.fHOutMultCL1),
 314   fHOutMultV0MvsZDN(ana.fHOutMultV0MvsZDN),
 315   fHOutMultZEMvsZDN(ana.fHOutMultZEMvsZDN),
 316   fHOutMultV0MvsZDC(ana.fHOutMultV0MvsZDC),
 317   fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC),
 318   fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1),
 319   fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK),
 320   fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1),
 321   fHOutCentV0M_qual1(ana.fHOutCentV0M_qual1),
 322   fHOutCentTRK_qual1(ana.fHOutCentTRK_qual1),
 323   fHOutCentCL1_qual1(ana.fHOutCentCL1_qual1),
 324   fHOutCentV0M_qual2(ana.fHOutCentV0M_qual2),
 325   fHOutCentTRK_qual2(ana.fHOutCentTRK_qual2),
 326   fHOutCentCL1_qual2(ana.fHOutCentCL1_qual2),
 327   fHOutQuality(ana.fHOutQuality),
 328   fHOutVertex(ana.fHOutVertex)
 329 {
 330   // Copy Constructor
 331     for (Int_t i=0; i<2667; i++) {
 332         V0MScaleFactor[i]=0.0;
 333         SPDScaleFactor[i]=0.0;
 334         TPCScaleFactor[i]=0.0;
 335         V0MScaleFactorMC[i]=0.0;
 336     }
 337
 338 }
 339
 340 //________________________________________________________________________
 341 AliCentralitySelectionTask::~AliCentralitySelectionTask()
 342 {
 343   // Destructor
 344   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;
 345   if (fTrackCuts) delete fTrackCuts;
 346 }
 347
 348 //________________________________________________________________________
 349 void AliCentralitySelectionTask::UserCreateOutputObjects()
 350 {
 351   // Create the output containers
 352   if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
 353   AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);
 354
 355   fOutputList = new TList();
 356   fOutputList->SetOwner();
 357   fHOutCentV0M     = new TH1F("fHOutCentV0M","fHOutCentV0M; Centrality V0",501,0,101);
 358   fHOutCentFMD     = new TH1F("fHOutCentFMD","fHOutCentFMD; Centrality FMD",501,0,101);
 359   fHOutCentTRK     = new TH1F("fHOutCentTRK","fHOutCentTRK; Centrality TPC",501,0,101);
 360   fHOutCentTKL     = new TH1F("fHOutCentTKL","fHOutCentTKL; Centrality tracklets",501,0,101);
 361   fHOutCentCL0     = new TH1F("fHOutCentCL0","fHOutCentCL0; Centrality SPD inner",501,0,101);
 362   fHOutCentCL1     = new TH1F("fHOutCentCL1","fHOutCentCL1; Centrality SPD outer",501,0,101);
 363   fHOutCentV0MvsFMD= new TH1F("fHOutCentV0MvsFMD","fHOutCentV0MvsFMD; Centrality V0 vs FMD",501,0,101);
 364   fHOutCentTKLvsV0M= new TH1F("fHOutCentTKLvsV0M","fHOutCentTKLvsV0M; Centrality tracklets vs V0",501,0,101);
 365   fHOutCentZEMvsZDC= new TH1F("fHOutCentZEMvsZDC","fHOutCentZEMvsZDC; Centrality ZEM vs ZDC",501,0,101);
 366   fHOutCentV0MvsCentCL1= new TH2F("fHOutCentV0MvsCentCL1","fHOutCentV0MvsCentCL1; Cent V0 vs Cent SPD",501,0,101,501,0,101);
 367   fHOutCentV0MvsCentTRK= new TH2F("fHOutCentV0MvsCentTRK","fHOutCentV0MvsCentTRK; Cent V0 vs Cent TPC",501,0,101,501,0,101);
 368   fHOutCentTRKvsCentCL1= new TH2F("fHOutCentTRKvsCentCL1","fHOutCentTRKvsCentCL1; Cent TPC vs Cent SPD",501,0,101,501,0,101);
 369
 370   fHOutMultV0M = new TH1F("fHOutMultV0M","fHOutMultV0M; Multiplicity V0",25000,0,25000);
 371   fHOutMultV0R = new TH1F("fHOutMultV0R","fHOutMultV0R; Multiplicity V0",30000,0,30000);
 372   fHOutMultFMD = new TH1F("fHOutMultFMD","fHOutMultFMD; Multiplicity FMD",24000,0,24000);
 373   fHOutMultTRK = new TH1F("fHOutMultTRK","fHOutMultTRK; Multiplicity TPC",4000,0,4000);
 374   fHOutMultTKL = new TH1F("fHOutMultTKL","fHOutMultTKL; Multiplicity tracklets",5000,0,5000);
 375   fHOutMultCL0 = new TH1F("fHOutMultCL0","fHOutMultCL0; Multiplicity SPD inner",7000,0,7000);
 376   fHOutMultCL1 = new TH1F("fHOutMultCL1","fHOutMultCL1; Multiplicity SPD outer",7000,0,7000);
 377   fHOutMultV0MvsZDN = new TH2F("fHOutMultV0MvsZDN","fHOutMultV0MvsZDN; Multiplicity V0; Energy ZDC-N",500,0,25000,500,0,180000);
 378   fHOutMultZEMvsZDN = new TH2F("fHOutMultZEMvsZDN","fHOutMultZEMvsZDN; Energy ZEM; Energy ZDC-N",500,0,2500,500,0,180000);
 379   fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,25000,500,0,200000);
 380   fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,200000);
 381   fHOutMultV0MvsCL1 = new TH2F("fHOutMultV0MvsCL1","fHOutMultV0MvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);
 382   fHOutMultV0MvsTRK = new TH2F("fHOutMultV0MvsTRK","fHOutMultV0MvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);
 383   fHOutMultTRKvsCL1 = new TH2F("fHOutMultTRKvsCL1","fHOutMultTRKvsCL1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);
 384
 385   fHOutCentV0M_qual1 = new TH1F("fHOutCentV0M_qual1","fHOutCentV0M_qual1; Centrality V0",501,0,100);
 386   fHOutCentTRK_qual1 = new TH1F("fHOutCentTRK_qual1","fHOutCentTRK_qual1; Centrality TPC",501,0,100);
 387   fHOutCentCL1_qual1 = new TH1F("fHOutCentCL1_qual1","fHOutCentCL1_qual1; Centrality SPD outer",501,0,100);
 388
 389   fHOutCentV0M_qual2 = new TH1F("fHOutCentV0M_qual2","fHOutCentV0M_qual2; Centrality V0",101,-0.1,100.1);
 390   fHOutCentTRK_qual2 = new TH1F("fHOutCentTRK_qual2","fHOutCentTRK_qual2; Centrality TPC",101,-0.1,100.1);
 391   fHOutCentCL1_qual2 = new TH1F("fHOutCentCL1_qual2","fHOutCentCL1_qual2; Centrality SPD outer",101,-0.1,100.1);
 392
 393   fHOutQuality = new TH1F("fHOutQuality", "fHOutQuality", 10,-0.5,9.5);
 394   fHOutVertex  = new TH1F("fHOutVertex", "fHOutVertex", 100,-20,20);
 395
 396   fOutputList->Add(  fHOutCentV0M     );
 397   fOutputList->Add(  fHOutCentFMD     );
 398   fOutputList->Add(  fHOutCentTRK     );
 399   fOutputList->Add(  fHOutCentTKL     );
 400   fOutputList->Add(  fHOutCentCL0     );
 401   fOutputList->Add(  fHOutCentCL1     );
 402   fOutputList->Add(  fHOutCentV0MvsFMD);
 403   fOutputList->Add(  fHOutCentTKLvsV0M);
 404   fOutputList->Add(  fHOutCentZEMvsZDC);
 405   fOutputList->Add(  fHOutCentV0MvsCentCL1);
 406   fOutputList->Add(  fHOutCentV0MvsCentTRK);
 407   fOutputList->Add(  fHOutCentTRKvsCentCL1);
 408   fOutputList->Add(  fHOutMultV0M);
 409   fOutputList->Add(  fHOutMultV0R);
 410   fOutputList->Add(  fHOutMultFMD);
 411   fOutputList->Add(  fHOutMultTRK);
 412   fOutputList->Add(  fHOutMultTKL);
 413   fOutputList->Add(  fHOutMultCL0);
 414   fOutputList->Add(  fHOutMultCL1);
 415   fOutputList->Add(  fHOutMultV0MvsZDN);
 416   fOutputList->Add(  fHOutMultZEMvsZDN);
 417   fOutputList->Add(  fHOutMultV0MvsZDC);
 418   fOutputList->Add(  fHOutMultZEMvsZDC);
 419   fOutputList->Add(  fHOutMultV0MvsCL1);
 420   fOutputList->Add(  fHOutMultV0MvsTRK);
 421   fOutputList->Add(  fHOutMultTRKvsCL1);
 422   fOutputList->Add(  fHOutCentV0M_qual1 );
 423   fOutputList->Add(  fHOutCentTRK_qual1 );
 424   fOutputList->Add(  fHOutCentCL1_qual1 );
 425   fOutputList->Add(  fHOutCentV0M_qual2 );
 426   fOutputList->Add(  fHOutCentTRK_qual2 );
 427   fOutputList->Add(  fHOutCentCL1_qual2 );
 428   fOutputList->Add(  fHOutQuality );
 429   fOutputList->Add(  fHOutVertex );
 430
 431
 432   fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 433
 434   PostData(1, fOutputList);
 435
 436   MyInitScaleFactor();
 437   if (fIsMCInput) MyInitScaleFactorMC();
 438 }
 439
 440 //________________________________________________________________________
 441 void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
 442 {
 443   // Execute analysis for current event:
 444   if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
 445
 446   Float_t  zncEnergy = 0.;          //  ZNC Energy
 447   Float_t  zpcEnergy = 0.;          //  ZPC Energy
 448   Float_t  znaEnergy = 0.;          //  ZNA Energy
 449   Float_t  zpaEnergy = 0.;          //  ZPA Energy
 450   Float_t  zem1Energy = 0.;         //  ZEM1 Energy
 451   Float_t  zem2Energy = 0.;         //  ZEM2 Energy
 452   Bool_t   zdcEnergyCal = kFALSE;   // if zdc is calibrated (in pass2)
 453
 454   Int_t    nTracks = 0;             //  no. tracks
 455   Int_t    nTracklets = 0;          //  no. tracklets
 456   Int_t    nClusters[6] = {0};      //  no. clusters on 6 ITS layers
 457   Int_t    nChips[2];               //  no. chips on 2 SPD layers
 458   Float_t  spdCorr =0;              //  corrected spd2 multiplicity
 459
 460   Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
 461   Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
 462   Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A
 463   Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C
 464
 465   Short_t v0Corr = 0;               // corrected V0 multiplicity
 466   Short_t v0CorrResc = 0;           // corrected and rescaled V0 multiplicity
 467
 468   Float_t zvtx =0;                  // z-vertex SPD
 469
 470   AliCentrality *esdCent = 0;
 471
 472   if(fAnalysisInput.CompareTo("ESD")==0){
 473
 474     AliVEvent* event = InputEvent();
 475     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 476     if (!esd) {
 477         AliError("No ESD Event");
 478         return;
 479     }
 480
 481     if (fRunNo<=0) {
 482       if (SetupRun(esd)<0)
 483          AliFatal("Centrality File not available for this run");
 484     }
 485
 486     esdCent = esd->GetCentrality();
 487
 488     // ***** V0 info
 489     AliESDVZERO* esdV0 = esd->GetVZEROData();
 490     multV0A=esdV0->GetMTotV0A();
 491     multV0C=esdV0->GetMTotV0C();
 492
 493     Float_t v0CorrR;
 494     v0Corr = (Short_t)AliESDUtils::GetCorrV0(esd,v0CorrR);
 495     v0CorrResc = (Short_t)v0CorrR;
 496
 497     // ***** Vertex Info
 498     const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
 499     zvtx        = vtxESD->GetZ();
 500
 501     // ***** CB info (tracklets, clusters, chips)
 502     //nTracks    = event->GetNumberOfTracks();
 503     nTracks    = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
 504
 505     const AliMultiplicity *mult = esd->GetMultiplicity();
 506
 507     nTracklets = mult->GetNumberOfTracklets();
 508
 509     for(Int_t ilay=0; ilay<6; ilay++){
 510       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
 511     }
 512
 513     for(Int_t ilay=0; ilay<2; ilay++){
 514       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
 515     }
 516
 517     spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);
 518
 519     // ***** FMD info
 520     AliESDFMD *fmd = esd->GetFMDData();
 521     Float_t totalMultA = 0;
 522     Float_t totalMultC = 0;
 523     const Float_t fFMDLowCut = 0.4;
 524
 525     for(UShort_t det=1;det<=3;det++) {
 526       Int_t nRings = (det==1 ? 1 : 2);
 527       for (UShort_t ir = 0; ir < nRings; ir++) {
 528         Char_t   ring = (ir == 0 ? 'I' : 'O');
 529         UShort_t nsec = (ir == 0 ? 20  : 40);
 530         UShort_t nstr = (ir == 0 ? 512 : 256);
 531         for(UShort_t sec =0; sec < nsec;  sec++)  {
 532           for(UShort_t strip = 0; strip < nstr; strip++) {
 533
 534             Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
 535             if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
 536
 537             Float_t nParticles=0;
 538
 539             if(FMDmult > fFMDLowCut) {
 540               nParticles = 1.;
 541             }
 542
 543             if (det<3) totalMultA = totalMultA + nParticles;
 544             else totalMultC = totalMultC + nParticles;
 545
 546           }
 547         }
 548       }
 549     }
 550     multFMDA = totalMultA;
 551     multFMDC = totalMultC;
 552
 553     // ***** ZDC info
 554     AliESDZDC *esdZDC = esd->GetESDZDC();
 555     zdcEnergyCal = esdZDC->AliESDZDC::TestBit(AliESDZDC::kEnergyCalibratedSignal);
 556     if (zdcEnergyCal) {
 557       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
 558       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
 559       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
 560       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
 561     } else {
 562       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;
 563       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;
 564       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;
 565       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;
 566     }
 567     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;
 568     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;
 569
 570   }
 571   else if(fAnalysisInput.CompareTo("AOD")==0){
 572     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
 573     // to be implemented
 574     printf("  AOD analysis not yet implemented!!!\n\n");
 575     return;
 576   }
 577
 578
 579   // ***** Scaling
 580   // ***** Scaling for MC
 581   if (fIsMCInput) {
 582     fUseScaling=kFALSE;
 583     Float_t temp_scalefactorV0M = MyGetScaleFactorMC(fCurrentRun);
 584     v0Corr  = Short_t((multV0A+multV0C)  * temp_scalefactorV0M);
 585   }
 586   // ***** Scaling for Data
 587   if (fUseScaling) {
 588     Float_t temp_scalefactorV0M = MyGetScaleFactor(fCurrentRun,0);
 589     Float_t temp_scalefactorSPD = MyGetScaleFactor(fCurrentRun,1);
 590     Float_t temp_scalefactorTPC = MyGetScaleFactor(fCurrentRun,2);
 591     v0Corr  = Short_t(v0Corr / temp_scalefactorV0M);
 592     spdCorr = spdCorr / temp_scalefactorSPD;
 593     nTracks = Int_t(nTracks / temp_scalefactorTPC);
 594   }
 595
 596   // ***** Centrality Selection
 597   if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
 598   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
 599   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
 600   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
 601   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
 602   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
 603
 604   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
 605   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
 606   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin(zem1Energy+zem2Energy,zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
 607
 608   // ***** Cleaning
 609   if (fUseCleaning) {
 610       fQuality=0;
 611       fZVCut=10;
 612       fOutliersCut=6;
 613
 614       // ***** vertex
 615       if (TMath::Abs(zvtx)>fZVCut) fQuality += 1;
 616
 617       // ***** outliers
 618       // **** V0 vs SPD
 619       if (IsOutlierV0MSPD(spdCorr, v0Corr, int(fCentV0M))) fQuality  += 2;
 620       // ***** V0 vs TPC
 621       if (IsOutlierV0MTPC(nTracks, v0Corr, int(fCentV0M))) fQuality  += 4;
 622       // ***** V0 vs ZDC
 623        if (IsOutlierV0MZDC((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
 624            (zdcEnergyCal==kFALSE) && !(fIsMCInput)) fQuality  += 8;
 625        if (IsOutlierV0MZDCECal((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
 626            ((zdcEnergyCal==kTRUE) || (fIsMCInput))) fQuality  += 8;
 627   } else {
 628       fQuality = 0;
 629   }
 630
 631
 632   if (esdCent) {
 633       esdCent->SetQuality(fQuality);
 634       esdCent->SetCentralityV0M(fCentV0M);
 635       esdCent->SetCentralityFMD(fCentFMD);
 636       esdCent->SetCentralityTRK(fCentTRK);
 637       esdCent->SetCentralityTKL(fCentTKL);
 638       esdCent->SetCentralityCL0(fCentCL0);
 639       esdCent->SetCentralityCL1(fCentCL1);
 640       esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
 641       esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
 642       esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
 643   }
 644
 645   fHOutQuality->Fill(fQuality);
 646   fHOutVertex->Fill(zvtx);
 647
 648   fHOutMultV0MvsZDN->Fill(v0Corr,(zncEnergy+znaEnergy));
 649   fHOutMultZEMvsZDN->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy));
 650   fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
 651   fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
 652   fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
 653   fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
 654
 655   if (fQuality==0) {
 656     fHOutCentV0M->Fill(fCentV0M);
 657     fHOutCentFMD->Fill(fCentFMD);
 658     fHOutCentTRK->Fill(fCentTRK);
 659     fHOutCentTKL->Fill(fCentTKL);
 660     fHOutCentCL0->Fill(fCentCL0);
 661     fHOutCentCL1->Fill(fCentCL1);
 662     fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);
 663     fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);
 664     fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);
 665     fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentCL1);
 666     fHOutCentV0MvsCentTRK->Fill(fCentV0M,fCentTRK);
 667     fHOutCentTRKvsCentCL1->Fill(fCentTRK,fCentCL1);
 668     fHOutMultV0M->Fill(v0Corr);
 669     fHOutMultV0R->Fill(multV0A+multV0C);
 670     fHOutMultFMD->Fill((multFMDA+multFMDC));
 671     fHOutMultTRK->Fill(nTracks);
 672     fHOutMultTKL->Fill(nTracklets);
 673     fHOutMultCL0->Fill(nClusters[0]);
 674     fHOutMultCL1->Fill(spdCorr);
 675     fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
 676   } else if (fQuality ==1) {
 677     fHOutCentV0M_qual1->Fill(fCentV0M);
 678     fHOutCentTRK_qual1->Fill(fCentTRK);
 679     fHOutCentCL1_qual1->Fill(fCentCL1);
 680   } else {
 681     fHOutCentV0M_qual2->Fill(fCentV0M);
 682     fHOutCentTRK_qual2->Fill(fCentTRK);
 683     fHOutCentCL1_qual2->Fill(fCentCL1);
 684   }
 685
 686   PostData(1, fOutputList);
 687 }
 688
 689 //________________________________________________________________________
 690 void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename)
 691 {
 692   //  Read centrality histograms
 693   TDirectory *owd = gDirectory;
 694   // Check if the file is present
 695   TString path = gSystem->ExpandPathName(fCentfilename.Data());
 696   if (gSystem->AccessPathName(path)) {
 697      AliError(Form("File %s does not exist", path.Data()));
 698      return;
 699   }
 700   fFile  = TFile::Open(fCentfilename);
 701   owd->cd();
 702   fHtempV0M  = (TH1F*) (fFile->Get("hmultV0_percentile"));
 703   fHtempFMD  = (TH1F*) (fFile->Get("hmultFMD_percentile"));
 704   fHtempTRK  = (TH1F*) (fFile->Get("hNtracks_percentile"));
 705   fHtempTKL  = (TH1F*) (fFile->Get("hNtracklets_percentile"));
 706   fHtempCL0  = (TH1F*) (fFile->Get("hNclusters0_percentile"));
 707   fHtempCL1  = (TH1F*) (fFile->Get("hNclusters1_percentile"));
 708
 709   if (!fHtempV0M) AliWarning(Form("Calibration for V0M does not exist in %s", path.Data()));
 710   if (!fHtempFMD) AliWarning(Form("Calibration for FMD does not exist in %s", path.Data()));
 711   if (!fHtempTRK) AliWarning(Form("Calibration for TRK does not exist in %s", path.Data()));
 712   if (!fHtempTKL) AliWarning(Form("Calibration for TKL does not exist in %s", path.Data()));
 713   if (!fHtempCL0) AliWarning(Form("Calibration for CL0 does not exist in %s", path.Data()));
 714   if (!fHtempCL1) AliWarning(Form("Calibration for CL1 does not exist in %s", path.Data()));
 715
 716   owd->cd();
 717 }
 718
 719 //________________________________________________________________________
 720 void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2)
 721 {
 722   //  Read centrality histograms
 723   TDirectory *owd = gDirectory;
 724   TString path = gSystem->ExpandPathName(fCentfilename2.Data());
 725   if (gSystem->AccessPathName(path)) {
 726      AliError(Form("File %s does not exist", path.Data()));
 727      return;
 728   }
 729   fFile2  = TFile::Open(fCentfilename2);
 730   owd->cd();
 731   fHtempV0MvsFMD =  (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
 732   fHtempTKLvsV0M  = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
 733   fHtempZEMvsZDC  = (TH2F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
 734
 735   if (!fHtempV0MvsFMD) AliWarning(Form("Calibration for V0MvsFMD does not exist in %s", path.Data()));
 736   if (!fHtempTKLvsV0M) AliWarning(Form("Calibration for TKLvsV0M does not exist in %s", path.Data()));
 737   if (!fHtempZEMvsZDC) AliWarning(Form("Calibration for ZEMvsZDC does not exist in %s", path.Data()));
 738
 739   owd->cd();
 740 }
 741
 742 //________________________________________________________________________
 743 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
 744 {
 745   // Terminate analysis
 746   if (fFile && fFile->IsOpen())
 747     fFile->Close();
 748   if (fFile2 && fFile2->IsOpen())
 749     fFile2->Close();
 750 }
 751 //________________________________________________________________________
 752 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* const esd)
 753 {
 754   // Setup files for run
 755
 756   if (!esd)
 757     return -1;
 758
 759   // check if something to be done
 760   if (fCurrentRun == esd->GetRunNumber())
 761     return 0;
 762   else
 763     fCurrentRun = esd->GetRunNumber();
 764
 765   AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
 766
 767   // CHANGE HERE FOR RUN RANGES
 768   if ( fCurrentRun <= 137165 ) fRunNo = 137161;
 769   else fRunNo = 137366;
 770   // CHANGE HERE FOR RUN RANGES
 771
 772   TString fileName(Form("%s/COMMON/CENTRALITY/data/AliCentralityBy1D_%d.root", AliAnalysisManager::GetOADBPath(), fRunNo));
 773   TString fileName2(Form("%s/COMMON/CENTRALITY/data/AliCentralityByFunction_%d.root", AliAnalysisManager::GetOADBPath(), fRunNo));
 774
 775   AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data()));
 776   ReadCentralityHistos(fileName.Data());
 777   ReadCentralityHistos2(fileName2.Data());
 778   if (!fFile && !fFile2) {
 779      AliFatal(Form("Run %d not known to centrality selection!", fCurrentRun));
 780      return -1;
 781   }
 782   return 0;
 783 }
 784
 785 //________________________________________________________________________
 786 Bool_t AliCentralitySelectionTask::IsOutlierV0MSPD(Float_t spd, Float_t v0, Int_t cent) const
 787 {
 788 // Clean outliers
 789   Float_t val= -0.143789 + 0.288874 * v0;
 790   Float_t SPDsigma[101]={231.483, 189.446, 183.359, 179.923, 174.229, 170.309, 165.021,
 791                          160.84, 159.33, 154.453, 151.644, 148.337, 145.215, 142.353,
 792                          139.351, 136, 133.838, 129.885, 127.36, 125.032, 122.21, 120.3,
 793                          117.766, 114.77, 113.1, 110.268, 107.463, 105.293, 102.845,
 794                          100.835, 98.9632, 97.3287, 93.6887, 92.1066, 89.3224, 87.8382,
 795                          86.04, 83.6431, 81.9655, 80.0491, 77.8208, 76.4716, 74.2165,
 796                          72.2752, 70.4875, 68.9414, 66.8622, 65.022, 63.5134, 61.8228,
 797                          59.7166, 58.5008, 56.2789, 54.9615, 53.386, 51.2165, 49.4842,
 798                          48.259, 47.1129, 45.3115, 43.8486, 41.9207, 40.5754, 39.3872,
 799                          38.1897, 36.5401, 35.1283, 33.9702, 32.6429, 31.3612, 29.5876,
 800                          28.9319, 27.564, 26.0443, 25.2836, 23.9753, 22.8936, 21.5665,
 801                          20.7048, 19.8016, 18.7095, 18.1144, 17.2095, 16.602, 16.3233,
 802                          15.7185, 15.3006, 14.7432, 14.4174, 14.0805, 13.7638, 13.7638,
 803                          13.7638, 13.7638, 13.7638, 13.7638, 13.7638, 13.7638, 13.7638, 18.0803, 18.0803};
 804
 805   if ( TMath::Abs(spd-val) > fOutliersCut*SPDsigma[cent] )
 806     return kTRUE;
 807   else
 808     return kFALSE;
 809 }
 810
 811 //________________________________________________________________________
 812 Bool_t AliCentralitySelectionTask::IsOutlierV0MTPC(Int_t tracks, Float_t v0, Int_t cent) const
 813 {
 814 // Clean outliers
 815   Float_t val = -0.540691 + 0.128358 * v0;
 816   Float_t TPCsigma[101]={106.439, 89.2834, 86.7568, 85.3641, 83.379, 81.6093, 79.3189,
 817                          78.0616, 77.2167, 75.0021, 73.9957, 72.0926, 71.0442, 69.8395,
 818                          68.1169, 66.6676, 66.0038, 64.2284, 63.3845, 61.7439, 60.642,
 819                          59.5383, 58.3696, 57.0227, 56.0619, 54.7108, 53.8382, 52.3398,
 820                          51.5297, 49.9488, 49.1126, 48.208, 46.8566, 45.7724, 44.7829,
 821                          43.8726, 42.7499, 41.9307, 40.6874, 39.9619, 38.5534, 38.0884,
 822                          36.6141, 35.7482, 34.8804, 34.1769, 33.1278, 32.3435, 31.4783,
 823                          30.2587, 29.4741, 28.8575, 27.9298, 26.9752, 26.1675, 25.1234,
 824                          24.4702, 23.6843, 22.9764, 21.8579, 21.2924, 20.3241, 19.8296,
 825                          19.2465, 18.4474, 17.7216, 16.8956, 16.342, 15.626, 15.0329,
 826                          14.3911, 13.9301, 13.254, 12.6745, 12.2436, 11.7776, 11.1795,
 827                          10.673, 10.27, 9.95646, 9.50939, 9.26162, 8.95315, 8.73439,
 828                          8.67375, 8.43029, 8.34818, 8.33484, 8.40709, 8.3974, 8.32814,
 829                          8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 8.32814, 12.351, 12.351};
 830
 831   if ( TMath::Abs(tracks-val) > fOutliersCut*TPCsigma[cent] )
 832     return kTRUE;
 833   else
 834     return kFALSE;
 835 }
 836
 837 //________________________________________________________________________
 838 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDC(Float_t zdc, Float_t v0) const
 839 {
 840 // Clean outliers
 841   Float_t val1 = 6350. - 0.26 * v0;
 842   Float_t val2 = 5580.;
 843   if ((zdc >  val1) || (zdc > val2))
 844     return kTRUE;
 845   else
 846     return kFALSE;
 847 }
 848
 849 //________________________________________________________________________
 850 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDCECal(Float_t /*zdc*/, Float_t /*v0*/) const
 851 {
 852 // Clean outliers
 853     return kFALSE;
 854 }
 855
 856 //________________________________________________________________________
 857 Float_t AliCentralitySelectionTask::MyGetScaleFactor(Int_t runnumber, Int_t flag) const
 858 {
 859 // Get scaling factor
 860   if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
 861     cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
 862     return 0.0;
 863   }
 864
 865   Float_t scalefactor=0.0;
 866   if (flag==0)
 867     scalefactor = V0MScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
 868   else if (flag==1)
 869     scalefactor = SPDScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
 870   else if (flag==2)
 871     scalefactor = TPCScaleFactor[runnumber - fLowRunN]; // subtracting reference offset index
 872
 873   return scalefactor;
 874
 875 }
 876
 877 //________________________________________________________________________
 878 Float_t AliCentralitySelectionTask::MyGetScaleFactorMC(Int_t runnumber) const
 879 {
 880 // Get MC scaling factor
 881   if (! (runnumber >= fLowRunN && runnumber <=fHighRunN)) {
 882     cout << "MyGetScaleFactor error in run number range " << runnumber << endl;
 883     return 0.0;
 884   }
 885
 886   Float_t scalefactor= V0MScaleFactorMC[runnumber - fLowRunN]; // subtracting reference offset index
 887   return scalefactor;
 888
 889 }
 890
 891 //________________________________________________________________________
 892 void AliCentralitySelectionTask::MyInitScaleFactor ()
 893 {
 894 // Initialize the scaling factors
 895   for (int i=0; i<(fHighRunN-fLowRunN); i++) V0MScaleFactor[i] = 0.0;
 896   for (int i=0; i<(fHighRunN-fLowRunN); i++) SPDScaleFactor[i] = 0.0;
 897   for (int i=0; i<(fHighRunN-fLowRunN); i++) TPCScaleFactor[i] = 0.0;
 898
 899   // scale factors determined from <V0 charge> on a run-by-run basis
 900    V0MScaleFactor[310] = 1.;
 901    V0MScaleFactor[311] = 1.;
 902    V0MScaleFactor[514] = 1.0046;
 903    V0MScaleFactor[515] = 0.983535;
 904    V0MScaleFactor[579] = 0.988185;
 905    V0MScaleFactor[580] = 0.983351;
 906    V0MScaleFactor[581] = 0.989013;
 907    V0MScaleFactor[583] = 0.990056;
 908    V0MScaleFactor[588] = 0.974438;
 909    V0MScaleFactor[589] = 0.981572;
 910    V0MScaleFactor[590] = 0.989316;
 911    V0MScaleFactor[592] = 0.98345;
 912    V0MScaleFactor[688] = 0.993647;
 913    V0MScaleFactor[690] = 0.994758;
 914    V0MScaleFactor[693] = 0.989569;
 915    V0MScaleFactor[698] = 0.993119;
 916    V0MScaleFactor[744] = 0.989583;
 917    V0MScaleFactor[757] = 0.990377;
 918    V0MScaleFactor[787] = 0.990176;
 919    V0MScaleFactor[788] = 0.98723;
 920    V0MScaleFactor[834] = 1.00403;
 921    V0MScaleFactor[835] = 0.994376;
 922    V0MScaleFactor[840] = 0.99759;
 923    V0MScaleFactor[842] = 1.01031;
 924    V0MScaleFactor[900] = 0.996216;
 925    V0MScaleFactor[901] = 0.994205;
 926    V0MScaleFactor[992] = 0.998479;
 927    V0MScaleFactor[997] = 1.00078;
 928    V0MScaleFactor[1299] = 1.00515;
 929    V0MScaleFactor[1303] = 1.00094;
 930    V0MScaleFactor[1339] = 0.986596;
 931    V0MScaleFactor[1346] = 0.972226;
 932    V0MScaleFactor[1349] = 0.960358;
 933    V0MScaleFactor[1350] = 0.970023;
 934    V0MScaleFactor[1374] = 1.00575;
 935    V0MScaleFactor[1545] = 1.00471;
 936    V0MScaleFactor[1587] = 1.00611;
 937    V0MScaleFactor[1588] = 1.00976;
 938    V0MScaleFactor[1618] = 1.00771;
 939    V0MScaleFactor[1682] = 1.01622;
 940    V0MScaleFactor[1727] = 1.01305;
 941    V0MScaleFactor[1728] = 1.00388;
 942    V0MScaleFactor[1731] = 1.00673;
 943    V0MScaleFactor[1732] = 1.00916;
 944    V0MScaleFactor[1770] = 1.0092;
 945    V0MScaleFactor[1773] = 1.00728;
 946    V0MScaleFactor[1786] = 1.01655;
 947    V0MScaleFactor[1787] = 1.00672;
 948    V0MScaleFactor[1801] = 0.983339;
 949    V0MScaleFactor[1802] = 1.00754;
 950    V0MScaleFactor[1811] = 1.00608;
 951    V0MScaleFactor[1815] = 1.01227;
 952    V0MScaleFactor[1879] = 0.99907;
 953    V0MScaleFactor[1881] = 0.995696;
 954    V0MScaleFactor[2186] = 1.00559;
 955    V0MScaleFactor[2187] = 1.00631;
 956    V0MScaleFactor[2254] = 1.01512;
 957    V0MScaleFactor[2256] = 0.998727;
 958    V0MScaleFactor[2321] = 1.00701;
 959    SPDScaleFactor[310] = 1.00211;
 960    SPDScaleFactor[311] = 1.00067;
 961    SPDScaleFactor[514] = 1.02268;
 962    SPDScaleFactor[515] = 0.994902;
 963    SPDScaleFactor[579] = 1.00215;
 964    SPDScaleFactor[580] = 0.993421;
 965    SPDScaleFactor[581] = 1.00129;
 966    SPDScaleFactor[583] = 1.00242;
 967    SPDScaleFactor[588] = 0.984762;
 968    SPDScaleFactor[589] = 0.994355;
 969    SPDScaleFactor[590] = 1.00073;
 970    SPDScaleFactor[592] = 0.995889;
 971    SPDScaleFactor[688] = 0.994532;
 972    SPDScaleFactor[690] = 0.998307;
 973    SPDScaleFactor[693] = 0.994052;
 974    SPDScaleFactor[698] = 0.993224;
 975    SPDScaleFactor[744] = 0.993279;
 976    SPDScaleFactor[757] = 0.992494;
 977    SPDScaleFactor[787] = 0.992678;
 978    SPDScaleFactor[788] = 0.996563;
 979    SPDScaleFactor[834] = 1.01116;
 980    SPDScaleFactor[835] = 0.993108;
 981    SPDScaleFactor[840] = 0.997574;
 982    SPDScaleFactor[842] = 1.01829;
 983    SPDScaleFactor[900] = 0.999438;
 984    SPDScaleFactor[901] = 0.995849;
 985    SPDScaleFactor[992] = 0.999227;
 986    SPDScaleFactor[997] = 1.00575;
 987    SPDScaleFactor[1299] = 0.99877;
 988    SPDScaleFactor[1303] = 0.999682;
 989    SPDScaleFactor[1339] = 0.978198;
 990    SPDScaleFactor[1346] = 0.964178;
 991    SPDScaleFactor[1349] = 0.959439;
 992    SPDScaleFactor[1350] = 0.956945;
 993    SPDScaleFactor[1374] = 0.994434;
 994    SPDScaleFactor[1545] = 1.0016;
 995    SPDScaleFactor[1587] = 1.00153;
 996    SPDScaleFactor[1588] = 1.00698;
 997    SPDScaleFactor[1618] = 1.00554;
 998    SPDScaleFactor[1682] = 1.0123;
 999    SPDScaleFactor[1727] = 1.011;
1000    SPDScaleFactor[1728] = 1.00143;
1001    SPDScaleFactor[1731] = 1.00486;
1002    SPDScaleFactor[1732] = 1.00646;
1003    SPDScaleFactor[1770] = 1.00515;
1004    SPDScaleFactor[1773] = 1.00485;
1005    SPDScaleFactor[1786] = 1.01215;
1006    SPDScaleFactor[1787] = 1.00419;
1007    SPDScaleFactor[1801] = 0.983327;
1008    SPDScaleFactor[1802] = 1.00529;
1009    SPDScaleFactor[1811] = 1.00367;
1010    SPDScaleFactor[1815] = 1.01045;
1011    SPDScaleFactor[1879] = 0.996374;
1012    SPDScaleFactor[1881] = 0.988827;
1013    SPDScaleFactor[2186] = 1.00354;
1014    SPDScaleFactor[2187] = 1.00397;
1015    SPDScaleFactor[2254] = 1.01138;
1016    SPDScaleFactor[2256] = 0.996641;
1017    SPDScaleFactor[2321] = 1.00357;
1018    TPCScaleFactor[310] = 1.00434;
1019    TPCScaleFactor[311] = 1.0056;
1020    TPCScaleFactor[514] = 1.02185;
1021    TPCScaleFactor[515] = 1.0036;
1022    TPCScaleFactor[579] = 1.00607;
1023    TPCScaleFactor[580] = 1.00183;
1024    TPCScaleFactor[581] = 1.00693;
1025    TPCScaleFactor[583] = 1.00746;
1026    TPCScaleFactor[588] = 0.990524;
1027    TPCScaleFactor[589] = 0.998582;
1028    TPCScaleFactor[590] = 1.00618;
1029    TPCScaleFactor[592] = 1.00088;
1030    TPCScaleFactor[688] = 1.00598;
1031    TPCScaleFactor[690] = 1.00658;
1032    TPCScaleFactor[693] = 1.00144;
1033    TPCScaleFactor[698] = 1.00279;
1034    TPCScaleFactor[744] = 1.00122;
1035    TPCScaleFactor[757] = 1.002;
1036    TPCScaleFactor[787] = 0.997818;
1037    TPCScaleFactor[788] = 0.994583;
1038    TPCScaleFactor[834] = 1.01508;
1039    TPCScaleFactor[835] = 1.00218;
1040    TPCScaleFactor[840] = 1.00569;
1041    TPCScaleFactor[842] = 1.01789;
1042    TPCScaleFactor[900] = 1.00739;
1043    TPCScaleFactor[901] = 1.00462;
1044    TPCScaleFactor[992] = 1.00502;
1045    TPCScaleFactor[997] = 1.00943;
1046    TPCScaleFactor[1299] = 1.00438;
1047    TPCScaleFactor[1303] = 0.996701;
1048    TPCScaleFactor[1339] = 0.978641;
1049    TPCScaleFactor[1346] = 0.968906;
1050    TPCScaleFactor[1349] = 0.954311;
1051    TPCScaleFactor[1350] = 0.958764;
1052    TPCScaleFactor[1374] = 0.997899;
1053    TPCScaleFactor[1545] = 0.992;
1054    TPCScaleFactor[1587] = 0.992635;
1055    TPCScaleFactor[1588] = 1.00207;
1056    TPCScaleFactor[1618] = 1.00126;
1057    TPCScaleFactor[1682] = 1.00324;
1058    TPCScaleFactor[1727] = 1.00042;
1059    TPCScaleFactor[1728] = 0.978881;
1060    TPCScaleFactor[1731] = 0.999818;
1061    TPCScaleFactor[1732] = 1.00109;
1062    TPCScaleFactor[1770] = 0.99935;
1063    TPCScaleFactor[1773] = 0.998531;
1064    TPCScaleFactor[1786] = 0.999125;
1065    TPCScaleFactor[1787] = 0.998479;
1066    TPCScaleFactor[1801] = 0.9775;
1067    TPCScaleFactor[1802] = 0.999095;
1068    TPCScaleFactor[1811] = 0.998197;
1069    TPCScaleFactor[1815] = 1.00413;
1070    TPCScaleFactor[1879] = 0.990916;
1071    TPCScaleFactor[1881] = 0.987241;
1072    TPCScaleFactor[2186] = 1.00048;
1073    TPCScaleFactor[2187] = 1.00057;
1074    TPCScaleFactor[2254] = 1.00588;
1075    TPCScaleFactor[2256] = 0.991503;
1076    TPCScaleFactor[2321] = 1.00057;
1077
1078
1079   // set all missing values to the value of the run before it ....
1080   for (int i=0; i<(fHighRunN-fLowRunN); i++) {
1081     if (V0MScaleFactor[i] == 0.0) {
1082       if (i==0) {
1083         V0MScaleFactor[i] = 1.00;
1084       } else {
1085         // search for last run number with non-zero value
1086         for (int j=i-1;j>=0;j--) {
1087           if (V0MScaleFactor[j] != 0.0) {
1088             V0MScaleFactor[i] = V0MScaleFactor[j];
1089             break;
1090           }
1091         }
1092       }
1093     }
1094   } // end loop over checking all run-numbers
1095
1096   for (int i=0; i<(fHighRunN-fLowRunN); i++) {
1097     if (SPDScaleFactor[i] == 0.0) {
1098       if (i==0) {
1099         SPDScaleFactor[i] = 1.00;
1100       } else {
1101         for (int j=i-1;j>=0;j--) {
1102           if (SPDScaleFactor[j] != 0.0) {
1103             SPDScaleFactor[i] = SPDScaleFactor[j];
1104             break;
1105           }
1106         }
1107       }
1108     }
1109   }
1110
1111   for (int i=0; i<(fHighRunN-fLowRunN); i++) {
1112     if (TPCScaleFactor[i] == 0.0) {
1113       if (i==0) {
1114         TPCScaleFactor[i] = 1.00;
1115       } else {
1116         for (int j=i-1;j>=0;j--) {
1117           if (TPCScaleFactor[j] != 0.0) {
1118             TPCScaleFactor[i] = TPCScaleFactor[j];
1119             break;
1120           }
1121         }
1122       }
1123     }
1124   }
1125
1126
1127   //    for (int i=0; i<(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << V0MScaleFactor[i] << " for Run " << i+fLowRunN << endl;
1128
1129   return;
1130
1131 }
1132
1133
1134 //________________________________________________________________________
1135 void AliCentralitySelectionTask::MyInitScaleFactorMC()
1136 {
1137 // Initialize the MC scaling factors
1138   for (int i=0; i<(fHighRunN-fLowRunN); i++) V0MScaleFactorMC[i] = 0.0;
1139   // scale factors determined from <V0 charge> on a run-by-run basis
1140   V0MScaleFactorMC[0] = 0.75108;
1141   // set all missing values to the value of the run before it ....
1142   for (int i=0; i<(fHighRunN-fLowRunN); i++) {
1143     if (V0MScaleFactorMC[i] == 0.0) {
1144       if (i==0) {
1145         V0MScaleFactorMC[i] = 1.00;
1146       } else {
1147         // search for last run number with non-zero value
1148         for (int j=i-1;j>=0;j--) {
1149           if (V0MScaleFactorMC[j] != 0.0) {
1150             V0MScaleFactorMC[i] = V0MScaleFactorMC[j];
1151             break;
1152           }
1153         }
1154       }
1155     }
1156   } // end loop over checking all run-numbers
1157
1158
1159   //    for (int i=0; i<(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << V0MScaleFactorMC[i] << " for Run " << i+fLowRunN << endl;
1160
1161   return;
1162
1163 }
1164