Coverity fixed
[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
321     if (fUseMCRP) {
322       AliAODMCHeader *headerH = dynamic_cast<AliAODMCHeader*>(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
323       if (headerH) fRP = headerH->GetReactionPlaneAngle();
324     }
325   
326     if (aod){
327       esdEP = aod->GetHeader()->GetEventplaneP();
328       esdEP->Reset(); 
329      
330       Int_t maxID = 0;
331       TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID);
332         
333     if (fSaveTrackContribution) {
334       esdEP->GetQContributionXArray()->Set(maxID+1);
335       esdEP->GetQContributionYArray()->Set(maxID+1);
336       esdEP->GetQContributionXArraysub1()->Set(maxID+1);
337       esdEP->GetQContributionYArraysub1()->Set(maxID+1);
338       esdEP->GetQContributionXArraysub2()->Set(maxID+1);
339       esdEP->GetQContributionYArraysub2()->Set(maxID+1);
340     }
341         
342     const int NT = tracklist->GetEntries();
343       
344   if (NT>4){
345     fQVector = new TVector2(GetQ(esdEP,tracklist));
346     fEventplaneQ = fQVector->Phi()/2.; 
347     GetQsub(qq1, qq2, tracklist, esdEP);
348     fQsub1 = new TVector2(qq1);
349     fQsub2 = new TVector2(qq2);
350     fQsubRes = (fQsub1->Phi()/2. - fQsub2->Phi()/2.);
351         
352     esdEP->SetQVector(fQVector);
353     esdEP->SetEventplaneQ(fEventplaneQ);
354     esdEP->SetQsub(fQsub1,fQsub2);
355     esdEP->SetQsubRes(fQsubRes);
356         
357     fHOutEventplaneQ->Fill(fEventplaneQ);
358     fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.);
359     fHOutNTEPRes->Fill(NT,fQsubRes);
360
361         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
362         
363         for (int iter = 0; iter<NT;iter++){
364           AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
365           if (track) {
366             float delta = track->Phi()-fEventplaneQ;
367             while (delta < 0) delta += TMath::Pi();
368             while (delta > TMath::Pi()) delta -= TMath::Pi();
369             fHOutPTPsi->Fill(track->Pt(),delta);
370             fHOutPhi->Fill(track->Phi());
371             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
372           }
373         }
374         
375         AliAODTrack* trmax = aod->GetTrack(0);
376         for (int iter = 1; iter<NT;iter++){
377           AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
378           if (track && (track->Pt() > trmax->Pt())) trmax = track;
379         }
380         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
381       }     
382       delete tracklist;
383       tracklist = 0;
384     }   
385         
386     
387   }  
388
389   
390   else {
391     printf(" Analysis Input not known!\n\n ");
392     return;
393   }
394   PostData(1, fOutputList); 
395 }
396
397 //________________________________________________________________________
398 void AliEPSelectionTask::Terminate(Option_t */*option*/)
399 {
400   // Terminate analysis
401 }
402
403 //__________________________________________________________________________
404 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
405 {
406 // Get the Q vector
407   TVector2 mQ;
408   float mQx=0, mQy=0;
409   AliVTrack* track;
410   Double_t weight;
411   Int_t idtemp = -1;
412   
413   int nt = tracklist->GetEntries();
414
415   for (int i=0; i<nt; i++){
416     weight = 1;
417     track = dynamic_cast<AliVTrack*> (tracklist->At(i));
418     if (track) {
419       weight = GetWeight(track);
420     if (fSaveTrackContribution){
421       idtemp = track->GetID(); 
422       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
423       EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),idtemp);
424       EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),idtemp);
425      }
426      mQx += (weight*cos(2*track->Phi()));
427      mQy += (weight*sin(2*track->Phi()));
428     }
429   }
430   mQ.Set(mQx,mQy);
431   return mQ;
432 }
433   
434   //________________________________________________________________________
435 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP)
436 {
437 // Get Qsub
438   TVector2 mQ[2];
439   float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
440   Double_t weight;
441
442   AliVTrack* track;
443   TRandom2 rn = 0;
444   
445   int nt = tracklist->GetEntries();
446   int trackcounter1=0, trackcounter2=0;
447   int idtemp = 0;
448
449   if (fSplitMethod == AliEPSelectionTask::kRandom){
450     
451     for (Int_t i = 0; i < nt; i++) {
452       weight = 1;
453       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
454       if (!track) continue;
455       weight = GetWeight(track);
456       idtemp = track->GetID(); 
457       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
458     
459       // This loop splits the track set into 2 random subsets
460       if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
461         float random = rn.Rndm();
462         if(random < .5){
463           mQx1 += (weight*cos(2*track->Phi()));
464           mQy1 += (weight*sin(2*track->Phi()));
465           if (fSaveTrackContribution){
466             EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
467             EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
468           }
469           trackcounter1++;
470         }
471         else {
472           mQx2 += (weight*cos(2*track->Phi()));
473           mQy2 += (weight*sin(2*track->Phi()));
474           if (fSaveTrackContribution){
475             EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
476             EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
477           }
478           trackcounter2++;
479         }
480       }
481       else if( trackcounter1 >= int(nt/2.)){
482         mQx2 += (weight*cos(2*track->Phi()));
483         mQy2 += (weight*sin(2*track->Phi()));
484         if (fSaveTrackContribution){
485           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
486           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
487         }
488         trackcounter2++;
489       }
490       else {
491         mQx1 += (weight*cos(2*track->Phi()));
492         mQy1 += (weight*sin(2*track->Phi()));
493         if (fSaveTrackContribution){
494           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
495           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
496         }
497         trackcounter1++;
498       }
499     }
500   } else if (fSplitMethod == AliEPSelectionTask::kEta) {
501      
502     for (Int_t i = 0; i < nt; i++) {
503       weight = 1;
504       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
505       if (!track) continue;
506       weight = GetWeight(track);
507       Double_t eta = track->Eta();
508       idtemp = track->GetID(); 
509       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
510
511       if (eta > fEtaGap/2.) {  
512         mQx1 += (weight*cos(2*track->Phi()));
513         mQy1 += (weight*sin(2*track->Phi()));
514         if (fSaveTrackContribution){
515           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
516           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
517         }
518       } else if (eta < -1.*fEtaGap/2.) {
519         mQx2 += (weight*cos(2*track->Phi()));
520         mQy2 += (weight*sin(2*track->Phi()));
521         if (fSaveTrackContribution){
522           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
523           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
524         }
525       }
526     }
527   } else {
528     printf("plane resolution determination method not available!\n\n ");
529     return;
530   }
531      
532   mQ[0].Set(mQx1,mQy1);
533   mQ[1].Set(mQx2,mQy2);
534   Q1 = mQ[0];
535   Q2 = mQ[1];
536 }
537
538 //________________________________________________________________________
539 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
540   
541   if(fESDtrackCuts){ 
542     delete fESDtrackCuts;
543     fESDtrackCuts = 0;
544   }
545   if (fAnalysisInput.CompareTo("AOD")==0){
546     AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");  
547     fUsercuts = kFALSE;
548     SetTrackType("TPC");
549     return;
550   } 
551   fUsercuts = kTRUE;
552   fESDtrackCuts = trackcuts;
553 }
554
555 //________________________________________________________________________
556 void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup){
557   
558   if(fESDtrackCuts){ 
559     delete fESDtrackCuts;
560     fESDtrackCuts = 0;
561   }
562   if (fAnalysisInput.CompareTo("ESD")==0){
563     AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");  
564     fUsercuts = kFALSE;
565     SetTrackType("TPC");
566     return;
567   }
568   fUsercuts = kTRUE;
569   fESDtrackCuts = new AliESDtrackCuts();
570   fESDtrackCuts->SetPtRange(ptlow,ptup);
571   fESDtrackCuts->SetEtaRange(etalow,etaup);
572   fAODfilterbit = filterbit;
573 }
574
575 //_____________________________________________________________________________
576
577 void AliEPSelectionTask::SetTrackType(TString tracktype){
578   fTrackType = tracktype;
579   if (!fUsercuts) {
580   if (fTrackType.CompareTo("GLOBAL")==0){ 
581     fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
582     fAODfilterbit = 32;
583   }     
584   if (fTrackType.CompareTo("TPC")==0){  
585     fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
586     fAODfilterbit = 128;
587   }
588   fESDtrackCuts->SetPtRange(0.15,20.);
589   fESDtrackCuts->SetEtaRange(-0.8,0.8);
590   }
591 }
592
593 //________________________________________________________________________
594 Double_t AliEPSelectionTask::GetWeight(TObject* track1)
595 {
596   Double_t ptweight=1;
597   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
598   if (fUsePtWeight && track) {      
599     if (track->Pt()<2) ptweight=track->Pt();
600     else ptweight=2;
601   }
602   return ptweight*GetPhiWeight(track);
603 }
604
605 //________________________________________________________________________
606 Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
607 {
608   Double_t phiweight=1;
609   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
610   
611   if (fUsePhiWeight && fPhiDist && track) {
612     Double_t nParticles = fPhiDist->Integral();
613     Double_t nPhibins = fPhiDist->GetNbinsX();
614   
615     Double_t Phi = track->Phi();
616     
617     while (Phi<0) Phi += TMath::TwoPi();
618     while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
619       
620     Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
621     
622     if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
623   }
624   return phiweight;
625 }
626
627 //__________________________________________________________________________
628 void AliEPSelectionTask::SetPhiDist() 
629 {
630   if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
631
632     fPhiDist = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");
633     if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
634
635   } 
636   else {
637     AliInfo("Using Custom Phi Distribution");
638   }
639     
640   Bool_t emptybins;
641
642   int iter = 0;  
643   while (iter<3){
644       emptybins = kFALSE;
645    
646       for (int i=1; i<fPhiDist->GetNbinsX(); i++){
647         if (!((fPhiDist->GetBinContent(i))>0)) {
648           emptybins = kTRUE;
649         }
650       }  
651       if (emptybins) {
652         cout << "empty bins - rebinning!" << endl;
653         fPhiDist->Rebin();
654         iter++;
655       }      
656       else iter = 3;
657   }
658   
659   if (emptybins) {
660     AliError("After Maximum of rebinning still empty Phi-bins!!!");
661   }
662 }
663
664 //__________________________________________________________________________
665 void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
666 {
667   
668   fUserphidist = kTRUE;
669   
670   TFile f(infilename);
671   TObject* list = f.Get(listname);
672   fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
673   if (!fPhiDist) AliFatal("Phi Distribution not found!!!");
674
675   f.Close();
676
677
678
679 //_________________________________________________________________________
680 TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
681 {
682   TObjArray *acctracks = new TObjArray();
683   
684   AliAODTrack *tr = 0;
685   Int_t maxid1 = 0;
686   Int_t maxidtemp = -1;
687   Float_t ptlow = 0;
688   Float_t ptup = 0;
689   Float_t etalow = 0;
690   Float_t etaup = 0;
691   fESDtrackCuts->GetPtRange(ptlow,ptup);
692   fESDtrackCuts->GetEtaRange(etalow,etaup);
693   
694   for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
695      tr = aod->GetTrack(i);
696      maxidtemp = tr->GetID(); 
697      if(maxidtemp < 0 && fAODfilterbit != 128) continue;
698      if(maxidtemp > -1 && fAODfilterbit == 128) continue;
699      if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
700      if (maxidtemp > maxid1) maxid1 = maxidtemp;
701      if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow){
702      acctracks->Add(tr);
703      }
704   }
705   
706   maxid = maxid1;
707   
708   return acctracks;
709   
710 }