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