Correct the not needed conversion of phi angle from radians to degrees
[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 }
77c17385 178 if (fPeriod.CompareTo("LHC11h")==0){
9663fa4b 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
77c17385 616 TH1F *phiDist = 0x0;
617 if(track) phiDist = SelectPhiDist(track);
ce7adfe9 618
9663fa4b 619 if (fUsePhiWeight && phiDist && track) {
620 Double_t nParticles = phiDist->Integral();
621 Double_t nPhibins = phiDist->GetNbinsX();
ce7adfe9 622
e51055a0 623 Double_t Phi = track->Phi();
ce7adfe9 624
e51055a0 625 while (Phi<0) Phi += TMath::TwoPi();
626 while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
ce7adfe9 627
9663fa4b 628 Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
ce7adfe9 629
e51055a0 630 if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
ce7adfe9 631 }
632 return phiweight;
633}
90267ca6 634
635//__________________________________________________________________________
636void AliEPSelectionTask::SetPhiDist()
637{
1e552991 638 if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
90267ca6 639
9663fa4b 640 if (fPeriod.CompareTo("LHC10h")==0)
641 {
642 fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");}
643 else if(fPeriod.CompareTo("LHC11h")==0){
644 Int_t runbin=fHruns->FindBin(fRunNumber);
645 if (fHruns->GetBinContent(runbin) > 1){
646 fSparseDist->GetAxis(0)->SetRange(runbin,runbin);
647 }
648 else if(fHruns->GetBinContent(runbin) < 2){
649 fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights
650 AliInfo("Using integrated Phi-weights for this run");
651 }
652 for (Int_t i = 0; i<4 ;i++)
653 {
654 if(fPhiDist[i]){
655 delete fPhiDist[i];
656 fPhiDist[i] = 0x0;
657 }
658 if(i == 0){
659 fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge
660 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
661 if(i == 1){
662 fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge
663 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
664 if(i == 2){
665 fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge
666 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
667 if(i == 3){
668 fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge
669 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
670 fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi
671 fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber));
672 fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes
673 fSparseDist->GetAxis(2)->SetRange(1,2);
674 }
675 fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis
676 }
677
678 if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
90267ca6 679
680 }
681 else {
682 AliInfo("Using Custom Phi Distribution");
683 }
684
9663fa4b 685 if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist){
686 Bool_t emptybins;
90267ca6 687
9663fa4b 688 int iter = 0;
689 while (iter<3){
690 emptybins = kFALSE;
90267ca6 691
9663fa4b 692 for (int i=1; i<fPhiDist[0]->GetNbinsX(); i++){
693 if (!((fPhiDist[0]->GetBinContent(i))>0)) {
694 emptybins = kTRUE;
695 }
696 }
697 if (emptybins) {
698 cout << "empty bins - rebinning!" << endl;
699 fPhiDist[0]->Rebin();
700 iter++;
701 }
702 else iter = 3;
703 }
90267ca6 704
9663fa4b 705 if (emptybins) {
706 AliError("After Maximum of rebinning still empty Phi-bins!!!");
707 }
90267ca6 708 }
709}
710
711//__________________________________________________________________________
71916547 712void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
90267ca6 713{
e51055a0 714
1e552991 715 fUserphidist = kTRUE;
90267ca6 716
717 TFile f(infilename);
718 TObject* list = f.Get(listname);
9663fa4b 719 fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi");
720 if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!");
90267ca6 721
722 f.Close();
723}
e51055a0 724
725
726//_________________________________________________________________________
727TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
728{
729 TObjArray *acctracks = new TObjArray();
730
731 AliAODTrack *tr = 0;
732 Int_t maxid1 = 0;
733 Int_t maxidtemp = -1;
734 Float_t ptlow = 0;
735 Float_t ptup = 0;
736 Float_t etalow = 0;
737 Float_t etaup = 0;
738 fESDtrackCuts->GetPtRange(ptlow,ptup);
739 fESDtrackCuts->GetEtaRange(etalow,etaup);
740
741 for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
742 tr = aod->GetTrack(i);
743 maxidtemp = tr->GetID();
744 if(maxidtemp < 0 && fAODfilterbit != 128) continue;
745 if(maxidtemp > -1 && fAODfilterbit == 128) continue;
746 if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
747 if (maxidtemp > maxid1) maxid1 = maxidtemp;
748 if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow){
749 acctracks->Add(tr);
750 }
751 }
752
753 maxid = maxid1;
754
755 return acctracks;
756
757}
9663fa4b 758
759//_________________________________________________________________________
760void AliEPSelectionTask::SetOADBandPeriod()
761{
762 TString oadbfilename;
763
764 if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h
765 {fPeriod = "LHC10h";
766 if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
767
768 if (fAnalysisInput.CompareTo("AOD")==0){
769 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
770 } else if (fAnalysisInput.CompareTo("ESD")==0){
771 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
772 }
773
774 TFile foadb(oadbfilename);
775 if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
776
777 AliInfo("Using Standard OADB");
778 fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
779 if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
780 foadb.Close();
781 }
782 }
783
784 if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h
785 {fPeriod = "LHC11h";
786 if (!fUserphidist) {
787 // if it's already set and custom class is required, we use the one provided by the user
788
789 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
790 TFile *foadb = TFile::Open(oadbfilename);
791 if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
792
793 AliInfo("Using Standard OADB");
794 fSparseDist = (THnSparse*) foadb->Get("Default");
795 if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
796 foadb->Close();
7cdafa1e 797 if(!fHruns){
798 fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
799 fHruns->SetName("runsHisto");
800 }
9663fa4b 801 }
802 }
803}
804
805//_________________________________________________________________________
806TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track)
807{
808 if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist) return fPhiDist[0];
809 else if(fPeriod.CompareTo("LHC11h")==0)
810 {
811 if (track->Charge() < 0)
812 {
813 if(track->Eta() < 0.) return fPhiDist[0];
814 else if (track->Eta() > 0.) return fPhiDist[2];
815 }
816 else if (track->Charge() > 0)
817 {
818 if(track->Eta() < 0.) return fPhiDist[1];
819 else if (track->Eta() > 0.) return fPhiDist[3];
820 }
821
822 }
823 return 0;
824}
825
826TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd)
827{
828 // 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
829 TObjArray *acctracks = new TObjArray();
830 acctracks->SetOwner(kTRUE);
831
832 const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
833
834 Float_t ptlow = 0;
835 Float_t ptup = 0;
836 Float_t etalow = 0;
837 Float_t etaup = 0;
838 fESDtrackCuts->GetPtRange(ptlow,ptup);
839 fESDtrackCuts->GetEtaRange(etalow,etaup);
840
841 Double_t pDCA[3] = { 0. }; // momentum at DCA
842 Double_t rDCA[3] = { 0. }; // position at DCA
843 Float_t dDCA[2] = {0.}; // DCA to the vertex d and z
844 Float_t cDCA[3] = {0.}; // covariance of impact parameters
845
846
847
848 for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){
849
850 AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise need to work with a copy
851 //
852 // Track selection
853 if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue;
854
855 // create a tpc only tracl
856 AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
857 if(!track) continue;
858
859 if(track->Pt()>0.)
860 {
861 // only constrain tracks above threshold
862 AliExternalTrackParam exParam;
863 // take the B-field from the ESD, no 3D fieldMap available at this point
864 Bool_t relate = false;
865 relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
866 if(!relate){
867 delete track;
868 continue;
869 }
870 // fetch the track parameters at the DCA (unconstraint)
871 if(track->GetTPCInnerParam()){
872 track->GetTPCInnerParam()->GetPxPyPz(pDCA);
873 track->GetTPCInnerParam()->GetXYZ(rDCA);
874 }
875 // get the DCA to the vertex:
876 track->GetImpactParametersTPC(dDCA,cDCA);
877 // set the constrained parameters to the track
878 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
879 }
880
881
882 Float_t eta = track->Eta();
883 Float_t pT = track->Pt();
884
885 if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){
886 delete track;
887 continue;
888 }
889
890 acctracks->Add(track);
891 }
892
893
894 return acctracks;
895
896}
897