cosmetics
[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(),
ce7adfe9 67 fAnalysisInput("ESD"),
90267ca6 68 fTrackType("TPC"),
ce7adfe9 69 fUseMCRP(kFALSE),
70 fUsePhiWeight(kFALSE),
71 fUsePtWeight(kFALSE),
72 fSaveTrackContribution(kFALSE),
90267ca6 73 fuserphidist(kFALSE),
74 fusercuts(kFALSE),
75 frunNumber(-15),
ce7adfe9 76 fESDtrackCuts(0),
77 ftracklist(0),
90267ca6 78 fEPContainer(0),
ce7adfe9 79 fPhiDist(0),
80 fQVector(0),
81 fQContributionX(0),
82 fQContributionY(0),
83 fEventplaneQ(0),
84 fQsub1(0),
85 fQsub2(0),
86 fQsubRes(0),
87 fOutputList(0),
88 fHOutEventplaneQ(0),
89 fHOutPhi(0),
90 fHOutPhiCorr(0),
91 fHOutsub1sub2(0),
92 fHOutNTEPRes(0),
93 fHOutPTPsi(0),
94 fHOutDiff(0),
95 fHOutleadPTPsi(0)
96{
97 // Default constructor
98 AliInfo("Event Plane Selection enabled.");
99}
100
101//________________________________________________________________________
102AliEPSelectionTask::AliEPSelectionTask(const char *name):
103 AliAnalysisTaskSE(name),
ce7adfe9 104 fAnalysisInput("ESD"),
90267ca6 105 fTrackType("TPC"),
ce7adfe9 106 fUseMCRP(kFALSE),
107 fUsePhiWeight(kFALSE),
108 fUsePtWeight(kFALSE),
109 fSaveTrackContribution(kFALSE),
90267ca6 110 fuserphidist(kFALSE),
111 fusercuts(kFALSE),
112 frunNumber(-15),
ce7adfe9 113 fESDtrackCuts(0),
114 ftracklist(0),
90267ca6 115 fEPContainer(0),
ce7adfe9 116 fPhiDist(0),
117 fQVector(0),
118 fQContributionX(0),
119 fQContributionY(0),
120 fEventplaneQ(0),
121 fQsub1(0),
122 fQsub2(0),
123 fQsubRes(0),
124 fOutputList(0),
125 fHOutEventplaneQ(0),
126 fHOutPhi(0),
127 fHOutPhiCorr(0),
128 fHOutsub1sub2(0),
129 fHOutNTEPRes(0),
130 fHOutPTPsi(0),
131 fHOutDiff(0),
132 fHOutleadPTPsi(0)
133{
134 // Default constructor
135 AliInfo("Event Plane Selection enabled.");
136 DefineOutput(1, TList::Class());
137}
138
139//________________________________________________________________________
140AliEPSelectionTask::~AliEPSelectionTask()
141{
142 // Destructor
143 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
144 delete fOutputList;
145 fOutputList = 0;
146 }
147 if (fESDtrackCuts){
148 delete fESDtrackCuts;
149 fESDtrackCuts = 0;
150 }
0b304ccd 151
90267ca6 152 if (fPhiDist){
153 delete fPhiDist;
154 fPhiDist = 0;
155 }
156 if (ftracklist){
157 delete ftracklist;
158 ftracklist = 0;
159 }
160 if (fEPContainer){
161 delete fEPContainer;
162 fEPContainer = 0;
163 }
ce7adfe9 164}
165
166//________________________________________________________________________
167void AliEPSelectionTask::UserCreateOutputObjects()
168{
169 // Create the output containers
170 if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
171 AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
172
173 fOutputList = new TList();
174 fOutputList->SetOwner();
175 fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
176 fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
177 fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
178 fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
179 fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
180 fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
181 fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
182 fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
183
184 fOutputList->Add(fHOutEventplaneQ);
185 fOutputList->Add(fHOutPhi);
186 fOutputList->Add(fHOutPhiCorr);
187 fOutputList->Add(fHOutsub1sub2);
188 fOutputList->Add(fHOutNTEPRes);
189 fOutputList->Add(fHOutPTPsi);
190 fOutputList->Add(fHOutleadPTPsi);
191 fOutputList->Add(fHOutDiff);
192
193 PostData(1, fOutputList);
90267ca6 194
195 if(!fuserphidist) { // if it's already set and custom class is required, we use the one provided by the user
196
197 TString oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
198
199 TFile foadb(oadbfilename);
200 if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
201
202 AliInfo("Using Standard OADB");
203 fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
204 if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
205 foadb.Close();
206 }
ce7adfe9 207}
208
209//________________________________________________________________________
210void AliEPSelectionTask::UserExec(Option_t */*option*/)
211{
212 // Execute analysis for current event:
213 if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
90267ca6 214
215// frunNumber = -15;
ce7adfe9 216
217 AliEventplane* esdEP = 0;
71916547 218 TVector2 qq1;
219 TVector2 qq2;
ce7adfe9 220 Double_t fRP = 0.; // the monte carlo reaction plane angle
221
222 if (fAnalysisInput.CompareTo("ESD")==0){
92955eca 223
ce7adfe9 224 AliVEvent* event = InputEvent();
225 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
92955eca 226 if (esd){
227 if (!(frunNumber == esd->GetRunNumber())) {
228 frunNumber = esd->GetRunNumber();
229 SetPhiDist();
ce7adfe9 230 }
92955eca 231
232
233 if (fUseMCRP) {
234 AliMCEvent* mcEvent = MCEvent();
235 if (mcEvent && mcEvent->GenEventHeader()) {
236 AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
237 if (headerH) fRP = headerH->ReactionPlaneAngle();
238 }
239 }
240
ce7adfe9 241 esdEP = esd->GetEventplane();
242 if (fSaveTrackContribution) {
243 esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
244 esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
245 }
246
90267ca6 247 if (fTrackType.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
248 if (fTrackType.CompareTo("TPC")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
71916547 249 const int nt = ftracklist->GetEntries();
ce7adfe9 250
71916547 251 if (nt>4){
ce7adfe9 252 fQVector = new TVector2(GetQ(esdEP,ftracklist));
253 fEventplaneQ = fQVector->Phi()/2;
71916547 254 GetQsub(qq1, qq2, ftracklist);
255 fQsub1 = new TVector2(qq1);
256 fQsub2 = new TVector2(qq2);
ce7adfe9 257 fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
90267ca6 258
ce7adfe9 259 esdEP->SetQVector(fQVector);
260 esdEP->SetEventplaneQ(fEventplaneQ);
261 esdEP->SetQsub(fQsub1,fQsub2);
262 esdEP->SetQsubRes(fQsubRes);
263
264 fHOutEventplaneQ->Fill(fEventplaneQ);
265 fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
71916547 266 fHOutNTEPRes->Fill(nt,fQsubRes);
ce7adfe9 267
268 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
269
71916547 270 for (int iter = 0; iter<nt;iter++){
ce7adfe9 271 AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
b1a983b4 272 if (track) {
273 float delta = track->Phi()-fEventplaneQ;
274 while (delta < 0) delta += TMath::Pi();
275 while (delta > TMath::Pi()) delta -= TMath::Pi();
276 fHOutPTPsi->Fill(track->Pt(),delta);
277 fHOutPhi->Fill(track->Phi());
278 fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
279 }
ce7adfe9 280 }
281
282 AliESDtrack* trmax = esd->GetTrack(0);
71916547 283 for (int iter = 1; iter<nt;iter++){
ce7adfe9 284 AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
b1a983b4 285 if (track && (track->Pt() > trmax->Pt())) trmax = track;
ce7adfe9 286 }
287 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
288 }
289 delete ftracklist;
290 }
291 }
292
293 else if (fAnalysisInput.CompareTo("AOD")==0){
294 //AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
295 // to be implemented
296 printf(" AOD analysis not yet implemented!!!\n\n");
297 return;
298 }
299 else {
300 printf(" Analysis Input not known!\n\n ");
301 return;
302 }
303 PostData(1, fOutputList);
304}
305
306//________________________________________________________________________
307void AliEPSelectionTask::Terminate(Option_t */*option*/)
308{
309 // Terminate analysis
310}
311
312//__________________________________________________________________________
313TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
314{
71916547 315// Get the Q vector
ce7adfe9 316 TVector2 mQ;
317 float mQx=0, mQy=0;
318 AliESDtrack* track;
319 Double_t weight;
320
71916547 321 int nt = tracklist->GetEntries();
ce7adfe9 322
71916547 323 for (int i=0; i<nt; i++){
ce7adfe9 324 weight = 1;
325 track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
b1a983b4 326 if (track) {
327 weight = GetWeight(track);
328 if (fSaveTrackContribution){
329 EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),track->GetID());
330 EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),track->GetID());
331 }
332 mQx += (weight*cos(2*track->Phi()));
333 mQy += (weight*sin(2*track->Phi()));
ce7adfe9 334 }
ce7adfe9 335 }
336 mQ.Set(mQx,mQy);
337 return mQ;
338}
339
340 //________________________________________________________________________
341void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
342{
71916547 343// Get Qsub
ce7adfe9 344 TVector2 mQ[2];
345 float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
346 Double_t weight;
347
348 AliESDtrack* track;
71916547 349 TRandom2 rn = 0;
ce7adfe9 350
71916547 351 int nt = tracklist->GetEntries();
ce7adfe9 352 int trackcounter1=0, trackcounter2=0;
353
71916547 354 for (Int_t i = 0; i < nt; i++) {
ce7adfe9 355 weight = 1;
356 track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
b1a983b4 357 if (!track) continue;
ce7adfe9 358 weight = GetWeight(track);
359
360 // This loop splits the track set into 2 random subsets
71916547 361 if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
362 float random = rn.Rndm();
ce7adfe9 363 if(random < .5){
364 mQx1 += (weight*cos(2*track->Phi()));
365 mQy1 += (weight*sin(2*track->Phi()));
366 trackcounter1++;
367 }
368 else {
369 mQx2 += (weight*cos(2*track->Phi()));
370 mQy2 += (weight*sin(2*track->Phi()));
371 trackcounter2++;
372 }
373 }
71916547 374 else if( trackcounter1 >= int(nt/2.)){
ce7adfe9 375 mQx2 += (weight*cos(2*track->Phi()));
376 mQy2 += (weight*sin(2*track->Phi()));
377 trackcounter2++;
378 }
379 else {
380 mQx1 += (weight*cos(2*track->Phi()));
381 mQy1 += (weight*sin(2*track->Phi()));
382 trackcounter1++;
383 }
384 }
385 mQ[0].Set(mQx1,mQy1);
386 mQ[1].Set(mQx2,mQy2);
387 Q1 = mQ[0];
388 Q2 = mQ[1];
389}
390
391//________________________________________________________________________
90267ca6 392void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
ce7adfe9 393
90267ca6 394 fusercuts = kTRUE;
395 fESDtrackCuts = trackcuts;
ce7adfe9 396}
397
90267ca6 398//_____________________________________________________________________________
399void AliEPSelectionTask::SetTrackType(TString tracktype){
71916547 400// Set the track type
90267ca6 401 fTrackType = tracktype;
402 if (!fusercuts) {
403 if (fTrackType.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
404 if (fTrackType.CompareTo("TPC")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
405 fESDtrackCuts->SetPtRange(0.15,20);
406 fESDtrackCuts->SetEtaRange(-0.8,0.8);
ce7adfe9 407 }
ce7adfe9 408}
409
410//________________________________________________________________________
411Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
412{
71916547 413// Get weight for track
ce7adfe9 414 Double_t ptweight=1;
415
416 if (fUsePtWeight) {
417 if (track->Pt()<2) ptweight=track->Pt();
418 else ptweight=2;
419 }
420 return ptweight*GetPhiWeight(track);
421}
422
423//________________________________________________________________________
424Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
425{
71916547 426// Get phi weight for track
ce7adfe9 427 Double_t phiweight=1;
428
e4f1ed0c 429 if (fUsePhiWeight && fPhiDist && track) {
ce7adfe9 430 Double_t nParticles = fPhiDist->Integral();
431 Double_t nPhibins = fPhiDist->GetNbinsX();
432
71916547 433 Double_t phi = track->Phi();
ce7adfe9 434
71916547 435 while (phi<0) phi += TMath::TwoPi();
436 while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
ce7adfe9 437
71916547 438 Double_t phiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
ce7adfe9 439
71916547 440 if (phiDistValue > 0) phiweight = nParticles/nPhibins/phiDistValue;
ce7adfe9 441 }
442 return phiweight;
443}
90267ca6 444
445//__________________________________________________________________________
446void AliEPSelectionTask::SetPhiDist()
447{
71916547 448// Set the phi distribution
90267ca6 449 if(!fuserphidist) { // if it's already set and custom class is required, we use the one provided by the user
450
451 fPhiDist = (TH1F*) fEPContainer->GetObject(frunNumber, "Default");
452 if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", frunNumber));
453
454 }
455 else {
456 AliInfo("Using Custom Phi Distribution");
457 }
458
459 Bool_t emptybins;
460
461 int iter = 0;
462 while (iter<3){
463 emptybins = kFALSE;
464
465 for (int i=1; i<fPhiDist->GetNbinsX(); i++){
466 if (!((fPhiDist->GetBinContent(i))>0)) {
467 emptybins = kTRUE;
468 }
469 }
470 if (emptybins) {
471 cout << "empty bins - rebinning!" << endl;
472 fPhiDist->Rebin();
473 iter++;
474 }
475 else iter = 3;
476 }
477
478 if (emptybins) {
479 AliError("After Maximum of rebinning still empty Phi-bins!!!");
480 }
481}
482
483//__________________________________________________________________________
71916547 484void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
90267ca6 485{
71916547 486 // Set a personal phi distribution
90267ca6 487 fuserphidist = kTRUE;
488
489 TFile f(infilename);
490 TObject* list = f.Get(listname);
491 fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
492 if (!fPhiDist) AliFatal("Phi Distribution not found!!!");
493
494 f.Close();
495}