[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliAnalysisTaskCentral.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.                  *
**************************************************************************/

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

#include "TChain.h"
#include "TTree.h"
#include "TH1D.h"
#include "TList.h"
#include "TObjArray.h"
#include "TString.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TLegend.h"

#include "AliAnalysisManager.h"
#include "AliMCEvent.h"
#include "AliMCEventHandler.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliESDtrackCuts.h"
#include "AliCFContainer.h"
#include "AliAnalysisDataContainer.h"
#include "AliAnalysisDataSlot.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;
    }

	InitCuts(); //initialize the analysis specific cuts	
    if (!fCutsList) {
		Printf("ERROR: fCutsList not available");
		return;
    }

}


//________________________________________________________________________
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");
	    fSim = kFALSE;
	}
	else{
	    fMC = mcH->MCEvent();
	    fSim = kTRUE;
	}
    }
}

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

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

	fNoEvt = new TH1D("TaskCentral_NoEvt", "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("TaskCentral_CFCont_Pi","container for tracks",nstep,nvar,iBin);
	fCFContainerK = new AliCFContainer("TaskCentral_CFCont_K","container for tracks",nstep,nvar,iBin);
	fCFContainerP = new AliCFContainer("TaskCentral_CFCont_P","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();
// 	printf("AliAnalysisTaskCentral::CheckCuts IN\n");

    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)){
// 			printf("AliAnalysisTaskCentral::CheckCuts OUT\n");
			return kFALSE;
		}
    }

// 	printf("AliAnalysisTaskCentral::CheckCuts OUT\n");
    return kTRUE;
}


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

// Main loop
// Called for each event

	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();


//----------- Plot the extrapolated spectra -----------------
	TCanvas *ccorrdata = new TCanvas();
	ccorrdata->Divide(3,2);
	
	ccorrdata->cd(1);
	ccorrdata->cd(1)->SetLogy();
	TH1D *extBoltzPi = dynamic_cast<TH1D*>(extOutListPi->FindObject("PtExtBoltzmann"));
	if(extBoltzPi){
		extBoltzPi->Draw("p e1");
	}

	ccorrdata->cd(4);
	ccorrdata->cd(4)->SetLogy();
	TH1D *extTsallisPi = dynamic_cast<TH1D*>(extOutListPi->FindObject("PtExtTsallis"));
	if(extTsallisPi){
		extTsallisPi->Draw("p e1");
	}


	ccorrdata->cd(2);
	ccorrdata->cd(2)->SetLogy();
	TH1D *extBoltzK = dynamic_cast<TH1D*>(extOutListK->FindObject("PtExtBoltzmann"));
	if(extBoltzK){
		extBoltzK->Draw("p e1");
	}

	ccorrdata->cd(5);
	ccorrdata->cd(5)->SetLogy();
	TH1D *extTsallisK = dynamic_cast<TH1D*>(extOutListK->FindObject("PtExtTsallis"));
	if(extTsallisK){
		extTsallisK->Draw("p e1");
	}


	ccorrdata->cd(3);
	ccorrdata->cd(3)->SetLogy();
	TH1D *extBoltzP = dynamic_cast<TH1D*>(extOutListP->FindObject("PtExtBoltzmann"));
	if(extBoltzP){
		extBoltzP->Draw("p e1");
	}

	ccorrdata->cd(6);
	ccorrdata->cd(6)->SetLogy();
	TH1D *extTsallisP = dynamic_cast<TH1D*>(extOutListP->FindObject("PtExtTsallis"));
	if(extTsallisP){
		extTsallisP->Draw("p e1");
	}

