ANALYSIS/AliEPSelectionTask.cxx

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