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