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       if(!esdEP) return; // protection against missing EP branch (nanoAODs)
 336       esdEP->Reset();
 337
 338       Int_t maxID = 0;
 339       TObjArray* tracklist = GetAODTracksAndMaxID(aod,maxID);
 340
 341       if (fSaveTrackContribution) {
 342         esdEP->GetQContributionXArray()->Set(maxID+1);
 343         esdEP->GetQContributionYArray()->Set(maxID+1);
 344         esdEP->GetQContributionXArraysub1()->Set(maxID+1);
 345         esdEP->GetQContributionYArraysub1()->Set(maxID+1);
 346         esdEP->GetQContributionXArraysub2()->Set(maxID+1);
 347         esdEP->GetQContributionYArraysub2()->Set(maxID+1);
 348       }
 349
 350       const int NT = tracklist->GetEntries();
 351
 352       if (NT>4){
 353         fQVector = new TVector2(GetQ(esdEP,tracklist));
 354         fEventplaneQ = fQVector->Phi()/2.;
 355         GetQsub(qq1, qq2, tracklist, esdEP);
 356         fQsub1 = new TVector2(qq1);
 357         fQsub2 = new TVector2(qq2);
 358         fQsubRes = (fQsub1->Phi()/2. - fQsub2->Phi()/2.);
 359
 360         esdEP->SetQVector(fQVector);
 361         esdEP->SetEventplaneQ(fEventplaneQ);
 362         esdEP->SetQsub(fQsub1,fQsub2);
 363         esdEP->SetQsubRes(fQsubRes);
 364
 365         fHOutEventplaneQ->Fill(fEventplaneQ);
 366         fHOutsub1sub2->Fill(fQsub1->Phi()/2.,fQsub2->Phi()/2.);
 367         fHOutNTEPRes->Fill(NT,fQsubRes);
 368
 369         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
 370
 371         for (int iter = 0; iter<NT;iter++){
 372           AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
 373           if (track) {
 374             float delta = track->Phi()-fEventplaneQ;
 375             while (delta < 0) delta += TMath::Pi();
 376             while (delta > TMath::Pi()) delta -= TMath::Pi();
 377             fHOutPTPsi->Fill(track->Pt(),delta);
 378             fHOutPhi->Fill(track->Phi());
 379             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
 380           }
 381         }
 382
 383         AliAODTrack* trmax = aod->GetTrack(0);
 384         for (int iter = 1; iter<NT;iter++){
 385           AliAODTrack* track = dynamic_cast<AliAODTrack*> (tracklist->At(iter));
 386           if (track && (track->Pt() > trmax->Pt())) trmax = track;
 387         }
 388         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
 389       }
 390       delete tracklist;
 391       tracklist = 0;
 392     }
 393
 394
 395   }
 396
 397
 398   else {
 399     printf(" Analysis Input not known!\n\n ");
 400     return;
 401   }
 402   PostData(1, fOutputList);
 403 }
 404
 405 //________________________________________________________________________
 406 void AliEPSelectionTask::Terminate(Option_t */*option*/)
 407 {
 408   // Terminate analysis
 409 }
 410
 411 //__________________________________________________________________________
 412 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
 413 {
 414 // Get the Q vector
 415   TVector2 mQ;
 416   float mQx=0, mQy=0;
 417   AliVTrack* track;
 418   Double_t weight;
 419   Int_t idtemp = -1;
 420
 421   int nt = tracklist->GetEntries();
 422
 423   for (int i=0; i<nt; i++){
 424     weight = 1;
 425     track = dynamic_cast<AliVTrack*> (tracklist->At(i));
 426     if (track) {
 427       weight = GetWeight(track);
 428     if (fSaveTrackContribution){
 429       idtemp = track->GetID();
 430       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
 431       EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),idtemp);
 432       EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),idtemp);
 433      }
 434      mQx += (weight*cos(2*track->Phi()));
 435      mQy += (weight*sin(2*track->Phi()));
 436     }
 437   }
 438   mQ.Set(mQx,mQy);
 439   return mQ;
 440 }
 441
 442   //________________________________________________________________________
 443 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist,AliEventplane* EP)
 444 {
 445 // Get Qsub
 446   TVector2 mQ[2];
 447   float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
 448   Double_t weight;
 449
 450   AliVTrack* track;
 451   TRandom2 rn = 0;
 452
 453   int nt = tracklist->GetEntries();
 454   int trackcounter1=0, trackcounter2=0;
 455   int idtemp = 0;
 456
 457   if (fSplitMethod == AliEPSelectionTask::kRandom){
 458
 459     for (Int_t i = 0; i < nt; i++) {
 460       weight = 1;
 461       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
 462       if (!track) continue;
 463       weight = GetWeight(track);
 464       idtemp = track->GetID();
 465       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
 466
 467       // This loop splits the track set into 2 random subsets
 468       if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
 469         float random = rn.Rndm();
 470         if(random < .5){
 471           mQx1 += (weight*cos(2*track->Phi()));
 472           mQy1 += (weight*sin(2*track->Phi()));
 473           if (fSaveTrackContribution){
 474             EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
 475             EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
 476           }
 477           trackcounter1++;
 478         }
 479         else {
 480           mQx2 += (weight*cos(2*track->Phi()));
 481           mQy2 += (weight*sin(2*track->Phi()));
 482           if (fSaveTrackContribution){
 483             EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
 484             EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
 485           }
 486           trackcounter2++;
 487         }
 488       }
 489       else if( trackcounter1 >= int(nt/2.)){
 490         mQx2 += (weight*cos(2*track->Phi()));
 491         mQy2 += (weight*sin(2*track->Phi()));
 492         if (fSaveTrackContribution){
 493           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
 494           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
 495         }
 496         trackcounter2++;
 497       }
 498       else {
 499         mQx1 += (weight*cos(2*track->Phi()));
 500         mQy1 += (weight*sin(2*track->Phi()));
 501         if (fSaveTrackContribution){
 502           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
 503           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
 504         }
 505         trackcounter1++;
 506       }
 507     }
 508   } else if (fSplitMethod == AliEPSelectionTask::kEta) {
 509
 510     for (Int_t i = 0; i < nt; i++) {
 511       weight = 1;
 512       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
 513       if (!track) continue;
 514       weight = GetWeight(track);
 515       Double_t eta = track->Eta();
 516       idtemp = track->GetID();
 517       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
 518
 519       if (eta > fEtaGap/2.) {
 520         mQx1 += (weight*cos(2*track->Phi()));
 521         mQy1 += (weight*sin(2*track->Phi()));
 522         if (fSaveTrackContribution){
 523           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
 524           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
 525         }
 526       } else if (eta < -1.*fEtaGap/2.) {
 527         mQx2 += (weight*cos(2*track->Phi()));
 528         mQy2 += (weight*sin(2*track->Phi()));
 529         if (fSaveTrackContribution){
 530           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
 531           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
 532         }
 533       }
 534     }
 535   } else if (fSplitMethod == AliEPSelectionTask::kCharge) {
 536
 537     for (Int_t i = 0; i < nt; i++) {
 538       weight = 1;
 539       track = dynamic_cast<AliVTrack*> (tracklist->At(i));
 540       if (!track) continue;
 541       weight = GetWeight(track);
 542       Short_t cha = track->Charge();
 543       idtemp = track->GetID();
 544       if ((fAnalysisInput.CompareTo("AOD")==0) && (fAODfilterbit == 128)) idtemp = idtemp*(-1) - 1;
 545
 546       if (cha > 0) {
 547         mQx1 += (weight*cos(2*track->Phi()));
 548         mQy1 += (weight*sin(2*track->Phi()));
 549         if (fSaveTrackContribution){
 550           EP->GetQContributionXArraysub1()->AddAt(weight*cos(2*track->Phi()),idtemp);
 551           EP->GetQContributionYArraysub1()->AddAt(weight*sin(2*track->Phi()),idtemp);
 552         }
 553       } else if (cha < 0) {
 554         mQx2 += (weight*cos(2*track->Phi()));
 555         mQy2 += (weight*sin(2*track->Phi()));
 556         if (fSaveTrackContribution){
 557           EP->GetQContributionXArraysub2()->AddAt(weight*cos(2*track->Phi()),idtemp);
 558           EP->GetQContributionYArraysub2()->AddAt(weight*sin(2*track->Phi()),idtemp);
 559         }
 560       }
 561     }
 562   } else {
 563     printf("plane resolution determination method not available!\n\n ");
 564     return;
 565   }
 566
 567   mQ[0].Set(mQx1,mQy1);
 568   mQ[1].Set(mQx2,mQy2);
 569   Q1 = mQ[0];
 570   Q2 = mQ[1];
 571 }
 572
 573 //________________________________________________________________________
 574 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
 575
 576   if(fESDtrackCuts){
 577     delete fESDtrackCuts;
 578     fESDtrackCuts = 0;
 579   }
 580   if (fAnalysisInput.CompareTo("AOD")==0){
 581     AliInfo("ESD track cuts not possible for AOD analysis; please use SetPersonalAODtrackCuts(); using TPC only track cuts");
 582     fUsercuts = kFALSE;
 583     SetTrackType("TPC");
 584     return;
 585   }
 586   fUsercuts = kTRUE;
 587   fESDtrackCuts = trackcuts;
 588 }
 589
 590 //________________________________________________________________________
 591 void AliEPSelectionTask::SetPersonalAODtrackCuts(UInt_t filterbit, Float_t etalow, Float_t etaup, Float_t ptlow, Float_t ptup, Int_t ntpc){
 592
 593   if(fESDtrackCuts){
 594     delete fESDtrackCuts;
 595     fESDtrackCuts = 0;
 596   }
 597   if (fAnalysisInput.CompareTo("ESD")==0){
 598     AliInfo("AOD track cuts not possible for ESD analysis; please use SetPersonalESDtrackCuts(); using TPC only track cuts");
 599     fUsercuts = kFALSE;
 600     SetTrackType("TPC");
 601     return;
 602   }
 603   fUsercuts = kTRUE;
 604   fESDtrackCuts = new AliESDtrackCuts();
 605   fESDtrackCuts->SetPtRange(ptlow,ptup);
 606   fESDtrackCuts->SetMinNClustersTPC(ntpc);
 607   fESDtrackCuts->SetEtaRange(etalow,etaup);
 608   fAODfilterbit = filterbit;
 609 }
 610
 611 //_____________________________________________________________________________
 612
 613 void AliEPSelectionTask::SetTrackType(TString tracktype){
 614   fTrackType = tracktype;
 615   if (!fUsercuts) {
 616   if (fTrackType.CompareTo("GLOBAL")==0){
 617     fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
 618     fAODfilterbit = 32;
 619   }
 620   if (fTrackType.CompareTo("TPC")==0){
 621     fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 622     fAODfilterbit = 128;
 623   }
 624   fESDtrackCuts->SetPtRange(0.15,20.);
 625   fESDtrackCuts->SetEtaRange(-0.8,0.8);
 626   }
 627 }
 628
 629 //________________________________________________________________________
 630 Double_t AliEPSelectionTask::GetWeight(TObject* track1)
 631 {
 632   Double_t ptweight=1;
 633   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
 634   if (fUsePtWeight && track) {
 635     if (track->Pt()<2) ptweight=track->Pt();
 636     else ptweight=2;
 637   }
 638   return ptweight*GetPhiWeight(track);
 639 }
 640
 641 //________________________________________________________________________
 642 Double_t AliEPSelectionTask::GetPhiWeight(TObject* track1)
 643 {
 644   Double_t phiweight=1;
 645   AliVTrack* track = dynamic_cast<AliVTrack*>(track1);
 646
 647   TH1F *phiDist = 0x0;
 648   if(track) phiDist = SelectPhiDist(track);
 649
 650   if (fUsePhiWeight && phiDist && track) {
 651     Double_t nParticles = phiDist->Integral();
 652     Double_t nPhibins = phiDist->GetNbinsX();
 653
 654     Double_t Phi = track->Phi();
 655
 656     while (Phi<0) Phi += TMath::TwoPi();
 657     while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
 658
 659     Double_t PhiDistValue = phiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
 660
 661     if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
 662   }
 663   return phiweight;
 664 }
 665
 666 //__________________________________________________________________________
 667 void AliEPSelectionTask::SetPhiDist()
 668 {
 669   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
 670
 671     if (fPeriod.CompareTo("LHC10h")==0)
 672        {
 673         fPhiDist[0] = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");}
 674         else if(fPeriod.CompareTo("LHC11h")==0){
 675             Int_t runbin=fHruns->FindBin(fRunNumber);
 676             if (fHruns->GetBinContent(runbin) > 1){
 677               fSparseDist->GetAxis(0)->SetRange(runbin,runbin);
 678               }
 679             else if(fHruns->GetBinContent(runbin) < 2){
 680                fSparseDist->GetAxis(0)->SetRange(1,2901); // not calibrated run, use integrated phi-weights
 681                AliInfo("Using integrated Phi-weights for this run");
 682                }
 683             for (Int_t i = 0; i<4 ;i++)
 684             {
 685              if(fPhiDist[i]){
 686                delete fPhiDist[i];
 687                fPhiDist[i] = 0x0;
 688                }
 689              if(i == 0){
 690                fSparseDist->GetAxis(1)->SetRange(1,1);  // neg charge
 691                fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
 692              if(i == 1){
 693                fSparseDist->GetAxis(1)->SetRange(2,2);  // pos charge
 694                fSparseDist->GetAxis(2)->SetRange(1,1);} // neg eta
 695              if(i == 2){
 696                fSparseDist->GetAxis(1)->SetRange(1,1);  // neg charge
 697                fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
 698              if(i == 3){
 699                fSparseDist->GetAxis(1)->SetRange(2,2);  // pos charge
 700                fSparseDist->GetAxis(2)->SetRange(2,2);} // pos eta
 701              fPhiDist[i] = (TH1F*)fSparseDist->Projection(3); // Projection on Phi
 702              fPhiDist[i]->SetName(Form("phidist%d%d",i,fRunNumber));
 703              fSparseDist->GetAxis(1)->SetRange(1,2); // reset axes
 704              fSparseDist->GetAxis(2)->SetRange(1,2);
 705              }
 706              fSparseDist->GetAxis(0)->SetRange(1,2901);// reset run axis
 707        }
 708
 709     if (!fPhiDist[0]) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
 710
 711   }
 712
 713
 714   if (fPeriod.CompareTo("LHC10h")==0 || fUserphidist){
 715      Bool_t emptybins;
 716
 717      int iter = 0;
 718      while (iter<3){
 719          emptybins = kFALSE;
 720
 721          for (int i=1; i<fPhiDist[0]->GetNbinsX(); i++){
 722            if (!((fPhiDist[0]->GetBinContent(i))>0)) {
 723              emptybins = kTRUE;
 724            }
 725          }
 726          if (emptybins) {
 727            cout << "empty bins - rebinning!" << endl;
 728            fPhiDist[0]->Rebin();
 729            iter++;
 730          }
 731          else iter = 3;
 732      }
 733
 734      if (emptybins) {
 735        AliError("After Maximum of rebinning still empty Phi-bins!!!");
 736      }
 737   }
 738   if (fPeriod.CompareTo("LHC10h") != 0 && fPeriod.CompareTo("LHC11h") != 0 && !fUserphidist){
 739   AliInfo("No Phi-weights available. All Phi weights set to 1");
 740   SetUsePhiWeight(kFALSE);
 741   }
 742 }
 743
 744 //__________________________________________________________________________
 745 void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
 746 {
 747
 748   fUserphidist = kTRUE;
 749
 750   TFile f(infilename);
 751   TObject* list = f.Get(listname);
 752   fPhiDist[0] = (TH1F*)list->FindObject("fHOutPhi");
 753   if (!fPhiDist[0]) AliFatal("Phi Distribution not found!!!");
 754
 755   f.Close();
 756 }
 757
 758
 759 //_________________________________________________________________________
 760 TObjArray* AliEPSelectionTask::GetAODTracksAndMaxID(AliAODEvent* aod, Int_t& maxid)
 761 {
 762   TObjArray *acctracks = new TObjArray();
 763
 764   AliAODTrack *tr = 0;
 765   Int_t maxid1 = 0;
 766   Int_t maxidtemp = -1;
 767   Float_t ptlow = 0;
 768   Float_t ptup = 0;
 769   Float_t etalow = 0;
 770   Float_t etaup = 0;
 771   fESDtrackCuts->GetPtRange(ptlow,ptup);
 772   fESDtrackCuts->GetEtaRange(etalow,etaup);
 773   Int_t ntpc = fESDtrackCuts->GetMinNClusterTPC();
 774
 775   for (Int_t i = 0; i < aod->GetNumberOfTracks() ; i++){
 776      tr = aod->GetTrack(i);
 777      maxidtemp = tr->GetID();
 778      if(maxidtemp < 0 && fAODfilterbit != 128) continue;
 779      if(maxidtemp > -1 && fAODfilterbit == 128) continue;
 780      if (fAODfilterbit == 128) maxidtemp = maxidtemp*(-1) - 1;
 781      if (maxidtemp > maxid1) maxid1 = maxidtemp;
 782      if(tr->TestFilterBit(fAODfilterbit) && tr->Pt() < ptup && tr->Pt() > ptlow && tr->Eta() < etaup && tr->Eta() > etalow && tr->GetTPCNcls() > ntpc){
 783      acctracks->Add(tr);
 784      }
 785   }
 786
 787   maxid = maxid1;
 788
 789   return acctracks;
 790
 791 }
 792
 793 //_________________________________________________________________________
 794 void AliEPSelectionTask::SetOADBandPeriod()
 795 {
 796   TString oadbfilename;
 797
 798   if (fRunNumber >= 136851 && fRunNumber <= 139515) // LHC10h
 799     {fPeriod = "LHC10h";
 800      if (!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
 801
 802         if (fAnalysisInput.CompareTo("AOD")==0){
 803            oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.aod.root", AliAnalysisManager::GetOADBPath()));
 804            } else if (fAnalysisInput.CompareTo("ESD")==0){
 805            oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
 806            }
 807
 808        TFile foadb(oadbfilename);
 809        if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
 810
 811        AliInfo("Using Standard OADB");
 812        fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
 813        if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
 814        foadb.Close();
 815        }
 816      }
 817
 818   if (fRunNumber >= 166529 && fRunNumber <= 170593) // LHC11h
 819      {fPeriod = "LHC11h";
 820       if (!fUserphidist) {
 821       // if it's already set and custom class is required, we use the one provided by the user
 822
 823       oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist2011.root", AliAnalysisManager::GetOADBPath()));
 824       TFile *foadb = TFile::Open(oadbfilename);
 825       if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
 826
 827       AliInfo("Using Standard OADB");
 828       fSparseDist = (THnSparse*) foadb->Get("Default");
 829       if (!fSparseDist) AliFatal("Cannot fetch OADB container for EP selection");
 830       foadb->Close();
 831       if(!fHruns){
 832            fHruns = (TH1F*)fSparseDist->Projection(0); //projection on run axis;
 833            fHruns->SetName("runsHisto");
 834       }
 835       }
 836      }
 837 }
 838
 839 //_________________________________________________________________________
 840 TH1F* AliEPSelectionTask::SelectPhiDist(AliVTrack *track)
 841 {
 842   if (fPeriod.CompareTo("LHC10h")==0  || fUserphidist) return fPhiDist[0];
 843   else if(fPeriod.CompareTo("LHC11h")==0)
 844     {
 845      if (track->Charge() < 0)
 846        {
 847         if(track->Eta() < 0.)       return fPhiDist[0];
 848         else if (track->Eta() > 0.) return fPhiDist[2];
 849        }
 850       else if (track->Charge() > 0)
 851        {
 852         if(track->Eta() < 0.)       return fPhiDist[1];
 853         else if (track->Eta() > 0.) return fPhiDist[3];
 854        }
 855
 856     }
 857   return 0;
 858 }
 859
 860 TObjArray* AliEPSelectionTask::GetTracksForLHC11h(AliESDEvent* esd)
 861 {
 862   // 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
 863   TObjArray *acctracks = new TObjArray();
 864   acctracks->SetOwner(kTRUE);
 865
 866   const AliESDVertex *vtxSPD = esd->GetPrimaryVertexSPD();
 867
 868   Float_t ptlow = 0;
 869   Float_t ptup = 0;
 870   Float_t etalow = 0;
 871   Float_t etaup = 0;
 872   fESDtrackCuts->GetPtRange(ptlow,ptup);
 873   fESDtrackCuts->GetEtaRange(etalow,etaup);
 874
 875   Double_t pDCA[3] = { 0. }; // momentum at DCA
 876   Double_t rDCA[3] = { 0. }; // position at DCA
 877   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
 878   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
 879
 880
 881
 882   for (Int_t i = 0; i < esd->GetNumberOfTracks() ; i++){
 883
 884     AliESDtrack* esdTrack = esd->GetTrack(i); //carefull do not modify it othwise  need to work with a copy
 885     //
 886     // Track selection
 887     if (!fESDtrackCuts->AcceptTrack(esdTrack)) continue;
 888
 889     // create a tpc only tracl
 890     AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(esd),esdTrack->GetID());
 891     if(!track) continue;
 892
 893     if(track->Pt()>0.)
 894     {
 895       // only constrain tracks above threshold
 896       AliExternalTrackParam exParam;
 897       // take the B-field from the ESD, no 3D fieldMap available at this point
 898       Bool_t relate = false;
 899       relate = track->RelateToVertexTPC(vtxSPD,esd->GetMagneticField(),kVeryBig,&exParam);
 900       if(!relate){
 901         delete track;
 902         continue;
 903       }
 904       // fetch the track parameters at the DCA (unconstraint)
 905       if(track->GetTPCInnerParam()){
 906         track->GetTPCInnerParam()->GetPxPyPz(pDCA);
 907         track->GetTPCInnerParam()->GetXYZ(rDCA);
 908       }
 909       // get the DCA to the vertex:
 910       track->GetImpactParametersTPC(dDCA,cDCA);
 911       // set the constrained parameters to the track
 912       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 913     }
 914
 915
 916     Float_t eta = track->Eta();
 917     Float_t pT = track->Pt();
 918
 919     if(pT < ptlow || pT > ptup || eta < etalow || eta > etaup){
 920       delete track;
 921       continue;
 922     }
 923
 924     acctracks->Add(track);
 925    }
 926
 927
 928   return acctracks;
 929
 930 }
 931