Coverity
[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         if (headerH) 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           if (track) {
243             float delta = track->Phi()-fEventplaneQ;
244             while (delta < 0) delta += TMath::Pi();
245             while (delta > TMath::Pi()) delta -= TMath::Pi();
246             fHOutPTPsi->Fill(track->Pt(),delta);
247             fHOutPhi->Fill(track->Phi());
248             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
249           }
250         }
251         
252         AliESDtrack* trmax = esd->GetTrack(0);
253         for (int iter = 1; iter<NT;iter++){
254           AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
255           if (track && (track->Pt() > trmax->Pt())) trmax = track;
256         }
257         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
258       }
259       delete ftracklist;
260     }
261   }
262   
263   else if (fAnalysisInput.CompareTo("AOD")==0){
264     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
265     // to be implemented
266     printf("  AOD analysis not yet implemented!!!\n\n");
267     return;
268   }  
269   else {
270     printf(" Analysis Input not known!\n\n ");
271     return;
272   }
273   PostData(1, fOutputList); 
274 }
275
276 //________________________________________________________________________
277 void AliEPSelectionTask::Terminate(Option_t */*option*/)
278 {
279   // Terminate analysis
280 }
281
282 //__________________________________________________________________________
283 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
284 {
285   TVector2 mQ;
286   float mQx=0, mQy=0;
287   AliESDtrack* track;
288   Double_t weight;
289   
290   int NT = tracklist->GetEntries();
291
292   for (int i=0; i<NT; i++){
293     weight = 1;
294     track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
295     if (track) {
296       weight = GetWeight(track);
297       if (fSaveTrackContribution){
298         EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),track->GetID());
299         EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),track->GetID());
300       }
301       mQx += (weight*cos(2*track->Phi()));
302       mQy += (weight*sin(2*track->Phi()));
303     }
304   }
305   mQ.Set(mQx,mQy);
306   return mQ;
307 }
308   
309   //________________________________________________________________________
310 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
311 {
312   TVector2 mQ[2];
313   float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
314   Double_t weight;
315
316   AliESDtrack* track;
317   TRandom2 RN = 0;
318   
319   int NT = tracklist->GetEntries();
320   int trackcounter1=0, trackcounter2=0;
321   
322   for (Int_t i = 0; i < NT; i++) {
323     weight = 1;
324     track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
325     if (!track) continue;
326     weight = GetWeight(track);
327     
328     // This loop splits the track set into 2 random subsets
329     if( trackcounter1 < int(NT/2.) && trackcounter2 < int(NT/2.)){
330       float random = RN.Rndm();
331       if(random < .5){
332         mQx1 += (weight*cos(2*track->Phi()));
333         mQy1 += (weight*sin(2*track->Phi()));
334         trackcounter1++;
335       }
336       else {
337         mQx2 += (weight*cos(2*track->Phi()));
338         mQy2 += (weight*sin(2*track->Phi()));
339         trackcounter2++;
340       }
341     }
342     else if( trackcounter1 >= int(NT/2.)){
343       mQx2 += (weight*cos(2*track->Phi()));
344       mQy2 += (weight*sin(2*track->Phi()));
345       trackcounter2++;
346     }
347     else {
348       mQx1 += (weight*cos(2*track->Phi()));
349       mQy1 += (weight*sin(2*track->Phi()));
350       trackcounter1++;
351     }
352   }
353   mQ[0].Set(mQx1,mQy1);
354   mQ[1].Set(mQx2,mQy2);
355   Q1 = mQ[0];
356   Q2 = mQ[1];
357 }
358
359 //________________________________________________________________________
360 void AliEPSelectionTask::SetESDtrackCuts(TString status){
361   
362   fStatus = status;
363   if (fStatus.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
364   if (fStatus.CompareTo("TPC")==0)    fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
365   fESDtrackCuts->SetPtRange(0.15,20);
366   fESDtrackCuts->SetEtaRange(-0.8,0.8);
367 }
368
369 //__________________________________________________________________________
370 void AliEPSelectionTask::SetPhiDistribution(char* infilename, char* listname)
371 {
372   TFile f(infilename);
373   TObject* list = f.Get(listname);
374   fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
375   if (!fPhiDist) {
376     cout << "Phi Distribution not found!!!" << endl;
377     return;
378   }
379   
380   Bool_t emptybins;
381
382   int iter = 0;  
383   while (iter<3){
384       emptybins = kFALSE;
385    
386       for (int i=1; i<fPhiDist->GetNbinsX(); i++){
387         if (!((fPhiDist->GetBinContent(i))>0)) {
388           emptybins = kTRUE;
389         }
390       }  
391       if (emptybins) {
392         cout << "empty bins - rebinning!" << endl;
393         fPhiDist->Rebin();
394         iter++;
395       }      
396       else iter = 3;
397   }
398   
399   if (emptybins) {
400     AliError("After Maximum of rebinning still empty Phi-bins!!!");
401   }
402   f.Close();
403 }
404
405 //________________________________________________________________________
406 Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
407 {
408   Double_t ptweight=1;
409
410   if (fUsePtWeight) {      
411     if (track->Pt()<2) ptweight=track->Pt();
412     else ptweight=2;
413   }
414   return ptweight*GetPhiWeight(track);
415 }
416
417 //________________________________________________________________________
418 Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
419 {
420   Double_t phiweight=1;
421   
422   if (fUsePhiWeight && fPhiDist && track) {
423     Double_t nParticles = fPhiDist->Integral();
424     Double_t nPhibins = fPhiDist->GetNbinsX();
425   
426     Double_t Phi = track->Phi();
427     
428     while (Phi<0) Phi += TMath::TwoPi();
429     while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
430       
431     Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
432     
433     if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
434   }
435   return phiweight;
436 }