[u/mrichter/AliRoot.git] / ANALYSIS / AliCentralitySelectionTask.cxx

/**************************************************************************
 * 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),
  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)
{   
  // 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),
  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)
{
  // 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),
  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)
{
  // 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]));
  fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(nClusters[1]));
  
  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->SetCentralityCL1(fCentCL1);
  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")); 
  fHtempCL1  = (TH1D*) (fFile->Get("hNclusters1_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();  
}
//________________________________________________________________________
Commit	Line	Data
e8472cd0	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 <TTree.h>
	23	#include <TList.h>
	24	#include <TH1F.h>
	25	#include <TH2F.h>
	26	#include <TProfile.h>
	27	#include <TFile.h>
	28	#include <TString.h>
	29	#include <TCanvas.h>
	30	#include <iostream>
	31
	32	#include "AliAnalysisManager.h"
	33	#include "AliVEvent.h"
	34	#include "AliESD.h"
	35	#include "AliESDEvent.h"
	36	#include "AliESDHeader.h"
	37	#include "AliESDInputHandler.h"
	38	#include "AliESDZDC.h"
	39	#include "AliESDFMD.h"
	40	#include "AliESDVZERO.h"
	41	#include "AliESDCentrality.h"
	42	#include "AliMultiplicity.h"
	43	#include "AliAODHandler.h"
	44	#include "AliAODEvent.h"
	45	#include "AliAODVertex.h"
	46	#include "AliAODMCHeader.h"
	47	#include "AliMCEvent.h"
	48	#include "AliMCEventHandler.h"
	49	#include "AliMCParticle.h"
	50	#include "AliStack.h"
	51	#include "AliHeader.h"
	52	#include "AliAODMCParticle.h"
	53	#include "AliAnalysisTaskSE.h"
	54	#include "AliGenEventHeader.h"
	55	#include "AliGenHijingEventHeader.h"
	56	#include "AliPhysicsSelectionTask.h"
	57	#include "AliPhysicsSelection.h"
	58	#include "AliBackgroundSelection.h"
	59	#include "AliCentralitySelectionTask.h"
	60
	61	ClassImp(AliCentralitySelectionTask)
	62
	63
	64	//________________________________________________________________________
65	AliCentralitySelectionTask::AliCentralitySelectionTask():
66	AliAnalysisTaskSE(),
67	fDebug(0),
68	fAnalysisInput("ESD"),
69	fIsMCInput(kFALSE),
70	fFile(0),
d15bf53f	71	fFile2(0),
e8472cd0	72	fCentfilename(0),
d15bf53f	73	fCentfilename2(0),
	74	fCentV0M(0),
	75	fCentFMD(0),
	76	fCentTRK(0),
	77	fCentTKL(0),
	78	fCentCL0(0),
be0d4e9b	79	fCentCL1(0),
d15bf53f	80	fCentV0MvsFMD(0),
	81	fCentTKLvsV0M(0),
	82	fCentZEMvsZDC(0),
	83	fHtempV0M(0),
	84	fHtempFMD(0),
	85	fHtempTRK(0),
	86	fHtempTKL(0),
	87	fHtempCL0(0),
be0d4e9b	88	fHtempCL1(0),
d15bf53f	89	fHtempV0MvsFMD(0),
	90	fHtempTKLvsV0M(0),
	91	fHtempZEMvsZDC(0)
e8472cd0	92	{
	93	// Default constructor
	94	}
	95
	96	//________________________________________________________________________
	97	AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
	98	AliAnalysisTaskSE(name),
	99	fDebug(0),
	100	fAnalysisInput("ESD"),
	101	fIsMCInput(kFALSE),
	102	fFile(0),
d15bf53f	103	fFile2(0),
e8472cd0	104	fCentfilename(0),
d15bf53f	105	fCentfilename2(0),
	106	fCentV0M(0),
	107	fCentFMD(0),
	108	fCentTRK(0),
	109	fCentTKL(0),
	110	fCentCL0(0),
be0d4e9b	111	fCentCL1(0),
d15bf53f	112	fCentV0MvsFMD(0),
	113	fCentTKLvsV0M(0),
	114	fCentZEMvsZDC(0),
	115	fHtempV0M(0),
	116	fHtempFMD(0),
	117	fHtempTRK(0),
	118	fHtempTKL(0),
	119	fHtempCL0(0),
be0d4e9b	120	fHtempCL1(0),
d15bf53f	121	fHtempV0MvsFMD(0),
	122	fHtempTKLvsV0M(0),
	123	fHtempZEMvsZDC(0)
e8472cd0	124	{
	125	// Default constructor
	126	}
	127
	128	//________________________________________________________________________
	129	AliCentralitySelectionTask& AliCentralitySelectionTask::operator=(const AliCentralitySelectionTask& c)
	130	{
	131	// Assignment operator
	132	if (this!=&c) {
	133	AliAnalysisTaskSE::operator=(c);
	134	}
	135	return *this;
	136	}
	137
	138	//________________________________________________________________________
	139	AliCentralitySelectionTask::AliCentralitySelectionTask(const AliCentralitySelectionTask& ana):
	140	AliAnalysisTaskSE(ana),
	141	fDebug(ana.fDebug),
	142	fAnalysisInput(ana.fDebug),
	143	fIsMCInput(ana.fIsMCInput),
	144	fFile(ana.fFile),
d15bf53f	145	fFile2(ana.fFile2),
e8472cd0	146	fCentfilename(ana.fCentfilename),
d15bf53f	147	fCentfilename2(ana.fCentfilename2),
	148	fCentV0M(ana.fCentV0M),
	149	fCentFMD(ana.fCentFMD),
	150	fCentTRK(ana.fCentTRK),
	151	fCentTKL(ana.fCentTKL),
	152	fCentCL0(ana.fCentCL0),
be0d4e9b	153	fCentCL1(ana.fCentCL1),
d15bf53f	154	fCentV0MvsFMD(ana.fCentV0MvsFMD),
	155	fCentTKLvsV0M(ana.fCentTKLvsV0M),
	156	fCentZEMvsZDC(ana.fCentZEMvsZDC),
	157	fHtempV0M(ana.fHtempV0M),
	158	fHtempFMD(ana.fHtempFMD),
	159	fHtempTRK(ana.fHtempTRK),
	160	fHtempTKL(ana.fHtempTKL),
	161	fHtempCL0(ana.fHtempCL0),
be0d4e9b	162	fHtempCL1(ana.fHtempCL1),
d15bf53f	163	fHtempV0MvsFMD(ana.fHtempV0MvsFMD),
	164	fHtempTKLvsV0M(ana.fHtempTKLvsV0M),
	165	fHtempZEMvsZDC(ana.fHtempZEMvsZDC)
e8472cd0	166	{
	167	// Copy Constructor
	168	}
	169
	170	//________________________________________________________________________
	171	AliCentralitySelectionTask::~AliCentralitySelectionTask()
	172	{
	173	// Destructor
	174	}
	175
	176	//________________________________________________________________________
	177	void AliCentralitySelectionTask::UserCreateOutputObjects()
	178	{
	179	// Create the output containers
	180	if(fDebug>1) printf("AnalysisCentralitySelectionTask::UserCreateOutputObjects() \n");
	181	}
	182
	183	//________________________________________________________________________
	184	void AliCentralitySelectionTask::UserExec(Option_t /option*/)
	185	{
	186	// Execute analysis for current event:
	187	if(fDebug>1) printf(" **** AliCentralitySelectionTask::UserExec() \n");
	188
be0d4e9b	189	Float_t zncEnergy; // ZNC Energy
	190	Float_t zpcEnergy; // ZPC Energy
	191	Float_t znaEnergy; // ZNA Energy
	192	Float_t zpaEnergy; // ZPA Energy
d15bf53f	193	Float_t zem1Energy = 0.; // ZEM1 Energy
d15bf53f	194	Float_t zem2Energy = 0.; // ZEM2 Energy
e8472cd0	195
be0d4e9b	196	Int_t nTracks = 0; // no. tracks
d15bf53f	197	Int_t nTracklets = 0; // no. tracklets
	198	Int_t nClusters[6]; // no. clusters on 6 ITS layers
	199	Int_t nChips[2]; // no. chips on 2 SPD layers
e8472cd0	200
d15bf53f	201	Float_t multV0A = 0; // multiplicity from V0 reco side A
	202	Float_t multV0C = 0; // multiplicity from V0 reco side C
	203	Float_t multFMDA = 0; // multiplicity from FMD on detector A
	204	Float_t multFMDC = 0; // multiplicity from FMD on detector C
e8472cd0	205
	206	AliESDCentrality *esdCent = 0;
	207
	208	if(fAnalysisInput.CompareTo("ESD")==0){
	209
	210	AliVEvent* event = InputEvent();
	211	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
d15bf53f	212
e8472cd0	213	esdCent = esd->GetCentrality();
	214
	215	// ***** V0 info
	216	AliESDVZERO* esdV0 = esd->GetVZEROData();
	217	multV0A=esdV0->GetMTotV0A();
	218	multV0C=esdV0->GetMTotV0C();
	219
	220	// ***** CB info (tracklets, clusters, chips)
	221	nTracks = event->GetNumberOfTracks();
	222
	223	const AliMultiplicity *mult = esd->GetMultiplicity();
	224
	225	nTracklets = mult->GetNumberOfTracklets();
	226
	227	for(Int_t ilay=0; ilay<6; ilay++){
	228	nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
	229	}
	230
	231	for(Int_t ilay=0; ilay<2; ilay++){
	232	nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
	233	}
	234
	235
	236	// ***** FMD info
	237	AliESDFMD *fmd = esd->GetFMDData();
	238	Float_t totalMultA = 0;
	239	Float_t totalMultC = 0;
	240	const Float_t fFMDLowCut = 0.4;
	241
	242	for(UShort_t det=1;det<=3;det++) {
	243	Int_t nRings = (det==1 ? 1 : 2);
	244	for (UShort_t ir = 0; ir < nRings; ir++) {
	245	Char_t ring = (ir == 0 ? 'I' : 'O');
	246	UShort_t nsec = (ir == 0 ? 20 : 40);
	247	UShort_t nstr = (ir == 0 ? 512 : 256);
	248	for(UShort_t sec =0; sec < nsec; sec++) {
	249	for(UShort_t strip = 0; strip < nstr; strip++) {
	250
	251	Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
	252	if(FMDmult == 0 \|\| FMDmult == AliESDFMD::kInvalidMult) continue;
	253
	254	Float_t nParticles=0;
	255
	256	if(FMDmult > fFMDLowCut) {
	257	nParticles = 1.;
	258	}
	259
	260	if (det<3) totalMultA = totalMultA + nParticles;
	261	else totalMultC = totalMultC + nParticles;
	262
	263	}
	264	}
	265	}
	266	}
	267	multFMDA = totalMultA;
	268	multFMDC = totalMultC;
	269
	270	// ***** ZDC info
	271	AliESDZDC *esdZDC = esd->GetESDZDC();
d15bf53f	272	zncEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
	273	zpcEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
	274	znaEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
	275	zpaEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
e8472cd0	276	zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
	277	zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
	278
	279	}
	280	else if(fAnalysisInput.CompareTo("AOD")==0){
	281	//AliAODEvent aod = dynamic_cast<AliAODEvent> (InputEvent());
	282	// to be implemented
	283	printf(" AOD analysis not yet implemented!!!\n\n");
	284	return;
	285	}
	286
d15bf53f	287	// ***** Centrality Selection
	288	fCentV0M = fHtempV0M->GetBinContent(fHtempV0M->FindBin((multV0A+multV0C)));
	289	fCentFMD = fHtempFMD->GetBinContent(fHtempFMD->FindBin((multFMDA+multFMDC)));
	290	fCentTRK = fHtempTRK->GetBinContent(fHtempTRK->FindBin(nTracks));
	291	fCentTKL = fHtempTKL->GetBinContent(fHtempTKL->FindBin(nTracklets));
	292	fCentCL0 = fHtempCL0->GetBinContent(fHtempCL0->FindBin(nClusters[0]));
be0d4e9b	293	fCentCL1 = fHtempCL1->GetBinContent(fHtempCL1->FindBin(nClusters[1]));
d15bf53f	294
	295	fCentV0MvsFMD = fHtempV0MvsFMD->GetBinContent(fHtempV0MvsFMD->FindBin((multV0A+multV0C)));
	296	fCentTKLvsV0M = fHtempTKLvsV0M->GetBinContent(fHtempTKLvsV0M->FindBin(nTracklets));
	297	fCentZEMvsZDC = fHtempZEMvsZDC->GetBinContent(fHtempZEMvsZDC->FindBin((zem1Energy+zem2Energy)/1000.));
e8472cd0	298
d15bf53f	299	esdCent->SetCentralityV0M(fCentV0M);
	300	esdCent->SetCentralityFMD(fCentFMD);
	301	esdCent->SetCentralityTRK(fCentTRK);
	302	esdCent->SetCentralityTKL(fCentTKL);
	303	esdCent->SetCentralityCL0(fCentCL0);
be0d4e9b	304	esdCent->SetCentralityCL1(fCentCL1);
d15bf53f	305	esdCent->SetCentralityV0MvsFMD(fCentV0MvsFMD);
	306	esdCent->SetCentralityTKLvsV0M(fCentTKLvsV0M);
	307	esdCent->SetCentralityZEMvsZDC(fCentZEMvsZDC);
e8472cd0	308	}
	309
	310	//________________________________________________________________________
