AliPhysics  a6017e1 (a6017e1)
AliFMDCorrector.cxx
Go to the documentation of this file.
1 //
2 // This class calculates the exclusive charged particle density
3 // in each for the 5 FMD rings.
4 //
5 #include "AliFMDCorrector.h"
6 #include <AliESDFMD.h>
7 #include <TAxis.h>
8 #include <TList.h>
9 #include <TMath.h>
11 #include "AliFMDCorrSecondaryMap.h"
12 #include "AliFMDCorrVertexBias.h"
14 #include "AliFMDCorrAcceptance.h"
15 #include "AliLog.h"
16 #include <TH2D.h>
17 #include <TROOT.h>
18 #include <THStack.h>
19 #include <iostream>
20 #include <iomanip>
21 
22 ClassImp(AliFMDCorrector)
23 #if 0
24 ; // For Emacs
25 #endif
26 
27 //____________________________________________________________________
29  : TNamed(),
30  fRingHistos(),
31  fUseSecondaryMap(true),
32  fUseVertexBias(false),
33  fUseAcceptance(false),
34  fUseMergingEfficiency(false),
35  fDebug(0)
36 {
37  // Constructor
38  DGUARD(fDebug, 3, "Default CTOR of AliFMDCorrector");
39 }
40 
41 //____________________________________________________________________
43  : TNamed("fmdCorrector", title),
44  fRingHistos(),
45  fUseSecondaryMap(true),
46  fUseVertexBias(false),
47  fUseAcceptance(false),
48  fUseMergingEfficiency(false),
49  fDebug(0)
50 {
51  // Constructor
52  //
53  // Parameters:
54  // title Title
55  DGUARD(fDebug, 3, "Named CTOR of AliFMDCorrector: %s", title);
56  fRingHistos.SetName(GetName());
57 }
58 
59 //____________________________________________________________________
61  : TNamed(o),
62  fRingHistos(),
67  fDebug(o.fDebug)
68 {
69  // Copy constructor
70  //
71  // Parameters:
72  // o Object to copy from
73  DGUARD(fDebug, 3, "Copy CTOR of AliFMDCorrector");
74  TIter next(&o.fRingHistos);
75  TObject* obj = 0;
76  while ((obj = next())) fRingHistos.Add(obj);
77 }
78 
79 //____________________________________________________________________
81 {
82  // Destructor
83  //
84  //
85  DGUARD(fDebug, 3, "DTOR of AliFMDCorrector");
86  fRingHistos.Delete();
87 }
88 
89 //____________________________________________________________________
92 {
93  // Assignment operator
94  //
95  // Parameters:
96  // o Object to assign from
97  DGUARD(fDebug, 3, "Assignment of AliFMDCorrector");
98  if (&o == this) return *this;
99  TNamed::operator=(o);
100 
101  fDebug = o.fDebug;
102  fRingHistos.Delete();
107  TIter next(&o.fRingHistos);
108  TObject* obj = 0;
109  while ((obj = next())) fRingHistos.Add(obj);
110 
111  return *this;
112 }
113 
114 //____________________________________________________________________
115 void
117 {
118  //
119  // Initialize this object
120  //
121  // Parameters:
122  // etaAxis Eta axis to use
123  //
124  DGUARD(fDebug, 1, "Initialization of AliFMDCorrector");
125  if (!fUseSecondaryMap)
126  AliWarning("Secondary maps not used - BE CAREFUL");
127  if (!fUseVertexBias)
128  AliWarning("Vertex bias not used");
129  if (!fUseAcceptance)
130  AliWarning("Acceptance from dead-channels not used");
131 }
132 
133 //____________________________________________________________________
136 {
137  //
138  // Get the ring histogram container
139  // Parameters:
140  // d Detector
141  // r Ring
142  //
143  // Return:
144  // Ring histogram container
145  //
146  Int_t idx = -1;
147  switch (d) {
148  case 1: idx = 0; break;
149  case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
150  case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
151  }
152  if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
153 
154  return static_cast<RingHistos*>(fRingHistos.At(idx));
155 }
156 
157 //____________________________________________________________________
158 void
159 AliFMDCorrector::DivideMap(TH2* num, const TH2* denom,
160  Bool_t alsoUnderOver) const
161 {
162  //
163  // Implement TH1::Divide but
164  // - Assume compatible histograms
165  // - Unless third argument is true, do not divide over/under flow bins
166  //
167  if (!num || !denom) return;
168 
169  Int_t first = (alsoUnderOver ? 0 : 1);
170  Int_t lastX = num->GetNbinsX() + (alsoUnderOver ? 1 : 0);
171  Int_t lastY = num->GetNbinsY() + (alsoUnderOver ? 1 : 0);
172 
173  for (Int_t ix = first; ix <= lastX; ix++) {
174  for (Int_t iy = first; iy <= lastY; iy++) {
175  Int_t bin = num->GetBin(ix,iy);
176  Double_t c0 = num->GetBinContent(bin);
177  Double_t c1 = denom->GetBinContent(bin);
178  if (!c1) {
179  num->SetBinContent(bin,0);
180  num->SetBinError(bin, 0);
181  continue;
182  }
183  Double_t w = c0 / c1;
184  Double_t e0 = num->GetBinError(bin);
185  Double_t e1 = denom->GetBinError(bin);
186  Double_t c12 = c1*c1;
187  Double_t e2 = (e0*e0*c1*c1 + e1*e1*c0*c0)/(c12*c12);
188 
189  num->SetBinContent(bin, w);
190  num->SetBinError(bin, TMath::Sqrt(e2));
191  }
192  }
193 }
194 //____________________________________________________________________
195 Bool_t
197  UShort_t vtxbin)
198 {
199  //
200  // Do the calculations
201  // Parameters:
202  // hists Cache of histograms
203  // vtxBin Vertex bin
204  //
205  // Return:
206  // true on successs
207  //
208  DGUARD(fDebug, 1, "Correct histograms of AliFMDCorrector");
210 
211  UShort_t uvb = vtxbin;
212  for (UShort_t d=1; d<=3; d++) {
213  UShort_t nr = (d == 1 ? 1 : 2);
214  for (UShort_t q=0; q<nr; q++) {
215  Char_t r = (q == 0 ? 'I' : 'O');
216  TH2D* h = hists.Get(d,r);
217  RingHistos* rh = GetRingHistos(d,r);
218 
219  if (fUseSecondaryMap) {
220  TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d, r, uvb);
221  if (!bg) {
222  AliWarning(Form("No secondary correction for FMDM%d%c "
223  "in vertex bin %d", d, r, uvb));
224  continue;
225  }
226  // Divide by primary/total ratio
227  DivideMap(h, bg, false);
228  }
229  if (fUseVertexBias) {
230  TH2D* ef = fcm.GetVertexBias()->GetCorrection(r, uvb);
231  if (!ef) {
232  AliWarning(Form("No event %s vertex bias correction in vertex bin %d",
233  (r == 'I' || r == 'i' ? "inner" : "outer"), uvb));
234  continue;
235  }
236  // Divide by the event selection efficiency
237  DivideMap(h, ef, false);
238  }
239  if (fUseAcceptance) {
240  TH2D* ac = fcm.GetAcceptance()->GetCorrection(d, r, uvb);
241  if (!ac) {
242  AliWarning(Form("No acceptance correction for FMD%d%c in "
243  "vertex bin %d", d, r, uvb));
244  continue;
245  }
246  // Fill overflow bin with ones
247  for (Int_t i = 1; i <= h->GetNbinsX(); i++)
248  h->SetBinContent(i, h->GetNbinsY()+1, 1);
249 
250  // Divide by the acceptance correction -
251  DivideMap(h, ac, fcm.GetAcceptance()->HasOverflow());
252  }
253 
254  if (fUseMergingEfficiency) {
255  if (!fcm.GetMergingEfficiency()) {
256  AliWarning("No merging efficiencies");
257  continue;
258  }
259  TH1D* sf = fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb);
260  if (!fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb)) {
261  AliWarning(Form("No merging efficiency for FMD%d%c at vertex bin %d",
262  d, r, uvb));
263  continue;
264  }
265 
266 
267  for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
268  Float_t c = sf->GetBinContent(ieta);
269  Float_t ec = sf->GetBinError(ieta);
270 
271  for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) {
272  if (c == 0) {
273  h->SetBinContent(ieta,iphi,0);
274  h->SetBinError(ieta,iphi,0);
275  continue;
276  }
277 
278  Double_t m = h->GetBinContent(ieta, iphi) / c;
279  Double_t em = h->GetBinError(ieta, iphi);
280 
281  Double_t e = TMath::Sqrt(em * em + (m * ec) * (m * ec));
282 
283  h->SetBinContent(ieta,iphi,m);
284  h->SetBinError(ieta,iphi,e);
285  }
286  }
287  }
288 
289  rh->fDensity->Add(h);
290  }
291  }
292 
293  return kTRUE;
294 }
295 
296 //____________________________________________________________________
297 void
299 {
300  //
301  // Scale the histograms to the total number of events
302  // Parameters:
303  // dir Where the output is stored
304  // nEvents Number of events
305  //
306  DGUARD(fDebug, 1, "Scale histograms of AliFMDCorrector");
307  if (nEvents <= 0) return;
308  TList* d = static_cast<TList*>(dir->FindObject(GetName()));
309  if (!d) return;
310 
311  TList* out = new TList;
312  out->SetName(d->GetName());
313  out->SetOwner();
314 
315  fRingHistos.Add(new RingHistos(1, 'I'));
316  fRingHistos.Add(new RingHistos(2, 'I'));
317  fRingHistos.Add(new RingHistos(2, 'O'));
318  fRingHistos.Add(new RingHistos(3, 'I'));
319  fRingHistos.Add(new RingHistos(3, 'O'));
320  TIter next(&fRingHistos);
321  RingHistos* o = 0;
322  THStack* sums = new THStack("sums", "Sums of ring results");
323  while ((o = static_cast<RingHistos*>(next()))) {
324  o->Terminate(d, nEvents);
325  if (!o->fDensity) {
326  Warning("Terminate", "No density from %s", o->GetName());
327  continue;
328  }
329  TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
330  o->fDensity->GetNbinsY(),"e");
331  sum->Scale(1., "width");
332  sum->SetTitle(o->GetName());
333  sum->SetDirectory(0);
334  sum->SetYTitle("#sum N_{ch,primary}");
335  sums->Add(sum);
336  }
337  out->Add(sums);
338  output->Add(out);
339 }
340 //____________________________________________________________________
341 void
343 {
344  //
345  // Output diagnostic histograms to directory
346  //
347  // Parameters:
348  // dir List to write in
349  //
350  DGUARD(fDebug, 1, "Define output of AliFMDCorrector");
351  TList* d = new TList;
352  d->SetOwner();
353  d->SetName(GetName());
354  dir->Add(d);
355 
356  d->Add(AliForwardUtil::MakeParameter("secondary", fUseSecondaryMap));
357  d->Add(AliForwardUtil::MakeParameter("vertexBias", fUseVertexBias));
358  d->Add(AliForwardUtil::MakeParameter("acceptance", fUseAcceptance));
360 
361 
362  fRingHistos.Add(new RingHistos(1, 'I'));
363  fRingHistos.Add(new RingHistos(2, 'I'));
364  fRingHistos.Add(new RingHistos(2, 'O'));
365  fRingHistos.Add(new RingHistos(3, 'I'));
366  fRingHistos.Add(new RingHistos(3, 'O'));
367  TIter next(&fRingHistos);
368  RingHistos* o = 0;
369  while ((o = static_cast<RingHistos*>(next()))) {
370  o->CreateOutputObjects(d);
371  }
372 }
373 
374 #define PF(N,V,...) \
375  AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
376 #define PFB(N,FLAG) \
377  do { \
378  AliForwardUtil::PrintName(N); \
379  std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
380  } while(false)
381 #define PFV(N,VALUE) \
382  do { \
383  AliForwardUtil::PrintName(N); \
384  std::cout << (VALUE) << std::endl; } while(false)
385 
386 //____________________________________________________________________
387 void
388 AliFMDCorrector::Print(Option_t* /* option */) const
389 {
390  //
391  // Print information
392  // Parameters:
393  // option Not used
394  //
396  gROOT->IncreaseDirLevel();
397  PFB("Use secondary maps", fUseSecondaryMap );
398  PFB("Use vertex bias", fUseVertexBias );
399  PFB("Use acceptance", fUseAcceptance );
400  PFB("Use merging efficiency", fUseMergingEfficiency);
401  gROOT->DecreaseDirLevel();
402 }
403 
404 //====================================================================
407  fDensity(0)
408 {
409  // Constructor
410  //
411  //
412 }
413 
414 //____________________________________________________________________
416  : AliForwardUtil::RingHistos(d,r),
417  fDensity(0)
418 {
419  //
420  // Constructor
421  // Parameters:
422  // d detector
423  // r ring
424  //
425  fDensity = new TH2D("primaryDensity",
426  "#sum N_{ch,primary}/(#Delta#eta#Delta#phi)",
427  200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
428  0, 2*TMath::Pi());
429  fDensity->SetDirectory(0);
430  fDensity->SetMarkerColor(Color());
431  fDensity->Sumw2();
432  fDensity->SetXTitle("#eta");
433  fDensity->SetYTitle("#phi [radians]");
434  fDensity->SetZTitle("Primary N_{ch} density");
435 }
436 //____________________________________________________________________
438  : AliForwardUtil::RingHistos(o),
439  fDensity(o.fDensity)
440 {
441  //
442  // Copy constructor
443  // Parameters:
444  // o Object to copy from
445  //
446 }
447 
448 //____________________________________________________________________
451 {
452  //
453  // Assignment operator
454  // Parameters:
455  // o Object to assign from
456  //
457  // Return:
458  // Reference to this
459  //
460  if (&o == this) return *this;
462 
463  if (fDensity) delete fDensity;
464 
465  fDensity = static_cast<TH2D*>(o.fDensity->Clone());
466 
467  return *this;
468 }
469 //____________________________________________________________________
471 {
472  //
473  // Destructor
474  //
475  // if (fDensity) delete fDensity;
476 }
477 
478 //____________________________________________________________________
479 void
481 {
482  //
483  // Make output
484  // Parameters:
485  // dir Where to put it
486  //
487  TList* d = DefineOutputList(dir);
488  d->Add(fDensity);
489 }
490 
491 //____________________________________________________________________
492 void
494 {
495  //
496  // Scale the histograms to the total number of events
497  // Parameters:
498  // dir where the output is stored
499  // nEvents Number of events
500  //
501  TList* l = GetOutputList(dir);
502  if (!l) return;
503 
504  TH2D* density = static_cast<TH2D*>(GetOutputHist(l,"primaryDensity"));
505  if (density) density->Scale(1./nEvents);
506  fDensity = density;
507 }
508 
509 //____________________________________________________________________
510 //
511 // EOF
512 //
513 
514 
515 
TH2D * GetCorrection(UShort_t d, Char_t r, Double_t v) const
RingHistos & operator=(const RingHistos &o)
double Double_t
Definition: External.C:58
TH2D * GetCorrection(Char_t r, Double_t v) const
const char * title
Definition: MakeQAPdf.C:27
RingHistos * GetRingHistos(UShort_t d, Char_t r) const
virtual ~AliFMDCorrector()
TH2D * GetCorrection(UShort_t d, Char_t r, Double_t v) const
char Char_t
Definition: External.C:18
virtual void CreateOutputObjects(TList *dir)
virtual void Print(Option_t *option="") const
const char * GetName() const
TCanvas * c
Definition: TestFitELoss.C:172
#define PFB(N, FLAG)
virtual void Terminate(const TList *dir, TList *output, Int_t nEvents)
void Terminate(TList *dir, Int_t nEvents)
virtual Bool_t Correct(AliForwardUtil::Histos &hists, UShort_t vtxBin)
void CreateOutputObjects(TList *dir)
TH1D * GetCorrection(UShort_t d, Char_t r, Double_t v) const
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
const AliFMDCorrMergingEfficiency * GetMergingEfficiency() const
Definition: External.C:228
Definition: External.C:212
AliFMDCorrector & operator=(const AliFMDCorrector &)
Double_t nEvents
plot quality messages
#define DGUARD(L, N, F,...)
static void PrintTask(const TObject &o)
TH2D * Get(UShort_t d, Char_t r) const
static TObject * MakeParameter(const char *name, UShort_t value)
TH1 * GetOutputHist(const TList *d, const char *name) const
RingHistos & operator=(const RingHistos &o)
void DivideMap(TH2 *num, const TH2 *denom, Bool_t alsoUnderOver=false) const
Definition: External.C:220
const AliFMDCorrAcceptance * GetAcceptance() const
const AliFMDCorrVertexBias * GetVertexBias() const
virtual void SetupForData(const TAxis &etaAxis)
unsigned short UShort_t
Definition: External.C:28
TList * DefineOutputList(TList *d) const
TList * ef
Definition: TestFitELoss.C:145
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
const AliFMDCorrSecondaryMap * GetSecondaryMap() const
Bool_t fUseMergingEfficiency
static AliForwardCorrectionManager & Instance()
TList * GetOutputList(const TList *d) const
TDirectoryFile * dir