AliPhysics  ebc57ae (ebc57ae)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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),
94  fUsePhiWeights(1)
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),
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, 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[]
option to what and if export to output file
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)
ClassImp(AliFMDEventPlaneFinder) AliFMDEventPlaneFinder
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