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 && (fPeriod.CompareTo("LHC10h") == 0 || fPeriod.CompareTo("LHC11h") == 0)) { // 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
 682
 683   if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist){
 684      Bool_t emptybins;
 685
 686      int iter = 0;
 687      while (iter<3){
 688          emptybins = kFALSE;
 689
 690          for (int i=1; i<fPhiDist[0]->GetNbinsX(); i++){
 691            if (!((fPhiDist[0]->GetBinContent(i))>0)) {
 692              emptybins = kTRUE;
 693            }
 694          }
 695          if (emptybins) {
 696            cout << "empty bins - rebinning!" << endl;
 697            fPhiDist[0]->Rebin();
 698            iter++;
 699          }
 700          else iter = 3;
 701      }
 702
 703      if (emptybins) {
 704        AliError("After Maximum of rebinning still empty Phi-bins!!!");
 705      }
 706   }
 707   if (fPeriod.CompareTo("LHC10h") != 0 && fPeriod.CompareTo("LHC11h") != 0 && !fUserphidist){
 708   AliInfo("No Phi-weights available. All Phi weights set to 1");
 709   SetUsePhiWeight(kFALSE);
 710   }
 711 }
 712
 713 //__________________________________________________________________________
 714 void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
 715 {
 716
 717   fUserphidist = kTRUE;
 718
 719   TFile f(infilename);
 720   TObject* list = f.Get(listname);
 721   fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi");
 722   if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!");
 723
 724   f.Close();
 725 }
 726
 727
 728 //_________________________________________________________________________
 729 TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
 730 {
 731   TObjArray *acctracks = new TObjArray();
 732
 733   AliAODTrack *tr = 0;
 734   Int_t maxid1 = 0;
 735   Int_t maxidtemp = -1;
 736   Float_t ptlow = 0;
 737   Float_t ptup = 0;
 738   Float_t etalow = 0;
 739   Float_t etaup = 0;
 740   fESDtrackCuts->GetPtRange(ptlow,ptup);
 741   fESDtrackCuts->GetEtaRange(etalow,etaup);
 742
 743   for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
 744      tr = aod->GetTrack(i);
 745      maxidtemp = tr->GetID();
 746      if(maxidtemp < 0 && fAODfilterbit != 128) continue;
 747      if(maxidtemp > -1 && fAODfilterbit == 128) continue;
 748      if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
 749      if (maxidtemp > maxid1) maxid1 = maxidtemp;
 750      if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow){
 751      acctracks->Add(tr);
 752      }
 753   }
 754
 755   maxid = maxid1;
 756
 757   return acctracks;
 758
 759 }
 760
 761 //_________________________________________________________________________
 762 void AliEPSelectionTask::SetOADBandPeriod()
 763 {
 764   TString oadbfilename;
 765
 766   if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h
 767     {fPeriod = "LHC10h";
 768      if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
 769
 770         if (fAnalysisInput.CompareTo("AOD")==0){
 771            oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
 772            } else if (fAnalysisInput.CompareTo("ESD")==0){
 773            oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
 774            }
 775
 776        TFile foadb(oadbfilename);
 777        if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
 778
 779        AliInfo("Using Standard OADB");
 780        fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
 781        if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
 782        foadb.Close();
 783        }
 784      }
 785
 786   if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h
 787      {fPeriod = "LHC11h";
 788       if (!fUserphidist) {
 789       // if it's already set and custom class is required, we use the one provided by the user
 790
 791       oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
 792       TFile *foadb = TFile::Open(oadbfilename);
 793       if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
 794
 795       AliInfo("Using Standard OADB");
 796       fSparseDist = (THnSparse*) foadb->Get("Default");
 797       if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
 798       foadb->Close();
 799       if(!fHruns){
 800            fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
 801            fHruns->SetName("runsHisto");
 802       }
 803       }
 804      }
 805 }
 806
 807 //_________________________________________________________________________
 808 TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track)
 809 {
 810   if (fPeriod.CompareTo("LHC10h")==0  || fUserphidist) return fPhiDist[0];
 811   else if(fPeriod.CompareTo("LHC11h")==0)
 812     {
 813      if (track->Charge() < 0)
 814        {
 815         if(track->Eta() < 0.)       return fPhiDist[0];
 816         else if (track->Eta() > 0.) return fPhiDist[2];
 817        }
 818       else if (track->Charge() > 0)
 819        {
 820         if(track->Eta() < 0.)       return fPhiDist[1];
 821         else if (track->Eta() > 0.) return fPhiDist[3];
 822        }
 823
 824     }
 825   return 0;
 826 }
 827
 828 TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd)
 829 {
 830   // 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
 831   TObjArray *acctracks = new TObjArray();
 832   acctracks->SetOwner(kTRUE);
 833
 834   const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
 835
 836   Float_t ptlow = 0;
 837   Float_t ptup = 0;
 838   Float_t etalow = 0;
 839   Float_t etaup = 0;
 840   fESDtrackCuts->GetPtRange(ptlow,ptup);
 841   fESDtrackCuts->GetEtaRange(etalow,etaup);
 842
 843   Double_t pDCA[3] = { 0. }; // momentum at DCA
 844   Double_t rDCA[3] = { 0. }; // position at DCA
 845   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
 846   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
 847
 848
 849
 850   for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){
 851
 852     AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise  need to work with a copy
 853     //
 854     // Track selection
 855     if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue;
 856
 857     // create a tpc only tracl
 858     AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
 859     if(!track) continue;
 860
 861     if(track->Pt()>0.)
 862     {
 863       // only constrain tracks above threshold
 864       AliExternalTrackParam exParam;
 865       // take the B-field from the ESD, no 3D fieldMap available at this point
 866       Bool_t relate = false;
 867       relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
 868       if(!relate){
 869         delete track;
 870         continue;
 871       }
 872       // fetch the track parameters at the DCA (unconstraint)
 873       if(track->GetTPCInnerParam()){
 874         track->GetTPCInnerParam()->GetPxPyPz(pDCA);
 875         track->GetTPCInnerParam()->GetXYZ(rDCA);
 876       }
 877       // get the DCA to the vertex:
 878       track->GetImpactParametersTPC(dDCA,cDCA);
 879       // set the constrained parameters to the track
 880       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 881     }
 882
 883
 884     Float_t eta = track->Eta();
 885     Float_t pT = track->Pt();
 886
 887     if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){
 888       delete track;
 889       continue;
 890     }
 891
 892     acctracks->Add(track);
 893    }
 894
 895
 896   return acctracks;
 897
 898 }
 899