AliPhysics  e6c8d43 (e6c8d43)
AliBaseMultTask.cxx
Go to the documentation of this file.
1 
12 #include <TH1D.h>
13 #include "AliBaseMultTask.h"
14 #include "AliAODForwardMult.h"
15 #include "AliAODCentralMult.h"
16 #include "AliAODEvent.h"
17 #include "AliAODMCHeader.h"
18 #include "AliAODMCParticle.h"
19 
20 
21 ClassImp(AliBaseMultTask)
22 #if 0
23 ; // This is for Emacs - do not delete
24 #endif
25 //_____________________________________________________________________
27  : AliBaseAODTask(),
28  fBins(),
29  fIsSelected(false),
30  fMaxMult(500),
31  fMCIsNSD(false)
32 {
33  //
34  // Default Constructor
35  //
36 }
37 
38 //_____________________________________________________________________
40  : AliBaseAODTask(name,"AliBaseMultTask"),
41  fBins(),
42  fIsSelected(false),
43  fMaxMult(500),
44  fMCIsNSD(false)
45 {
46  //
47  // Constructor
48  //
49 }
50 
51 //_____________________________________________________________________
52 void
54 {
55  //Add Full eta-ranges
56  AddBin(-3.4,5.1);
57 
58  //Add Symmetric eta bins.
59  Double_t limits[] = { 3.4, 3.0, 2.5, 2.4, 2.0, 1.5, 1.4, 1.0, 0.5, 0. };
60  Double_t* limit = limits;
61  while ((*limit) > 0.1) {
62  AddBin(-(*limit), +(*limit));
63  limit++;
64  }
65 }
66 //_____________________________________________________________________
67 Bool_t
69 {
71  // We always return true, so that we can process MC truth;
72  return true;
73 }
74 //_____________________________________________________________________
76 {
77  //
78  // Create Output Objects
79  //
80  //Loop over all individual eta bins, and define their hisograms
81  TIter next(&fBins);
82  Bin * bin = 0;
83  while ((bin = static_cast<Bin*>(next()))) {
85  }
86  return true;
87 }
88 
89 //_____________________________________________________________________
91 {
93 }
94 //_____________________________________________________________________
96 {
97  // if (fMCIsNSD) Printf("Returning true because we assume it");
98  return (fMCIsNSD ? true : m->IsTriggerBits(AliAODForwardMult::kMCNSD));
99 }
100 
101 //_____________________________________________________________________
103 {
104  //
105  // User Exec
106  //
107  // Here, we used to get ForwardMC!
108  AliAODForwardMult* aodForward = GetForward(aod);
109  if (!aodForward) return false;
110 
111  // Here, we used to get CentralClusters!
112  AliAODCentralMult* aodCentral = GetCentral(aod);
113  if (!aodCentral) return false;
114 
115  TH2D& forward = aodForward->GetHistogram();
116  TH2D& central = aodCentral->GetHistogram();
117 
118  //check if event is NSD according to either MC truth or analysis (ESD)
119  Bool_t isMCClass = IsMCClass(aodForward);
120  Bool_t isESDClass = IsESDClass(aodForward);
121  Bool_t isPileUp = aodForward->IsTriggerBits(AliAODForwardMult::kPileUp);
122  //primary dndeta dist
123  TH2D* primHist = GetPrimary(aod);
124  TH1D* dndetaForward = forward.ProjectionX("dndetaForward",1,
125  forward.GetNbinsY(),"");
126  TH1D* dndetaCentral = central.ProjectionX("dndetaCentral",1,
127  central.GetNbinsY(),"");
128  TH1D* dndetaMC = (primHist ?
129  primHist->ProjectionX("dndetaMC",1,
130  primHist->GetNbinsY(),"") :
131  0);
132 
133  // underflow eta bin of forward/central carry information on whether
134  // there is acceptance. 1= acceptance, 0= no acceptance
135  TH1D* normForward = forward.ProjectionX("dndetaForwardNorm",0,0,"");
136  TH1D* normCentral = central.ProjectionX("dndetaCentralNorm",0,0,"");
137 
138  // Interaction point
139  Double_t ipZ = aodForward->GetIpZ();
140 
141  Process(dndetaForward,
142  dndetaCentral,
143  normForward,
144  normCentral,
145  dndetaMC,
146  ipZ,
147  isPileUp,
148  fIsSelected,
149  isMCClass,
150  isESDClass,
151  aod);
152  // Clean up after projections
153  if (dndetaForward) delete dndetaForward;
154  if (dndetaCentral) delete dndetaCentral;
155  if (dndetaMC) delete dndetaMC;
156 
157  return true;
158 }
159 
160 //_____________________________________________________________________
161 void AliBaseMultTask::Process(TH1D* dndetaForward,
162  TH1D* dndetaCentral,
163  TH1D* normForward,
164  TH1D* normCentral,
165  TH1D* dndetaMC,
166  Double_t ipZ,
167  Bool_t pileup,
168  Bool_t selectedTrigger,
169  Bool_t isMCClass,
170  Bool_t isESDClass,
171  const AliAODEvent& aodevent)
172 {
173  // loop over all eta bins,
174  TIter next(&fBins);
175  Bin * bin = 0;
176  while ((bin = static_cast<Bin*>(next()))) {
177  bin->Process(dndetaForward,
178  dndetaCentral,
179  normForward,
180  normCentral,
181  dndetaMC,
182  ipZ,
183  pileup,
184  fIsSelected,
185  isMCClass,
186  isESDClass,
187  aodevent,
188  fIPzAxis.GetXmin(),
189  fIPzAxis.GetXmax());
190  }
191 }
192 
193 //=====================================================================
195  : TNamed("", ""),
196  fEtaLow(etaLow), // low eta limit
197  fEtaHigh(etaHigh), // high eta limit
198  fHist(), // multiplicity histogram
199  fHistMC(), // multiplicity histogram MC truth primaries
200  fAcceptance(), // histogram showing the 'holes' in
201  // acceptance. BinContent of 1 shows a hole,
202  // and BinContent of 10 shows data coverage
203  fVtxZvsNdataBins() // VtxZ vs. number of data acceptance bins
204 {
205  //
206  // Constructor
207  //
208  const char* name = FormBinName(etaLow,etaHigh);
209  SetName(name);
210  SetTitle(Form("%+4.1f < #eta < %+4.1f", fEtaLow, fEtaHigh));
211 
212 }
213 //_____________________________________________________________________
215  : TNamed(),
216  fEtaLow(), // low eta limit
217  fEtaHigh(), // high eta limit
218  fHist(), // multiplicity histogram
219  fHistMC(), // multiplicity histogram MC truth primaries
220  fAcceptance(), // histogram showing the 'holes' in
221  // acceptance. BinContent of 1 shows a hole,
222  // and BinContent of 10 shows data coverage
223  fVtxZvsNdataBins() // VtxZ vs. number of data acceptance bins
224  // (normalised to the eta range)
225 {
226  //
227  // Default Constructor
228  //
229 }
230 //_____________________________________________________________________
232 {
233  //
234  // Define eta bin output histos
235  //
236  TList* out = new TList;
237  out->SetName(GetName());
238  out->SetOwner();
239  cont->Add(out);
240 
241  fHist = new TH1D("mult", GetTitle(), max, -0.5, max-.5);
242  fHistMC = new TH1D("multTruth", GetTitle(), max, -0.5, max-.5);
243  fVtxZvsNdataBins = new TH2D("VtxZvsNdataBins",
244  "VtxZ vs dataAcceptance/etaRange;"
245  "z-vtz;dataAcceptance/etaRange",
246  20, -10,10, 130,0,1.3);
247  fAcceptance = new TH2D("Acceptance","Acceptance;#eta;z-vtx",
248  200,-4, 6 , 20,-10,10);
249 
250  out->Add(fHist);
251  out->Add(fHistMC);
252  out->Add(fVtxZvsNdataBins);
253  out->Add(fAcceptance);
254 }
255 
256 
257 //_____________________________________________________________________
258 Double_t
260  TH1D* dndetaCentral,
261  TH1D* normForward,
262  TH1D* normCentral,
263  TH1D* mc,
264  Double_t ipZ,
265  Double_t& statErr,
266  Double_t& sysErr,
267  Double_t& mcMult,
268  Double_t& mcErr)
269 {
270  // Find bin bounds
271  Int_t first = dndetaForward->GetXaxis()->FindBin(fEtaLow);
272  Int_t last = dndetaForward->GetXaxis()->FindBin(fEtaHigh-.0001);
273 
274  // Fill acceptance histogram
275  Double_t acceptanceBins=0;
276  for(Int_t n= first;n<=last;n++){
277  if(normForward->GetBinContent(n)>0||normCentral->GetBinContent(n)>0){
278  acceptanceBins++;
279  }
280  fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(ipZ), 1);
281  if(normForward->GetBinContent(n)>0||normCentral->GetBinContent(n)>0)
282  fAcceptance->SetBinContent(n, fAcceptance->GetYaxis()->FindBin(ipZ),10);
283  }
284  fVtxZvsNdataBins->Fill(ipZ, (Double_t)acceptanceBins/(last-first+1));
285 
286  // Zero output
287  statErr = 0;
288  sysErr = 0;
289  mcMult = -1;
290  mcErr = -1;
291  // Sum up measurements
292  Double_t c = 0;
293  Double_t e2 = 0;
294  Double_t cPrimary = 0;
295  Double_t e2Primary= 0;
296  Double_t tCentral = 0;
297  Double_t tForward = 0;
298  for (Int_t i = first; i <= last; i++){
299  Double_t cForward = 0;
300  Double_t cCentral = 0;
301  Double_t e2Forward = 0;
302  Double_t e2Central = 0;
303  Bool_t aForward = normForward->GetBinContent(i) > 0;
304  Bool_t aCentral = normCentral->GetBinContent(i) > 0;
305  if (aForward) {
306  cForward = dndetaForward->GetBinContent(i);
307  e2Forward = TMath::Power(dndetaForward->GetBinError(i),2);
308  tForward += cForward;
309  }
310  if (aCentral) {
311  cCentral = dndetaCentral->GetBinContent(i);
312  e2Central = TMath::Power(dndetaCentral->GetBinError(i),2);
313  tCentral += cCentral;
314  }
315  Double_t cc = 0;
316  Double_t ee2 = 0;
317  if (aCentral && aForward) {
318  cc = 0.5 * (cForward + cCentral);
319  ee2 = 0.5 * (e2Forward + e2Central);
320  }
321  else if (aCentral) {
322  cc = cCentral;
323  ee2 = e2Central;
324  }
325  else if (aForward) {
326  cc = cForward;
327  ee2 = e2Forward;
328  }
329  cPrimary += (mc ? mc->GetBinContent(i) : 0);
330  e2Primary += (mc ? TMath::Power(mc->GetBinError(i),2) : 0);
331  c += cc;
332  e2 += ee2;
333  }
334  // Filter out zero-SPD events
335 #if 0
336  if (tCentral < 1e-4) {
337  Printf("Total forward: %f total central: %f -> Argh!",
338  tForward, tCentral);
339  cPrimary = c = -1;
340  }
341 #endif
342 
343  // Systematic errors from here
344  Int_t fmd=0;
345  Int_t spd=0;
346  Int_t overlap=0;
347  // number of eta bins in fmd, spd and overlap respectively
348  for(Int_t i = first;i<=last;i++){
349  if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)<1)
350  fmd++;
351  if(normForward->GetBinContent(i)>0&&normCentral->GetBinContent(i)>0)
352  overlap++;
353  if(normCentral->GetBinContent(i)>0&&normForward->GetBinContent(i)<1)
354  spd++;
355  }
356 
357  Double_t sysErrorSquared = 0;
358  // estimate of systematic uncertainty on the event multiplicity -
359  // hardcoded :(. estimates taken from Hans Hjersing Dalsgaard or
360  // Casper Nygaard phd theses.
361  Double_t fmdSysError= 0.08;
362  Double_t spdSysError= 0.04;
363  Double_t total = 0;
364  total= fmd + spd + overlap;
365  if(total){
366  // Combined systematc event uncertainty, by weighting with the
367  // number of eta-bins of fmd-only, spd-only and the overlap.
368  sysErrorSquared= (fmd*TMath::Power(fmdSysError,2)+
369  spd*TMath::Power(spdSysError,2)+
370  0.5*overlap*TMath::Power(fmdSysError,2)+
371  0.5*overlap*TMath::Power(spdSysError,2))/total;
372  }
373 
374  // Return values
375  statErr = TMath::Sqrt(e2);
376  sysErr = TMath::Sqrt(sysErrorSquared);
377  mcMult = mc ? cPrimary : -1;
378  mcErr = mc ? TMath::Sqrt(e2Primary) : -1;
379 
380  return c;
381 }
382 
383 //_____________________________________________________________________
384 void
386  TH1D* dndetaCentral,
387  TH1D* normForward,
388  TH1D* normCentral,
389  TH1D* mc,
390  Double_t ipZ,
391  Bool_t pileup,
392  Bool_t selectedTrigger,
393  Bool_t isMCClass,
394  Bool_t isESDClass,
395  const AliAODEvent& aodevent,
396  Double_t minIPz,
397  Double_t maxIPz)
398 {
399  //
400  // Process a single eta bin
401  //
402  if (pileup) return;
403  if (!selectedTrigger) return;
404  if (ipZ < minIPz || ipZ > maxIPz) return;
405 
406  Double_t mcMult, mcErr, statErr, sysErr;
407  Double_t mult = CalcMult(dndetaForward,
408  dndetaCentral,
409  normForward,
410  normCentral,
411  mc,
412  ipZ,
413  statErr,
414  sysErr,
415  mcMult,
416  mcErr);
417  fHist->Fill(mult);
418  if (mc) fHistMC->Fill(mcMult);
419 }
420 //______________________________________________________________________
421 const Char_t*
423 {
424  //
425  // Form name of eta bin
426  //
427  TString sLow(TString::Format("%+03d",int(10*etaLow)));
428  TString sHigh(TString::Format("%+03d",int(10*etaHigh)));
429  sLow.ReplaceAll("+", "plus");
430  sLow.ReplaceAll("-", "minus");
431  sHigh.ReplaceAll("+", "plus");
432  sHigh.ReplaceAll("-", "minus");
433  static TString tmp;
434  tmp = TString::Format("%s_%s", sLow.Data(), sHigh.Data());
435  return tmp.Data();
436 }
437 //______________________________________________________________________
438 //
439 // EOF
440 //
virtual Bool_t CheckEvent(const AliAODForwardMult &forward)
const Color_t cc[]
Definition: DrawKs.C:1
double Double_t
Definition: External.C:58
virtual void Process(TH1D *dndetaForward, TH1D *dndetaCentral, TH1D *normForward, TH1D *normCentral, TH1D *mc, Double_t ipZ, Bool_t pileup, Bool_t selectedTrigger, Bool_t isMCClass, Bool_t isESDClass, const AliAODEvent &aodevent)
char Char_t
Definition: External.C:18
TCanvas * c
Definition: TestFitELoss.C:172
virtual Bool_t IsMCClass(AliAODForwardMult *m) const
Base task for multiplicity distribution tasks.
Float_t GetIpZ() const
virtual Bool_t Event(AliAODEvent &aod)
Per-event per bin.
int Int_t
Definition: External.C:63
virtual Bool_t Book()
static Bool_t IsTriggerBits(UInt_t bits, UInt_t trg)
virtual void CreateOutputObjects(TList *cont, Int_t max)
Definition: External.C:228
Definition: External.C:212
virtual Bool_t IsESDClass(AliAODForwardMult *m) const
static const Char_t * FormBinName(Double_t etaLow, Double_t etaHigh)
virtual void DefaultBins()
const TH2D & GetHistogram() const
const char * fwd
void AddBin(Double_t etaLow, Double_t etaHigh)
virtual Double_t CalcMult(TH1D *dndetaForward, TH1D *dndetaCentral, TH1D *normForward, TH1D *normCentral, TH1D *mc, Double_t ipZ, Double_t &statErr, Double_t &sysErr, Double_t &mcMult, Double_t &mcErr)
virtual void Process(TH1D *dndetaForward, TH1D *dndetaCentral, TH1D *normForward, TH1D *normCentral, TH1D *mc, Double_t ipZ, Bool_t pileup, Bool_t selectedTrigger, Bool_t isMCClass, Bool_t isESDClass, const AliAODEvent &aodevent, Double_t minIPz, Double_t maxIPz)
AliAODForwardMult * GetForward(const AliAODEvent &aod, Bool_t mc=false, Bool_t verb=true)
TH2D * GetPrimary(const AliAODEvent &aod)
AliAODCentralMult * GetCentral(const AliAODEvent &aod, Bool_t mc=false, Bool_t verb=true)
bool Bool_t
Definition: External.C:53
virtual Bool_t CheckEvent(const AliAODForwardMult &fwd)
const TH2D & GetHistogram() const