git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
TPC cluster cut on AODs
[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);
c39c35d1 531 }
532 }
533 }
534 } else if (fSplitMethod == AliEPSelectionTask::kCharge) {
535
536 for (Int_t i = 0; i < nt; i++) {
537 weight = 1;
538 track = dynamic_cast<AliVTrack*> (tracklist->At(i));
539 if (!track) continue;
540 weight = GetWeight(track);
541 Short_t cha = track->Charge();
542 idtemp = track->GetID();
543 if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
544
545 if (cha > 0) {
546 mQx1 += (weight*cos(2*track->Phi()));
547 mQy1 += (weight*sin(2*track->Phi()));
548 if (fSaveTrackContribution){
549 EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
550 EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
551 }
552 } else if (cha < 0) {
553 mQx2 += (weight*cos(2*track->Phi()));
554 mQy2 += (weight*sin(2*track->Phi()));
555 if (fSaveTrackContribution){
556 EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
557 EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
e51055a0 558 }
ce7adfe9 559 }
560 }
e51055a0 561 } else {
562 printf("plane resolution determination method not available!\n\n ");
563 return;
ce7adfe9 564 }
e51055a0 565
ce7adfe9 566 mQ[0].Set(mQx1,mQy1);
567 mQ[1].Set(mQx2,mQy2);
568 Q1 = mQ[0];
569 Q2 = mQ[1];
570}
571
572//________________________________________________________________________
90267ca6 573void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
ce7adfe9 574
1e552991 575 if(fESDtrackCuts){
576 delete fESDtrackCuts;
577 fESDtrackCuts = 0;
578 }
e51055a0 579 if (fAnalysisInput.CompareTo("AOD")==0){
580 AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");
581 fUsercuts = kFALSE;
582 SetTrackType("TPC");
583 return;
584 }
1e552991 585 fUsercuts = kTRUE;
90267ca6 586 fESDtrackCuts = trackcuts;
ce7adfe9 587}
588
e51055a0 589//________________________________________________________________________
e1ffd90d 590void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup, Int_t ntpc){
e51055a0 591
592 if(fESDtrackCuts){
593 delete fESDtrackCuts;
594 fESDtrackCuts = 0;
595 }
596 if (fAnalysisInput.CompareTo("ESD")==0){
597 AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");
598 fUsercuts = kFALSE;
599 SetTrackType("TPC");
600 return;
601 }
602 fUsercuts = kTRUE;
603 fESDtrackCuts = new AliESDtrackCuts();
604 fESDtrackCuts->SetPtRange(ptlow,ptup);
e1ffd90d 605 fESDtrackCuts->SetMinNClustersTPC(ntpc);
e51055a0 606 fESDtrackCuts->SetEtaRange(etalow,etaup);
607 fAODfilterbit = filterbit;
608}
609
90267ca6 610//_____________________________________________________________________________
e51055a0 611
90267ca6 612void AliEPSelectionTask::SetTrackType(TString tracktype){
613 fTrackType = tracktype;
1e552991 614 if (!fUsercuts) {
e51055a0 615 if (fTrackType.CompareTo("GLOBAL")==0){
616 fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
617 fAODfilterbit = 32;
618 }
619 if (fTrackType.CompareTo("TPC")==0){
620 fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
621 fAODfilterbit = 128;
622 }
623 fESDtrackCuts->SetPtRange(0.15,20.);
90267ca6 624 fESDtrackCuts->SetEtaRange(-0.8,0.8);
ce7adfe9 625 }
ce7adfe9 626}
627
628//________________________________________________________________________
e51055a0 629Double_t AliEPSelectionTask::GetWeight(TObject* track1)
ce7adfe9 630{
631 Double_t ptweight=1;
e51055a0 632 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
5b53b816 633 if (fUsePtWeight && track) {
ce7adfe9 634 if (track->Pt()<2) ptweight=track->Pt();
635 else ptweight=2;
636 }
637 return ptweight*GetPhiWeight(track);
638}
639
640//________________________________________________________________________
e51055a0 641Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
ce7adfe9 642{
643 Double_t phiweight=1;
e51055a0 644 AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
9663fa4b 645
77c17385 646 TH1F *phiDist = 0x0;
647 if(track) phiDist = SelectPhiDist(track);
ce7adfe9 648
9663fa4b 649 if (fUsePhiWeight && phiDist && track) {
650 Double_t nParticles = phiDist->Integral();
651 Double_t nPhibins = phiDist->GetNbinsX();
ce7adfe9 652
e51055a0 653 Double_t Phi = track->Phi();
ce7adfe9 654
e51055a0 655 while (Phi<0) Phi += TMath::TwoPi();
656 while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
ce7adfe9 657
9663fa4b 658 Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
ce7adfe9 659
e51055a0 660 if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
ce7adfe9 661 }
662 return phiweight;
663}
90267ca6 664
665//__________________________________________________________________________
666void AliEPSelectionTask::SetPhiDist()
667{
9593a21e 668 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 669
9663fa4b 670 if (fPeriod.CompareTo("LHC10h")==0)
671 {
672 fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");}
673 else if(fPeriod.CompareTo("LHC11h")==0){
674 Int_t runbin=fHruns->FindBin(fRunNumber);
675 if (fHruns->GetBinContent(runbin) > 1){
676 fSparseDist->GetAxis(0)->SetRange(runbin,runbin);
677 }
678 else if(fHruns->GetBinContent(runbin) < 2){
679 fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights
680 AliInfo("Using integrated Phi-weights for this run");
681 }
682 for (Int_t i = 0; i<4 ;i++)
683 {
684 if(fPhiDist[i]){
685 delete fPhiDist[i];
686 fPhiDist[i] = 0x0;
687 }
688 if(i == 0){
689 fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge
690 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
691 if(i == 1){
692 fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge
693 fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
694 if(i == 2){
695 fSparseDist->GetAxis(1)->SetRange(1,1); // neg charge
696 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
697 if(i == 3){
698 fSparseDist->GetAxis(1)->SetRange(2,2); // pos charge
699 fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
700 fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi
701 fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber));
702 fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes
703 fSparseDist->GetAxis(2)->SetRange(1,2);
704 }
705 fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis
706 }
707
708 if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
90267ca6 709
710 }
9593a21e 711
90267ca6 712
9663fa4b 713 if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist){
714 Bool_t emptybins;
90267ca6 715
9663fa4b 716 int iter = 0;
717 while (iter<3){
718 emptybins = kFALSE;
90267ca6 719
9663fa4b 720 for (int i=1; i<fPhiDist[0]->GetNbinsX(); i++){
721 if (!((fPhiDist[0]->GetBinContent(i))>0)) {
722 emptybins = kTRUE;
723 }
724 }
725 if (emptybins) {
726 cout << "empty bins - rebinning!" << endl;
727 fPhiDist[0]->Rebin();
728 iter++;
729 }
730 else iter = 3;
731 }
90267ca6 732
9663fa4b 733 if (emptybins) {
734 AliError("After Maximum of rebinning still empty Phi-bins!!!");
735 }
90267ca6 736 }
9593a21e 737 if (fPeriod.CompareTo("LHC10h") != 0 && fPeriod.CompareTo("LHC11h") != 0 && !fUserphidist){
738 AliInfo("No Phi-weights available. All Phi weights set to 1");
739 SetUsePhiWeight(kFALSE);
740 }
90267ca6 741}
742
743//__________________________________________________________________________
71916547 744void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
90267ca6 745{
e51055a0 746
1e552991 747 fUserphidist = kTRUE;
90267ca6 748
749 TFile f(infilename);
750 TObject* list = f.Get(listname);
9663fa4b 751 fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi");
752 if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!");
90267ca6 753
754 f.Close();
755}
e51055a0 756
757
758//_________________________________________________________________________
759TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
760{
761 TObjArray *acctracks = new TObjArray();
762
763 AliAODTrack *tr = 0;
764 Int_t maxid1 = 0;
765 Int_t maxidtemp = -1;
766 Float_t ptlow = 0;
767 Float_t ptup = 0;
768 Float_t etalow = 0;
769 Float_t etaup = 0;
770 fESDtrackCuts->GetPtRange(ptlow,ptup);
771 fESDtrackCuts->GetEtaRange(etalow,etaup);
e1ffd90d 772 Int_t ntpc = fESDtrackCuts->GetMinNClusterTPC();
e51055a0 773
774 for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
775 tr = aod->GetTrack(i);
776 maxidtemp = tr->GetID();
777 if(maxidtemp < 0 && fAODfilterbit != 128) continue;
778 if(maxidtemp > -1 && fAODfilterbit == 128) continue;
779 if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
780 if (maxidtemp > maxid1) maxid1 = maxidtemp;
e1ffd90d 781 if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow && tr->GetTPCNcls() > ntpc){
e51055a0 782 acctracks->Add(tr);
783 }
784 }
785
786 maxid = maxid1;
787
788 return acctracks;
789
790}
9663fa4b 791
792//_________________________________________________________________________
793void AliEPSelectionTask::SetOADBandPeriod()
794{
795 TString oadbfilename;
796
797 if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h
798 {fPeriod = "LHC10h";
799 if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
800
801 if (fAnalysisInput.CompareTo("AOD")==0){
802 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
803 } else if (fAnalysisInput.CompareTo("ESD")==0){
804 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
805 }
806
807 TFile foadb(oadbfilename);
808 if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
809
810 AliInfo("Using Standard OADB");
811 fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
812 if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
813 foadb.Close();
814 }
815 }
816
817 if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h
818 {fPeriod = "LHC11h";
819 if (!fUserphidist) {
820 // if it's already set and custom class is required, we use the one provided by the user
821
822 oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
823 TFile *foadb = TFile::Open(oadbfilename);
824 if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
825
826 AliInfo("Using Standard OADB");
827 fSparseDist = (THnSparse*) foadb->Get("Default");
828 if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
829 foadb->Close();
7cdafa1e 830 if(!fHruns){
831 fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
832 fHruns->SetName("runsHisto");
833 }
9663fa4b 834 }
835 }
836}
837
838//_________________________________________________________________________
839TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track)
840{
841 if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist) return fPhiDist[0];
842 else if(fPeriod.CompareTo("LHC11h")==0)
843 {
844 if (track->Charge() < 0)
845 {
846 if(track->Eta() < 0.) return fPhiDist[0];
847 else if (track->Eta() > 0.) return fPhiDist[2];
848 }
849 else if (track->Charge() > 0)
850 {
851 if(track->Eta() < 0.) return fPhiDist[1];
852 else if (track->Eta() > 0.) return fPhiDist[3];
853 }
854
855 }
856 return 0;
857}
858
859TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd)
860{
861 // 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
862 TObjArray *acctracks = new TObjArray();
863 acctracks->SetOwner(kTRUE);
864
865 const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
866
867 Float_t ptlow = 0;
868 Float_t ptup = 0;
869 Float_t etalow = 0;
870 Float_t etaup = 0;
871 fESDtrackCuts->GetPtRange(ptlow,ptup);
872 fESDtrackCuts->GetEtaRange(etalow,etaup);
873
874 Double_t pDCA[3] = { 0. }; // momentum at DCA
875 Double_t rDCA[3] = { 0. }; // position at DCA
876 Float_t dDCA[2] = {0.}; // DCA to the vertex d and z
877 Float_t cDCA[3] = {0.}; // covariance of impact parameters
878
879
880
881 for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){
882
883 AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise need to work with a copy
884 //
885 // Track selection
886 if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue;
887
888 // create a tpc only tracl
889 AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
890 if(!track) continue;
891
892 if(track->Pt()>0.)
893 {
894 // only constrain tracks above threshold
895 AliExternalTrackParam exParam;
896 // take the B-field from the ESD, no 3D fieldMap available at this point
897 Bool_t relate = false;
898 relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
899 if(!relate){
900 delete track;
901 continue;
902 }
903 // fetch the track parameters at the DCA (unconstraint)
904 if(track->GetTPCInnerParam()){
905 track->GetTPCInnerParam()->GetPxPyPz(pDCA);
906 track->GetTPCInnerParam()->GetXYZ(rDCA);
907 }
908 // get the DCA to the vertex:
909 track->GetImpactParametersTPC(dDCA,cDCA);
910 // set the constrained parameters to the track
911 track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
912 }
913
914
915 Float_t eta = track->Eta();
916 Float_t pT = track->Pt();
917
918 if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){
919 delete track;
920 continue;
921 }
922
923 acctracks->Add(track);
924 }
925
926
927 return acctracks;
928
929}
930