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