coverity fix
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / MultEvShape / AliAnalysisTaskCentral.cxx
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 }