Coverity 16550
[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
59 #include "AliEventplane.h"
60
61 ClassImp(AliEPSelectionTask)
62
63 //________________________________________________________________________
64 AliEPSelectionTask::AliEPSelectionTask():
65 AliAnalysisTaskSE(),
66   fDebug(0),
67   fAnalysisInput("ESD"),
68   fStatus("TPC"),
69   fUseMCRP(kFALSE),
70   fUsePhiWeight(kFALSE),
71   fUsePtWeight(kFALSE),
72   fSaveTrackContribution(kFALSE),
73   fESDtrackCuts(0),
74   ftracklist(0),
75   fPhiDist(0),
76   fQVector(0),
77   fQContributionX(0),
78   fQContributionY(0),
79   fEventplaneQ(0),
80   fQsub1(0),
81   fQsub2(0),
82   fQsubRes(0),  
83   fOutputList(0),
84   fHOutEventplaneQ(0),
85   fHOutPhi(0),
86   fHOutPhiCorr(0),
87   fHOutsub1sub2(0),
88   fHOutNTEPRes(0),
89   fHOutPTPsi(0),
90   fHOutDiff(0),
91   fHOutleadPTPsi(0)
92 {   
93   // Default constructor
94   AliInfo("Event Plane Selection enabled.");
95 }   
96
97 //________________________________________________________________________
98 AliEPSelectionTask::AliEPSelectionTask(const char *name):
99   AliAnalysisTaskSE(name),
100   fDebug(0),
101   fAnalysisInput("ESD"),
102   fStatus("TPC"),
103   fUseMCRP(kFALSE),
104   fUsePhiWeight(kFALSE),
105   fUsePtWeight(kFALSE),  
106   fSaveTrackContribution(kFALSE),
107   fESDtrackCuts(0),
108   ftracklist(0),
109   fPhiDist(0),
110   fQVector(0),
111   fQContributionX(0),
112   fQContributionY(0),
113   fEventplaneQ(0),
114   fQsub1(0),
115   fQsub2(0),
116   fQsubRes(0),
117   fOutputList(0),
118   fHOutEventplaneQ(0),
119   fHOutPhi(0),
120   fHOutPhiCorr(0),
121   fHOutsub1sub2(0),
122   fHOutNTEPRes(0),
123   fHOutPTPsi(0),
124   fHOutDiff(0),
125   fHOutleadPTPsi(0)
126 {
127   // Default constructor
128   AliInfo("Event Plane Selection enabled.");
129   DefineOutput(1, TList::Class());
130 }
131  
132 //________________________________________________________________________
133 AliEPSelectionTask::~AliEPSelectionTask()
134 {
135   // Destructor  
136   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
137       delete fOutputList;
138       fOutputList = 0;
139   }
140   if (fESDtrackCuts){
141       delete fESDtrackCuts;
142       fESDtrackCuts = 0;
143   }
144   if (fQVector){
145       delete fQVector;
146       fQVector = 0;
147   }
148   if (fQsub1){
149       delete fQsub1;
150       fQsub1 = 0;
151   }
152     if (fQsub2){
153       delete fQsub2;
154       fQsub2 = 0;
155   }
156 }  
157
158 //________________________________________________________________________
159 void AliEPSelectionTask::UserCreateOutputObjects()
160 {  
161   // Create the output containers
162   if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
163   AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
164
165   fOutputList = new TList();
166   fOutputList->SetOwner();
167   fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
168   fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
169   fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
170   fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
171   fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
172   fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
173   fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
174   fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
175
176   fOutputList->Add(fHOutEventplaneQ);
177   fOutputList->Add(fHOutPhi);
178   fOutputList->Add(fHOutPhiCorr);
179   fOutputList->Add(fHOutsub1sub2);
180   fOutputList->Add(fHOutNTEPRes);
181   fOutputList->Add(fHOutPTPsi);
182   fOutputList->Add(fHOutleadPTPsi);
183   fOutputList->Add(fHOutDiff);
184   
185   PostData(1, fOutputList); 
186 }
187
188 //________________________________________________________________________
189 void AliEPSelectionTask::UserExec(Option_t */*option*/)
190
191   // Execute analysis for current event:
192   if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
193  
194   AliEventplane* esdEP = 0;
195   TVector2 QQ1;
196   TVector2 QQ2;
197   Double_t fRP = 0.; // the monte carlo reaction plane angle
198     
199   if (fAnalysisInput.CompareTo("ESD")==0){
200     
201     AliVEvent* event = InputEvent();
202     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
203     
204     if (fUseMCRP) {
205       AliMCEvent* mcEvent  = MCEvent();      
206       if (mcEvent && mcEvent->GenEventHeader()) {
207         AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
208         fRP = headerH->ReactionPlaneAngle();
209       }
210     }
211     if (esd){
212       esdEP = esd->GetEventplane();
213       if (fSaveTrackContribution) {
214         esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
215         esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
216       }
217       
218       if (fStatus.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
219       if (fStatus.CompareTo("TPC")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
220       const int NT = ftracklist->GetEntries();
221       
222       if (NT>4){
223         fQVector = new TVector2(GetQ(esdEP,ftracklist));
224         fEventplaneQ = fQVector->Phi()/2; 
225         GetQsub(QQ1, QQ2, ftracklist);
226         fQsub1 = new TVector2(QQ1);
227         fQsub2 = new TVector2(QQ2);
228         fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
229         esdEP->SetQVector(fQVector);
230         esdEP->SetEventplaneQ(fEventplaneQ);
231         esdEP->SetQsub(fQsub1,fQsub2);
232         esdEP->SetQsubRes(fQsubRes);
233         
234         fHOutEventplaneQ->Fill(fEventplaneQ);
235         fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
236         fHOutNTEPRes->Fill(NT,fQsubRes);
237
238         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
239         
240         for (int iter = 0; iter<NT;iter++){
241           AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
242           float delta = track->Phi()-fEventplaneQ;
243           while (delta < 0) delta += TMath::Pi();
244           while (delta > TMath::Pi()) delta -= TMath::Pi();
245           fHOutPTPsi->Fill(track->Pt(),delta);
246           fHOutPhi->Fill(track->Phi());
247           fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
248         }
249         
250         AliESDtrack* trmax = esd->GetTrack(0);
251         for (int iter = 1; iter<NT;iter++){
252           AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
253           if (track->Pt() > trmax->Pt()) trmax = track;
254         }
255         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
256       }
257       delete ftracklist;
258     }
259   }
260   
261   else if (fAnalysisInput.CompareTo("AOD")==0){
262     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
263     // to be implemented
264     printf("  AOD analysis not yet implemented!!!\n\n");
265     return;
266   }  
267   else {
268     printf(" Analysis Input not known!\n\n ");
269     return;
270   }
271   PostData(1, fOutputList); 
272 }
273
274 //________________________________________________________________________
275 void AliEPSelectionTask::Terminate(Option_t */*option*/)
276 {
277   // Terminate analysis
278 }
279
280 //__________________________________________________________________________
281 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
282 {
283   TVector2 mQ;
284   float mQx=0, mQy=0;
285   AliESDtrack* track;
286   Double_t weight;
287   
288   int NT = tracklist->GetEntries();
289
290   for (int i=0; i<NT; i++){
291     weight = 1;
292     track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
293     weight = GetWeight(track);
294     if (fSaveTrackContribution){
295       EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),track->GetID());
296       EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),track->GetID());
297     }
298     mQx += (weight*cos(2*track->Phi()));
299     mQy += (weight*sin(2*track->Phi()));
300   }
301   mQ.Set(mQx,mQy);
302   return mQ;
303 }
304   
305   //________________________________________________________________________
306 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
307 {
308   TVector2 mQ[2];
309   float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
310   Double_t weight;
311
312   AliESDtrack* track;
313   TRandom2 RN = 0;
314   
315   int NT = tracklist->GetEntries();
316   int trackcounter1=0, trackcounter2=0;
317   
318   for (Int_t i = 0; i < NT; i++) {
319     weight = 1;
320     track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
321     weight = GetWeight(track);
322     
323     // This loop splits the track set into 2 random subsets
324     if( trackcounter1 < int(NT/2.) && trackcounter2 < int(NT/2.)){
325       float random = RN.Rndm();
326       if(random < .5){
327         mQx1 += (weight*cos(2*track->Phi()));
328         mQy1 += (weight*sin(2*track->Phi()));
329         trackcounter1++;
330       }
331       else {
332         mQx2 += (weight*cos(2*track->Phi()));
333         mQy2 += (weight*sin(2*track->Phi()));
334         trackcounter2++;
335       }
336     }
337     else if( trackcounter1 >= int(NT/2.)){
338       mQx2 += (weight*cos(2*track->Phi()));
339       mQy2 += (weight*sin(2*track->Phi()));
340       trackcounter2++;
341     }
342     else {
343       mQx1 += (weight*cos(2*track->Phi()));
344       mQy1 += (weight*sin(2*track->Phi()));
345       trackcounter1++;
346     }
347   }
348   mQ[0].Set(mQx1,mQy1);
349   mQ[1].Set(mQx2,mQy2);
350   Q1 = mQ[0];
351   Q2 = mQ[1];
352 }
353
354 //________________________________________________________________________
355 void AliEPSelectionTask::SetESDtrackCuts(TString status){
356   
357   fStatus = status;
358   if (fStatus.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
359   if (fStatus.CompareTo("TPC")==0)    fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
360   fESDtrackCuts->SetPtRange(0.15,20);
361   fESDtrackCuts->SetEtaRange(-0.8,0.8);
362 }
363
364 //__________________________________________________________________________
365 void AliEPSelectionTask::SetPhiDistribution(char* infilename, char* listname)
366 {
367   TFile f(infilename);
368   TObject* list = f.Get(listname);
369   fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
370   if (!fPhiDist) {
371     cout << "Phi Distribution not found!!!" << endl;
372     return;
373   }
374   
375   Bool_t emptybins;
376
377   int iter = 0;  
378   while (iter<3){
379       emptybins = kFALSE;
380    
381       for (int i=1; i<fPhiDist->GetNbinsX(); i++){
382         if (!((fPhiDist->GetBinContent(i))>0)) {
383           emptybins = kTRUE;
384         }
385       }  
386       if (emptybins) {
387         cout << "empty bins - rebinning!" << endl;
388         fPhiDist->Rebin();
389         iter++;
390       }      
391       else iter = 3;
392   }
393   
394   if (emptybins) {
395     AliError("After Maximum of rebinning still empty Phi-bins!!!");
396   }
397   f.Close();
398 }
399
400 //________________________________________________________________________
401 Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
402 {
403   Double_t ptweight=1;
404
405   if (fUsePtWeight) {      
406     if (track->Pt()<2) ptweight=track->Pt();
407     else ptweight=2;
408   }
409   return ptweight*GetPhiWeight(track);
410 }
411
412 //________________________________________________________________________
413 Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
414 {
415   Double_t phiweight=1;
416   
417   if (fUsePhiWeight && fPhiDist && track) {
418     Double_t nParticles = fPhiDist->Integral();
419     Double_t nPhibins = fPhiDist->GetNbinsX();
420   
421     Double_t Phi = track->Phi();
422     
423     while (Phi<0) Phi += TMath::TwoPi();
424     while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
425       
426     Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
427     
428     if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
429   }
430   return phiweight;
431 }