AliPhysics  32b88a8 (32b88a8)
 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 [%d,%d]",
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  Double_t aVtx = TMath::Abs(vtx);
103  if (aVtx < 37.5 || aVtx > 400) return;
104 
105  TH2* hists[] = { data, dataMC };
106 
107  // In satellite vertices FMD2i is cut away manually at this point
108  // for certain vertices. It could be done in the ESDs, but as of
109  // this writing not for specific vertices.
110  //
111  // cholm comment: It would be difficult to setup the filter in the
112  // reconstruction pass, but it could perhaps be done in the AOD
113  // filtering.
114  //
115  // This is what was done for
116  // the Pb-Pb paper (arXiv:1304.0347).
117  for (Int_t iX = 0; iX<=data->GetNbinsX(); iX++) {
118  // Do all checks up front - as soon as we can - branching is
119  // expensive!
120  Double_t x = data->GetXaxis()->GetBinCenter(iX);
121  Bool_t zero = false;
122  if (((vtx > 60 && vtx < 90) && x < 3) ||
123  ((vtx > 330 && vtx < 350) && x > -2.5) ||
124  ((vtx < 100 || vtx > 305) && TMath::Abs(x) < 4.5) ||
125  (vtx < 50 && TMath::Abs(x) < 4.75))
126  zero = true;
127  if (!zero) continue;
128 
129  for (Int_t iH = 0; iH < 2; iH++) {
130  if (!hists[iH]) continue;
131  // if (iX > hists[iH]->GetNbinsX()+1) continue;
132  // Also zero coverage and phi acceptance for this
133  for (Int_t iY = 0; iY<=hists[iH]->GetNbinsY()+1; iY++) {
134  hists[iH]->SetBinContent(iX, iY, 0);
135  hists[iH]->SetBinError(iX, iY, 0);
136  }
137  }
138  }
139 
140  if (fCorrEmpty) {
141  // Now, since we have some dead areas in FMD2i (sectors 16 and
142  // 17), we need to remove the corresponding bins from the
143  // histogram. However, it is not obvious which bins (in eta) to
144  // remove, so remove everything starting from the most negative to
145  // the middle of the histogram.
146  //
147  // This hack was first introduced by HHD, but was done at the end of
148  // the event processing (CentralityBin::MakeResults). That is,
149  // however, not very practical, as we'd like to normalize to the phi
150  // acceptance rather than the eta coverage and then correct for
151  // empty bins. Since the only way to really update the phi
152  // acceptance stored in the overflow bin is on the event level, we
153  // should really do it here.
154  const Int_t phiBin1 = 17; // Sector 16
155  const Int_t phiBin2 = 18; // Sector 17
156  for (Int_t iH = 0; iH < 2; iH++) {
157  if (!hists[iH]) continue;
158 
159  Int_t midX = hists[iH]->GetNbinsX() / 2;
160  // Int_t nY = hists[iH]->GetNbinsY();
161  for (Int_t i = 1; i <= midX; i++) {
162  hists[iH]->SetBinContent(i, phiBin1, 0);
163  hists[iH]->SetBinContent(i, phiBin2, 0);
164  hists[iH]->SetBinError(i, phiBin1, 0);
165  hists[iH]->SetBinError(i, phiBin2, 0);
166 
167  // Here, we should also modify the overflow bin to reflect the
168  // new phi acceptance. First get the old phi acceptance -
169  // then multiply this on the number of bins. This gives us -
170  // roughly - the number of sectors we had. Then take out two
171  // from that number, and then calculate the new phi
172  // Acceptance. Note, if the sectors where already taken out in
173  // the AOD production, we _will_ end up with a wrong number,
174  // so we should _not_ do that in the AOD production. This is
175  // tricky and may not work at all. For now, we should rely on
176  // the old way of correcting to the eta coverage and
177  // correcting for empty bins.
178  }
179  }
180  }
181 }
182 //____________________________________________________________________
183 Bool_t
185 {
186  TString path(prx);
187  if (gSystem->ExpandPathName(path)) {
188  // Expand with TString argument return 0 on success, 1 on failure
189  return false;
190  }
191  if (!path.Contains("empirical"))
192  path = gSystem->ConcatFileName(path.Data(), "empirical.root");
193  TUrl empUrl(path);
194  if (!empUrl.GetAnchor() || empUrl.GetAnchor()[0] == '\0')
195  empUrl.SetAnchor("default");
196  TFile* empFile = TFile::Open(empUrl.GetUrl());
197  if (!empFile) return false;
198 
199  TString base(GetName()); base.ReplaceAll("dNdeta", "");
200  TString empAnch = empUrl.GetAnchor();
201  TObject* empObj = empFile->Get(Form("%s/%s",base.Data(),empAnch.Data()));
202  if (!(empObj &&
203  (empObj->IsA()->InheritsFrom(TH1::Class()) ||
204  empObj->IsA()->InheritsFrom(TGraphAsymmErrors::Class())))) {
205  Warning("LoadEmpirical", "Didn't get Forward/%s from %s",
206  empAnch.Data(), empUrl.GetUrl());
207  return false;
208  }
209  // Store correction in output list
210  static_cast<TNamed*>(empObj)->SetName("empirical");
211  fResults->Add(empObj);
212 
213  if (!empObj->IsA()->InheritsFrom(TH1::Class()))
214  return true;
215 
216  // Release from directory
217  TH1* h = static_cast<TH1*>(empObj);
218  h->SetDirectory(0);
219 
220  // Check for IP delta
221  TH1* xy = static_cast<TH2*>(fSums->FindObject("vertexAccXY"));
222  Double_t delta = 0;
223  if (xy) {
224  Double_t meanIpx = xy->GetMean(1);
225  Double_t meanIpy = xy->GetMean(2);
226  const Double_t refX = -0.004;
227  const Double_t refY = 0.184;
228  Double_t dx = (meanIpx - refX);
229  Double_t dy = (meanIpy - refY);
230  Info("LoadEmpirical","Shifts (%f-%f)=%f, (%f-%f)=%f",
231  meanIpx, refX, dx,
232  meanIpy, refY, dy);
233  if (TMath::Abs(dx) > 1e-3 || TMath::Abs(dy) > 1e-3)
234  delta = TMath::Sqrt(dx*dx+dy*dy);
235  // Store delta in output list
236  fResults->Add(AliForwardUtil::MakeParameter("deltaIP", delta));
237  }
238  if (delta > 0.2) {
239  // Only correct if delta is larger than 2mm (2% correction)
240  TF1* f = new TF1("deltaCorr", "1+[2]*([0]+(x<[1])*pow([0]*(x-[1]),2))");
241  f->SetParNames("\\delta","\\eta_{0}","a");
242  f->SetParameter(0,delta);
243  f->SetParameter(1,-2.0);
244  f->SetParameter(2,.10); //TMath::Sqrt(2)); // 0.5);
245  Info("LoadEmpirical","Appying correction for IP delta=%f", delta);
246  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
247  Double_t c = h->GetBinContent(i);
248  if (c < 1e-6) continue;
249 
250  Double_t e = h->GetBinError(i);
251  Double_t eta = h->GetXaxis()->GetBinCenter(i);
252  Double_t cor = f->Eval(eta);
253  // Info("", "%5.2f -> %7.4f", eta, cor);
254  h->SetBinContent(i, c*cor);
255  h->SetBinError(i, e*cor);
256  }
257  // Adding correction to result list
258  fResults->Add(f);
259  }
260  return true;
261 }
262 
263 //____________________________________________________________________
264 Bool_t
266 {
267  // See if we can find the empirical correction so that the bins may
268  // apply it
269  TString oadb(gSystem->ConcatFileName(AliAnalysisManager::GetOADBPath(),
270  "PWGLF/FORWARD/CORRECTIONS/data"));
271  const char* dirs[] = {
272  "${PWD}",
273  "${FWD}",
274  oadb.Data(),
275  "${OADB_PATH}/PWGLF/FORWARD/EMPIRICAL",
276  "${ALICE_PHYSICS}/OADB/PWGLF/FORWARD/EMPIRICAL",
277  0
278  };
279  const char** pdir = dirs;
280  Bool_t ok = false;
281  while (*pdir) {
282  const char* fns[] = { "", "empirical_000138190.root", 0 };
283  const char** pfn = fns;
284  while (*pfn) {
285  TString path(gSystem->ConcatFileName(*pdir, *pfn));
286  if ((ok = LoadEmpirical(path))) break;
287  pfn++;
288  }
289  if (ok) break;
290  pdir++;
291  }
293 }
294 
295 //========================================================================
296 TH1*
298 {
299  TList* out = fOutput;
300 
301  TString base(GetName()); base.ReplaceAll("dNdeta", "");
302  TString hName(Form("dndeta%s",base.Data()));
303  TH1* h = static_cast<TH1*>(out->FindObject(hName));
304  if (!h) {
305  Warning("End", "%s not found in %s",
306  hName.Data(), out->GetName());
307  out->ls();
308  return 0;
309  }
310 
311  TObject* o = static_cast<TH1*>(results->FindObject("empirical"));
312  if (!o) {
313  Info("EmpiricalCorrection", "Empirical not found in %s",
314  results->GetName());
315  return 0;
316  }
317 
318  Info("EmpiricalCorrection", "Correcting %s/%s with %s",
319  out->GetName(), h->GetName(), o->GetName());
320 
321  // Make a clone
322  h = static_cast<TH1*>(h->Clone(Form("%sEmp", h->GetName())));
323  h->SetDirectory(0);
324 
325  if (o->IsA()->InheritsFrom(TGraphAsymmErrors::Class())) {
326  TGraphAsymmErrors* empCorr = static_cast<TGraphAsymmErrors*>(o);
327  TAxis* xAxis = h->GetXaxis();
328  for (Int_t i = 1; i <= xAxis->GetNbins(); i++) {
329  Double_t x = xAxis->GetBinCenter(i);
330  Double_t y = h->GetBinContent(i);
331  Double_t c = empCorr->Eval(x);
332  h->SetBinContent(i, y / c);
333  }
334  }
335  else if (o->IsA()->InheritsFrom(TH1::Class())) {
336  TH1* empCorr = static_cast<TH1*>(o);
337  h->Divide(empCorr);
338  }
339  else {
340  Warning("CorrectEmpirical",
341  "Don't know how to apply a %s as an empirical correction",
342  o->IsA()->GetName());
343  delete h;
344  return 0;
345  }
346  // Adding the corrected histogram to output
347  out->Add(h);
348 
349  return h;
350 }
351 //____________________________________________________________________
352 void
354  TList* results,
355  UShort_t scheme,
356  Double_t trigEff,
357  Double_t trigEff0,
358  Bool_t rootProj,
359  Bool_t corrEmpty,
360  Int_t triggerMask,
361  Int_t marker,
362  Int_t color,
363  TList* mclist,
364  TList* truthlist )
365 {
366  DGUARD(fDebug, 1,"In End of %s with corrEmpty=%d, rootProj=%d",
367  GetName(), corrEmpty, rootProj);
368  AliBasedNdetaTask::CentralityBin::End(sums, results, scheme, trigEff,
369  trigEff0, rootProj, corrEmpty,
370  triggerMask, marker, color, mclist,
371  truthlist);
372 
373  // TH1* h = EmpiricalCorrection(results);
374 
375  if (!IsAllBin()) return;
376  TFile* file = TFile::Open("forward.root", "READ");
377  if (!file) return;
378 
379  TList* forward = static_cast<TList*>(file->Get("ForwardSums"));
380  if (!forward) {
381  AliError("List Forward not found in forward.root");
382  return;
383  }
384  TList* rings = static_cast<TList*>(forward->FindObject("ringResults"));
385  if (!rings) {
386  AliError("List ringResults not found in forward.root");
387  return;
388  }
389  THStack* res = static_cast<THStack*>(rings->FindObject("all"));
390  if (!res) {
391  AliError(Form("Stack all not found in %s", rings->GetName()));
392  return;
393  }
394  if (!fTriggers) {
395  AliError("Triggers histogram not set");
396  return;
397  }
398 
399  Double_t ntotal = 0;
400  Double_t epsilonT = trigEff;
401 #if 0
402  // TEMPORARY FIX
403  if (triggerMask == AliAODForwardMult::kNSD) {
404  // This is a local change
405  epsilonT = 0.92;
406  AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
407  }
408 #endif
409  AliInfo("Adding per-ring histograms to output");
410  TString text;
411  Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal, &text);
412  TIter next(res->GetHists());
413  TH1* hist = 0;
414  while ((hist = static_cast<TH1*>(next()))) hist->Scale(scaler);
415  res->SetName("dndetaRings");
416  fOutput->Add(res);
417  fOutput->Add(new TNamed("normCalc", text.Data()));
418 }
419 
420 //________________________________________________________________________
421 //
422 // EOF
423 //
Int_t color[]
double Double_t
Definition: External.C:58
virtual void CheckEventData(Double_t vtx, TH2 *data, TH2 *mcData)
TSystem * gSystem
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