[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	

}


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

	if(fESD) delete fESD;
	if(fMC) delete fMC;

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


	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

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