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