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