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