/*************************************************************************
* 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.                  *
**************************************************************************/

//  ******************************************
//  * analysis task for azimuthal isotropic  *
//  * expansion in highly central collisions *
//  * author: Cristian Andrei                *
//  *         acristian@niham.nipne.ro       *
//  ******************************************

#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TList.h"
#include "TObjArray.h"
#include "TString.h"
// #include "TCanvas.h"

#include "AliAnalysisManager.h"
#include "AliMCEvent.h"
#include "AliMCEventHandler.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliESDtrackCuts.h"
#include "AliCFContainer.h"

#include "AliAnalysisCentralCutMC.h"
#include "AliAnalysisCentralCutESD.h"
#include "AliAnalysisCentralCutEvtMC.h"
#include "AliAnalysisCentralCutEvtESD.h"
#include "AliAnalysisCentralExtrapolate.h"
#include "AliAnalysisTaskCentral.h"


ClassImp(AliAnalysisTaskCentral)


//________________________________________________________________________
AliAnalysisTaskCentral::AliAnalysisTaskCentral(const char *name) 
  : AliAnalysisTask(name, "")
  ,fESD(0), fMC(0)
  ,fNoEvt(0)
  ,fCFContainerPi(0)
  ,fCFContainerK(0)
  ,fCFContainerP(0)
  ,fSim(kFALSE)
  ,fOutList(NULL)

{
  // Constructor
	printf("AliAnalysisTaskCentral::AliAnalysisTaskCentral(const char *name)\n");

  // Define input and output slots here
	DefineInput(0, TChain::Class());

	DefineOutput(0, TList::Class()); 


	for(Int_t i=0; i<10; i++){
		fCutsList[i] = 0;
    } 

}


//________________________________________________________________________
AliAnalysisTaskCentral::~AliAnalysisTaskCentral() 
{
// Destructor
// Delete the created objects

	if(fESD) delete fESD;
	if(fMC) delete fMC;
	
	if(fCutsList){ 
    	delete [] fCutsList;
  	}
	
	if(fNoEvt) delete fNoEvt;
	
	if(fCFContainerPi) delete fCFContainerPi;
	if(fCFContainerPi) delete fCFContainerK;
	if(fCFContainerPi) delete fCFContainerP;
	
	if(fOutList) delete fOutList;

}

//________________________________________________________________________
void AliAnalysisTaskCentral::ConnectInputData(Option_t *) {
// get the event from the input chain

    TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
	printf("tree(%p)\n", (void*)tree);

    if (!tree) {
	Printf("ERROR: Could not read chain from input slot 0");
    } else {
	AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	if (!esdH) {
    	    Printf("ERROR: Could not get ESDInputHandler");
	} 
	else fESD = esdH->GetEvent();
  
	AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
	if (!mcH) {
	    Printf("ERROR: Could not get MCInputHandler");
	}
	else fMC = mcH->MCEvent();
    }	
}

//________________________________________________________________________
void AliAnalysisTaskCentral::CreateOutputObjects(){
//Creates the output list

	fNoEvt = new TH1F(Form("%s_NoEvt",GetName()), "Number of processed events", 1, 0.0, 1.0);

//set here the:     
    //min and max variable values (to do: move these to a separate cuts class)
    const Double_t ptmin  = 0.0 ; //bins in pt
    const Double_t ptmax  = 5.0 ; //bins in y 
    const Double_t etamin  = -0.9 ; //bins in pt
    const Double_t etamax  = 0.9 ; //bins in y 
    //apropiate number of bins for histos and CF Container
    const Int_t nbinpt  = 25 ; //bins in pt
    const Int_t nbineta  = 25 ; //bins in y 

//Correction Framework CONTAINER DEFINITION
    //the sensitive variables, their indices
	UInt_t ipt = 0;
	UInt_t ieta  = 1;
    //Setting up the container grid... 
    const UInt_t nstep = 2 ; //number of selection steps MC & ESD 
    const Int_t nvar = 2 ; //number of variables on the grid:pt,eta
	
    //arrays for the number of bins in each dimension
    //nbin is defined above =
    Int_t iBin[nvar];
    iBin[0]=nbinpt;
    iBin[1]=nbineta;

    //arrays for lower bounds :
    Double_t *binLim1=new Double_t[nbinpt+1];
    Double_t *binLim2=new Double_t[nbineta+1];

    //values for bin lower bounds
    for(Int_t i=0; i<=nbinpt; i++) binLim1[i]=(Double_t)ptmin + (ptmax-ptmin)/nbinpt*(Double_t)i ; 
    for(Int_t i=0; i<=nbineta; i++) binLim2[i]=(Double_t)etamin  + (etamax-etamin)  /nbineta*(Double_t)i ;

	//container creation
	fCFContainerPi = new AliCFContainer(Form("%s_CFCont_Pi",GetName()),"container for tracks",nstep,nvar,iBin);
	fCFContainerK = new AliCFContainer(Form("%s_CFCont_K",GetName()),"container for tracks",nstep,nvar,iBin);
	fCFContainerP = new AliCFContainer(Form("%s_CFCont_P",GetName()),"container for tracks",nstep,nvar,iBin);
	
    //setting the bin limits
    fCFContainerPi -> SetBinLimits(ipt,binLim1);
    fCFContainerPi -> SetBinLimits(ieta,binLim2);
	
	fCFContainerK -> SetBinLimits(ipt,binLim1);
	fCFContainerK -> SetBinLimits(ieta,binLim2);
	
	fCFContainerP -> SetBinLimits(ipt,binLim1);
	fCFContainerP -> SetBinLimits(ieta,binLim2);
	
//outout list creation
	fOutList = new TList();
	fOutList->Add(fNoEvt);
	fOutList->Add(fCFContainerPi);
	fOutList->Add(fCFContainerK);
	fOutList->Add(fCFContainerP);
}


