AliPhysics  97344c9 (97344c9)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliForwarddNdetaTask.cxx
Go to the documentation of this file.
1 //====================================================================
2 #include "AliForwarddNdetaTask.h"
3 #include <TMath.h>
4 #include <TH2D.h>
5 #include <TH1D.h>
6 #include <THStack.h>
7 #include <TList.h>
8 #include <TFile.h>
9 #include <TF1.h>
10 #include <TSystem.h>
11 #include <TGraphAsymmErrors.h>
12 #include <AliAnalysisManager.h>
13 #include <AliAODEvent.h>
14 #include <AliAODHandler.h>
15 #include <AliAODInputHandler.h>
16 #include "AliForwardUtil.h"
17 #include "AliAODForwardMult.h"
18 
19 //____________________________________________________________________
22 {
23  //
24  // Constructor
25  //
26  DGUARD(fDebug, 3, "Default CTOR of AliForwarddNdetaTask");
27 }
28 
29 //____________________________________________________________________
31  : AliBasedNdetaTask("Forward")
32 {
33  //
34  // Constructor
35  //
36  // Paramters
37  // name Name of task
38  // SetTitle("FMD");
39  DGUARD(fDebug, 3, "Named CTOR of AliForwarddNdetaTask");
40 }
41 
42 //____________________________________________________________________
45 {
46  //
47  // Copy constructor
48  //
49  DGUARD(fDebug, 3, "Copy CTOR of AliForwarddNdetaTask");
50 }
51 
52 //____________________________________________________________________
55  const
56 {
57  //
58  // Make a new centrality bin
59  //
60  // Parameters:
61  // name Histogram names
62  // l Lower cut
63  // h Upper cut
64  //
65  // Return:
66  // Newly allocated object (of our type)
67  //
68  DGUARD(fDebug, 3,
69  "Make a centrality bin for AliForwarddNdetaTask: %s [%5.1f%%,%5.1f%%]",
70  name, l, h);
71  return new AliForwarddNdetaTask::CentralityBin(name, l, h);
72 }
73 
74 
75 //____________________________________________________________________
76 TH2D*
78 {
79  //
80  // Retrieve the histogram
81  //
82  // Parameters:
83  // aod AOD event
84  // mc Whether to get the MC histogram or not
85  //
86  // Return:
87  // Retrieved histogram or null
88  //
89  // We should have a forward object at least
90  AliAODForwardMult* forward = GetForward(aod, mc, !mc);
91  if (!forward) return 0;
92  return &(forward->GetHistogram());
93 }
94 //____________________________________________________________________
95 void
97  TH2* data,
98  TH2* dataMC)
99 {
100  // Check if this is satellite
101  // if (!fSatelliteVertices) return;
102 
103  // If we don't care about satellites, get out
104  if (vtx < fIPzAxis.GetXmin() || vtx > fIPzAxis.GetXmax()) return;
105 
106  // Check if this is a satelllite
107  Double_t aVtx = TMath::Abs(vtx);
108  if (aVtx < 37.5 || aVtx > 400) return;
109 
110  DMSG(fDebug,0,"Got satelitte vertex %f", vtx);
111 
112  TH2* hists[] = { data, dataMC };
113 
114  // In satellite vertices FMD2i is cut away manually at this point
115  // for certain vertices. It could be done in the ESDs, but as of
116  // this writing not for specific vertices.
117  //
118  // cholm comment: It would be difficult to setup the filter in the
119  // reconstruction pass, but it could perhaps be done in the AOD
120  // filtering.
121  //
122  // This is what was done for
123  // the Pb-Pb paper (arXiv:1304.0347).
124  for (Int_t iX = 0; iX<=data->GetNbinsX(); iX++) {
125  // Do all checks up front - as soon as we can - branching is
126  // expensive!
127  Double_t x = data->GetXaxis()->GetBinCenter(iX);
128  Bool_t zero = false;
129  if (((vtx > 60 && vtx < 90) && x < 3) ||
130  ((vtx > 330 && vtx < 350) && x > -2.5) ||
131  ((vtx < 100 || vtx > 305) && TMath::Abs(x) < 4.5) ||
132  (vtx < 50 && TMath::Abs(x) < 4.75))
133  zero = true;
134  if (!zero) continue;
135  for (Int_t iH = 0; iH < 2; iH++) {
136  if (!hists[iH]) continue;
137  // if (iX > hists[iH]->GetNbinsX()+1) continue;
138  // Also zero coverage and phi acceptance for this
139  for (Int_t iY = 0; iY<=hists[iH]->GetNbinsY()+1; iY++) {
140  hists[iH]->SetBinContent(iX, iY, 0);
141  hists[iH]->SetBinError(iX, iY, 0);
142  }
143  }
144  }
145 
146  if (fCorrEmpty) {
147  DMSG(fDebug,1,"Correcting with corrEmpty=true");
148  // Now, since we have some dead areas in FMD2i (sectors 16 and
149  // 17), we need to remove the corresponding bins from the
150  // histogram. However, it is not obvious which bins (in eta) to
151  // remove, so remove everything starting from the most negative to
152  // the middle of the histogram.
153  //
154  // This hack was first introduced by HHD, but was done at the end of
155  // the event processing (CentralityBin::MakeResults). That is,
156  // however, not very practical, as we'd like to normalize to the phi
157  // acceptance rather than the eta coverage and then correct for
158  // empty bins. Since the only way to really update the phi
159  // acceptance stored in the overflow bin is on the event level, we
160  // should really do it here.
161  const Int_t phiBin1 = 17; // Sector 16
162  const Int_t phiBin2 = 18; // Sector 17
163  for (Int_t iH = 0; iH < 2; iH++) {
164  if (!hists[iH]) continue;
165 
166  Int_t midX = hists[iH]->GetNbinsX() / 2;
167  // Int_t nY = hists[iH]->GetNbinsY();
168  for (Int_t i = 1; i <= midX; i++) {
169  hists[iH]->SetBinContent(i, phiBin1, 0);
170  hists[iH]->SetBinContent(i, phiBin2, 0);
171  hists[iH]->SetBinError(i, phiBin1, 0);
172  hists[iH]->SetBinError(i, phiBin2, 0);
173 
174  // Here, we should also modify the overflow bin to reflect the
175  // new phi acceptance. First get the old phi acceptance -
176  // then multiply this on the number of bins. This gives us -
177  // roughly - the number of sectors we had. Then take out two
178  // from that number, and then calculate the new phi
179  // Acceptance. Note, if the sectors where already taken out in
180  // the AOD production, we _will_ end up with a wrong number,
181  // so we should _not_ do that in the AOD production. This is
182  // tricky and may not work at all. For now, we should rely on
183  // the old way of correcting to the eta coverage and
184  // correcting for empty bins.
185  }
186  }
187  }
188 }
189 //____________________________________________________________________
190 Bool_t
192 {
193  TString path(prx);
194  TUrl empUrl;
195  TFile* empFile = 0;
196 
197  if (gSystem->ExpandPathName(path)) {
198  // Expand with TString argument return 0 on success, 1 on failure
199  return false;
200  }
201  if (!path.Contains("empirical"))
202  path = gSystem->ConcatFileName(path.Data(), "empirical.root");
203  empUrl.SetUrl(path);
204  if (!empUrl.GetAnchor() || empUrl.GetAnchor()[0] == '\0')
205  empUrl.SetAnchor("default");
206  empFile = TFile::Open(empUrl.GetUrl());
207  if (!empFile) {
208  DMSG(fDebug,1,"%s not found", empUrl.GetUrl());
209  return false;
210  }
211  DMSG(fDebug,0,"Got empirical file %s", empUrl.GetUrl());
212 
213  TString base(GetName()); base.ReplaceAll("dNdeta", "");
214  TString empAnch = empUrl.GetAnchor();
215  TObject* empObj = empFile->Get(Form("%s/%s",base.Data(),empAnch.Data()));
216  if (!(empObj &&
217  (empObj->IsA()->InheritsFrom(TH1::Class()) ||
218  empObj->IsA()->InheritsFrom(TGraphAsymmErrors::Class())))) {
219  Warning("LoadEmpirical", "Didn't get Forward/%s from %s",
220  empAnch.Data(), empUrl.GetUrl());
221  return false;
222  }
223  DMSG(fDebug,0,"Got empirical correction %s [%s]", empObj->GetName(),
224  empObj->ClassName());
225  // Store correction in output list
226  static_cast<TNamed*>(empObj)->SetName("empirical");
227  fResults->Add(empObj);
228 
229  if (!empObj->IsA()->InheritsFrom(TH1::Class()))
230  return true;
231 
232  // Release from directory
233  TH1* h = static_cast<TH1*>(empObj);
234  h->SetDirectory(0);
235 
236  // Check for IP delta
237  TH1* xy = static_cast<TH2*>(fSums->FindObject("vertexAccXY"));
238  Double_t delta = 0;
239  if (xy) {
240  Double_t meanIpx = xy->GetMean(1);
241  Double_t meanIpy = xy->GetMean(2);
242  const Double_t refX = -0.004;
243  const Double_t refY = 0.184;
244  Double_t dx = (meanIpx - refX);
245  Double_t dy = (meanIpy - refY);
246  Info("LoadEmpirical","Shifts (%f-%f)=%f, (%f-%f)=%f",
247  meanIpx, refX, dx,
248  meanIpy, refY, dy);
249  if (TMath::Abs(dx) > 1e-3 || TMath::Abs(dy) > 1e-3)
250  delta = TMath::Sqrt(dx*dx+dy*dy);
251  // Store delta in output list
252  fResults->Add(AliForwardUtil::MakeParameter("deltaIP", delta));
253  }
254  if (delta > 0.2) {
255  // Only correct if delta is larger than 2mm (2% correction)
256  TF1* f = new TF1("deltaCorr", "1+[2]*([0]+(x<[1])*pow([0]*(x-[1]),2))");
257  f->SetParNames("\\delta","\\eta_{0}","a");
258  f->SetParameter(0,delta);
259  f->SetParameter(1,-2.0);
260  f->SetParameter(2,.10); //TMath::Sqrt(2)); // 0.5);
261  Info("LoadEmpirical","Appying correction for IP delta=%f", delta);
262  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
263  Double_t c = h->GetBinContent(i);
264  if (c < 1e-6) continue;
265 
266  Double_t e = h->GetBinError(i);
267  Double_t eta = h->GetXaxis()->GetBinCenter(i);
268  Double_t cor = f->Eval(eta);
269  // Info("", "%5.2f -> %7.4f", eta, cor);
270  h->SetBinContent(i, c*cor);
271  h->SetBinError(i, e*cor);
272  }
273  // Adding correction to result list
274  fResults->Add(f);
275  }
276  return true;
277 }
278 
279 //____________________________________________________________________
280 Bool_t
282 {
283  // See if we can find the empirical correction so that the bins may
284  // apply it
285  TString oadb(gSystem->ConcatFileName(AliAnalysisManager::GetOADBPath(),
286  "PWGLF/FORWARD/CORRECTIONS/data"));
287  const char* dirs[] = {
288  "${PWD}",
289  "${FWD}",
290  oadb.Data(),
291  "${OADB_PATH}/PWGLF/FORWARD/EMPIRICAL",
292  "${ALICE_PHYSICS}/OADB/PWGLF/FORWARD/EMPIRICAL",
293  0
294  };
295  const char** pdir = dirs;
296  Bool_t ok = false;
297  while (*pdir) {
298  const char* fns[] = { "", "empirical_000138190.root", 0 };
299  const char** pfn = fns;
300  while (*pfn) {
302  TString path(gSystem->ConcatFileName(*pdir, *pfn));
303  if ((ok = LoadEmpirical(path))) break;
304  pfn++;
305  }
306  if (ok) break;
307  pdir++;
308  }
310 }
311 
312 //========================================================================
313 TH1*
315 {
316  TList* out = fOutput;
317 
318  TString base(GetName()); base.ReplaceAll("dNdeta", "");
319  TString hName(Form("dndeta%s",base.Data()));
320  TH1* h = static_cast<TH1*>(out->FindObject(hName));
321  if (!h) {
322  Warning("End", "%s not found in %s",
323  hName.Data(), out->GetName());
324  out->ls();
325  return 0;
326  }
327 
328  TObject* o = static_cast<TH1*>(results->FindObject("empirical"));
329  if (!o) {
330  Info("EmpiricalCorrection", "Empirical not found in %s",
331  results->GetName());
332  return 0;
333  }
334 
335  Info("EmpiricalCorrection", "Correcting %s/%s with %s",
336  out->GetName(), h->GetName(), o->GetName());
337 
338  // Make a clone
339  h = static_cast<TH1*>(h->Clone(Form("%sEmp", h->GetName())));
340  h->SetDirectory(0);
341 
342  if (o->IsA()->InheritsFrom(TGraphAsymmErrors::Class())) {
343  TGraphAsymmErrors* empCorr = static_cast<TGraphAsymmErrors*>(o);
344  TAxis* xAxis = h->GetXaxis();
345  for (Int_t i = 1; i <= xAxis->GetNbins(); i++) {
346  Double_t x = xAxis->GetBinCenter(i);
347  Double_t y = h->GetBinContent(i);
348  Double_t c = empCorr->Eval(x);
349  h->SetBinContent(i, y / c);
350  }
351  }
352  else if (o->IsA()->InheritsFrom(TH1::Class())) {
353  TH1* empCorr = static_cast<TH1*>(o);
354  h->Divide(empCorr);
355  }
356  else {
357  Warning("CorrectEmpirical",
358  "Don't know how to apply a %s as an empirical correction",
359  o->IsA()->GetName());
360  delete h;
361  return 0;
362  }
363  // Adding the corrected histogram to output
364  out->Add(h);
365 
366  return h;
367 }
368 //____________________________________________________________________
369 void
371  TList* results,
372  UShort_t scheme,
373  Double_t trigEff,
374  Double_t trigEff0,
375  Bool_t rootProj,
376  Bool_t corrEmpty,
377  Int_t triggerMask,
378  Int_t marker,
379  Int_t color,
380  TList* mclist,
381  TList* truthlist )
382 {
383  DGUARD(fDebug, 1,"In End of %s with corrEmpty=%d, rootProj=%d",
384  GetName(), corrEmpty, rootProj);
385  AliBasedNdetaTask::CentralityBin::End(sums, results, scheme, trigEff,
386  trigEff0, rootProj, corrEmpty,
387  triggerMask, marker, color, mclist,
388  truthlist);
389 
390  TH1* h = EmpiricalCorrection(results);
391  Info("End", "Applied empirical correction: %p (%s)",
392  h, h ? h->GetName() : "");
393 
394  if (!IsAllBin()) return;
395 
396  THStack* res = 0;
397  {
399  TFile* file = TFile::Open("forward.root", "READ");
400  if (!file) return;
401 
402  TList* forward = static_cast<TList*>(file->Get("ForwardSums"));
403  if (!forward) {
404  AliError("List Forward not found in forward.root");
405  return;
406  }
407  TList* rings = static_cast<TList*>(forward->FindObject("ringResults"));
408  if (!rings) {
409  AliError("List ringResults not found in forward.root");
410  return;
411  }
412  res = static_cast<THStack*>(rings->FindObject("all"));
413  if (!res) {
414  AliError(Form("Stack all not found in %s", rings->GetName()));
415  return;
416  }
417  }
418  if (!fTriggers) {
419  AliError("Triggers histogram not set");
420  return;
421  }
422 
423  Double_t ntotal = 0;
424  Double_t epsilonT = trigEff;
425 #if 0
426  // TEMPORARY FIX
427  if (triggerMask == AliAODForwardMult::kNSD) {
428  // This is a local change
429  epsilonT = 0.92;
430  AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
431  }
432 #endif
433  AliInfo("Adding per-ring histograms to output");
434  TString text;
435  Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal, &text);
436  TIter next(res->GetHists());
437  TH1* hist = 0;
438  while ((hist = static_cast<TH1*>(next()))) hist->Scale(scaler);
439  res->SetName("dndetaRings");
440  fOutput->Add(res);
441  fOutput->Add(new TNamed("normCalc", text.Data()));
442 }
443 
444 //________________________________________________________________________
445 //
446 // EOF
447 //
Int_t color[]
double Double_t
Definition: External.C:58
virtual void CheckEventData(Double_t vtx, TH2 *data, TH2 *mcData)
TSystem * gSystem
#define DMSG(L, N, F,...)
TCanvas * c
Definition: TestFitELoss.C:172
Bool_t LoadEmpirical(const char *path)
Definition: External.C:92
virtual void End(TList *sums, TList *results, UShort_t scheme, Double_t trigEff, Double_t trigEff0, Bool_t rootProj, Bool_t corrEmpty, Int_t triggerMask, Int_t marker, Int_t color, TList *mclist, TList *truthlist)
virtual Bool_t Finalize()
Per-event per bin.
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Various utilities used in PWGLF/FORWARD.
Definition: External.C:228
#define DGUARD(L, N, F,...)
static TObject * MakeParameter(const char *name, UShort_t value)
TH2D * GetHistogram(const AliAODEvent &aod, Bool_t mc)
Definition: External.C:220
AliAODForwardMult * GetForward(const AliAODEvent &aod, Bool_t mc=false, Bool_t verb=true)
AliBasedNdetaTask::CentralityBin * MakeCentralityBin(const char *name, Float_t l, Float_t h) const
TFile * file
virtual void End(TList *sums, TList *results, UShort_t scheme, Double_t trigEff, Double_t trigEff0, Bool_t rootProj, Bool_t corrEmpty, Int_t triggerMask, Int_t marker, Int_t color, TList *mclist, TList *truthlist)
unsigned short UShort_t
Definition: External.C:28
bool Bool_t
Definition: External.C:53
Definition: External.C:196
const TH2D & GetHistogram() const