// ------------------ Save the results -----------------------
    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
	if(!mgr){
		printf("Unable to get AnalysisManager! \n");
        return;
	}

	TString anaType;
	mgr->GetAnalysisTypeString(anaType);

	if(anaType.Contains("local")){
		fOutList->Add(extOutListPi);
		fOutList->Add(extOutListK);
		fOutList->Add(extOutListP);
	}
	else{
		AliAnalysisDataContainer *cont = GetOutputSlot(0)->GetContainer();
		if(!cont){
			printf("Unable to get DataContainer! \n");
			return;
		}

		printf("file name = %s\n", cont->GetFileName());
		TFile file(cont->GetFileName(),"update");

		file.cd("PWG2Central");

		gFile->WriteObject(extOutListPi,"pion_list","SingleKey");
		gFile->WriteObject(extOutListK,"kaon_list","SingleKey");
		gFile->WriteObject(extOutListP,"proton_list","SingleKey");
		file.Close();
	}

}
Commit	Line	Data
bde49c8a	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
410ff8cc	16	// ------------------------------------------
	17	// analysis task for azimuthal isotropic
	18	// expansion in highly central collisions
	19	// author: Cristian Andrei
	20	// acristian@niham.nipne.ro
	21	// ------------------------------------------
bde49c8a	22
	23	#include "TChain.h"
	24	#include "TTree.h"
0fb1f0cf	25	#include "TH1D.h"
bde49c8a	26	#include "TList.h"
	27	#include "TObjArray.h"
	28	#include "TString.h"
410ff8cc	29	#include "TFile.h"
0fb1f0cf	30	#include "TCanvas.h"
0fb1f0cf	31	#include "TLegend.h"
bde49c8a	32
	33	#include "AliAnalysisManager.h"
	34	#include "AliMCEvent.h"
	35	#include "AliMCEventHandler.h"
	36	#include "AliESDEvent.h"
	37	#include "AliESDInputHandler.h"
	38	#include "AliESDtrackCuts.h"
	39	#include "AliCFContainer.h"