//________________________________________________________________
void AliAnalysisTaskCentral::InitCuts(){

//Create and set cuts

//------------------EVENT LEVEL CUTS--------------------
//-----------------------MC-----------------------------
	AliAnalysisCentralCutEvtMC *evMCCuts = new AliAnalysisCentralCutEvtMC();
	evMCCuts->SetMultiplicityRange(1,100);
	evMCCuts->SetDirectivityRange(0.0, 1.0);
    evMCCuts->SetDirUnitRange(0.0, 1.0);

//-----------------------ESD----------------------------
	AliAnalysisCentralCutEvtESD *evESDCuts = new AliAnalysisCentralCutEvtESD();
	evESDCuts->SetMultiplicityRange(1,100);
	evESDCuts->SetDirectivityRange(0.0, 1.0);
    evESDCuts->SetSPDMultiplicityRange(1,20000);
    evESDCuts->SetSPDDirectivityRange(0.0, 1.0);

	TObjArray* mcEventCuts = new TObjArray();
	mcEventCuts->AddLast(evMCCuts);

	TObjArray* esdEventCuts = new TObjArray();
	esdEventCuts->AddLast(evESDCuts);


//------------------PARTICLE LEVEL CUTS-----------------
//-------------------General MC Cuts--------------------
	AliAnalysisCentralCutMC *mcCutsGen = new AliAnalysisCentralCutMC();
	mcCutsGen->SetOnlyPrimaries(kTRUE);
	mcCutsGen->SetPtRange(0.2,4.0);
	mcCutsGen->SetEtaRange(-0.2,0.2);

//-------------------Specific MC Cuts-------------------
	AliAnalysisCentralCutMC *mcCutsPi = new AliAnalysisCentralCutMC();
	mcCutsPi->SetPDGCode(211); //211 pion;  321 kaon; 2212 proton

	AliAnalysisCentralCutMC *mcCutsK = new AliAnalysisCentralCutMC();
	mcCutsK->SetPDGCode(321); //211 pion;  321 kaon; 2212 proton

	AliAnalysisCentralCutMC *mcCutsP = new AliAnalysisCentralCutMC();
	mcCutsP->SetPDGCode(2212); //211 pion;  321 kaon; 2212 proton

//--------each task has its own mcList of cuts----------
	TObjArray *mcListGen = new TObjArray(); //general MC cuts
	mcListGen->AddLast(mcCutsGen);
	
	TObjArray *mcListPi = new TObjArray(); //task pt pions
	mcListPi->AddLast(mcCutsPi);

	TObjArray *mcListK = new TObjArray(); //task pt kaons
	mcListK->AddLast(mcCutsK);

	TObjArray *mcListP = new TObjArray(); //task pt protons
	mcListP->AddLast(mcCutsP);


//-------------------General ESD Cuts-------------------
	AliESDtrackCuts *esdCutsGen = new AliESDtrackCuts("AliESDtrackCuts", "Loose");
	esdCutsGen->SetMinNClustersTPC(50);
	esdCutsGen->SetMaxChi2PerClusterTPC(2.2);
	esdCutsGen->SetMaxCovDiagonalElements(0.5,0.5,0.5,0.5,0.5);
	esdCutsGen->SetRequireTPCRefit(kTRUE);
	esdCutsGen->SetAcceptKinkDaughters(kFALSE);
	esdCutsGen->SetMaxNsigmaToVertex(2.0);
	esdCutsGen->SetRequireSigmaToVertex(kTRUE);
	esdCutsGen->SetPtRange(0.2,4.0);
	esdCutsGen->SetEtaRange(-0.2,0.2);

//-------------------Specific ESD Cuts------------------
	AliAnalysisCentralCutESD *esdCutsPi = new AliAnalysisCentralCutESD("AliAnalysisCentralCutESD","NIHAM");
	esdCutsPi->SetPIDtype("Custom");
	esdCutsPi->SetPriorFunctions(kFALSE);
	esdCutsPi->SetPartType(kPiPlus);

	AliAnalysisCentralCutESD *esdCutsK = new AliAnalysisCentralCutESD("AliAnalysisCentralCutESD","NIHAM");
	esdCutsK->SetPIDtype("Custom");
	esdCutsK->SetPriorFunctions(kFALSE);
	esdCutsK->SetPartType(kKPlus);

	AliAnalysisCentralCutESD *esdCutsP = new AliAnalysisCentralCutESD("AliAnalysisCentralCutESD","NIHAM");
	esdCutsP->SetPIDtype("Custom");
	esdCutsP->SetPriorFunctions(kFALSE);
	esdCutsP->SetPartType(kProton);

//--------each task has its own esdList of cuts---------
	TObjArray* esdListGen = new TObjArray(); //general ESD track cuts
	esdListGen->AddLast(esdCutsGen);
	
	TObjArray* esdListPi = new TObjArray();
	esdListPi->AddLast(esdCutsPi);

	TObjArray* esdListK = new TObjArray();
	esdListK->AddLast(esdCutsK);

	TObjArray* esdListP = new TObjArray();
	esdListP->AddLast(esdCutsP);

//------set the cuts to the RIGHT! fCutsList slots-------
// // event level cuts
	SetCuts(0, mcEventCuts);
	SetCuts(1, esdEventCuts);
 
 //particle level cuts
	SetCuts(2, mcListGen);
	SetCuts(3, mcListPi);
	SetCuts(4, mcListK);
	SetCuts(5, mcListP);
 
	SetCuts(6, esdListGen);
	SetCuts(7, esdListPi);
	SetCuts(8, esdListK);
	SetCuts(9, esdListP);

}


