Update mini-task with big output flag and add current status macros
[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){
236
237 AliVEvent* event = InputEvent();
238 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
90267ca6 239 if (!(frunNumber == esd->GetRunNumber())) {
240 frunNumber = esd->GetRunNumber();
241 SetPhiDist();
242 }
243
ce7adfe9 244
245 if (fUseMCRP) {
246 AliMCEvent* mcEvent = MCEvent();
247 if (mcEvent && mcEvent->GenEventHeader()) {
248 AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
b1a983b4 249 if (headerH) fRP = headerH->ReactionPlaneAngle();
ce7adfe9 250 }
251 }
252 if (esd){
253 esdEP = esd->GetEventplane();
254 if (fSaveTrackContribution) {
255 esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
256 esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
257 }
258
90267ca6 259 if (fTrackType.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
260 if (fTrackType.CompareTo("TPC")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
ce7adfe9 261 const int NT = ftracklist->GetEntries();
262
263 if (NT>4){
264 fQVector = new TVector2(GetQ(esdEP,ftracklist));
265 fEventplaneQ = fQVector->Phi()/2;
266 GetQsub(QQ1, QQ2, ftracklist);
267 fQsub1 = new TVector2(QQ1);
268 fQsub2 = new TVector2(QQ2);
269 fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
90267ca6 270
ce7adfe9 271 esdEP->SetQVector(fQVector);
272 esdEP->SetEventplaneQ(fEventplaneQ);
273 esdEP->SetQsub(fQsub1,fQsub2);
274 esdEP->SetQsubRes(fQsubRes);
275
276 fHOutEventplaneQ->Fill(fEventplaneQ);
277 fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
278 fHOutNTEPRes->Fill(NT,fQsubRes);
279
280 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
281
282 for (int iter = 0; iter<NT;iter++){
283 AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
b1a983b4 284 if (track) {
285 float delta = track->Phi()-fEventplaneQ;
286 while (delta < 0) delta += TMath::Pi();
287 while (delta > TMath::Pi()) delta -= TMath::Pi();
288 fHOutPTPsi->Fill(track->Pt(),delta);
289 fHOutPhi->Fill(track->Phi());
290 fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
291 }
ce7adfe9 292 }
293
294 AliESDtrack* trmax = esd->GetTrack(0);
295 for (int iter = 1; iter<NT;iter++){
296 AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
b1a983b4 297 if (track && (track->Pt() > trmax->Pt())) trmax = track;
ce7adfe9 298 }
299 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
300 }
301 delete ftracklist;
302 }
303 }
304
305 else if (fAnalysisInput.CompareTo("AOD")==0){
306 //AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
307 // to be implemented
308 printf(" AOD analysis not yet implemented!!!\n\n");
309 return;
310 }
311 else {
312 printf(" Analysis Input not known!\n\n ");
313 return;
314 }
315 PostData(1, fOutputList);
316}
317
318//________________________________________________________________________
319void AliEPSelectionTask::Terminate(Option_t */*option*/)
320{
321 // Terminate analysis
322}
323
324//__________________________________________________________________________
325TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
326{
327 TVector2 mQ;
328 float mQx=0, mQy=0;
329 AliESDtrack* track;
330 Double_t weight;
331
332 int NT = tracklist->GetEntries();
333
334 for (int i=0; i<NT; i++){
335 weight = 1;
336 track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
b1a983b4 337 if (track) {
338 weight = GetWeight(track);
339 if (fSaveTrackContribution){
340 EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),track->GetID());
341 EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),track->GetID());
342 }
343 mQx += (weight*cos(2*track->Phi()));
344 mQy += (weight*sin(2*track->Phi()));
ce7adfe9 345 }
ce7adfe9 346 }
347 mQ.Set(mQx,mQy);
348 return mQ;
349}
350
351 //________________________________________________________________________
352void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
353{
354 TVector2 mQ[2];
355 float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
356 Double_t weight;
357
358 AliESDtrack* track;
359 TRandom2 RN = 0;
360
361 int NT = tracklist->GetEntries();
362 int trackcounter1=0, trackcounter2=0;
363
364 for (Int_t i = 0; i < NT; i++) {
365 weight = 1;
366 track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
b1a983b4 367 if (!track) continue;
ce7adfe9 368 weight = GetWeight(track);
369
370 // This loop splits the track set into 2 random subsets
371 if( trackcounter1 < int(NT/2.) && trackcounter2 < int(NT/2.)){
372 float random = RN.Rndm();
373 if(random < .5){
374 mQx1 += (weight*cos(2*track->Phi()));
375 mQy1 += (weight*sin(2*track->Phi()));
376 trackcounter1++;
377 }
378 else {
379 mQx2 += (weight*cos(2*track->Phi()));
380 mQy2 += (weight*sin(2*track->Phi()));
381 trackcounter2++;
382 }
383 }
384 else if( trackcounter1 >= int(NT/2.)){
385 mQx2 += (weight*cos(2*track->Phi()));
386 mQy2 += (weight*sin(2*track->Phi()));
387 trackcounter2++;
388 }
389 else {
390 mQx1 += (weight*cos(2*track->Phi()));
391 mQy1 += (weight*sin(2*track->Phi()));
392 trackcounter1++;
393 }
394 }
395 mQ[0].Set(mQx1,mQy1);
396 mQ[1].Set(mQx2,mQy2);
397 Q1 = mQ[0];
398 Q2 = mQ[1];
399}
400
401//________________________________________________________________________
90267ca6 402void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
ce7adfe9 403
90267ca6 404 fusercuts = kTRUE;
405 fESDtrackCuts = trackcuts;
ce7adfe9 406}
407
90267ca6 408//_____________________________________________________________________________
409void AliEPSelectionTask::SetTrackType(TString tracktype){
410 fTrackType = tracktype;
411 if (!fusercuts) {
412 if (fTrackType.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
413 if (fTrackType.CompareTo("TPC")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
414 fESDtrackCuts->SetPtRange(0.15,20);
415 fESDtrackCuts->SetEtaRange(-0.8,0.8);
ce7adfe9 416 }
ce7adfe9 417}
418
419//________________________________________________________________________
420Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
421{
422 Double_t ptweight=1;
423
424 if (fUsePtWeight) {
425 if (track->Pt()<2) ptweight=track->Pt();
426 else ptweight=2;
427 }
428 return ptweight*GetPhiWeight(track);
429}
430
431//________________________________________________________________________
432Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
433{
434 Double_t phiweight=1;
435
e4f1ed0c 436 if (fUsePhiWeight && fPhiDist && track) {
ce7adfe9 437 Double_t nParticles = fPhiDist->Integral();
438 Double_t nPhibins = fPhiDist->GetNbinsX();
439
440 Double_t Phi = track->Phi();
441
442 while (Phi<0) Phi += TMath::TwoPi();
443 while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
444
08c34a44 445 Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
ce7adfe9 446
447 if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
448 }
449 return phiweight;
450}
90267ca6 451
452//__________________________________________________________________________
453void AliEPSelectionTask::SetPhiDist()
454{
455 if(!fuserphidist) { // if it's already set and custom class is required, we use the one provided by the user
456
457 fPhiDist = (TH1F*) fEPContainer->GetObject(frunNumber, "Default");
458 if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", frunNumber));
459
460 }
461 else {
462 AliInfo("Using Custom Phi Distribution");
463 }
464
465 Bool_t emptybins;
466
467 int iter = 0;
468 while (iter<3){
469 emptybins = kFALSE;
470
471 for (int i=1; i<fPhiDist->GetNbinsX(); i++){
472 if (!((fPhiDist->GetBinContent(i))>0)) {
473 emptybins = kTRUE;
474 }
475 }
476 if (emptybins) {
477 cout << "empty bins - rebinning!" << endl;
478 fPhiDist->Rebin();
479 iter++;
480 }
481 else iter = 3;
482 }
483
484 if (emptybins) {
485 AliError("After Maximum of rebinning still empty Phi-bins!!!");
486 }
487}
488
489//__________________________________________________________________________
490void AliEPSelectionTask::SetPersonalPhiDistribution(char* infilename, char* listname)
491{
492
493 fuserphidist = kTRUE;
494
495 TFile f(infilename);
496 TObject* list = f.Get(listname);
497 fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
498 if (!fPhiDist) AliFatal("Phi Distribution not found!!!");
499
500 f.Close();
501}