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