/**************************************************************************
 * 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),
  fCentfilename(0),
  fMethod(0),
  fCent(0),
  fHtemp(0)
{   
  // Default constructor
}   

//________________________________________________________________________
AliCentralitySelectionTask::AliCentralitySelectionTask(const char *name):
  AliAnalysisTaskSE(name),
  fDebug(0),
  fAnalysisInput("ESD"),
  fIsMCInput(kFALSE),
  fFile(0),
  fCentfilename(0),
  fMethod(0),
  fCent(0),
  fHtemp(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),
  fCentfilename(ana.fCentfilename),
  fMethod(ana.fMethod),
  fCent(ana.fCent),
  fHtemp(ana.fHtemp)
{
  // 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;                //  ZEM1 Energy
  Float_t  zem2Energy;                //  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
  if(fMethod.CompareTo("V0")==0){
    fCent = fHtemp->GetBinContent(fHtemp->FindBin((multV0A+multV0C)));
  }
  if(fMethod.CompareTo("FMD")==0){
    fCent = fHtemp->GetBinContent(fHtemp->FindBin((multFMDA+multFMDC)));
  }
  if(fMethod.CompareTo("TRACKS")==0) {
    fCent = fHtemp->GetBinContent(fHtemp->FindBin(nTracks));
  }
  if(fMethod.CompareTo("TRACKLETS")==0){
    fCent = fHtemp->GetBinContent(fHtemp->FindBin(nTracklets));
  }
  if(fMethod.CompareTo("CLUSTERS")==0) {
    fCent = fHtemp->GetBinContent(fHtemp->FindBin(nClusters[0]));
  }
  printf(" **** centrality is %3.2f \n", fCent);
  
  esdCent->SetCentrality(fCent);
}

//________________________________________________________________________
void AliCentralitySelectionTask::SetCentralityMethod(const char* x) 
{
  fMethod = x;

  fFile  = new TFile(fCentfilename);
  
  if(fMethod.CompareTo("V0")==0)
    fHtemp  = (TH1D*) (fFile->Get("hmultV0_percentile")); 
  
  if(fMethod.CompareTo("FMD")==0)
    fHtemp  = (TH1D*) (fFile->Get("hmultFMD_percentile")); 
  
  if(fMethod.CompareTo("TRACKS")==0) 
    fHtemp  = (TH1D*) (fFile->Get("hNtracks_percentile")); 
  
  if(fMethod.CompareTo("TRACKLETS")==0)
    fHtemp  = (TH1D*) (fFile->Get("hNtracklets_percentile")); 
  
  if(fMethod.CompareTo("CLUSTERS")==0) 
    fHtemp  = (TH1D*) (fFile->Get("hNclusters0_percentile")); 
  
}

//________________________________________________________________________
void AliCentralitySelectionTask::Terminate(Option_t */*option*/)
{
  // Terminate analysis
  fFile->Close();  
}
//________________________________________________________________________