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