//_______________________________________________________________________
void AliAnalysisTaskCentral::SendEvent(TObject *obj) const{

// Some cuts (ie MC IsPrimary) need the MC Event info

    if (!fCutsList) {
		printf("No particle cut list found!\n\n");
		return;
    }
    else {
		for(Int_t isel=0;isel< 10; isel++){
			if(!fCutsList[isel]) continue;
			TObjArrayIter iter(fCutsList[isel]);
			AliAnalysisCuts *cut = 0;
	
			while ((cut = (AliAnalysisCuts*)iter.Next())) {
				
				TString cutName=cut->GetName();
				if(!cutName){
					printf("No cutname!\n");
					return;
				}
		
				Bool_t checkCut=cutName.Contains("AliAnalysisCentralCutMC");
					
				if(checkCut){
					AliAnalysisCentralCutMC *newcut = dynamic_cast<AliAnalysisCentralCutMC *>(cut);
					newcut->ReceiveEvt(obj);
				}
			}
		}
    }

}


//________________________________________________________________________
Bool_t AliAnalysisTaskCentral::CheckCuts(Int_t no, TObject *obj) const{ 

// For each cut run IsSelected();

    if(no > 10){
		printf("\nAliAnalysisTaskCentral::CheckCuts -> Cut number is not ok! \n");
		return kFALSE;
    }

    if(!fCutsList[no]){
		printf("AliAnalysisTaskCentral::CheckCuts -> cuts list problem! \n");
		return kFALSE;
    }

    TObjArrayIter iter(fCutsList[no]);
    AliAnalysisCuts *cut = 0;

    while((cut = (AliAnalysisCuts*)iter.Next())){
		if(!cut->IsSelected(obj)) return kFALSE;
    }

    return kTRUE;
}


