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