410ff8cc	40	#include "AliAnalysisDataContainer.h"
410ff8cc	41	#include "AliAnalysisDataSlot.h"
bde49c8a	42
	43	#include "AliAnalysisCentralCutMC.h"
	44	#include "AliAnalysisCentralCutESD.h"
	45	#include "AliAnalysisCentralCutEvtMC.h"
	46	#include "AliAnalysisCentralCutEvtESD.h"
	47	#include "AliAnalysisCentralExtrapolate.h"
	48	#include "AliAnalysisTaskCentral.h"
	49
	50
	51	ClassImp(AliAnalysisTaskCentral)
	52
	53
	54	//________________________________________________________________________
	55	AliAnalysisTaskCentral::AliAnalysisTaskCentral(const char *name)
	56	: AliAnalysisTask(name, "")
	57	,fESD(0), fMC(0)
	58	,fNoEvt(0)
	59	,fCFContainerPi(0)
	60	,fCFContainerK(0)
	61	,fCFContainerP(0)
	62	,fSim(kFALSE)
	63	,fOutList(NULL)
	64
	65	{
	66	// Constructor
	67	printf("AliAnalysisTaskCentral::AliAnalysisTaskCentral(const char *name)\n");
	68
	69	// Define input and output slots here
	70	DefineInput(0, TChain::Class());
	71
	72	DefineOutput(0, TList::Class());
	73
	74
	75	for(Int_t i=0; i<10; i++){
	76	fCutsList[i] = 0;
0fb1f0cf	77	}
bde49c8a	78
410ff8cc	79	InitCuts(); //initialize the analysis specific cuts
	80	if (!fCutsList) {
	81	Printf("ERROR: fCutsList not available");
	82	return;
	83	}
0fb1f0cf	84
bde49c8a	85	}
	86
	87
	88	//________________________________________________________________________
	89	AliAnalysisTaskCentral::~AliAnalysisTaskCentral()
	90	{
	91	// Destructor
	92	// Delete the created objects
	93
	94	if(fESD) delete fESD;
	95	if(fMC) delete fMC;
0fb1f0cf	96
bde49c8a	97	if(fCutsList){
	98	delete [] fCutsList;
	99	}
0fb1f0cf	100
bde49c8a	101	if(fNoEvt) delete fNoEvt;
0fb1f0cf	102
bde49c8a	103	if(fCFContainerPi) delete fCFContainerPi;
	104	if(fCFContainerPi) delete fCFContainerK;
	105	if(fCFContainerPi) delete fCFContainerP;
0fb1f0cf	106
bde49c8a	107	if(fOutList) delete fOutList;
bde49c8a	108
0fb1f0cf	109
bde49c8a	110	}
	111
	112	//________________________________________________________________________
	113	void AliAnalysisTaskCentral::ConnectInputData(Option_t *) {
	114	// get the event from the input chain
	115
	116	TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
	117	printf("tree(%p)\n", (void*)tree);
	118
	119	if (!tree) {
	120	Printf("ERROR: Could not read chain from input slot 0");
	121	} else {
	122	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	123	if (!esdH) {
	124	Printf("ERROR: Could not get ESDInputHandler");
0fb1f0cf	125	}
bde49c8a	126	else fESD = esdH->GetEvent();
0fb1f0cf	127
bde49c8a	128	AliMCEventHandler mcH = (AliMCEventHandler) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
	129	if (!mcH) {
	130	Printf("ERROR: Could not get MCInputHandler");
7d3a1a1c	131	fSim = kFALSE;
	132	}
	133	else{
	134	fMC = mcH->MCEvent();
	135	fSim = kTRUE;
bde49c8a	136	}
0fb1f0cf	137	}
bde49c8a	138	}
	139
	140	//________________________________________________________________________
	141	void AliAnalysisTaskCentral::CreateOutputObjects(){
	142	//Creates the output list
	143
410ff8cc	144	// fNoEvt = new TH1F(Form("%s_NoEvt",GetName()), "Number of processed events", 1, 0.0, 1.0);
bde49c8a	145
0fb1f0cf	146	fNoEvt = new TH1D("TaskCentral_NoEvt", "Number of processed events", 1, 0.0, 1.0);
0fb1f0cf	147	//set here the min and max variable values (to do: move these to a separate cuts class)
bde49c8a	148	const Double_t ptmin = 0.0 ; //bins in pt
0fb1f0cf	149	const Double_t ptmax = 5.0 ; //bins in y
bde49c8a	150	const Double_t etamin = -0.9 ; //bins in pt
0fb1f0cf	151	const Double_t etamax = 0.9 ; //bins in y
bde49c8a	152	//apropiate number of bins for histos and CF Container
bde49c8a	153	const Int_t nbinpt = 25 ; //bins in pt
0fb1f0cf	154	const Int_t nbineta = 25 ; //bins in y
bde49c8a	155
	156	//Correction Framework CONTAINER DEFINITION
	157	//the sensitive variables, their indices
	158	UInt_t ipt = 0;
	159	UInt_t ieta = 1;
0fb1f0cf	160	//Setting up the container grid...
bde49c8a	161	const UInt_t nstep = 2 ; //number of selection steps MC & ESD
bde49c8a	162	const Int_t nvar = 2 ; //number of variables on the grid:pt,eta
0fb1f0cf	163
bde49c8a	164	//arrays for the number of bins in each dimension
0fb1f0cf	165	//nbin is defined above
bde49c8a	166	Int_t iBin[nvar];
	167	iBin[0]=nbinpt;
	168	iBin[1]=nbineta;
	169
	170	//arrays for lower bounds :
	171	Double_t *binLim1=new Double_t[nbinpt+1];
	172	Double_t *binLim2=new Double_t[nbineta+1];
	173
	174	//values for bin lower bounds
0fb1f0cf	175	for(Int_t i=0; i<=nbinpt; i++) binLim1[i]=(Double_t)ptmin + (ptmax-ptmin)/nbinpt*(Double_t)i ;
bde49c8a	176	for(Int_t i=0; i<=nbineta; i++) binLim2[i]=(Double_t)etamin + (etamax-etamin) /nbineta*(Double_t)i ;
	177
	178	//container creation
410ff8cc	179	fCFContainerPi = new AliCFContainer("TaskCentral_CFCont_Pi","container for tracks",nstep,nvar,iBin);
	180	fCFContainerK = new AliCFContainer("TaskCentral_CFCont_K","container for tracks",nstep,nvar,iBin);
	181	fCFContainerP = new AliCFContainer("TaskCentral_CFCont_P","container for tracks",nstep,nvar,iBin);
0fb1f0cf	182
bde49c8a	183	//setting the bin limits
	184	fCFContainerPi -> SetBinLimits(ipt,binLim1);
	185	fCFContainerPi -> SetBinLimits(ieta,binLim2);
0fb1f0cf	186
bde49c8a	187	fCFContainerK -> SetBinLimits(ipt,binLim1);
bde49c8a	188	fCFContainerK -> SetBinLimits(ieta,binLim2);
0fb1f0cf	189
bde49c8a	190	fCFContainerP -> SetBinLimits(ipt,binLim1);
bde49c8a	191	fCFContainerP -> SetBinLimits(ieta,binLim2);
0fb1f0cf	192
bde49c8a	193	//outout list creation
	194	fOutList = new TList();
	195	fOutList->Add(fNoEvt);
	196	fOutList->Add(fCFContainerPi);
	197	fOutList->Add(fCFContainerK);
	198	fOutList->Add(fCFContainerP);
	199	}
	200
	201
	202	//________________________________________________________________
	203	void AliAnalysisTaskCentral::InitCuts(){
	204
	205	//Create and set cuts
	206
	207	//------------------EVENT LEVEL CUTS--------------------
	208	//-----------------------MC-----------------------------
	209	AliAnalysisCentralCutEvtMC *evMCCuts = new AliAnalysisCentralCutEvtMC();
	210	evMCCuts->SetMultiplicityRange(1,100);
	211	evMCCuts->SetDirectivityRange(0.0, 1.0);
	212	evMCCuts->SetDirUnitRange(0.0, 1.0);
	213
	214	//-----------------------ESD----------------------------
	215	AliAnalysisCentralCutEvtESD *evESDCuts = new AliAnalysisCentralCutEvtESD();
	216	evESDCuts->SetMultiplicityRange(1,100);
	217	evESDCuts->SetDirectivityRange(0.0, 1.0);
	218	evESDCuts->SetSPDMultiplicityRange(1,20000);
	219	evESDCuts->SetSPDDirectivityRange(0.0, 1.0);
	220
	221	TObjArray* mcEventCuts = new TObjArray();
	222	mcEventCuts->AddLast(evMCCuts);
	223
	224	TObjArray* esdEventCuts = new TObjArray();
	225	esdEventCuts->AddLast(evESDCuts);
	226
	227
	228	//------------------PARTICLE LEVEL CUTS-----------------
	229	//-------------------General MC Cuts--------------------
	230	AliAnalysisCentralCutMC *mcCutsGen = new AliAnalysisCentralCutMC();
	231	mcCutsGen->SetOnlyPrimaries(kTRUE);
	232	mcCutsGen->SetPtRange(0.2,4.0);
	233	mcCutsGen->SetEtaRange(-0.2,0.2);
	234
	235	//-------------------Specific MC Cuts-------------------
	236	AliAnalysisCentralCutMC *mcCutsPi = new AliAnalysisCentralCutMC();
	237	mcCutsPi->SetPDGCode(211); //211 pion; 321 kaon; 2212 proton
	238
	239	AliAnalysisCentralCutMC *mcCutsK = new AliAnalysisCentralCutMC();
	240	mcCutsK->SetPDGCode(321); //211 pion; 321 kaon; 2212 proton
	241
	242	AliAnalysisCentralCutMC *mcCutsP = new AliAnalysisCentralCutMC();
	243	mcCutsP->SetPDGCode(2212); //211 pion; 321 kaon; 2212 proton
	244
	245	//--------each task has its own mcList of cuts----------
	246	TObjArray *mcListGen = new TObjArray(); //general MC cuts
	247	mcListGen->AddLast(mcCutsGen);
	248
	249	TObjArray *mcListPi = new TObjArray(); //task pt pions
	250	mcListPi->AddLast(mcCutsPi);
	251
	252	TObjArray *mcListK = new TObjArray(); //task pt kaons
	253	mcListK->AddLast(mcCutsK);
	254
	255	TObjArray *mcListP = new TObjArray(); //task pt protons
	256	mcListP->AddLast(mcCutsP);
257
258
259	//-------------------General ESD Cuts-------------------
260	AliESDtrackCuts *esdCutsGen = new AliESDtrackCuts("AliESDtrackCuts", "Loose");
261	esdCutsGen->SetMinNClustersTPC(50);
262	esdCutsGen->SetMaxChi2PerClusterTPC(2.2);
263	esdCutsGen->SetMaxCovDiagonalElements(0.5,0.5,0.5,0.5,0.5);
264	esdCutsGen->SetRequireTPCRefit(kTRUE);
265	esdCutsGen->SetAcceptKinkDaughters(kFALSE);
266	esdCutsGen->SetMaxNsigmaToVertex(2.0);
267	esdCutsGen->SetRequireSigmaToVertex(kTRUE);
268	esdCutsGen->SetPtRange(0.2,4.0);
269	esdCutsGen->SetEtaRange(-0.2,0.2);
270
271	//-------------------Specific ESD Cuts------------------
272	AliAnalysisCentralCutESD *esdCutsPi = new AliAnalysisCentralCutESD("AliAnalysisCentralCutESD","NIHAM");
273	esdCutsPi->SetPIDtype("Custom");
274	esdCutsPi->SetPriorFunctions(kFALSE);
275	esdCutsPi->SetPartType(kPiPlus);
276
277	AliAnalysisCentralCutESD *esdCutsK = new AliAnalysisCentralCutESD("AliAnalysisCentralCutESD","NIHAM");
278	esdCutsK->SetPIDtype("Custom");
279	esdCutsK->SetPriorFunctions(kFALSE);
280	esdCutsK->SetPartType(kKPlus);
281
282	AliAnalysisCentralCutESD *esdCutsP = new AliAnalysisCentralCutESD("AliAnalysisCentralCutESD","NIHAM");
283	esdCutsP->SetPIDtype("Custom");
284	esdCutsP->SetPriorFunctions(kFALSE);
285	esdCutsP->SetPartType(kProton);
286
287	//--------each task has its own esdList of cuts---------
288	TObjArray* esdListGen = new TObjArray(); //general ESD track cuts
289	esdListGen->AddLast(esdCutsGen);
290
291	TObjArray* esdListPi = new TObjArray();
292	esdListPi->AddLast(esdCutsPi);
293
294	TObjArray* esdListK = new TObjArray();
295	esdListK->AddLast(esdCutsK);
296
297	TObjArray* esdListP = new TObjArray();
298	esdListP->AddLast(esdCutsP);
299
300	//------set the cuts to the RIGHT! fCutsList slots-------
0fb1f0cf	301	// event level cuts
bde49c8a	302	SetCuts(0, mcEventCuts);
bde49c8a	303	SetCuts(1, esdEventCuts);
0fb1f0cf	304
0fb1f0cf	305	// particle level cuts
bde49c8a	306	SetCuts(2, mcListGen);
	307	SetCuts(3, mcListPi);
	308	SetCuts(4, mcListK);
	309	SetCuts(5, mcListP);
0fb1f0cf	310
bde49c8a	311	SetCuts(6, esdListGen);
	312	SetCuts(7, esdListPi);
	313	SetCuts(8, esdListK);
	314	SetCuts(9, esdListP);
	315
	316	}
	317
	318
	319	//_______________________________________________________________________
	320	void AliAnalysisTaskCentral::SendEvent(TObject *obj) const{
	321
	322	// Some cuts (ie MC IsPrimary) need the MC Event info
	323
	324	if (!fCutsList) {
8a993eac	325	printf("No particle cut list found!\n\n");
8a993eac	326	return;
bde49c8a	327	}
	328	else {
	329	for(Int_t isel=0;isel< 10; isel++){
8a993eac	330	if(!fCutsList[isel]) continue;
bde49c8a	331	TObjArrayIter iter(fCutsList[isel]);
	332	AliAnalysisCuts *cut = 0;
	333
	334	while ((cut = (AliAnalysisCuts*)iter.Next())) {
	335
	336	TString cutName=cut->GetName();
	337	if(!cutName){
	338	printf("No cutname!\n");
8a993eac	339	return;
bde49c8a	340	}
	341
	342	Bool_t checkCut=cutName.Contains("AliAnalysisCentralCutMC");
	343
	344	if(checkCut){
	345	AliAnalysisCentralCutMC newcut = dynamic_cast<AliAnalysisCentralCutMC >(cut);
	346	newcut->ReceiveEvt(obj);
	347	}
	348	}
	349	}
	350	}
	351
	352	}
	353
	354
	355	//________________________________________________________________________
	356	Bool_t AliAnalysisTaskCentral::CheckCuts(Int_t no, TObject *obj) const{
	357
	358	// For each cut run IsSelected();
410ff8cc	359	// printf("AliAnalysisTaskCentral::CheckCuts IN\n");
bde49c8a	360
	361	if(no > 10){
	362	printf("\nAliAnalysisTaskCentral::CheckCuts -> Cut number is not ok! \n");
	363	return kFALSE;
	364	}
	365
	366	if(!fCutsList[no]){
	367	printf("AliAnalysisTaskCentral::CheckCuts -> cuts list problem! \n");
	368	return kFALSE;
	369	}
	370
	371	TObjArrayIter iter(fCutsList[no]);
	372	AliAnalysisCuts *cut = 0;
	373
	374	while((cut = (AliAnalysisCuts*)iter.Next())){
410ff8cc	375	if(!cut->IsSelected(obj)){
	376	// printf("AliAnalysisTaskCentral::CheckCuts OUT\n");
	377	return kFALSE;
	378	}
bde49c8a	379	}
bde49c8a	380
410ff8cc	381	// printf("AliAnalysisTaskCentral::CheckCuts OUT\n");
bde49c8a	382	return kTRUE;
	383	}
	384
	385
	386	//________________________________________________________________________
	387	void AliAnalysisTaskCentral::Exec(Option_t *) {
	388
	389	// Main loop
	390	// Called for each event
0fb1f0cf	391
bde49c8a	392	Double_t pt, eta;
0fb1f0cf	393	Double_t w = 1.0;
bde49c8a	394	const Int_t nvar=2; //number of variables on the grid:pt,vtx
bde49c8a	395	Double_t value[nvar];//used to fill the CF container
0fb1f0cf	396
bde49c8a	397	if(fSim){ // if running on simulations -> look at MC Truth
	398
	399	if (!fMC) {
	400	Printf("ERROR: fMC not available");
	401	return;
	402	}
0fb1f0cf	403
bde49c8a	404	if(CheckCuts(0, fMC)){ //check event level cuts
	405
	406	SendEvent(fMC);
	407
	408	// MC loop
	409	for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++) {
	410
	411	AliMCParticle particle = dynamic_cast<AliMCParticle>(fMC->GetTrack(ipart));
	412
	413	if(!particle){
	414	printf("\nMCParticle pointer is null!!!\n");
	415	continue;
	416	}
	417
	418	if (!CheckCuts(2, particle)) continue; //check the MC general particle cuts
	419
	420	pt = particle->Pt();
	421	eta = particle->Eta();
	422
	423	if(pt>0) w = 1.0/pt; //invariant distribution
	424
	425	value[0]=pt;
	426	value[1]=eta;
	427
	428	if(CheckCuts(3, particle)){ //fill the right container for each particle
	429	fCFContainerPi->Fill(value,0,w);
	430	}
	431	else if(CheckCuts(4, particle)){
	432	fCFContainerK->Fill(value,0,w);
	433	}
	434	else if(CheckCuts(5, particle)){
	435	fCFContainerP->Fill(value,0,w);
	436	}
	437	} //end MC particle loop
	438	}
	439	}
