Coverity 17281, 17282
[u/mrichter/AliRoot.git] / ANALYSIS / AliEPSelectionTask.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 //   Class AliEPSelectionTask
18 //   Class to determine event plane            
19 //   author: Alberica Toia, Johanna Gramling
20 //*****************************************************
21
22 #include "AliEPSelectionTask.h"
23
24 #include <TTree.h>
25 #include <TList.h>
26 #include <TH1F.h>
27 #include <TH2F.h>
28 #include <TProfile.h>
29 #include <TFile.h>
30 #include <TObjString.h>
31 #include <TString.h>
32 #include <TCanvas.h>
33 #include <TROOT.h>
34 #include <TDirectory.h>
35 #include <TSystem.h>
36 #include <iostream>
37 #include <TRandom2.h>
38 #include <TArrayF.h>
39
40 #include "AliAnalysisManager.h"
41 #include "AliVEvent.h"
42 #include "AliESD.h"
43 #include "AliESDEvent.h"
44 #include "AliMCEvent.h"
45 #include "AliESDtrack.h"
46 #include "AliESDtrackCuts.h"
47 #include "AliESDHeader.h"
48 #include "AliESDInputHandler.h"
49 #include "AliAODHandler.h"
50 #include "AliAODEvent.h"
51 #include "AliHeader.h"
52 #include "AliGenHijingEventHeader.h"
53 #include "AliAnalysisTaskSE.h"
54 #include "AliPhysicsSelectionTask.h"
55 #include "AliPhysicsSelection.h"
56 #include "AliBackgroundSelection.h"
57 #include "AliESDUtils.h"
58 #include "AliOADBContainer.h"
59
60 #include "AliEventplane.h"
61
62 ClassImp(AliEPSelectionTask)
63
64 //________________________________________________________________________
65 AliEPSelectionTask::AliEPSelectionTask():
66 AliAnalysisTaskSE(),
67   fDebug(0),
68   fAnalysisInput("ESD"),
69   fTrackType("TPC"),
70   fUseMCRP(kFALSE),
71   fUsePhiWeight(kFALSE),
72   fUsePtWeight(kFALSE),
73   fSaveTrackContribution(kFALSE),
74   fuserphidist(kFALSE),
75   fusercuts(kFALSE),
76   frunNumber(-15),
77   fESDtrackCuts(0),
78   ftracklist(0),
79   fEPContainer(0),
80   fPhiDist(0),
81   fQVector(0),
82   fQContributionX(0),
83   fQContributionY(0),
84   fEventplaneQ(0),
85   fQsub1(0),
86   fQsub2(0),
87   fQsubRes(0),  
88   fOutputList(0),
89   fHOutEventplaneQ(0),
90   fHOutPhi(0),
91   fHOutPhiCorr(0),
92   fHOutsub1sub2(0),
93   fHOutNTEPRes(0),
94   fHOutPTPsi(0),
95   fHOutDiff(0),
96   fHOutleadPTPsi(0)
97 {   
98   // Default constructor
99   AliInfo("Event Plane Selection enabled.");
100 }   
101
102 //________________________________________________________________________
103 AliEPSelectionTask::AliEPSelectionTask(const char *name):
104   AliAnalysisTaskSE(name),
105   fDebug(0),
106   fAnalysisInput("ESD"),
107   fTrackType("TPC"),
108   fUseMCRP(kFALSE),
109   fUsePhiWeight(kFALSE),
110   fUsePtWeight(kFALSE),  
111   fSaveTrackContribution(kFALSE),
112   fuserphidist(kFALSE),
113   fusercuts(kFALSE),
114   frunNumber(-15),
115   fESDtrackCuts(0),
116   ftracklist(0),
117   fEPContainer(0),
118   fPhiDist(0),
119   fQVector(0),
120   fQContributionX(0),
121   fQContributionY(0),
122   fEventplaneQ(0),
123   fQsub1(0),
124   fQsub2(0),
125   fQsubRes(0),
126   fOutputList(0),
127   fHOutEventplaneQ(0),
128   fHOutPhi(0),
129   fHOutPhiCorr(0),
130   fHOutsub1sub2(0),
131   fHOutNTEPRes(0),
132   fHOutPTPsi(0),
133   fHOutDiff(0),
134   fHOutleadPTPsi(0)
135 {
136   // Default constructor
137   AliInfo("Event Plane Selection enabled.");
138   DefineOutput(1, TList::Class());
139 }
140  
141 //________________________________________________________________________
142 AliEPSelectionTask::~AliEPSelectionTask()
143 {
144   // Destructor  
145   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
146       delete fOutputList;
147       fOutputList = 0;
148   }
149   if (fESDtrackCuts){
150       delete fESDtrackCuts;
151       fESDtrackCuts = 0;
152   }
153   if (fQVector){
154       delete fQVector;
155       fQVector = 0;
156   }
157   if (fQsub1){
158       delete fQsub1;
159       fQsub1 = 0;
160   }
161   if (fQsub2){
162       delete fQsub2;
163       fQsub2 = 0;
164   }
165   if (fPhiDist){
166       delete fPhiDist;
167       fPhiDist = 0;
168   }
169   if (ftracklist){
170       delete ftracklist;
171       ftracklist = 0;
172   }
173   if (fEPContainer){
174       delete fEPContainer;
175       fEPContainer = 0;
176   }
177 }  
178
179 //________________________________________________________________________
180 void AliEPSelectionTask::UserCreateOutputObjects()
181 {  
182   // Create the output containers
183   if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
184   AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
185
186   fOutputList = new TList();
187   fOutputList->SetOwner();
188   fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
189   fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
190   fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
191   fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
192   fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
193   fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
194   fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
195   fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
196
197   fOutputList->Add(fHOutEventplaneQ);
198   fOutputList->Add(fHOutPhi);
199   fOutputList->Add(fHOutPhiCorr);
200   fOutputList->Add(fHOutsub1sub2);
201   fOutputList->Add(fHOutNTEPRes);
202   fOutputList->Add(fHOutPTPsi);
203   fOutputList->Add(fHOutleadPTPsi);
204   fOutputList->Add(fHOutDiff);
205   
206   PostData(1, fOutputList); 
207   
208     if(!fuserphidist) { // if it's already set and custom class is required, we use the one provided by the user
209
210     TString oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
211
212     TFile foadb(oadbfilename); 
213     if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
214
215     AliInfo("Using Standard OADB");
216     fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");    
217     if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
218     foadb.Close();
219     }
220 }
221
222 //________________________________________________________________________
223 void AliEPSelectionTask::UserExec(Option_t */*option*/)
224
225   // Execute analysis for current event:
226   if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
227   
228 //   frunNumber = -15;
229  
230   AliEventplane* esdEP = 0;
231   TVector2 QQ1;
232   TVector2 QQ2;
233   Double_t fRP = 0.; // the monte carlo reaction plane angle
234     
235   if (fAnalysisInput.CompareTo("ESD")==0){
236
237     AliVEvent* event = InputEvent();
238     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
239     if (esd){    
240       if (!(frunNumber == esd->GetRunNumber())) {
241           frunNumber = esd->GetRunNumber();
242             SetPhiDist();      
243       }
244       
245       
246       if (fUseMCRP) {
247           AliMCEvent* mcEvent  = MCEvent();      
248           if (mcEvent && mcEvent->GenEventHeader()) {
249               AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
250               if (headerH) fRP = headerH->ReactionPlaneAngle();
251           }
252       }
253       
254       esdEP = esd->GetEventplane();
255       if (fSaveTrackContribution) {
256         esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
257         esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
258       }
259       
260       if (fTrackType.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
261       if (fTrackType.CompareTo("TPC")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
262       const int NT = ftracklist->GetEntries();
263       
264       if (NT>4){
265         fQVector = new TVector2(GetQ(esdEP,ftracklist));
266         fEventplaneQ = fQVector->Phi()/2; 
267         GetQsub(QQ1, QQ2, ftracklist);
268         fQsub1 = new TVector2(QQ1);
269         fQsub2 = new TVector2(QQ2);
270         fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
271         
272         esdEP->SetQVector(fQVector);
273         esdEP->SetEventplaneQ(fEventplaneQ);
274         esdEP->SetQsub(fQsub1,fQsub2);
275         esdEP->SetQsubRes(fQsubRes);
276         
277         fHOutEventplaneQ->Fill(fEventplaneQ);
278         fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
279         fHOutNTEPRes->Fill(NT,fQsubRes);
280
281         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
282         
283         for (int iter = 0; iter<NT;iter++){
284           AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
285           if (track) {
286             float delta = track->Phi()-fEventplaneQ;
287             while (delta < 0) delta += TMath::Pi();
288             while (delta > TMath::Pi()) delta -= TMath::Pi();
289             fHOutPTPsi->Fill(track->Pt(),delta);
290             fHOutPhi->Fill(track->Phi());
291             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
292           }
293         }
294         
295         AliESDtrack* trmax = esd->GetTrack(0);
296         for (int iter = 1; iter<NT;iter++){
297           AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
298           if (track && (track->Pt() > trmax->Pt())) trmax = track;
299         }
300         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
301       }
302       delete ftracklist;
303     }
304   }
305   
306   else if (fAnalysisInput.CompareTo("AOD")==0){
307     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
308     // to be implemented
309     printf("  AOD analysis not yet implemented!!!\n\n");
310     return;
311   }  
312   else {
313     printf(" Analysis Input not known!\n\n ");
314     return;
315   }
316   PostData(1, fOutputList); 
317 }
318
319 //________________________________________________________________________
320 void AliEPSelectionTask::Terminate(Option_t */*option*/)
321 {
322   // Terminate analysis
323 }
324
325 //__________________________________________________________________________
326 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
327 {
328   TVector2 mQ;
329   float mQx=0, mQy=0;
330   AliESDtrack* track;
331   Double_t weight;
332   
333   int NT = tracklist->GetEntries();
334
335   for (int i=0; i<NT; i++){
336     weight = 1;
337     track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
338     if (track) {
339       weight = GetWeight(track);
340       if (fSaveTrackContribution){
341         EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),track->GetID());
342         EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),track->GetID());
343       }
344       mQx += (weight*cos(2*track->Phi()));
345       mQy += (weight*sin(2*track->Phi()));
346     }
347   }
348   mQ.Set(mQx,mQy);
349   return mQ;
350 }
351   
352   //________________________________________________________________________
353 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
354 {
355   TVector2 mQ[2];
356   float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
357   Double_t weight;
358
359   AliESDtrack* track;
360   TRandom2 RN = 0;
361   
362   int NT = tracklist->GetEntries();
363   int trackcounter1=0, trackcounter2=0;
364   
365   for (Int_t i = 0; i < NT; i++) {
366     weight = 1;
367     track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
368     if (!track) continue;
369     weight = GetWeight(track);
370     
371     // This loop splits the track set into 2 random subsets
372     if( trackcounter1 < int(NT/2.) && trackcounter2 < int(NT/2.)){
373       float random = RN.Rndm();
374       if(random < .5){
375         mQx1 += (weight*cos(2*track->Phi()));
376         mQy1 += (weight*sin(2*track->Phi()));
377         trackcounter1++;
378       }
379       else {
380         mQx2 += (weight*cos(2*track->Phi()));
381         mQy2 += (weight*sin(2*track->Phi()));
382         trackcounter2++;
383       }
384     }
385     else if( trackcounter1 >= int(NT/2.)){
386       mQx2 += (weight*cos(2*track->Phi()));
387       mQy2 += (weight*sin(2*track->Phi()));
388       trackcounter2++;
389     }
390     else {
391       mQx1 += (weight*cos(2*track->Phi()));
392       mQy1 += (weight*sin(2*track->Phi()));
393       trackcounter1++;
394     }
395   }
396   mQ[0].Set(mQx1,mQy1);
397   mQ[1].Set(mQx2,mQy2);
398   Q1 = mQ[0];
399   Q2 = mQ[1];
400 }
401
402 //________________________________________________________________________
403 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
404   
405   fusercuts = kTRUE;
406   fESDtrackCuts = trackcuts;
407 }
408
409 //_____________________________________________________________________________
410 void AliEPSelectionTask::SetTrackType(TString tracktype){
411   fTrackType = tracktype;
412   if (!fusercuts) {
413   if (fTrackType.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
414   if (fTrackType.CompareTo("TPC")==0)    fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
415   fESDtrackCuts->SetPtRange(0.15,20);
416   fESDtrackCuts->SetEtaRange(-0.8,0.8);
417   }
418 }
419
420 //________________________________________________________________________
421 Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
422 {
423   Double_t ptweight=1;
424
425   if (fUsePtWeight) {      
426     if (track->Pt()<2) ptweight=track->Pt();
427     else ptweight=2;
428   }
429   return ptweight*GetPhiWeight(track);
430 }
431
432 //________________________________________________________________________
433 Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
434 {
435   Double_t phiweight=1;
436   
437   if (fUsePhiWeight && fPhiDist && track) {
438     Double_t nParticles = fPhiDist->Integral();
439     Double_t nPhibins = fPhiDist->GetNbinsX();
440   
441     Double_t Phi = track->Phi();
442     
443     while (Phi<0) Phi += TMath::TwoPi();
444     while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
445       
446     Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
447     
448     if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
449   }
450   return phiweight;
451 }
452
453 //__________________________________________________________________________
454 void AliEPSelectionTask::SetPhiDist() 
455 {
456   if(!fuserphidist) { // if it's already set and custom class is required, we use the one provided by the user
457
458     fPhiDist = (TH1F*) fEPContainer->GetObject(frunNumber, "Default");
459     if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", frunNumber));
460
461   } 
462   else {
463     AliInfo("Using Custom Phi Distribution");
464   }
465     
466   Bool_t emptybins;
467
468   int iter = 0;  
469   while (iter<3){
470       emptybins = kFALSE;
471    
472       for (int i=1; i<fPhiDist->GetNbinsX(); i++){
473         if (!((fPhiDist->GetBinContent(i))>0)) {
474           emptybins = kTRUE;
475         }
476       }  
477       if (emptybins) {
478         cout << "empty bins - rebinning!" << endl;
479         fPhiDist->Rebin();
480         iter++;
481       }      
482       else iter = 3;
483   }
484   
485   if (emptybins) {
486     AliError("After Maximum of rebinning still empty Phi-bins!!!");
487   }
488 }
489
490 //__________________________________________________________________________
491 void AliEPSelectionTask::SetPersonalPhiDistribution(char* infilename, char* listname)
492 {
493   
494   fuserphidist = kTRUE;
495   
496   TFile f(infilename);
497   TObject* list = f.Get(listname);
498   fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
499   if (!fPhiDist) AliFatal("Phi Distribution not found!!!");
500
501   f.Close();
502