ANALYSIS/AliCentralitySelectionTask.cxx

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