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 <THnSparse.h>
  29 #include <TProfile.h>
  30 #include <TFile.h>
  31 #include <TObjString.h>
  32 #include <TString.h>
  33 #include <TCanvas.h>
  34 #include <TROOT.h>
  35 #include <TDirectory.h>
  36 #include <TSystem.h>
  37 #include <iostream>
  38 #include <TRandom2.h>
  39 #include <TArrayF.h>
  40
  41 #include "AliAnalysisManager.h"
  42 #include "AliVEvent.h"
  43 #include "AliESD.h"
  44 #include "AliESDEvent.h"
  45 #include "AliMCEvent.h"
  46 #include "AliESDtrack.h"
  47 #include "AliESDtrackCuts.h"
  48 #include "AliESDHeader.h"
  49 #include "AliESDInputHandler.h"
  50 #include "AliAODHandler.h"
  51 #include "AliAODEvent.h"
  52 #include "AliHeader.h"
  53 #include "AliGenHijingEventHeader.h"
  54 #include "AliAnalysisTaskSE.h"
  55 #include "AliPhysicsSelectionTask.h"
  56 #include "AliPhysicsSelection.h"
  57 #include "AliBackgroundSelection.h"
  58 #include "AliESDUtils.h"
  59 #include "AliOADBContainer.h"
  60 #include "AliAODMCHeader.h"
  61 #include "AliAODTrack.h"
  62 #include "AliVTrack.h"
  63 #include "AliEventplane.h"
  64
  65 ClassImp(AliEPSelectionTask)
  66
  67 //________________________________________________________________________
  68 AliEPSelectionTask::AliEPSelectionTask():
  69 AliAnalysisTaskSE(),
  70   fAnalysisInput("ESD"),
  71   fTrackType("TPC"),
  72   fPeriod(""),
  73   fUseMCRP(kFALSE),
  74   fUsePhiWeight(kFALSE),
  75   fUsePtWeight(kFALSE),
  76   fSaveTrackContribution(kFALSE),
  77   fUserphidist(kFALSE),
  78   fUsercuts(kFALSE),
  79   fRunNumber(-15),
  80   fAODfilterbit(1),
  81   fEtaGap(0.),
  82   fSplitMethod(0),
  83   fESDtrackCuts(0),
  84   fEPContainer(0),
  85   fSparseDist(0),
  86   fHruns(0),
  87   fQVector(0),
  88   fQContributionX(0),
  89   fQContributionY(0),
  90   fEventplaneQ(0),
  91   fQsub1(0),
  92   fQsub2(0),
  93   fQsubRes(0),
  94   fOutputList(0),
  95   fHOutEventplaneQ(0),
  96   fHOutPhi(0),
  97   fHOutPhiCorr(0),
  98   fHOutsub1sub2(0),
  99   fHOutNTEPRes(0),
 100   fHOutPTPsi(0),
 101   fHOutDiff(0),
 102   fHOutleadPTPsi(0)
 103 {
 104   // Default constructor
 105   AliInfo("Event Plane Selection enabled.");
 106   for(Int_t i = 0; i < 4; ++i) {
 107      fPhiDist[i] = 0;
 108   }
 109 }
 110
 111 //________________________________________________________________________
 112 AliEPSelectionTask::AliEPSelectionTask(const char *name):
 113   AliAnalysisTaskSE(name),
 114   fAnalysisInput("ESD"),
 115   fTrackType("TPC"),
 116   fPeriod(""),
 117   fUseMCRP(kFALSE),
 118   fUsePhiWeight(kFALSE),
 119   fUsePtWeight(kFALSE),
 120   fSaveTrackContribution(kFALSE),
 121   fUserphidist(kFALSE),
 122   fUsercuts(kFALSE),
 123   fRunNumber(-15),
 124   fAODfilterbit(1),
 125   fEtaGap(0.),
 126   fSplitMethod(0),
 127   fESDtrackCuts(0),
 128   fEPContainer(0),
 129   fSparseDist(0),
 130   fHruns(0),
 131   fQVector(0),
 132   fQContributionX(0),
 133   fQContributionY(0),
 134   fEventplaneQ(0),
 135   fQsub1(0),
 136   fQsub2(0),
 137   fQsubRes(0),
 138   fOutputList(0),
 139   fHOutEventplaneQ(0),
 140   fHOutPhi(0),
 141   fHOutPhiCorr(0),
 142   fHOutsub1sub2(0),
 143   fHOutNTEPRes(0),
 144   fHOutPTPsi(0),
 145   fHOutDiff(0),
 146   fHOutleadPTPsi(0)
 147 {
 148   // Default constructor
 149   AliInfo("Event Plane Selection enabled.");
 150   DefineOutput(1, TList::Class());
 151   for(Int_t i = 0; i < 4; i++) {
 152      fPhiDist[i] = 0;
 153   }
 154 }
 155
 156 //________________________________________________________________________
 157 AliEPSelectionTask::~AliEPSelectionTask()
 158 {
 159   // Destructor
 160   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
 161       delete fOutputList;
 162       fOutputList = 0;
 163   }
 164   if (fESDtrackCuts){
 165       delete fESDtrackCuts;
 166       fESDtrackCuts = 0;
 167   }
 168   if (fUserphidist) {
 169     if (fPhiDist[0]) {
 170       delete fPhiDist[0];
 171       fPhiDist[0] = 0;
 172     }
 173   }
 174   if (fEPContainer){
 175       delete fEPContainer;
 176       fEPContainer = 0;
 177   }
 178   if (fPeriod.CompareTo("LHC11h")==0){
 179       for(Int_t i = 0; i < 4; i++) {
 180         if(fPhiDist[i]){
 181           delete fPhiDist[i];
 182           fPhiDist[i] = 0;
 183         }
 184       }
 185       if(fHruns) delete fHruns;
 186   }
 187 }
 188
 189 //________________________________________________________________________
 190 void AliEPSelectionTask::UserCreateOutputObjects()
 191 {
 192   // Create the output containers
 193   if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
 194   AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
 195
 196   fOutputList = new TList();
 197   fOutputList->SetOwner();
 198   fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
 199   fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
 200   fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
 201   fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
 202   fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
 203   fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
 204   fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
 205   fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
 206
 207   fOutputList->Add(fHOutEventplaneQ);
 208   fOutputList->Add(fHOutPhi);
 209   fOutputList->Add(fHOutPhiCorr);
 210   fOutputList->Add(fHOutsub1sub2);
 211   fOutputList->Add(fHOutNTEPRes);
 212   fOutputList->Add(fHOutPTPsi);
 213   fOutputList->Add(fHOutleadPTPsi);
 214   fOutputList->Add(fHOutDiff);
 215
 216   PostData(1, fOutputList);
 217
 218
 219 }
 220
 221 //________________________________________________________________________
 222 void AliEPSelectionTask::UserExec(Option_t */*option*/)
 223 {
 224   // Execute analysis for current event:
 225   if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
 226
 227 //   fRunNumber = -15;
 228
 229   AliEventplane *esdEP;
 230   TVector2 qq1;
 231   TVector2 qq2;
 232   Double_t fRP = 0.; // monte carlo reaction plane angle
 233
 234   if (fAnalysisInput.CompareTo("ESD")==0){
 235
 236     AliVEvent* event = InputEvent();
 237     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 238     if (esd){
 239       if (!(fRunNumber == esd->GetRunNumber())) {
 240           fRunNumber = esd->GetRunNumber();
 241           AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber));
 242           SetOADBandPeriod();
 243           SetPhiDist();
 244       }
 245
 246
 247       if (fUseMCRP) {
 248           AliMCEvent* mcEvent  = MCEvent();
 249           if (mcEvent && mcEvent->GenEventHeader()) {
 250               AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
 251               if (headerH) fRP = headerH->ReactionPlaneAngle();
 252           }
 253       }
 254
 255       esdEP = esd->GetEventplane();
 256       if (fSaveTrackContribution) {
 257         esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
 258         esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
 259         esdEP->GetQContributionXArraysub1()->Set(esd->GetNumberOfTracks());
 260         esdEP->GetQContributionYArraysub1()->Set(esd->GetNumberOfTracks());
 261         esdEP->GetQContributionXArraysub2()->Set(esd->GetNumberOfTracks());
 262         esdEP->GetQContributionYArraysub2()->Set(esd->GetNumberOfTracks());
 263       }
 264
 265       TObjArray* tracklist = new TObjArray;
 266       if (fTrackType.CompareTo("GLOBAL")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
 267       if (fTrackType.CompareTo("TPC")==0 && fPeriod.CompareTo("LHC10h")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
 268       else if (fTrackType.CompareTo("TPC")==0 && fPeriod.CompareTo("LHC11h")==0) tracklist = GetTracksForLHC11h(esd);
 269       const int nt = tracklist->GetEntries();
 270
 271       if (nt>4){
 272         fQVector = new TVector2(GetQ(esdEP,tracklist));
 273         fEventplaneQ = fQVector->Phi()/2;
 274         GetQsub(qq1, qq2, tracklist, esdEP);
 275         fQsub1 = new TVector2(qq1);
 276         fQsub2 = new TVector2(qq2);
 277         fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
 278
 279         esdEP->SetQVector(fQVector);
 280         esdEP->SetEventplaneQ(fEventplaneQ);
 281         esdEP->SetQsub(fQsub1,fQsub2);
 282         esdEP->SetQsubRes(fQsubRes);
 283
 284         fHOutEventplaneQ->Fill(fEventplaneQ);
 285         fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
 286         fHOutNTEPRes->Fill(nt,fQsubRes);
 287
 288         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
 289
 290         for (int iter = 0; iter<nt;iter++){
 291           AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
 292           if (track) {
 293             float delta = track->Phi()-fEventplaneQ;
 294             while (delta < 0) delta += TMath::Pi();
 295             while (delta > TMath::Pi()) delta -= TMath::Pi();
 296             fHOutPTPsi->Fill(track->Pt(),delta);
 297             fHOutPhi->Fill(track->Phi());
 298             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
 299           }
 300         }
 301
 302         AliESDtrack* trmax = esd->GetTrack(0);
 303         for (int iter = 1; iter<nt;iter++){
 304           AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
 305           if (track && (track->Pt() > trmax->Pt())) trmax = track;
 306         }
 307         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
 308       }
 309       tracklist->Clear();
 310       delete tracklist;
 311       tracklist = 0;
 312     }
 313   }
 314
 315     else if (fAnalysisInput.CompareTo("AOD")==0){
 316     AliVEvent* event = InputEvent();
 317     AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
 318
 319     if (aod){
 320       if (!(fRunNumber == aod->GetRunNumber())) {
 321         fRunNumber = aod->GetRunNumber();
 322         AliInfo(Form("Changing Phi-distribution to run %d",fRunNumber));
 323         SetOADBandPeriod();
 324         SetPhiDist();
 325       }
 326
 327       if (fUseMCRP) {
 328         AliAODMCHeader *headerH = dynamic_cast<AliAODMCHeader*>(aod->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 329         if (headerH) fRP = headerH->GetReactionPlaneAngle();
 330       }
 331
 332       esdEP = aod->GetHeader()->GetEventplaneP();
 333       esdEP->Reset();
 334
 335       Int_t maxID = 0;
 336       TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID);
 337
 338       if (fSaveTrackContribution) {
 339         esdEP->GetQContributionXArray()->Set(maxID+1);
 340         esdEP->GetQContributionYArray()->Set(maxID+1);
 341         esdEP->GetQContributionXArraysub1()->Set(maxID+1);
 342         esdEP->GetQContributionYArraysub1()->Set(maxID+1);
 343         esdEP->GetQContributionXArraysub2()->Set(maxID+1);
 344         esdEP->GetQContributionYArraysub2()->Set(maxID+1);
 345       }
 346
 347       const int NT = tracklist->GetEntries();
 348
 349       if (NT>4){
 350         fQVector = new TVector2(GetQ(esdEP,tracklist));
 351         fEventplaneQ = fQVector->Phi()/2.;
 352         GetQsub(qq1, qq2, tracklist, esdEP);
 353         fQsub1 = new TVector2(qq1);
 354         fQsub2 = new TVector2(qq2);
 355         fQsubRes = (fQsub1->Phi()/2. - fQsub2->Phi()/2.);
 356
 357         esdEP->SetQVector(fQVector);
 358         esdEP->SetEventplaneQ(fEventplaneQ);
 359         esdEP->SetQsub(fQsub1,fQsub2);
 360         esdEP->SetQsubRes(fQsubRes);
 361
 362         fHOutEventplaneQ->Fill(fEventplaneQ);
 363         fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.);
 364         fHOutNTEPRes->Fill(NT,fQsubRes);
 365
 366         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
 367
 368         for (int iter = 0; iter<NT;iter++){
 369           AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
 370           if (track) {
 371             float delta = track->Phi()-fEventplaneQ;
 372             while (delta < 0) delta += TMath::Pi();
 373             while (delta > TMath::Pi()) delta -= TMath::Pi();
 374             fHOutPTPsi->Fill(track->Pt(),delta);
 375             fHOutPhi->Fill(track->Phi());
 376             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
 377           }
 378         }
 379
 380         AliAODTrack* trmax = aod->GetTrack(0);
 381         for (int iter = 1; iter<NT;iter++){
 382           AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
 383           if (track && (track->Pt() > trmax->Pt())) trmax = track;
 384         }
 385         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
 386       }
 387       delete tracklist;
 388       tracklist = 0;
 389     }
 390
 391
 392   }
 393
 394
 395   else {
 396     printf(" Analysis Input not known!\n\n ");
 397     return;
 398   }
 399   PostData(1, fOutputList);
 400 }
 401
 402 //________________________________________________________________________
 403 void AliEPSelectionTask::Terminate(Option_t */*option*/)
 404 {
 405   // Terminate analysis
 406 }
 407
 408 //__________________________________________________________________________
 409 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
 410 {
 411 // Get the Q vector
 412   TVector2 mQ;
 413   float mQx=0, mQy=0;
 414   AliVTrack* track;
 415   Double_t weight;
 416   Int_t idtemp = -1;
 417
 418   int nt = tracklist->GetEntries();
 419
 420   for (int i=0; i<nt; i++){
 421     weight = 1;
 422     track = dynamic_cast<AliVTrack*> (tracklist->At(i));
 423     if (track) {
 424       weight = GetWeight(track);
 425     if (fSaveTrackContribution){
 426       idtemp = track->GetID();
 427       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
 428       EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),idtemp);
 429       EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),idtemp);
 430      }
 431      mQx += (weight*cos(2*track->Phi()));
 432      mQy += (weight*sin(2*track->Phi()));
 433     }
 434   }
 435   mQ.Set(mQx,mQy);
 436   return mQ;
 437 }
 438
 439   //________________________________________________________________________
 440 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP)
 441 {
 442 // Get Qsub
 443   TVector2 mQ[2];
 444   float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
 445   Double_t weight;
 446
 447   AliVTrack* track;
 448   TRandom2 rn = 0;
 449
 450   int nt = tracklist->GetEntries();
 451   int trackcounter1=0, trackcounter2=0;
 452   int idtemp = 0;
 453
 454   if (fSplitMethod == AliEPSelectionTask::kRandom){
 455
 456     for (Int_t i = 0; i < nt; i++) {
 457       weight = 1;
 458       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
 459       if (!track) continue;
 460       weight = GetWeight(track);
 461       idtemp = track->GetID();
 462       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
 463
 464       // This loop splits the track set into 2 random subsets
 465       if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
 466         float random = rn.Rndm();
 467         if(random < .5){
 468           mQx1 += (weight*cos(2*track->Phi()));
 469           mQy1 += (weight*sin(2*track->Phi()));
 470           if (fSaveTrackContribution){
 471             EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
 472             EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
 473           }
 474           trackcounter1++;
 475         }
 476         else {
 477           mQx2 += (weight*cos(2*track->Phi()));
 478           mQy2 += (weight*sin(2*track->Phi()));
 479           if (fSaveTrackContribution){
 480             EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
 481             EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
 482           }
 483           trackcounter2++;
 484         }
 485       }
 486       else if( trackcounter1 >= int(nt/2.)){
 487         mQx2 += (weight*cos(2*track->Phi()));
 488         mQy2 += (weight*sin(2*track->Phi()));
 489         if (fSaveTrackContribution){
 490           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
 491           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
 492         }
 493         trackcounter2++;
 494       }
 495       else {
 496         mQx1 += (weight*cos(2*track->Phi()));
 497         mQy1 += (weight*sin(2*track->Phi()));
 498         if (fSaveTrackContribution){
 499           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
 500           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
 501         }
 502         trackcounter1++;
 503       }
 504     }
 505   } else if (fSplitMethod == AliEPSelectionTask::kEta) {
 506
 507     for (Int_t i = 0; i < nt; i++) {
 508       weight = 1;
 509       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
 510       if (!track) continue;
 511       weight = GetWeight(track);
 512       Double_t eta = track->Eta();
 513       idtemp = track->GetID();
 514       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
 515
 516       if (eta > fEtaGap/2.) {
 517         mQx1 += (weight*cos(2*track->Phi()));
 518         mQy1 += (weight*sin(2*track->Phi()));
 519         if (fSaveTrackContribution){
 520           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
 521           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
 522         }
 523       } else if (eta < -1.*fEtaGap/2.) {
 524         mQx2 += (weight*cos(2*track->Phi()));
 525         mQy2 += (weight*sin(2*track->Phi()));
 526         if (fSaveTrackContribution){
 527           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
 528           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
 529         }
 530       }
 531     }
 532   } else {
 533     printf("plane resolution determination method not available!\n\n ");
 534     return;
 535   }
 536
 537   mQ[0].Set(mQx1,mQy1);
 538   mQ[1].Set(mQx2,mQy2);
 539   Q1 = mQ[0];
 540   Q2 = mQ[1];
 541 }
 542
 543 //________________________________________________________________________
 544 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
 545
 546   if(fESDtrackCuts){
 547     delete fESDtrackCuts;
 548     fESDtrackCuts = 0;
 549   }
 550   if (fAnalysisInput.CompareTo("AOD")==0){
 551     AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");
 552     fUsercuts = kFALSE;
 553     SetTrackType("TPC");
 554     return;
 555   }
 556   fUsercuts = kTRUE;
 557   fESDtrackCuts = trackcuts;
 558 }
 559
 560 //________________________________________________________________________
 561 void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup){
 562
 563   if(fESDtrackCuts){
 564     delete fESDtrackCuts;
 565     fESDtrackCuts = 0;
 566   }
 567   if (fAnalysisInput.CompareTo("ESD")==0){
 568     AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");
 569     fUsercuts = kFALSE;
 570     SetTrackType("TPC");
 571     return;
 572   }
 573   fUsercuts = kTRUE;
 574   fESDtrackCuts = new AliESDtrackCuts();
 575   fESDtrackCuts->SetPtRange(ptlow,ptup);
 576   fESDtrackCuts->SetEtaRange(etalow,etaup);
 577   fAODfilterbit = filterbit;
 578 }
 579
 580 //_____________________________________________________________________________
 581
 582 void AliEPSelectionTask::SetTrackType(TString tracktype){
 583   fTrackType = tracktype;
 584   if (!fUsercuts) {
 585   if (fTrackType.CompareTo("GLOBAL")==0){
 586     fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
 587     fAODfilterbit = 32;
 588   }
 589   if (fTrackType.CompareTo("TPC")==0){
 590     fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 591     fAODfilterbit = 128;
 592   }
 593   fESDtrackCuts->SetPtRange(0.15,20.);
 594   fESDtrackCuts->SetEtaRange(-0.8,0.8);
 595   }
 596 }
 597
 598 //________________________________________________________________________
 599 Double_t AliEPSelectionTask::GetWeight(TObject* track1)
 600 {
 601   Double_t ptweight=1;
 602   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
 603   if (fUsePtWeight && track) {
 604     if (track->Pt()<2) ptweight=track->Pt();
 605     else ptweight=2;
 606   }
 607   return ptweight*GetPhiWeight(track);
 608 }
 609
 610 //________________________________________________________________________
 611 Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
 612 {
 613   Double_t phiweight=1;
 614   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
 615
 616   TH1F *phiDist = 0x0;
 617   if(track) phiDist = SelectPhiDist(track);
 618
 619   if (fUsePhiWeight && phiDist && track) {
 620     Double_t nParticles = phiDist->Integral();
 621     Double_t nPhibins = phiDist->GetNbinsX();
 622
 623     Double_t Phi = track->Phi();
 624
 625     while (Phi<0) Phi += TMath::TwoPi();
 626     while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
 627
 628     Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
 629
 630     if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
 631   }
 632   return phiweight;
 633 }
 634
 635 //__________________________________________________________________________
 636 void AliEPSelectionTask::SetPhiDist()
 637 {
 638   if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
 639
 640     if (fPeriod.CompareTo("LHC10h")==0)
 641        {
 642         fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");}
 643         else if(fPeriod.CompareTo("LHC11h")==0){
 644             Int_t runbin=fHruns->FindBin(fRunNumber);
 645             if (fHruns->GetBinContent(runbin) > 1){
 646               fSparseDist->GetAxis(0)->SetRange(runbin,runbin);
 647               }
 648             else if(fHruns->GetBinContent(runbin) < 2){
 649                fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights
 650                AliInfo("Using integrated Phi-weights for this run");
 651                }
 652             for (Int_t i = 0; i<4 ;i++)
 653             {
 654              if(fPhiDist[i]){
 655                delete fPhiDist[i];
 656                fPhiDist[i] = 0x0;
 657                }
 658              if(i == 0){
 659                fSparseDist->GetAxis(1)->SetRange(1,1);  // neg charge
 660                fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
 661              if(i == 1){
 662                fSparseDist->GetAxis(1)->SetRange(2,2);  // pos charge
 663                fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
 664              if(i == 2){
 665                fSparseDist->GetAxis(1)->SetRange(1,1);  // neg charge
 666                fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
 667              if(i == 3){
 668                fSparseDist->GetAxis(1)->SetRange(2,2);  // pos charge
 669                fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
 670              fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi
 671              fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber));
 672              fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes
 673              fSparseDist->GetAxis(2)->SetRange(1,2);
 674              }
 675              fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis
 676        }
 677
 678     if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
 679
 680   }
 681   else {
 682     AliInfo("Using Custom Phi Distribution");
 683   }
 684
 685   if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist){
 686      Bool_t emptybins;
 687
 688      int iter = 0;
 689      while (iter<3){
 690          emptybins = kFALSE;
 691
 692          for (int i=1; i<fPhiDist[0]->GetNbinsX(); i++){
 693            if (!((fPhiDist[0]->GetBinContent(i))>0)) {
 694              emptybins = kTRUE;
 695            }
 696          }
 697          if (emptybins) {
 698            cout << "empty bins - rebinning!" << endl;
 699            fPhiDist[0]->Rebin();
 700            iter++;
 701          }
 702          else iter = 3;
 703      }
 704
 705      if (emptybins) {
 706        AliError("After Maximum of rebinning still empty Phi-bins!!!");
 707      }
 708   }
 709 }
 710
 711 //__________________________________________________________________________
 712 void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
 713 {
 714
 715   fUserphidist = kTRUE;
 716
 717   TFile f(infilename);
 718   TObject* list = f.Get(listname);
 719   fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi");
 720   if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!");
 721
 722   f.Close();
 723 }
 724
 725
 726 //_________________________________________________________________________
 727 TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
 728 {
 729   TObjArray *acctracks = new TObjArray();
 730
 731   AliAODTrack *tr = 0;
 732   Int_t maxid1 = 0;
 733   Int_t maxidtemp = -1;
 734   Float_t ptlow = 0;
 735   Float_t ptup = 0;
 736   Float_t etalow = 0;
 737   Float_t etaup = 0;
 738   fESDtrackCuts->GetPtRange(ptlow,ptup);
 739   fESDtrackCuts->GetEtaRange(etalow,etaup);
 740
 741   for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
 742      tr = aod->GetTrack(i);
 743      maxidtemp = tr->GetID();
 744      if(maxidtemp < 0 && fAODfilterbit != 128) continue;
 745      if(maxidtemp > -1 && fAODfilterbit == 128) continue;
 746      if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
 747      if (maxidtemp > maxid1) maxid1 = maxidtemp;
 748      if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow){
 749      acctracks->Add(tr);
 750      }
 751   }
 752
 753   maxid = maxid1;
 754
 755   return acctracks;
 756
 757 }
 758
 759 //_________________________________________________________________________
 760 void AliEPSelectionTask::SetOADBandPeriod()
 761 {
 762   TString oadbfilename;
 763
 764   if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h
 765     {fPeriod = "LHC10h";
 766      if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
 767
 768         if (fAnalysisInput.CompareTo("AOD")==0){
 769            oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
 770            } else if (fAnalysisInput.CompareTo("ESD")==0){
 771            oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
 772            }
 773
 774        TFile foadb(oadbfilename);
 775        if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
 776
 777        AliInfo("Using Standard OADB");
 778        fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
 779        if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
 780        foadb.Close();
 781        }
 782      }
 783
 784   if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h
 785      {fPeriod = "LHC11h";
 786       if (!fUserphidist) {
 787       // if it's already set and custom class is required, we use the one provided by the user
 788
 789       oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
 790       TFile *foadb = TFile::Open(oadbfilename);
 791       if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
 792
 793       AliInfo("Using Standard OADB");
 794       fSparseDist = (THnSparse*) foadb->Get("Default");
 795       if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
 796       foadb->Close();
 797       if(!fHruns) fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
 798       }
 799      }
 800 }
 801
 802 //_________________________________________________________________________
 803 TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track)
 804 {
 805   if (fPeriod.CompareTo("LHC10h")==0  || fUserphidist) return fPhiDist[0];
 806   else if(fPeriod.CompareTo("LHC11h")==0)
 807     {
 808      if (track->Charge() < 0)
 809        {
 810         if(track->Eta() < 0.)       return fPhiDist[0];
 811         else if (track->Eta() > 0.) return fPhiDist[2];
 812        }
 813       else if (track->Charge() > 0)
 814        {
 815         if(track->Eta() < 0.)       return fPhiDist[1];
 816         else if (track->Eta() > 0.) return fPhiDist[3];
 817        }
 818
 819     }
 820   return 0;
 821 }
 822
 823 TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd)
 824 {
 825   // Need to do this hack beacuse only this type of TPC only tracks in AOD is available and one type of Phi-weights is used
 826   TObjArray *acctracks = new TObjArray();
 827   acctracks->SetOwner(kTRUE);
 828
 829   const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
 830
 831   Float_t ptlow = 0;
 832   Float_t ptup = 0;
 833   Float_t etalow = 0;
 834   Float_t etaup = 0;
 835   fESDtrackCuts->GetPtRange(ptlow,ptup);
 836   fESDtrackCuts->GetEtaRange(etalow,etaup);
 837
 838   Double_t pDCA[3] = { 0. }; // momentum at DCA
 839   Double_t rDCA[3] = { 0. }; // position at DCA
 840   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
 841   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
 842
 843
 844
 845   for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){
 846
 847     AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise  need to work with a copy
 848     //
 849     // Track selection
 850     if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue;
 851
 852     // create a tpc only tracl
 853     AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
 854     if(!track) continue;
 855
 856     if(track->Pt()>0.)
 857     {
 858       // only constrain tracks above threshold
 859       AliExternalTrackParam exParam;
 860       // take the B-field from the ESD, no 3D fieldMap available at this point
 861       Bool_t relate = false;
 862       relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
 863       if(!relate){
 864         delete track;
 865         continue;
 866       }
 867       // fetch the track parameters at the DCA (unconstraint)
 868       if(track->GetTPCInnerParam()){
 869         track->GetTPCInnerParam()->GetPxPyPz(pDCA);
 870         track->GetTPCInnerParam()->GetXYZ(rDCA);
 871       }
 872       // get the DCA to the vertex:
 873       track->GetImpactParametersTPC(dDCA,cDCA);
 874       // set the constrained parameters to the track
 875       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 876     }
 877
 878
 879     Float_t eta = track->Eta();
 880     Float_t pT = track->Pt();
 881
 882     if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){
 883       delete track;
 884       continue;
 885     }
 886
 887     acctracks->Add(track);
 888    }
 889
 890
 891   return acctracks;
 892
 893 }
 894