track iso fix
[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"
e51055a0 59#include "AliAODMCHeader.h"
60#include "AliAODTrack.h"
61#include "AliVTrack.h"
ce7adfe9 62#include "AliEventplane.h"
63
64ClassImp(AliEPSelectionTask)
65
66//________________________________________________________________________
67AliEPSelectionTask::AliEPSelectionTask():
68AliAnalysisTaskSE(),
ce7adfe9 69 fAnalysisInput("ESD"),
90267ca6 70 fTrackType("TPC"),
ce7adfe9 71 fUseMCRP(kFALSE),
72 fUsePhiWeight(kFALSE),
73 fUsePtWeight(kFALSE),
74 fSaveTrackContribution(kFALSE),
1e552991 75 fUserphidist(kFALSE),
76 fUsercuts(kFALSE),
77 fRunNumber(-15),
e51055a0 78 fAODfilterbit(1),
79 fEtaGap(0.),
80 fSplitMethod(0),
ce7adfe9 81 fESDtrackCuts(0),
90267ca6 82 fEPContainer(0),
ce7adfe9 83 fPhiDist(0),
84 fQVector(0),
85 fQContributionX(0),
86 fQContributionY(0),
87 fEventplaneQ(0),
88 fQsub1(0),
89 fQsub2(0),
90 fQsubRes(0),
91 fOutputList(0),
92 fHOutEventplaneQ(0),
93 fHOutPhi(0),
94 fHOutPhiCorr(0),
95 fHOutsub1sub2(0),
96 fHOutNTEPRes(0),
97 fHOutPTPsi(0),
98 fHOutDiff(0),
99 fHOutleadPTPsi(0)
100{
101 // Default constructor
102 AliInfo("Event Plane Selection enabled.");
103}
104
105//________________________________________________________________________
106AliEPSelectionTask::AliEPSelectionTask(const char *name):
107 AliAnalysisTaskSE(name),
ce7adfe9 108 fAnalysisInput("ESD"),
90267ca6 109 fTrackType("TPC"),
ce7adfe9 110 fUseMCRP(kFALSE),
111 fUsePhiWeight(kFALSE),
112 fUsePtWeight(kFALSE),
113 fSaveTrackContribution(kFALSE),
1e552991 114 fUserphidist(kFALSE),
115 fUsercuts(kFALSE),
116 fRunNumber(-15),
e51055a0 117 fAODfilterbit(1),
118 fEtaGap(0.),
119 fSplitMethod(0),
ce7adfe9 120 fESDtrackCuts(0),
90267ca6 121 fEPContainer(0),
ce7adfe9 122 fPhiDist(0),
123 fQVector(0),
124 fQContributionX(0),
125 fQContributionY(0),
126 fEventplaneQ(0),
127 fQsub1(0),
128 fQsub2(0),
129 fQsubRes(0),
130 fOutputList(0),
131 fHOutEventplaneQ(0),
132 fHOutPhi(0),
133 fHOutPhiCorr(0),
134 fHOutsub1sub2(0),
135 fHOutNTEPRes(0),
136 fHOutPTPsi(0),
137 fHOutDiff(0),
138 fHOutleadPTPsi(0)
139{
140 // Default constructor
141 AliInfo("Event Plane Selection enabled.");
142 DefineOutput(1, TList::Class());
143}
144
145//________________________________________________________________________
146AliEPSelectionTask::~AliEPSelectionTask()
147{
148 // Destructor
149 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
150 delete fOutputList;
151 fOutputList = 0;
152 }
153 if (fESDtrackCuts){
154 delete fESDtrackCuts;
155 fESDtrackCuts = 0;
156 }
1e552991 157 if (fUserphidist) {
158 if (fPhiDist) {
90267ca6 159 delete fPhiDist;
160 fPhiDist = 0;
1e552991 161 }
90267ca6 162 }
163 if (fEPContainer){
164 delete fEPContainer;
165 fEPContainer = 0;
166 }
ce7adfe9 167}
168
169//________________________________________________________________________
170void AliEPSelectionTask::UserCreateOutputObjects()
171{
172 // Create the output containers
173 if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
174 AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
175
176 fOutputList = new TList();
177 fOutputList->SetOwner();
178 fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
179 fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
180 fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
181 fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
182 fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
183 fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
184 fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
185 fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
186
187 fOutputList->Add(fHOutEventplaneQ);
188 fOutputList->Add(fHOutPhi);
189 fOutputList->Add(fHOutPhiCorr);
190 fOutputList->Add(fHOutsub1sub2);
191 fOutputList->Add(fHOutNTEPRes);
192 fOutputList->Add(fHOutPTPsi);
193 fOutputList->Add(fHOutleadPTPsi);
194 fOutputList->Add(fHOutDiff);
195
196 PostData(1, fOutputList);
90267ca6 197
1e552991 198 if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
90267ca6 199
e51055a0 200
201 TString oadbfilename;
90267ca6 202
e51055a0 203 if (fAnalysisInput.CompareTo("AOD")==0){
204 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
205 } else if (fAnalysisInput.CompareTo("ESD")==0){
206 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
207 }
208
90267ca6 209 TFile foadb(oadbfilename);
210 if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
211
212 AliInfo("Using Standard OADB");
213 fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
214 if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
215 foadb.Close();
216 }
e51055a0 217
ce7adfe9 218}
219
220//________________________________________________________________________
221void AliEPSelectionTask::UserExec(Option_t */*option*/)
222{
223 // Execute analysis for current event:
224 if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
90267ca6 225
1e552991 226// fRunNumber = -15;
ce7adfe9 227
e51055a0 228 AliEventplane *esdEP;
71916547 229 TVector2 qq1;
230 TVector2 qq2;
ce7adfe9 231 Double_t fRP = 0.; // the monte carlo reaction plane angle
232
233 if (fAnalysisInput.CompareTo("ESD")==0){
92955eca 234
ce7adfe9 235 AliVEvent* event = InputEvent();
236 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
92955eca 237 if (esd){
1e552991 238 if (!(fRunNumber == esd->GetRunNumber())) {
239 fRunNumber = esd->GetRunNumber();
92955eca 240 SetPhiDist();
ce7adfe9 241 }
92955eca 242
243
244 if (fUseMCRP) {
245 AliMCEvent* mcEvent = MCEvent();
246 if (mcEvent && mcEvent->GenEventHeader()) {
247 AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
248 if (headerH) fRP = headerH->ReactionPlaneAngle();
249 }
250 }
251
ce7adfe9 252 esdEP = esd->GetEventplane();
253 if (fSaveTrackContribution) {
254 esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
255 esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
e51055a0 256 esdEP->GetQContributionXArraysub1()->Set(esd->GetNumberOfTracks());
257 esdEP->GetQContributionYArraysub1()->Set(esd->GetNumberOfTracks());
258 esdEP->GetQContributionXArraysub2()->Set(esd->GetNumberOfTracks());
259 esdEP->GetQContributionYArraysub2()->Set(esd->GetNumberOfTracks());
ce7adfe9 260 }
261
1e552991 262 TObjArray* tracklist = new TObjArray;
263 if (fTrackType.CompareTo("GLOBAL")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
264 if (fTrackType.CompareTo("TPC")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
265 const int nt = tracklist->GetEntries();
ce7adfe9 266
71916547 267 if (nt>4){
1e552991 268 fQVector = new TVector2(GetQ(esdEP,tracklist));
ce7adfe9 269 fEventplaneQ = fQVector->Phi()/2;
e51055a0 270 GetQsub(qq1, qq2, tracklist, esdEP);
71916547 271 fQsub1 = new TVector2(qq1);
272 fQsub2 = new TVector2(qq2);
ce7adfe9 273 fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
90267ca6 274
ce7adfe9 275 esdEP->SetQVector(fQVector);
276 esdEP->SetEventplaneQ(fEventplaneQ);
277 esdEP->SetQsub(fQsub1,fQsub2);
278 esdEP->SetQsubRes(fQsubRes);
279
280 fHOutEventplaneQ->Fill(fEventplaneQ);
281 fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
71916547 282 fHOutNTEPRes->Fill(nt,fQsubRes);
ce7adfe9 283
284 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
285
71916547 286 for (int iter = 0; iter<nt;iter++){
1e552991 287 AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
b1a983b4 288 if (track) {
289 float delta = track->Phi()-fEventplaneQ;
290 while (delta < 0) delta += TMath::Pi();
291 while (delta > TMath::Pi()) delta -= TMath::Pi();
292 fHOutPTPsi->Fill(track->Pt(),delta);
293 fHOutPhi->Fill(track->Phi());
294 fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
295 }
ce7adfe9 296 }
297
298 AliESDtrack* trmax = esd->GetTrack(0);
71916547 299 for (int iter = 1; iter<nt;iter++){
1e552991 300 AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
b1a983b4 301 if (track && (track->Pt() > trmax->Pt())) trmax = track;
ce7adfe9 302 }
303 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
304 }
1e552991 305 delete tracklist;
306 tracklist = 0;
ce7adfe9 307 }
308 }
309
81e3aec0 310 else if (fAnalysisInput.CompareTo("AOD")==0){
5b53b816 311 AliVEvent* event = InputEvent();
312 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
e51055a0 313
dce05057 314 if (aod){
315 if (!(fRunNumber == aod->GetRunNumber())) {
81e3aec0 316 fRunNumber = aod->GetRunNumber();
317 SetPhiDist();
dce05057 318 }
81e3aec0 319 }
320
321 if (fUseMCRP) {
322 AliAODMCHeader *headerH = dynamic_cast<AliAODMCHeader*>(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
323 if (headerH) fRP = headerH->GetReactionPlaneAngle();
324 }
e51055a0 325
81e3aec0 326 if (aod){
e51055a0 327 esdEP = aod->GetHeader()->GetEventplaneP();
5b53b816 328 esdEP->Reset();
e51055a0 329
5b53b816 330 Int_t maxID = 0;
331 TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID);
e51055a0 332
81e3aec0 333 if (fSaveTrackContribution) {
334 esdEP->GetQContributionXArray()->Set(maxID+1);
335 esdEP->GetQContributionYArray()->Set(maxID+1);
336 esdEP->GetQContributionXArraysub1()->Set(maxID+1);
337 esdEP->GetQContributionYArraysub1()->Set(maxID+1);
338 esdEP->GetQContributionXArraysub2()->Set(maxID+1);
339 esdEP->GetQContributionYArraysub2()->Set(maxID+1);
340 }
e51055a0 341
81e3aec0 342 const int NT = tracklist->GetEntries();
e51055a0 343
81e3aec0 344 if (NT>4){
345 fQVector = new TVector2(GetQ(esdEP,tracklist));
346 fEventplaneQ = fQVector->Phi()/2.;
347 GetQsub(qq1, qq2, tracklist, esdEP);
348 fQsub1 = new TVector2(qq1);
349 fQsub2 = new TVector2(qq2);
350 fQsubRes = (fQsub1->Phi()/2. - fQsub2->Phi()/2.);
e51055a0 351
81e3aec0 352 esdEP->SetQVector(fQVector);
353 esdEP->SetEventplaneQ(fEventplaneQ);
354 esdEP->SetQsub(fQsub1,fQsub2);
355 esdEP->SetQsubRes(fQsubRes);
e51055a0 356
81e3aec0 357 fHOutEventplaneQ->Fill(fEventplaneQ);
358 fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.);
359 fHOutNTEPRes->Fill(NT,fQsubRes);
e51055a0 360
361 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
362
363 for (int iter = 0; iter<NT;iter++){
364 AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
365 if (track) {
366 float delta = track->Phi()-fEventplaneQ;
367 while (delta < 0) delta += TMath::Pi();
368 while (delta > TMath::Pi()) delta -= TMath::Pi();
369 fHOutPTPsi->Fill(track->Pt(),delta);
370 fHOutPhi->Fill(track->Phi());
371 fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
372 }
373 }
374
375 AliAODTrack* trmax = aod->GetTrack(0);
376 for (int iter = 1; iter<NT;iter++){
377 AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
378 if (track && (track->Pt() > trmax->Pt())) trmax = track;
379 }
380 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
381 }
382 delete tracklist;
383 tracklist = 0;
384 }
385
386
ce7adfe9 387 }
e51055a0 388
389
ce7adfe9 390 else {
391 printf(" Analysis Input not known!\n\n ");
392 return;
393 }
394 PostData(1, fOutputList);
395}
396
397//________________________________________________________________________
398void AliEPSelectionTask::Terminate(Option_t */*option*/)
399{
400 // Terminate analysis
401}
402
403//__________________________________________________________________________
404TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
405{
71916547 406// Get the Q vector
ce7adfe9 407 TVector2 mQ;
408 float mQx=0, mQy=0;
e51055a0 409 AliVTrack* track;
ce7adfe9 410 Double_t weight;
e51055a0 411 Int_t idtemp = -1;
ce7adfe9 412
71916547 413 int nt = tracklist->GetEntries();
ce7adfe9 414
71916547 415 for (int i=0; i<nt; i++){
ce7adfe9 416 weight = 1;
e51055a0 417 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
b1a983b4 418 if (track) {
419 weight = GetWeight(track);
e51055a0 420 if (fSaveTrackContribution){
421 idtemp = track->GetID();
422 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
423 EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),idtemp);
424 EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),idtemp);
425 }
426 mQx += (weight*cos(2*track->Phi()));
427 mQy += (weight*sin(2*track->Phi()));
ce7adfe9 428 }
ce7adfe9 429 }
430 mQ.Set(mQx,mQy);
431 return mQ;
432}
433
434 //________________________________________________________________________
e51055a0 435void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP)
ce7adfe9 436{
71916547 437// Get Qsub
ce7adfe9 438 TVector2 mQ[2];
439 float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
440 Double_t weight;
441
e51055a0 442 AliVTrack* track;
71916547 443 TRandom2 rn = 0;
ce7adfe9 444
71916547 445 int nt = tracklist->GetEntries();
ce7adfe9 446 int trackcounter1=0, trackcounter2=0;
e51055a0 447 int idtemp = 0;
448
449 if (fSplitMethod == AliEPSelectionTask::kRandom){
450
451 for (Int_t i = 0; i < nt; i++) {
452 weight = 1;
453 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
454 if (!track) continue;
455 weight = GetWeight(track);
456 idtemp = track->GetID();
457 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
ce7adfe9 458
e51055a0 459 // This loop splits the track set into 2 random subsets
460 if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
461 float random = rn.Rndm();
462 if(random < .5){
463 mQx1 += (weight*cos(2*track->Phi()));
464 mQy1 += (weight*sin(2*track->Phi()));
465 if (fSaveTrackContribution){
466 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
467 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
468 }
469 trackcounter1++;
470 }
471 else {
472 mQx2 += (weight*cos(2*track->Phi()));
473 mQy2 += (weight*sin(2*track->Phi()));
474 if (fSaveTrackContribution){
475 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
476 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
477 }
478 trackcounter2++;
479 }
480 }
481 else if( trackcounter1 >= int(nt/2.)){
482 mQx2 += (weight*cos(2*track->Phi()));
483 mQy2 += (weight*sin(2*track->Phi()));
484 if (fSaveTrackContribution){
485 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
486 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
487 }
488 trackcounter2++;
489 }
490 else {
ce7adfe9 491 mQx1 += (weight*cos(2*track->Phi()));
492 mQy1 += (weight*sin(2*track->Phi()));
e51055a0 493 if (fSaveTrackContribution){
494 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
495 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
496 }
ce7adfe9 497 trackcounter1++;
498 }
e51055a0 499 }
500 } else if (fSplitMethod == AliEPSelectionTask::kEta) {
501
502 for (Int_t i = 0; i < nt; i++) {
503 weight = 1;
504 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
505 if (!track) continue;
506 weight = GetWeight(track);
507 Double_t eta = track->Eta();
508 idtemp = track->GetID();
509 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
510
511 if (eta > fEtaGap/2.) {
512 mQx1 += (weight*cos(2*track->Phi()));
513 mQy1 += (weight*sin(2*track->Phi()));
514 if (fSaveTrackContribution){
515 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
516 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
517 }
518 } else if (eta < -1.*fEtaGap/2.) {
ce7adfe9 519 mQx2 += (weight*cos(2*track->Phi()));
520 mQy2 += (weight*sin(2*track->Phi()));
e51055a0 521 if (fSaveTrackContribution){
522 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
523 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
524 }
ce7adfe9 525 }
526 }
e51055a0 527 } else {
528 printf("plane resolution determination method not available!\n\n ");
529 return;
ce7adfe9 530 }
e51055a0 531
ce7adfe9 532 mQ[0].Set(mQx1,mQy1);
533 mQ[1].Set(mQx2,mQy2);
534 Q1 = mQ[0];
535 Q2 = mQ[1];
536}
537
538//________________________________________________________________________
90267ca6 539void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
ce7adfe9 540
1e552991 541 if(fESDtrackCuts){
542 delete fESDtrackCuts;
543 fESDtrackCuts = 0;
544 }
e51055a0 545 if (fAnalysisInput.CompareTo("AOD")==0){
546 AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");
547 fUsercuts = kFALSE;
548 SetTrackType("TPC");
549 return;
550 }
1e552991 551 fUsercuts = kTRUE;
90267ca6 552 fESDtrackCuts = trackcuts;
ce7adfe9 553}
554
e51055a0 555//________________________________________________________________________
556void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup){
557
558 if(fESDtrackCuts){
559 delete fESDtrackCuts;
560 fESDtrackCuts = 0;
561 }
562 if (fAnalysisInput.CompareTo("ESD")==0){
563 AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");
564 fUsercuts = kFALSE;
565 SetTrackType("TPC");
566 return;
567 }
568 fUsercuts = kTRUE;
569 fESDtrackCuts = new AliESDtrackCuts();
570 fESDtrackCuts->SetPtRange(ptlow,ptup);
571 fESDtrackCuts->SetEtaRange(etalow,etaup);
572 fAODfilterbit = filterbit;
573}
574
90267ca6 575//_____________________________________________________________________________
e51055a0 576
90267ca6 577void AliEPSelectionTask::SetTrackType(TString tracktype){
578 fTrackType = tracktype;
1e552991 579 if (!fUsercuts) {
e51055a0 580 if (fTrackType.CompareTo("GLOBAL")==0){
581 fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
582 fAODfilterbit = 32;
583 }
584 if (fTrackType.CompareTo("TPC")==0){
585 fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
586 fAODfilterbit = 128;
587 }
588 fESDtrackCuts->SetPtRange(0.15,20.);
90267ca6 589 fESDtrackCuts->SetEtaRange(-0.8,0.8);
ce7adfe9 590 }
ce7adfe9 591}
592
593//________________________________________________________________________
e51055a0 594Double_t AliEPSelectionTask::GetWeight(TObject* track1)
ce7adfe9 595{
596 Double_t ptweight=1;
e51055a0 597 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
5b53b816 598 if (fUsePtWeight && track) {
ce7adfe9 599 if (track->Pt()<2) ptweight=track->Pt();
600 else ptweight=2;
601 }
602 return ptweight*GetPhiWeight(track);
603}
604
605//________________________________________________________________________
e51055a0 606Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
ce7adfe9 607{
608 Double_t phiweight=1;
e51055a0 609 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
ce7adfe9 610
e4f1ed0c 611 if (fUsePhiWeight && fPhiDist && track) {
ce7adfe9 612 Double_t nParticles = fPhiDist->Integral();
613 Double_t nPhibins = fPhiDist->GetNbinsX();
614
e51055a0 615 Double_t Phi = track->Phi();
ce7adfe9 616
e51055a0 617 while (Phi<0) Phi += TMath::TwoPi();
618 while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
ce7adfe9 619
e51055a0 620 Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
ce7adfe9 621
e51055a0 622 if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
ce7adfe9 623 }
624 return phiweight;
625}
90267ca6 626
627//__________________________________________________________________________
628void AliEPSelectionTask::SetPhiDist()
629{
1e552991 630 if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
90267ca6 631
1e552991 632 fPhiDist = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");
633 if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
90267ca6 634
635 }
636 else {
637 AliInfo("Using Custom Phi Distribution");
638 }
639
640 Bool_t emptybins;
641
642 int iter = 0;
643 while (iter<3){
644 emptybins = kFALSE;
645
646 for (int i=1; i<fPhiDist->GetNbinsX(); i++){
647 if (!((fPhiDist->GetBinContent(i))>0)) {
648 emptybins = kTRUE;
649 }
650 }
651 if (emptybins) {
652 cout << "empty bins - rebinning!" << endl;
653 fPhiDist->Rebin();
654 iter++;
655 }
656 else iter = 3;
657 }
658
659 if (emptybins) {
660 AliError("After Maximum of rebinning still empty Phi-bins!!!");
661 }
662}
663
664//__________________________________________________________________________
71916547 665void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
90267ca6 666{
e51055a0 667
1e552991 668 fUserphidist = kTRUE;
90267ca6 669
670 TFile f(infilename);
671 TObject* list = f.Get(listname);
672 fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
673 if (!fPhiDist) AliFatal("Phi Distribution not found!!!");
674
675 f.Close();
676}
e51055a0 677
678
679//_________________________________________________________________________
680TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
681{
682 TObjArray *acctracks = new TObjArray();
683
684 AliAODTrack *tr = 0;
685 Int_t maxid1 = 0;
686 Int_t maxidtemp = -1;
687 Float_t ptlow = 0;
688 Float_t ptup = 0;
689 Float_t etalow = 0;
690 Float_t etaup = 0;
691 fESDtrackCuts->GetPtRange(ptlow,ptup);
692 fESDtrackCuts->GetEtaRange(etalow,etaup);
693
694 for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
695 tr = aod->GetTrack(i);
696 maxidtemp = tr->GetID();
697 if(maxidtemp < 0 && fAODfilterbit != 128) continue;
698 if(maxidtemp > -1 && fAODfilterbit == 128) continue;
699 if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
700 if (maxidtemp > maxid1) maxid1 = maxidtemp;
701 if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow){
702 acctracks->Add(tr);
703 }
704 }
705
706 maxid = maxid1;
707
708 return acctracks;
709
710}