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 <TProfile.h>
  29 #include <TFile.h>
  30 #include <TObjString.h>
  31 #include <TString.h>
  32 #include <TCanvas.h>
  33 #include <TROOT.h>
  34 #include <TDirectory.h>
  35 #include <TSystem.h>
  36 #include <iostream>
  37 #include <TRandom2.h>
  38 #include <TArrayF.h>
  39
  40 #include "AliAnalysisManager.h"
  41 #include "AliVEvent.h"
  42 #include "AliESD.h"
  43 #include "AliESDEvent.h"
  44 #include "AliMCEvent.h"
  45 #include "AliESDtrack.h"
  46 #include "AliESDtrackCuts.h"
  47 #include "AliESDHeader.h"
  48 #include "AliESDInputHandler.h"
  49 #include "AliAODHandler.h"
  50 #include "AliAODEvent.h"
  51 #include "AliHeader.h"
  52 #include "AliGenHijingEventHeader.h"
  53 #include "AliAnalysisTaskSE.h"
  54 #include "AliPhysicsSelectionTask.h"
  55 #include "AliPhysicsSelection.h"
  56 #include "AliBackgroundSelection.h"
  57 #include "AliESDUtils.h"
  58 #include "AliOADBContainer.h"
  59
  60 #include "AliEventplane.h"
  61
  62 ClassImp(AliEPSelectionTask)
  63
  64 //________________________________________________________________________
  65 AliEPSelectionTask::AliEPSelectionTask():
  66 AliAnalysisTaskSE(),
  67   fDebug(0),
  68   fAnalysisInput("ESD"),
  69   fTrackType("TPC"),
  70   fUseMCRP(kFALSE),
  71   fUsePhiWeight(kFALSE),
  72   fUsePtWeight(kFALSE),
  73   fSaveTrackContribution(kFALSE),
  74   fuserphidist(kFALSE),
  75   fusercuts(kFALSE),
  76   frunNumber(-15),
  77   fESDtrackCuts(0),
  78   ftracklist(0),
  79   fEPContainer(0),
  80   fPhiDist(0),
  81   fQVector(0),
  82   fQContributionX(0),
  83   fQContributionY(0),
  84   fEventplaneQ(0),
  85   fQsub1(0),
  86   fQsub2(0),
  87   fQsubRes(0),
  88   fOutputList(0),
  89   fHOutEventplaneQ(0),
  90   fHOutPhi(0),
  91   fHOutPhiCorr(0),
  92   fHOutsub1sub2(0),
  93   fHOutNTEPRes(0),
  94   fHOutPTPsi(0),
  95   fHOutDiff(0),
  96   fHOutleadPTPsi(0)
  97 {
  98   // Default constructor
  99   AliInfo("Event Plane Selection enabled.");
 100 }
 101
 102 //________________________________________________________________________
 103 AliEPSelectionTask::AliEPSelectionTask(const char *name):
 104   AliAnalysisTaskSE(name),
 105   fDebug(0),
 106   fAnalysisInput("ESD"),
 107   fTrackType("TPC"),
 108   fUseMCRP(kFALSE),
 109   fUsePhiWeight(kFALSE),
 110   fUsePtWeight(kFALSE),
 111   fSaveTrackContribution(kFALSE),
 112   fuserphidist(kFALSE),
 113   fusercuts(kFALSE),
 114   frunNumber(-15),
 115   fESDtrackCuts(0),
 116   ftracklist(0),
 117   fEPContainer(0),
 118   fPhiDist(0),
 119   fQVector(0),
 120   fQContributionX(0),
 121   fQContributionY(0),
 122   fEventplaneQ(0),
 123   fQsub1(0),
 124   fQsub2(0),
 125   fQsubRes(0),
 126   fOutputList(0),
 127   fHOutEventplaneQ(0),
 128   fHOutPhi(0),
 129   fHOutPhiCorr(0),
 130   fHOutsub1sub2(0),
 131   fHOutNTEPRes(0),
 132   fHOutPTPsi(0),
 133   fHOutDiff(0),
 134   fHOutleadPTPsi(0)
 135 {
 136   // Default constructor
 137   AliInfo("Event Plane Selection enabled.");
 138   DefineOutput(1, TList::Class());
 139 }
 140
 141 //________________________________________________________________________
 142 AliEPSelectionTask::~AliEPSelectionTask()
 143 {
 144   // Destructor
 145   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
 146       delete fOutputList;
 147       fOutputList = 0;
 148   }
 149   if (fESDtrackCuts){
 150       delete fESDtrackCuts;
 151       fESDtrackCuts = 0;
 152   }
 153   if (fQVector){
 154       delete fQVector;
 155       fQVector = 0;
 156   }
 157   if (fQsub1){
 158       delete fQsub1;
 159       fQsub1 = 0;
 160   }
 161   if (fQsub2){
 162       delete fQsub2;
 163       fQsub2 = 0;
 164   }
 165   if (fPhiDist){
 166       delete fPhiDist;
 167       fPhiDist = 0;
 168   }
 169   if (ftracklist){
 170       delete ftracklist;
 171       ftracklist = 0;
 172   }
 173   if (fEPContainer){
 174       delete fEPContainer;
 175       fEPContainer = 0;
 176   }
 177 }
 178
 179 //________________________________________________________________________
 180 void AliEPSelectionTask::UserCreateOutputObjects()
 181 {
 182   // Create the output containers
 183   if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
 184   AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
 185
 186   fOutputList = new TList();
 187   fOutputList->SetOwner();
 188   fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
 189   fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
 190   fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
 191   fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
 192   fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
 193   fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
 194   fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
 195   fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
 196
 197   fOutputList->Add(fHOutEventplaneQ);
 198   fOutputList->Add(fHOutPhi);
 199   fOutputList->Add(fHOutPhiCorr);
 200   fOutputList->Add(fHOutsub1sub2);
 201   fOutputList->Add(fHOutNTEPRes);
 202   fOutputList->Add(fHOutPTPsi);
 203   fOutputList->Add(fHOutleadPTPsi);
 204   fOutputList->Add(fHOutDiff);
 205
 206   PostData(1, fOutputList);
 207
 208     if(!fuserphidist) { // if it's already set and custom class is required, we use the one provided by the user
 209
 210     TString oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
 211
 212     TFile foadb(oadbfilename);
 213     if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
 214
 215     AliInfo("Using Standard OADB");
 216     fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
 217     if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
 218     foadb.Close();
 219     }
 220 }
 221
 222 //________________________________________________________________________
 223 void AliEPSelectionTask::UserExec(Option_t */*option*/)
 224 {
 225   // Execute analysis for current event:
 226   if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
 227
 228 //   frunNumber = -15;
 229
 230   AliEventplane* esdEP = 0;
 231   TVector2 QQ1;
 232   TVector2 QQ2;
 233   Double_t fRP = 0.; // the monte carlo reaction plane angle
 234
 235   if (fAnalysisInput.CompareTo("ESD")==0){
 236
 237     AliVEvent* event = InputEvent();
 238     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 239     if (esd){
 240       if (!(frunNumber == esd->GetRunNumber())) {
 241           frunNumber = esd->GetRunNumber();
 242             SetPhiDist();
 243       }
 244
 245
 246       if (fUseMCRP) {
 247           AliMCEvent* mcEvent  = MCEvent();
 248           if (mcEvent && mcEvent->GenEventHeader()) {
 249               AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
 250               if (headerH) fRP = headerH->ReactionPlaneAngle();
 251           }
 252       }
 253
 254       esdEP = esd->GetEventplane();
 255       if (fSaveTrackContribution) {
 256         esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
 257         esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
 258       }
 259
 260       if (fTrackType.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
 261       if (fTrackType.CompareTo("TPC")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
 262       const int NT = ftracklist->GetEntries();
 263
 264       if (NT>4){
 265         fQVector = new TVector2(GetQ(esdEP,ftracklist));
 266         fEventplaneQ = fQVector->Phi()/2;
 267         GetQsub(QQ1, QQ2, ftracklist);
 268         fQsub1 = new TVector2(QQ1);
 269         fQsub2 = new TVector2(QQ2);
 270         fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
 271
 272         esdEP->SetQVector(fQVector);
 273         esdEP->SetEventplaneQ(fEventplaneQ);
 274         esdEP->SetQsub(fQsub1,fQsub2);
 275         esdEP->SetQsubRes(fQsubRes);
 276
 277         fHOutEventplaneQ->Fill(fEventplaneQ);
 278         fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
 279         fHOutNTEPRes->Fill(NT,fQsubRes);
 280
 281         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
 282
 283         for (int iter = 0; iter<NT;iter++){
 284           AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
 285           if (track) {
 286             float delta = track->Phi()-fEventplaneQ;
 287             while (delta < 0) delta += TMath::Pi();
 288             while (delta > TMath::Pi()) delta -= TMath::Pi();
 289             fHOutPTPsi->Fill(track->Pt(),delta);
 290             fHOutPhi->Fill(track->Phi());
 291             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
 292           }
 293         }
 294
 295         AliESDtrack* trmax = esd->GetTrack(0);
 296         for (int iter = 1; iter<NT;iter++){
 297           AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
 298           if (track && (track->Pt() > trmax->Pt())) trmax = track;
 299         }
 300         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
 301       }
 302       delete ftracklist;
 303     }
 304   }
 305
 306   else if (fAnalysisInput.CompareTo("AOD")==0){
 307     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
 308     // to be implemented
 309     printf("  AOD analysis not yet implemented!!!\n\n");
 310     return;
 311   }
 312   else {
 313     printf(" Analysis Input not known!\n\n ");
 314     return;
 315   }
 316   PostData(1, fOutputList);
 317 }
 318
 319 //________________________________________________________________________
 320 void AliEPSelectionTask::Terminate(Option_t */*option*/)
 321 {
 322   // Terminate analysis
 323 }
 324
 325 //__________________________________________________________________________
 326 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
 327 {
 328   TVector2 mQ;
 329   float mQx=0, mQy=0;
 330   AliESDtrack* track;
 331   Double_t weight;
 332
 333   int NT = tracklist->GetEntries();
 334
 335   for (int i=0; i<NT; i++){
 336     weight = 1;
 337     track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
 338     if (track) {
 339       weight = GetWeight(track);
 340       if (fSaveTrackContribution){
 341         EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),track->GetID());
 342         EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),track->GetID());
 343       }
 344       mQx += (weight*cos(2*track->Phi()));
 345       mQy += (weight*sin(2*track->Phi()));
 346     }
 347   }
 348   mQ.Set(mQx,mQy);
 349   return mQ;
 350 }
 351
 352   //________________________________________________________________________
 353 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
 354 {
 355   TVector2 mQ[2];
 356   float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
 357   Double_t weight;
 358
 359   AliESDtrack* track;
 360   TRandom2 RN = 0;
 361
 362   int NT = tracklist->GetEntries();
 363   int trackcounter1=0, trackcounter2=0;
 364
 365   for (Int_t i = 0; i < NT; i++) {
 366     weight = 1;
 367     track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
 368     if (!track) continue;
 369     weight = GetWeight(track);
 370
 371     // This loop splits the track set into 2 random subsets
 372     if( trackcounter1 < int(NT/2.) && trackcounter2 < int(NT/2.)){
 373       float random = RN.Rndm();
 374       if(random < .5){
 375         mQx1 += (weight*cos(2*track->Phi()));
 376         mQy1 += (weight*sin(2*track->Phi()));
 377         trackcounter1++;
 378       }
 379       else {
 380         mQx2 += (weight*cos(2*track->Phi()));
 381         mQy2 += (weight*sin(2*track->Phi()));
 382         trackcounter2++;
 383       }
 384     }
 385     else if( trackcounter1 >= int(NT/2.)){
 386       mQx2 += (weight*cos(2*track->Phi()));
 387       mQy2 += (weight*sin(2*track->Phi()));
 388       trackcounter2++;
 389     }
 390     else {
 391       mQx1 += (weight*cos(2*track->Phi()));
 392       mQy1 += (weight*sin(2*track->Phi()));
 393       trackcounter1++;
 394     }
 395   }
 396   mQ[0].Set(mQx1,mQy1);
 397   mQ[1].Set(mQx2,mQy2);
 398   Q1 = mQ[0];
 399   Q2 = mQ[1];
 400 }
 401
 402 //________________________________________________________________________
 403 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
 404
 405   fusercuts = kTRUE;
 406   fESDtrackCuts = trackcuts;
 407 }
 408
 409 //_____________________________________________________________________________
 410 void AliEPSelectionTask::SetTrackType(TString tracktype){
 411   fTrackType = tracktype;
 412   if (!fusercuts) {
 413   if (fTrackType.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
 414   if (fTrackType.CompareTo("TPC")==0)    fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 415   fESDtrackCuts->SetPtRange(0.15,20);
 416   fESDtrackCuts->SetEtaRange(-0.8,0.8);
 417   }
 418 }
 419
 420 //________________________________________________________________________
 421 Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
 422 {
 423   Double_t ptweight=1;
 424
 425   if (fUsePtWeight) {
 426     if (track->Pt()<2) ptweight=track->Pt();
 427     else ptweight=2;
 428   }
 429   return ptweight*GetPhiWeight(track);
 430 }
 431
 432 //________________________________________________________________________
 433 Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
 434 {
 435   Double_t phiweight=1;
 436
 437   if (fUsePhiWeight && fPhiDist && track) {
 438     Double_t nParticles = fPhiDist->Integral();
 439     Double_t nPhibins = fPhiDist->GetNbinsX();
 440
 441     Double_t Phi = track->Phi();
 442
 443     while (Phi<0) Phi += TMath::TwoPi();
 444     while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
 445
 446     Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
 447
 448     if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
 449   }
 450   return phiweight;
 451 }
 452
 453 //__________________________________________________________________________
 454 void AliEPSelectionTask::SetPhiDist()
 455 {
 456   if(!fuserphidist) { // if it's already set and custom class is required, we use the one provided by the user
 457
 458     fPhiDist = (TH1F*) fEPContainer->GetObject(frunNumber, "Default");
 459     if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", frunNumber));
 460
 461   }
 462   else {
 463     AliInfo("Using Custom Phi Distribution");
 464   }
 465
 466   Bool_t emptybins;
 467
 468   int iter = 0;
 469   while (iter<3){
 470       emptybins = kFALSE;
 471
 472       for (int i=1; i<fPhiDist->GetNbinsX(); i++){
 473         if (!((fPhiDist->GetBinContent(i))>0)) {
 474           emptybins = kTRUE;
 475         }
 476       }
 477       if (emptybins) {
 478         cout << "empty bins - rebinning!" << endl;
 479         fPhiDist->Rebin();
 480         iter++;
 481       }
 482       else iter = 3;
 483   }
 484
 485   if (emptybins) {
 486     AliError("After Maximum of rebinning still empty Phi-bins!!!");
 487   }
 488 }
 489
 490 //__________________________________________________________________________
 491 void AliEPSelectionTask::SetPersonalPhiDistribution(char* infilename, char* listname)
 492 {
 493
 494   fuserphidist = kTRUE;
 495
 496   TFile f(infilename);
 497   TObject* list = f.Get(listname);
 498   fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
 499   if (!fPhiDist) AliFatal("Phi Distribution not found!!!");
 500
 501   f.Close();
 502 }