[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 "TCanvas.h"

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

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


ClassImp(AliAnalysisTaskCentral)


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

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

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

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


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

}


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

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

}

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

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

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

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

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

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

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

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

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

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


//________________________________________________________________
void AliAnalysisTaskCentral::InitCuts(){

//Create and set cuts

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

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

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

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


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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

}


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

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

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

}


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

// For each cut run IsSelected();

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

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

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

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

    return kTRUE;
}


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

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

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

    if(CheckCuts(1, fESD)){

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

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

} 

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

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

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

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

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