]> 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
Calculation Ecross resored
[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
319e6239 350 AliAODHeader *header=dynamic_cast<AliAODHeader*>(aod->GetHeader());
351 if(!header) AliFatal("Not a standard AOD");
6b7c66c9 352 AliCentrality *centrality=header->GetCentralityP();
353 if(!centrality) AliError(Form("No AliCentrality attached to AOD header"));
354 fCentrality = centrality->GetCentralityPercentile("V0M");
355
dce05057 356 if (!(fRunNumber == aod->GetRunNumber())) {
81e3aec0 357 fRunNumber = aod->GetRunNumber();
9663fa4b 358 AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber));
359 SetOADBandPeriod();
6b7c66c9 360 SetPhiDist();
361 SetQvectorDist();
dce05057 362 }
81e3aec0 363
e24b883a 364 if (fUseMCRP) {
365 AliAODMCHeader *headerH = dynamic_cast<AliAODMCHeader*>(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
366 if (headerH) fRP = headerH->GetReactionPlaneAngle();
367 }
e51055a0 368
0a918d8d 369 esdEP = header->GetEventplaneP();
960edfd3 370 if(!esdEP) return; // protection against missing EP branch (nanoAODs)
5b53b816 371 esdEP->Reset();
e51055a0 372
5b53b816 373 Int_t maxID = 0;
374 TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID);
e51055a0 375
e24b883a 376 if (fSaveTrackContribution) {
377 esdEP->GetQContributionXArray()->Set(maxID+1);
378 esdEP->GetQContributionYArray()->Set(maxID+1);
379 esdEP->GetQContributionXArraysub1()->Set(maxID+1);
380 esdEP->GetQContributionYArraysub1()->Set(maxID+1);
381 esdEP->GetQContributionXArraysub2()->Set(maxID+1);
382 esdEP->GetQContributionYArraysub2()->Set(maxID+1);
383 }
e51055a0 384
e24b883a 385 const int NT = tracklist->GetEntries();
e51055a0 386
e24b883a 387 if (NT>4){
6b7c66c9 388
389 // qvector full event
e24b883a 390 fQVector = new TVector2(GetQ(esdEP,tracklist));
6b7c66c9 391 fEventplaneQ = fQVector->Phi()/2;
392 // q vector subevents
e24b883a 393 GetQsub(qq1, qq2, tracklist, esdEP);
394 fQsub1 = new TVector2(qq1);
395 fQsub2 = new TVector2(qq2);
6b7c66c9 396 fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
397
e24b883a 398 esdEP->SetQVector(fQVector);
399 esdEP->SetEventplaneQ(fEventplaneQ);
400 esdEP->SetQsub(fQsub1,fQsub2);
401 esdEP->SetQsubRes(fQsubRes);
402
403 fHOutEventplaneQ->Fill(fEventplaneQ);
404 fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.);
405 fHOutNTEPRes->Fill(NT,fQsubRes);
e51055a0 406
e51055a0 407 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
408
409 for (int iter = 0; iter<NT;iter++){
410 AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
411 if (track) {
412 float delta = track->Phi()-fEventplaneQ;
413 while (delta < 0) delta += TMath::Pi();
414 while (delta > TMath::Pi()) delta -= TMath::Pi();
415 fHOutPTPsi->Fill(track->Pt(),delta);
416 fHOutPhi->Fill(track->Phi());
417 fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
418 }
419 }
420
f15c1f69 421 AliAODTrack* trmax = dynamic_cast<AliAODTrack*>(aod->GetTrack(0));
422 if(!trmax) AliFatal("Not a standard AOD");
e51055a0 423 for (int iter = 1; iter<NT;iter++){
424 AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
425 if (track && (track->Pt() > trmax->Pt())) trmax = track;
426 }
427 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
428 }
429 delete tracklist;
430 tracklist = 0;
431 }
432
433
ce7adfe9 434 }
e51055a0 435
436
ce7adfe9 437 else {
438 printf(" Analysis Input not known!\n\n ");
439 return;
440 }
441 PostData(1, fOutputList);
442}
443
444//________________________________________________________________________
445void AliEPSelectionTask::Terminate(Option_t */*option*/)
446{
447 // Terminate analysis
448}
449
450//__________________________________________________________________________
451TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
452{
6b7c66c9 453 // Get the Q vector
ce7adfe9 454 TVector2 mQ;
455 float mQx=0, mQy=0;
e51055a0 456 AliVTrack* track;
ce7adfe9 457 Double_t weight;
e51055a0 458 Int_t idtemp = -1;
6b7c66c9 459 // get recentering values
460 Double_t mean[2], rms[2];
461 Recenter(0, mean);
462 Recenter(1, rms);
463
71916547 464 int nt = tracklist->GetEntries();
ce7adfe9 465
71916547 466 for (int i=0; i<nt; i++){
ce7adfe9 467 weight = 1;
e51055a0 468 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
b1a983b4 469 if (track) {
470 weight = GetWeight(track);
e51055a0 471 if (fSaveTrackContribution){
472 idtemp = track->GetID();
473 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
6b7c66c9 474 EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
475 EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
e51055a0 476 }
6b7c66c9 477 mQx += (weight*cos(2*track->Phi())/rms[0]);
478 mQy += (weight*sin(2*track->Phi())/rms[1]);
ce7adfe9 479 }
ce7adfe9 480 }
6b7c66c9 481 mQ.Set(mQx-(mean[0]/rms[0]), mQy-(mean[1]/rms[1]));
ce7adfe9 482 return mQ;
483}
484
485 //________________________________________________________________________
e51055a0 486void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP)
ce7adfe9 487{
6b7c66c9 488 // Get Qsub
ce7adfe9 489 TVector2 mQ[2];
490 float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
491 Double_t weight;
6b7c66c9 492 // get recentering values
493 Double_t mean[2], rms[2];
494 Recenter(0, mean);
495 Recenter(1, rms);
ce7adfe9 496
e51055a0 497 AliVTrack* track;
71916547 498 TRandom2 rn = 0;
ce7adfe9 499
71916547 500 int nt = tracklist->GetEntries();
ce7adfe9 501 int trackcounter1=0, trackcounter2=0;
e51055a0 502 int idtemp = 0;
503
504 if (fSplitMethod == AliEPSelectionTask::kRandom){
505
506 for (Int_t i = 0; i < nt; i++) {
507 weight = 1;
508 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
509 if (!track) continue;
510 weight = GetWeight(track);
511 idtemp = track->GetID();
512 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
ce7adfe9 513
e51055a0 514 // This loop splits the track set into 2 random subsets
515 if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
516 float random = rn.Rndm();
517 if(random < .5){
6b7c66c9 518 mQx1 += (weight*cos(2*track->Phi())/rms[0]);
519 mQy1 += (weight*sin(2*track->Phi())/rms[1]);
e51055a0 520 if (fSaveTrackContribution){
6b7c66c9 521 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
522 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
e51055a0 523 }
524 trackcounter1++;
525 }
526 else {
6b7c66c9 527 mQx2 += (weight*cos(2*track->Phi())/rms[0]);
528 mQy2 += (weight*sin(2*track->Phi())/rms[1]);
e51055a0 529 if (fSaveTrackContribution){
6b7c66c9 530 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
531 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
e51055a0 532 }
533 trackcounter2++;
534 }
535 }
536 else if( trackcounter1 >= int(nt/2.)){
6b7c66c9 537 mQx2 += (weight*cos(2*track->Phi())/rms[0]);
538 mQy2 += (weight*sin(2*track->Phi())/rms[1]);
e51055a0 539 if (fSaveTrackContribution){
6b7c66c9 540 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
541 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
e51055a0 542 }
543 trackcounter2++;
544 }
545 else {
6b7c66c9 546 mQx1 += (weight*cos(2*track->Phi())/rms[0]);
547 mQy1 += (weight*sin(2*track->Phi())/rms[1]);
e51055a0 548 if (fSaveTrackContribution){
6b7c66c9 549 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
550 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
e51055a0 551 }
ce7adfe9 552 trackcounter1++;
553 }
e51055a0 554 }
555 } else if (fSplitMethod == AliEPSelectionTask::kEta) {
556
557 for (Int_t i = 0; i < nt; i++) {
558 weight = 1;
559 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
560 if (!track) continue;
561 weight = GetWeight(track);
562 Double_t eta = track->Eta();
563 idtemp = track->GetID();
564 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
565
566 if (eta > fEtaGap/2.) {
6b7c66c9 567 mQx1 += (weight*cos(2*track->Phi())/rms[0]);
568 mQy1 += (weight*sin(2*track->Phi())/rms[1]);
e51055a0 569 if (fSaveTrackContribution){
6b7c66c9 570 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
571 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
e51055a0 572 }
573 } else if (eta < -1.*fEtaGap/2.) {
6b7c66c9 574 mQx2 += (weight*cos(2*track->Phi())/rms[0]);
575 mQy2 += (weight*sin(2*track->Phi())/rms[1]);
e51055a0 576 if (fSaveTrackContribution){
6b7c66c9 577 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
578 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
c39c35d1 579 }
580 }
581 }
582 } else if (fSplitMethod == AliEPSelectionTask::kCharge) {
583
584 for (Int_t i = 0; i < nt; i++) {
585 weight = 1;
586 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
587 if (!track) continue;
588 weight = GetWeight(track);
589 Short_t cha = track->Charge();
590 idtemp = track->GetID();
591 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
592
593 if (cha > 0) {
6b7c66c9 594 mQx1 += (weight*cos(2*track->Phi())/rms[0]);
595 mQy1 += (weight*sin(2*track->Phi())/rms[1]);
c39c35d1 596 if (fSaveTrackContribution){
6b7c66c9 597 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
598 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
c39c35d1 599 }
600 } else if (cha < 0) {
6b7c66c9 601 mQx2 += (weight*cos(2*track->Phi())/rms[0]);
602 mQy2 += (weight*sin(2*track->Phi())/rms[1]);
c39c35d1 603 if (fSaveTrackContribution){
6b7c66c9 604 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi())/rms[0],idtemp);
605 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi())/rms[1],idtemp);
e51055a0 606 }
ce7adfe9 607 }
608 }
e51055a0 609 } else {
610 printf("plane resolution determination method not available!\n\n ");
611 return;
ce7adfe9 612 }
6b7c66c9 613 // apply recenetering
614 mQ[0].Set(mQx1-(mean[0]/rms[0]), mQy1-(mean[1]/rms[1]));
615 mQ[1].Set(mQx2-(mean[0]/rms[0]), mQy2-(mean[1]/rms[1]));
ce7adfe9 616 Q1 = mQ[0];
617 Q2 = mQ[1];
618}
619
620//________________________________________________________________________
90267ca6 621void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
ce7adfe9 622
1e552991 623 if(fESDtrackCuts){
624 delete fESDtrackCuts;
625 fESDtrackCuts = 0;
626 }
e51055a0 627 if (fAnalysisInput.CompareTo("AOD")==0){
628 AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");
629 fUsercuts = kFALSE;
630 SetTrackType("TPC");
631 return;
632 }
1e552991 633 fUsercuts = kTRUE;
90267ca6 634 fESDtrackCuts = trackcuts;
ce7adfe9 635}
636
e51055a0 637//________________________________________________________________________
e1ffd90d 638void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup, Int_t ntpc){
e51055a0 639
640 if(fESDtrackCuts){
641 delete fESDtrackCuts;
642 fESDtrackCuts = 0;
643 }
644 if (fAnalysisInput.CompareTo("ESD")==0){
645 AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");
646 fUsercuts = kFALSE;
647 SetTrackType("TPC");
648 return;
649 }
650 fUsercuts = kTRUE;
651 fESDtrackCuts = new AliESDtrackCuts();
652 fESDtrackCuts->SetPtRange(ptlow,ptup);
e1ffd90d 653 fESDtrackCuts->SetMinNClustersTPC(ntpc);
e51055a0 654 fESDtrackCuts->SetEtaRange(etalow,etaup);
655 fAODfilterbit = filterbit;
656}
657
90267ca6 658//_____________________________________________________________________________
e51055a0 659
90267ca6 660void AliEPSelectionTask::SetTrackType(TString tracktype){
661 fTrackType = tracktype;
1e552991 662 if (!fUsercuts) {
e51055a0 663 if (fTrackType.CompareTo("GLOBAL")==0){
664 fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
665 fAODfilterbit = 32;
666 }
667 if (fTrackType.CompareTo("TPC")==0){
668 fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
669 fAODfilterbit = 128;
670 }
671 fESDtrackCuts->SetPtRange(0.15,20.);
90267ca6 672 fESDtrackCuts->SetEtaRange(-0.8,0.8);
ce7adfe9 673 }
ce7adfe9 674}
675
676//________________________________________________________________________
e51055a0 677Double_t AliEPSelectionTask::GetWeight(TObject* track1)
ce7adfe9 678{
679 Double_t ptweight=1;
e51055a0 680 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
5b53b816 681 if (fUsePtWeight && track) {
ce7adfe9 682 if (track->Pt()<2) ptweight=track->Pt();
683 else ptweight=2;
684 }
685 return ptweight*GetPhiWeight(track);
686}
687
688//________________________________________________________________________
e51055a0 689Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
ce7adfe9 690{
691 Double_t phiweight=1;
e51055a0 692 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
9663fa4b 693
77c17385 694 TH1F *phiDist = 0x0;
695 if(track) phiDist = SelectPhiDist(track);
ce7adfe9 696
9663fa4b 697 if (fUsePhiWeight && phiDist && track) {
698 Double_t nParticles = phiDist->Integral();
699 Double_t nPhibins = phiDist->GetNbinsX();
ce7adfe9 700
e51055a0 701 Double_t Phi = track->Phi();
ce7adfe9 702
e51055a0 703 while (Phi<0) Phi += TMath::TwoPi();
704 while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
ce7adfe9 705
9663fa4b 706 Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
ce7adfe9 707
e51055a0 708 if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
ce7adfe9 709 }
710 return phiweight;
711}
90267ca6 712
6b7c66c9 713//________________________________________________________________________
714void AliEPSelectionTask::Recenter(Int_t var, Double_t * values)
715{
716
717 if (fUseRecentering && fQDist[0] && fQDist[1] && fCentrality!=-1.) {
718 Int_t centbin = fQDist[0]->FindBin(fCentrality);
719
720 if(var==0) { // fill mean
721 values[0] = fQDist[0]->GetBinContent(centbin);
722 values[1] = fQDist[1]->GetBinContent(centbin);
723 }
724 else if(var==1) { // fill rms
725 values[0] = fQDist[0]->GetBinError(centbin);
726 values[1] = fQDist[1]->GetBinError(centbin);
727 // protection against division by zero
728 if(values[0]==0.0) values[0]=1.0;
729 if(values[1]==0.0) values[1]=1.0;
730 }
731 }
732 else { //default (no recentering)
733 values[0] = var;
734 values[1] = var;
735 }
736 return;
737}
738
90267ca6 739//__________________________________________________________________________
740void AliEPSelectionTask::SetPhiDist()
741{
9593a21e 742 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 743
9663fa4b 744 if (fPeriod.CompareTo("LHC10h")==0)
745 {
746 fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");}
747 else if(fPeriod.CompareTo("LHC11h")==0){
748 Int_t runbin=fHruns->FindBin(fRunNumber);
749 if (fHruns->GetBinContent(runbin) > 1){
750 fSparseDist->GetAxis(0)->SetRange(runbin,runbin);
751 }
752 else if(fHruns->GetBinContent(runbin) < 2){
753 fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights
754 AliInfo("Using integrated Phi-weights for this run");
755 }
756 for (Int_t i = 0; i<4 ;i++)
757 {
758 if(fPhiDist[i]){
759 delete fPhiDist[i];
760 fPhiDist[i] = 0x0;
761 }
762 if(i == 0){
763 fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge
764 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
765 if(i == 1){
766 fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge
767 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
768 if(i == 2){
769 fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge
770 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
771 if(i == 3){
772 fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge
773 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
774 fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi
775 fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber));
776 fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes
777 fSparseDist->GetAxis(2)->SetRange(1,2);
778 }
779 fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis
780 }
781
782 if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
90267ca6 783
784 }
9593a21e 785
90267ca6 786
9663fa4b 787 if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist){
788 Bool_t emptybins;
90267ca6 789
9663fa4b 790 int iter = 0;
791 while (iter<3){
792 emptybins = kFALSE;
90267ca6 793
9663fa4b 794 for (int i=1; i<fPhiDist[0]->GetNbinsX(); i++){
795 if (!((fPhiDist[0]->GetBinContent(i))>0)) {
796 emptybins = kTRUE;
797 }
798 }
799 if (emptybins) {
800 cout << "empty bins - rebinning!" << endl;
801 fPhiDist[0]->Rebin();
802 iter++;
803 }
804 else iter = 3;
805 }
90267ca6 806
9663fa4b 807 if (emptybins) {
808 AliError("After Maximum of rebinning still empty Phi-bins!!!");
809 }
90267ca6 810 }
9593a21e 811 if (fPeriod.CompareTo("LHC10h") != 0 && fPeriod.CompareTo("LHC11h") != 0 && !fUserphidist){
812 AliInfo("No Phi-weights available. All Phi weights set to 1");
813 SetUsePhiWeight(kFALSE);
814 }
90267ca6 815}
816
6b7c66c9 817//__________________________________________________________________________
818void AliEPSelectionTask::SetQvectorDist()
819{
820 if(!fUseRecentering) return;
821 AliInfo(Form("Setting q vector distributions"));
822 fQDist[0] = (TProfile*) fQxContainer->GetObject(fRunNumber, "Default");
823 fQDist[1] = (TProfile*) fQyContainer->GetObject(fRunNumber, "Default");
824
825 if (!fQDist[0] || !fQDist[1]) {
826 AliError(Form("Cannot find OADB q-vector distributions for run %d. Using default values (mean=0,rms=1).", fRunNumber));
827 return;
828 }
829
830 Bool_t emptybins;
831
832 int iter = 0;
833 while (iter<3){
834 emptybins = kFALSE;
835
836 for (int i=1; i<=fQDist[0]->GetNbinsX()-2; i++){
837 if (!((fQDist[0]->GetBinEntries(i))>0)) {
838 emptybins = kTRUE;
839 }
840 }
841 if (emptybins) {
842 cout << "empty bins - rebinning!" << endl;
843 fQDist[0]->Rebin();
844 fQDist[1]->Rebin();
845 iter++;
846 }
847 else iter = 3;
848 }
849
850 if (emptybins) {
851 AliError("After Maximum of rebinning still empty Qxy-bins!!!");
852 }
853}
854
90267ca6 855//__________________________________________________________________________
71916547 856void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
90267ca6 857{
e51055a0 858
1e552991 859 fUserphidist = kTRUE;
90267ca6 860
861 TFile f(infilename);
862 TObject* list = f.Get(listname);
9663fa4b 863 fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi");
864 if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!");
90267ca6 865
866 f.Close();
867}
e51055a0 868
e51055a0 869//_________________________________________________________________________
870TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
871{
872 TObjArray *acctracks = new TObjArray();
873
874 AliAODTrack *tr = 0;
875 Int_t maxid1 = 0;
876 Int_t maxidtemp = -1;
877 Float_t ptlow = 0;
878 Float_t ptup = 0;
879 Float_t etalow = 0;
880 Float_t etaup = 0;
881 fESDtrackCuts->GetPtRange(ptlow,ptup);
882 fESDtrackCuts->GetEtaRange(etalow,etaup);
e1ffd90d 883 Int_t ntpc = fESDtrackCuts->GetMinNClusterTPC();
e51055a0 884
885 for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
f15c1f69 886 tr = dynamic_cast<AliAODTrack*>(aod->GetTrack(i));
887 if(!tr) AliFatal("Not a standard AOD");
e51055a0 888 maxidtemp = tr->GetID();
889 if(maxidtemp < 0 && fAODfilterbit != 128) continue;
890 if(maxidtemp > -1 && fAODfilterbit == 128) continue;
891 if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
892 if (maxidtemp > maxid1) maxid1 = maxidtemp;
e1ffd90d 893 if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow && tr->GetTPCNcls() > ntpc){
e51055a0 894 acctracks->Add(tr);
895 }
896 }
897
898 maxid = maxid1;
899
900 return acctracks;
901
902}
9663fa4b 903
904//_________________________________________________________________________
905void AliEPSelectionTask::SetOADBandPeriod()
906{
907 TString oadbfilename;
908
909 if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h
910 {fPeriod = "LHC10h";
911 if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
912
913 if (fAnalysisInput.CompareTo("AOD")==0){
914 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
915 } else if (fAnalysisInput.CompareTo("ESD")==0){
916 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
917 }
918
919 TFile foadb(oadbfilename);
920 if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
921
922 AliInfo("Using Standard OADB");
923 fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
924 if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
925 foadb.Close();
926 }
927 }
928
929 if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h
930 {fPeriod = "LHC11h";
931 if (!fUserphidist) {
932 // if it's already set and custom class is required, we use the one provided by the user
933
934 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
935 TFile *foadb = TFile::Open(oadbfilename);
936 if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
937
938 AliInfo("Using Standard OADB");
939 fSparseDist = (THnSparse*) foadb->Get("Default");
940 if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
941 foadb->Close();
7cdafa1e 942 if(!fHruns){
943 fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
944 fHruns->SetName("runsHisto");
945 }
6b7c66c9 946 }
947
948 if(fUseRecentering) {
949 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/eprecentering.root", AliAnalysisManager::GetOADBPath()));
950 TFile *foadb = TFile::Open(oadbfilename);
951 if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
952
953 AliInfo("Using Standard OADB");
954 fQxContainer = (AliOADBContainer*) foadb->Get("eprecentering.Qx");
955 fQyContainer = (AliOADBContainer*) foadb->Get("eprecentering.Qy");
956 if (!fQxContainer || !fQyContainer) AliFatal("Cannot fetch OADB container for EP recentering");
957 foadb->Close();
958 }
959
9663fa4b 960 }
961}
962
963//_________________________________________________________________________
964TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track)
965{
966 if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist) return fPhiDist[0];
967 else if(fPeriod.CompareTo("LHC11h")==0)
968 {
969 if (track->Charge() < 0)
970 {
971 if(track->Eta() < 0.) return fPhiDist[0];
972 else if (track->Eta() > 0.) return fPhiDist[2];
973 }
974 else if (track->Charge() > 0)
975 {
976 if(track->Eta() < 0.) return fPhiDist[1];
977 else if (track->Eta() > 0.) return fPhiDist[3];
978 }
979
980 }
981 return 0;
982}
983
984TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd)
985{
986 // 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
987 TObjArray *acctracks = new TObjArray();
988 acctracks->SetOwner(kTRUE);
989
990 const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
991
992 Float_t ptlow = 0;
993 Float_t ptup = 0;
994 Float_t etalow = 0;
995 Float_t etaup = 0;
996 fESDtrackCuts->GetPtRange(ptlow,ptup);
997 fESDtrackCuts->GetEtaRange(etalow,etaup);
998
999 Double_t pDCA[3] = { 0. }; // momentum at DCA
1000 Double_t rDCA[3] = { 0. }; // position at DCA
1001 Float_t dDCA[2] = {0.}; // DCA to the vertex d and z
1002 Float_t cDCA[3] = {0.}; // covariance of impact parameters
1003
1004
1005
1006 for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){
1007
1008 AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise need to work with a copy
1009 //
1010 // Track selection
1011 if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue;
1012
1013 // create a tpc only tracl
1014 AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
1015 if(!track) continue;
1016
1017 if(track->Pt()>0.)
1018 {
1019 // only constrain tracks above threshold
1020 AliExternalTrackParam exParam;
1021 // take the B-field from the ESD, no 3D fieldMap available at this point
1022 Bool_t relate = false;
1023 relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
1024 if(!relate){
1025 delete track;
1026 continue;
1027 }
1028 // fetch the track parameters at the DCA (unconstraint)
1029 if(track->GetTPCInnerParam()){
1030 track->GetTPCInnerParam()->GetPxPyPz(pDCA);
1031 track->GetTPCInnerParam()->GetXYZ(rDCA);
1032 }
1033 // get the DCA to the vertex:
1034 track->GetImpactParametersTPC(dDCA,cDCA);
1035 // set the constrained parameters to the track
1036 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1037 }
1038
1039
1040 Float_t eta = track->Eta();
1041 Float_t pT = track->Pt();
1042
1043 if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){
1044 delete track;
1045 continue;
1046 }
1047
1048 acctracks->Add(track);
1049 }
1050
1051
1052 return acctracks;
1053
1054}
1055