updated track cut systematic study
[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"
90267ca6 58#include "AliOADBContainer.h"
ce7adfe9 59
60#include "AliEventplane.h"
61
62ClassImp(AliEPSelectionTask)
63
64//________________________________________________________________________
65AliEPSelectionTask::AliEPSelectionTask():
66AliAnalysisTaskSE(),
67 fDebug(0),
68 fAnalysisInput("ESD"),
90267ca6 69 fTrackType("TPC"),
ce7adfe9 70 fUseMCRP(kFALSE),
71 fUsePhiWeight(kFALSE),
72 fUsePtWeight(kFALSE),
73 fSaveTrackContribution(kFALSE),
90267ca6 74 fuserphidist(kFALSE),
75 fusercuts(kFALSE),
76 frunNumber(-15),
ce7adfe9 77 fESDtrackCuts(0),
78 ftracklist(0),
90267ca6 79 fEPContainer(0),
ce7adfe9 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//________________________________________________________________________
103AliEPSelectionTask::AliEPSelectionTask(const char *name):
104 AliAnalysisTaskSE(name),
105 fDebug(0),
106 fAnalysisInput("ESD"),
90267ca6 107 fTrackType("TPC"),
ce7adfe9 108 fUseMCRP(kFALSE),
109 fUsePhiWeight(kFALSE),
110 fUsePtWeight(kFALSE),
111 fSaveTrackContribution(kFALSE),
90267ca6 112 fuserphidist(kFALSE),
113 fusercuts(kFALSE),
114 frunNumber(-15),
ce7adfe9 115 fESDtrackCuts(0),
116 ftracklist(0),
90267ca6 117 fEPContainer(0),
ce7adfe9 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//________________________________________________________________________
142AliEPSelectionTask::~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 }
90267ca6 161 if (fQsub2){
ce7adfe9 162 delete fQsub2;
163 fQsub2 = 0;
164 }
90267ca6 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 }
ce7adfe9 177}
178
179//________________________________________________________________________
180void 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);
90267ca6 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 }
ce7adfe9 220}
221
222//________________________________________________________________________
223void AliEPSelectionTask::UserExec(Option_t */*option*/)
224{
225 // Execute analysis for current event:
226 if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
90267ca6 227
228// frunNumber = -15;
ce7adfe9 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){
92955eca 236
ce7adfe9 237 AliVEvent* event = InputEvent();
238 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
92955eca 239 if (esd){
240 if (!(frunNumber == esd->GetRunNumber())) {
241 frunNumber = esd->GetRunNumber();
242 SetPhiDist();
ce7adfe9 243 }
92955eca 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
ce7adfe9 254 esdEP = esd->GetEventplane();
255 if (fSaveTrackContribution) {
256 esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
257 esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
258 }
259
90267ca6 260 if (fTrackType.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
261 if (fTrackType.CompareTo("TPC")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
ce7adfe9 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);
90267ca6 271
ce7adfe9 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));
b1a983b4 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 }
ce7adfe9 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));
b1a983b4 298 if (track && (track->Pt() > trmax->Pt())) trmax = track;
ce7adfe9 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//________________________________________________________________________
320void AliEPSelectionTask::Terminate(Option_t */*option*/)
321{
322 // Terminate analysis
323}
324
325//__________________________________________________________________________
326TVector2 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));
b1a983b4 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()));
ce7adfe9 346 }
ce7adfe9 347 }
348 mQ.Set(mQx,mQy);
349 return mQ;
350}
351
352 //________________________________________________________________________
353void 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));
b1a983b4 368 if (!track) continue;
ce7adfe9 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//________________________________________________________________________
90267ca6 403void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
ce7adfe9 404
90267ca6 405 fusercuts = kTRUE;
406 fESDtrackCuts = trackcuts;
ce7adfe9 407}
408
90267ca6 409//_____________________________________________________________________________
410void 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);
ce7adfe9 417 }
ce7adfe9 418}
419
420//________________________________________________________________________
421Double_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//________________________________________________________________________
433Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
434{
435 Double_t phiweight=1;
436
e4f1ed0c 437 if (fUsePhiWeight && fPhiDist && track) {
ce7adfe9 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
08c34a44 446 Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
ce7adfe9 447
448 if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
449 }
450 return phiweight;
451}
90267ca6 452
453//__________________________________________________________________________
454void 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//__________________________________________________________________________
491void 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}