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