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