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