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