PWGLF/FORWARD/analysis2/AliFMDEventPlaneFinder.cxx

   1 // This task finds the eventplane
   2 // using the FMD
   3 //
   4 #include "AliFMDEventPlaneFinder.h"
   5 #include <TList.h>
   6 #include <TMath.h>
   7 #include "AliLog.h"
   8 #include <TH2D.h>
   9 #include <TH3D.h>
  10 #include "TROOT.h"
  11 #include <iostream>
  12 #include <iomanip>
  13 #include <TString.h>
  14 #include <TFile.h>
  15 #include "AliOADBContainer.h"
  16 #include "AliAnalysisManager.h"
  17 #include "AliVEvent.h"
  18 #include "AliEventplane.h"
  19 #include "AliCentrality.h"
  20 #include "AliAODForwardEP.h"
  21 #include "AliAODForwardMult.h"
  22 #include "AliAODEvent.h"
  23
  24 ClassImp(AliFMDEventPlaneFinder)
  25 #if 0
  26 ; // For Emacs
  27 #endif
  28
  29 //____________________________________________________________________
  30 AliFMDEventPlaneFinder::AliFMDEventPlaneFinder()
  31   : TNamed(),
  32     fList(0),
  33     fEvent(0),
  34     fQt(),
  35     fQa(),
  36     fQc(),
  37     fQ1(),
  38     fQ2(),
  39     fQeta(),
  40     fHepFMD(0),
  41     fHepFMDA(0),
  42     fHepFMDC(0),
  43     fHepFMDQC1(0),
  44     fHepFMDQC2(0),
  45     fHdiffFMDAC(0),
  46     fHdiffFMDTPC(0),
  47     fHdiffFMDVZERO(0),
  48     fHcorrFMDAC(0),
  49     fHcorrFMDTPC(0),
  50     fHcorrFMDVZERO(0),
  51     fHPhi(0),
  52     fDebug(0),
  53     fOADBFileName(0),
  54     fOADBContainer(0),
  55     fPhiDist(0),
  56     fRunNumber(0),
  57     fUsePhiWeights(0)
  58 {
  59   //
  60   // Constructor
  61   //
  62   DGUARD(fDebug,0,"Default CTOR of AliFMDEventPlaneFinder");
  63
  64 }
  65
  66 //____________________________________________________________________
  67 AliFMDEventPlaneFinder::AliFMDEventPlaneFinder(const char* title)
  68   : TNamed("fmdEventPlaneFinder", title),
  69     fList(0),
  70     fEvent(0),
  71     fQt(),
  72     fQa(),
  73     fQc(),
  74     fQ1(),
  75     fQ2(),
  76     fQeta(),
  77     fHepFMD(0),
  78     fHepFMDA(0),
  79     fHepFMDC(0),
  80     fHepFMDQC1(0),
  81     fHepFMDQC2(0),
  82     fHdiffFMDAC(0),
  83     fHdiffFMDTPC(0),
  84     fHdiffFMDVZERO(0),
  85     fHcorrFMDAC(0),
  86     fHcorrFMDTPC(0),
  87     fHcorrFMDVZERO(0),
  88     fHPhi(0),
  89     fDebug(0),
  90     fOADBFileName(0),
  91     fOADBContainer(0),
  92     fPhiDist(0),
  93     fRunNumber(0),
  94     fUsePhiWeights(1)
  95 {
  96   //
  97   // Constructor
  98   //
  99   // Parameters:
 100   //    name Name of object
 101   //
 102   DGUARD(fDebug,0,"Named CTOR of AliFMDEventPlaneFinder: %s", title);
 103 }
 104
 105 //____________________________________________________________________
 106 AliFMDEventPlaneFinder::AliFMDEventPlaneFinder(const
 107                                                  AliFMDEventPlaneFinder& o)
 108   : TNamed(o),
 109     fList(o.fList),
 110     fEvent(o.fEvent),
 111     fQt(o.fQt),
 112     fQa(o.fQa),
 113     fQc(o.fQc),
 114     fQ1(o.fQ1),
 115     fQ2(o.fQ2),
 116     fQeta(o.fQeta),
 117     fHepFMD(o.fHepFMD),
 118     fHepFMDA(o.fHepFMDA),
 119     fHepFMDC(o.fHepFMDC),
 120     fHepFMDQC1(o.fHepFMDQC1),
 121     fHepFMDQC2(o.fHepFMDQC2),
 122     fHdiffFMDAC(o.fHdiffFMDAC),
 123     fHdiffFMDTPC(o.fHdiffFMDTPC),
 124     fHdiffFMDVZERO(o.fHdiffFMDVZERO),
 125     fHcorrFMDAC(o.fHcorrFMDAC),
 126     fHcorrFMDTPC(o.fHcorrFMDTPC),
 127     fHcorrFMDVZERO(o.fHcorrFMDVZERO),
 128     fHPhi(o.fHPhi),
 129     fDebug(o.fDebug),
 130     fOADBFileName(o.fOADBFileName),
 131     fOADBContainer(o.fOADBContainer),
 132     fPhiDist(o.fPhiDist),
 133     fRunNumber(o.fRunNumber),
 134     fUsePhiWeights(o.fUsePhiWeights)
 135 {
 136   //
 137   // Copy constructor
 138   //
 139   // Parameters:
 140   //    o Object to copy from
 141   //
 142   DGUARD(fDebug,0,"Copy CTOR of AliFMDEventPlaneFinder");
 143 }
 144
 145 //____________________________________________________________________
 146 AliFMDEventPlaneFinder::~AliFMDEventPlaneFinder()
 147 {
 148   //
 149   // Destructor
 150   //
 151 }
 152
 153 //____________________________________________________________________
 154 AliFMDEventPlaneFinder&
 155 AliFMDEventPlaneFinder::operator=(const AliFMDEventPlaneFinder& o)
 156 {
 157   //
 158   // Assignement operator
 159   //
 160   // Parameters:
 161   //    o Object to assign from
 162   //
 163   // Return:
 164   //    Reference to this object
 165   //
 166   DGUARD(fDebug,3,"Assignment of AliFMDEventPlaneFinder");
 167   if (&o == this) return *this;
 168   TNamed::operator=(o);
 169
 170   fList                = o.fList;
 171   fEvent               = o.fEvent;
 172   fQt                  = o.fQt;
 173   fQa                  = o.fQa;
 174   fQc                  = o.fQc;
 175   fQ1                  = o.fQ1;
 176   fQ2                  = o.fQ2;
 177   fQeta                = o.fQeta;
 178   fHepFMD              = o.fHepFMD;
 179   fHepFMDA             = o.fHepFMDA;
 180   fHepFMDC             = o.fHepFMDC;
 181   fHepFMDQC1           = o.fHepFMDQC1;
 182   fHepFMDQC2           = o.fHepFMDQC2;
 183   fHdiffFMDAC          = o.fHdiffFMDAC;
 184   fHdiffFMDTPC         = o.fHdiffFMDTPC;
 185   fHdiffFMDVZERO       = o.fHdiffFMDVZERO;
 186   fHcorrFMDAC          = o.fHcorrFMDAC;
 187   fHcorrFMDTPC         = o.fHcorrFMDTPC;
 188   fHcorrFMDVZERO       = o.fHcorrFMDVZERO;
 189   fHPhi                = o.fHPhi;
 190   fDebug               = o.fDebug;
 191   fOADBFileName        = o.fOADBFileName;
 192   fOADBContainer       = o.fOADBContainer;
 193   fPhiDist             = o.fPhiDist;
 194   fRunNumber           = o.fRunNumber;
 195   fUsePhiWeights       = o.fUsePhiWeights;
 196
 197   return *this;
 198 }
 199
 200 //____________________________________________________________________
 201 TH1D*
 202 AliFMDEventPlaneFinder::MakePsiRHist(const char* name,
 203                                      const char* title,
 204                                      Int_t color)
 205 {
 206   // Generate a 1D histogram of Psi_R.
 207   TH1D* ret = new TH1D(name, Form("#Psi_{R} - %s",title), 100, 0, TMath::Pi());
 208   ret->SetDirectory(0);
 209   ret->SetXTitle("#Psi_{R} [radians]");
 210   ret->SetYTitle("Events");
 211   ret->SetLineColor(color);
 212   ret->SetFillColor(color);
 213   ret->SetFillStyle(3001);
 214   ret->Sumw2();
 215   fList->Add(ret);
 216   return ret;
 217 }
 218 //____________________________________________________________________
 219 TH1D*
 220 AliFMDEventPlaneFinder::MakeDiffHist(const char* name,
 221                                      const char* first,
 222                                      const char* second,
 223                                      Int_t color)
 224 {
 225   TH1D* ret = new TH1D(name, Form("#Delta#Psi_{R} - %s minus %s",first,second),
 226                        100, -TMath::Pi()/2, +TMath::Pi()/2);
 227   ret->SetDirectory(0);
 228   ret->SetXTitle(Form("#Psi_{R,%s}-#Psi_{R,%s} [radians]", first, second));
 229   ret->SetYTitle("Events");
 230   ret->SetLineColor(color);
 231   ret->SetFillColor(color);
 232   ret->SetFillStyle(3001);
 233   ret->Sumw2();
 234   fList->Add(ret);
 235   return ret;
 236
 237 }
 238 //____________________________________________________________________
 239 TH2F*
 240 AliFMDEventPlaneFinder::MakeCorrHist(const char* name,
 241                                      const char* first,
 242                                      const char* second)
 243 {
 244   TH2F* ret = new TH2F(name, Form("#Psi_{R} - %s vs %s", first, second),
 245                        100, 0, TMath::Pi(), 100, 0, TMath::Pi());
 246   ret->SetDirectory(0);
 247   ret->Sumw2();
 248   ret->SetXTitle(Form("#Psi_{R,%s} [radians]", first));
 249   ret->SetYTitle(Form("#Psi_{R,%s} [radians]", second));
 250   ret->SetZTitle("Events");
 251   fList->Add(ret);
 252   return ret;
 253 }
 254
 255 //____________________________________________________________________
 256 void
 257 AliFMDEventPlaneFinder::Init(const TAxis& etaAxis)
 258 {
 259   // Intialize this sub-algorithm
 260   //
 261   // Parameters:
 262   //   etaAxis   fmd eta axis binning
 263   //
 264   DGUARD(fDebug,1,"Initalization of AliFMDEventPlaneFinder");
 265
 266   fHepFMD    = MakePsiRHist("epFMD", "FMD", kRed+1);
 267   fHepFMDA   = MakePsiRHist("epFMDA","FMD - A side", kGreen+1);
 268   fHepFMDC   = MakePsiRHist("epFMDC","FMD - C side", kBlue+1);
 269   fHepFMDQC1 = MakePsiRHist("epFMDQC1", "FMD QC{1}", kMagenta+1);
 270   fHepFMDQC2 = MakePsiRHist("epEMDQC2", "FMD QC{2}", kCyan+1);
 271
 272   fHdiffFMDAC     = MakeDiffHist("diffFMDAC",    "FMD-A", "FMD-C", kRed+1);
 273   fHdiffFMDTPC    = MakeDiffHist("diffFMDTPC",   "FMD",   "TPC",   kGreen+1);
 274   fHdiffFMDVZERO  = MakeDiffHist("diffFMDVZERO", "FMD",   "VZERO", kBlue+1);
 275   fHcorrFMDAC     = MakeCorrHist("corrFMDAC",    "FMD-A", "FMD-C");
 276   fHcorrFMDTPC    = MakeCorrHist("corrFMDTPC",   "FMD",   "TPC");
 277   fHcorrFMDVZERO  = MakeCorrHist("corrFMDVZERO", "FMD",   "VZERO");
 278
 279   //
 280   fHPhi = new TH2D("hPhi", "Phi distribution in FMD",
 281                    etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(),
 282                    20, 0., TMath::TwoPi());
 283   fHPhi->Sumw2();
 284   fHPhi->SetXTitle("#eta");
 285   fHPhi->SetYTitle("#phi [radians]");
 286   fHPhi->SetZTitle("Events");
 287   fHPhi->SetDirectory(0);
 288   fList->Add(fHPhi);
 289
 290   //
 291
 292   if (!fUsePhiWeights) return;
 293
 294   if (fOADBFileName.Length()==0)
 295     fOADBFileName = Form("%s/PWGLF/FORWARD/FMDEVENTPLANE/data/fmdEPoadb.root",
 296                         AliAnalysisManager::GetOADBPath());
 297
 298   TFile foadb(fOADBFileName);
 299   if(!foadb.IsOpen()) {
 300     AliError(Form("Cannot open OADB file %s", fOADBFileName.Data()));
 301     return;
 302   }
 303   fOADBContainer = (AliOADBContainer*)foadb.Get("FMDphidist");
 304   if (!fOADBContainer)
 305     AliError(Form("No OADB container found in %s", fOADBFileName.Data()));
 306   foadb.Close();
 307
 308   return;
 309 }
 310
 311 //_____________________________________________________________________
 312 void
 313 AliFMDEventPlaneFinder::DefineOutput(TList* dir)
 314 {
 315   //
 316   // Output diagnostic histograms to directory
 317   //
 318   // Parameters:
 319   //    dir List to write in
 320   //
 321   DGUARD(fDebug,1,"Define output of AliFMDEventPlaneFinder");
 322   fList = new TList();
 323   fList->SetOwner();
 324   fList->SetName(GetName());
 325   dir->Add(fList);
 326
 327   return;
 328 }
 329
 330 //____________________________________________________________________
 331 Bool_t
 332 AliFMDEventPlaneFinder::FindEventplane(AliVEvent* ev,
 333                                        AliAODForwardEP& aodep,
 334                                        TH2D* h,
 335                                        AliForwardUtil::Histos* hists)
 336 {
 337   //
 338   // Do the calculations
 339   //
 340   // Parameters:
 341   //    hists    Histogram cache
 342   //    ep       calculated eventplane
 343   //
 344   // Return:
 345   //    true on successs
 346   DGUARD(fDebug,1,"Find the event plane in AliFMDEventPlaneFinder");
 347
 348   fQt.Set(0., 0.);
 349   fQa.Set(0., 0.);
 350   fQc.Set(0., 0.);
 351   fQ1.Set(0., 0.);
 352   fQ2.Set(0., 0.);
 353
 354   TH1D epEtaHist = aodep.GetHistogram();
 355
 356   fEvent = ev;
 357   if (hists) {
 358     for (UShort_t d=1; d<=3; d++) {
 359       UShort_t nr = (d == 1 ? 1 : 2);
 360       for (UShort_t q=0; q<nr; q++) {
 361         Char_t r = (q == 0 ? 'I' : 'O');
 362         CalcQVectors(hists->Get(d,r), &epEtaHist);
 363       } // for q
 364     } // for d
 365   }
 366   else if (h) {
 367     CalcQVectors(h, &epEtaHist);
 368   }
 369
 370   aodep.SetEventplane(CalcEventplane(fQt));
 371   aodep.SetEventplaneA(CalcEventplane(fQa));
 372   aodep.SetEventplaneC(CalcEventplane(fQc));
 373   // aodep.SetEventplane1(CalcEventplane(fQ1));
 374   // aodep.SetEventplane2(CalcEventplane(fQ2));
 375
 376   FillHists(&aodep);
 377
 378   return kTRUE;
 379 }
 380
 381 //_____________________________________________________________________
 382 void
 383 AliFMDEventPlaneFinder::CalcQVectors(TH2D* h, TH1D* eHist)
 384 {
 385   //
 386   // Calculate the Q vectors
 387   //
 388   DGUARD(fDebug,2,"Calculate Q-vectors in AliFMDEventPlaneFinder");
 389   Double_t phi = 0, eta = 0, weight = 0;
 390   for (Int_t e = 1; e <= h->GetNbinsX(); e++) {
 391     Double_t qx = 0, qy = 0;
 392     eta = h->GetXaxis()->GetBinCenter(e);
 393     for (Int_t p = 1; p <= h->GetNbinsY(); p++) {
 394       phi = h->GetYaxis()->GetBinCenter(p);
 395       weight = h->GetBinContent(e, p);
 396       if (fUsePhiWeights) weight *= GetPhiWeight(e, p);
 397       // fix for FMD1 hole
 398       if (e > 168 && p == 14) {
 399         weight = h->GetBinContent(e, 4);
 400         if (fUsePhiWeights) weight *= GetPhiWeight(e, 4);
 401       }
 402       fHPhi->Fill(eta, phi, weight);
 403       // for underflowbin total Nch/eta
 404       fHPhi->Fill(eta, -1., weight);
 405
 406       // increment Q vectors
 407       qx += weight*TMath::Cos(2.*phi);
 408       qy += weight*TMath::Sin(2.*phi);
 409     }
 410     TVector2 qVec = TVector2(qx, qy);
 411     fQt += qVec;
 412     if (eta < 0) fQa += qVec;
 413     if (eta > 0) fQc += qVec;
 414     // TODO: Add random ep increments
 415     fQeta += qVec;
 416     if (e % 10 == 0) {
 417       eHist->Fill(eta, CalcEventplane(fQeta));
 418       fQeta.Set(0., 0.);
 419     }
 420   }
 421
 422   return;
 423 }
 424
 425 //_____________________________________________________________________
 426 Double_t
 427 AliFMDEventPlaneFinder::CalcEventplane(const TVector2& v) const
 428 {
 429   //
 430   // Calculate the eventplane
 431   //
 432   DGUARD(fDebug,2,"Calculate Event plane in AliFMDEventPlaneFinder");
 433   Double_t ep = -1;
 434
 435   if (v.Mod() == 0.) return ep;
 436
 437   ep = v.Phi();
 438   ep /= 2.;
 439
 440   if (fDebug > 0)
 441     AliInfo(Form("Eventplane found to be: %f", ep));
 442
 443   return ep;
 444 }
 445
 446 //_____________________________________________________________________
 447 Double_t
 448 AliFMDEventPlaneFinder::CalcDifference(Double_t a1, Double_t a2) const
 449 {
 450   Double_t diff = a1 - a2;
 451   if (diff <  TMath::Pi()/2) diff = TMath::Pi() - diff;
 452   if (diff >= TMath::Pi()/2) diff = TMath::Pi() - diff;
 453   return diff;
 454 }
 455
 456 //_____________________________________________________________________
 457 void
 458 AliFMDEventPlaneFinder::FillHists(AliAODForwardEP* fmdEP)
 459 {
 460   //
 461   // Fill diagnostics histograms
 462   //
 463   DGUARD(fDebug,2,"Fill histograms in AliFMDEventPlaneFinder");
 464   Double_t fmdEPt = fmdEP->GetEventplane();
 465   Double_t fmdEPa = fmdEP->GetEventplaneA();
 466   Double_t fmdEPc = fmdEP->GetEventplaneC();
 467   // Double_t fmdEP1 = fmdEP->GetEventplane1();
 468   // Double_t fmdEP2 = fmdEP->GetEventplane2();
 469
 470   // FMD hists
 471   fHepFMD->Fill(fmdEPt);
 472   fHepFMDA->Fill(fmdEPa);
 473   fHepFMDC->Fill(fmdEPc);
 474
 475   fHdiffFMDAC->Fill(CalcDifference(fmdEPa,fmdEPc));
 476   fHcorrFMDAC->Fill(fmdEPa, fmdEPc);
 477   // fHepFMDQC1->Fill(fmdEP1);
 478   // fHepFMDQC2->Fill(fmdEP2);
 479
 480   // TPC comparison
 481   AliEventplane* ep = fEvent->GetEventplane();
 482   Double_t tpcEP   = (ep ? ep->GetEventplane("Q") : -1);
 483   fHdiffFMDTPC->Fill(CalcDifference(fmdEPt, tpcEP));
 484   fHcorrFMDTPC->Fill(fmdEPt, tpcEP);
 485
 486   // VZERO comparison
 487   Double_t vzeroEP = ep ? ep->GetEventplane("V0", fEvent, 2) : -1;
 488   fHdiffFMDVZERO->Fill(CalcDifference(fmdEPt, vzeroEP));
 489   fHcorrFMDVZERO->Fill(fmdEPt, vzeroEP);
 490
 491   return;
 492 }
 493
 494 //_____________________________________________________________________
 495 Double_t
 496 AliFMDEventPlaneFinder::GetPhiWeight(Int_t etaBin, Int_t phiBin) const
 497 {
 498   //
 499   // Get phi weight for flattening
 500   //
 501   Double_t phiWeight = 1;
 502
 503   if (!fPhiDist) return phiWeight;
 504
 505   Double_t nParticles   = fPhiDist->GetBinContent(etaBin, 0);
 506   Double_t nPhiBins     = fPhiDist->GetNbinsY();
 507   Double_t phiDistValue = fPhiDist->GetBinContent(etaBin, phiBin);
 508
 509   if (phiDistValue > 0) phiWeight = nParticles/nPhiBins/phiDistValue;
 510
 511   return phiWeight;
 512 }
 513
 514 //____________________________________________________________________
 515 void
 516 AliFMDEventPlaneFinder::SetRunNumber(Int_t run)
 517 {
 518   //
 519   // Set run number
 520   //
 521   if (fRunNumber == run) return;
 522
 523   fRunNumber = run;
 524   if (fUsePhiWeights) GetPhiDist();
 525
 526   return;
 527 }
 528
 529 //____________________________________________________________________
 530 void
 531 AliFMDEventPlaneFinder::GetPhiDist()
 532 {
 533   //
 534   // Get phi dist from OADB
 535   //
 536   if (!fOADBContainer) return;
 537
 538   fPhiDist = static_cast<TH2D*>(fOADBContainer
 539                                 ->GetObject(fRunNumber, "Default")
 540                                 ->Clone(Form("fPhiDist_%d", fRunNumber)));
 541   fList->Add(fPhiDist);
 542
 543   return;
 544 }
 545
 546 //____________________________________________________________________
 547 void
 548 AliFMDEventPlaneFinder::Print(Option_t* /*option*/) const
 549 {
 550   //
 551   // Print information
 552   //
 553   // Parameters:
 554   //    option Not used
 555   //
 556   char ind[gROOT->GetDirLevel()+3];
 557   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 558   ind[gROOT->GetDirLevel()] = '\0';
 559   std::cout << ind << ClassName() << ": " << GetName() << '\n'
 560             << std::boolalpha
 561             << ind << "Eventplane finder active!" << '\n'
 562             << ind << "Loading OADB object for run number: "
 563             << fRunNumber << '\n'
 564             << ind << "Using phi weights: " << fUsePhiWeights  << '\n'
 565             << std::noboolalpha
 566             << std::endl;
 567 }
 568 //____________________________________________________________________
 569 //
 570 // EOF
 571 //
 572
 573
 574