0fb1f0cf	440	else{ // if we DONT run in simulated data we fill the MC step of the CFCont with 0
bde49c8a	441	value[0]=0;
	442	value[1]=0;
	443
	444	fCFContainerPi->Fill(value,0);
	445	fCFContainerK->Fill(value,0);
	446	fCFContainerP->Fill(value,0);
	447	}
	448
	449
	450	if (!fESD) {
	451	Printf("ERROR: fESD not available");
	452	return;
	453	}
	454
	455	if(CheckCuts(1, fESD)){
	456
	457	Printf("There are %d ESD tracks in this event\n", fESD->GetNumberOfTracks());
	458
	459	// ESD loop
	460	for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
	461
	462	AliESDtrack* track = fESD->GetTrack(iTracks);
	463
	464	if (!track) {
	465	Printf("ERROR: Could not receive track %d", iTracks);
	466	continue;
	467	}
	468
	469	if(!CheckCuts(6, track)) continue; //check general ESD track cuts
	470
	471	pt = track->Pt();
	472	eta = track->Eta();
	473
	474	if(pt>0) w = 1.0/pt; // invariant
	475
	476	value[0]=pt;
	477	value[1]=eta;
	478
	479	if(CheckCuts(7, track)){
	480	fCFContainerPi->Fill(value,1,w);
	481	}
	482	else if(CheckCuts(8, track)){
	483	fCFContainerK->Fill(value,1,w);
	484	}
	485	else if(CheckCuts(9, track)){
	486	fCFContainerP->Fill(value,1,w);
	487	}
	488
	489	} //end ESD track loop
	490
	491	fNoEvt->Fill(0); //get the number of analyzed events
