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
  59 #include "AliEventplane.h"
  60
  61 ClassImp(AliEPSelectionTask)
  62
  63 //________________________________________________________________________
  64 AliEPSelectionTask::AliEPSelectionTask():
  65 AliAnalysisTaskSE(),
  66   fDebug(0),
  67   fAnalysisInput("ESD"),
  68   fStatus("TPC"),
  69   fUseMCRP(kFALSE),
  70   fUsePhiWeight(kFALSE),
  71   fUsePtWeight(kFALSE),
  72   fSaveTrackContribution(kFALSE),
  73   fESDtrackCuts(0),
  74   ftracklist(0),
  75   fPhiDist(0),
  76   fQVector(0),
  77   fQContributionX(0),
  78   fQContributionY(0),
  79   fEventplaneQ(0),
  80   fQsub1(0),
  81   fQsub2(0),
  82   fQsubRes(0),
  83   fOutputList(0),
  84   fHOutEventplaneQ(0),
  85   fHOutPhi(0),
  86   fHOutPhiCorr(0),
  87   fHOutsub1sub2(0),
  88   fHOutNTEPRes(0),
  89   fHOutPTPsi(0),
  90   fHOutDiff(0),
  91   fHOutleadPTPsi(0)
  92 {
  93   // Default constructor
  94   AliInfo("Event Plane Selection enabled.");
  95 }
  96
  97 //________________________________________________________________________
  98 AliEPSelectionTask::AliEPSelectionTask(const char *name):
  99   AliAnalysisTaskSE(name),
 100   fDebug(0),
 101   fAnalysisInput("ESD"),
 102   fStatus("TPC"),
 103   fUseMCRP(kFALSE),
 104   fUsePhiWeight(kFALSE),
 105   fUsePtWeight(kFALSE),
 106   fSaveTrackContribution(kFALSE),
 107   fESDtrackCuts(0),
 108   ftracklist(0),
 109   fPhiDist(0),
 110   fQVector(0),
 111   fQContributionX(0),
 112   fQContributionY(0),
 113   fEventplaneQ(0),
 114   fQsub1(0),
 115   fQsub2(0),
 116   fQsubRes(0),
 117   fOutputList(0),
 118   fHOutEventplaneQ(0),
 119   fHOutPhi(0),
 120   fHOutPhiCorr(0),
 121   fHOutsub1sub2(0),
 122   fHOutNTEPRes(0),
 123   fHOutPTPsi(0),
 124   fHOutDiff(0),
 125   fHOutleadPTPsi(0)
 126 {
 127   // Default constructor
 128   AliInfo("Event Plane Selection enabled.");
 129   DefineOutput(1, TList::Class());
 130 }
 131
 132 //________________________________________________________________________
 133 AliEPSelectionTask::~AliEPSelectionTask()
 134 {
 135   // Destructor
 136   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
 137       delete fOutputList;
 138       fOutputList = 0;
 139   }
 140   if (fESDtrackCuts){
 141       delete fESDtrackCuts;
 142       fESDtrackCuts = 0;
 143   }
 144   if (fQVector){
 145       delete fQVector;
 146       fQVector = 0;
 147   }
 148   if (fQsub1){
 149       delete fQsub1;
 150       fQsub1 = 0;
 151   }
 152     if (fQsub2){
 153       delete fQsub2;
 154       fQsub2 = 0;
 155   }
 156 }
 157
 158 //________________________________________________________________________
 159 void AliEPSelectionTask::UserCreateOutputObjects()
 160 {
 161   // Create the output containers
 162   if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
 163   AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
 164
 165   fOutputList = new TList();
 166   fOutputList->SetOwner();
 167   fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
 168   fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
 169   fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
 170   fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
 171   fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
 172   fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
 173   fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
 174   fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
 175
 176   fOutputList->Add(fHOutEventplaneQ);
 177   fOutputList->Add(fHOutPhi);
 178   fOutputList->Add(fHOutPhiCorr);
 179   fOutputList->Add(fHOutsub1sub2);
 180   fOutputList->Add(fHOutNTEPRes);
 181   fOutputList->Add(fHOutPTPsi);
 182   fOutputList->Add(fHOutleadPTPsi);
 183   fOutputList->Add(fHOutDiff);
 184
 185   PostData(1, fOutputList);
 186 }
 187
 188 //________________________________________________________________________
 189 void AliEPSelectionTask::UserExec(Option_t */*option*/)
 190 {
 191   // Execute analysis for current event:
 192   if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
 193
 194   AliEventplane* esdEP = 0;
 195   TVector2 QQ1;
 196   TVector2 QQ2;
 197   Double_t fRP = 0.; // the monte carlo reaction plane angle
 198
 199   if (fAnalysisInput.CompareTo("ESD")==0){
 200
 201     AliVEvent* event = InputEvent();
 202     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
 203
 204     if (fUseMCRP) {
 205       AliMCEvent* mcEvent  = MCEvent();
 206       if (mcEvent && mcEvent->GenEventHeader()) {
 207         AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
 208         if (headerH) fRP = headerH->ReactionPlaneAngle();
 209       }
 210     }
 211     if (esd){
 212       esdEP = esd->GetEventplane();
 213       if (fSaveTrackContribution) {
 214         esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
 215         esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
 216       }
 217
 218       if (fStatus.CompareTo("GLOBAL")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
 219       if (fStatus.CompareTo("TPC")==0) ftracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
 220       const int NT = ftracklist->GetEntries();
 221
 222       if (NT>4){
 223         fQVector = new TVector2(GetQ(esdEP,ftracklist));
 224         fEventplaneQ = fQVector->Phi()/2;
 225         GetQsub(QQ1, QQ2, ftracklist);
 226         fQsub1 = new TVector2(QQ1);
 227         fQsub2 = new TVector2(QQ2);
 228         fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
 229         esdEP->SetQVector(fQVector);
 230         esdEP->SetEventplaneQ(fEventplaneQ);
 231         esdEP->SetQsub(fQsub1,fQsub2);
 232         esdEP->SetQsubRes(fQsubRes);
 233
 234         fHOutEventplaneQ->Fill(fEventplaneQ);
 235         fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
 236         fHOutNTEPRes->Fill(NT,fQsubRes);
 237
 238         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
 239
 240         for (int iter = 0; iter<NT;iter++){
 241           AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
 242           if (track) {
 243             float delta = track->Phi()-fEventplaneQ;
 244             while (delta < 0) delta += TMath::Pi();
 245             while (delta > TMath::Pi()) delta -= TMath::Pi();
 246             fHOutPTPsi->Fill(track->Pt(),delta);
 247             fHOutPhi->Fill(track->Phi());
 248             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
 249           }
 250         }
 251
 252         AliESDtrack* trmax = esd->GetTrack(0);
 253         for (int iter = 1; iter<NT;iter++){
 254           AliESDtrack* track = dynamic_cast<AliESDtrack*> (ftracklist->At(iter));
 255           if (track && (track->Pt() > trmax->Pt())) trmax = track;
 256         }
 257         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
 258       }
 259       delete ftracklist;
 260     }
 261   }
 262
 263   else if (fAnalysisInput.CompareTo("AOD")==0){
 264     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
 265     // to be implemented
 266     printf("  AOD analysis not yet implemented!!!\n\n");
 267     return;
 268   }
 269   else {
 270     printf(" Analysis Input not known!\n\n ");
 271     return;
 272   }
 273   PostData(1, fOutputList);
 274 }
 275
 276 //________________________________________________________________________
 277 void AliEPSelectionTask::Terminate(Option_t */*option*/)
 278 {
 279   // Terminate analysis
 280 }
 281
 282 //__________________________________________________________________________
 283 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
 284 {
 285   TVector2 mQ;
 286   float mQx=0, mQy=0;
 287   AliESDtrack* track;
 288   Double_t weight;
 289
 290   int NT = tracklist->GetEntries();
 291
 292   for (int i=0; i<NT; i++){
 293     weight = 1;
 294     track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
 295     if (track) {
 296       weight = GetWeight(track);
 297       if (fSaveTrackContribution){
 298         EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),track->GetID());
 299         EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),track->GetID());
 300       }
 301       mQx += (weight*cos(2*track->Phi()));
 302       mQy += (weight*sin(2*track->Phi()));
 303     }
 304   }
 305   mQ.Set(mQx,mQy);
 306   return mQ;
 307 }
 308
 309   //________________________________________________________________________
 310 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
 311 {
 312   TVector2 mQ[2];
 313   float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
 314   Double_t weight;
 315
 316   AliESDtrack* track;
 317   TRandom2 RN = 0;
 318
 319   int NT = tracklist->GetEntries();
 320   int trackcounter1=0, trackcounter2=0;
 321
 322   for (Int_t i = 0; i < NT; i++) {
 323     weight = 1;
 324     track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
 325     if (!track) continue;
 326     weight = GetWeight(track);
 327
 328     // This loop splits the track set into 2 random subsets
 329     if( trackcounter1 < int(NT/2.) && trackcounter2 < int(NT/2.)){
 330       float random = RN.Rndm();
 331       if(random < .5){
 332         mQx1 += (weight*cos(2*track->Phi()));
 333         mQy1 += (weight*sin(2*track->Phi()));
 334         trackcounter1++;
 335       }
 336       else {
 337         mQx2 += (weight*cos(2*track->Phi()));
 338         mQy2 += (weight*sin(2*track->Phi()));
 339         trackcounter2++;
 340       }
 341     }
 342     else if( trackcounter1 >= int(NT/2.)){
 343       mQx2 += (weight*cos(2*track->Phi()));
 344       mQy2 += (weight*sin(2*track->Phi()));
 345       trackcounter2++;
 346     }
 347     else {
 348       mQx1 += (weight*cos(2*track->Phi()));
 349       mQy1 += (weight*sin(2*track->Phi()));
 350       trackcounter1++;
 351     }
 352   }
 353   mQ[0].Set(mQx1,mQy1);
 354   mQ[1].Set(mQx2,mQy2);
 355   Q1 = mQ[0];
 356   Q2 = mQ[1];
 357 }
 358
 359 //________________________________________________________________________
 360 void AliEPSelectionTask::SetESDtrackCuts(TString status){
 361
 362   fStatus = status;
 363   if (fStatus.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
 364   if (fStatus.CompareTo("TPC")==0)    fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 365   fESDtrackCuts->SetPtRange(0.15,20);
 366   fESDtrackCuts->SetEtaRange(-0.8,0.8);
 367 }
 368
 369 //__________________________________________________________________________
 370 void AliEPSelectionTask::SetPhiDistribution(char* infilename, char* listname)
 371 {
 372   TFile f(infilename);
 373   TObject* list = f.Get(listname);
 374   fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
 375   if (!fPhiDist) {
 376     cout << "Phi Distribution not found!!!" << endl;
 377     return;
 378   }
 379
 380   Bool_t emptybins;
 381
 382   int iter = 0;
 383   while (iter<3){
 384       emptybins = kFALSE;
 385
 386       for (int i=1; i<fPhiDist->GetNbinsX(); i++){
 387         if (!((fPhiDist->GetBinContent(i))>0)) {
 388           emptybins = kTRUE;
 389         }
 390       }
 391       if (emptybins) {
 392         cout << "empty bins - rebinning!" << endl;
 393         fPhiDist->Rebin();
 394         iter++;
 395       }
 396       else iter = 3;
 397   }
 398
 399   if (emptybins) {
 400     AliError("After Maximum of rebinning still empty Phi-bins!!!");
 401   }
 402   f.Close();
 403 }
 404
 405 //________________________________________________________________________
 406 Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
 407 {
 408   Double_t ptweight=1;
 409
 410   if (fUsePtWeight) {
 411     if (track->Pt()<2) ptweight=track->Pt();
 412     else ptweight=2;
 413   }
 414   return ptweight*GetPhiWeight(track);
 415 }
 416
 417 //________________________________________________________________________
 418 Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
 419 {
 420   Double_t phiweight=1;
 421
 422   if (fUsePhiWeight && fPhiDist && track) {
 423     Double_t nParticles = fPhiDist->Integral();
 424     Double_t nPhibins = fPhiDist->GetNbinsX();
 425
 426     Double_t Phi = track->Phi();
 427
 428     while (Phi<0) Phi += TMath::TwoPi();
 429     while (Phi>TMath::TwoPi()) Phi -= TMath::TwoPi();
 430
 431     Double_t PhiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
 432
 433     if (PhiDistValue > 0) phiweight = nParticles/nPhibins/PhiDistValue;
 434   }
 435   return phiweight;
 436 }