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 }