AliPhysics  3337bb0 (3337bb0)
 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  {
199  if (gSystem->ExpandPathName(path)) {
200  // Expand with TString argument return 0 on success, 1 on failure
201  return false;
202  }
203  if (!path.Contains("empirical"))
204  path = gSystem->ConcatFileName(path.Data(), "empirical.root");
205  empUrl.SetUrl(path);
206  if (!empUrl.GetAnchor() || empUrl.GetAnchor()[0] == '\0')
207  empUrl.SetAnchor("default");
208  empFile = TFile::Open(empUrl.GetUrl());
209  if (!empFile) {
210  DMSG(fDebug,1,"%s not found", empUrl.GetUrl());
211  return false;
212  }
213  }
214  DMSG(fDebug,0,"Got empirical file %s", empUrl.GetUrl());
215 
216  TString base(GetName()); base.ReplaceAll("dNdeta", "");
217  TString empAnch = empUrl.GetAnchor();
218  TObject* empObj = empFile->Get(Form("%s/%s",base.Data(),empAnch.Data()));
219  if (!(empObj &&
220  (empObj->IsA()->InheritsFrom(TH1::Class()) ||
221  empObj->IsA()->InheritsFrom(TF1::Class()) ||
222  empObj->IsA()->InheritsFrom(TGraphAsymmErrors::Class())))) {
223  Warning("LoadEmpirical", "Didn't get Forward/%s from %s",
224  empAnch.Data(), empUrl.GetUrl());
225  return false;
226  }
227  DMSG(fDebug,0,"Got empirical correction %s [%s]", empObj->GetName(),
228  empObj->ClassName());
229  // Store correction in output list
230  static_cast<TNamed*>(empObj)->SetName("empirical");
231  fResults->Add(empObj);
232 
233  if (!empObj->IsA()->InheritsFrom(TH1::Class()))
234  return true;
235 
236  // Release from directory
237  TH1* h = static_cast<TH1*>(empObj);
238  h->SetDirectory(0);
239 
240  // Check for IP delta
241  TH1* xy = static_cast<TH2*>(fSums->FindObject("vertexAccXY"));
242  Double_t delta = 0;
243  if (xy) {
244  Double_t meanIpx = xy->GetMean(1);
245  Double_t meanIpy = xy->GetMean(2);
246  const Double_t refX = -0.004;
247  const Double_t refY = 0.184;
248  Double_t dx = (meanIpx - refX);
249  Double_t dy = (meanIpy - refY);
250  Info("LoadEmpirical","Shifts (%f-%f)=%f, (%f-%f)=%f",
251  meanIpx, refX, dx,
252  meanIpy, refY, dy);
253  if (TMath::Abs(dx) > 1e-3 || TMath::Abs(dy) > 1e-3)
254  delta = TMath::Sqrt(dx*dx+dy*dy);
255  // Store delta in output list
256  fResults->Add(AliForwardUtil::MakeParameter("deltaIP", delta));
257  }
258  if (delta > 0.2) {
259  // Only correct if delta is larger than 2mm (2% correction)
260  TF1* f = new TF1("deltaCorr", "1+[2]*([0]+(x<[1])*pow([0]*(x-[1]),2))");
261  f->SetParNames("\\delta","\\eta_{0}","a");
262  f->SetParameter(0,delta);
263  f->SetParameter(1,-2.0);
264  f->SetParameter(2,.10); //TMath::Sqrt(2)); // 0.5);
265  Info("LoadEmpirical","Appying correction for IP delta=%f", delta);
266  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
267  Double_t c = h->GetBinContent(i);
268  if (c < 1e-6) continue;
269 
270  Double_t e = h->GetBinError(i);
271  Double_t eta = h->GetXaxis()->GetBinCenter(i);
272  Double_t cor = f->Eval(eta);
273  // Info("", "%5.2f -> %7.4f", eta, cor);
274  h->SetBinContent(i, c*cor);
275  h->SetBinError(i, e*cor);
276  }
277  // Adding correction to result list
278  fResults->Add(f);
279  }
280  return true;
281 }
282 
283 //____________________________________________________________________
284 Bool_t
286 {
287  // See if we can find the empirical correction so that the bins may
288  // apply it
289  TString oadb(gSystem->ConcatFileName(AliAnalysisManager::GetOADBPath(),
290  "PWGLF/FORWARD/CORRECTIONS/data"));
291  const char* dirs[] = {
292  "${PWD}",
293  "${FWD}",
294  oadb.Data(),
295  "${OADB_PATH}/PWGLF/FORWARD/EMPIRICAL",
296  "${ALICE_PHYSICS}/OADB/PWGLF/FORWARD/EMPIRICAL",
297  0
298  };
299  const char** pdir = dirs;
300  Bool_t ok = false;
301  while (*pdir) {
302  const char* fns[] = { "", "empirical_000138190.root", 0 };
303  const char** pfn = fns;
304  const char* dir = *pdir;
305  pdir++;
306  while (*pfn) {
308  TString path(gSystem->ConcatFileName(dir, *pfn));
309  pfn++;
310  const char* ancs[] = { "param", "default", 0 };
311  const char** pan = ancs;
312  while (*pan) {
313  const char* anch = *pan;
314  pan++;
315  TString u(path); u.Append("#"); u.Append(anch);
316  if ((ok = LoadEmpirical(u))) break;
317  }
318  if (ok) break;
319  }
320  if (ok) break;
321  }
323 }
324 
325 //========================================================================
326 TH1*
328 {
329  TList* out = fOutput;
330 
331  TString base(GetName()); base.ReplaceAll("dNdeta", "");
332  TString hName(Form("dndeta%s",base.Data()));
333  TH1* h = static_cast<TH1*>(out->FindObject(hName));
334  if (!h) {
335  Warning("End", "%s not found in %s",
336  hName.Data(), out->GetName());
337  out->ls();
338  return 0;
339  }
340 
341  TObject* o = static_cast<TH1*>(results->FindObject("empirical"));
342  if (!o) {
343  Info("EmpiricalCorrection", "Empirical not found in %s",
344  results->GetName());
345  return 0;
346  }
347 
348  Info("EmpiricalCorrection", "Correcting %s/%s with %s [%s]",
349  out->GetName(), h->GetName(), o->GetName(), o->ClassName());
350 
351  // Make a clone
352  h = static_cast<TH1*>(h->Clone(Form("%sEmp", h->GetName())));
353  h->SetDirectory(0);
354 
355  if (o->IsA()->InheritsFrom(TGraphAsymmErrors::Class())) {
356  TGraphAsymmErrors* empCorr = static_cast<TGraphAsymmErrors*>(o);
357  TAxis* xAxis = h->GetXaxis();
358  for (Int_t i = 1; i <= xAxis->GetNbins(); i++) {
359  Double_t x = xAxis->GetBinCenter(i);
360  Double_t y = h->GetBinContent(i);
361  Double_t c = empCorr->Eval(x);
362  h->SetBinContent(i, y / c);
363  }
364  }
365  else if (o->IsA()->InheritsFrom(TF1::Class())) {
366  TF1* empCorr = static_cast<TF1*>(o);
367  h->Divide(empCorr);
368  }
369  else if (o->IsA()->InheritsFrom(TH1::Class())) {
370  TH1* empCorr = static_cast<TH1*>(o);
371  h->Divide(empCorr);
372  }
373  else {
374  Warning("CorrectEmpirical",
375  "Don't know how to apply a %s as an empirical correction",
376  o->IsA()->GetName());
377  delete h;
378  return 0;
379  }
380  // Adding the corrected histogram to output
381  out->Add(h);
382 
383  return h;
384 }
385 //____________________________________________________________________
386 void
388  TList* results,
389  UShort_t scheme,
390  Double_t trigEff,
391  Double_t trigEff0,
392  Bool_t rootProj,
393  Bool_t corrEmpty,
394  Int_t triggerMask,
395  Int_t marker,
396  Int_t color,
397  TList* mclist,
398  TList* truthlist )
399 {
400  DGUARD(fDebug, 1,"In End of %s with corrEmpty=%d, rootProj=%d",
401  GetName(), corrEmpty, rootProj);
402  AliBasedNdetaTask::CentralityBin::End(sums, results, scheme, trigEff,
403  trigEff0, rootProj, corrEmpty,
404  triggerMask, marker, color, mclist,
405  truthlist);
406 
407  TH1* h = EmpiricalCorrection(results);
408  Info("End", "Applied empirical correction: %p (%s)",
409  h, h ? h->GetName() : "");
410 
411  if (!IsAllBin()) return;
412 
413  THStack* res = 0;
414  {
415  if (gSystem->AccessPathName("forward.root")) return;
416 
417  TFile* file = TFile::Open("forward.root", "READ");
418  if (!file) return;
419 
420  TList* forward = static_cast<TList*>(file->Get("ForwardSums"));
421  if (!forward) {
422  AliError("List Forward not found in forward.root");
423  return;
424  }
425  TList* rings = static_cast<TList*>(forward->FindObject("ringResults"));
426  if (!rings) {
427  AliError("List ringResults not found in forward.root");
428  return;
429  }
430  res = static_cast<THStack*>(rings->FindObject("all"));
431  if (!res) {
432  AliError(Form("Stack all not found in %s", rings->GetName()));
433  return;
434  }
435  }
436  if (!fTriggers) {
437  AliError("Triggers histogram not set");
438  return;
439  }
440 
441  Double_t ntotal = 0;
442  Double_t epsilonT = trigEff;
443 #if 0
444  // TEMPORARY FIX
445  if (triggerMask == AliAODForwardMult::kNSD) {
446  // This is a local change
447  epsilonT = 0.92;
448  AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
449  }
450 #endif
451  AliInfo("Adding per-ring histograms to output");
452  TString text;
453  Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal, &text);
454  TIter next(res->GetHists());
455  TH1* hist = 0;
456  while ((hist = static_cast<TH1*>(next()))) hist->Scale(scaler);
457  res->SetName("dndetaRings");
458  fOutput->Add(res);
459  fOutput->Add(new TNamed("normCalc", text.Data()));
460 }
461 
462 //________________________________________________________________________
463 //
464 // EOF
465 //
Int_t color[]
option to what and if export to output file
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
TList with histograms for a given trigger.
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
TDirectoryFile * dir