Fix for coverity 17562
[u/mrichter/AliRoot.git] / ANALYSIS / AliEPSelectionTask.cxx
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 <TProfile.h>
29 #include <TFile.h>
30 #include <TObjString.h>
31 #include <TString.h>
32 #include <TCanvas.h>
33 #include <TROOT.h>
34 #include <TDirectory.h>
35 #include <TSystem.h>
36 #include <iostream>
37 #include <TRandom2.h>
38 #include <TArrayF.h>
39
40 #include "AliAnalysisManager.h"
41 #include "AliVEvent.h"
42 #include "AliESD.h"
43 #include "AliESDEvent.h"
44 #include "AliMCEvent.h"
45 #include "AliESDtrack.h"
46 #include "AliESDtrackCuts.h"
47 #include "AliESDHeader.h"
48 #include "AliESDInputHandler.h"
49 #include "AliAODHandler.h"
50 #include "AliAODEvent.h"
51 #include "AliHeader.h"
52 #include "AliGenHijingEventHeader.h"
53 #include "AliAnalysisTaskSE.h"
54 #include "AliPhysicsSelectionTask.h"
55 #include "AliPhysicsSelection.h"
56 #include "AliBackgroundSelection.h"
57 #include "AliESDUtils.h"
58 #include "AliOADBContainer.h"
59 #include "AliAODMCHeader.h"
60 #include "AliAODTrack.h"
61 #include "AliVTrack.h"
62 #include "AliEventplane.h"
63
64 ClassImp(AliEPSelectionTask)
65
66 //________________________________________________________________________
67 AliEPSelectionTask::AliEPSelectionTask():
68 AliAnalysisTaskSE(),
69   fAnalysisInput("ESD"),
70   fTrackType("TPC"),
71   fUseMCRP(kFALSE),
72   fUsePhiWeight(kFALSE),
73   fUsePtWeight(kFALSE),
74   fSaveTrackContribution(kFALSE),
75   fUserphidist(kFALSE),
76   fUsercuts(kFALSE),
77   fRunNumber(-15),
78   fAODfilterbit(1),
79   fEtaGap(0.),
80   fSplitMethod(0),
81   fESDtrackCuts(0),
82   fEPContainer(0),
83   fPhiDist(0),
84   fQVector(0),
85   fQContributionX(0),
86   fQContributionY(0),
87   fEventplaneQ(0),
88   fQsub1(0),
89   fQsub2(0),
90   fQsubRes(0),  
91   fOutputList(0),
92   fHOutEventplaneQ(0),
93   fHOutPhi(0),
94   fHOutPhiCorr(0),
95   fHOutsub1sub2(0),
96   fHOutNTEPRes(0),
97   fHOutPTPsi(0),
98   fHOutDiff(0),
99   fHOutleadPTPsi(0)
100 {   
101   // Default constructor
102   AliInfo("Event Plane Selection enabled.");
103 }   
104
105 //________________________________________________________________________
106 AliEPSelectionTask::AliEPSelectionTask(const char *name):
107   AliAnalysisTaskSE(name),
108   fAnalysisInput("ESD"),
109   fTrackType("TPC"),
110   fUseMCRP(kFALSE),
111   fUsePhiWeight(kFALSE),
112   fUsePtWeight(kFALSE),  
113   fSaveTrackContribution(kFALSE),
114   fUserphidist(kFALSE),
115   fUsercuts(kFALSE),
116   fRunNumber(-15),
117   fAODfilterbit(1),
118   fEtaGap(0.),
119   fSplitMethod(0),
120   fESDtrackCuts(0),
121   fEPContainer(0),
122   fPhiDist(0),
123   fQVector(0),
124   fQContributionX(0),
125   fQContributionY(0),
126   fEventplaneQ(0),
127   fQsub1(0),
128   fQsub2(0),
129   fQsubRes(0),
130   fOutputList(0),
131   fHOutEventplaneQ(0),
132   fHOutPhi(0),
133   fHOutPhiCorr(0),
134   fHOutsub1sub2(0),
135   fHOutNTEPRes(0),
136   fHOutPTPsi(0),
137   fHOutDiff(0),
138   fHOutleadPTPsi(0)
139 {
140   // Default constructor
141   AliInfo("Event Plane Selection enabled.");
142   DefineOutput(1, TList::Class());
143 }
144  
145 //________________________________________________________________________
146 AliEPSelectionTask::~AliEPSelectionTask()
147 {
148   // Destructor  
149   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
150       delete fOutputList;
151       fOutputList = 0;
152   }
153   if (fESDtrackCuts){
154       delete fESDtrackCuts;
155       fESDtrackCuts = 0;
156   }
157   if (fUserphidist) {
158     if (fPhiDist) {
159       delete fPhiDist;
160       fPhiDist = 0;
161     }
162   }
163   if (fEPContainer){
164       delete fEPContainer;
165       fEPContainer = 0;
166   }
167 }  
168
169 //________________________________________________________________________
170 void AliEPSelectionTask::UserCreateOutputObjects()
171 {  
172   // Create the output containers
173   if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
174   AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
175
176   fOutputList = new TList();
177   fOutputList->SetOwner();
178   fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
179   fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
180   fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
181   fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
182   fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
183   fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
184   fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
185   fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
186
187   fOutputList->Add(fHOutEventplaneQ);
188   fOutputList->Add(fHOutPhi);
189   fOutputList->Add(fHOutPhiCorr);
190   fOutputList->Add(fHOutsub1sub2);
191   fOutputList->Add(fHOutNTEPRes);
192   fOutputList->Add(fHOutPTPsi);
193   fOutputList->Add(fHOutleadPTPsi);
194   fOutputList->Add(fHOutDiff);
195   
196   PostData(1, fOutputList); 
197   
198     if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
199
200     
201     TString oadbfilename; 
202
203     if (fAnalysisInput.CompareTo("AOD")==0){
204       oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
205     } else if (fAnalysisInput.CompareTo("ESD")==0){
206       oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
207     }
208  
209     TFile foadb(oadbfilename); 
210     if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
211
212     AliInfo("Using Standard OADB");
213     fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");    
214     if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
215     foadb.Close();
216     }
217
218 }
219
220 //________________________________________________________________________
221 void AliEPSelectionTask::UserExec(Option_t */*option*/)
222
223   // Execute analysis for current event:
224   if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
225   
226 //   fRunNumber = -15;
227  
228   AliEventplane *esdEP;
229   TVector2 qq1;
230   TVector2 qq2;
231   Double_t fRP = 0.; // the monte carlo reaction plane angle
232     
233   if (fAnalysisInput.CompareTo("ESD")==0){
234
235     AliVEvent* event = InputEvent();
236     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
237     if (esd){    
238       if (!(fRunNumber == esd->GetRunNumber())) {
239           fRunNumber = esd->GetRunNumber();
240             SetPhiDist();      
241       }
242       
243       
244       if (fUseMCRP) {
245           AliMCEvent* mcEvent  = MCEvent();      
246           if (mcEvent && mcEvent->GenEventHeader()) {
247               AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
248               if (headerH) fRP = headerH->ReactionPlaneAngle();
249           }
250       }
251       
252       esdEP = esd->GetEventplane();
253       if (fSaveTrackContribution) {
254         esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
255         esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
256         esdEP->GetQContributionXArraysub1()->Set(esd->GetNumberOfTracks());
257         esdEP->GetQContributionYArraysub1()->Set(esd->GetNumberOfTracks());
258         esdEP->GetQContributionXArraysub2()->Set(esd->GetNumberOfTracks());
259         esdEP->GetQContributionYArraysub2()->Set(esd->GetNumberOfTracks());
260       }
261       
262       TObjArray* tracklist = new TObjArray;
263       if (fTrackType.CompareTo("GLOBAL")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
264       if (fTrackType.CompareTo("TPC")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
265       const int nt = tracklist->GetEntries();
266       
267       if (nt>4){
268         fQVector = new TVector2(GetQ(esdEP,tracklist));
269         fEventplaneQ = fQVector->Phi()/2; 
270         GetQsub(qq1, qq2, tracklist, esdEP);
271         fQsub1 = new TVector2(qq1);
272         fQsub2 = new TVector2(qq2);
273         fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
274         
275         esdEP->SetQVector(fQVector);
276         esdEP->SetEventplaneQ(fEventplaneQ);
277         esdEP->SetQsub(fQsub1,fQsub2);
278         esdEP->SetQsubRes(fQsubRes);
279         
280         fHOutEventplaneQ->Fill(fEventplaneQ);
281         fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
282         fHOutNTEPRes->Fill(nt,fQsubRes);
283
284         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
285         
286         for (int iter = 0; iter<nt;iter++){
287           AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
288           if (track) {
289             float delta = track->Phi()-fEventplaneQ;
290             while (delta < 0) delta += TMath::Pi();
291             while (delta > TMath::Pi()) delta -= TMath::Pi();
292             fHOutPTPsi->Fill(track->Pt(),delta);
293             fHOutPhi->Fill(track->Phi());
294             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
295           }
296         }
297         
298         AliESDtrack* trmax = esd->GetTrack(0);
299         for (int iter = 1; iter<nt;iter++){
300           AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
301           if (track && (track->Pt() > trmax->Pt())) trmax = track;
302         }
303         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
304       }
305       delete tracklist;
306       tracklist = 0;
307     }
308   }
309   
310   else if (fAnalysisInput.CompareTo("AOD")==0){
311     AliVEvent* event = InputEvent();
312     AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
313
314     if (aod){
315       if (!(fRunNumber == aod->GetRunNumber())) {
316         fRunNumber = aod->GetRunNumber();
317         SetPhiDist();      
318       }
319   
320       if (fUseMCRP) {
321         AliAODMCHeader *headerH = dynamic_cast<AliAODMCHeader*>(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
322         if (headerH) fRP = headerH->GetReactionPlaneAngle();
323       }
324   
325       esdEP = aod->GetHeader()->GetEventplaneP();
326       esdEP->Reset(); 
327      
328       Int_t maxID = 0;
329       TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID);
330         
331       if (fSaveTrackContribution) {
332         esdEP->GetQContributionXArray()->Set(maxID+1);
333         esdEP->GetQContributionYArray()->Set(maxID+1);
334         esdEP->GetQContributionXArraysub1()->Set(maxID+1);
335         esdEP->GetQContributionYArraysub1()->Set(maxID+1);
336         esdEP->GetQContributionXArraysub2()->Set(maxID+1);
337         esdEP->GetQContributionYArraysub2()->Set(maxID+1);
338       }
339         
340       const int NT = tracklist->GetEntries();
341       
342       if (NT>4){
343         fQVector = new TVector2(GetQ(esdEP,tracklist));
344         fEventplaneQ = fQVector->Phi()/2.; 
345         GetQsub(qq1, qq2, tracklist, esdEP);
346         fQsub1 = new TVector2(qq1);
347         fQsub2 = new TVector2(qq2);
348         fQsubRes = (fQsub1->Phi()/2. - fQsub2->Phi()/2.);
349         
350         esdEP->SetQVector(fQVector);
351         esdEP->SetEventplaneQ(fEventplaneQ);
352         esdEP->SetQsub(fQsub1,fQsub2);
353         esdEP->SetQsubRes(fQsubRes);
354         
355         fHOutEventplaneQ->Fill(fEventplaneQ);
356         fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.);
357         fHOutNTEPRes->Fill(NT,fQsubRes);
358
359         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
360         
361         for (int iter = 0; iter<NT;iter++){
362           AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
363           if (track) {
364             float delta = track->Phi()-fEventplaneQ;
365             while (delta < 0) delta += TMath::Pi();
366             while (delta > TMath::Pi()) delta -= TMath::Pi();
367             fHOutPTPsi->Fill(track->Pt(),delta);
368             fHOutPhi->Fill(track->Phi());
369             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
370           }
371         }
372         
373         AliAODTrack* trmax = aod->GetTrack(0);
374         for (int iter = 1; iter<NT;iter++){
375           AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
376           if (track && (track->Pt() > trmax->Pt())) trmax = track;
377         }
378         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
379       }     
380       delete tracklist;
381       tracklist = 0;
382     }   
383         
384     
385   }  
386
387   
388   else {
389     printf(" Analysis Input not known!\n\n ");
390     return;
391   }
392   PostData(1, fOutputList); 
393 }
394
395 //________________________________________________________________________
396 void AliEPSelectionTask::Terminate(Option_t */*option*/)
397 {
398   // Terminate analysis
399 }
400
401 //__________________________________________________________________________
402 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
403 {
404 // Get the Q vector
405   TVector2 mQ;
406   float mQx=0, mQy=0;
407   AliVTrack* track;
408   Double_t weight;
409   Int_t idtemp = -1;
410   
411   int nt = tracklist->GetEntries();
412
413   for (int i=0; i<nt; i++){
414     weight = 1;
415     track = dynamic_cast<AliVTrack*> (tracklist->At(i));
416     if (track) {
417       weight = GetWeight(track);
418     if (fSaveTrackContribution){
419       idtemp = track->GetID(); 
420       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
421       EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),idtemp);
422       EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),idtemp);
423      }
424      mQx += (weight*cos(2*track->Phi()));
425      mQy += (weight*sin(2*track->Phi()));
426     }
427   }
428   mQ.Set(mQx,mQy);
429   return mQ;
430 }
431   
432   //________________________________________________________________________
433 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP)
434 {
435 // Get Qsub
436   TVector2 mQ[2];
437   float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
438   Double_t weight;
439
440   AliVTrack* track;
441   TRandom2 rn = 0;
442   
443   int nt = tracklist->GetEntries();
444   int trackcounter1=0, trackcounter2=0;
445   int idtemp = 0;
446
447   if (fSplitMethod == AliEPSelectionTask::kRandom){
448     
449     for (Int_t i = 0; i < nt; i++) {
450       weight = 1;
451       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
452       if (!track) continue;
453       weight = GetWeight(track);
454       idtemp = track->GetID(); 
455       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
456     
457       // This loop splits the track set into 2 random subsets
458       if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
459         float random = rn.Rndm();
460         if(random < .5){
461           mQx1 += (weight*cos(2*track->Phi()));
462           mQy1 += (weight*sin(2*track->Phi()));
463           if (fSaveTrackContribution){
464             EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
465             EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
466           }
467           trackcounter1++;
468         }
469         else {
470           mQx2 += (weight*cos(2*track->Phi()));
471           mQy2 += (weight*sin(2*track->Phi()));
472           if (fSaveTrackContribution){
473             EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
474             EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
475           }
476           trackcounter2++;
477         }
478       }
479       else if( trackcounter1 >= int(nt/2.)){
480         mQx2 += (weight*cos(2*track->Phi()));
481         mQy2 += (weight*sin(2*track->Phi()));
482         if (fSaveTrackContribution){
483           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
484           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
485         }
486         trackcounter2++;
487       }
488       else {
489         mQx1 += (weight*cos(2*track->Phi()));
490         mQy1 += (weight*sin(2*track->Phi()));
491         if (fSaveTrackContribution){
492           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
493           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
494         }
495         trackcounter1++;
496       }
497     }
498   } else if (fSplitMethod == AliEPSelectionTask::kEta) {
499      
500     for (Int_t i = 0; i < nt; i++) {
501       weight = 1;
502       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
503       if (!track) continue;
504       weight = GetWeight(track);
505       Double_t eta = track->Eta();
506       idtemp = track->GetID(); 
507       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
508
509       if (eta > fEtaGap/2.) {  
510         mQx1 += (weight*cos(2*track->Phi()));
511         mQy1 += (weight*sin(2*track->Phi()));
512         if (fSaveTrackContribution){
513           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
514           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
515         }
516       } else if (eta < -1.*fEtaGap/2.) {
517         mQx2 += (weight*cos(2*track->Phi()));
518         mQy2 += (weight*sin(2*track->Phi()));
519         if (fSaveTrackContribution){
520           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
521           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
522         }
523       }
524     }
525   } else {
526     printf("plane resolution determination method not available!\n\n ");
527     return;
528   }
529      
530   mQ[0].Set(mQx1,mQy1);
531   mQ[1].Set(mQx2,mQy2);
532   Q1 = mQ[0];
533   Q2 = mQ[1];
534 }
535
536 //________________________________________________________________________
537 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
538   
539   if(fESDtrackCuts){ 
540     delete fESDtrackCuts;
541     fESDtrackCuts = 0;
542   }
543   if (fAnalysisInput.CompareTo("AOD")==0){
544     AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");  
545     fUsercuts = kFALSE;
546     SetTrackType("TPC");
547     return;
548   } 
549   fUsercuts = kTRUE;
550   fESDtrackCuts = trackcuts;
551 }
552
553 //________________________________________________________________________
554 void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup){
555   
556   if(fESDtrackCuts){ 
557     delete fESDtrackCuts;
558     fESDtrackCuts = 0;
559   }
560   if (fAnalysisInput.CompareTo("ESD")==0){
561     AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");  
562     fUsercuts = kFALSE;
563     SetTrackType("TPC");
564     return;
565   }
566   fUsercuts = kTRUE;
567   fESDtrackCuts = new AliESDtrackCuts();
568   fESDtrackCuts->SetPtRange(ptlow,ptup);
569   fESDtrackCuts->SetEtaRange(etalow,etaup);
570   fAODfilterbit = filterbit;
571 }
572
573 //_____________________________________________________________________________
574
575 void AliEPSelectionTask::SetTrackType(TString tracktype){
576   fTrackType = tracktype;
577   if (!fUsercuts) {
578   if (fTrackType.CompareTo("GLOBAL")==0){ 
579     fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
580     fAODfilterbit = 32;
581   }     
582   if (fTrackType.CompareTo("TPC")==0){  
583     fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
584     fAODfilterbit = 128;
585   }
586   fESDtrackCuts->SetPtRange(0.15,20.);
587   fESDtrackCuts->SetEtaRange(-0.8,0.8);
588   }
589 }
590
591 //________________________________________________________________________
592 Double_t AliEPSelectionTask::GetWeight(TObject* track1)
593 {
594   Double_t ptweight=1;
595   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
596   if (fUsePtWeight && track) {      
597     if (track->Pt()<2) ptweight=track->Pt();
598     else ptweight=2;
599   }
600   return ptweight*GetPhiWeight(track);
601 }
602
603 //________________________________________________________________________
604 Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
605 {
606   Double_t phiweight=1;
607   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
608   
609   if (fUsePhiWeight && fPhiDist && track) {
610     Double_t nParticles = fPhiDist->Integral();
611     Double_t nPhibins = fPhiDist->GetNbinsX();
612   
613     Double_t Phi = track->Phi();
614     
615     while (Phi<0) Phi += TMath::TwoPi();
616     while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
617       
618     Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
619     
620     if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
621   }
622   return phiweight;
623 }
624
625 //__________________________________________________________________________
626 void AliEPSelectionTask::SetPhiDist() 
627 {
628   if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
629
630     fPhiDist = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");
631     if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
632
633   } 
634   else {
635     AliInfo("Using Custom Phi Distribution");
636   }
637     
638   Bool_t emptybins;
639
640   int iter = 0;  
641   while (iter<3){
642       emptybins = kFALSE;
643    
644       for (int i=1; i<fPhiDist->GetNbinsX(); i++){
645         if (!((fPhiDist->GetBinContent(i))>0)) {
646           emptybins = kTRUE;
647         }
648       }  
649       if (emptybins) {
650         cout << "empty bins - rebinning!" << endl;
651         fPhiDist->Rebin();
652         iter++;
653       }      
654       else iter = 3;
655   }
656   
657   if (emptybins) {
658     AliError("After Maximum of rebinning still empty Phi-bins!!!");
659   }
660 }
661
662 //__________________________________________________________________________
663 void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
664 {
665   
666   fUserphidist = kTRUE;
667   
668   TFile f(infilename);
669   TObject* list = f.Get(listname);
670   fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
671   if (!fPhiDist) AliFatal("Phi Distribution not found!!!");
672
673   f.Close();
674
675
676
677 //_________________________________________________________________________
678 TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
679 {
680   TObjArray *acctracks = new TObjArray();
681   
682   AliAODTrack *tr = 0;
683   Int_t maxid1 = 0;
684   Int_t maxidtemp = -1;
685   Float_t ptlow = 0;
686   Float_t ptup = 0;
687   Float_t etalow = 0;
688   Float_t etaup = 0;
689   fESDtrackCuts->GetPtRange(ptlow,ptup);
690   fESDtrackCuts->GetEtaRange(etalow,etaup);
691   
692   for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
693      tr = aod->GetTrack(i);
694      maxidtemp = tr->GetID(); 
695      if(maxidtemp < 0 && fAODfilterbit != 128) continue;
696      if(maxidtemp > -1 && fAODfilterbit == 128) continue;
697      if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
698      if (maxidtemp > maxid1) maxid1 = maxidtemp;
699      if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow){
700      acctracks->Add(tr);
701      }
702   }
703   
704   maxid = maxid1;
705   
706   return acctracks;
707   
708 }