[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 "TH1F.h"
#include "TList.h"
#include "TObjArray.h"
#include "TString.h"
#include "TFile.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;
    }
	
	printf("AliAnalysisTaskCentral::Constructor\n");
}


//________________________________________________________________________
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;

	printf("AliAnalysisTaskCentral::Destructor\n");
}

//________________________________________________________________________
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);

	fNoEvt = new TH1F("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();


    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"
	25	#include "TH1F.h"
	26	#include "TList.h"
	27	#include "TObjArray.h"
	28	#include "TString.h"
410ff8cc	29	#include "TFile.h"
bde49c8a	30
	31	#include "AliAnalysisManager.h"
	32	#include "AliMCEvent.h"
	33	#include "AliMCEventHandler.h"
	34	#include "AliESDEvent.h"
	35	#include "AliESDInputHandler.h"
	36	#include "AliESDtrackCuts.h"
	37	#include "AliCFContainer.h"
410ff8cc	38	#include "AliAnalysisDataContainer.h"
410ff8cc	39	#include "AliAnalysisDataSlot.h"
bde49c8a	40
	41	#include "AliAnalysisCentralCutMC.h"
	42	#include "AliAnalysisCentralCutESD.h"
	43	#include "AliAnalysisCentralCutEvtMC.h"
	44	#include "AliAnalysisCentralCutEvtESD.h"
	45	#include "AliAnalysisCentralExtrapolate.h"
	46	#include "AliAnalysisTaskCentral.h"
	47
	48
	49	ClassImp(AliAnalysisTaskCentral)
	50
	51
	52	//________________________________________________________________________
	53	AliAnalysisTaskCentral::AliAnalysisTaskCentral(const char *name)
	54	: AliAnalysisTask(name, "")
	55	,fESD(0), fMC(0)
	56	,fNoEvt(0)
	57	,fCFContainerPi(0)
	58	,fCFContainerK(0)
	59	,fCFContainerP(0)
	60	,fSim(kFALSE)
	61	,fOutList(NULL)
	62
	63	{
	64	// Constructor
	65	printf("AliAnalysisTaskCentral::AliAnalysisTaskCentral(const char *name)\n");
	66
	67	// Define input and output slots here
	68	DefineInput(0, TChain::Class());
	69
	70	DefineOutput(0, TList::Class());
	71
	72
	73	for(Int_t i=0; i<10; i++){
	74	fCutsList[i] = 0;
	75	}
	76
410ff8cc	77	InitCuts(); //initialize the analysis specific cuts
	78	if (!fCutsList) {
	79	Printf("ERROR: fCutsList not available");
	80	return;
	81	}
	82
	83	printf("AliAnalysisTaskCentral::Constructor\n");
bde49c8a	84	}
	85
	86
	87	//________________________________________________________________________
	88	AliAnalysisTaskCentral::~AliAnalysisTaskCentral()
	89	{
	90	// Destructor
	91	// Delete the created objects
	92
	93	if(fESD) delete fESD;
	94	if(fMC) delete fMC;
	95
	96	if(fCutsList){
	97	delete [] fCutsList;
	98	}
	99
	100	if(fNoEvt) delete fNoEvt;
	101
	102	if(fCFContainerPi) delete fCFContainerPi;
	103	if(fCFContainerPi) delete fCFContainerK;
	104	if(fCFContainerPi) delete fCFContainerP;
	105
	106	if(fOutList) delete fOutList;
	107
410ff8cc	108	printf("AliAnalysisTaskCentral::Destructor\n");
bde49c8a	109	}
	110
	111	//________________________________________________________________________
	112	void AliAnalysisTaskCentral::ConnectInputData(Option_t *) {
	113	// get the event from the input chain
	114
	115	TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
	116	printf("tree(%p)\n", (void*)tree);
	117
	118	if (!tree) {
	119	Printf("ERROR: Could not read chain from input slot 0");
	120	} else {
	121	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	122	if (!esdH) {
	123	Printf("ERROR: Could not get ESDInputHandler");
	124	}
	125	else fESD = esdH->GetEvent();
	126
	127	AliMCEventHandler mcH = (AliMCEventHandler) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
	128	if (!mcH) {
	129	Printf("ERROR: Could not get MCInputHandler");
	130	}
	131	else fMC = mcH->MCEvent();
	132	}
	133	}
	134
	135	//________________________________________________________________________
	136	void AliAnalysisTaskCentral::CreateOutputObjects(){
	137	//Creates the output list
	138
410ff8cc	139	// fNoEvt = new TH1F(Form("%s_NoEvt",GetName()), "Number of processed events", 1, 0.0, 1.0);
bde49c8a	140
410ff8cc	141	fNoEvt = new TH1F("TaskCentral_NoEvt", "Number of processed events", 1, 0.0, 1.0);
bde49c8a	142	//set here the:
	143	//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
410ff8cc	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);
bde49c8a	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	if (!fCutsList) {
8a993eac	321	printf("No particle cut list found!\n\n");
8a993eac	322	return;
bde49c8a	323	}
	324	else {
	325	for(Int_t isel=0;isel< 10; isel++){
8a993eac	326	if(!fCutsList[isel]) continue;
bde49c8a	327	TObjArrayIter iter(fCutsList[isel]);
	328	AliAnalysisCuts *cut = 0;
	329
	330	while ((cut = (AliAnalysisCuts*)iter.Next())) {
	331
	332	TString cutName=cut->GetName();
	333	if(!cutName){
	334	printf("No cutname!\n");
8a993eac	335	return;
bde49c8a	336	}
	337
	338	Bool_t checkCut=cutName.Contains("AliAnalysisCentralCutMC");
	339
	340	if(checkCut){
	341	AliAnalysisCentralCutMC newcut = dynamic_cast<AliAnalysisCentralCutMC >(cut);
	342	newcut->ReceiveEvt(obj);
	343	}
	344	}
	345	}
	346	}
	347
	348	}
	349
	350
	351	//________________________________________________________________________
	352	Bool_t AliAnalysisTaskCentral::CheckCuts(Int_t no, TObject *obj) const{
	353
	354	// For each cut run IsSelected();
410ff8cc	355	// printf("AliAnalysisTaskCentral::CheckCuts IN\n");
bde49c8a	356
	357	if(no > 10){
	358	printf("\nAliAnalysisTaskCentral::CheckCuts -> Cut number is not ok! \n");
	359	return kFALSE;
	360	}
	361
	362	if(!fCutsList[no]){
	363	printf("AliAnalysisTaskCentral::CheckCuts -> cuts list problem! \n");
	364	return kFALSE;
	365	}
	366
	367	TObjArrayIter iter(fCutsList[no]);
	368	AliAnalysisCuts *cut = 0;
	369
	370	while((cut = (AliAnalysisCuts*)iter.Next())){
410ff8cc	371	if(!cut->IsSelected(obj)){
	372	// printf("AliAnalysisTaskCentral::CheckCuts OUT\n");
	373	return kFALSE;
	374	}
bde49c8a	375	}
bde49c8a	376
410ff8cc	377	// printf("AliAnalysisTaskCentral::CheckCuts OUT\n");
bde49c8a	378	return kTRUE;
	379	}
	380
	381
	382	//________________________________________________________________________
	383	void AliAnalysisTaskCentral::Exec(Option_t *) {
	384
	385	// Main loop
	386	// Called for each event
	387
bde49c8a	388	Double_t pt, eta;
	389	Double_t w = 1.0;
	390	const Int_t nvar=2; //number of variables on the grid:pt,vtx
	391	Double_t value[nvar];//used to fill the CF container
	392
	393	if(fSim){ // if running on simulations -> look at MC Truth
	394
	395	if (!fMC) {
	396	Printf("ERROR: fMC not available");
	397	return;
	398	}
	399
	400	if(CheckCuts(0, fMC)){ //check event level cuts
	401
	402	SendEvent(fMC);
	403
	404	// MC loop
	405	for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++) {
	406
	407	AliMCParticle particle = dynamic_cast<AliMCParticle>(fMC->GetTrack(ipart));
	408
	409	if(!particle){
	410	printf("\nMCParticle pointer is null!!!\n");
	411	continue;
	412	}
	413
	414	if (!CheckCuts(2, particle)) continue; //check the MC general particle cuts
	415
	416	pt = particle->Pt();
	417	eta = particle->Eta();
	418
	419	if(pt>0) w = 1.0/pt; //invariant distribution
	420
	421	value[0]=pt;
	422	value[1]=eta;
	423
	424	if(CheckCuts(3, particle)){ //fill the right container for each particle
	425	fCFContainerPi->Fill(value,0,w);
	426	}
	427	else if(CheckCuts(4, particle)){
	428	fCFContainerK->Fill(value,0,w);
	429	}
	430	else if(CheckCuts(5, particle)){
	431	fCFContainerP->Fill(value,0,w);
	432	}
	433	} //end MC particle loop
	434	}
	435	}
	436	else{ // if we DONT run in simulated data we fill the MC step of the CFCont with 0
	437	value[0]=0;
	438	value[1]=0;
	439
	440	fCFContainerPi->Fill(value,0);
	441	fCFContainerK->Fill(value,0);
	442	fCFContainerP->Fill(value,0);
	443	}
	444
	445
	446	if (!fESD) {
	447	Printf("ERROR: fESD not available");
	448	return;
	449	}
	450
	451	if(CheckCuts(1, fESD)){
452
453	Printf("There are %d ESD tracks in this event\n", fESD->GetNumberOfTracks());
454
455	// ESD loop
456	for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
457
458	AliESDtrack* track = fESD->GetTrack(iTracks);
459
460	if (!track) {
461	Printf("ERROR: Could not receive track %d", iTracks);
462	continue;
463	}
464
465	if(!CheckCuts(6, track)) continue; //check general ESD track cuts
466
467	pt = track->Pt();
468	eta = track->Eta();
469
470	if(pt>0) w = 1.0/pt; // invariant
471
472	value[0]=pt;
473	value[1]=eta;
474
475	if(CheckCuts(7, track)){
476	fCFContainerPi->Fill(value,1,w);
477	}
478	else if(CheckCuts(8, track)){
479	fCFContainerK->Fill(value,1,w);
480	}
481	else if(CheckCuts(9, track)){
482	fCFContainerP->Fill(value,1,w);
483	}
484
485	} //end ESD track loop
486
487	fNoEvt->Fill(0); //get the number of analyzed events
488	}
489
490	// Post output data.
491	PostData(0, fOutList);
492
493	}
494
495	//________________________________________________________________________
496	void AliAnalysisTaskCentral::Terminate(Option_t *) {
497
498	// Called once at the end of the query
499	printf("\n\n****************************************\n");
500	printf("\tAliAnalysisCentralExtrapolate Terminate... \n");
501
502	TList outList = dynamic_cast<TList>(GetOutputData(0));
503	if(!outList){
504	printf("Unable to get output list! \n");
8a993eac	505	return;
bde49c8a	506	}
	507
	508	AliAnalysisCentralExtrapolate *extPi = new AliAnalysisCentralExtrapolate("extrapolationpi");
	509	extPi->SetInputList(outList); //send the outlist to the extrapolation class
	510	extPi->SetParticle("kPiPlus"); //set the particle type
	511	extPi->ApplyEff(); //correct the pt distribution !!HAS TO RUN BEFORE extrapolation!!
	512	extPi->BoltzmannFit(); //fit and extrapolate using Boltzmann-Gibbs Blast wave model
	513	extPi->TsallisFit(); //fit and extrapolate using Tsallis Blast wave model
	514	TList *extOutListPi = extPi->GetOutputList();
	515
	516	AliAnalysisCentralExtrapolate *extK = new AliAnalysisCentralExtrapolate("extrapolationK");
	517	extK->SetInputList(outList);
	518	extK->SetParticle("kKPlus");
	519	extK->ApplyEff();
	520	extK->BoltzmannFit();
410ff8cc	521	extK->TsallisFit();
bde49c8a	522	TList *extOutListK = extK->GetOutputList();
	523
	524	AliAnalysisCentralExtrapolate *extP = new AliAnalysisCentralExtrapolate("extrapolationP");
	525	extP->SetInputList(outList);
	526	extP->SetParticle("kProton");
	527	extP->ApplyEff();
	528	extP->BoltzmannFit();
410ff8cc	529	extP->TsallisFit();
bde49c8a	530	TList *extOutListP = extP->GetOutputList();
bde49c8a	531
410ff8cc	532
	533	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
	534	if(!mgr){
	535	printf("Unable to get AnalysisManager! \n");
	536	return;
	537	}
	538
	539	TString anaType;
	540	mgr->GetAnalysisTypeString(anaType);
	541
	542	if(anaType.Contains("local")){
	543	fOutList->Add(extOutListPi);
	544	fOutList->Add(extOutListK);
	545	fOutList->Add(extOutListP);
	546	}
	547	else{
	548	AliAnalysisDataContainer *cont = GetOutputSlot(0)->GetContainer();
	549	if(!cont){
	550	printf("Unable to get DataContainer! \n");
	551	return;
	552	}
	553
	554	printf("file name = %s\n", cont->GetFileName());
	555	TFile file(cont->GetFileName(),"update");
	556
	557	file.cd("PWG2Central");
	558
	559	gFile->WriteObject(extOutListPi,"pion_list","SingleKey");
	560	gFile->WriteObject(extOutListK,"kaon_list","SingleKey");
	561	gFile->WriteObject(extOutListP,"proton_list","SingleKey");
	562	file.Close();
	563	}
	564
bde49c8a	565
bde49c8a	566	}