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