write only the complement in the case of hybrid tracks
[u/mrichter/AliRoot.git] / ANALYSIS / AliEPSelectionTask.cxx
index 1830b51..2d6007b 100644 (file)
 #include "AliPhysicsSelection.h"
 #include "AliBackgroundSelection.h"
 #include "AliESDUtils.h"
-
+#include "AliOADBContainer.h"
+#include "AliAODMCHeader.h"
+#include "AliAODTrack.h"
+#include "AliVTrack.h"
 #include "AliEventplane.h"
 
 ClassImp(AliEPSelectionTask)
@@ -63,15 +66,20 @@ ClassImp(AliEPSelectionTask)
 //________________________________________________________________________
 AliEPSelectionTask::AliEPSelectionTask():
 AliAnalysisTaskSE(),
-  fDebug(0),
   fAnalysisInput("ESD"),
-  fStatus("TPC"),
+  fTrackType("TPC"),
   fUseMCRP(kFALSE),
   fUsePhiWeight(kFALSE),
   fUsePtWeight(kFALSE),
   fSaveTrackContribution(kFALSE),
+  fUserphidist(kFALSE),
+  fUsercuts(kFALSE),
+  fRunNumber(-15),
+  fAODfilterbit(1),
+  fEtaGap(0.),
+  fSplitMethod(0),
   fESDtrackCuts(0),
-  ftracklist(0),
+  fEPContainer(0),
   fPhiDist(0),
   fQVector(0),
   fQContributionX(0),
@@ -97,15 +105,20 @@ AliAnalysisTaskSE(),
 //________________________________________________________________________
 AliEPSelectionTask::AliEPSelectionTask(const char *name):
   AliAnalysisTaskSE(name),
-  fDebug(0),
   fAnalysisInput("ESD"),
-  fStatus("TPC"),
+  fTrackType("TPC"),
   fUseMCRP(kFALSE),
   fUsePhiWeight(kFALSE),
   fUsePtWeight(kFALSE),  
   fSaveTrackContribution(kFALSE),
+  fUserphidist(kFALSE),
+  fUsercuts(kFALSE),
+  fRunNumber(-15),
+  fAODfilterbit(1),
+  fEtaGap(0.),
+  fSplitMethod(0),
   fESDtrackCuts(0),
-  ftracklist(0),
+  fEPContainer(0),
   fPhiDist(0),
   fQVector(0),
   fQContributionX(0),
@@ -141,17 +154,15 @@ AliEPSelectionTask::~AliEPSelectionTask()
       delete fESDtrackCuts;
       fESDtrackCuts = 0;
   }
-  if (fQVector){
-      delete fQVector;
-      fQVector = 0;
-  }
-  if (fQsub1){
-      delete fQsub1;
-      fQsub1 = 0;
+  if (fUserphidist) {
+    if (fPhiDist) {
+      delete fPhiDist;
+      fPhiDist = 0;
+    }
   }
-    if (fQsub2){
-      delete fQsub2;
-      fQsub2 = 0;
+  if (fEPContainer){
+      delete fEPContainer;
+      fEPContainer = 0;
   }
 }  
 
@@ -183,6 +194,27 @@ void AliEPSelectionTask::UserCreateOutputObjects()
   fOutputList->Add(fHOutDiff);
   
   PostData(1, fOutputList); 
+  
+    if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
+
+    
+    TString oadbfilename; 
+
+    if (fAnalysisInput.CompareTo("AOD")==0){
+      oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
+    } else if (fAnalysisInput.CompareTo("ESD")==0){
+      oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
+    }
+    TFile foadb(oadbfilename); 
+    if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
+
+    AliInfo("Using Standard OADB");
+    fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");    
+    if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
+    foadb.Close();
+    }
+
 }
 
 //________________________________________________________________________
