AliPhysics  857879d (857879d)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliCentralMCCorrectionsTask.cxx
Go to the documentation of this file.
1 //
2 // Calculate the corrections in the central regions
3 //
4 // Inputs:
5 // - AliESDEvent
6 //
7 // Outputs:
8 // - AliAODCentralMult
9 //
10 // Histograms
11 //
12 // Corrections used
13 //
16 #include "AliTriggerAnalysis.h"
17 #include "AliPhysicsSelection.h"
18 #include "AliLog.h"
19 #include "AliHeader.h"
20 #include "AliGenEventHeader.h"
21 #include "AliESDEvent.h"
22 #include "AliAODHandler.h"
23 #include "AliMultiplicity.h"
24 #include "AliInputEventHandler.h"
25 #include "AliStack.h"
26 #include "AliMCEvent.h"
27 #include "AliAODForwardMult.h"
30 #include "AliForwardUtil.h"
31 #include "AliMultiplicity.h"
32 #include <TVector3.h>
33 #include <TH1.h>
34 #include <TH2D.h>
35 #include <TDirectory.h>
36 #include <TList.h>
37 #include <TROOT.h>
38 #include <iostream>
39 
40 //====================================================================
43  fTrackDensity(),
44  fSecCorr(0),
45  fAccCorr(0),
46  fNPhiBins(20),
47  fEffectiveCorr(true),
48  fEtaCut(1.9),
49  fCorrCut(0.6)
50 {
51  //
52  // Constructor
53  //
54  // Parameters:
55  // name Name of task
56  //
57  DGUARD(fDebug, 3,"Default CTOR of AliCentralMCCorrectionsTask");
58 }
59 
60 //____________________________________________________________________
63  fTrackDensity("trackDensity"),
64  fSecCorr(0),
65  fAccCorr(0),
66  fNPhiBins(20),
67  fEffectiveCorr(true),
68  fEtaCut(1.9),
69  fCorrCut(0.6)
70 {
71  //
72  // Constructor
73  //
74  // Parameters:
75  // name Name of task
76  //
77  DGUARD(fDebug, 3,"Named CTOR of AliCentralMCCorrectionsTask: %s",name);
78 }
79 
80 //____________________________________________________________________
83 {
85 }
86 
87 //____________________________________________________________________
88 Bool_t
90  const AliMCEvent& mc,
92  const TVector3& ip)
93 {
95  static_cast<AliCentralMCCorrectionsTask::VtxBin&>(bin);
96 
97  // Now process our input data and store in own ESD object
98  fTrackDensity.Calculate(mc, ip, *vb.fHits, bin.fPrimary);
99 
100  // Get the ESD object
101  const AliMultiplicity* spdmult = esd.GetMultiplicity();
102 
103  // Count number of tracklets per bin
104  for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
105  vb.fClusters->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
106  //...and then the unused clusters in layer 1
107  for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
108  Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
109  vb.fClusters->Fill(eta, spdmult->GetPhiSingle(j));
110  }
111 
112  // Count events
113  bin.fCounts->Fill(0.5);
114 
115  return true;
116 }
117 //____________________________________________________________________
118 void
120 {
123 
126 
127  results->Add(fSecCorr);
128  results->Add(fAccCorr);
129 }
130 
131 //____________________________________________________________________
132 Bool_t
134  bin, UShort_t iVz)
135 {
136 
138  static_cast<AliCentralMCCorrectionsTask::VtxBin*>(bin);
141  return true;
142 }
143 
144 
145 //____________________________________________________________________
146 void
148 {
150  std::cout << " # of phi bins: " << fNPhiBins << "\n"
151  << " Effective corr.: " << fEffectiveCorr << "\n"
152  << " Eta cut-off: " << fEtaCut << "\n"
153  << " Acceptance cut: " << fCorrCut
154  << std::endl;
155  gROOT->IncreaseDirLevel();
156  fTrackDensity.Print(option);
157  gROOT->DecreaseDirLevel();
158 }
159 
160 //====================================================================
163  fHits(0),
164  fClusters(0)
165 {
166 }
167 //____________________________________________________________________
169  Double_t high,
170  const TAxis& axis,
171  UShort_t nPhi)
172  : AliBaseMCCorrectionsTask::VtxBin(low, high, axis, nPhi),
173  fHits(0),
174  fClusters(0)
175 {
176  fHits = static_cast<TH2D*>(fPrimary->Clone("hits"));
177  fHits->SetTitle("Hits");
178  fHits->SetDirectory(0);
179 
180  fClusters = static_cast<TH2D*>(fPrimary->Clone("clusters"));
181  fClusters->SetTitle("Clusters");
182  fClusters->SetDirectory(0);
183 }
184 
185 
186 //____________________________________________________________________
187 TList*
189 {
191 
192  d->Add(fHits);
193  d->Add(fClusters);
194 
195  return d;
196 }
197 //____________________________________________________________________
198 void
200  TList* output,
201  UShort_t iVz,
202  Bool_t effectiveCorr,
203  Double_t etaCut,
204  Double_t accCut,
207 {
208  TList* out = new TList;
209  out->SetName(GetName());
210  out->SetOwner();
211  output->Add(out);
212 
213  TList* l = static_cast<TList*>(input->FindObject(GetName()));
214  if (!l) {
215  AliError(Form("List %s not found in %s", GetName(), input->GetName()));
216  return;
217  }
218 
219  // Get the sums
220  TH2D* hits = static_cast<TH2D*>(l->FindObject("hits"));
221  TH2D* clus = static_cast<TH2D*>(l->FindObject("clusters"));
222  TH2D* prim = static_cast<TH2D*>(l->FindObject("primary"));
223  if (!hits || !prim) {
224  AliError(Form("Missing histograms: %p, %p", hits, prim));
225  return;
226  }
227 
228  // Clone cluster and hit map
229  TH2D* secMapEff = static_cast<TH2D*>(clus->Clone("secMapEff"));
230  TH2D* secMapHit = static_cast<TH2D*>(hits->Clone("secMapHit"));
231  secMapEff->SetTitle("2^{nd} map from clusters");
232  secMapEff->SetDirectory(0);
233  secMapHit->SetTitle("2^{nd} map from MC hits");
234  secMapHit->SetDirectory(0);
235 
236  // Divide cluster and hit map with number of primaries
237  secMapEff->Divide(prim);
238  secMapHit->Divide(prim);
239 
240  // Create acceptance histograms
241  TH1D* accEff = new TH1D("accEff",
242  "Acceptance correction for SPD (from 2^{nd} map)" ,
243  fPrimary->GetXaxis()->GetNbins(),
244  fPrimary->GetXaxis()->GetXmin(),
245  fPrimary->GetXaxis()->GetXmax());
246  TH1D* accHit = static_cast<TH1D*>(accEff->Clone("accHit"));
247  accHit->SetTitle("Acceptance correction for SPD (from clusters)");
248 
249  // Diagnostics histogra,
250  TH2* dia = static_cast<TH2D*>(clus->Clone("diagnostics"));
251  dia->SetTitle("Scaled cluster density");
252 
253  // Total number of channels along phi and # of eta bins
254  Int_t nTotal = secMapHit->GetNbinsY();
255  Int_t nEta = secMapHit->GetNbinsX();
256 
257  for(Int_t xx = 1; xx <= nEta; xx++) {
258  Double_t eta = secMapHit->GetXaxis()->GetBinCenter(xx);
259  Bool_t ins = TMath::Abs(eta) <= etaCut;
260  Double_t mm = 0;
261  if (ins) {
262  // Find the maximum cluster signal in this phi range
263  for (Int_t yy = 1; yy <= nTotal; yy++) {
264  Double_t c = clus->GetBinContent(xx,yy);
265  mm = TMath::Max(mm, c);
266  }
267  }
268  // Count number of phi bins with enough clusters or high enough
269  // correction.
270  Int_t nOKEff = 0;
271  Int_t nOKHit = 0;
272  for(Int_t yy = 1; yy <=nTotal; yy++) {
273  if (!ins) { // Not inside Eta cut
274  secMapEff->SetBinContent(xx,yy,0.);
275  secMapEff->SetBinError(xx,yy,0.);
276  secMapHit->SetBinContent(xx,yy,0.);
277  secMapHit->SetBinError(xx,yy,0.);
278  dia->SetBinContent(xx,yy,0);
279  continue;
280  }
281 
282  // Check if the background correction is large enough, or zero map
283  if(secMapEff->GetBinContent(xx,yy) > 0.9) {
284  // acc->Fill(h->GetXaxis()->GetBinCenter(xx));
285  nOKEff++;
286  }
287  else {
288  secMapEff->SetBinContent(xx,yy,0.);
289  secMapEff->SetBinError(xx,yy,0.);
290  }
291 
292  // Check if the number of cluster is large enough, or zero map
293  Double_t c = clus->GetBinContent(xx,yy);
294  Double_t s = (mm < 1e-6) ? 0 : c / mm;
295  dia->SetBinContent(xx,yy,s);
296  if (s >= accCut) {
297  nOKHit++;
298  }
299  else {
300  secMapHit->SetBinContent(xx,yy,0);
301  secMapHit->SetBinError(xx,yy,0);
302  }
303  }
304 
305  // Calculate acceptance as ratio of bins with enough clusters and
306  // total number of phi bins.
307  Double_t accXX = float(nOKHit) / nTotal;
308  if (accXX < 0.2) accXX = 0;
309  accHit->SetBinContent(xx, accXX);
310 
311  // Calculate acceptance as ratio of bins with large enough
312  // correction and total number of phi bins.
313  accXX = float(nOKEff) / nTotal;
314  if (accXX < 0.2) accXX = 0;
315  accEff->SetBinContent(xx, accXX);
316  }
317 
318  TH2D* secMap = (effectiveCorr ? secMapEff : secMapHit);
319  TH2D* secMapAlt = (effectiveCorr ? secMapHit : secMapEff);
320  TH1D* acc = (effectiveCorr ? accEff : accHit);
321  TH1D* accAlt = (effectiveCorr ? accHit : accEff);
322  out->Add(secMap->Clone("secMap"));
323  out->Add(secMapAlt->Clone());
324  out->Add(acc->Clone("acc"));
325  out->Add(accAlt->Clone());
326  out->Add(dia->Clone());
327 
328  map->SetCorrection(iVz, secMap);
329  acorr->SetCorrection(iVz, acc);
330 }
331 
332 //
333 // EOF
334 //
double Double_t
Definition: External.C:58
void Print(Option_t *option="") const
AliCentralCorrAcceptance * fAccCorr
void SetVertexAxis(const TAxis &axis)
AliBaseMCCorrectionsTask::VtxBin * CreateVtxBin(Double_t low, Double_t high)
TCanvas * c
Definition: TestFitELoss.C:172
Bool_t Calculate(const AliMCEvent &event, const TVector3 &ip, TH2D &output, TH2D *primary)
Per-event per bin.
int Int_t
Definition: External.C:63
virtual void CreateCorrections(TList *results)
Various utilities used in PWGLF/FORWARD.
Definition: External.C:228
Definition: External.C:212
virtual TList * CreateOutputObjects(TList *list)
void SetVertexAxis(const TAxis &axis)
#define DGUARD(L, N, F,...)
void Terminate(const TList *i, TList *o, UShort_t iVz, Bool_t effective, Double_t etaCut, Double_t accCut, AliCentralCorrSecondaryMap *map, AliCentralCorrAcceptance *acorr)
Bool_t SetCorrection(Double_t v, TH2D *h)
virtual Bool_t FinalizeVtxBin(AliBaseMCCorrectionsTask::VtxBin *bin, UShort_t iVz)
AliCentralCorrSecondaryMap * fSecCorr
Definition: External.C:220
Bool_t SetCorrection(Double_t v, TH1D *h)
void Print(Option_t *option="") const
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
void Print(Option_t *option="") const
bool Bool_t
Definition: External.C:53
Bool_t ProcessESD(const AliESDEvent &esd, const AliMCEvent &mc, AliBaseMCCorrectionsTask::VtxBin &bin, const TVector3 &ip)