//________________________________________________________________________
void AliAnalysisTaskCentral::Exec(Option_t *) {

// Main loop
// Called for each event
	
	InitCuts(); //initialize the analysis specific cuts	
    if (!fCutsList) {
		Printf("ERROR: fCutsList not available");
		return;
    }

	Double_t pt, eta;
	Double_t w = 1.0;    
	const Int_t nvar=2; //number of variables on the grid:pt,vtx
	Double_t value[nvar];//used to fill the CF container
	
	if(fSim){  // if running on simulations -> look at MC Truth
		
		if (!fMC) {
			Printf("ERROR: fMC not available");
			return;
		}
    
		if(CheckCuts(0, fMC)){ //check event level cuts
	
			SendEvent(fMC);
		
			// MC loop
			for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++) {
		
				AliMCParticle *particle  = dynamic_cast<AliMCParticle*>(fMC->GetTrack(ipart));
						
				if(!particle){
						printf("\nMCParticle pointer is null!!!\n");
						continue;
				}
				
				if (!CheckCuts(2, particle)) continue; //check the MC general particle cuts 
				
				pt = particle->Pt();
				eta = particle->Eta();
				
				if(pt>0) w = 1.0/pt; //invariant distribution
				
				value[0]=pt;
				value[1]=eta;
				
				if(CheckCuts(3, particle)){  //fill the right container for each particle
					fCFContainerPi->Fill(value,0,w);
				}
				else if(CheckCuts(4, particle)){
					fCFContainerK->Fill(value,0,w);
				}
				else if(CheckCuts(5, particle)){
					fCFContainerP->Fill(value,0,w);
				}
			} //end MC particle loop 
		}
	}
	else{  // if we DONT run in simulated data we fill the MC step of the CFCont with 0 
		value[0]=0;
		value[1]=0;
			
		fCFContainerPi->Fill(value,0);
		fCFContainerK->Fill(value,0);
		fCFContainerP->Fill(value,0);
	}
	
	
    if (!fESD) {
	Printf("ERROR: fESD not available");
	return;
    }

    if(CheckCuts(1, fESD)){

		Printf("There are %d ESD tracks in this event\n", fESD->GetNumberOfTracks());
	
		// ESD loop 
		for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
	
			AliESDtrack* track = fESD->GetTrack(iTracks);
			
			if (!track) {
					Printf("ERROR: Could not receive track %d", iTracks);
				continue;
			}
		
			if(!CheckCuts(6, track)) continue; //check general ESD track cuts
		
			pt = track->Pt();
			eta = track->Eta();
		
			if(pt>0) w = 1.0/pt; // invariant
			
			value[0]=pt;
			value[1]=eta;
		
			if(CheckCuts(7, track)){
				fCFContainerPi->Fill(value,1,w);
			}
			else if(CheckCuts(8, track)){
				fCFContainerK->Fill(value,1,w);
			}
			else if(CheckCuts(9, track)){
				fCFContainerP->Fill(value,1,w);
			}
			
		} //end ESD track loop 
		
		fNoEvt->Fill(0); //get the number of analyzed events
	}	

  // Post output data.
	PostData(0, fOutList);

} 

//________________________________________________________________________
void AliAnalysisTaskCentral::Terminate(Option_t *) {

// Called once at the end of the query
	printf("\n\n****************************************\n");
    printf("\tAliAnalysisCentralExtrapolate Terminate... \n");

	TList *outList = dynamic_cast<TList*>(GetOutputData(0));
	if(!outList){
		printf("Unable to get output list! \n");
		return;
	}
 
	AliAnalysisCentralExtrapolate *extPi = new AliAnalysisCentralExtrapolate("extrapolationpi");
	extPi->SetInputList(outList);  //send the outlist to the extrapolation class
	extPi->SetParticle("kPiPlus"); //set the particle type
	extPi->ApplyEff();             //correct the pt distribution !!HAS TO RUN BEFORE extrapolation!!
	extPi->BoltzmannFit();         //fit and extrapolate using Boltzmann-Gibbs Blast wave model
	extPi->TsallisFit();           //fit and extrapolate using Tsallis Blast wave model
	TList *extOutListPi = extPi->GetOutputList();
	
	AliAnalysisCentralExtrapolate *extK = new AliAnalysisCentralExtrapolate("extrapolationK");
	extK->SetInputList(outList);
	extK->SetParticle("kKPlus");
	extK->ApplyEff();
	extK->BoltzmannFit();
	extK->TsallisFit();
	TList *extOutListK = extK->GetOutputList();
	
	AliAnalysisCentralExtrapolate *extP = new AliAnalysisCentralExtrapolate("extrapolationP");
	extP->SetInputList(outList);
	extP->SetParticle("kProton");
	extP->ApplyEff();
	extP->BoltzmannFit();
	extP->TsallisFit();
	TList *extOutListP = extP->GetOutputList();

	fOutList->Add(extOutListPi);
	fOutList->Add(extOutListK);
	fOutList->Add(extOutListP);

}