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