]> 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
Reconstructor for digits and fix of some bugs
[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>
9663fa4b 28#include <THnSparse.h>
ce7adfe9 29#include <TProfile.h>
30#include <TFile.h>
31#include <TObjString.h>
32#include <TString.h>
33#include <TCanvas.h>
34#include <TROOT.h>
35#include <TDirectory.h>
36#include <TSystem.h>
37#include <iostream>
38#include <TRandom2.h>
39#include <TArrayF.h>
40
41#include "AliAnalysisManager.h"
42#include "AliVEvent.h"
43#include "AliESD.h"
44#include "AliESDEvent.h"
45#include "AliMCEvent.h"
46#include "AliESDtrack.h"
47#include "AliESDtrackCuts.h"
48#include "AliESDHeader.h"
49#include "AliESDInputHandler.h"
50#include "AliAODHandler.h"
51#include "AliAODEvent.h"
52#include "AliHeader.h"
53#include "AliGenHijingEventHeader.h"
54#include "AliAnalysisTaskSE.h"
55#include "AliPhysicsSelectionTask.h"
56#include "AliPhysicsSelection.h"
57#include "AliBackgroundSelection.h"
58#include "AliESDUtils.h"
90267ca6 59#include "AliOADBContainer.h"
e51055a0 60#include "AliAODMCHeader.h"
61#include "AliAODTrack.h"
62#include "AliVTrack.h"
ce7adfe9 63#include "AliEventplane.h"
64
c82bb898 65using std::cout;
66using std::endl;
ce7adfe9 67ClassImp(AliEPSelectionTask)
68
69//________________________________________________________________________
70AliEPSelectionTask::AliEPSelectionTask():
71AliAnalysisTaskSE(),
ce7adfe9 72 fAnalysisInput("ESD"),
90267ca6 73 fTrackType("TPC"),
9663fa4b 74 fPeriod(""),
6b7c66c9 75 fCentrality(-1.),
ce7adfe9 76 fUseMCRP(kFALSE),
77 fUsePhiWeight(kFALSE),
78 fUsePtWeight(kFALSE),
6b7c66c9 79 fUseRecentering(kFALSE),
ce7adfe9 80 fSaveTrackContribution(kFALSE),
1e552991 81 fUserphidist(kFALSE),
82 fUsercuts(kFALSE),
83 fRunNumber(-15),
e51055a0 84 fAODfilterbit(1),
85 fEtaGap(0.),
86 fSplitMethod(0),
ce7adfe9 87 fESDtrackCuts(0),
90267ca6 88 fEPContainer(0),
6b7c66c9 89 fQxContainer(0),
90 fQyContainer(0),
9663fa4b 91 fSparseDist(0),
92 fHruns(0),
ce7adfe9 93 fQVector(0),
94 fQContributionX(0),
95 fQContributionY(0),
96 fEventplaneQ(0),
97 fQsub1(0),
98 fQsub2(0),
99 fQsubRes(0),
100 fOutputList(0),
101 fHOutEventplaneQ(0),
102 fHOutPhi(0),
103 fHOutPhiCorr(0),
104 fHOutsub1sub2(0),
105 fHOutNTEPRes(0),
106 fHOutPTPsi(0),
107 fHOutDiff(0),
108 fHOutleadPTPsi(0)
109{
110 // Default constructor
111 AliInfo("Event Plane Selection enabled.");
9663fa4b 112 for(Int_t i = 0; i < 4; ++i) {
113 fPhiDist[i] = 0;
114 }
6b7c66c9 115 for(Int_t i = 0; i < 2; ++i) {
116 fQDist[i] = 0;
117 }
118}
ce7adfe9 119
120//________________________________________________________________________
121AliEPSelectionTask::AliEPSelectionTask(const char *name):
122 AliAnalysisTaskSE(name),
ce7adfe9 123 fAnalysisInput("ESD"),
90267ca6 124 fTrackType("TPC"),
9663fa4b 125 fPeriod(""),
6b7c66c9 126 fCentrality(-1.),
ce7adfe9 127 fUseMCRP(kFALSE),
128 fUsePhiWeight(kFALSE),
129 fUsePtWeight(kFALSE),
6b7c66c9 130 fUseRecentering(kFALSE),
ce7adfe9 131 fSaveTrackContribution(kFALSE),
1e552991 132 fUserphidist(kFALSE),
133 fUsercuts(kFALSE),
134 fRunNumber(-15),
e51055a0 135 fAODfilterbit(1),
136 fEtaGap(0.),
137 fSplitMethod(0),
ce7adfe9 138 fESDtrackCuts(0),
90267ca6 139 fEPContainer(0),
6b7c66c9 140 fQxContainer(0),
141 fQyContainer(0),
9663fa4b 142 fSparseDist(0),
143 fHruns(0),
ce7adfe9 144 fQVector(0),
145 fQContributionX(0),
146 fQContributionY(0),
147 fEventplaneQ(0),
148 fQsub1(0),
149 fQsub2(0),
150 fQsubRes(0),
151 fOutputList(0),
152 fHOutEventplaneQ(0),
153 fHOutPhi(0),
154 fHOutPhiCorr(0),
155 fHOutsub1sub2(0),
156 fHOutNTEPRes(0),
157 fHOutPTPsi(0),
158 fHOutDiff(0),
159 fHOutleadPTPsi(0)
160{
161 // Default constructor
162 AliInfo("Event Plane Selection enabled.");
163 DefineOutput(1, TList::Class());
9663fa4b 164 for(Int_t i = 0; i < 4; i++) {
165 fPhiDist[i] = 0;
166 }
6b7c66c9 167 for(Int_t i = 0; i < 2; ++i) {
168 fQDist[i] = 0;
169 }
ce7adfe9 170}
171
172//________________________________________________________________________
173AliEPSelectionTask::~AliEPSelectionTask()
174{
175 // Destructor
176 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
177 delete fOutputList;
178 fOutputList = 0;
179 }
180 if (fESDtrackCuts){
181 delete fESDtrackCuts;
182 fESDtrackCuts = 0;
183 }
1e552991 184 if (fUserphidist) {
9663fa4b 185 if (fPhiDist[0]) {
186 delete fPhiDist[0];
187 fPhiDist[0] = 0;
1e552991 188 }
90267ca6 189 }
190 if (fEPContainer){
191 delete fEPContainer;
192 fEPContainer = 0;
193 }
77c17385 194 if (fPeriod.CompareTo("LHC11h")==0){
9663fa4b 195 for(Int_t i = 0; i < 4; i++) {
196 if(fPhiDist[i]){
197 delete fPhiDist[i];
198 fPhiDist[i] = 0;
199 }
200 }
201 if(fHruns) delete fHruns;
202 }
6b7c66c9 203 if(fQDist[0] && fQDist[1]) {
204 for(Int_t i = 0; i < 2; i++) {
205 if(fQDist[i]){
206 delete fQDist[i];
207 fQDist[i] = 0;
208 }
209 }
210 }
ce7adfe9 211}
212
213//________________________________________________________________________
214void AliEPSelectionTask::UserCreateOutputObjects()
215{
216 // Create the output containers
217 if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
218 AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
219
220 fOutputList = new TList();
221 fOutputList->SetOwner();
222 fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
223 fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
224 fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
225 fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
226 fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
227 fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
228 fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
229 fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
230
231 fOutputList->Add(fHOutEventplaneQ);
232 fOutputList->Add(fHOutPhi);
233 fOutputList->Add(fHOutPhiCorr);
234 fOutputList->Add(fHOutsub1sub2);
235 fOutputList->Add(fHOutNTEPRes);
236 fOutputList->Add(fHOutPTPsi);
237 fOutputList->Add(fHOutleadPTPsi);
238 fOutputList->Add(fHOutDiff);
239
240 PostData(1, fOutputList);
90267ca6 241
e51055a0 242
ce7adfe9 243}
244
245//________________________________________________________________________
246void AliEPSelectionTask::UserExec(Option_t */*option*/)
247{
248 // Execute analysis for current event:
249 if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
90267ca6 250
1e552991 251// fRunNumber = -15;
ce7adfe9 252
e51055a0 253 AliEventplane *esdEP;
71916547 254 TVector2 qq1;
255 TVector2 qq2;
9663fa4b 256 Double_t fRP = 0.; // monte carlo reaction plane angle
ce7adfe9 257
258 if (fAnalysisInput.CompareTo("ESD")==0){
92955eca 259
ce7adfe9 260 AliVEvent* event = InputEvent();
261 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
92955eca 262 if (esd){
1e552991 263 if (!(fRunNumber == esd->GetRunNumber())) {
264 fRunNumber = esd->GetRunNumber();
9663fa4b 265 AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber));
266 SetOADBandPeriod();
267 SetPhiDist();
6b7c66c9 268 SetQvectorDist(); // recentring objects
ce7adfe9 269 }
92955eca 270
271
272 if (fUseMCRP) {
273 AliMCEvent* mcEvent = MCEvent();
274 if (mcEvent && mcEvent->GenEventHeader()) {
275 AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
276 if (headerH) fRP = headerH->ReactionPlaneAngle();
277 }
278 }
279
ce7adfe9 280 esdEP = esd->GetEventplane();
281 if (fSaveTrackContribution) {
282 esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
283 esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
e51055a0 284 esdEP->GetQContributionXArraysub1()->Set(esd->GetNumberOfTracks());
285 esdEP->GetQContributionYArraysub1()->Set(esd->GetNumberOfTracks());
286 esdEP->GetQContributionXArraysub2()->Set(esd->GetNumberOfTracks());
287 esdEP->GetQContributionYArraysub2()->Set(esd->GetNumberOfTracks());
ce7adfe9 288 }
289
1e552991 290 TObjArray* tracklist = new TObjArray;
291 if (fTrackType.CompareTo("GLOBAL")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
9663fa4b 292 if (fTrackType.CompareTo("TPC")==0 && fPeriod.CompareTo("LHC10h")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
293 else if (fTrackType.CompareTo("TPC")==0 && fPeriod.CompareTo("LHC11h")==0) tracklist = GetTracksForLHC11h(esd);
1e552991 294 const int nt = tracklist->GetEntries();
ce7adfe9 295
71916547 296 if (nt>4){
6b7c66c9 297
298 // qvector full event
1e552991 299 fQVector = new TVector2(GetQ(esdEP,tracklist));
ce7adfe9 300 fEventplaneQ = fQVector->Phi()/2;
6b7c66c9 301 // q vector subevents
e51055a0 302 GetQsub(qq1, qq2, tracklist, esdEP);
71916547 303 fQsub1 = new TVector2(qq1);
304 fQsub2 = new TVector2(qq2);
ce7adfe9 305 fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
90267ca6 306
ce7adfe9 307 esdEP->SetQVector(fQVector);
308 esdEP->SetEventplaneQ(fEventplaneQ);
309 esdEP->SetQsub(fQsub1,fQsub2);
310 esdEP->SetQsubRes(fQsubRes);
311
312 fHOutEventplaneQ->Fill(fEventplaneQ);
313 fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
71916547 314 fHOutNTEPRes->Fill(nt,fQsubRes);
ce7adfe9 315
316 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
317
71916547 318 for (int iter = 0; iter<nt;iter++){
1e552991 319 AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
b1a983b4 320 if (track) {
321 float delta = track->Phi()-fEventplaneQ;
322 while (delta < 0) delta += TMath::Pi();
323 while (delta > TMath::Pi()) delta -= TMath::Pi();
324 fHOutPTPsi->Fill(track->Pt(),delta);
325 fHOutPhi->Fill(track->Phi());
9663fa4b 326 fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
327 }
ce7adfe9 328 }
329
330 AliESDtrack* trmax = esd->GetTrack(0);
71916547 331 for (int iter = 1; iter<nt;iter++){
1e552991 332 AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
b1a983b4 333 if (track && (track->Pt() > trmax->Pt())) trmax = track;
ce7adfe9 334 }
335 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
336 }
9663fa4b 337 tracklist->Clear();
1e552991 338 delete tracklist;
339 tracklist = 0;
ce7adfe9 340 }
341 }
342
81e3aec0 343 else if (fAnalysisInput.CompareTo("AOD")==0){
5b53b816 344 AliVEvent* event = InputEvent();
345 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
e51055a0 346
dce05057 347 if (aod){
6b7c66c9 348
349 // get centrality of the event
350 AliAODHeader *header=aod->GetHeader();
351 AliCentrality *centrality=header->GetCentralityP();
352 if(!centrality) AliError(Form("No AliCentrality attached to AOD header"));
353 fCentrality = centrality->GetCentralityPercentile("V0M");
354
dce05057 355 if (!(fRunNumber == aod->GetRunNumber())) {
81e3aec0 356 fRunNumber = aod->GetRunNumber();
9663fa4b 357 AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber));
358 SetOADBandPeriod();
6b7c66c9 359 SetPhiDist();
360 SetQvectorDist();
dce05057 361 }
81e3aec0 362
e24b883a 363 if (fUseMCRP) {
364 AliAODMCHeader *headerH = dynamic_cast<AliAODMCHeader*>(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
365 if (headerH) fRP = headerH->GetReactionPlaneAngle();
366 }
e51055a0 367
e51055a0 368 esdEP = aod->GetHeader()->GetEventplaneP();
960edfd3 369 if(!esdEP) return; // protection against missing EP branch (nanoAODs)
5b53b816 370 esdEP->Reset();
e51055a0 371
5b53b816 372 Int_t maxID = 0;
373 TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID);
e51055a0 374
e24b883a 375 if (fSaveTrackContribution) {
376 esdEP->GetQContributionXArray()->Set(maxID+1);
377 esdEP->GetQContributionYArray()->Set(maxID+1);
378 esdEP->GetQContributionXArraysub1()->Set(maxID+1);
379 esdEP->GetQContributionYArraysub1()->Set(maxID+1);
380 esdEP->GetQContributionXArraysub2()->Set(maxID+1);
381 esdEP->GetQContributionYArraysub2()->Set(maxID+1);
382 }
e51055a0 383
e24b883a 384 const int NT = tracklist->GetEntries();
e51055a0 385
e24b883a 386 if (NT>4){
6b7c66c9 387
388 // qvector full event
e24b883a 389 fQVector = new TVector2(GetQ(esdEP,tracklist));
6b7c66c9 390 fEventplaneQ = fQVector->Phi()/2;
391 // q vector subevents
e24b883a 392 GetQsub(qq1, qq2, tracklist, esdEP);
393 fQsub1 = new TVector2(qq1);
394 fQsub2 = new TVector2(qq2);
6b7c66c9 395 fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
396
e24b883a 397 esdEP->SetQVector(fQVector);
398 esdEP->SetEventplaneQ(fEventplaneQ);
399 esdEP->SetQsub(fQsub1,fQsub2);
400 esdEP->SetQsubRes(fQsubRes);
401
402 fHOutEventplaneQ->Fill(fEventplaneQ);
403 fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.);
404 fHOutNTEPRes->Fill(NT,fQsubRes);
e51055a0 405
e51055a0 406 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
407
408 for (int iter = 0; iter<NT;iter++){
409 AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
410 if (track) {
411 float delta = track->Phi()-fEventplaneQ;
412 while (delta < 0) delta += TMath::Pi();
413 while (delta > TMath::Pi()) delta -= TMath::Pi();
414 fHOutPTPsi->Fill(track->Pt(),delta);
415 fHOutPhi->Fill(track->Phi());
416 fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
417 }
418 }
419
420 AliAODTrack* trmax = aod->GetTrack(0);
421 for (int iter = 1; iter<NT;iter++){
422 AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
423 if (track && (track->Pt() > trmax->Pt())) trmax = track;
424 }
425 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
426 }
427 delete tracklist;
428 tracklist = 0;
429 }
430
431
ce7adfe9 432 }
e51055a0 433
434
ce7adfe9 435 else {
436 printf(" Analysis Input not known!\n\n ");
437 return;
438 }
439 PostData(1, fOutputList);
440}
441
442//________________________________________________________________________
443void AliEPSelectionTask::Terminate(Option_t */*option*/)
444{
445 // Terminate analysis
446}
447
448//__________________________________________________________________________
449TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
450{
6b7c66c9 451 // Get the Q vector
ce7adfe9 452 TVector2 mQ;
453 float mQx=0, mQy=0;
e51055a0 454 AliVTrack* track;
ce7adfe9 455 Double_t weight;
e51055a0 456 Int_t idtemp = -1;
6b7c66c9 457 // get recentering values
458 Double_t mean[2], rms[2];
459 Recenter(0, mean);
460 Recenter(1, rms);
461
71916547 462 int nt = tracklist->GetEntries();
ce7adfe9 463
71916547 464 for (int i=0; i<nt; i++){
ce7adfe9 465 weight = 1;
e51055a0 466 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
b1a983b4 467 if (track) {
468 weight = GetWeight(track);
e51055a0 469 if (fSaveTrackContribution){
470 idtemp = track->GetID();
471 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
6b7c66c9 472 EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
473 EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
e51055a0 474 }
6b7c66c9 475 mQx += (weight*cos(2*track->Phi())/rms[0]);
476 mQy += (weight*sin(2*track->Phi())/rms[1]);
ce7adfe9 477 }
ce7adfe9 478 }
6b7c66c9 479 mQ.Set(mQx-(mean[0]/rms[0]), mQy-(mean[1]/rms[1]));
ce7adfe9 480 return mQ;
481}
482
483 //________________________________________________________________________
e51055a0 484void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP)
ce7adfe9 485{
6b7c66c9 486 // Get Qsub
ce7adfe9 487 TVector2 mQ[2];
488 float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
489 Double_t weight;
6b7c66c9 490 // get recentering values
491 Double_t mean[2], rms[2];
492 Recenter(0, mean);
493 Recenter(1, rms);
ce7adfe9 494
e51055a0 495 AliVTrack* track;
71916547 496 TRandom2 rn = 0;
ce7adfe9 497
71916547 498 int nt = tracklist->GetEntries();
ce7adfe9 499 int trackcounter1=0, trackcounter2=0;
e51055a0 500 int idtemp = 0;
501
502 if (fSplitMethod == AliEPSelectionTask::kRandom){
503
504 for (Int_t i = 0; i < nt; i++) {
505 weight = 1;
506 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
507 if (!track) continue;
508 weight = GetWeight(track);
509 idtemp = track->GetID();
510 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
ce7adfe9 511
e51055a0 512 // This loop splits the track set into 2 random subsets
513 if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
514 float random = rn.Rndm();
515 if(random < .5){
6b7c66c9 516 mQx1 += (weight*cos(2*track->Phi())/rms[0]);
517 mQy1 += (weight*sin(2*track->Phi())/rms[1]);
e51055a0 518 if (fSaveTrackContribution){
6b7c66c9 519 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
520 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
e51055a0 521 }
522 trackcounter1++;
523 }
524 else {
6b7c66c9 525 mQx2 += (weight*cos(2*track->Phi())/rms[0]);
526 mQy2 += (weight*sin(2*track->Phi())/rms[1]);
e51055a0 527 if (fSaveTrackContribution){
6b7c66c9 528 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
529 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
e51055a0 530 }
531 trackcounter2++;
532 }
533 }
534 else if( trackcounter1 >= int(nt/2.)){
6b7c66c9 535 mQx2 += (weight*cos(2*track->Phi())/rms[0]);
536 mQy2 += (weight*sin(2*track->Phi())/rms[1]);
e51055a0 537 if (fSaveTrackContribution){
6b7c66c9 538 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
539 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
e51055a0 540 }
541 trackcounter2++;
542 }
543 else {
6b7c66c9 544 mQx1 += (weight*cos(2*track->Phi())/rms[0]);
545 mQy1 += (weight*sin(2*track->Phi())/rms[1]);
e51055a0 546 if (fSaveTrackContribution){
6b7c66c9 547 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
548 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
e51055a0 549 }
ce7adfe9 550 trackcounter1++;
551 }
e51055a0 552 }
553 } else if (fSplitMethod == AliEPSelectionTask::kEta) {
554
555 for (Int_t i = 0; i < nt; i++) {
556 weight = 1;
557 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
558 if (!track) continue;
559 weight = GetWeight(track);
560 Double_t eta = track->Eta();
561 idtemp = track->GetID();
562 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
563
564 if (eta > fEtaGap/2.) {
6b7c66c9 565 mQx1 += (weight*cos(2*track->Phi())/rms[0]);
566 mQy1 += (weight*sin(2*track->Phi())/rms[1]);
e51055a0 567 if (fSaveTrackContribution){
6b7c66c9 568 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
569 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
e51055a0 570 }
571 } else if (eta < -1.*fEtaGap/2.) {
6b7c66c9 572 mQx2 += (weight*cos(2*track->Phi())/rms[0]);
573 mQy2 += (weight*sin(2*track->Phi())/rms[1]);
e51055a0 574 if (fSaveTrackContribution){
6b7c66c9 575 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
576 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
c39c35d1 577 }
578 }
579 }
580 } else if (fSplitMethod == AliEPSelectionTask::kCharge) {
581
582 for (Int_t i = 0; i < nt; i++) {
583 weight = 1;
584 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
585 if (!track) continue;
586 weight = GetWeight(track);
587 Short_t cha = track->Charge();
588 idtemp = track->GetID();
589 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
590
591 if (cha > 0) {
6b7c66c9 592 mQx1 += (weight*cos(2*track->Phi())/rms[0]);
593 mQy1 += (weight*sin(2*track->Phi())/rms[1]);
c39c35d1 594 if (fSaveTrackContribution){
6b7c66c9 595 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
596 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
c39c35d1 597 }
598 } else if (cha < 0) {
6b7c66c9 599 mQx2 += (weight*cos(2*track->Phi())/rms[0]);
600 mQy2 += (weight*sin(2*track->Phi())/rms[1]);
c39c35d1 601 if (fSaveTrackContribution){
6b7c66c9 602 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
603 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
e51055a0 604 }
ce7adfe9 605 }
606 }
e51055a0 607 } else {
608 printf("plane resolution determination method not available!\n\n ");
609 return;
ce7adfe9 610 }
6b7c66c9 611 // apply recenetering
612 mQ[0].Set(mQx1-(mean[0]/rms[0]), mQy1-(mean[1]/rms[1]));
613 mQ[1].Set(mQx2-(mean[0]/rms[0]), mQy2-(mean[1]/rms[1]));
ce7adfe9 614 Q1 = mQ[0];
615 Q2 = mQ[1];
616}
617
618//________________________________________________________________________
90267ca6 619void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
ce7adfe9 620
1e552991 621 if(fESDtrackCuts){
622 delete fESDtrackCuts;
623 fESDtrackCuts = 0;
624 }
e51055a0 625 if (fAnalysisInput.CompareTo("AOD")==0){
626 AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");
627 fUsercuts = kFALSE;
628 SetTrackType("TPC");
629 return;
630 }
1e552991 631 fUsercuts = kTRUE;
90267ca6 632 fESDtrackCuts = trackcuts;
ce7adfe9 633}
634
e51055a0 635//________________________________________________________________________
e1ffd90d 636void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup, Int_t ntpc){
e51055a0 637
638 if(fESDtrackCuts){
639 delete fESDtrackCuts;
640 fESDtrackCuts = 0;
641 }
642 if (fAnalysisInput.CompareTo("ESD")==0){
643 AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");
644 fUsercuts = kFALSE;
645 SetTrackType("TPC");
646 return;
647 }
648 fUsercuts = kTRUE;
649 fESDtrackCuts = new AliESDtrackCuts();
650 fESDtrackCuts->SetPtRange(ptlow,ptup);
e1ffd90d 651 fESDtrackCuts->SetMinNClustersTPC(ntpc);
e51055a0 652 fESDtrackCuts->SetEtaRange(etalow,etaup);
653 fAODfilterbit = filterbit;
654}
655
90267ca6 656//_____________________________________________________________________________
e51055a0 657
90267ca6 658void AliEPSelectionTask::SetTrackType(TString tracktype){
659 fTrackType = tracktype;
1e552991 660 if (!fUsercuts) {
e51055a0 661 if (fTrackType.CompareTo("GLOBAL")==0){
662 fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
663 fAODfilterbit = 32;
664 }
665 if (fTrackType.CompareTo("TPC")==0){
666 fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
667 fAODfilterbit = 128;
668 }
669 fESDtrackCuts->SetPtRange(0.15,20.);
90267ca6 670 fESDtrackCuts->SetEtaRange(-0.8,0.8);
ce7adfe9 671 }
ce7adfe9 672}
673
674//________________________________________________________________________
e51055a0 675Double_t AliEPSelectionTask::GetWeight(TObject* track1)
ce7adfe9 676{
677 Double_t ptweight=1;
e51055a0 678 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
5b53b816 679 if (fUsePtWeight && track) {
ce7adfe9 680 if (track->Pt()<2) ptweight=track->Pt();
681 else ptweight=2;
682 }
683 return ptweight*GetPhiWeight(track);
684}
685
686//________________________________________________________________________
e51055a0 687Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
ce7adfe9 688{
689 Double_t phiweight=1;
e51055a0 690 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
9663fa4b 691
77c17385 692 TH1F *phiDist = 0x0;
693 if(track) phiDist = SelectPhiDist(track);
ce7adfe9 694
9663fa4b 695 if (fUsePhiWeight && phiDist && track) {
696 Double_t nParticles = phiDist->Integral();
697 Double_t nPhibins = phiDist->GetNbinsX();
ce7adfe9 698
e51055a0 699 Double_t Phi = track->Phi();
ce7adfe9 700
e51055a0 701 while (Phi<0) Phi += TMath::TwoPi();
702 while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
ce7adfe9 703
9663fa4b 704 Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
ce7adfe9 705
e51055a0 706 if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
ce7adfe9 707 }
708 return phiweight;
709}
90267ca6 710
6b7c66c9 711//________________________________________________________________________
712void AliEPSelectionTask::Recenter(Int_t var, Double_t * values)
713{
714
715 if (fUseRecentering && fQDist[0] && fQDist[1] && fCentrality!=-1.) {
716 Int_t centbin = fQDist[0]->FindBin(fCentrality);
717
718 if(var==0) { // fill mean
719 values[0] = fQDist[0]->GetBinContent(centbin);
720 values[1] = fQDist[1]->GetBinContent(centbin);
721 }
722 else if(var==1) { // fill rms
723 values[0] = fQDist[0]->GetBinError(centbin);
724 values[1] = fQDist[1]->GetBinError(centbin);
725 // protection against division by zero
726 if(values[0]==0.0) values[0]=1.0;
727 if(values[1]==0.0) values[1]=1.0;
728 }
729 }
730 else { //default (no recentering)
731 values[0] = var;
732 values[1] = var;
733 }
734 return;
735}
736
90267ca6 737//__________________________________________________________________________
738void AliEPSelectionTask::SetPhiDist()
739{
9593a21e 740 if(!fUserphidist && (fPeriod.CompareTo("LHC10h") == 0 || fPeriod.CompareTo("LHC11h") == 0)) { // if it's already set and custom class is required, we use the one provided by the user
90267ca6 741
9663fa4b 742 if (fPeriod.CompareTo("LHC10h")==0)
743 {
744 fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");}
745 else if(fPeriod.CompareTo("LHC11h")==0){
746 Int_t runbin=fHruns->FindBin(fRunNumber);
747 if (fHruns->GetBinContent(runbin) > 1){
748 fSparseDist->GetAxis(0)->SetRange(runbin,runbin);
749 }
750 else if(fHruns->GetBinContent(runbin) < 2){
751 fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights
752 AliInfo("Using integrated Phi-weights for this run");
753 }
754 for (Int_t i = 0; i<4 ;i++)
755 {
756 if(fPhiDist[i]){
757 delete fPhiDist[i];
758 fPhiDist[i] = 0x0;
759 }
760 if(i == 0){
761 fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge
762 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
763 if(i == 1){
764 fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge
765 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
766 if(i == 2){
767 fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge
768 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
769 if(i == 3){
770 fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge
771 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
772 fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi
773 fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber));
774 fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes
775 fSparseDist->GetAxis(2)->SetRange(1,2);
776 }
777 fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis
778 }
779
780 if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
90267ca6 781
782 }
9593a21e 783
90267ca6 784
9663fa4b 785 if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist){
786 Bool_t emptybins;
90267ca6 787
9663fa4b 788 int iter = 0;
789 while (iter<3){
790 emptybins = kFALSE;
90267ca6 791
9663fa4b 792 for (int i=1; i<fPhiDist[0]->GetNbinsX(); i++){
793 if (!((fPhiDist[0]->GetBinContent(i))>0)) {
794 emptybins = kTRUE;
795 }
796 }
797 if (emptybins) {
798 cout << "empty bins - rebinning!" << endl;
799 fPhiDist[0]->Rebin();
800 iter++;
801 }
802 else iter = 3;
803 }
90267ca6 804
9663fa4b 805 if (emptybins) {
806 AliError("After Maximum of rebinning still empty Phi-bins!!!");
807 }
90267ca6 808 }
9593a21e 809 if (fPeriod.CompareTo("LHC10h") != 0 && fPeriod.CompareTo("LHC11h") != 0 && !fUserphidist){
810 AliInfo("No Phi-weights available. All Phi weights set to 1");
811 SetUsePhiWeight(kFALSE);
812 }
90267ca6 813}
814
6b7c66c9 815//__________________________________________________________________________
816void AliEPSelectionTask::SetQvectorDist()
817{
818 if(!fUseRecentering) return;
819 AliInfo(Form("Setting q vector distributions"));
820 fQDist[0] = (TProfile*) fQxContainer->GetObject(fRunNumber, "Default");
821 fQDist[1] = (TProfile*) fQyContainer->GetObject(fRunNumber, "Default");
822
823 if (!fQDist[0] || !fQDist[1]) {
824 AliError(Form("Cannot find OADB q-vector distributions for run %d. Using default values (mean=0,rms=1).", fRunNumber));
825 return;
826 }
827
828 Bool_t emptybins;
829
830 int iter = 0;
831 while (iter<3){
832 emptybins = kFALSE;
833
834 for (int i=1; i<=fQDist[0]->GetNbinsX()-2; i++){
835 if (!((fQDist[0]->GetBinEntries(i))>0)) {
836 emptybins = kTRUE;
837 }
838 }
839 if (emptybins) {
840 cout << "empty bins - rebinning!" << endl;
841 fQDist[0]->Rebin();
842 fQDist[1]->Rebin();
843 iter++;
844 }
845 else iter = 3;
846 }
847
848 if (emptybins) {
849 AliError("After Maximum of rebinning still empty Qxy-bins!!!");
850 }
851}
852
90267ca6 853//__________________________________________________________________________
71916547 854void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
90267ca6 855{
e51055a0 856
1e552991 857 fUserphidist = kTRUE;
90267ca6 858
859 TFile f(infilename);
860 TObject* list = f.Get(listname);
9663fa4b 861 fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi");
862 if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!");
90267ca6 863
864 f.Close();
865}
e51055a0 866
e51055a0 867//_________________________________________________________________________
868TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
869{
870 TObjArray *acctracks = new TObjArray();
871
872 AliAODTrack *tr = 0;
873 Int_t maxid1 = 0;
874 Int_t maxidtemp = -1;
875 Float_t ptlow = 0;
876 Float_t ptup = 0;
877 Float_t etalow = 0;
878 Float_t etaup = 0;
879 fESDtrackCuts->GetPtRange(ptlow,ptup);
880 fESDtrackCuts->GetEtaRange(etalow,etaup);
e1ffd90d 881 Int_t ntpc = fESDtrackCuts->GetMinNClusterTPC();
e51055a0 882
883 for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
884 tr = aod->GetTrack(i);
885 maxidtemp = tr->GetID();
886 if(maxidtemp < 0 && fAODfilterbit != 128) continue;
887 if(maxidtemp > -1 && fAODfilterbit == 128) continue;
888 if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
889 if (maxidtemp > maxid1) maxid1 = maxidtemp;
e1ffd90d 890 if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow && tr->GetTPCNcls() > ntpc){
e51055a0 891 acctracks->Add(tr);
892 }
893 }
894
895 maxid = maxid1;
896
897 return acctracks;
898
899}
9663fa4b 900
901//_________________________________________________________________________
902void AliEPSelectionTask::SetOADBandPeriod()
903{
904 TString oadbfilename;
905
906 if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h
907 {fPeriod = "LHC10h";
908 if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
909
910 if (fAnalysisInput.CompareTo("AOD")==0){
911 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
912 } else if (fAnalysisInput.CompareTo("ESD")==0){
913 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
914 }
915
916 TFile foadb(oadbfilename);
917 if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
918
919 AliInfo("Using Standard OADB");
920 fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
921 if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
922 foadb.Close();
923 }
924 }
925
926 if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h
927 {fPeriod = "LHC11h";
928 if (!fUserphidist) {
929 // if it's already set and custom class is required, we use the one provided by the user
930
931 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
932 TFile *foadb = TFile::Open(oadbfilename);
933 if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
934
935 AliInfo("Using Standard OADB");
936 fSparseDist = (THnSparse*) foadb->Get("Default");
937 if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
938 foadb->Close();
7cdafa1e 939 if(!fHruns){
940 fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
941 fHruns->SetName("runsHisto");
942 }
6b7c66c9 943 }
944
945 if(fUseRecentering) {
946 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/eprecentering.root", AliAnalysisManager::GetOADBPath()));
947 TFile *foadb = TFile::Open(oadbfilename);
948 if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
949
950 AliInfo("Using Standard OADB");
951 fQxContainer = (AliOADBContainer*) foadb->Get("eprecentering.Qx");
952 fQyContainer = (AliOADBContainer*) foadb->Get("eprecentering.Qy");
953 if (!fQxContainer || !fQyContainer) AliFatal("Cannot fetch OADB container for EP recentering");
954 foadb->Close();
955 }
956
9663fa4b 957 }
958}
959
960//_________________________________________________________________________
961TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track)
962{
963 if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist) return fPhiDist[0];
964 else if(fPeriod.CompareTo("LHC11h")==0)
965 {
966 if (track->Charge() < 0)
967 {
968 if(track->Eta() < 0.) return fPhiDist[0];
969 else if (track->Eta() > 0.) return fPhiDist[2];
970 }
971 else if (track->Charge() > 0)
972 {
973 if(track->Eta() < 0.) return fPhiDist[1];
974 else if (track->Eta() > 0.) return fPhiDist[3];
975 }
976
977 }
978 return 0;
979}
980
981TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd)
982{
983 // Need to do this hack beacuse only this type of TPC only tracks in AOD is available and one type of Phi-weights is used
984 TObjArray *acctracks = new TObjArray();
985 acctracks->SetOwner(kTRUE);
986
987 const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
988
989 Float_t ptlow = 0;
990 Float_t ptup = 0;
991 Float_t etalow = 0;
992 Float_t etaup = 0;
993 fESDtrackCuts->GetPtRange(ptlow,ptup);
994 fESDtrackCuts->GetEtaRange(etalow,etaup);
995
996 Double_t pDCA[3] = { 0. }; // momentum at DCA
997 Double_t rDCA[3] = { 0. }; // position at DCA
998 Float_t dDCA[2] = {0.}; // DCA to the vertex d and z
999 Float_t cDCA[3] = {0.}; // covariance of impact parameters
1000
1001
1002
1003 for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){
1004
1005 AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise need to work with a copy
1006 //
1007 // Track selection
1008 if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue;
1009
1010 // create a tpc only tracl
1011 AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
1012 if(!track) continue;
1013
1014 if(track->Pt()>0.)
1015 {
1016 // only constrain tracks above threshold
1017 AliExternalTrackParam exParam;
1018 // take the B-field from the ESD, no 3D fieldMap available at this point
1019 Bool_t relate = false;
1020 relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
1021 if(!relate){
1022 delete track;
1023 continue;
1024 }
1025 // fetch the track parameters at the DCA (unconstraint)
1026 if(track->GetTPCInnerParam()){
1027 track->GetTPCInnerParam()->GetPxPyPz(pDCA);
1028 track->GetTPCInnerParam()->GetXYZ(rDCA);
1029 }
1030 // get the DCA to the vertex:
1031 track->GetImpactParametersTPC(dDCA,cDCA);
1032 // set the constrained parameters to the track
1033 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1034 }
1035
1036
1037 Float_t eta = track->Eta();
1038 Float_t pT = track->Pt();
1039
1040 if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){
1041 delete track;
1042 continue;
1043 }
1044
1045 acctracks->Add(track);
1046 }
1047
1048
1049 return acctracks;
1050
1051}
1052