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 if (fSplitMethod == AliEPSelectionTask::kCharge) {
 535
 536     for (Int_t i = 0; i < nt; i++) {
 537       weight = 1;
 538       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
 539       if (!track) continue;
 540       weight = GetWeight(track);
 541       Short_t cha = track->Charge();
 542       idtemp = track->GetID();
 543       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
 544
 545       if (cha > 0) {
 546         mQx1 += (weight*cos(2*track->Phi()));
 547         mQy1 += (weight*sin(2*track->Phi()));
 548         if (fSaveTrackContribution){
 549           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
 550           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
 551         }
 552       } else if (cha < 0) {
 553         mQx2 += (weight*cos(2*track->Phi()));
 554         mQy2 += (weight*sin(2*track->Phi()));
 555         if (fSaveTrackContribution){
 556           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
 557           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
 558         }
 559       }
 560     }
 561   } else {
 562     printf("plane resolution determination method not available!\n\n ");
 563     return;
 564   }
 565
 566   mQ[0].Set(mQx1,mQy1);
 567   mQ[1].Set(mQx2,mQy2);
 568   Q1 = mQ[0];
 569   Q2 = mQ[1];
 570 }
 571
 572 //________________________________________________________________________
 573 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
 574
 575   if(fESDtrackCuts){
 576     delete fESDtrackCuts;
 577     fESDtrackCuts = 0;
 578   }
 579   if (fAnalysisInput.CompareTo("AOD")==0){
 580     AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");
 581     fUsercuts = kFALSE;
 582     SetTrackType("TPC");
 583     return;
 584   }
 585   fUsercuts = kTRUE;
 586   fESDtrackCuts = trackcuts;
 587 }
 588
 589 //________________________________________________________________________
 590 void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup){
 591
 592   if(fESDtrackCuts){
 593     delete fESDtrackCuts;
 594     fESDtrackCuts = 0;
 595   }
 596   if (fAnalysisInput.CompareTo("ESD")==0){
 597     AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");
 598     fUsercuts = kFALSE;
 599     SetTrackType("TPC");
 600     return;
 601   }
 602   fUsercuts = kTRUE;
 603   fESDtrackCuts = new AliESDtrackCuts();
 604   fESDtrackCuts->SetPtRange(ptlow,ptup);
 605   fESDtrackCuts->SetEtaRange(etalow,etaup);
 606   fAODfilterbit = filterbit;
 607 }
 608
 609 //_____________________________________________________________________________
 610
 611 void AliEPSelectionTask::SetTrackType(TString tracktype){
 612   fTrackType = tracktype;
 613   if (!fUsercuts) {
 614   if (fTrackType.CompareTo("GLOBAL")==0){
 615     fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
 616     fAODfilterbit = 32;
 617   }
 618   if (fTrackType.CompareTo("TPC")==0){
 619     fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 620     fAODfilterbit = 128;
 621   }
 622   fESDtrackCuts->SetPtRange(0.15,20.);
 623   fESDtrackCuts->SetEtaRange(-0.8,0.8);
 624   }
 625 }
 626
 627 //________________________________________________________________________
 628 Double_t AliEPSelectionTask::GetWeight(TObject* track1)
 629 {
 630   Double_t ptweight=1;
 631   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
 632   if (fUsePtWeight && track) {
 633     if (track->Pt()<2) ptweight=track->Pt();
 634     else ptweight=2;
 635   }
 636   return ptweight*GetPhiWeight(track);
 637 }
 638
 639 //________________________________________________________________________
 640 Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
 641 {
 642   Double_t phiweight=1;
 643   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
 644
 645   TH1F *phiDist = 0x0;
 646   if(track) phiDist = SelectPhiDist(track);
 647
 648   if (fUsePhiWeight && phiDist && track) {
 649     Double_t nParticles = phiDist->Integral();
 650     Double_t nPhibins = phiDist->GetNbinsX();
 651
 652     Double_t Phi = track->Phi();
 653
 654     while (Phi<0) Phi += TMath::TwoPi();
 655     while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
 656
 657     Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
 658
 659     if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
 660   }
 661   return phiweight;
 662 }
 663
 664 //__________________________________________________________________________
 665 void AliEPSelectionTask::SetPhiDist()
 666 {
 667   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
 668
 669     if (fPeriod.CompareTo("LHC10h")==0)
 670        {
 671         fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");}
 672         else if(fPeriod.CompareTo("LHC11h")==0){
 673             Int_t runbin=fHruns->FindBin(fRunNumber);
 674             if (fHruns->GetBinContent(runbin) > 1){
 675               fSparseDist->GetAxis(0)->SetRange(runbin,runbin);
 676               }
 677             else if(fHruns->GetBinContent(runbin) < 2){
 678                fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights
 679                AliInfo("Using integrated Phi-weights for this run");
 680                }
 681             for (Int_t i = 0; i<4 ;i++)
 682             {
 683              if(fPhiDist[i]){
 684                delete fPhiDist[i];
 685                fPhiDist[i] = 0x0;
 686                }
 687              if(i == 0){
 688                fSparseDist->GetAxis(1)->SetRange(1,1);  // neg charge
 689                fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
 690              if(i == 1){
 691                fSparseDist->GetAxis(1)->SetRange(2,2);  // pos charge
 692                fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
 693              if(i == 2){
 694                fSparseDist->GetAxis(1)->SetRange(1,1);  // neg charge
 695                fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
 696              if(i == 3){
 697                fSparseDist->GetAxis(1)->SetRange(2,2);  // pos charge
 698                fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
 699              fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi
 700              fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber));
 701              fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes
 702              fSparseDist->GetAxis(2)->SetRange(1,2);
 703              }
 704              fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis
 705        }
 706
 707     if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
 708
 709   }
 710
 711
 712   if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist){
 713      Bool_t emptybins;
 714
 715      int iter = 0;
 716      while (iter<3){
 717          emptybins = kFALSE;
 718
 719          for (int i=1; i<fPhiDist[0]->GetNbinsX(); i++){
 720            if (!((fPhiDist[0]->GetBinContent(i))>0)) {
 721              emptybins = kTRUE;
 722            }
 723          }
 724          if (emptybins) {
 725            cout << "empty bins - rebinning!" << endl;
 726            fPhiDist[0]->Rebin();
 727            iter++;
 728          }
 729          else iter = 3;
 730      }
 731
 732      if (emptybins) {
 733        AliError("After Maximum of rebinning still empty Phi-bins!!!");
 734      }
 735   }
 736   if (fPeriod.CompareTo("LHC10h") != 0 && fPeriod.CompareTo("LHC11h") != 0 && !fUserphidist){
 737   AliInfo("No Phi-weights available. All Phi weights set to 1");
 738   SetUsePhiWeight(kFALSE);
 739   }
 740 }
 741
 742 //__________________________________________________________________________
 743 void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
 744 {
 745
 746   fUserphidist = kTRUE;
 747
 748   TFile f(infilename);
 749   TObject* list = f.Get(listname);
 750   fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi");
 751   if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!");
 752
 753   f.Close();
 754 }
 755
 756
 757 //_________________________________________________________________________
 758 TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
 759 {
 760   TObjArray *acctracks = new TObjArray();
 761
 762   AliAODTrack *tr = 0;
 763   Int_t maxid1 = 0;
 764   Int_t maxidtemp = -1;
 765   Float_t ptlow = 0;
 766   Float_t ptup = 0;
 767   Float_t etalow = 0;
 768   Float_t etaup = 0;
 769   fESDtrackCuts->GetPtRange(ptlow,ptup);
 770   fESDtrackCuts->GetEtaRange(etalow,etaup);
 771
 772   for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
 773      tr = aod->GetTrack(i);
 774      maxidtemp = tr->GetID();
 775      if(maxidtemp < 0 && fAODfilterbit != 128) continue;
 776      if(maxidtemp > -1 && fAODfilterbit == 128) continue;
 777      if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
 778      if (maxidtemp > maxid1) maxid1 = maxidtemp;
 779      if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow){
 780      acctracks->Add(tr);
 781      }
 782   }
 783
 784   maxid = maxid1;
 785
 786   return acctracks;
 787
 788 }
 789
 790 //_________________________________________________________________________
 791 void AliEPSelectionTask::SetOADBandPeriod()
 792 {
 793   TString oadbfilename;
 794
 795   if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h
 796     {fPeriod = "LHC10h";
 797      if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
 798
 799         if (fAnalysisInput.CompareTo("AOD")==0){
 800            oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
 801            } else if (fAnalysisInput.CompareTo("ESD")==0){
 802            oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
 803            }
 804
 805        TFile foadb(oadbfilename);
 806        if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
 807
 808        AliInfo("Using Standard OADB");
 809        fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
 810        if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
 811        foadb.Close();
 812        }
 813      }
 814
 815   if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h
 816      {fPeriod = "LHC11h";
 817       if (!fUserphidist) {
 818       // if it's already set and custom class is required, we use the one provided by the user
 819
 820       oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
 821       TFile *foadb = TFile::Open(oadbfilename);
 822       if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
 823
 824       AliInfo("Using Standard OADB");
 825       fSparseDist = (THnSparse*) foadb->Get("Default");
 826       if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
 827       foadb->Close();
 828       if(!fHruns){
 829            fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
 830            fHruns->SetName("runsHisto");
 831       }
 832       }
 833      }
 834 }
 835
 836 //_________________________________________________________________________
 837 TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track)
 838 {
 839   if (fPeriod.CompareTo("LHC10h")==0  || fUserphidist) return fPhiDist[0];
 840   else if(fPeriod.CompareTo("LHC11h")==0)
 841     {
 842      if (track->Charge() < 0)
 843        {
 844         if(track->Eta() < 0.)       return fPhiDist[0];
 845         else if (track->Eta() > 0.) return fPhiDist[2];
 846        }
 847       else if (track->Charge() > 0)
 848        {
 849         if(track->Eta() < 0.)       return fPhiDist[1];
 850         else if (track->Eta() > 0.) return fPhiDist[3];
 851        }
 852
 853     }
 854   return 0;
 855 }
 856
 857 TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd)
 858 {
 859   // 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
 860   TObjArray *acctracks = new TObjArray();
 861   acctracks->SetOwner(kTRUE);
 862
 863   const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
 864
 865   Float_t ptlow = 0;
 866   Float_t ptup = 0;
 867   Float_t etalow = 0;
 868   Float_t etaup = 0;
 869   fESDtrackCuts->GetPtRange(ptlow,ptup);
 870   fESDtrackCuts->GetEtaRange(etalow,etaup);
 871
 872   Double_t pDCA[3] = { 0. }; // momentum at DCA
 873   Double_t rDCA[3] = { 0. }; // position at DCA
 874   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
 875   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
 876
 877
 878
 879   for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){
 880
 881     AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise  need to work with a copy
 882     //
 883     // Track selection
 884     if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue;
 885
 886     // create a tpc only tracl
 887     AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
 888     if(!track) continue;
 889
 890     if(track->Pt()>0.)
 891     {
 892       // only constrain tracks above threshold
 893       AliExternalTrackParam exParam;
 894       // take the B-field from the ESD, no 3D fieldMap available at this point
 895       Bool_t relate = false;
 896       relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
 897       if(!relate){
 898         delete track;
 899         continue;
 900       }
 901       // fetch the track parameters at the DCA (unconstraint)
 902       if(track->GetTPCInnerParam()){
 903         track->GetTPCInnerParam()->GetPxPyPz(pDCA);
 904         track->GetTPCInnerParam()->GetXYZ(rDCA);
 905       }
 906       // get the DCA to the vertex:
 907       track->GetImpactParametersTPC(dDCA,cDCA);
 908       // set the constrained parameters to the track
 909       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 910     }
 911
 912
 913     Float_t eta = track->Eta();
 914     Float_t pT = track->Pt();
 915
 916     if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){
 917       delete track;
 918       continue;
 919     }
 920
 921     acctracks->Add(track);
 922    }
 923
 924
 925   return acctracks;
 926
 927 }
 928