]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ANALYSIS/AliEPSelectionTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Extension for AOD
[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
e51055a0 310 else if (fAnalysisInput.CompareTo("AOD")==0){
311 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
312
313 if (!(fRunNumber == aod->GetRunNumber())) {
314 fRunNumber = aod->GetRunNumber();
315 SetPhiDist();
316 }
317
318 if (fUseMCRP) {
319 AliAODMCHeader *headerH = dynamic_cast<AliAODMCHeader*>(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
320 if (headerH) fRP = headerH->GetReactionPlaneAngle();
321 }
322
323 if (aod){
324 esdEP = aod->GetHeader()->GetEventplaneP();
325 if(esdEP) {esdEP->Reset();} // reset eventplane if not NULL
326
327 Int_t maxID = 0;
328 TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID);
329
330 if (fSaveTrackContribution) {
331 esdEP->GetQContributionXArray()->Set(maxID+1);
332 esdEP->GetQContributionYArray()->Set(maxID+1);
333 esdEP->GetQContributionXArraysub1()->Set(maxID+1);
334 esdEP->GetQContributionYArraysub1()->Set(maxID+1);
335 esdEP->GetQContributionXArraysub2()->Set(maxID+1);
336 esdEP->GetQContributionYArraysub2()->Set(maxID+1);
337 }
338
339 const int NT = tracklist->GetEntries();
340
341 if (NT>4){
342 fQVector = new TVector2(GetQ(esdEP,tracklist));
343 fEventplaneQ = fQVector->Phi()/2.;
344 GetQsub(qq1, qq2, tracklist, esdEP);
345 fQsub1 = new TVector2(qq1);
346 fQsub2 = new TVector2(qq2);
347 fQsubRes = (fQsub1->Phi()/2. - fQsub2->Phi()/2.);
348
349 esdEP->SetQVector(fQVector);
350 esdEP->SetEventplaneQ(fEventplaneQ);
351 esdEP->SetQsub(fQsub1,fQsub2);
352 esdEP->SetQsubRes(fQsubRes);
353
354 fHOutEventplaneQ->Fill(fEventplaneQ);
355 fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.);
356 fHOutNTEPRes->Fill(NT,fQsubRes);
357
358 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
359
360 for (int iter = 0; iter<NT;iter++){
361 AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
362 if (track) {
363 float delta = track->Phi()-fEventplaneQ;
364 while (delta < 0) delta += TMath::Pi();
365 while (delta > TMath::Pi()) delta -= TMath::Pi();
366 fHOutPTPsi->Fill(track->Pt(),delta);
367 fHOutPhi->Fill(track->Phi());
368 fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
369 }
370 }
371
372 AliAODTrack* trmax = aod->GetTrack(0);
373 for (int iter = 1; iter<NT;iter++){
374 AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
375 if (track && (track->Pt() > trmax->Pt())) trmax = track;
376 }
377 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
378 }
379 delete tracklist;
380 tracklist = 0;
381 }
382
383
ce7adfe9 384 }
e51055a0 385
386
ce7adfe9 387 else {
388 printf(" Analysis Input not known!\n\n ");
389 return;
390 }
391 PostData(1, fOutputList);
392}
393
394//________________________________________________________________________
395void AliEPSelectionTask::Terminate(Option_t */*option*/)
396{
397 // Terminate analysis
398}
399
400//__________________________________________________________________________
401TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
402{
71916547 403// Get the Q vector
ce7adfe9 404 TVector2 mQ;
405 float mQx=0, mQy=0;
e51055a0 406 AliVTrack* track;
ce7adfe9 407 Double_t weight;
e51055a0 408 Int_t idtemp = -1;
ce7adfe9 409
71916547 410 int nt = tracklist->GetEntries();
ce7adfe9 411
71916547 412 for (int i=0; i<nt; i++){
ce7adfe9 413 weight = 1;
e51055a0 414 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
b1a983b4 415 if (track) {
416 weight = GetWeight(track);
e51055a0 417 if (fSaveTrackContribution){
418 idtemp = track->GetID();
419 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
420 EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),idtemp);
421 EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),idtemp);
422 }
423 mQx += (weight*cos(2*track->Phi()));
424 mQy += (weight*sin(2*track->Phi()));
ce7adfe9 425 }
ce7adfe9 426 }
427 mQ.Set(mQx,mQy);
428 return mQ;
429}
430
431 //________________________________________________________________________
e51055a0 432void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP)
ce7adfe9 433{
71916547 434// Get Qsub
ce7adfe9 435 TVector2 mQ[2];
436 float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
437 Double_t weight;
438
e51055a0 439 AliVTrack* track;
71916547 440 TRandom2 rn = 0;
ce7adfe9 441
71916547 442 int nt = tracklist->GetEntries();
ce7adfe9 443 int trackcounter1=0, trackcounter2=0;
e51055a0 444 int idtemp = 0;
445
446 if (fSplitMethod == AliEPSelectionTask::kRandom){
447
448 for (Int_t i = 0; i < nt; i++) {
449 weight = 1;
450 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
451 if (!track) continue;
452 weight = GetWeight(track);
453 idtemp = track->GetID();
454 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
ce7adfe9 455
e51055a0 456 // This loop splits the track set into 2 random subsets
457 if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
458 float random = rn.Rndm();
459 if(random < .5){
460 mQx1 += (weight*cos(2*track->Phi()));
461 mQy1 += (weight*sin(2*track->Phi()));
462 if (fSaveTrackContribution){
463 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
464 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
465 }
466 trackcounter1++;
467 }
468 else {
469 mQx2 += (weight*cos(2*track->Phi()));
470 mQy2 += (weight*sin(2*track->Phi()));
471 if (fSaveTrackContribution){
472 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
473 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
474 }
475 trackcounter2++;
476 }
477 }
478 else if( trackcounter1 >= int(nt/2.)){
479 mQx2 += (weight*cos(2*track->Phi()));
480 mQy2 += (weight*sin(2*track->Phi()));
481 if (fSaveTrackContribution){
482 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
483 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
484 }
485 trackcounter2++;
486 }
487 else {
ce7adfe9 488 mQx1 += (weight*cos(2*track->Phi()));
489 mQy1 += (weight*sin(2*track->Phi()));
e51055a0 490 if (fSaveTrackContribution){
491 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
492 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
493 }
ce7adfe9 494 trackcounter1++;
495 }
e51055a0 496 }
497 } else if (fSplitMethod == AliEPSelectionTask::kEta) {
498
499 for (Int_t i = 0; i < nt; i++) {
500 weight = 1;
501 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
502 if (!track) continue;
503 weight = GetWeight(track);
504 Double_t eta = track->Eta();
505 idtemp = track->GetID();
506 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
507
508 if (eta > fEtaGap/2.) {
509 mQx1 += (weight*cos(2*track->Phi()));
510 mQy1 += (weight*sin(2*track->Phi()));
511 if (fSaveTrackContribution){
512 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
513 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
514 }
515 } else if (eta < -1.*fEtaGap/2.) {
ce7adfe9 516 mQx2 += (weight*cos(2*track->Phi()));
517 mQy2 += (weight*sin(2*track->Phi()));
e51055a0 518 if (fSaveTrackContribution){
519 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
520 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
521 }
ce7adfe9 522 }
523 }
e51055a0 524 } else {
525 printf("plane resolution determination method not available!\n\n ");
526 return;
ce7adfe9 527 }
e51055a0 528
ce7adfe9 529 mQ[0].Set(mQx1,mQy1);
530 mQ[1].Set(mQx2,mQy2);
531 Q1 = mQ[0];
532 Q2 = mQ[1];
533}
534
535//________________________________________________________________________
90267ca6 536void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
ce7adfe9 537
1e552991 538 if(fESDtrackCuts){
539 delete fESDtrackCuts;
540 fESDtrackCuts = 0;
541 }
e51055a0 542 if (fAnalysisInput.CompareTo("AOD")==0){
543 AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");
544 fUsercuts = kFALSE;
545 SetTrackType("TPC");
546 return;
547 }
1e552991 548 fUsercuts = kTRUE;
90267ca6 549 fESDtrackCuts = trackcuts;
ce7adfe9 550}
551
e51055a0 552//________________________________________________________________________
553void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup){
554
555 if(fESDtrackCuts){
556 delete fESDtrackCuts;
557 fESDtrackCuts = 0;
558 }
559 if (fAnalysisInput.CompareTo("ESD")==0){
560 AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");
561 fUsercuts = kFALSE;
562 SetTrackType("TPC");
563 return;
564 }
565 fUsercuts = kTRUE;
566 fESDtrackCuts = new AliESDtrackCuts();
567 fESDtrackCuts->SetPtRange(ptlow,ptup);
568 fESDtrackCuts->SetEtaRange(etalow,etaup);
569 fAODfilterbit = filterbit;
570}
571
90267ca6 572//_____________________________________________________________________________
e51055a0 573
90267ca6 574void AliEPSelectionTask::SetTrackType(TString tracktype){
575 fTrackType = tracktype;
1e552991 576 if (!fUsercuts) {
e51055a0 577 if (fTrackType.CompareTo("GLOBAL")==0){
578 fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
579 fAODfilterbit = 32;
580 }
581 if (fTrackType.CompareTo("TPC")==0){
582 fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
583 fAODfilterbit = 128;
584 }
585 fESDtrackCuts->SetPtRange(0.15,20.);
90267ca6 586 fESDtrackCuts->SetEtaRange(-0.8,0.8);
ce7adfe9 587 }
ce7adfe9 588}
589
590//________________________________________________________________________
e51055a0 591Double_t AliEPSelectionTask::GetWeight(TObject* track1)
ce7adfe9 592{
593 Double_t ptweight=1;
e51055a0 594 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
ce7adfe9 595 if (fUsePtWeight) {
596 if (track->Pt()<2) ptweight=track->Pt();
597 else ptweight=2;
598 }
599 return ptweight*GetPhiWeight(track);
600}
601
602//________________________________________________________________________
e51055a0 603Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
ce7adfe9 604{
605 Double_t phiweight=1;
e51055a0 606 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
ce7adfe9 607
e4f1ed0c 608 if (fUsePhiWeight && fPhiDist && track) {
ce7adfe9 609 Double_t nParticles = fPhiDist->Integral();
610 Double_t nPhibins = fPhiDist->GetNbinsX();
611
e51055a0 612 Double_t Phi = track->Phi();
ce7adfe9 613
e51055a0 614 while (Phi<0) Phi += TMath::TwoPi();
615 while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
ce7adfe9 616
e51055a0 617 Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
ce7adfe9 618
e51055a0 619 if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
ce7adfe9 620 }
621 return phiweight;
622}
90267ca6 623
624//__________________________________________________________________________
625void AliEPSelectionTask::SetPhiDist()
626{
1e552991 627 if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
90267ca6 628
1e552991 629 fPhiDist = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");
630 if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
90267ca6 631
632 }
633 else {
634 AliInfo("Using Custom Phi Distribution");
635 }
636
637 Bool_t emptybins;
638
639 int iter = 0;
640 while (iter<3){
641 emptybins = kFALSE;
642
643 for (int i=1; i<fPhiDist->GetNbinsX(); i++){
644 if (!((fPhiDist->GetBinContent(i))>0)) {
645 emptybins = kTRUE;
646 }
647 }
648 if (emptybins) {
649 cout << "empty bins - rebinning!" << endl;
650 fPhiDist->Rebin();
651 iter++;
652 }
653 else iter = 3;
654 }
655
656 if (emptybins) {
657 AliError("After Maximum of rebinning still empty Phi-bins!!!");
658 }
659}
660
661//__________________________________________________________________________
71916547 662void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
90267ca6 663{
e51055a0 664
1e552991 665 fUserphidist = kTRUE;
90267ca6 666
667 TFile f(infilename);
668 TObject* list = f.Get(listname);
669 fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
670 if (!fPhiDist) AliFatal("Phi Distribution not found!!!");
671
672 f.Close();
673}
e51055a0 674
675
676//_________________________________________________________________________
677TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
678{
679 TObjArray *acctracks = new TObjArray();
680
681 AliAODTrack *tr = 0;
682 Int_t maxid1 = 0;
683 Int_t maxidtemp = -1;
684 Float_t ptlow = 0;
685 Float_t ptup = 0;
686 Float_t etalow = 0;
687 Float_t etaup = 0;
688 fESDtrackCuts->GetPtRange(ptlow,ptup);
689 fESDtrackCuts->GetEtaRange(etalow,etaup);
690
691 for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
692 tr = aod->GetTrack(i);
693 maxidtemp = tr->GetID();
694 if(maxidtemp < 0 && fAODfilterbit != 128) continue;
695 if(maxidtemp > -1 && fAODfilterbit == 128) continue;
696 if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
697 if (maxidtemp > maxid1) maxid1 = maxidtemp;
698 if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow){
699 acctracks->Add(tr);
700 }
701 }
702
703 maxid = maxid1;
704
705 return acctracks;
706
707}