@@ -190,43 +222,56 @@ void AliEPSelectionTask::UserExec(Option_t */*option*/)
 { 
   // Execute analysis for current event:
   if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
+  
+//   fRunNumber = -15;
  
-  AliEventplane* esdEP = 0;
-  TVector2 QQ1;
-  TVector2 QQ2;
+  AliEventplane *esdEP;
+  TVector2 qq1;
+  TVector2 qq2;
   Double_t fRP = 0.; // the monte carlo reaction plane angle
     
   if (fAnalysisInput.CompareTo("ESD")==0){
-    
+
     AliVEvent* event = InputEvent();
     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
-    
-    if (fUseMCRP) {
-      AliMCEvent* mcEvent  = MCEvent();      
-      if (mcEvent && mcEvent->GenEventHeader()) {
-       AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
-       fRP = headerH->ReactionPlaneAngle();
+    if (esd){    
+      if (!(fRunNumber == esd->GetRunNumber())) {
+         fRunNumber = esd->GetRunNumber();
+           SetPhiDist();      
       }
-    }
-    if (esd){
+      
+      
+      if (fUseMCRP) {
+         AliMCEvent* mcEvent  = MCEvent();      
+         if (mcEvent && mcEvent->GenEventHeader()) {
+             AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
+             if (headerH) fRP = headerH->ReactionPlaneAngle();
+         }
+      }
+      
       esdEP = esd->GetEventplane();
       if (fSaveTrackContribution) {
        esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
        esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
+        esdEP->GetQContributionXArraysub1()->Set(esd->GetNumberOfTracks());
+       esdEP->GetQContributionYArraysub1()->Set(esd->GetNumberOfTracks());
+        esdEP->GetQContributionXArraysub2()->Set(esd->GetNumberOfTracks());
+       esdEP->GetQContributionYArraysub2()->Set(esd->GetNumberOfTracks());
       }
       
-      if (fStatus.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
-      if (fStatus.CompareTo("TPC")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
-      const int NT = ftracklist->GetEntries();
+      TObjArray* tracklist = new TObjArray;
+      if (fTrackType.CompareTo("GLOBAL")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
+      if (fTrackType.CompareTo("TPC")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
+      const int nt = tracklist->GetEntries();
       
-      if (NT>4){
-       fQVector = new TVector2(GetQ(esdEP,ftracklist));
+      if (nt>4){
+       fQVector = new TVector2(GetQ(esdEP,tracklist));
        fEventplaneQ = fQVector->Phi()/2; 
-       GetQsub(QQ1, QQ2, ftracklist);
-       fQsub1 = new TVector2(QQ1);
-       fQsub2 = new TVector2(QQ2);
+       GetQsub(qq1, qq2, tracklist, esdEP);
+       fQsub1 = new TVector2(qq1);
+       fQsub2 = new TVector2(qq2);
        fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
-  
+       
        esdEP->SetQVector(fQVector);
        esdEP->SetEventplaneQ(fEventplaneQ);
        esdEP->SetQsub(fQsub1,fQsub2);
@@ -234,37 +279,112 @@ void AliEPSelectionTask::UserExec(Option_t */*option*/)
        
        fHOutEventplaneQ->Fill(fEventplaneQ);
        fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
-       fHOutNTEPRes->Fill(NT,fQsubRes);
+       fHOutNTEPRes->Fill(nt,fQsubRes);
 
        if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
        
-       for (int iter = 0; iter<NT;iter++){
-         AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
-         float delta = track->Phi()-fEventplaneQ;
-         while (delta < 0) delta += TMath::Pi();
-         while (delta > TMath::Pi()) delta -= TMath::Pi();
-         fHOutPTPsi->Fill(track->Pt(),delta);
-         fHOutPhi->Fill(track->Phi());
-         fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
+       for (int iter = 0; iter<nt;iter++){
+         AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
+         if (track) {
+           float delta = track->Phi()-fEventplaneQ;
+           while (delta < 0) delta += TMath::Pi();
+           while (delta > TMath::Pi()) delta -= TMath::Pi();
+           fHOutPTPsi->Fill(track->Pt(),delta);
+           fHOutPhi->Fill(track->Phi());
+           fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
+         }
        }
        
        AliESDtrack* trmax = esd->GetTrack(0);
-       for (int iter = 1; iter<NT;iter++){
-         AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
-         if (track->Pt() > trmax->Pt()) trmax = track;
+       for (int iter = 1; iter<nt;iter++){
+         AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
+         if (track && (track->Pt() > trmax->Pt())) trmax = track;
        }
        fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
       }
-      delete ftracklist;
+      delete tracklist;
+      tracklist = 0;
     }
   }
   
-  else if (fAnalysisInput.CompareTo("AOD")==0){
-    //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
-    // to be implemented
-    printf("  AOD analysis not yet implemented!!!\n\n");
-    return;
+    else if (fAnalysisInput.CompareTo("AOD")==0){
+    AliVEvent* event = InputEvent();
+    AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
+
+    if (aod){
+      if (!(fRunNumber == aod->GetRunNumber())) {
+        fRunNumber = aod->GetRunNumber();
+        SetPhiDist();      
+      }
+
+      if (fUseMCRP) {
+       AliAODMCHeader *headerH = dynamic_cast<AliAODMCHeader*>(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+       if (headerH) fRP = headerH->GetReactionPlaneAngle();
+      }
+  
+      esdEP = aod->GetHeader()->GetEventplaneP();
+      esdEP->Reset(); 
+     
+      Int_t maxID = 0;
+      TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID);
+       
+      if (fSaveTrackContribution) {
+       esdEP->GetQContributionXArray()->Set(maxID+1);
+       esdEP->GetQContributionYArray()->Set(maxID+1);
+       esdEP->GetQContributionXArraysub1()->Set(maxID+1);
+       esdEP->GetQContributionYArraysub1()->Set(maxID+1);
+       esdEP->GetQContributionXArraysub2()->Set(maxID+1);
+       esdEP->GetQContributionYArraysub2()->Set(maxID+1);
+      }
+       
+      const int NT = tracklist->GetEntries();
+      
+      if (NT>4){
+       fQVector = new TVector2(GetQ(esdEP,tracklist));
+       fEventplaneQ = fQVector->Phi()/2.; 
+       GetQsub(qq1, qq2, tracklist, esdEP);
+       fQsub1 = new TVector2(qq1);
+       fQsub2 = new TVector2(qq2);
+       fQsubRes = (fQsub1->Phi()/2. - fQsub2->Phi()/2.);
+       
+       esdEP->SetQVector(fQVector);
+       esdEP->SetEventplaneQ(fEventplaneQ);
+       esdEP->SetQsub(fQsub1,fQsub2);
+       esdEP->SetQsubRes(fQsubRes);
+       
+       fHOutEventplaneQ->Fill(fEventplaneQ);
+       fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.);
+       fHOutNTEPRes->Fill(NT,fQsubRes);
+       
+       if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
+       
+       for (int iter = 0; iter<NT;iter++){
+         AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
+         if (track) {
+           float delta = track->Phi()-fEventplaneQ;
+           while (delta < 0) delta += TMath::Pi();
+           while (delta > TMath::Pi()) delta -= TMath::Pi();
+           fHOutPTPsi->Fill(track->Pt(),delta);
+           fHOutPhi->Fill(track->Phi());
+           fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
+         }
+       }
+       
+       AliAODTrack* trmax = aod->GetTrack(0);
+       for (int iter = 1; iter<NT;iter++){
+         AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
+         if (track && (track->Pt() > trmax->Pt())) trmax = track;
+       }
+       fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
+      }     
+      delete tracklist;
+      tracklist = 0;
+    }  
+       
+    
   }  
+
+  
   else {
     printf(" Analysis Input not known!\n\n ");
     return;
@@ -281,71 +401,132 @@ void AliEPSelectionTask::Terminate(Option_t */*option*/)
 //__________________________________________________________________________
 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
 {
+// Get the Q vector
   TVector2 mQ;
   float mQx=0, mQy=0;
-  AliESDtrack* track;
+  AliVTrack* track;
   Double_t weight;
+  Int_t idtemp = -1;
   
-  int NT = tracklist->GetEntries();
+  int nt = tracklist->GetEntries();
 
-  for (int i=0; i<NT; i++){
+  for (int i=0; i<nt; i++){
     weight = 1;
-    track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
-    weight = GetWeight(track);
+    track = dynamic_cast<AliVTrack*> (tracklist->At(i));
+    if (track) {
+      weight = GetWeight(track);
     if (fSaveTrackContribution){
-      EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),track->GetID());
-      EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),track->GetID());
+      idtemp = track->GetID(); 
+      if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
+      EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),idtemp);
+      EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),idtemp);
+     }
+     mQx += (weight*cos(2*track->Phi()));
+     mQy += (weight*sin(2*track->Phi()));
     }
-    mQx += (weight*cos(2*track->Phi()));
-    mQy += (weight*sin(2*track->Phi()));
   }
   mQ.Set(mQx,mQy);
   return mQ;
 }
   
   //________________________________________________________________________
-void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
+void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP)
 {
+// Get Qsub
   TVector2 mQ[2];
   float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
   Double_t weight;
 
-  AliESDtrack* track;
-  TRandom2 RN = 0;
+  AliVTrack* track;
+  TRandom2 rn = 0;
   
-  int NT = tracklist->GetEntries();
+  int nt = tracklist->GetEntries();
   int trackcounter1=0, trackcounter2=0;
-  
-  for (Int_t i = 0; i < NT; i++) {
-    weight = 1;
-    track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
-    weight = GetWeight(track);
+  int idtemp = 0;
+
+  if (fSplitMethod == AliEPSelectionTask::kRandom){
+    
+    for (Int_t i = 0; i < nt; i++) {
+      weight = 1;
+      track = dynamic_cast<AliVTrack*> (tracklist->At(i));
+      if (!track) continue;
+      weight = GetWeight(track);
+      idtemp = track->GetID(); 
+      if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
     
-    // This loop splits the track set into 2 random subsets
-    if( trackcounter1 < int(NT/2.) && trackcounter2 < int(NT/2.)){
-      float random = RN.Rndm();
-      if(random < .5){
+      // This loop splits the track set into 2 random subsets
+      if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
+        float random = rn.Rndm();
+        if(random < .5){
+          mQx1 += (weight*cos(2*track->Phi()));
+          mQy1 += (weight*sin(2*track->Phi()));
+          if (fSaveTrackContribution){
+            EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
+            EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
+          }
+          trackcounter1++;
+        }
+        else {
+          mQx2 += (weight*cos(2*track->Phi()));
+          mQy2 += (weight*sin(2*track->Phi()));
+          if (fSaveTrackContribution){
+            EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
+            EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
+          }
+          trackcounter2++;
+        }
+      }
+      else if( trackcounter1 >= int(nt/2.)){
+        mQx2 += (weight*cos(2*track->Phi()));
+        mQy2 += (weight*sin(2*track->Phi()));
+        if (fSaveTrackContribution){
+          EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
+          EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
+        }
+        trackcounter2++;
+      }
+      else {
         mQx1 += (weight*cos(2*track->Phi()));
         mQy1 += (weight*sin(2*track->Phi()));
+        if (fSaveTrackContribution){
+          EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
+          EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
+        }
         trackcounter1++;
       }
-      else {
+    }
+  } else if (fSplitMethod == AliEPSelectionTask::kEta) {
+     
+    for (Int_t i = 0; i < nt; i++) {
+      weight = 1;
+      track = dynamic_cast<AliVTrack*> (tracklist->At(i));
+      if (!track) continue;
+      weight = GetWeight(track);
+      Double_t eta = track->Eta();
+      idtemp = track->GetID(); 
+      if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
+
+      if (eta > fEtaGap/2.) {  
+        mQx1 += (weight*cos(2*track->Phi()));
+        mQy1 += (weight*sin(2*track->Phi()));
+        if (fSaveTrackContribution){
+          EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
+          EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
+        }
+      } else if (eta < -1.*fEtaGap/2.) {
         mQx2 += (weight*cos(2*track->Phi()));
         mQy2 += (weight*sin(2*track->Phi()));
-        trackcounter2++;
+        if (fSaveTrackContribution){
+          EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
+          EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
+        }
       }
     }
-    else if( trackcounter1 >= int(NT/2.)){
-      mQx2 += (weight*cos(2*track->Phi()));
-      mQy2 += (weight*sin(2*track->Phi()));
-      trackcounter2++;
-    }
-    else {
-      mQx1 += (weight*cos(2*track->Phi()));
-      mQy1 += (weight*sin(2*track->Phi()));
-      trackcounter1++;
-    }
+  } else {
+    printf("plane resolution determination method not available!\n\n ");
+    return;
   }
+     
   mQ[0].Set(mQx1,mQy1);
   mQ[1].Set(mQx2,mQy2);
   Q1 = mQ[0];
@@ -353,22 +534,107 @@ void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklis
 }
 
 //________________________________________________________________________
-void AliEPSelectionTask::SetESDtrackCuts(TString status){
+void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
   
-  fStatus = status;
-  if (fStatus.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
-  if (fStatus.CompareTo("TPC")==0)    fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
-  fESDtrackCuts->SetPtRange(0.15,20);
+  if(fESDtrackCuts){ 
+    delete fESDtrackCuts;
+    fESDtrackCuts = 0;
+  }
+  if (fAnalysisInput.CompareTo("AOD")==0){
+    AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");  
+    fUsercuts = kFALSE;
+    SetTrackType("TPC");
+    return;
+  } 
+  fUsercuts = kTRUE;
+  fESDtrackCuts = trackcuts;
+}
+
+//________________________________________________________________________
+void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup){
+  
+  if(fESDtrackCuts){ 
+    delete fESDtrackCuts;
+    fESDtrackCuts = 0;
+  }
+  if (fAnalysisInput.CompareTo("ESD")==0){
+    AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");  
+    fUsercuts = kFALSE;
+    SetTrackType("TPC");
+    return;
+  }
+  fUsercuts = kTRUE;
+  fESDtrackCuts = new AliESDtrackCuts();
+  fESDtrackCuts->SetPtRange(ptlow,ptup);
+  fESDtrackCuts->SetEtaRange(etalow,etaup);
+  fAODfilterbit = filterbit;
+}
+
+//_____________________________________________________________________________
+
+void AliEPSelectionTask::SetTrackType(TString tracktype){
+  fTrackType = tracktype;
+  if (!fUsercuts) {
+  if (fTrackType.CompareTo("GLOBAL")==0){ 
+    fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
+    fAODfilterbit = 32;
+  }    
+  if (fTrackType.CompareTo("TPC")==0){  
+    fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
+    fAODfilterbit = 128;
+  }
+  fESDtrackCuts->SetPtRange(0.15,20.);
   fESDtrackCuts->SetEtaRange(-0.8,0.8);
+  }
 }
 
-//__________________________________________________________________________
-void AliEPSelectionTask::SetPhiDistribution(char* infilename, char* listname){
-  TFile f(infilename);
-  TObject* list = f.Get(listname);
-  fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
-  if (!fPhiDist) cout << "Phi Distribution not found!!!" << endl;
+//________________________________________________________________________
+Double_t AliEPSelectionTask::GetWeight(TObject* track1)
+{
+  Double_t ptweight=1;
+  AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
+  if (fUsePtWeight && track) {      
+    if (track->Pt()<2) ptweight=track->Pt();
+    else ptweight=2;
+  }
+  return ptweight*GetPhiWeight(track);
+}
+
+//________________________________________________________________________
+Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
+{
+  Double_t phiweight=1;
+  AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
+  
+  if (fUsePhiWeight && fPhiDist && track) {
+    Double_t nParticles = fPhiDist->Integral();
+    Double_t nPhibins = fPhiDist->GetNbinsX();
   
+    Double_t Phi = track->Phi();
+    
+    while (Phi<0) Phi += TMath::TwoPi();
+    while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
+      
+    Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
+    
+    if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
+  }
+  return phiweight;
+}
+
+//__________________________________________________________________________
+void AliEPSelectionTask::SetPhiDist() 
+{
+  if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
+
+    fPhiDist = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");
+    if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
+
+  } 
+  else {
+    AliInfo("Using Custom Phi Distribution");
+  }
+    
   Bool_t emptybins;
 
   int iter = 0;  
@@ -391,38 +657,52 @@ void AliEPSelectionTask::SetPhiDistribution(char* infilename, char* listname){
   if (emptybins) {
     AliError("After Maximum of rebinning still empty Phi-bins!!!");
   }
-  f.Close();
 }
 
-//________________________________________________________________________
-Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
+//__________________________________________________________________________
+void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
 {
-  Double_t ptweight=1;
+  
+  fUserphidist = kTRUE;
+  
+  TFile f(infilename);
+  TObject* list = f.Get(listname);
+  fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
+  if (!fPhiDist) AliFatal("Phi Distribution not found!!!");
 
-  if (fUsePtWeight) {      
-    if (track->Pt()<2) ptweight=track->Pt();
-    else ptweight=2;
-  }
-  return ptweight*GetPhiWeight(track);
-}
+  f.Close();
+} 
 
-//________________________________________________________________________
-Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
+
+//_________________________________________________________________________
+TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
 {
-  Double_t phiweight=1;
+  TObjArray *acctracks = new TObjArray();
   
-  if (fUsePhiWeight) {
-    Double_t nParticles = fPhiDist->Integral();
-    Double_t nPhibins = fPhiDist->GetNbinsX();
+  AliAODTrack *tr = 0;
+  Int_t maxid1 = 0;
+  Int_t maxidtemp = -1;
+  Float_t ptlow = 0;
+  Float_t ptup = 0;
+  Float_t etalow = 0;
+  Float_t etaup = 0;
+  fESDtrackCuts->GetPtRange(ptlow,ptup);
+  fESDtrackCuts->GetEtaRange(etalow,etaup);
   
-    Double_t Phi = track->Phi();
-    
-    while (Phi<0) Phi += TMath::TwoPi();
-    while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
-      
-    Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::Floor((track->Phi())*nPhibins/TMath::TwoPi()));
-    
-    if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
+  for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
+     tr = aod->GetTrack(i);
+     maxidtemp = tr->GetID(); 
+     if(maxidtemp < 0 && fAODfilterbit != 128) continue;
+     if(maxidtemp > -1 && fAODfilterbit == 128) continue;
+     if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
+     if (maxidtemp > maxid1) maxid1 = maxidtemp;
+     if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow){
+     acctracks->Add(tr);
+     }
   }
-  return phiweight;
+  
+  maxid = maxid1;
+  
+  return acctracks;
+  
 }