Coverity 16550
[u/mrichter/AliRoot.git] / ANALYSIS / AliEPSelectionTask.cxx
CommitLineData
ce7adfe9 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
61ClassImp(AliEPSelectionTask)
62
63//________________________________________________________________________
64AliEPSelectionTask::AliEPSelectionTask():
65AliAnalysisTaskSE(),
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//________________________________________________________________________
98AliEPSelectionTask::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//________________________________________________________________________
133AliEPSelectionTask::~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//________________________________________________________________________
159void 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//________________________________________________________________________
189void 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);
ce7adfe9 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//________________________________________________________________________
275void AliEPSelectionTask::Terminate(Option_t */*option*/)
276{
277 // Terminate analysis
278}
279
280//__________________________________________________________________________
281TVector2 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 //________________________________________________________________________
306void 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//________________________________________________________________________
355void 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//__________________________________________________________________________
08c34a44 365void AliEPSelectionTask::SetPhiDistribution(char* infilename, char* listname)
366{
ce7adfe9 367 TFile f(infilename);
368 TObject* list = f.Get(listname);
369 fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
e4f1ed0c 370 if (!fPhiDist) {
371 cout << "Phi Distribution not found!!!" << endl;
372 return;
373 }
ce7adfe9 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//________________________________________________________________________
401Double_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//________________________________________________________________________
413Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
414{
415 Double_t phiweight=1;
416
e4f1ed0c 417 if (fUsePhiWeight && fPhiDist && track) {
ce7adfe9 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
08c34a44 426 Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
ce7adfe9 427
428 if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
429 }
430 return phiweight;
431}