ANALYSIS/AliCentralitySelectionTask.cxx

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