1 /*************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // ------------------------------------------
17 // analysis task for azimuthal isotropic
18 // expansion in highly central collisions
19 // author: Cristian Andrei
20 // acristian@niham.nipne.ro
21 // ------------------------------------------
27 #include "TObjArray.h"
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"
43 #include "AliAnalysisCentralCutMC.h"
44 #include "AliAnalysisCentralCutESD.h"
45 #include "AliAnalysisCentralCutEvtMC.h"
46 #include "AliAnalysisCentralCutEvtESD.h"
47 #include "AliAnalysisCentralExtrapolate.h"
48 #include "AliAnalysisTaskCentral.h"
51 ClassImp(AliAnalysisTaskCentral)
54 //________________________________________________________________________
55 AliAnalysisTaskCentral::AliAnalysisTaskCentral(const char *name)
56 : AliAnalysisTask(name, "")
67 printf("AliAnalysisTaskCentral::AliAnalysisTaskCentral(const char *name)\n");
69 // Define input and output slots here
70 DefineInput(0, TChain::Class());
72 DefineOutput(0, TList::Class());
75 for(Int_t i=0; i<10; i++){
79 InitCuts(); //initialize the analysis specific cuts
81 Printf("ERROR: fCutsList not available");
88 //________________________________________________________________________
89 AliAnalysisTaskCentral::~AliAnalysisTaskCentral()
92 // Delete the created objects
101 if(fNoEvt) delete fNoEvt;
103 if(fCFContainerPi) delete fCFContainerPi;
104 if(fCFContainerPi) delete fCFContainerK;
105 if(fCFContainerPi) delete fCFContainerP;
107 if(fOutList) delete fOutList;
112 //________________________________________________________________________
113 void AliAnalysisTaskCentral::ConnectInputData(Option_t *) {
114 // get the event from the input chain
116 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
117 printf("tree(%p)\n", (void*)tree);
120 Printf("ERROR: Could not read chain from input slot 0");
122 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
124 Printf("ERROR: Could not get ESDInputHandler");
126 else fESD = esdH->GetEvent();
128 AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
130 Printf("ERROR: Could not get MCInputHandler");
134 fMC = mcH->MCEvent();
140 //________________________________________________________________________
141 void AliAnalysisTaskCentral::CreateOutputObjects(){
142 //Creates the output list
144 // fNoEvt = new TH1F(Form("%s_NoEvt",GetName()), "Number of processed events", 1, 0.0, 1.0);
146 fNoEvt = new TH1D("TaskCentral_NoEvt", "Number of processed events", 1, 0.0, 1.0);
147 //set here the min and max variable values (to do: move these to a separate cuts class)
148 const Double_t ptmin = 0.0 ; //bins in pt
149 const Double_t ptmax = 5.0 ; //bins in y
150 const Double_t etamin = -0.9 ; //bins in pt
151 const Double_t etamax = 0.9 ; //bins in y
152 //apropiate number of bins for histos and CF Container
153 const Int_t nbinpt = 25 ; //bins in pt
154 const Int_t nbineta = 25 ; //bins in y
156 //Correction Framework CONTAINER DEFINITION
157 //the sensitive variables, their indices
160 //Setting up the container grid...
161 const UInt_t nstep = 2 ; //number of selection steps MC & ESD
162 const Int_t nvar = 2 ; //number of variables on the grid:pt,eta
164 //arrays for the number of bins in each dimension
165 //nbin is defined above
170 //arrays for lower bounds :
171 Double_t *binLim1=new Double_t[nbinpt+1];
172 Double_t *binLim2=new Double_t[nbineta+1];
174 //values for bin lower bounds
175 for(Int_t i=0; i<=nbinpt; i++) binLim1[i]=(Double_t)ptmin + (ptmax-ptmin)/nbinpt*(Double_t)i ;
176 for(Int_t i=0; i<=nbineta; i++) binLim2[i]=(Double_t)etamin + (etamax-etamin) /nbineta*(Double_t)i ;
179 fCFContainerPi = new AliCFContainer("TaskCentral_CFCont_Pi","container for tracks",nstep,nvar,iBin);
180 fCFContainerK = new AliCFContainer("TaskCentral_CFCont_K","container for tracks",nstep,nvar,iBin);
181 fCFContainerP = new AliCFContainer("TaskCentral_CFCont_P","container for tracks",nstep,nvar,iBin);
183 //setting the bin limits
184 fCFContainerPi -> SetBinLimits(ipt,binLim1);
185 fCFContainerPi -> SetBinLimits(ieta,binLim2);
187 fCFContainerK -> SetBinLimits(ipt,binLim1);
188 fCFContainerK -> SetBinLimits(ieta,binLim2);
190 fCFContainerP -> SetBinLimits(ipt,binLim1);
191 fCFContainerP -> SetBinLimits(ieta,binLim2);
193 //outout list creation
194 fOutList = new TList();
195 fOutList->Add(fNoEvt);
196 fOutList->Add(fCFContainerPi);
197 fOutList->Add(fCFContainerK);
198 fOutList->Add(fCFContainerP);
202 //________________________________________________________________
203 void AliAnalysisTaskCentral::InitCuts(){
205 //Create and set cuts
207 //------------------EVENT LEVEL CUTS--------------------
208 //-----------------------MC-----------------------------
209 AliAnalysisCentralCutEvtMC *evMCCuts = new AliAnalysisCentralCutEvtMC();
210 evMCCuts->SetMultiplicityRange(1,100);
211 evMCCuts->SetDirectivityRange(0.0, 1.0);
212 evMCCuts->SetDirUnitRange(0.0, 1.0);
214 //-----------------------ESD----------------------------
215 AliAnalysisCentralCutEvtESD *evESDCuts = new AliAnalysisCentralCutEvtESD();
216 evESDCuts->SetMultiplicityRange(1,100);
217 evESDCuts->SetDirectivityRange(0.0, 1.0);
218 evESDCuts->SetSPDMultiplicityRange(1,20000);
219 evESDCuts->SetSPDDirectivityRange(0.0, 1.0);
221 TObjArray* mcEventCuts = new TObjArray();
222 mcEventCuts->AddLast(evMCCuts);
224 TObjArray* esdEventCuts = new TObjArray();
225 esdEventCuts->AddLast(evESDCuts);
228 //------------------PARTICLE LEVEL CUTS-----------------
229 //-------------------General MC Cuts--------------------
230 AliAnalysisCentralCutMC *mcCutsGen = new AliAnalysisCentralCutMC();
231 mcCutsGen->SetOnlyPrimaries(kTRUE);
232 mcCutsGen->SetPtRange(0.2,4.0);
233 mcCutsGen->SetEtaRange(-0.2,0.2);
235 //-------------------Specific MC Cuts-------------------
236 AliAnalysisCentralCutMC *mcCutsPi = new AliAnalysisCentralCutMC();
237 mcCutsPi->SetPDGCode(211); //211 pion; 321 kaon; 2212 proton
239 AliAnalysisCentralCutMC *mcCutsK = new AliAnalysisCentralCutMC();
240 mcCutsK->SetPDGCode(321); //211 pion; 321 kaon; 2212 proton
242 AliAnalysisCentralCutMC *mcCutsP = new AliAnalysisCentralCutMC();
243 mcCutsP->SetPDGCode(2212); //211 pion; 321 kaon; 2212 proton
245 //--------each task has its own mcList of cuts----------
246 TObjArray *mcListGen = new TObjArray(); //general MC cuts
247 mcListGen->AddLast(mcCutsGen);
249 TObjArray *mcListPi = new TObjArray(); //task pt pions
250 mcListPi->AddLast(mcCutsPi);
252 TObjArray *mcListK = new TObjArray(); //task pt kaons
253 mcListK->AddLast(mcCutsK);
255 TObjArray *mcListP = new TObjArray(); //task pt protons
256 mcListP->AddLast(mcCutsP);
259 //-------------------General ESD Cuts-------------------
260 AliESDtrackCuts *esdCutsGen = new AliESDtrackCuts("AliESDtrackCuts", "Loose");
261 esdCutsGen->SetMinNClustersTPC(50);
262 esdCutsGen->SetMaxChi2PerClusterTPC(2.2);
263 esdCutsGen->SetMaxCovDiagonalElements(0.5,0.5,0.5,0.5,0.5);
264 esdCutsGen->SetRequireTPCRefit(kTRUE);
265 esdCutsGen->SetAcceptKinkDaughters(kFALSE);
266 esdCutsGen->SetMaxNsigmaToVertex(2.0);
267 esdCutsGen->SetRequireSigmaToVertex(kTRUE);
268 esdCutsGen->SetPtRange(0.2,4.0);
269 esdCutsGen->SetEtaRange(-0.2,0.2);
271 //-------------------Specific ESD Cuts------------------
272 AliAnalysisCentralCutESD *esdCutsPi = new AliAnalysisCentralCutESD("AliAnalysisCentralCutESD","NIHAM");
273 esdCutsPi->SetPIDtype("Custom");
274 esdCutsPi->SetPriorFunctions(kFALSE);
275 esdCutsPi->SetPartType(kPiPlus);
277 AliAnalysisCentralCutESD *esdCutsK = new AliAnalysisCentralCutESD("AliAnalysisCentralCutESD","NIHAM");
278 esdCutsK->SetPIDtype("Custom");
279 esdCutsK->SetPriorFunctions(kFALSE);
280 esdCutsK->SetPartType(kKPlus);
282 AliAnalysisCentralCutESD *esdCutsP = new AliAnalysisCentralCutESD("AliAnalysisCentralCutESD","NIHAM");
283 esdCutsP->SetPIDtype("Custom");
284 esdCutsP->SetPriorFunctions(kFALSE);
285 esdCutsP->SetPartType(kProton);
287 //--------each task has its own esdList of cuts---------
288 TObjArray* esdListGen = new TObjArray(); //general ESD track cuts
289 esdListGen->AddLast(esdCutsGen);
291 TObjArray* esdListPi = new TObjArray();
292 esdListPi->AddLast(esdCutsPi);
294 TObjArray* esdListK = new TObjArray();
295 esdListK->AddLast(esdCutsK);
297 TObjArray* esdListP = new TObjArray();
298 esdListP->AddLast(esdCutsP);
300 //------set the cuts to the RIGHT! fCutsList slots-------
302 SetCuts(0, mcEventCuts);
303 SetCuts(1, esdEventCuts);
305 // particle level cuts
306 SetCuts(2, mcListGen);
307 SetCuts(3, mcListPi);
311 SetCuts(6, esdListGen);
312 SetCuts(7, esdListPi);
313 SetCuts(8, esdListK);
314 SetCuts(9, esdListP);
319 //_______________________________________________________________________
320 void AliAnalysisTaskCentral::SendEvent(TObject *obj) const{
322 // Some cuts (ie MC IsPrimary) need the MC Event info
325 printf("No particle cut list found!\n\n");
329 for(Int_t isel=0;isel< 10; isel++){
330 if(!fCutsList[isel]) continue;
331 TObjArrayIter iter(fCutsList[isel]);
332 AliAnalysisCuts *cut = 0;
334 while ((cut = (AliAnalysisCuts*)iter.Next())) {
336 TString cutName=cut->GetName();
338 printf("No cutname!\n");
342 Bool_t checkCut=cutName.Contains("AliAnalysisCentralCutMC");
345 AliAnalysisCentralCutMC *newcut = dynamic_cast<AliAnalysisCentralCutMC *>(cut);
346 newcut->ReceiveEvt(obj);
355 //________________________________________________________________________
356 Bool_t AliAnalysisTaskCentral::CheckCuts(Int_t no, TObject *obj) const{
358 // For each cut run IsSelected();
359 // printf("AliAnalysisTaskCentral::CheckCuts IN\n");
362 printf("\nAliAnalysisTaskCentral::CheckCuts -> Cut number is not ok! \n");
367 printf("AliAnalysisTaskCentral::CheckCuts -> cuts list problem! \n");
371 TObjArrayIter iter(fCutsList[no]);
372 AliAnalysisCuts *cut = 0;
374 while((cut = (AliAnalysisCuts*)iter.Next())){
375 if(!cut->IsSelected(obj)){
376 // printf("AliAnalysisTaskCentral::CheckCuts OUT\n");
381 // printf("AliAnalysisTaskCentral::CheckCuts OUT\n");
386 //________________________________________________________________________
387 void AliAnalysisTaskCentral::Exec(Option_t *) {
390 // Called for each event
394 const Int_t nvar=2; //number of variables on the grid:pt,vtx
395 Double_t value[nvar];//used to fill the CF container
397 if(fSim){ // if running on simulations -> look at MC Truth
400 Printf("ERROR: fMC not available");
404 if(CheckCuts(0, fMC)){ //check event level cuts
409 for (Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++) {
411 AliMCParticle *particle = dynamic_cast<AliMCParticle*>(fMC->GetTrack(ipart));
414 printf("\nMCParticle pointer is null!!!\n");
418 if (!CheckCuts(2, particle)) continue; //check the MC general particle cuts
421 eta = particle->Eta();
423 if(pt>0) w = 1.0/pt; //invariant distribution
428 if(CheckCuts(3, particle)){ //fill the right container for each particle
429 fCFContainerPi->Fill(value,0,w);
431 else if(CheckCuts(4, particle)){
432 fCFContainerK->Fill(value,0,w);
434 else if(CheckCuts(5, particle)){
435 fCFContainerP->Fill(value,0,w);
437 } //end MC particle loop
440 else{ // if we DONT run in simulated data we fill the MC step of the CFCont with 0
444 fCFContainerPi->Fill(value,0);
445 fCFContainerK->Fill(value,0);
446 fCFContainerP->Fill(value,0);
451 Printf("ERROR: fESD not available");
455 if(CheckCuts(1, fESD)){
457 Printf("There are %d ESD tracks in this event\n", fESD->GetNumberOfTracks());
460 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
462 AliESDtrack* track = fESD->GetTrack(iTracks);
465 Printf("ERROR: Could not receive track %d", iTracks);
469 if(!CheckCuts(6, track)) continue; //check general ESD track cuts
474 if(pt>0) w = 1.0/pt; // invariant
479 if(CheckCuts(7, track)){
480 fCFContainerPi->Fill(value,1,w);
482 else if(CheckCuts(8, track)){
483 fCFContainerK->Fill(value,1,w);
485 else if(CheckCuts(9, track)){
486 fCFContainerP->Fill(value,1,w);
489 } //end ESD track loop
491 fNoEvt->Fill(0); //get the number of analyzed events
495 PostData(0, fOutList);
499 //________________________________________________________________________
500 void AliAnalysisTaskCentral::Terminate(Option_t *) {
502 // Called once at the end of the query
503 printf("\n\n****************************************\n");
504 printf("\tAliAnalysisCentralExtrapolate Terminate... \n");
506 TList *outList = dynamic_cast<TList*>(GetOutputData(0));
508 printf("Unable to get output list! \n");
512 AliAnalysisCentralExtrapolate *extPi = new AliAnalysisCentralExtrapolate("extrapolationpi");
513 extPi->SetInputList(outList); //send the outlist to the extrapolation class
514 extPi->SetParticle("kPiPlus"); //set the particle type
515 extPi->ApplyEff(); //correct the pt distribution !!HAS TO RUN BEFORE extrapolation!!
516 extPi->BoltzmannFit(); //fit and extrapolate using Boltzmann-Gibbs Blast wave model
517 // extPi->TsallisFit(); //fit and extrapolate using Tsallis Blast wave model
518 TList *extOutListPi = extPi->GetOutputList();
520 AliAnalysisCentralExtrapolate *extK = new AliAnalysisCentralExtrapolate("extrapolationK");
521 extK->SetInputList(outList);
522 extK->SetParticle("kKPlus");
524 extK->BoltzmannFit();
525 // extK->TsallisFit();
526 TList *extOutListK = extK->GetOutputList();
528 AliAnalysisCentralExtrapolate *extP = new AliAnalysisCentralExtrapolate("extrapolationP");
529 extP->SetInputList(outList);
530 extP->SetParticle("kProton");
532 extP->BoltzmannFit();
533 // extP->TsallisFit();
534 TList *extOutListP = extP->GetOutputList();
537 //----------- Plot the extrapolated spectra -----------------
538 TCanvas *ccorrdata = new TCanvas();
539 ccorrdata->Divide(3,2);
542 ccorrdata->cd(1)->SetLogy();
543 TH1D *extBoltzPi = dynamic_cast<TH1D*>(extOutListPi->FindObject("PtExtBoltzmann"));
545 extBoltzPi->Draw("p e1");
549 ccorrdata->cd(4)->SetLogy();
550 TH1D *extTsallisPi = dynamic_cast<TH1D*>(extOutListPi->FindObject("PtExtTsallis"));
552 extTsallisPi->Draw("p e1");
557 ccorrdata->cd(2)->SetLogy();
558 TH1D *extBoltzK = dynamic_cast<TH1D*>(extOutListK->FindObject("PtExtBoltzmann"));
560 extBoltzK->Draw("p e1");
564 ccorrdata->cd(5)->SetLogy();
565 TH1D *extTsallisK = dynamic_cast<TH1D*>(extOutListK->FindObject("PtExtTsallis"));
567 extTsallisK->Draw("p e1");
572 ccorrdata->cd(3)->SetLogy();
573 TH1D *extBoltzP = dynamic_cast<TH1D*>(extOutListP->FindObject("PtExtBoltzmann"));
575 extBoltzP->Draw("p e1");
579 ccorrdata->cd(6)->SetLogy();
580 TH1D *extTsallisP = dynamic_cast<TH1D*>(extOutListP->FindObject("PtExtTsallis"));
582 extTsallisP->Draw("p e1");
585 // ------------------ Save the results -----------------------
586 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
588 printf("Unable to get AnalysisManager! \n");
593 mgr->GetAnalysisTypeString(anaType);
595 if(anaType.Contains("local")){
596 fOutList->Add(extOutListPi);
597 fOutList->Add(extOutListK);
598 fOutList->Add(extOutListP);
601 AliAnalysisDataContainer *cont = GetOutputSlot(0)->GetContainer();
603 printf("Unable to get DataContainer! \n");
607 printf("file name = %s\n", cont->GetFileName());
608 TFile file(cont->GetFileName(),"update");
610 file.cd("PWG2Central");
612 gFile->WriteObject(extOutListPi,"pion_list","SingleKey");
613 gFile->WriteObject(extOutListK,"kaon_list","SingleKey");
614 gFile->WriteObject(extOutListP,"proton_list","SingleKey");