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