X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ANALYSIS%2FAliCentralitySelectionTask.cxx;h=966a6969a1d702d33033813361ddc8b399b5c4fa;hb=81858e84f574d65dfadef25985ddc29a32e1d311;hp=f3f157fb408b47acc859e15094fd21394e165963;hpb=d15bf53fad4a4bb28f94c1a59f99e8f5a6e35ac2;p=u%2Fmrichter%2FAliRoot.git diff --git a/ANALYSIS/AliCentralitySelectionTask.cxx b/ANALYSIS/AliCentralitySelectionTask.cxx index f3f157fb408..966a6969a1d 100644 --- a/ANALYSIS/AliCentralitySelectionTask.cxx +++ b/ANALYSIS/AliCentralitySelectionTask.cxx @@ -1,341 +1,583 @@ -/************************************************************************** - * 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 -#include -#include -#include -#include -#include -#include -#include -#include - -#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(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 (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(); -} -//________________________________________________________________________ - +/************************************************************************** + * 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 "AliCentralitySelectionTask.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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 "AliCentrality.h" +#include "AliESDtrackCuts.h" +#include "AliMultiplicity.h" +#include "AliAODHandler.h" +#include "AliAODEvent.h" +#include "AliESDVertex.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 "AliESDUtils.h" + +ClassImp(AliCentralitySelectionTask) + + +//________________________________________________________________________ +AliCentralitySelectionTask::AliCentralitySelectionTask(): +AliAnalysisTaskSE(), + fDebug(0), + fAnalysisInput("ESD"), + fIsMCInput(kFALSE), + fDirName("$ALICE_ROOT/OADB/COMMON/CENTRALITY/data"), + fFile(0), + fFile2(0), + fCurrentRun(-1), + fRunNo(-1), + fTrackCuts(0), + fCentV0M(0), + fCentFMD(0), + fCentTRK(0), + fCentTKL(0), + fCentCL0(0), + fCentCL1(0), + fCentV0MvsFMD(0), + fCentTKLvsV0M(0), + fCentZEMvsZDC(0), + fHtempV0M(0), + fHtempFMD(0), + fHtempTRK(0), + fHtempTKL(0), + fHtempCL0(0), + fHtempCL1(0), + fHtempV0MvsFMD(0), + fHtempTKLvsV0M(0), + fHtempZEMvsZDC(0), + fOutputList(0), + fHOutCentV0M (0), + fHOutCentFMD (0), + fHOutCentTRK (0), + fHOutCentTKL (0), + fHOutCentCL0 (0), + fHOutCentCL1 (0), + fHOutCentV0MvsFMD(0), + fHOutCentTKLvsV0M(0), + fHOutCentZEMvsZDC(0), + fHOutMultV0M(0), + fHOutMultV0R(0), + fHOutMultFMD(0), + fHOutMultTRK(0), + fHOutMultTKL(0), + fHOutMultCL0(0), + fHOutMultCL1(0), + fHOutMultV0MvsZDC(0), + fHOutMultZEMvsZDC(0), + fHOutMultV0MvsCL1(0), + fHOutMultV0MvsTRK(0), + fHOutMultTRKvsCL1(0) +{ + // Default constructor + AliInfo("Centrality Selection enabled."); +} + +//________________________________________________________________________ +AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name): + AliAnalysisTaskSE(name), + fDebug(0), + fAnalysisInput("ESD"), + fIsMCInput(kFALSE), + fDirName("$ALICE_ROOT/OADB/COMMON/CENTRALITY/data"), + fFile(0), + fFile2(0), + fCurrentRun(-1), + fRunNo(-1), + fTrackCuts(0), + fCentV0M(0), + fCentFMD(0), + fCentTRK(0), + fCentTKL(0), + fCentCL0(0), + fCentCL1(0), + fCentV0MvsFMD(0), + fCentTKLvsV0M(0), + fCentZEMvsZDC(0), + fHtempV0M(0), + fHtempFMD(0), + fHtempTRK(0), + fHtempTKL(0), + fHtempCL0(0), + fHtempCL1(0), + fHtempV0MvsFMD(0), + fHtempTKLvsV0M(0), + fHtempZEMvsZDC(0), + fOutputList(0), + fHOutCentV0M (0), + fHOutCentFMD (0), + fHOutCentTRK (0), + fHOutCentTKL (0), + fHOutCentCL0 (0), + fHOutCentCL1 (0), + fHOutCentV0MvsFMD(0), + fHOutCentTKLvsV0M(0), + fHOutCentZEMvsZDC(0), + fHOutMultV0M(0), + fHOutMultV0R(0), + fHOutMultFMD(0), + fHOutMultTRK(0), + fHOutMultTKL(0), + fHOutMultCL0(0), + fHOutMultCL1(0), + fHOutMultV0MvsZDC(0), + fHOutMultZEMvsZDC(0), + fHOutMultV0MvsCL1(0), + fHOutMultV0MvsTRK(0), + fHOutMultTRKvsCL1(0) +{ + // Default constructor + AliInfo("Centrality Selection enabled."); + DefineOutput(1, TList::Class()); +} + +//________________________________________________________________________ +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), + fDirName(ana.fDirName), + fFile(ana.fFile), + fFile2(ana.fFile2), + fCurrentRun(ana.fCurrentRun), + fRunNo(ana.fRunNo), + fTrackCuts(ana.fTrackCuts), + fCentV0M(ana.fCentV0M), + fCentFMD(ana.fCentFMD), + fCentTRK(ana.fCentTRK), + fCentTKL(ana.fCentTKL), + fCentCL0(ana.fCentCL0), + fCentCL1(ana.fCentCL1), + fCentV0MvsFMD(ana.fCentV0MvsFMD), + fCentTKLvsV0M(ana.fCentTKLvsV0M), + fCentZEMvsZDC(ana.fCentZEMvsZDC), + fHtempV0M(ana.fHtempV0M), + fHtempFMD(ana.fHtempFMD), + fHtempTRK(ana.fHtempTRK), + fHtempTKL(ana.fHtempTKL), + fHtempCL0(ana.fHtempCL0), + fHtempCL1(ana.fHtempCL1), + fHtempV0MvsFMD(ana.fHtempV0MvsFMD), + fHtempTKLvsV0M(ana.fHtempTKLvsV0M), + fHtempZEMvsZDC(ana.fHtempZEMvsZDC), + fOutputList(ana.fOutputList), + fHOutCentV0M (ana.fHOutCentV0M ), + fHOutCentFMD (ana.fHOutCentFMD ), + fHOutCentTRK (ana.fHOutCentTRK ), + fHOutCentTKL (ana.fHOutCentTKL ), + fHOutCentCL0 (ana.fHOutCentCL0 ), + fHOutCentCL1 (ana.fHOutCentCL1 ), + fHOutCentV0MvsFMD(ana.fHOutCentV0MvsFMD), + fHOutCentTKLvsV0M(ana.fHOutCentTKLvsV0M), + fHOutCentZEMvsZDC(ana.fHOutCentZEMvsZDC), + fHOutMultV0M(ana.fHOutMultV0M), + fHOutMultV0R(ana.fHOutMultV0R), + fHOutMultFMD(ana.fHOutMultFMD), + fHOutMultTRK(ana.fHOutMultTRK), + fHOutMultTKL(ana.fHOutMultTKL), + fHOutMultCL0(ana.fHOutMultCL0), + fHOutMultCL1(ana.fHOutMultCL1), + fHOutMultV0MvsZDC(ana.fHOutMultV0MvsZDC), + fHOutMultZEMvsZDC(ana.fHOutMultZEMvsZDC), + fHOutMultV0MvsCL1(ana.fHOutMultV0MvsCL1), + fHOutMultV0MvsTRK(ana.fHOutMultV0MvsTRK), + fHOutMultTRKvsCL1(ana.fHOutMultTRKvsCL1) +{ + // Copy Constructor +} + +//________________________________________________________________________ +AliCentralitySelectionTask::~AliCentralitySelectionTask() +{ + // Destructor + if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList; + if (fTrackCuts) delete fTrackCuts; +} + +//________________________________________________________________________ +void AliCentralitySelectionTask::UserCreateOutputObjects() +{ + // Create the output containers + if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n"); + AliLog::SetClassDebugLevel("AliCentralitySelectionTask", AliLog::kInfo); + + fOutputList = new TList(); + fOutputList->SetOwner(); + fHOutCentV0M = new TH1F("fHOutCentV0M","fHOutCentV0M; Centrality V0",101,-0.5,100.5); + fHOutCentFMD = new TH1F("fHOutCentFMD","fHOutCentFMD; Centrality FMD",101,-0.5,100.5); + fHOutCentTRK = new TH1F("fHOutCentTRK","fHOutCentTRK; Centrality TPC",101,-0.5,100.5); + fHOutCentTKL = new TH1F("fHOutCentTKL","fHOutCentTKL; Centrality tracklets",101,-0.5,100.5); + fHOutCentCL0 = new TH1F("fHOutCentCL0","fHOutCentCL0; Centrality SPD inner",101,-0.5,100.5); + fHOutCentCL1 = new TH1F("fHOutCentCL1","fHOutCentCL1; Centrality SPD outer",101,-0.5,100.5); + fHOutCentV0MvsFMD= new TH1F("fHOutCentV0MvsFMD","fHOutCentV0MvsFMD; Centrality V0 vs FMD",101,-0.5,100.5); + fHOutCentTKLvsV0M= new TH1F("fHOutCentTKLvsV0M","fHOutCentTKLvsV0M; Centrality tracklets vs V0",101,-0.5,100.5); + fHOutCentZEMvsZDC= new TH1F("fHOutCentZEMvsZDC","fHOutCentZEMvsZDC; Centrality ZEM vs ZDC",101,-0.5,100.5); + + fHOutMultV0M = new TH1F("fHOutMultV0M","fHOutMultV0M; Multiplicity V0",25000,0,25000); + fHOutMultV0R = new TH1F("fHOutMultV0R","fHOutMultV0R; Multiplicity V0",25000,0,25000); + fHOutMultFMD = new TH1F("fHOutMultFMD","fHOutMultFMD; Multiplicity FMD",24000,0,24000); + fHOutMultTRK = new TH1F("fHOutMultTRK","fHOutMultTRK; Multiplicity TPC",4000,0,4000); + fHOutMultTKL = new TH1F("fHOutMultTKL","fHOutMultTKL; Multiplicity tracklets",5000,0,5000); + fHOutMultCL0 = new TH1F("fHOutMultCL0","fHOutMultCL0; Multiplicity SPD inner",7000,0,7000); + fHOutMultCL1 = new TH1F("fHOutMultCL1","fHOutMultCL1; Multiplicity SPD outer",7000,0,7000); + fHOutMultV0MvsZDC = new TH2F("fHOutMultV0MvsZDC","fHOutMultV0MvsZDC; Multiplicity V0; Energy ZDC",500,0,25000,500,0,6000); + fHOutMultZEMvsZDC = new TH2F("fHOutMultZEMvsZDC","fHOutMultZEMvsZDC; Energy ZEM; Energy ZDC",500,0,2500,500,0,6000); + fHOutMultV0MvsCL1 = new TH2F("fHOutMultV0MvsCL1","fHOutMultV0MvsCL1; Multiplicity V0; Multiplicity SPD outer",2500,0,25000,700,0,7000); + fHOutMultV0MvsTRK = new TH2F("fHOutMultV0MvsTRK","fHOutMultV0MvsTRK; Multiplicity V0; Multiplicity TPC",2500,0,25000,400,0,4000); + fHOutMultTRKvsCL1 = new TH2F("fHOutMultTRKvsCL1","fHOutMultTRKvsCL1; Multiplicity TPC; Multiplicity SPD outer",400,0,4000,700,0,7000); + + fOutputList->Add( fHOutCentV0M ); + fOutputList->Add( fHOutCentFMD ); + fOutputList->Add( fHOutCentTRK ); + fOutputList->Add( fHOutCentTKL ); + fOutputList->Add( fHOutCentCL0 ); + fOutputList->Add( fHOutCentCL1 ); + fOutputList->Add( fHOutCentV0MvsFMD); + fOutputList->Add( fHOutCentTKLvsV0M); + fOutputList->Add( fHOutCentZEMvsZDC); + fOutputList->Add( fHOutMultV0M); + fOutputList->Add( fHOutMultV0R); + fOutputList->Add( fHOutMultFMD); + fOutputList->Add( fHOutMultTRK); + fOutputList->Add( fHOutMultTKL); + fOutputList->Add( fHOutMultCL0); + fOutputList->Add( fHOutMultCL1); + fOutputList->Add( fHOutMultV0MvsZDC); + fOutputList->Add( fHOutMultZEMvsZDC); + fOutputList->Add( fHOutMultV0MvsCL1); + fOutputList->Add( fHOutMultV0MvsTRK); + fOutputList->Add( fHOutMultTRKvsCL1); + + + fTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); + + PostData(1, fOutputList); +} + +//________________________________________________________________________ +void AliCentralitySelectionTask::UserExec(Option_t */*option*/) +{ + // Execute analysis for current event: + if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n"); + + Float_t zncEnergy = 0.; // ZNC Energy + Float_t zpcEnergy = 0.; // ZPC Energy + Float_t znaEnergy = 0.; // ZNA Energy + Float_t zpaEnergy = 0.; // 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 spdCorr =0; // corrected spd2 multiplicity + + 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 + + Short_t v0Corr = 0; // corrected V0 multiplicity + Short_t v0CorrResc = 0; // corrected and rescaled V0 multiplicity + + Float_t zvtx =0; // z-vertex SPD + + AliCentrality *esdCent = 0; + + if(fAnalysisInput.CompareTo("ESD")==0){ + + AliVEvent* event = InputEvent(); + AliESDEvent* esd = dynamic_cast(event); + + if (fRunNo<=0) { + if (SetupRun(esd)<0) + AliFatal("Centrality File not available for this run"); + } + + esdCent = esd->GetCentrality(); + + // ***** V0 info + AliESDVZERO* esdV0 = esd->GetVZEROData(); + multV0A=esdV0->GetMTotV0A(); + multV0C=esdV0->GetMTotV0C(); + + float v0CorrR; + v0Corr = (Short_t)AliESDUtils::GetCorrV0(esd,v0CorrR); + v0CorrResc = (Short_t)v0CorrR; + + // ***** Vertex Info + const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD(); + zvtx = vtxESD->GetZ(); + + // ***** CB info (tracklets, clusters, chips) + //nTracks = event->GetNumberOfTracks(); + nTracks = fTrackCuts ? (Short_t)fTrackCuts->GetReferenceMultiplicity(esd,kTRUE):-1; + + 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); + } + + spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],zvtx); + + // ***** 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())/8.; + zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy())/8.; + znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy())/8.; + zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy())/8.; + zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0))/8.; + zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1))/8.; + + } + else if(fAnalysisInput.CompareTo("AOD")==0){ + //AliAODEvent *aod = dynamic_cast (InputEvent()); + // to be implemented + printf(" AOD analysis not yet implemented!!!\n\n"); + return; + } + + // ***** Centrality Selection + if(fHtempV0M) fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((v0Corr))); + /// else printf(" Centrality by V0 not available!!!\n\n"); + if(fHtempFMD) fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC))); + // else printf(" Centrality by FMD not available!!!\n\n"); + if(fHtempTRK) fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks)); + // else printf(" Centrality by TRK not available!!!\n\n"); + if(fHtempTKL) fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets)); + // else printf(" Centrality by TKL not available!!!\n\n"); + if(fHtempCL0) fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0])); + // else printf(" Centrality by CL0 not available!!!\n\n"); + if(fHtempCL1) fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(spdCorr)); + /// else printf(" Centrality by CL1 not available!!!\n\n"); + + if(fHtempV0MvsFMD) fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C))); + // else printf(" Centrality by V0 vs FMD not available!!!\n\n"); + if(fHtempTKLvsV0M) fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets)); + // else printf(" Centrality by V0 vs TKL not available!!!\n\n"); + if(fHtempZEMvsZDC) fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.)); + // else printf(" Centrality by ZEM vs ZDC not available!!!\n\n"); + + esdCent->SetCentralityV0M(fCentV0M); + esdCent->SetCentralityFMD(fCentFMD); + esdCent->SetCentralityTRK(fCentTRK); + esdCent->SetCentralityTKL(fCentTKL); + esdCent->SetCentralityCL0(fCentCL0); + esdCent->SetCentralityCL1(fCentCL1); + esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD); + esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M); + esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC); + + fHOutCentV0M->Fill(fCentV0M); + fHOutCentFMD->Fill(fCentFMD); + fHOutCentTRK->Fill(fCentTRK); + fHOutCentTKL->Fill(fCentTKL); + fHOutCentCL0->Fill(fCentCL0); + fHOutCentCL1->Fill(fCentCL1); + fHOutCentV0MvsFMD->Fill(fCentV0MvsFMD); + fHOutCentTKLvsV0M->Fill(fCentTKLvsV0M); + fHOutCentZEMvsZDC->Fill(fCentZEMvsZDC); + fHOutMultV0M->Fill(v0Corr); + fHOutMultV0R->Fill(multV0A+multV0C); + fHOutMultFMD->Fill((multFMDA+multFMDC)); + fHOutMultTRK->Fill(nTracks); + fHOutMultTKL->Fill(nTracklets); + fHOutMultCL0->Fill(nClusters[0]); + fHOutMultCL1->Fill(spdCorr); + fHOutMultV0MvsZDC->Fill(v0Corr,(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy)); + fHOutMultZEMvsZDC->Fill((zem1Energy+zem2Energy),(zncEnergy+znaEnergy+zpcEnergy+zpaEnergy)); + fHOutMultV0MvsCL1->Fill(v0Corr,spdCorr); + fHOutMultV0MvsTRK->Fill(v0Corr,nTracks); + fHOutMultTRKvsCL1->Fill(nTracks,spdCorr); + + PostData(1, fOutputList); +} + +//________________________________________________________________________ +void AliCentralitySelectionTask::ReadCentralityHistos(TString fCentfilename) +{ + // Read centrality histograms + TDirectory *owd = gDirectory; + // Check if the file is present + TString path = gSystem->ExpandPathName(fCentfilename.Data()); + if (gSystem->AccessPathName(path)) { + AliError(Form("File %s does not exist", path.Data())); + return; + } + fFile = TFile::Open(fCentfilename); + owd->cd(); + fHtempV0M = (TH1F*) (fFile->Get("hmultV0_percentile")); + fHtempFMD = (TH1F*) (fFile->Get("hmultFMD_percentile")); + fHtempTRK = (TH1F*) (fFile->Get("hNtracks_percentile")); + fHtempTKL = (TH1F*) (fFile->Get("hNtracklets_percentile")); + fHtempCL0 = (TH1F*) (fFile->Get("hNclusters0_percentile")); + fHtempCL1 = (TH1F*) (fFile->Get("hNclusters1_percentile")); + owd->cd(); +} + +//________________________________________________________________________ +void AliCentralitySelectionTask::ReadCentralityHistos2(TString fCentfilename2) +{ + // Read centrality histograms + TDirectory *owd = gDirectory; + TString path = gSystem->ExpandPathName(fCentfilename2.Data()); + if (gSystem->AccessPathName(path)) { + AliError(Form("File %s does not exist", path.Data())); + return; + } + fFile2 = TFile::Open(fCentfilename2); + owd->cd(); + fHtempV0MvsFMD = (TH1F*) (fFile2->Get("hmultV0vsmultFMD_all_percentile")); + fHtempTKLvsV0M = (TH1F*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile")); + fHtempZEMvsZDC = (TH1F*) (fFile2->Get("hEzemvsEzdc_all_percentile")); + owd->cd(); +} + +//________________________________________________________________________ +void AliCentralitySelectionTask::Terminate(Option_t */*option*/) +{ + // Terminate analysis + if (fFile && fFile->IsOpen()) + fFile->Close(); + if (fFile2 && fFile2->IsOpen()) + fFile2->Close(); +} +//________________________________________________________________________ +Int_t AliCentralitySelectionTask::SetupRun(AliESDEvent* esd) +{ + // Setup files for run + + if (!esd) + return -1; + + // check if something to be done + if (fCurrentRun == esd->GetRunNumber()) + return 0; + else + fCurrentRun = esd->GetRunNumber(); + + AliInfo(Form("Setup Centrality Selection for run %d\n",fCurrentRun)); + + // CHANGE HERE FOR RUN RANGES + if ( fCurrentRun <= 137165 ) fRunNo = 137161; + else fRunNo = 137366; + // CHANGE HERE FOR RUN RANGES + + TString fileName(Form("%s/AliCentralityBy1D_%d.root", fDirName.Data(), fRunNo)); + TString fileName2(Form("%s/AliCentralityByFunction_%d.root", fDirName.Data(), fRunNo)); + + AliInfo(Form("Centrality Selection for run %d is initialized with %s", fCurrentRun, fileName.Data())); + ReadCentralityHistos(fileName.Data()); + ReadCentralityHistos2(fileName2.Data()); + if (!fFile && !fFile2) { + AliFatal(Form("Run %d not known to centrality selection!", fCurrentRun)); + return -1; + } + return 0; +}