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