AliPhysics  aaf9c62 (aaf9c62)
AliFMDEventPlaneFinder.cxx
Go to the documentation of this file.
1 // This task finds the eventplane
2 // using the FMD
3 //
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 //____________________________________________________________________
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, 3,"Default CTOR of AliFMDEventPlaneFinder");
63 
64 }
65 
66 //____________________________________________________________________
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),
95 {
96  //
97  // Constructor
98  //
99  // Parameters:
100  // name Name of object
101  //
102  DGUARD(fDebug, 3,"Named CTOR of AliFMDEventPlaneFinder: %s", title);
103 }
104 
105 //____________________________________________________________________
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),
128  fHPhi(o.fHPhi),
129  fDebug(o.fDebug),
132  fPhiDist(o.fPhiDist),
135 {
136  //
137  // Copy constructor
138  //
139  // Parameters:
140  // o Object to copy from
141  //
142  DGUARD(fDebug, 3,"Copy CTOR of AliFMDEventPlaneFinder");
143 }
144 
145 //____________________________________________________________________
147 {
148  //
149  // Destructor
150  //
151 }
152 
153 //____________________________________________________________________
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;
189  fHPhi = o.fHPhi;
190  fDebug = o.fDebug;
193  fPhiDist = o.fPhiDist;
196 
197  return *this;
198 }
199 
200 //____________________________________________________________________
201 TH1D*
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*
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*
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
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
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
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 
373  // aodep.SetEventplane1(CalcEventplane(fQ1));
374  // aodep.SetEventplane2(CalcEventplane(fQ2));
375 
376  FillHists(&aodep);
377 
378  return kTRUE;
379 }
380 
381 //_____________________________________________________________________
382 void
384 {
385  //
386  // Calculate the Q vectors
387  //
388  DGUARD(fDebug,5,"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(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
428 {
429  //
430  // Calculate the eventplane
431  //
432  DGUARD(fDebug,6,"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
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
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
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
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
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 #define PFV(N,VALUE) \
547  do { \
548  AliForwardUtil::PrintName(N); \
549  std::cout << (VALUE) << std::endl; } while(false)
550 #define PFB(N,FLAG) \
551  do { \
552  AliForwardUtil::PrintName(N); \
553  std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
554  } while(false)
555 
556 //____________________________________________________________________
557 void
559 {
560  //
561  // Print information
562  //
563  // Parameters:
564  // option Not used
565  //
567  gROOT->IncreaseDirLevel();
568  PFB("Is active", true);
569  PFV("Run number for OADB query", fRunNumber);
570  PFB("Use phi weights", fUsePhiWeights);
571  gROOT->DecreaseDirLevel();
572 }
573 //____________________________________________________________________
574 //
575 // EOF
576 //
577 
578 
579 
Int_t color[]
print message on plot with ok/not ok
void FillHists(AliAODForwardEP *fmdEP)
double Double_t
Definition: External.C:58
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:27
virtual void SetupForData(const TAxis &etaAxis)
#define PFB(N, FLAG)
Double_t CalcDifference(Double_t a1, Double_t a2) const
AliFMDEventPlaneFinder & operator=(const AliFMDEventPlaneFinder &o)
char Char_t
Definition: External.C:18
const TH1D & GetHistogram() const
#define PFV(N, VALUE)
Double_t CalcEventplane(const TVector2 &v) const
TH2F * MakeCorrHist(const char *name, const char *first, const char *second)
void SetEventplaneA(Double_t epA)
virtual void CreateOutputObjects(TList *dir)
AliOADBContainer * fOADBContainer
void SetEventplaneC(Double_t epC)
void CalcQVectors(TH2D *h, TH1D *eHist)
Double_t GetEventplaneA()
Per-event per bin.
int Int_t
Definition: External.C:63
Double_t GetEventplaneC()
void Print(Option_t *option="") const
Definition: External.C:228
Definition: External.C:212
TH1D * MakePsiRHist(const char *name, const char *title, Int_t color)
#define DGUARD(L, N, F,...)
TH1D * MakeDiffHist(const char *name, const char *first, const char *second, Int_t color)
static void PrintTask(const TObject &o)
void SetEventplane(Double_t ep)
TH2D * Get(UShort_t d, Char_t r) const
Double_t GetEventplane()
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
Bool_t FindEventplane(AliVEvent *esd, AliAODForwardEP &aodEp, TH2D *h, AliForwardUtil::Histos *hists)
Double_t GetPhiWeight(Int_t etaBin, Int_t phiBin) const
TDirectoryFile * dir