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