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
 505     esdCent = esd->GetCentrality();
 506
 507     // ***** V0 info
 508     AliESDVZERO* esdV0 = esd->GetVZEROData();
 509     multV0A=esdV0->GetMTotV0A();
 510     multV0C=esdV0->GetMTotV0C();
 511
 512     Float_t v0CorrR;
 513     v0Corr = (Short_t)AliESDUtils::GetCorrV0(esd,v0CorrR);
 514     v0CorrResc = (Short_t)v0CorrR;
 515
 516     // ***** Vertex Info
 517     const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
 518     zvtx        = vtxESD->GetZ();
 519
 520     // ***** CB info (tracklets, clusters, chips)
 521     //nTracks    = event->GetNumberOfTracks();
 522     nTracks    = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;
 523
 524     const AliMultiplicity *mult = esd->GetMultiplicity();
 525
 526     nTracklets = mult->GetNumberOfTracklets();
 527
 528     for(Int_t ilay=0; ilay<6; ilay++){
 529       nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
 530     }
 531
 532     for(Int_t ilay=0; ilay<2; ilay++){
 533       nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
 534     }
 535
 536     spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);
 537
 538     // ***** FMD info
 539     AliESDFMD *fmd = esd->GetFMDData();
 540     Float_t totalMultA = 0;
 541     Float_t totalMultC = 0;
 542     const Float_t fFMDLowCut = 0.4;
 543
 544     for(UShort_t det=1;det<=3;det++) {
 545       Int_t nRings = (det==1 ? 1 : 2);
 546       for (UShort_t ir = 0; ir < nRings; ir++) {
 547         Char_t   ring = (ir == 0 ? 'I' : 'O');
 548         UShort_t nsec = (ir == 0 ? 20  : 40);
 549         UShort_t nstr = (ir == 0 ? 512 : 256);
 550         for(UShort_t sec =0; sec < nsec;  sec++)  {
 551           for(UShort_t strip = 0; strip < nstr; strip++) {
 552
 553             Float_t fmdMult = fmd->Multiplicity(det,ring,sec,strip);
 554             if(fmdMult == 0 || fmdMult == AliESDFMD::kInvalidMult) continue;
 555
 556             Float_t nParticles=0;
 557
 558             if(fmdMult > fFMDLowCut) {
 559               nParticles = 1.;
 560             }
 561
 562             if (det<3) totalMultA = totalMultA + nParticles;
 563             else totalMultC = totalMultC + nParticles;
 564
 565           }
 566         }
 567       }
 568     }
 569     multFMDA = totalMultA;
 570     multFMDC = totalMultC;
 571
 572     // ***** ZDC info
 573     AliESDZDC *esdZDC = esd->GetESDZDC();
 574     zdcEnergyCal = esdZDC->AliESDZDC::TestBit(AliESDZDC::kEnergyCalibratedSignal);
 575     if (zdcEnergyCal) {
 576       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
 577       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
 578       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
 579       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
 580     } else {
 581       zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;
 582       zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;
 583       znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;
 584       zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;
 585     }
 586     zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;
 587     zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;
 588
 589   }
 590
 591   else if(fAnalysisInput.CompareTo("AOD")==0){
 592     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
 593     // to be implemented
 594     printf("  AOD analysis not yet implemented!!!\n\n");
 595     return;
 596   }
 597
 598   // ***** Scaling
 599   // ***** Scaling for pass2
 600   if (fPass==2) {
 601     fUseScaling=kFALSE;
 602   }
 603   // ***** Scaling for MC
 604   if (fIsMCInput) {
 605     fUseScaling=kFALSE;
 606     fUseCleaning=kFALSE;
 607     Float_t tempScalefactorV0M = MyGetScaleFactorMC(fCurrentRun);
 608     v0Corr  = Short_t((multV0A+multV0C)  * tempScalefactorV0M);
 609   }
 610   // ***** Scaling for Data
 611   if (fUseScaling) {
 612     Float_t tempScalefactorV0M = MyGetScaleFactor(fCurrentRun,0);
 613     Float_t tempScalefactorSPD = MyGetScaleFactor(fCurrentRun,1);
 614     Float_t tempScalefactorTPC = MyGetScaleFactor(fCurrentRun,2);
 615     v0Corr  = Short_t(v0Corr / tempScalefactorV0M);
 616     spdCorr = spdCorr / tempScalefactorSPD;
 617     nTracks = Int_t(nTracks / tempScalefactorTPC);
 618   }
 619
 620   // ***** Centrality Selection
 621   if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));
 622   if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
 623   if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
 624   if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
 625   if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
 626   if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));
 627
 628   if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
 629   if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
 630   if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin(zem1Energy+zem2Energy,zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
 631
 632   // ***** Cleaning
 633   if (fUseCleaning) {
 634       fQuality=0;
 635       fZVCut=10;
 636       fOutliersCut=6;
 637
 638       // ***** vertex
 639       if (TMath::Abs(zvtx)>fZVCut) fQuality += 1;
 640
 641       // ***** outliers
 642       // **** V0 vs SPD
 643       if (IsOutlierV0MSPD(spdCorr, v0Corr, int(fCentV0M))) fQuality  += 2;
 644       // ***** V0 vs TPC
 645       if (IsOutlierV0MTPC(nTracks, v0Corr, int(fCentV0M))) fQuality  += 4;
 646       // ***** V0 vs ZDC
 647        if (IsOutlierV0MZDC((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
 648            (zdcEnergyCal==kFALSE) && !(fIsMCInput)) fQuality  += 8;
 649        if (IsOutlierV0MZDCECal((zncEnergy+znaEnergy+zpcEnergy+zpaEnergy), v0Corr) &&
 650            ((zdcEnergyCal==kTRUE) || (fIsMCInput))) fQuality  += 8;
 651   } else {
 652       fQuality = 0;
 653   }
 654
 655
 656   if (esdCent) {
 657       esdCent->SetQuality(fQuality);
 658       esdCent->SetCentralityV0M(fCentV0M);
 659       esdCent->SetCentralityFMD(fCentFMD);
 660       esdCent->SetCentralityTRK(fCentTRK);
 661       esdCent->SetCentralityTKL(fCentTKL);
 662       esdCent->SetCentralityCL0(fCentCL0);
 663       esdCent->SetCentralityCL1(fCentCL1);
 664       esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
 665       esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
 666       esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
 667   }
 668
 669   fHOutQuality->Fill(fQuality);
 670   fHOutVertex->Fill(zvtx);
 671
 672   fHOutMultV0MvsZDN->Fill(v0Corr,(zncEnergy+znaEnergy));
 673   fHOutMultZEMvsZDN->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy));
 674   fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
 675   fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));
 676   fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);
 677   fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);
 678
 679   if (fQuality==0) {
 680     fHOutCentV0M->Fill(fCentV0M);
 681     fHOutCentFMD->Fill(fCentFMD);
 682     fHOutCentTRK->Fill(fCentTRK);
 683     fHOutCentTKL->Fill(fCentTKL);
 684     fHOutCentCL0->Fill(fCentCL0);
 685     fHOutCentCL1->Fill(fCentCL1);
 686     fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);
 687     fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);
 688     fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);
 689     fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentCL1);
 690     fHOutCentV0MvsCentTRK->Fill(fCentV0M,fCentTRK);
 691     fHOutCentTRKvsCentCL1->Fill(fCentTRK,fCentCL1);
 692     fHOutCentV0MvsCentCL1->Fill(fCentV0M,fCentZEMvsZDC);
 693     fHOutMultV0M->Fill(v0Corr);
 694     fHOutMultV0R->Fill(multV0A+multV0C);
 695     fHOutMultFMD->Fill((multFMDA+multFMDC));
 696     fHOutMultTRK->Fill(nTracks);
 697     fHOutMultTKL->Fill(nTracklets);
 698     fHOutMultCL0->Fill(nClusters[0]);
 699     fHOutMultCL1->Fill(spdCorr);
 700     fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);
 701   } else if (fQuality ==1) {
 702     fHOutCentV0Mqual1->Fill(fCentV0M);
 703     fHOutCentTRKqual1->Fill(fCentTRK);
 704     fHOutCentCL1qual1->Fill(fCentCL1);
 705   } else {
 706     fHOutCentV0Mqual2->Fill(fCentV0M);
 707     fHOutCentTRKqual2->Fill(fCentTRK);
 708     fHOutCentCL1qual2->Fill(fCentCL1);
 709   }
 710
 711   PostData(1, fOutputList);
 712 }
 713
 714 //________________________________________________________________________
 715 void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename)
 716 {
 717   //  Read centrality histograms
 718   TDirectory *owd = gDirectory;
 719   // Check if the file is present
 720   TString path = gSystem->ExpandPathName(fCentfilename.Data());
 721   if (gSystem->AccessPathName(path)) {
 722      AliError(Form("File %s does not exist", path.Data()));
 723      return;
 724   }
 725   fFile  = TFile::Open(fCentfilename);
 726   owd->cd();
 727   fHtempV0M  = (TH1F*) (fFile->Get("hmultV0_percentile"));
 728   fHtempFMD  = (TH1F*) (fFile->Get("hmultFMD_percentile"));
 729   fHtempTRK  = (TH1F*) (fFile->Get("hNtracks_percentile"));
 730   fHtempTKL  = (TH1F*) (fFile->Get("hNtracklets_percentile"));
 731   fHtempCL0  = (TH1F*) (fFile->Get("hNclusters0_percentile"));
 732   fHtempCL1  = (TH1F*) (fFile->Get("hNclusters1_percentile"));
 733
 734   if (!fHtempV0M) AliWarning(Form("Calibration for V0M does not exist in %s", path.Data()));
 735   if (!fHtempFMD) AliWarning(Form("Calibration for FMD does not exist in %s", path.Data()));
 736   if (!fHtempTRK) AliWarning(Form("Calibration for TRK does not exist in %s", path.Data()));
 737   if (!fHtempTKL) AliWarning(Form("Calibration for TKL does not exist in %s", path.Data()));
 738   if (!fHtempCL0) AliWarning(Form("Calibration for CL0 does not exist in %s", path.Data()));
 739   if (!fHtempCL1) AliWarning(Form("Calibration for CL1 does not exist in %s", path.Data()));
 740
 741   owd->cd();
 742 }
 743
 744 //________________________________________________________________________
 745 void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2)
 746 {
 747   //  Read centrality histograms
 748   TDirectory *owd = gDirectory;
 749   TString path = gSystem->ExpandPathName(fCentfilename2.Data());
 750   if (gSystem->AccessPathName(path)) {
 751      AliError(Form("File %s does not exist", path.Data()));
 752      return;
 753   }
 754   fFile2  = TFile::Open(fCentfilename2);
 755   owd->cd();
 756   fHtempV0MvsFMD =  (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
 757   fHtempTKLvsV0M  = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
 758   fHtempZEMvsZDC  = (TH2F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
 759
 760   if (!fHtempV0MvsFMD) AliWarning(Form("Calibration for V0MvsFMD does not exist in %s", path.Data()));
 761   if (!fHtempTKLvsV0M) AliWarning(Form("Calibration for TKLvsV0M does not exist in %s", path.Data()));
 762   if (!fHtempZEMvsZDC) AliWarning(Form("Calibration for ZEMvsZDC does not exist in %s", path.Data()));
 763
 764   owd->cd();
 765 }
 766
 767 //________________________________________________________________________
 768 void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
 769 {
 770   // Terminate analysis
 771   if (fFile && fFile->IsOpen())
 772     fFile->Close();
 773   if (fFile2 && fFile2->IsOpen())
 774     fFile2->Close();
 775 }
 776 //________________________________________________________________________
 777 Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* const esd)
 778 {
 779   // Setup files for run
 780
 781   if (!esd)
 782     return -1;
 783
 784   // check if something to be done
 785   if (fCurrentRun == esd->GetRunNumber())
 786     return 0;
 787   else
 788     fCurrentRun = esd->GetRunNumber();
 789
 790   AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));
 791
 792   // CHANGE HERE FOR RUN RANGES
 793   switch (fPass) {
 794     case 1:
 795       if ( fCurrentRun <= 137165 ) fRunNo = 137161;
 796       else fRunNo = 137366;
 797       break;
 798     case 2:
 799       if ( fCurrentRun >= 136851  && fCurrentRun <= 137165 ) fRunNo = 137161;
 800       else if ( ( fCurrentRun >= 137230  && fCurrentRun <= 137848 ) ||
 801                 ( fCurrentRun >= 138190  && fCurrentRun <= 138275 )) fRunNo = 137722;
 802       else if ( fCurrentRun >= 138125  && fCurrentRun <= 138154 ) fRunNo = 138150;
 803       else fRunNo = 139172;
 804       break;
 805     default:
 806       AliError(Form("Run %d not known to centrality selection!", fCurrentRun));
 807       return -1;
 808   }
 809   TString fileName =(Form("%s/COMMON/CENTRALITY/data/pass%d/AliCentralityBy1D_%d.root", AliAnalysisManager::GetOADBPath(), fPass, fRunNo));
 810   TString fileName2=(Form("%s/COMMON/CENTRALITY/data/pass%d/AliCentralityByFunction_%d.root", AliAnalysisManager::GetOADBPath(), fPass, fRunNo));
 811
 812   AliInfo(Form("Centrality Selection for run %d and pass %d is initialized with %s", fCurrentRun, fPass, fileName.Data()));
 813   ReadCentralityHistos(fileName.Data());
 814   ReadCentralityHistos2(fileName2.Data());
 815   if (!fFile && !fFile2) {
 816      AliError(Form("Run %d not known to centrality selection!", fCurrentRun));
 817      return -1;
 818   }
 819   return 0;
 820 }
 821
 822 //________________________________________________________________________
 823 Bool_t AliCentralitySelectionTask::IsOutlierV0MSPD(Float_t spd, Float_t v0, Int_t cent) const
 824 {
 825 // Clean outliers
 826   Float_t val= -0.579579 +  0.273949 * v0;
 827   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};
 828
 829   if ( TMath::Abs(spd-val) > fOutliersCut*spdSigma[cent] )
 830     return kTRUE;
 831   else
 832     return kFALSE;
 833 }
 834
 835 //________________________________________________________________________
 836 Bool_t AliCentralitySelectionTask::IsOutlierV0MTPC(Int_t tracks, Float_t v0, Int_t cent) const
 837 {
 838 // Clean outliers
 839   Float_t val = -1.03873 +  0.125691 * v0;
 840   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};
 841
 842   if ( TMath::Abs(tracks-val) > fOutliersCut*tpcSigma[cent] )
 843     return kTRUE;
 844   else
 845     return kFALSE;
 846 }
 847
 848 //________________________________________________________________________
 849 Bool_t AliCentralitySelectionTask::IsOutlierV0MZDC(Float_t zdc, Float_t v0) const
 850 {
 851 // Clean outliers
 852   Float_t val = 235000. - 9.5 * v0;
 853   if (zdc >  val)
 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   switch (fPass) {
1151   case 1:
1152     fV0MScaleFactorMC[0] = 0.75108;
1153     break;
1154   case 2:
1155     fV0MScaleFactorMC[0] = 0.8;
1156     break;
1157   default:
1158     AliError(Form("Unknown reconstruction pass (%d)! Setting MC scale in V0 to 1", fPass));
1159     fV0MScaleFactorMC[0] = 1.0;
1160     break;
1161   }
1162
1163   // set all missing values to the value of the run before it ....
1164   for (int i=0; i<=(fHighRunN-fLowRunN); i++) {
1165     if (fV0MScaleFactorMC[i] == 0.0) {
1166       if (i==0) {
1167         fV0MScaleFactorMC[i] = 1.00;
1168       } else {
1169         // search for last run number with non-zero value
1170         for (int j=i-1;j>=0;j--) {
1171           if (fV0MScaleFactorMC[j] != 0.0) {
1172             fV0MScaleFactorMC[i] = fV0MScaleFactorMC[j];
1173             break;
1174           }
1175         }
1176       }
1177     }
1178   } // end loop over checking all run-numbers
1179
1180
1181   //    for (int i=0; i<=(fHighRunN-fLowRunN); i++) cout << "Scale Factor = " << fV0MScaleFactorMC[i] << " for Run " << i+fLowRunN << endl;
1182
1183   return;
1184
1185 }
1186