Allow to pass trigger mask value as parameter to AddTaskQAsym
[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());
b1a983b4 208 if (headerH) fRP = headerH->ReactionPlaneAngle();
ce7adfe9 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));
b1a983b4 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 }
ce7adfe9 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));
b1a983b4 255 if (track && (track->Pt() > trmax->Pt())) trmax = track;
ce7adfe9 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//________________________________________________________________________
277void AliEPSelectionTask::Terminate(Option_t */*option*/)
278{
279 // Terminate analysis
280}
281
282//__________________________________________________________________________
283TVector2 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));
b1a983b4 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()));
ce7adfe9 303 }
ce7adfe9 304 }
305 mQ.Set(mQx,mQy);
306 return mQ;
307}
308
309 //________________________________________________________________________
310void 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));
b1a983b4 325 if (!track) continue;
ce7adfe9 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//________________________________________________________________________
360void 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//__________________________________________________________________________
08c34a44 370void AliEPSelectionTask::SetPhiDistribution(char* infilename, char* listname)
371{
ce7adfe9 372 TFile f(infilename);
373 TObject* list = f.Get(listname);
374 fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
e4f1ed0c 375 if (!fPhiDist) {
376 cout << "Phi Distribution not found!!!" << endl;
377 return;
378 }
ce7adfe9 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//________________________________________________________________________
406Double_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//________________________________________________________________________
418Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
419{
420 Double_t phiweight=1;
421
e4f1ed0c 422 if (fUsePhiWeight && fPhiDist && track) {
ce7adfe9 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
08c34a44 431 Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
ce7adfe9 432
433 if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
434 }
435 return phiweight;
436}