d15bf53f	311	void AliCentralitySelectionTask::ReadCentralityHistos()
e8472cd0	312	{
e8472cd0	313	fFile = new TFile(fCentfilename);
d15bf53f	314	fHtempV0M = (TH1D*) (fFile->Get("hmultV0_percentile"));
	315	fHtempFMD = (TH1D*) (fFile->Get("hmultFMD_percentile"));
	316	fHtempTRK = (TH1D*) (fFile->Get("hNtracks_percentile"));
	317	fHtempTKL = (TH1D*) (fFile->Get("hNtracklets_percentile"));
	318	fHtempCL0 = (TH1D*) (fFile->Get("hNclusters0_percentile"));
be0d4e9b	319	fHtempCL1 = (TH1D*) (fFile->Get("hNclusters1_percentile"));
d15bf53f	320	}
e8472cd0	321
d15bf53f	322	void AliCentralitySelectionTask::ReadCentralityHistos2()
	323	{
	324	fFile2 = new TFile(fCentfilename2);
	325	fHtempV0MvsFMD = (TH1D*) (fFile2->Get("hmultV0vsmultFMD_all_percentile"));
	326	fHtempTKLvsV0M = (TH1D*) (fFile2->Get("hNtrackletsvsmultV0_all_percentile"));
	327	fHtempZEMvsZDC = (TH1D*) (fFile2->Get("hEzemvsEzdc_all_percentile"));
	328	}
	329
	330	void AliCentralitySelectionTask::SetPercentileFile(TString filename)
	331	{
	332	fCentfilename = filename;
	333	ReadCentralityHistos();
	334	}
	335
	336	void AliCentralitySelectionTask::SetPercentileFile2(TString filename)
	337	{
	338	fCentfilename2 = filename;
	339	ReadCentralityHistos2();
e8472cd0	340	}
	341
	342	//________________________________________________________________________
	343	void AliCentralitySelectionTask::Terminate(Option_t /option*/)
	344	{
	345	// Terminate analysis
	346	fFile->Close();
d15bf53f	347	fFile2->Close();
e8472cd0	348	}
	349	//________________________________________________________________________
	350