0fb1f0cf	492	}
bde49c8a	493
	494	// Post output data.
	495	PostData(0, fOutList);
	496
0fb1f0cf	497	}
bde49c8a	498
	499	//________________________________________________________________________
	500	void AliAnalysisTaskCentral::Terminate(Option_t *) {
	501
	502	// Called once at the end of the query
	503	printf("\n\n****************************************\n");
	504	printf("\tAliAnalysisCentralExtrapolate Terminate... \n");
	505
	506	TList outList = dynamic_cast<TList>(GetOutputData(0));
	507	if(!outList){
	508	printf("Unable to get output list! \n");
8a993eac	509	return;
bde49c8a	510	}
0fb1f0cf	511
bde49c8a	512	AliAnalysisCentralExtrapolate *extPi = new AliAnalysisCentralExtrapolate("extrapolationpi");
	513	extPi->SetInputList(outList); //send the outlist to the extrapolation class
	514	extPi->SetParticle("kPiPlus"); //set the particle type
	515	extPi->ApplyEff(); //correct the pt distribution !!HAS TO RUN BEFORE extrapolation!!
	516	extPi->BoltzmannFit(); //fit and extrapolate using Boltzmann-Gibbs Blast wave model
32692aec	517	extPi->TsallisFit(); //fit and extrapolate using Tsallis Blast wave model
bde49c8a	518	TList *extOutListPi = extPi->GetOutputList();
0fb1f0cf	519
bde49c8a	520	AliAnalysisCentralExtrapolate *extK = new AliAnalysisCentralExtrapolate("extrapolationK");
	521	extK->SetInputList(outList);
	522	extK->SetParticle("kKPlus");
	523	extK->ApplyEff();
	524	extK->BoltzmannFit();
32692aec	525	extK->TsallisFit();
bde49c8a	526	TList *extOutListK = extK->GetOutputList();
0fb1f0cf	527
bde49c8a	528	AliAnalysisCentralExtrapolate *extP = new AliAnalysisCentralExtrapolate("extrapolationP");
	529	extP->SetInputList(outList);
	530	extP->SetParticle("kProton");
	531	extP->ApplyEff();
	532	extP->BoltzmannFit();
32692aec	533	extP->TsallisFit();
bde49c8a	534	TList *extOutListP = extP->GetOutputList();
bde49c8a	535
410ff8cc	536
0fb1f0cf	537	//----------- Plot the extrapolated spectra -----------------
	538	TCanvas *ccorrdata = new TCanvas();
	539	ccorrdata->Divide(3,2);
	540
	541	ccorrdata->cd(1);
	542	ccorrdata->cd(1)->SetLogy();
	543	TH1D extBoltzPi = dynamic_cast<TH1D>(extOutListPi->FindObject("PtExtBoltzmann"));
	544	if(extBoltzPi){
	545	extBoltzPi->Draw("p e1");
	546	}
	547
	548	ccorrdata->cd(4);
	549	ccorrdata->cd(4)->SetLogy();
	550	TH1D extTsallisPi = dynamic_cast<TH1D>(extOutListPi->FindObject("PtExtTsallis"));
	551	if(extTsallisPi){
	552	extTsallisPi->Draw("p e1");
	553	}
	554
	555
	556	ccorrdata->cd(2);
	557	ccorrdata->cd(2)->SetLogy();
	558	TH1D extBoltzK = dynamic_cast<TH1D>(extOutListK->FindObject("PtExtBoltzmann"));
	559	if(extBoltzK){
	560	extBoltzK->Draw("p e1");
	561	}
	562
	563	ccorrdata->cd(5);
	564	ccorrdata->cd(5)->SetLogy();
	565	TH1D extTsallisK = dynamic_cast<TH1D>(extOutListK->FindObject("PtExtTsallis"));
	566	if(extTsallisK){
	567	extTsallisK->Draw("p e1");
	568	}
	569
	570
	571	ccorrdata->cd(3);
	572	ccorrdata->cd(3)->SetLogy();
	573	TH1D extBoltzP = dynamic_cast<TH1D>(extOutListP->FindObject("PtExtBoltzmann"));
	574	if(extBoltzP){
	575	extBoltzP->Draw("p e1");
	576	}
	577
	578	ccorrdata->cd(6);
	579	ccorrdata->cd(6)->SetLogy();
	580	TH1D extTsallisP = dynamic_cast<TH1D>(extOutListP->FindObject("PtExtTsallis"));
	581	if(extTsallisP){
	582	extTsallisP->Draw("p e1");
	583	}
	584
	585	// ------------------ Save the results -----------------------
410ff8cc	586	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
	587	if(!mgr){
	588	printf("Unable to get AnalysisManager! \n");
	589	return;
	590	}
	591
	592	TString anaType;
	593	mgr->GetAnalysisTypeString(anaType);
	594
	595	if(anaType.Contains("local")){
	596	fOutList->Add(extOutListPi);
	597	fOutList->Add(extOutListK);
	598	fOutList->Add(extOutListP);
	599	}
	600	else{
	601	AliAnalysisDataContainer *cont = GetOutputSlot(0)->GetContainer();
	602	if(!cont){
	603	printf("Unable to get DataContainer! \n");
	604	return;
	605	}
0fb1f0cf	606
410ff8cc	607	printf("file name = %s\n", cont->GetFileName());
410ff8cc	608	TFile file(cont->GetFileName(),"update");
0fb1f0cf	609
410ff8cc	610	file.cd("PWG2Central");
0fb1f0cf	611
410ff8cc	612	gFile->WriteObject(extOutListPi,"pion_list","SingleKey");
	613	gFile->WriteObject(extOutListK,"kaon_list","SingleKey");
	614	gFile->WriteObject(extOutListP,"proton_list","SingleKey");
	615	file.Close();
	616	}
	617
bde49c8a	618	}