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