1 // This task finds the eventplane
4 #include "AliFMDEventPlaneFinder.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"
24 ClassImp(AliFMDEventPlaneFinder)
29 //____________________________________________________________________
30 AliFMDEventPlaneFinder::AliFMDEventPlaneFinder()
62 DGUARD(fDebug, 3,"Default CTOR of AliFMDEventPlaneFinder");
66 //____________________________________________________________________
67 AliFMDEventPlaneFinder::AliFMDEventPlaneFinder(const char* title)
68 : TNamed("fmdEventPlaneFinder", title),
100 // name Name of object
102 DGUARD(fDebug, 3,"Named CTOR of AliFMDEventPlaneFinder: %s", title);
105 //____________________________________________________________________
106 AliFMDEventPlaneFinder::AliFMDEventPlaneFinder(const
107 AliFMDEventPlaneFinder& o)
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),
130 fOADBFileName(o.fOADBFileName),
131 fOADBContainer(o.fOADBContainer),
132 fPhiDist(o.fPhiDist),
133 fRunNumber(o.fRunNumber),
134 fUsePhiWeights(o.fUsePhiWeights)
140 // o Object to copy from
142 DGUARD(fDebug, 3,"Copy CTOR of AliFMDEventPlaneFinder");
145 //____________________________________________________________________
146 AliFMDEventPlaneFinder::~AliFMDEventPlaneFinder()
153 //____________________________________________________________________
154 AliFMDEventPlaneFinder&
155 AliFMDEventPlaneFinder::operator=(const AliFMDEventPlaneFinder& o)
158 // Assignement operator
161 // o Object to assign from
164 // Reference to this object
166 DGUARD(fDebug,3,"Assignment of AliFMDEventPlaneFinder");
167 if (&o == this) return *this;
168 TNamed::operator=(o);
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;
191 fOADBFileName = o.fOADBFileName;
192 fOADBContainer = o.fOADBContainer;
193 fPhiDist = o.fPhiDist;
194 fRunNumber = o.fRunNumber;
195 fUsePhiWeights = o.fUsePhiWeights;
200 //____________________________________________________________________
202 AliFMDEventPlaneFinder::MakePsiRHist(const char* name,
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);
218 //____________________________________________________________________
220 AliFMDEventPlaneFinder::MakeDiffHist(const char* name,
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);
238 //____________________________________________________________________
240 AliFMDEventPlaneFinder::MakeCorrHist(const char* name,
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);
248 ret->SetXTitle(Form("#Psi_{R,%s} [radians]", first));
249 ret->SetYTitle(Form("#Psi_{R,%s} [radians]", second));
250 ret->SetZTitle("Events");
255 //____________________________________________________________________
257 AliFMDEventPlaneFinder::SetupForData(const TAxis& etaAxis)
259 // Intialize this sub-algorithm
262 // etaAxis fmd eta axis binning
264 DGUARD(fDebug,1,"Initalization of AliFMDEventPlaneFinder");
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);
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");
280 fHPhi = new TH2D("hPhi", "Phi distribution in FMD",
281 etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(),
282 20, 0., TMath::TwoPi());
284 fHPhi->SetXTitle("#eta");
285 fHPhi->SetYTitle("#phi [radians]");
286 fHPhi->SetZTitle("Events");
287 fHPhi->SetDirectory(0);
292 if (!fUsePhiWeights) return;
294 if (fOADBFileName.Length()==0)
295 fOADBFileName = Form("%s/PWGLF/FORWARD/FMDEVENTPLANE/data/fmdEPoadb.root",
296 AliAnalysisManager::GetOADBPath());
298 TFile foadb(fOADBFileName);
299 if(!foadb.IsOpen()) {
300 AliError(Form("Cannot open OADB file %s", fOADBFileName.Data()));
303 fOADBContainer = (AliOADBContainer*)foadb.Get("FMDphidist");
305 AliError(Form("No OADB container found in %s", fOADBFileName.Data()));
311 //_____________________________________________________________________
313 AliFMDEventPlaneFinder::CreateOutputObjects(TList* dir)
316 // Output diagnostic histograms to directory
319 // dir List to write in
321 DGUARD(fDebug,1,"Define output of AliFMDEventPlaneFinder");
324 fList->SetName(GetName());
330 //____________________________________________________________________
332 AliFMDEventPlaneFinder::FindEventplane(AliVEvent* ev,
333 AliAODForwardEP& aodep,
335 AliForwardUtil::Histos* hists)
338 // Do the calculations
341 // hists Histogram cache
342 // ep calculated eventplane
346 DGUARD(fDebug,1,"Find the event plane in AliFMDEventPlaneFinder");
354 TH1D epEtaHist = aodep.GetHistogram();
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);
367 CalcQVectors(h, &epEtaHist);
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));
381 //_____________________________________________________________________
383 AliFMDEventPlaneFinder::CalcQVectors(TH2D* h, TH1D* eHist)
386 // Calculate the Q vectors
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);
398 if (e > 168 && p == 14) {
399 weight = h->GetBinContent(e, 4);
400 if (fUsePhiWeights) weight *= GetPhiWeight(e, 4);
402 fHPhi->Fill(eta, phi, weight);
403 // for underflowbin total Nch/eta
404 fHPhi->Fill(eta, -1., weight);
406 // increment Q vectors
407 qx += weight*TMath::Cos(2.*phi);
408 qy += weight*TMath::Sin(2.*phi);
410 TVector2 qVec(qx, qy);
412 if (eta < 0) fQa += qVec;
413 if (eta > 0) fQc += qVec;
414 // TODO: Add random ep increments
417 eHist->Fill(eta, CalcEventplane(fQeta));
425 //_____________________________________________________________________
427 AliFMDEventPlaneFinder::CalcEventplane(const TVector2& v) const
430 // Calculate the eventplane
432 DGUARD(fDebug,2,"Calculate Event plane in AliFMDEventPlaneFinder");
435 if (v.Mod() == 0.) return ep;
441 AliInfo(Form("Eventplane found to be: %f", ep));
446 //_____________________________________________________________________
448 AliFMDEventPlaneFinder::CalcDifference(Double_t a1, Double_t a2) const
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;
456 //_____________________________________________________________________
458 AliFMDEventPlaneFinder::FillHists(AliAODForwardEP* fmdEP)
461 // Fill diagnostics histograms
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();
471 fHepFMD->Fill(fmdEPt);
472 fHepFMDA->Fill(fmdEPa);
473 fHepFMDC->Fill(fmdEPc);
475 fHdiffFMDAC->Fill(CalcDifference(fmdEPa,fmdEPc));
476 fHcorrFMDAC->Fill(fmdEPa, fmdEPc);
477 // fHepFMDQC1->Fill(fmdEP1);
478 // fHepFMDQC2->Fill(fmdEP2);
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);
487 Double_t vzeroEP = ep ? ep->GetEventplane("V0", fEvent, 2) : -1;
488 fHdiffFMDVZERO->Fill(CalcDifference(fmdEPt, vzeroEP));
489 fHcorrFMDVZERO->Fill(fmdEPt, vzeroEP);
494 //_____________________________________________________________________
496 AliFMDEventPlaneFinder::GetPhiWeight(Int_t etaBin, Int_t phiBin) const
499 // Get phi weight for flattening
501 Double_t phiWeight = 1;
503 if (!fPhiDist) return phiWeight;
505 Double_t nParticles = fPhiDist->GetBinContent(etaBin, 0);
506 Double_t nPhiBins = fPhiDist->GetNbinsY();
507 Double_t phiDistValue = fPhiDist->GetBinContent(etaBin, phiBin);
509 if (phiDistValue > 0) phiWeight = nParticles/nPhiBins/phiDistValue;
514 //____________________________________________________________________
516 AliFMDEventPlaneFinder::SetRunNumber(Int_t run)
521 if (fRunNumber == run) return;
524 if (fUsePhiWeights) GetPhiDist();
529 //____________________________________________________________________
531 AliFMDEventPlaneFinder::GetPhiDist()
534 // Get phi dist from OADB
536 if (!fOADBContainer) return;
538 fPhiDist = static_cast<TH2D*>(fOADBContainer
539 ->GetObject(fRunNumber, "Default")
540 ->Clone(Form("fPhiDist_%d", fRunNumber)));
541 fList->Add(fPhiDist);
546 //____________________________________________________________________
548 AliFMDEventPlaneFinder::Print(Option_t* /*option*/) const
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'
561 << ind << "Eventplane finder active!" << '\n'
562 << ind << "Loading OADB object for run number: "
563 << fRunNumber << '\n'
564 << ind << "Using phi weights: " << fUsePhiWeights << '\n'
568 //____________________________________________________________________