]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ANALYSIS/AliCentralitySelectionTask.cxx
Added pt-dep min d0 cut
[u/mrichter/AliRoot.git] / ANALYSIS / AliCentralitySelectionTask.cxx
index f3f157fb408b47acc859e15094fd21394e165963..966a6969a1d702d33033813361ddc8b399b5c4fa 100644 (file)
-/**************************************************************************
- * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-
-//*****************************************************
-//   Class AliCentralitySelectionTask
-//   Class to analyze determine centrality            
-//   author: Alberica Toia
-//*****************************************************
-
-#include <TTree.h>
-#include <TList.h>
-#include <TH1F.h>
-#include <TH2F.h>
-#include <TProfile.h>
-#include <TFile.h>
-#include <TString.h>
-#include <TCanvas.h>
-#include <iostream>
-
-#include "AliAnalysisManager.h"
-#include "AliVEvent.h"
-#include "AliESD.h"
-#include "AliESDEvent.h"
-#include "AliESDHeader.h"
-#include "AliESDInputHandler.h"
-#include "AliESDZDC.h"
-#include "AliESDFMD.h"
-#include "AliESDVZERO.h"
-#include "AliESDCentrality.h"
-#include "AliMultiplicity.h"
-#include "AliAODHandler.h"
-#include "AliAODEvent.h"
-#include "AliAODVertex.h"
-#include "AliAODMCHeader.h"
-#include "AliMCEvent.h"
-#include "AliMCEventHandler.h"
-#include "AliMCParticle.h"
-#include "AliStack.h"
-#include "AliHeader.h"
-#include "AliAODMCParticle.h"
-#include "AliAnalysisTaskSE.h"
-#include "AliGenEventHeader.h"
-#include "AliGenHijingEventHeader.h"
-#include "AliPhysicsSelectionTask.h"
-#include "AliPhysicsSelection.h"
-#include "AliBackgroundSelection.h"
-#include "AliCentralitySelectionTask.h"
-
-ClassImp(AliCentralitySelectionTask)
-
-
-//________________________________________________________________________
-AliCentralitySelectionTask::AliCentralitySelectionTask():
-AliAnalysisTaskSE(),
-  fDebug(0),
-  fAnalysisInput("ESD"),
-  fIsMCInput(kFALSE),
-  fFile(0),
-  fFile2(0),
-  fCentfilename(0),
-  fCentfilename2(0),
-  fCentV0M(0),
-  fCentFMD(0),
-  fCentTRK(0),
-  fCentTKL(0),
-  fCentCL0(0),
-  fCentV0MvsFMD(0),
-  fCentTKLvsV0M(0),
-  fCentZEMvsZDC(0),
-  fHtempV0M(0),
-  fHtempFMD(0),
-  fHtempTRK(0),
-  fHtempTKL(0),
-  fHtempCL0(0),
-  fHtempV0MvsFMD(0),
-  fHtempTKLvsV0M(0),
-  fHtempZEMvsZDC(0)
-{   
-  // Default constructor
-}   
-
-//________________________________________________________________________
-AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
-  AliAnalysisTaskSE(name),
-  fDebug(0),
-  fAnalysisInput("ESD"),
-  fIsMCInput(kFALSE),
-  fFile(0),
-  fFile2(0),
-  fCentfilename(0),
-  fCentfilename2(0),
-  fCentV0M(0),
-  fCentFMD(0),
-  fCentTRK(0),
-  fCentTKL(0),
-  fCentCL0(0),
-  fCentV0MvsFMD(0),
-  fCentTKLvsV0M(0),
-  fCentZEMvsZDC(0),
-  fHtempV0M(0),
-  fHtempFMD(0),
-  fHtempTRK(0),
-  fHtempTKL(0),
-  fHtempCL0(0),
-  fHtempV0MvsFMD(0),
-  fHtempTKLvsV0M(0),
-  fHtempZEMvsZDC(0)
-{
-  // Default constructor
-}
-
-//________________________________________________________________________
-AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
-{
-  // Assignment operator
-  if (this!=&c) {
-    AliAnalysisTaskSE::operator=(c);
-  }
-  return *this;
-}
-
-//________________________________________________________________________
-AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
-  AliAnalysisTaskSE(ana),
-  fDebug(ana.fDebug),    
-  fAnalysisInput(ana.fDebug),
-  fIsMCInput(ana.fIsMCInput),
-  fFile(ana.fFile),
-  fFile2(ana.fFile2),
-  fCentfilename(ana.fCentfilename),
-  fCentfilename2(ana.fCentfilename2),
-  fCentV0M(ana.fCentV0M),
-  fCentFMD(ana.fCentFMD),
-  fCentTRK(ana.fCentTRK),
-  fCentTKL(ana.fCentTKL),
-  fCentCL0(ana.fCentCL0),
-  fCentV0MvsFMD(ana.fCentV0MvsFMD),
-  fCentTKLvsV0M(ana.fCentTKLvsV0M),
-  fCentZEMvsZDC(ana.fCentZEMvsZDC),
-  fHtempV0M(ana.fHtempV0M),
-  fHtempFMD(ana.fHtempFMD),
-  fHtempTRK(ana.fHtempTRK),
-  fHtempTKL(ana.fHtempTKL),
-  fHtempCL0(ana.fHtempCL0),
-  fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
-  fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
-  fHtempZEMvsZDC(ana.fHtempZEMvsZDC)
-{
-  // Copy Constructor  
-}
-//________________________________________________________________________
- AliCentralitySelectionTask::~AliCentralitySelectionTask()
- {
-   // Destructor
- }  
-
-//________________________________________________________________________
-void AliCentralitySelectionTask::UserCreateOutputObjects()
-{  
-  // Create the output containers
-  if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
-}
-
-//________________________________________________________________________
-void AliCentralitySelectionTask::UserExec(Option_t */*option*/)
-{ 
-  // Execute analysis for current event:
-  if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
-  
-  Float_t  zncEnergy;               //  ZNC Energy
-  Float_t  zpcEnergy;               //  ZPC Energy
-  Float_t  znaEnergy;               //  ZNA Energy
-  Float_t  zpaEnergy;               //  ZPA Energy
-  Float_t  zem1Energy = 0.;         //  ZEM1 Energy
-  Float_t  zem2Energy = 0.;         //  ZEM2 Energy
-  
-  Int_t    nTracks    = 0;          //  no. tracks
-  Int_t    nTracklets = 0;          //  no. tracklets
-  Int_t    nClusters[6];            //  no. clusters on 6 ITS layers
-  Int_t    nChips[2];               //  no. chips on 2 SPD layers
-
-  Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A
-  Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C
-  Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A
-  Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C
-
-  AliESDCentrality *esdCent = 0;
-
-  if(fAnalysisInput.CompareTo("ESD")==0){
-
-    AliVEvent* event = InputEvent();
-    AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
-
-    esdCent = esd->GetCentrality();
-
-    // ***** V0 info    
-    AliESDVZERO* esdV0 = esd->GetVZEROData();
-    multV0A=esdV0->GetMTotV0A();
-    multV0C=esdV0->GetMTotV0C();
-    
-    // ***** CB info (tracklets, clusters, chips)
-    nTracks    = event->GetNumberOfTracks();     
-
-    const AliMultiplicity *mult = esd->GetMultiplicity();
-
-    nTracklets = mult->GetNumberOfTracklets();
-
-    for(Int_t ilay=0; ilay<6; ilay++){
-      nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
-    }
-
-    for(Int_t ilay=0; ilay<2; ilay++){
-      nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
-    }
-    
-
-    // ***** FMD info
-    AliESDFMD *fmd = esd->GetFMDData();
-    Float_t totalMultA = 0;
-    Float_t totalMultC = 0;
-    const Float_t fFMDLowCut = 0.4;
-    
-    for(UShort_t det=1;det<=3;det++) {
-      Int_t nRings = (det==1 ? 1 : 2);
-      for (UShort_t ir = 0; ir < nRings; ir++) {         
-       Char_t   ring = (ir == 0 ? 'I' : 'O');
-       UShort_t nsec = (ir == 0 ? 20  : 40);
-       UShort_t nstr = (ir == 0 ? 512 : 256);
-       for(UShort_t sec =0; sec < nsec;  sec++)  {
-         for(UShort_t strip = 0; strip < nstr; strip++) {
-           
-           Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
-           if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
-           
-           Float_t nParticles=0;
-           
-           if(FMDmult > fFMDLowCut) {
-             nParticles = 1.;
-           }
-           
-           if (det<3) totalMultA = totalMultA + nParticles;
-           else totalMultC = totalMultC + nParticles;
-           
-         }
-       }
-      }
-    }
-    multFMDA = totalMultA;
-    multFMDC = totalMultC;
-    
-    // ***** ZDC info
-    AliESDZDC *esdZDC = esd->GetESDZDC();
-    zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
-    zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
-    znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
-    zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
-    zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
-    zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
-    
-  }   
-  else if(fAnalysisInput.CompareTo("AOD")==0){
-    //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
-    // to be implemented
-    printf("  AOD analysis not yet implemented!!!\n\n");
-    return;
-  }
-
-  // ***** Centrality Selection
-  fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((multV0A+multV0C)));
-  fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
-  fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
-  fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
-  fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
-  
-  fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
-  fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
-  fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));
-  
-  esdCent->SetCentralityV0M(fCentV0M);
-  esdCent->SetCentralityFMD(fCentFMD);
-  esdCent->SetCentralityTRK(fCentTRK);
-  esdCent->SetCentralityTKL(fCentTKL);
-  esdCent->SetCentralityCL0(fCentCL0);
-  esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
-  esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
-  esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
-}
-
-//________________________________________________________________________
-void AliCentralitySelectionTask::ReadCentralityHistos() 
-{
-  fFile  = new TFile(fCentfilename);
-  fHtempV0M  = (TH1D*) (fFile->Get("hmultV0_percentile")); 
-  fHtempFMD  = (TH1D*) (fFile->Get("hmultFMD_percentile")); 
-  fHtempTRK  = (TH1D*) (fFile->Get("hNtracks_percentile")); 
-  fHtempTKL  = (TH1D*) (fFile->Get("hNtracklets_percentile")); 
-  fHtempCL0  = (TH1D*) (fFile->Get("hNclusters0_percentile")); 
-}  
-  
-void AliCentralitySelectionTask::ReadCentralityHistos2() 
-{
-  fFile2  = new TFile(fCentfilename2);
-  fHtempV0MvsFMD  = (TH1D*) (fFile2->Get("hmultV0vsmultFMD_all_percentile")); 
-  fHtempTKLvsV0M  = (TH1D*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile")); 
-  fHtempZEMvsZDC  = (TH1D*) (fFile2->Get("hEzemvsEzdc_all_percentile"));  
-}
-
-void AliCentralitySelectionTask::SetPercentileFile(TString filename) 
-{
-  fCentfilename = filename;
-  ReadCentralityHistos();
-}
-
-void AliCentralitySelectionTask::SetPercentileFile2(TString filename) 
-{
-  fCentfilename2 = filename;
-  ReadCentralityHistos2();
-}
-
-//________________________________________________________________________
-void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
-{
-  // Terminate analysis
-  fFile->Close();  
-  fFile2->Close();  
-}
-//________________________________________________________________________
-
+/**************************************************************************\r
+ * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *\r
+ *                                                                        *\r
+ * Author: The ALICE Off-line Project.                                    *\r
+ * Contributors are mentioned in the code where appropriate.              *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+\r
+//*****************************************************\r
+//   Class AliCentralitySelectionTask\r
+//   Class to analyze determine centrality            \r
+//   author: Alberica Toia\r
+//*****************************************************\r
+\r
+#include "AliCentralitySelectionTask.h"\r
+\r
+#include <TTree.h>\r
+#include <TList.h>\r
+#include <TH1F.h>\r
+#include <TH2F.h>\r
+#include <TProfile.h>\r
+#include <TFile.h>\r
+#include <TObjString.h>\r
+#include <TString.h>\r
+#include <TCanvas.h>\r
+#include <TROOT.h>\r
+#include <TDirectory.h>\r
+#include <TSystem.h>\r
+#include <iostream>\r
+\r
+#include "AliAnalysisManager.h"\r
+#include "AliVEvent.h"\r
+#include "AliESD.h"\r
+#include "AliESDEvent.h"\r
+#include "AliESDHeader.h"\r
+#include "AliESDInputHandler.h"\r
+#include "AliESDZDC.h"\r
+#include "AliESDFMD.h"\r
+#include "AliESDVZERO.h"\r
+#include "AliCentrality.h"\r
+#include "AliESDtrackCuts.h"\r
+#include "AliMultiplicity.h"\r
+#include "AliAODHandler.h"\r
+#include "AliAODEvent.h"\r
+#include "AliESDVertex.h"\r
+#include "AliAODVertex.h"\r
+#include "AliAODMCHeader.h"\r
+#include "AliMCEvent.h"\r
+#include "AliMCEventHandler.h"\r
+#include "AliMCParticle.h"\r
+#include "AliStack.h"\r
+#include "AliHeader.h"\r
+#include "AliAODMCParticle.h"\r
+#include "AliAnalysisTaskSE.h"\r
+#include "AliGenEventHeader.h"\r
+#include "AliGenHijingEventHeader.h"\r
+#include "AliPhysicsSelectionTask.h"\r
+#include "AliPhysicsSelection.h"\r
+#include "AliBackgroundSelection.h"\r
+#include "AliESDUtils.h"\r
+\r
+ClassImp(AliCentralitySelectionTask)\r
+\r
+\r
+//________________________________________________________________________\r
+AliCentralitySelectionTask::AliCentralitySelectionTask():\r
+AliAnalysisTaskSE(),\r
+  fDebug(0),\r
+  fAnalysisInput("ESD"),\r
+  fIsMCInput(kFALSE),\r
+  fDirName("$ALICE_ROOT/OADB/COMMON/CENTRALITY/data"),\r
+  fFile(0),\r
+  fFile2(0),\r
+  fCurrentRun(-1),\r
+  fRunNo(-1),\r
+  fTrackCuts(0),\r
+  fCentV0M(0),\r
+  fCentFMD(0),\r
+  fCentTRK(0),\r
+  fCentTKL(0),\r
+  fCentCL0(0),\r
+  fCentCL1(0),\r
+  fCentV0MvsFMD(0),\r
+  fCentTKLvsV0M(0),\r
+  fCentZEMvsZDC(0),\r
+  fHtempV0M(0),\r
+  fHtempFMD(0),\r
+  fHtempTRK(0),\r
+  fHtempTKL(0),\r
+  fHtempCL0(0),\r
+  fHtempCL1(0),\r
+  fHtempV0MvsFMD(0),\r
+  fHtempTKLvsV0M(0),\r
+  fHtempZEMvsZDC(0),\r
+  fOutputList(0),\r
+  fHOutCentV0M     (0),\r
+  fHOutCentFMD     (0),\r
+  fHOutCentTRK     (0),\r
+  fHOutCentTKL     (0),\r
+  fHOutCentCL0     (0),\r
+  fHOutCentCL1     (0),\r
+  fHOutCentV0MvsFMD(0),\r
+  fHOutCentTKLvsV0M(0),\r
+  fHOutCentZEMvsZDC(0),\r
+  fHOutMultV0M(0),\r
+  fHOutMultV0R(0),\r
+  fHOutMultFMD(0),\r
+  fHOutMultTRK(0),\r
+  fHOutMultTKL(0),\r
+  fHOutMultCL0(0),\r
+  fHOutMultCL1(0),\r
+  fHOutMultV0MvsZDC(0),\r
+  fHOutMultZEMvsZDC(0),\r
+  fHOutMultV0MvsCL1(0),\r
+  fHOutMultV0MvsTRK(0),\r
+  fHOutMultTRKvsCL1(0)\r
+{   \r
+  // Default constructor\r
+  AliInfo("Centrality Selection enabled.");\r
+}   \r
+\r
+//________________________________________________________________________\r
+AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):\r
+  AliAnalysisTaskSE(name),\r
+  fDebug(0),\r
+  fAnalysisInput("ESD"),\r
+  fIsMCInput(kFALSE),\r
+  fDirName("$ALICE_ROOT/OADB/COMMON/CENTRALITY/data"),\r
+  fFile(0),\r
+  fFile2(0),\r
+  fCurrentRun(-1),\r
+  fRunNo(-1),\r
+  fTrackCuts(0),\r
+  fCentV0M(0),\r
+  fCentFMD(0),\r
+  fCentTRK(0),\r
+  fCentTKL(0),\r
+  fCentCL0(0),\r
+  fCentCL1(0),\r
+  fCentV0MvsFMD(0),\r
+  fCentTKLvsV0M(0),\r
+  fCentZEMvsZDC(0),\r
+  fHtempV0M(0),\r
+  fHtempFMD(0),\r
+  fHtempTRK(0),\r
+  fHtempTKL(0),\r
+  fHtempCL0(0),\r
+  fHtempCL1(0),\r
+  fHtempV0MvsFMD(0),\r
+  fHtempTKLvsV0M(0),\r
+  fHtempZEMvsZDC(0),\r
+  fOutputList(0),\r
+  fHOutCentV0M     (0),\r
+  fHOutCentFMD     (0),\r
+  fHOutCentTRK     (0),\r
+  fHOutCentTKL     (0),\r
+  fHOutCentCL0     (0),\r
+  fHOutCentCL1     (0),\r
+  fHOutCentV0MvsFMD(0),\r
+  fHOutCentTKLvsV0M(0),\r
+  fHOutCentZEMvsZDC(0),\r
+  fHOutMultV0M(0),\r
+  fHOutMultV0R(0),\r
+  fHOutMultFMD(0),\r
+  fHOutMultTRK(0),\r
+  fHOutMultTKL(0),\r
+  fHOutMultCL0(0),\r
+  fHOutMultCL1(0),\r
+  fHOutMultV0MvsZDC(0),\r
+  fHOutMultZEMvsZDC(0),\r
+  fHOutMultV0MvsCL1(0),\r
+  fHOutMultV0MvsTRK(0),\r
+  fHOutMultTRKvsCL1(0)\r
+{\r
+  // Default constructor\r
+  AliInfo("Centrality Selection enabled.");\r
+  DefineOutput(1, TList::Class());\r
+}\r
+\r
+//________________________________________________________________________\r
+AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)\r
+{\r
+  // Assignment operator\r
+  if (this!=&c) {\r
+    AliAnalysisTaskSE::operator=(c);\r
+  }\r
+  return *this;\r
+}\r
+\r
+//________________________________________________________________________\r
+AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):\r
+  AliAnalysisTaskSE(ana),\r
+  fDebug(ana.fDebug),    \r
+  fAnalysisInput(ana.fDebug),\r
+  fIsMCInput(ana.fIsMCInput),\r
+  fDirName(ana.fDirName),\r
+  fFile(ana.fFile),\r
+  fFile2(ana.fFile2),\r
+  fCurrentRun(ana.fCurrentRun),\r
+  fRunNo(ana.fRunNo),\r
+  fTrackCuts(ana.fTrackCuts),\r
+  fCentV0M(ana.fCentV0M),\r
+  fCentFMD(ana.fCentFMD),\r
+  fCentTRK(ana.fCentTRK),\r
+  fCentTKL(ana.fCentTKL),\r
+  fCentCL0(ana.fCentCL0),\r
+  fCentCL1(ana.fCentCL1),\r
+  fCentV0MvsFMD(ana.fCentV0MvsFMD),\r
+  fCentTKLvsV0M(ana.fCentTKLvsV0M),\r
+  fCentZEMvsZDC(ana.fCentZEMvsZDC),\r
+  fHtempV0M(ana.fHtempV0M),\r
+  fHtempFMD(ana.fHtempFMD),\r
+  fHtempTRK(ana.fHtempTRK),\r
+  fHtempTKL(ana.fHtempTKL),\r
+  fHtempCL0(ana.fHtempCL0),\r
+  fHtempCL1(ana.fHtempCL1),\r
+  fHtempV0MvsFMD(ana.fHtempV0MvsFMD),\r
+  fHtempTKLvsV0M(ana.fHtempTKLvsV0M),\r
+  fHtempZEMvsZDC(ana.fHtempZEMvsZDC),\r
+  fOutputList(ana.fOutputList),\r
+  fHOutCentV0M     (ana.fHOutCentV0M     ),\r
+  fHOutCentFMD     (ana.fHOutCentFMD     ),\r
+  fHOutCentTRK     (ana.fHOutCentTRK     ),\r
+  fHOutCentTKL     (ana.fHOutCentTKL     ),\r
+  fHOutCentCL0     (ana.fHOutCentCL0     ),\r
+  fHOutCentCL1     (ana.fHOutCentCL1     ),\r
+  fHOutCentV0MvsFMD(ana.fHOutCentV0MvsFMD),\r
+  fHOutCentTKLvsV0M(ana.fHOutCentTKLvsV0M),\r
+  fHOutCentZEMvsZDC(ana.fHOutCentZEMvsZDC),\r
+  fHOutMultV0M(ana.fHOutMultV0M),\r
+  fHOutMultV0R(ana.fHOutMultV0R),\r
+  fHOutMultFMD(ana.fHOutMultFMD),\r
+  fHOutMultTRK(ana.fHOutMultTRK),\r
+  fHOutMultTKL(ana.fHOutMultTKL),\r
+  fHOutMultCL0(ana.fHOutMultCL0),\r
+  fHOutMultCL1(ana.fHOutMultCL1),\r
+  fHOutMultV0MvsZDC(ana.fHOutMultV0MvsZDC),\r
+  fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC),\r
+  fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1),\r
+  fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK),\r
+  fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1)\r
+{\r
+  // Copy Constructor  \r
+}\r
\r
+//________________________________________________________________________\r
+AliCentralitySelectionTask::~AliCentralitySelectionTask()\r
+{\r
+  // Destructor  \r
+  if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList;\r
+  if (fTrackCuts) delete fTrackCuts;\r
+}  \r
+\r
+//________________________________________________________________________\r
+void AliCentralitySelectionTask::UserCreateOutputObjects()\r
+{  \r
+  // Create the output containers\r
+  if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");\r
+  AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo);\r
+\r
+  fOutputList = new TList();\r
+  fOutputList->SetOwner();\r
+  fHOutCentV0M     = new TH1F("fHOutCentV0M","fHOutCentV0M; Centrality V0",101,-0.5,100.5);\r
+  fHOutCentFMD     = new TH1F("fHOutCentFMD","fHOutCentFMD; Centrality FMD",101,-0.5,100.5);\r
+  fHOutCentTRK     = new TH1F("fHOutCentTRK","fHOutCentTRK; Centrality TPC",101,-0.5,100.5);\r
+  fHOutCentTKL     = new TH1F("fHOutCentTKL","fHOutCentTKL; Centrality tracklets",101,-0.5,100.5);\r
+  fHOutCentCL0     = new TH1F("fHOutCentCL0","fHOutCentCL0; Centrality SPD inner",101,-0.5,100.5);\r
+  fHOutCentCL1     = new TH1F("fHOutCentCL1","fHOutCentCL1; Centrality SPD outer",101,-0.5,100.5);\r
+  fHOutCentV0MvsFMD= new TH1F("fHOutCentV0MvsFMD","fHOutCentV0MvsFMD; Centrality V0 vs FMD",101,-0.5,100.5);\r
+  fHOutCentTKLvsV0M= new TH1F("fHOutCentTKLvsV0M","fHOutCentTKLvsV0M; Centrality tracklets vs V0",101,-0.5,100.5);\r
+  fHOutCentZEMvsZDC= new TH1F("fHOutCentZEMvsZDC","fHOutCentZEMvsZDC; Centrality ZEM vs ZDC",101,-0.5,100.5);\r
+\r
+  fHOutMultV0M = new TH1F("fHOutMultV0M","fHOutMultV0M; Multiplicity V0",25000,0,25000);\r
+  fHOutMultV0R = new TH1F("fHOutMultV0R","fHOutMultV0R; Multiplicity V0",25000,0,25000);\r
+  fHOutMultFMD = new TH1F("fHOutMultFMD","fHOutMultFMD; Multiplicity FMD",24000,0,24000);\r
+  fHOutMultTRK = new TH1F("fHOutMultTRK","fHOutMultTRK; Multiplicity TPC",4000,0,4000);\r
+  fHOutMultTKL = new TH1F("fHOutMultTKL","fHOutMultTKL; Multiplicity tracklets",5000,0,5000);\r
+  fHOutMultCL0 = new TH1F("fHOutMultCL0","fHOutMultCL0; Multiplicity SPD inner",7000,0,7000);\r
+  fHOutMultCL1 = new TH1F("fHOutMultCL1","fHOutMultCL1; Multiplicity SPD outer",7000,0,7000);\r
+  fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,25000,500,0,6000);\r
+  fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,6000);\r
+  fHOutMultV0MvsCL1 = new TH2F("fHOutMultV0MvsCL1","fHOutMultV0MvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000);\r
+  fHOutMultV0MvsTRK = new TH2F("fHOutMultV0MvsTRK","fHOutMultV0MvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000);\r
+  fHOutMultTRKvsCL1 = new TH2F("fHOutMultTRKvsCL1","fHOutMultTRKvsCL1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000);\r
+\r
+  fOutputList->Add(  fHOutCentV0M     );\r
+  fOutputList->Add(  fHOutCentFMD     );\r
+  fOutputList->Add(  fHOutCentTRK     );\r
+  fOutputList->Add(  fHOutCentTKL     );\r
+  fOutputList->Add(  fHOutCentCL0     );\r
+  fOutputList->Add(  fHOutCentCL1     );\r
+  fOutputList->Add(  fHOutCentV0MvsFMD);\r
+  fOutputList->Add(  fHOutCentTKLvsV0M);\r
+  fOutputList->Add(  fHOutCentZEMvsZDC);\r
+  fOutputList->Add(  fHOutMultV0M); \r
+  fOutputList->Add(  fHOutMultV0R); \r
+  fOutputList->Add(  fHOutMultFMD); \r
+  fOutputList->Add(  fHOutMultTRK); \r
+  fOutputList->Add(  fHOutMultTKL); \r
+  fOutputList->Add(  fHOutMultCL0); \r
+  fOutputList->Add(  fHOutMultCL1); \r
+  fOutputList->Add(  fHOutMultV0MvsZDC);\r
+  fOutputList->Add(  fHOutMultZEMvsZDC);\r
+  fOutputList->Add(  fHOutMultV0MvsCL1);\r
+  fOutputList->Add(  fHOutMultV0MvsTRK);\r
+  fOutputList->Add(  fHOutMultTRKvsCL1);\r
+\r
+\r
+  fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();\r
+\r
+  PostData(1, fOutputList); \r
+}\r
+\r
+//________________________________________________________________________\r
+void AliCentralitySelectionTask::UserExec(Option_t */*option*/)\r
+{ \r
+  // Execute analysis for current event:\r
+  if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");\r
+  \r
+  Float_t  zncEnergy = 0.;          //  ZNC Energy\r
+  Float_t  zpcEnergy = 0.;          //  ZPC Energy\r
+  Float_t  znaEnergy = 0.;          //  ZNA Energy\r
+  Float_t  zpaEnergy = 0.;          //  ZPA Energy\r
+  Float_t  zem1Energy = 0.;         //  ZEM1 Energy\r
+  Float_t  zem2Energy = 0.;         //  ZEM2 Energy\r
+  \r
+  Int_t    nTracks = 0;             //  no. tracks\r
+  Int_t    nTracklets = 0;          //  no. tracklets\r
+  Int_t    nClusters[6];            //  no. clusters on 6 ITS layers\r
+  Int_t    nChips[2];               //  no. chips on 2 SPD layers\r
+  Float_t  spdCorr =0;              //  corrected spd2 multiplicity\r
+\r
+  Float_t  multV0A  = 0;            //  multiplicity from V0 reco side A\r
+  Float_t  multV0C  = 0;            //  multiplicity from V0 reco side C\r
+  Float_t  multFMDA = 0;            //  multiplicity from FMD on detector A\r
+  Float_t  multFMDC = 0;            //  multiplicity from FMD on detector C\r
+\r
+  Short_t v0Corr = 0;               // corrected V0 multiplicity\r
+  Short_t v0CorrResc = 0;           // corrected and rescaled V0 multiplicity\r
+\r
+  Float_t zvtx =0;                  // z-vertex SPD\r
\r
+  AliCentrality *esdCent = 0;\r
+\r
+  if(fAnalysisInput.CompareTo("ESD")==0){\r
+\r
+    AliVEvent* event = InputEvent();\r
+    AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);\r
+  \r
+    if (fRunNo<=0) {\r
+      if (SetupRun(esd)<0)\r
+         AliFatal("Centrality File not available for this run");\r
+    }\r
+\r
+    esdCent = esd->GetCentrality();\r
+\r
+    // ***** V0 info    \r
+    AliESDVZERO* esdV0 = esd->GetVZEROData();\r
+    multV0A=esdV0->GetMTotV0A();\r
+    multV0C=esdV0->GetMTotV0C();\r
+\r
+    float v0CorrR;\r
+    v0Corr = (Short_t)AliESDUtils::GetCorrV0(esd,v0CorrR);\r
+    v0CorrResc = (Short_t)v0CorrR;\r
+\r
+    // ***** Vertex Info\r
+    const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();\r
+    zvtx        = vtxESD->GetZ(); \r
+\r
+    // ***** CB info (tracklets, clusters, chips)\r
+    //nTracks    = event->GetNumberOfTracks();     \r
+    nTracks    = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1;\r
+\r
+    const AliMultiplicity *mult = esd->GetMultiplicity();\r
+\r
+    nTracklets = mult->GetNumberOfTracklets();\r
+\r
+    for(Int_t ilay=0; ilay<6; ilay++){\r
+      nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);\r
+    }\r
+\r
+    for(Int_t ilay=0; ilay<2; ilay++){\r
+      nChips[ilay] = mult->GetNumberOfFiredChips(ilay);\r
+    }\r
+\r
+    spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx);    \r
+\r
+    // ***** FMD info\r
+    AliESDFMD *fmd = esd->GetFMDData();\r
+    Float_t totalMultA = 0;\r
+    Float_t totalMultC = 0;\r
+    const Float_t fFMDLowCut = 0.4;\r
+    \r
+    for(UShort_t det=1;det<=3;det++) {\r
+      Int_t nRings = (det==1 ? 1 : 2);\r
+      for (UShort_t ir = 0; ir < nRings; ir++) {         \r
+       Char_t   ring = (ir == 0 ? 'I' : 'O');\r
+       UShort_t nsec = (ir == 0 ? 20  : 40);\r
+       UShort_t nstr = (ir == 0 ? 512 : 256);\r
+       for(UShort_t sec =0; sec < nsec;  sec++)  {\r
+         for(UShort_t strip = 0; strip < nstr; strip++) {\r
+           \r
+           Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);\r
+           if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;\r
+           \r
+           Float_t nParticles=0;\r
+           \r
+           if(FMDmult > fFMDLowCut) {\r
+             nParticles = 1.;\r
+           }\r
+           \r
+           if (det<3) totalMultA = totalMultA + nParticles;\r
+           else totalMultC = totalMultC + nParticles;\r
+           \r
+         }\r
+       }\r
+      }\r
+    }\r
+    multFMDA = totalMultA;\r
+    multFMDC = totalMultC;\r
+    \r
+    // ***** ZDC info\r
+    AliESDZDC *esdZDC = esd->GetESDZDC();\r
+    zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy())/8.;\r
+    zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.;\r
+    znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.;\r
+    zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.;\r
+    zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.;\r
+    zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.;\r
+    \r
+  }   \r
+  else if(fAnalysisInput.CompareTo("AOD")==0){\r
+    //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());\r
+    // to be implemented\r
+    printf("  AOD analysis not yet implemented!!!\n\n");\r
+    return;\r
+  }\r
+\r
+  // ***** Centrality Selection\r
+  if(fHtempV0M)  fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr)));\r
+  ///  else     printf("  Centrality by V0 not available!!!\n\n");\r
+  if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));\r
+  //  else     printf("  Centrality by FMD not available!!!\n\n");\r
+  if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));\r
+  //  else     printf("  Centrality by TRK not available!!!\n\n");\r
+  if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));\r
+  //  else     printf("  Centrality by TKL not available!!!\n\n");\r
+  if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));\r
+  //  else     printf("  Centrality by CL0 not available!!!\n\n");\r
+  if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr));\r
+  ///  else     printf("  Centrality by CL1 not available!!!\n\n");\r
+  \r
+  if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));\r
+  //  else     printf("  Centrality by V0 vs FMD not available!!!\n\n");\r
+  if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));\r
+  //  else     printf("  Centrality by V0 vs TKL not available!!!\n\n");\r
+  if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));\r
+  //  else     printf("  Centrality by ZEM vs ZDC not available!!!\n\n");\r
+\r
+  esdCent->SetCentralityV0M(fCentV0M);\r
+  esdCent->SetCentralityFMD(fCentFMD);\r
+  esdCent->SetCentralityTRK(fCentTRK);\r
+  esdCent->SetCentralityTKL(fCentTKL);\r
+  esdCent->SetCentralityCL0(fCentCL0);\r
+  esdCent->SetCentralityCL1(fCentCL1);\r
+  esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);\r
+  esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);\r
+  esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);\r
+\r
+  fHOutCentV0M->Fill(fCentV0M);\r
+  fHOutCentFMD->Fill(fCentFMD);\r
+  fHOutCentTRK->Fill(fCentTRK);\r
+  fHOutCentTKL->Fill(fCentTKL);\r
+  fHOutCentCL0->Fill(fCentCL0);\r
+  fHOutCentCL1->Fill(fCentCL1);\r
+  fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD);\r
+  fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M);\r
+  fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC);\r
+  fHOutMultV0M->Fill(v0Corr);\r
+  fHOutMultV0R->Fill(multV0A+multV0C);\r
+  fHOutMultFMD->Fill((multFMDA+multFMDC));\r
+  fHOutMultTRK->Fill(nTracks);\r
+  fHOutMultTKL->Fill(nTracklets);\r
+  fHOutMultCL0->Fill(nClusters[0]);\r
+  fHOutMultCL1->Fill(spdCorr);\r
+  fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));\r
+  fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy));\r
+  fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr);\r
+  fHOutMultV0MvsTRK->Fill(v0Corr,nTracks);\r
+  fHOutMultTRKvsCL1->Fill(nTracks,spdCorr);\r
+\r
+  PostData(1, fOutputList); \r
+}\r
+\r
+//________________________________________________________________________\r
+void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename) \r
+{\r
+  //  Read centrality histograms\r
+  TDirectory *owd = gDirectory;\r
+  // Check if the file is present\r
+  TString path = gSystem->ExpandPathName(fCentfilename.Data());\r
+  if (gSystem->AccessPathName(path)) {\r
+     AliError(Form("File %s does not exist", path.Data()));\r
+     return;\r
+  }\r
+  fFile  = TFile::Open(fCentfilename);\r
+  owd->cd();\r
+  fHtempV0M  = (TH1F*) (fFile->Get("hmultV0_percentile"));\r
+  fHtempFMD  = (TH1F*) (fFile->Get("hmultFMD_percentile"));\r
+  fHtempTRK  = (TH1F*) (fFile->Get("hNtracks_percentile"));\r
+  fHtempTKL  = (TH1F*) (fFile->Get("hNtracklets_percentile"));\r
+  fHtempCL0  = (TH1F*) (fFile->Get("hNclusters0_percentile"));\r
+  fHtempCL1  = (TH1F*) (fFile->Get("hNclusters1_percentile"));\r
+  owd->cd();\r
+}  \r
+\r
+//________________________________________________________________________\r
+void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2) \r
+{\r
+  //  Read centrality histograms\r
+  TDirectory *owd = gDirectory;\r
+  TString path = gSystem->ExpandPathName(fCentfilename2.Data());\r
+  if (gSystem->AccessPathName(path)) {\r
+     AliError(Form("File %s does not exist", path.Data()));\r
+     return;\r
+  }   \r
+  fFile2  = TFile::Open(fCentfilename2);\r
+  owd->cd();\r
+  fHtempV0MvsFMD =  (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));\r
+  fHtempTKLvsV0M  = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));\r
+  fHtempZEMvsZDC  = (TH1F*) (fFile2->Get("hEzemvsEzdc_all_percentile"));\r
+  owd->cd();\r
+}\r
+\r
+//________________________________________________________________________\r
+void AliCentralitySelectionTask::Terminate(Option_t */*option*/)\r
+{\r
+  // Terminate analysis\r
+  if (fFile && fFile->IsOpen())\r
+    fFile->Close();  \r
+  if (fFile2 && fFile2->IsOpen())\r
+    fFile2->Close();  \r
+}\r
+//________________________________________________________________________\r
+Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* esd)\r
+{\r
+  // Setup files for run\r
+\r
+  if (!esd)\r
+    return -1;\r
+\r
+  // check if something to be done\r
+  if (fCurrentRun == esd->GetRunNumber())\r
+    return 0;\r
+  else\r
+    fCurrentRun = esd->GetRunNumber();\r
+  \r
+  AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun));\r
+\r
+  // CHANGE HERE FOR RUN RANGES\r
+  if ( fCurrentRun <= 137165 ) fRunNo = 137161;\r
+  else fRunNo = 137366;\r
+  // CHANGE HERE FOR RUN RANGES\r
+  \r
+  TString fileName(Form("%s/AliCentralityBy1D_%d.root", fDirName.Data(), fRunNo));\r
+  TString fileName2(Form("%s/AliCentralityByFunction_%d.root", fDirName.Data(), fRunNo));\r
+  \r
+  AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data()));\r
+  ReadCentralityHistos(fileName.Data());\r
+  ReadCentralityHistos2(fileName2.Data());\r
+  if (!fFile && !fFile2) {\r
+     AliFatal(Form("Run %d not known to centrality selection!", fCurrentRun));       \r
+     return -1;\r
+  }   \r
+  return 0;\r
+}\r