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