ANALYSIS/AliEPSelectionTask.cxx

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