]> 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
Added job information to fileinfo.log: total duration, CPU time reported by ROOT...
[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();
5b53b816 335 esdEP->Reset();
e51055a0 336
5b53b816 337 Int_t maxID = 0;
338 TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID);
e51055a0 339
e24b883a 340 if (fSaveTrackContribution) {
341 esdEP->GetQContributionXArray()->Set(maxID+1);
342 esdEP->GetQContributionYArray()->Set(maxID+1);
343 esdEP->GetQContributionXArraysub1()->Set(maxID+1);
344 esdEP->GetQContributionYArraysub1()->Set(maxID+1);
345 esdEP->GetQContributionXArraysub2()->Set(maxID+1);
346 esdEP->GetQContributionYArraysub2()->Set(maxID+1);
347 }
e51055a0 348
e24b883a 349 const int NT = tracklist->GetEntries();
e51055a0 350
e24b883a 351 if (NT>4){
352 fQVector = new TVector2(GetQ(esdEP,tracklist));
353 fEventplaneQ = fQVector->Phi()/2.;
354 GetQsub(qq1, qq2, tracklist, esdEP);
355 fQsub1 = new TVector2(qq1);
356 fQsub2 = new TVector2(qq2);
357 fQsubRes = (fQsub1->Phi()/2. - fQsub2->Phi()/2.);
e51055a0 358
e24b883a 359 esdEP->SetQVector(fQVector);
360 esdEP->SetEventplaneQ(fEventplaneQ);
361 esdEP->SetQsub(fQsub1,fQsub2);
362 esdEP->SetQsubRes(fQsubRes);
363
364 fHOutEventplaneQ->Fill(fEventplaneQ);
365 fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.);
366 fHOutNTEPRes->Fill(NT,fQsubRes);
e51055a0 367
e51055a0 368 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
369
370 for (int iter = 0; iter<NT;iter++){
371 AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
372 if (track) {
373 float delta = track->Phi()-fEventplaneQ;
374 while (delta < 0) delta += TMath::Pi();
375 while (delta > TMath::Pi()) delta -= TMath::Pi();
376 fHOutPTPsi->Fill(track->Pt(),delta);
377 fHOutPhi->Fill(track->Phi());
378 fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
379 }
380 }
381
382 AliAODTrack* trmax = aod->GetTrack(0);
383 for (int iter = 1; iter<NT;iter++){
384 AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
385 if (track && (track->Pt() > trmax->Pt())) trmax = track;
386 }
387 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
388 }
389 delete tracklist;
390 tracklist = 0;
391 }
392
393
ce7adfe9 394 }
e51055a0 395
396
ce7adfe9 397 else {
398 printf(" Analysis Input not known!\n\n ");
399 return;
400 }
401 PostData(1, fOutputList);
402}
403
404//________________________________________________________________________
405void AliEPSelectionTask::Terminate(Option_t */*option*/)
406{
407 // Terminate analysis
408}
409
410//__________________________________________________________________________
411TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
412{
71916547 413// Get the Q vector
ce7adfe9 414 TVector2 mQ;
415 float mQx=0, mQy=0;
e51055a0 416 AliVTrack* track;
ce7adfe9 417 Double_t weight;
e51055a0 418 Int_t idtemp = -1;
ce7adfe9 419
71916547 420 int nt = tracklist->GetEntries();
ce7adfe9 421
71916547 422 for (int i=0; i<nt; i++){
ce7adfe9 423 weight = 1;
e51055a0 424 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
b1a983b4 425 if (track) {
426 weight = GetWeight(track);
e51055a0 427 if (fSaveTrackContribution){
428 idtemp = track->GetID();
429 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
430 EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),idtemp);
431 EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),idtemp);
432 }
433 mQx += (weight*cos(2*track->Phi()));
434 mQy += (weight*sin(2*track->Phi()));
ce7adfe9 435 }
ce7adfe9 436 }
437 mQ.Set(mQx,mQy);
438 return mQ;
439}
440
441 //________________________________________________________________________
e51055a0 442void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP)
ce7adfe9 443{
71916547 444// Get Qsub
ce7adfe9 445 TVector2 mQ[2];
446 float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
447 Double_t weight;
448
e51055a0 449 AliVTrack* track;
71916547 450 TRandom2 rn = 0;
ce7adfe9 451
71916547 452 int nt = tracklist->GetEntries();
ce7adfe9 453 int trackcounter1=0, trackcounter2=0;
e51055a0 454 int idtemp = 0;
455
456 if (fSplitMethod == AliEPSelectionTask::kRandom){
457
458 for (Int_t i = 0; i < nt; i++) {
459 weight = 1;
460 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
461 if (!track) continue;
462 weight = GetWeight(track);
463 idtemp = track->GetID();
464 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
ce7adfe9 465
e51055a0 466 // This loop splits the track set into 2 random subsets
467 if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
468 float random = rn.Rndm();
469 if(random < .5){
470 mQx1 += (weight*cos(2*track->Phi()));
471 mQy1 += (weight*sin(2*track->Phi()));
472 if (fSaveTrackContribution){
473 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
474 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
475 }
476 trackcounter1++;
477 }
478 else {
479 mQx2 += (weight*cos(2*track->Phi()));
480 mQy2 += (weight*sin(2*track->Phi()));
481 if (fSaveTrackContribution){
482 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
483 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
484 }
485 trackcounter2++;
486 }
487 }
488 else if( trackcounter1 >= int(nt/2.)){
489 mQx2 += (weight*cos(2*track->Phi()));
490 mQy2 += (weight*sin(2*track->Phi()));
491 if (fSaveTrackContribution){
492 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
493 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
494 }
495 trackcounter2++;
496 }
497 else {
ce7adfe9 498 mQx1 += (weight*cos(2*track->Phi()));
499 mQy1 += (weight*sin(2*track->Phi()));
e51055a0 500 if (fSaveTrackContribution){
501 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
502 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
503 }
ce7adfe9 504 trackcounter1++;
505 }
e51055a0 506 }
507 } else if (fSplitMethod == AliEPSelectionTask::kEta) {
508
509 for (Int_t i = 0; i < nt; i++) {
510 weight = 1;
511 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
512 if (!track) continue;
513 weight = GetWeight(track);
514 Double_t eta = track->Eta();
515 idtemp = track->GetID();
516 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
517
518 if (eta > fEtaGap/2.) {
519 mQx1 += (weight*cos(2*track->Phi()));
520 mQy1 += (weight*sin(2*track->Phi()));
521 if (fSaveTrackContribution){
522 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
523 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
524 }
525 } else if (eta < -1.*fEtaGap/2.) {
ce7adfe9 526 mQx2 += (weight*cos(2*track->Phi()));
527 mQy2 += (weight*sin(2*track->Phi()));
e51055a0 528 if (fSaveTrackContribution){
529 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
530 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
531 }
ce7adfe9 532 }
533 }
e51055a0 534 } else {
535 printf("plane resolution determination method not available!\n\n ");
536 return;
ce7adfe9 537 }
e51055a0 538
ce7adfe9 539 mQ[0].Set(mQx1,mQy1);
540 mQ[1].Set(mQx2,mQy2);
541 Q1 = mQ[0];
542 Q2 = mQ[1];
543}
544
545//________________________________________________________________________
90267ca6 546void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
ce7adfe9 547
1e552991 548 if(fESDtrackCuts){
549 delete fESDtrackCuts;
550 fESDtrackCuts = 0;
551 }
e51055a0 552 if (fAnalysisInput.CompareTo("AOD")==0){
553 AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");
554 fUsercuts = kFALSE;
555 SetTrackType("TPC");
556 return;
557 }
1e552991 558 fUsercuts = kTRUE;
90267ca6 559 fESDtrackCuts = trackcuts;
ce7adfe9 560}
561
e51055a0 562//________________________________________________________________________
563void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup){
564
565 if(fESDtrackCuts){
566 delete fESDtrackCuts;
567 fESDtrackCuts = 0;
568 }
569 if (fAnalysisInput.CompareTo("ESD")==0){
570 AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");
571 fUsercuts = kFALSE;
572 SetTrackType("TPC");
573 return;
574 }
575 fUsercuts = kTRUE;
576 fESDtrackCuts = new AliESDtrackCuts();
577 fESDtrackCuts->SetPtRange(ptlow,ptup);
578 fESDtrackCuts->SetEtaRange(etalow,etaup);
579 fAODfilterbit = filterbit;
580}
581
90267ca6 582//_____________________________________________________________________________
e51055a0 583
90267ca6 584void AliEPSelectionTask::SetTrackType(TString tracktype){
585 fTrackType = tracktype;
1e552991 586 if (!fUsercuts) {
e51055a0 587 if (fTrackType.CompareTo("GLOBAL")==0){
588 fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
589 fAODfilterbit = 32;
590 }
591 if (fTrackType.CompareTo("TPC")==0){
592 fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
593 fAODfilterbit = 128;
594 }
595 fESDtrackCuts->SetPtRange(0.15,20.);
90267ca6 596 fESDtrackCuts->SetEtaRange(-0.8,0.8);
ce7adfe9 597 }
ce7adfe9 598}
599
600//________________________________________________________________________
e51055a0 601Double_t AliEPSelectionTask::GetWeight(TObject* track1)
ce7adfe9 602{
603 Double_t ptweight=1;
e51055a0 604 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
5b53b816 605 if (fUsePtWeight && track) {
ce7adfe9 606 if (track->Pt()<2) ptweight=track->Pt();
607 else ptweight=2;
608 }
609 return ptweight*GetPhiWeight(track);
610}
611
612//________________________________________________________________________
e51055a0 613Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
ce7adfe9 614{
615 Double_t phiweight=1;
e51055a0 616 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
9663fa4b 617
77c17385 618 TH1F *phiDist = 0x0;
619 if(track) phiDist = SelectPhiDist(track);
ce7adfe9 620
9663fa4b 621 if (fUsePhiWeight && phiDist && track) {
622 Double_t nParticles = phiDist->Integral();
623 Double_t nPhibins = phiDist->GetNbinsX();
ce7adfe9 624
e51055a0 625 Double_t Phi = track->Phi();
ce7adfe9 626
e51055a0 627 while (Phi<0) Phi += TMath::TwoPi();
628 while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
ce7adfe9 629
9663fa4b 630 Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
ce7adfe9 631
e51055a0 632 if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
ce7adfe9 633 }
634 return phiweight;
635}
90267ca6 636
637//__________________________________________________________________________
638void AliEPSelectionTask::SetPhiDist()
639{
9593a21e 640 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 641
9663fa4b 642 if (fPeriod.CompareTo("LHC10h")==0)
643 {
644 fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");}
645 else if(fPeriod.CompareTo("LHC11h")==0){
646 Int_t runbin=fHruns->FindBin(fRunNumber);
647 if (fHruns->GetBinContent(runbin) > 1){
648 fSparseDist->GetAxis(0)->SetRange(runbin,runbin);
649 }
650 else if(fHruns->GetBinContent(runbin) < 2){
651 fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights
652 AliInfo("Using integrated Phi-weights for this run");
653 }
654 for (Int_t i = 0; i<4 ;i++)
655 {
656 if(fPhiDist[i]){
657 delete fPhiDist[i];
658 fPhiDist[i] = 0x0;
659 }
660 if(i == 0){
661 fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge
662 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
663 if(i == 1){
664 fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge
665 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
666 if(i == 2){
667 fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge
668 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
669 if(i == 3){
670 fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge
671 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
672 fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi
673 fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber));
674 fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes
675 fSparseDist->GetAxis(2)->SetRange(1,2);
676 }
677 fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis
678 }
679
680 if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
90267ca6 681
682 }
9593a21e 683
90267ca6 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 }
9593a21e 709 if (fPeriod.CompareTo("LHC10h") != 0 && fPeriod.CompareTo("LHC11h") != 0 && !fUserphidist){
710 AliInfo("No Phi-weights available. All Phi weights set to 1");
711 SetUsePhiWeight(kFALSE);
712 }
90267ca6 713}
714
715//__________________________________________________________________________
71916547 716void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
90267ca6 717{
e51055a0 718
1e552991 719 fUserphidist = kTRUE;
90267ca6 720
721 TFile f(infilename);
722 TObject* list = f.Get(listname);
9663fa4b 723 fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi");
724 if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!");
90267ca6 725
726 f.Close();
727}
e51055a0 728
729
730//_________________________________________________________________________
731TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
732{
733 TObjArray *acctracks = new TObjArray();
734
735 AliAODTrack *tr = 0;
736 Int_t maxid1 = 0;
737 Int_t maxidtemp = -1;
738 Float_t ptlow = 0;
739 Float_t ptup = 0;
740 Float_t etalow = 0;
741 Float_t etaup = 0;
742 fESDtrackCuts->GetPtRange(ptlow,ptup);
743 fESDtrackCuts->GetEtaRange(etalow,etaup);
744
745 for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
746 tr = aod->GetTrack(i);
747 maxidtemp = tr->GetID();
748 if(maxidtemp < 0 && fAODfilterbit != 128) continue;
749 if(maxidtemp > -1 && fAODfilterbit == 128) continue;
750 if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
751 if (maxidtemp > maxid1) maxid1 = maxidtemp;
752 if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow){
753 acctracks->Add(tr);
754 }
755 }
756
757 maxid = maxid1;
758
759 return acctracks;
760
761}
9663fa4b 762
763//_________________________________________________________________________
764void AliEPSelectionTask::SetOADBandPeriod()
765{
766 TString oadbfilename;
767
768 if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h
769 {fPeriod = "LHC10h";
770 if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
771
772 if (fAnalysisInput.CompareTo("AOD")==0){
773 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
774 } else if (fAnalysisInput.CompareTo("ESD")==0){
775 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
776 }
777
778 TFile foadb(oadbfilename);
779 if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
780
781 AliInfo("Using Standard OADB");
782 fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
783 if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
784 foadb.Close();
785 }
786 }
787
788 if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h
789 {fPeriod = "LHC11h";
790 if (!fUserphidist) {
791 // if it's already set and custom class is required, we use the one provided by the user
792
793 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
794 TFile *foadb = TFile::Open(oadbfilename);
795 if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
796
797 AliInfo("Using Standard OADB");
798 fSparseDist = (THnSparse*) foadb->Get("Default");
799 if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
800 foadb->Close();
7cdafa1e 801 if(!fHruns){
802 fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
803 fHruns->SetName("runsHisto");
804 }
9663fa4b 805 }
806 }
807}
808
809//_________________________________________________________________________
810TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track)
811{
812 if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist) return fPhiDist[0];
813 else if(fPeriod.CompareTo("LHC11h")==0)
814 {
815 if (track->Charge() < 0)
816 {
817 if(track->Eta() < 0.) return fPhiDist[0];
818 else if (track->Eta() > 0.) return fPhiDist[2];
819 }
820 else if (track->Charge() > 0)
821 {
822 if(track->Eta() < 0.) return fPhiDist[1];
823 else if (track->Eta() > 0.) return fPhiDist[3];
824 }
825
826 }
827 return 0;
828}
829
830TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd)
831{
832 // 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
833 TObjArray *acctracks = new TObjArray();
834 acctracks->SetOwner(kTRUE);
835
836 const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
837
838 Float_t ptlow = 0;
839 Float_t ptup = 0;
840 Float_t etalow = 0;
841 Float_t etaup = 0;
842 fESDtrackCuts->GetPtRange(ptlow,ptup);
843 fESDtrackCuts->GetEtaRange(etalow,etaup);
844
845 Double_t pDCA[3] = { 0. }; // momentum at DCA
846 Double_t rDCA[3] = { 0. }; // position at DCA
847 Float_t dDCA[2] = {0.}; // DCA to the vertex d and z
848 Float_t cDCA[3] = {0.}; // covariance of impact parameters
849
850
851
852 for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){
853
854 AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise need to work with a copy
855 //
856 // Track selection
857 if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue;
858
859 // create a tpc only tracl
860 AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
861 if(!track) continue;
862
863 if(track->Pt()>0.)
864 {
865 // only constrain tracks above threshold
866 AliExternalTrackParam exParam;
867 // take the B-field from the ESD, no 3D fieldMap available at this point
868 Bool_t relate = false;
869 relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
870 if(!relate){
871 delete track;
872 continue;
873 }
874 // fetch the track parameters at the DCA (unconstraint)
875 if(track->GetTPCInnerParam()){
876 track->GetTPCInnerParam()->GetPxPyPz(pDCA);
877 track->GetTPCInnerParam()->GetXYZ(rDCA);
878 }
879 // get the DCA to the vertex:
880 track->GetImpactParametersTPC(dDCA,cDCA);
881 // set the constrained parameters to the track
882 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
883 }
884
885
886 Float_t eta = track->Eta();
887 Float_t pT = track->Pt();
888
889 if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){
890 delete track;
891 continue;
892 }
893
894 acctracks->Add(track);
895 }
896
897
898 return acctracks;
899
900}
901