1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 /* $Id: AliUEHist.cxx 20164 2007-08-14 15:31:50Z morsch $ */
20 // encapsulate histogram and corrections for one underlying event histogram
23 // Author: Jan Fiete Grosse-Oetringhaus, Sara Vallero
25 #include "AliUEHist.h"
26 #include "AliCFContainer.h"
27 #include "THnSparse.h"
30 #include "TCollection.h"
42 const Int_t AliUEHist::fgkCFSteps = 11;
44 AliUEHist::AliUEHist(const char* reqHist) :
48 fTrackHistEfficiency(0),
58 fContaminationEnhancement(0),
62 fGetMultCacheOn(kFALSE),
64 fHistogramType(reqHist)
68 for (UInt_t i=0; i<fkRegions; i++)
71 if (strlen(reqHist) == 0)
74 Printf("Creating AliUEHist with %s", reqHist);
76 AliLog::SetClassDebugLevel("AliCFContainer", -1);
77 AliLog::SetClassDebugLevel("AliCFGridSparse", -3);
79 const char* title = "";
82 Int_t nTrackVars = 4; // eta vs pT vs pT,lead (vs delta phi) vs multiplicity
84 Double_t* trackBins[6];
85 const char* trackAxisTitle[6];
88 const Int_t kNEtaBins = 20;
89 Double_t etaBins[20+1];
90 for (Int_t i=0; i<=kNEtaBins; i++)
91 etaBins[i] = -1.0 + 0.1 * i;
92 const char* etaTitle = "#eta";
93 iTrackBin[0] = kNEtaBins;
94 trackBins[0] = etaBins;
95 trackAxisTitle[0] = etaTitle;
98 const Int_t kNDeltaEtaBins = 48;
99 Double_t deltaEtaBins[kNDeltaEtaBins+1] = { -2.4, -2.3, -2.2, -2.1,
100 -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1,
101 -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1,
103 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
104 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0,
105 2.1, 2.2, 2.3, 2.4 };
107 const Int_t kNDeltaEtaBinsTTR = 40+4;
108 Double_t deltaEtaBinsTTR[kNDeltaEtaBinsTTR+1] = { -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1,
109 -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1,
110 -0.05, -0.025, 0, 0.025, 0.05,
111 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
112 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0 };
116 Double_t pTBins[] = {0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 8.0 };
117 trackBins[1] = pTBins;
118 trackAxisTitle[1] = "p_{T} (GeV/c)";
121 const Int_t kNpTBinsFine = 22;
122 Double_t pTBinsFine[kNpTBinsFine+1] = {0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0};
125 //const Int_t kNpTBinsWideRange = 20;
126 //Double_t pTBinsWideRange[kNpTBinsWideRange+1] = {0.15, 0.2, 0.3, 0.4, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 10.0, 12.0, 15.0};
129 const Int_t kNLeadingpTBins = 100;
130 Double_t leadingpTBins[kNLeadingpTBins+1];
131 for (Int_t i=0; i<=kNLeadingpTBins; i++)
132 leadingpTBins[i] = 0.5 * i;
135 //const Int_t kNLeadingpTBins2 = 8;
136 // Double_t leadingpTBins2[] = { 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0 };
137 const Int_t kNLeadingpTBins2 = 6;
138 Double_t leadingpTBins2[] = { 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0 };
140 // phi,lead; this binning starts at -pi/2 and is modulo 3
141 const Int_t kNLeadingPhiSpacing = 72;
142 const Int_t kNLeadingPhiBins = kNLeadingPhiSpacing;
143 Double_t leadingPhiBins[kNLeadingPhiBins+1];
144 for (Int_t i=0; i<=kNLeadingPhiSpacing; i++)
145 leadingPhiBins[i] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * i;
147 const Int_t kNLeadingPhiBinsTTR = kNLeadingPhiSpacing+4;
148 Double_t leadingPhiBinsTTR[kNLeadingPhiBinsTTR+1];
149 for (Int_t i=0; i<=17; i++)
150 leadingPhiBinsTTR[i] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * i;
151 leadingPhiBinsTTR[18] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * 17 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() / 2;
152 leadingPhiBinsTTR[19] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * 17 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() / 2 * 1.5;
153 leadingPhiBinsTTR[20] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * 17 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() / 2 * 2; // = 0
154 leadingPhiBinsTTR[21] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * 17 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() / 2 * 2.5;
155 leadingPhiBinsTTR[22] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * 17 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() / 2 * 3;
156 for (Int_t i=19; i<=kNLeadingPhiSpacing; i++)
157 leadingPhiBinsTTR[i+4] = -TMath::Pi() / 2 + 1.0 / kNLeadingPhiSpacing * TMath::TwoPi() * i;
160 const Int_t kNMultiplicityBins = 15;
161 Double_t multiplicityBins[kNMultiplicityBins+1];
162 for (Int_t i=0; i<=kNMultiplicityBins; i++)
163 multiplicityBins[i] = -0.5 + i;
164 multiplicityBins[kNMultiplicityBins] = 200;
166 trackBins[3] = multiplicityBins;
167 iTrackBin[3] = kNMultiplicityBins;
168 trackAxisTitle[3] = "multiplicity";
171 const Int_t kNCentralityBins = 5 + 1 + 9;
172 Double_t centralityBins[] = { 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1 };
174 const Int_t kNCentralityBinsCourse = 5;
175 Double_t centralityBinsCourse[] = { 0, 20, 40, 60, 80, 100.1 };
178 const Int_t kNSpeciesBins = 4; // pi, K, p, rest
179 Double_t speciesBins[] = { -0.5, 0.5, 1.5, 2.5, 3.5 };
182 const char* vertexTitle = "z-vtx (cm)";
183 const Int_t kNVertexBins = 7;
184 Double_t vertexBins[] = { -7, -5, -3, -1, 1, 3, 5, 7 };
185 const Int_t kNVertexBins2 = 10;
186 Double_t vertexBins2[] = { -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10 };
188 Int_t useVtxAxis = 0;
189 Bool_t useTTRBinning = kFALSE;
190 Bool_t useCourseCentralityBinning = kFALSE;
192 // selection depending on requested histogram
193 Int_t axis = -1; // 0 = pT,lead, 1 = phi,lead
194 if (strcmp(reqHist, "NumberDensitypT") == 0)
197 title = "d^{2}N_{ch}/d#varphid#eta";
199 else if (strcmp(reqHist, "NumberDensityPhi") == 0)
202 title = "d^{2}N_{ch}/d#varphid#eta";
204 else if (TString(reqHist).BeginsWith("NumberDensityPhiCentrality"))
206 if (TString(reqHist).Contains("Vtx"))
208 if (TString(reqHist).Contains("Vtx10"))
210 if (TString(reqHist).Contains("TTR"))
211 useTTRBinning = kTRUE;
212 if (TString(reqHist).Contains("Course"))
213 useCourseCentralityBinning = kTRUE;
215 reqHist = "NumberDensityPhiCentrality";
216 fHistogramType = reqHist;
218 title = "d^{2}N_{ch}/d#varphid#eta";
220 else if (strcmp(reqHist, "SumpT") == 0)
223 title = "d^{2}#Sigma p_{T}/d#varphid#eta";
226 AliFatal(Form("Invalid histogram requested: %s", reqHist));
228 UInt_t initRegions = fkRegions;
232 trackBins[2] = leadingpTBins;
233 iTrackBin[2] = kNLeadingpTBins;
234 trackAxisTitle[2] = "leading p_{T} (GeV/c)";
242 iTrackBin[2] = kNLeadingpTBins2;
243 trackBins[2] = leadingpTBins2;
244 trackAxisTitle[2] = "leading p_{T} (GeV/c)";
246 iTrackBin[4] = kNLeadingPhiBins;
247 trackBins[4] = leadingPhiBins;
248 trackAxisTitle[4] = "#Delta #varphi w.r.t. leading track";
255 iTrackBin[0] = (useTTRBinning) ? kNDeltaEtaBinsTTR : kNDeltaEtaBins;
256 trackBins[0] = (useTTRBinning) ? deltaEtaBinsTTR : deltaEtaBins;
257 trackAxisTitle[0] = "#Delta#eta";
259 iTrackBin[2] = kNLeadingpTBins2;
260 trackBins[2] = leadingpTBins2;
261 trackAxisTitle[2] = "leading p_{T} (GeV/c)";
263 trackBins[3] = (useCourseCentralityBinning) ? centralityBinsCourse : centralityBins;
264 iTrackBin[3] = (useCourseCentralityBinning) ? kNCentralityBinsCourse : kNCentralityBins;
265 trackAxisTitle[3] = "centrality";
267 iTrackBin[4] = (useTTRBinning) ? kNLeadingPhiBinsTTR : kNLeadingPhiBins;
268 trackBins[4] = (useTTRBinning) ? leadingPhiBinsTTR : leadingPhiBins;
269 trackAxisTitle[4] = "#Delta#varphi (rad.)";
274 iTrackBin[5] = (useVtxAxis == 1) ? kNVertexBins : kNVertexBins2;
275 trackBins[5] = (useVtxAxis == 1) ? vertexBins : vertexBins2;
276 trackAxisTitle[5] = vertexTitle;
280 for (UInt_t i=0; i<initRegions; i++)
282 if (strcmp(reqHist, "NumberDensityPhiCentrality") == 0)
283 fTrackHist[i] = new AliTHn(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
285 fTrackHist[i] = new AliCFContainer(Form("fTrackHist_%d", i), title, fgkCFSteps, nTrackVars, iTrackBin);
287 for (Int_t j=0; j<nTrackVars; j++)
289 fTrackHist[i]->SetBinLimits(j, trackBins[j]);
290 fTrackHist[i]->SetVarTitle(j, trackAxisTitle[j]);
293 SetStepNames(fTrackHist[i]);
297 Int_t nEventVars = 2;
300 // track 3rd and 4th axis --> event 1st and 2nd axis
301 iEventBin[0] = iTrackBin[2];
302 iEventBin[1] = iTrackBin[3];
304 // plus track 5th axis (in certain cases)
305 if (axis == 2 && useVtxAxis)
308 iEventBin[2] = iTrackBin[5];
311 fEventHist = new AliCFContainer("fEventHist", title, fgkCFSteps, nEventVars, iEventBin);
313 fEventHist->SetBinLimits(0, trackBins[2]);
314 fEventHist->SetVarTitle(0, trackAxisTitle[2]);
316 fEventHist->SetBinLimits(1, trackBins[3]);
317 fEventHist->SetVarTitle(1, trackAxisTitle[3]);
319 if (axis == 2 && useVtxAxis)
321 fEventHist->SetBinLimits(2, trackBins[5]);
322 fEventHist->SetVarTitle(2, trackAxisTitle[5]);
325 SetStepNames(fEventHist);
327 iTrackBin[1] = kNpTBinsFine;
328 iTrackBin[2] = kNSpeciesBins;
329 iTrackBin[4] = kNVertexBins2;
331 fTrackHistEfficiency = new AliCFContainer("fTrackHistEfficiency", "Tracking efficiency", 4, 5, iTrackBin);
332 fTrackHistEfficiency->SetBinLimits(0, trackBins[0]);
333 fTrackHistEfficiency->SetVarTitle(0, trackAxisTitle[0]);
334 fTrackHistEfficiency->SetBinLimits(1, pTBinsFine);
335 fTrackHistEfficiency->SetVarTitle(1, trackAxisTitle[1]);
336 fTrackHistEfficiency->SetBinLimits(2, speciesBins);
337 fTrackHistEfficiency->SetVarTitle(2, "particle species");
338 fTrackHistEfficiency->SetBinLimits(3, trackBins[3]);
339 fTrackHistEfficiency->SetVarTitle(3, trackAxisTitle[3]);
340 fTrackHistEfficiency->SetBinLimits(4, vertexBins2);
341 fTrackHistEfficiency->SetVarTitle(4, vertexTitle);
343 fFakePt = new TH3F("fFakePt","fFakePt;p_{T,rec};p_{T};centrality", 200, 0, 20, 200, 0, 20, 20, 0, 100);
346 //_____________________________________________________________________________
347 AliUEHist::AliUEHist(const AliUEHist &c) :
351 fTrackHistEfficiency(0),
361 fContaminationEnhancement(0),
365 fGetMultCacheOn(kFALSE),
370 // AliUEHist copy constructor
373 ((AliUEHist &) c).Copy(*this);
376 //____________________________________________________________________
377 void AliUEHist::SetStepNames(AliCFContainer* container)
379 // sets the names of the correction steps
381 for (Int_t i=0; i<fgkCFSteps; i++)
382 container->SetStepTitle(i, GetStepTitle((CFStep) i));
385 //____________________________________________________________________
386 AliUEHist::~AliUEHist()
390 for (UInt_t i=0; i<fkRegions; i++)
394 delete fTrackHist[i];
405 if (fTrackHistEfficiency)
407 delete fTrackHistEfficiency;
408 fTrackHistEfficiency = 0;
424 //____________________________________________________________________
425 AliUEHist &AliUEHist::operator=(const AliUEHist &c)
427 // assigment operator
430 ((AliUEHist &) c).Copy(*this);
435 //____________________________________________________________________
436 void AliUEHist::Copy(TObject& c) const
440 AliUEHist& target = (AliUEHist &) c;
442 for (UInt_t i=0; i<fkRegions; i++)
444 target.fTrackHist[i] = dynamic_cast<AliCFContainer*> (fTrackHist[i]->Clone());
447 target.fEventHist = dynamic_cast<AliCFContainer*> (fEventHist->Clone());
449 if (fTrackHistEfficiency)
450 target.fTrackHistEfficiency = dynamic_cast<AliCFContainer*> (fTrackHistEfficiency->Clone());
453 target.fFakePt = dynamic_cast<TH3F*> (fFakePt->Clone());
455 target.fEtaMin = fEtaMin;
456 target.fEtaMax = fEtaMax;
457 target.fPtMin = fPtMin;
458 target.fPtMax = fPtMax;
459 target.fCentralityMin = fCentralityMin;
460 target.fCentralityMax = fCentralityMax;
461 target.fZVtxMin = fZVtxMin;
462 target.fZVtxMax = fZVtxMax;
464 if (fContaminationEnhancement)
465 target.fContaminationEnhancement = dynamic_cast<TH1F*> (fContaminationEnhancement->Clone());
467 target.fCombineMinMax = fCombineMinMax;
468 target.fTrackEtaCut = fTrackEtaCut;
469 target.fHistogramType = fHistogramType;
472 //____________________________________________________________________
473 Long64_t AliUEHist::Merge(TCollection* list)
475 // Merge a list of AliUEHist objects with this (needed for
477 // Returns the number of merged objects (including this).
485 TIterator* iter = list->MakeIterator();
488 // collections of objects
489 const UInt_t kMaxLists = fkRegions+3;
490 TList** lists = new TList*[kMaxLists];
492 for (UInt_t i=0; i<kMaxLists; i++)
493 lists[i] = new TList;
496 while ((obj = iter->Next())) {
498 AliUEHist* entry = dynamic_cast<AliUEHist*> (obj);
502 for (UInt_t i=0; i<fkRegions; i++)
503 if (entry->fTrackHist[i])
504 lists[i]->Add(entry->fTrackHist[i]);
506 lists[fkRegions]->Add(entry->fEventHist);
507 lists[fkRegions+1]->Add(entry->fTrackHistEfficiency);
509 lists[fkRegions+2]->Add(entry->fFakePt);
513 for (UInt_t i=0; i<fkRegions; i++)
515 fTrackHist[i]->Merge(lists[i]);
517 fEventHist->Merge(lists[fkRegions]);
518 fTrackHistEfficiency->Merge(lists[fkRegions+1]);
520 fFakePt->Merge(lists[fkRegions+2]);
522 for (UInt_t i=0; i<kMaxLists; i++)
530 //____________________________________________________________________
531 void AliUEHist::SetBinLimits(THnBase* grid)
533 // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
535 if (fEtaMax > fEtaMin)
536 grid->GetAxis(0)->SetRangeUser(fEtaMin, fEtaMax);
538 grid->GetAxis(1)->SetRangeUser(fPtMin, fPtMax);
541 //____________________________________________________________________
542 void AliUEHist::SetBinLimits(AliCFGridSparse* grid)
544 // sets the bin limits in eta and pT defined by fEtaMin/Max, fPtMin/Max
546 SetBinLimits(grid->GetGrid());
549 //____________________________________________________________________
550 void AliUEHist::ResetBinLimits(THnBase* grid)
552 // resets all bin limits
554 for (Int_t i=0; i<grid->GetNdimensions(); i++)
555 if (grid->GetAxis(i)->TestBit(TAxis::kAxisRange))
556 grid->GetAxis(i)->SetRangeUser(0, -1);
559 //____________________________________________________________________
560 void AliUEHist::ResetBinLimits(AliCFGridSparse* grid)
562 // resets all bin limits
564 ResetBinLimits(grid->GetGrid());
567 //____________________________________________________________________
568 void AliUEHist::CountEmptyBins(AliUEHist::CFStep step, Float_t ptLeadMin, Float_t ptLeadMax)
570 // prints the number of empty bins in the track end event histograms in the given step
575 for (Int_t i=0; i<4; i++)
578 binEnd[i] = fTrackHist[0]->GetNBins(i);
581 if (fEtaMax > fEtaMin)
583 binBegin[0] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMin);
584 binEnd[0] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(0)->FindBin(fEtaMax);
589 binBegin[1] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMin);
590 binEnd[1] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(1)->FindBin(fPtMax);
593 if (ptLeadMax > ptLeadMin)
595 binBegin[2] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMin);
596 binEnd[2] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(2)->FindBin(ptLeadMax);
599 // start from multiplicity 1
600 binBegin[3] = fTrackHist[0]->GetGrid(step)->GetGrid()->GetAxis(3)->FindBin(1);
602 for (UInt_t region=0; region<fkRegions; region++)
608 for (Int_t i=0; i<4; i++)
609 vars[i] = binBegin[i];
611 AliCFGridSparse* grid = fTrackHist[region]->GetGrid(step);
614 if (grid->GetElement(vars) == 0)
616 Printf("Empty bin at eta=%.2f pt=%.2f pt_lead=%.2f mult=%.1f",
617 grid->GetBinCenter(0, vars[0]),
618 grid->GetBinCenter(1, vars[1]),
619 grid->GetBinCenter(2, vars[2]),
620 grid->GetBinCenter(3, vars[3])
626 for (Int_t i=3; i>0; i--)
628 if (vars[i] == binEnd[i]+1)
630 vars[i] = binBegin[i];
635 if (vars[0] == binEnd[0]+1)
640 Printf("Region %s has %d empty bins (out of %d bins)", GetRegionTitle((Region) region), count, total);
644 //____________________________________________________________________
645 TH1* AliUEHist::GetUEHist(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Int_t twoD, Bool_t etaNorm, Int_t* normEvents)
647 // Extracts the UE histogram at the given step and in the given region by projection and dividing tracks by events
649 // ptLeadMin, ptLeadMax: Only meaningful for vs. delta phi plot (third axis is ptLead)
650 // Histogram has to be deleted by the caller of the function
652 // twoD: 0: 1D histogram as function of phi
653 // 1: 2D histogram, deltaphi, deltaeta
654 // 10: 1D histogram, within |deltaeta| < 0.8
655 // 11: 1D histogram, within 0.8 < |deltaeta| < 1.6
657 // etaNorm: if kTRUE (default), the distributions are divided by the area in delta eta
659 // normEvents: if non-0 the number of events/trigger particles for the normalization is filled
662 ResetBinLimits(fTrackHist[region]->GetGrid(step));
663 ResetBinLimits(fEventHist->GetGrid(step));
665 SetBinLimits(fTrackHist[region]->GetGrid(step));
668 if (fZVtxMax > fZVtxMin)
670 Printf("Using z-vtx range %f --> %f", fZVtxMin, fZVtxMax);
671 fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(5)->SetRangeUser(fZVtxMin, fZVtxMax);
672 fEventHist->GetGrid(step)->GetGrid()->GetAxis(2)->SetRangeUser(fZVtxMin, fZVtxMax);
679 tracks = fTrackHist[region]->ShowProjection(2, step);
680 tracks->GetYaxis()->SetTitle(fTrackHist[region]->GetTitle());
681 if (fCombineMinMax && region == kMin)
683 ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
684 SetBinLimits(fTrackHist[kMax]->GetGrid(step));
686 TH1D* tracks2 = fTrackHist[kMax]->ShowProjection(2, step);
687 tracks->Add(tracks2);
689 ResetBinLimits(fTrackHist[kMax]->GetGrid(step));
692 // normalize to get a density (deta dphi)
693 TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
694 Float_t phiRegion = TMath::TwoPi() / 3;
695 if (!fCombineMinMax && region == kMin)
697 Float_t normalization = phiRegion;
699 normalization *= (axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst()));
700 //Printf("Normalization: %f", normalization);
701 tracks->Scale(1.0 / normalization);
703 TH1D* events = fEventHist->ShowProjection(0, step);
704 tracks->Divide(events);
708 if (multBinEnd >= multBinBegin)
710 Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
711 fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
712 fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
715 Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
716 Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
718 Printf("Using leading pT range %d --> %d", firstBin, lastBin);
720 fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(firstBin, lastBin);
723 tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4);
725 tracks = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(4, 0);
727 Printf("Calculated histogram --> %f tracks", tracks->Integral());
728 fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
730 if (twoD == 10 || twoD == 11)
732 Float_t etaLimit = 1.0;
735 tracks = (TH1D*) ((TH2*) tracks)->ProjectionX("proj", tracks->GetYaxis()->FindBin(-etaLimit + 0.01), tracks->GetYaxis()->FindBin(etaLimit - 0.01))->Clone();
737 // TODO calculate acc with 2 * (deta - 0.5 * deta*deta / 1.6)
738 tracks->Scale(1. / 0.75);
742 TH1D* tracksTmp1 = (TH1D*) ((TH2*) tracks)->ProjectionX("proj1", tracks->GetYaxis()->FindBin(-etaLimit * 2 + 0.01), tracks->GetYaxis()->FindBin(-etaLimit - 0.01))->Clone();
743 TH1D* tracksTmp2 = ((TH2*) tracks)->ProjectionX("proj2", tracks->GetYaxis()->FindBin(etaLimit + 0.01), tracks->GetYaxis()->FindBin(2 * etaLimit - 0.01));
745 tracksTmp1->Add(tracksTmp2);
751 tracks->Scale(1. / 0.25);
755 // normalize to get a density (deta dphi)
756 // TODO normalization may be off for 2d histogram
757 Float_t normalization = fTrackHist[region]->GetGrid(step)->GetAxis(4)->GetBinWidth(1);
761 TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
762 if (strcmp(axis->GetTitle(), "#eta") == 0)
764 Printf("Normalizing using eta-axis range");
765 normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
769 Printf("Normalizing assuming accepted range of |eta| < 0.8");
770 normalization *= 0.8 * 2;
774 //Printf("Normalization: %f", normalization);
775 tracks->Scale(1.0 / normalization);
777 // NOTE fEventHist contains the number of events for the underlying event analysis and the number of trigger particles for the azimuthal correlation analysis. In the latter case the naming is therefore somewhat misleading!
778 TH1D* events = fEventHist->ShowProjection(0, step);
779 Int_t nEvents = (Int_t) events->Integral(firstBin, lastBin);
780 Printf("Calculated histogram --> %d events", nEvents);
782 *normEvents = nEvents;
785 tracks->Scale(1.0 / nEvents);
790 ResetBinLimits(fTrackHist[region]->GetGrid(step));
791 ResetBinLimits(fEventHist->GetGrid(step));
796 //____________________________________________________________________
797 void AliUEHist::GetHistsZVtx(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, TH3** trackHist, TH1** eventHist)
799 // Calculates a 3d histogram with deltaphi, deltaeta, zvtx on track level and a 1d histogram on event level (as fct of zvtx)
800 // Histogram has to be deleted by the caller of the function
802 if (!fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5))
803 AliFatal("Histogram without vertex axis provided");
806 ResetBinLimits(fTrackHist[region]->GetGrid(step));
807 ResetBinLimits(fEventHist->GetGrid(step));
809 SetBinLimits(fTrackHist[region]->GetGrid(step));
811 if (multBinEnd >= multBinBegin)
813 Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
814 fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
815 fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
818 Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
819 Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
820 Printf("Using leading pT range %d --> %d", firstBin, lastBin);
821 fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(firstBin, lastBin);
822 fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
824 *trackHist = (TH3*) fTrackHist[region]->GetGrid(step)->Project(4, 0, 5);
825 *eventHist = (TH1*) fEventHist->GetGrid(step)->Project(2);
827 ResetBinLimits(fTrackHist[region]->GetGrid(step));
828 ResetBinLimits(fEventHist->GetGrid(step));
831 //____________________________________________________________________
832 void AliUEHist::GetHistsZVtxMult(AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, THnBase** trackHist, TH2** eventHist)
834 // Calculates a 4d histogram with deltaphi, deltaeta, zvtx, multiplicity on track level and
835 // a 2d histogram on event level (as fct of zvtx, multiplicity)
836 // Histograms has to be deleted by the caller of the function
838 if (!fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5))
839 AliFatal("Histogram without vertex axis provided");
841 THnBase* sparse = fTrackHist[region]->GetGrid(step)->GetGrid();
846 fGetMultCache = ChangeToThn(sparse);
847 // should work but causes SEGV in ProjectionND below
850 sparse = fGetMultCache;
854 ResetBinLimits(sparse);
855 ResetBinLimits(fEventHist->GetGrid(step));
857 SetBinLimits(sparse);
859 Int_t firstBin = sparse->GetAxis(2)->FindBin(ptLeadMin);
860 Int_t lastBin = sparse->GetAxis(2)->FindBin(ptLeadMax);
861 Printf("Using leading pT range %d --> %d", firstBin, lastBin);
862 sparse->GetAxis(2)->SetRange(firstBin, lastBin);
863 fEventHist->GetGrid(step)->GetGrid()->GetAxis(0)->SetRange(firstBin, lastBin);
865 Int_t dimensions[] = { 4, 0, 5, 3 };
866 THnBase* tmpTrackHist = sparse->ProjectionND(4, dimensions, "E");
867 *eventHist = (TH2*) fEventHist->GetGrid(step)->Project(2, 1);
869 ResetBinLimits(sparse);
870 ResetBinLimits(fEventHist->GetGrid(step));
873 *trackHist = ChangeToThn(tmpTrackHist);
876 //____________________________________________________________________
877 TH2* AliUEHist::GetSumOfRatios2(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd)
879 // Calls GetUEHist(...) for *each* vertex bin and multiplicity bin and performs a sum of ratios:
880 // 1_N [ (same/mixed)_1 + (same/mixed)_2 + (same/mixed)_3 + ... ]
881 // where N is the total number of events/trigger particles and the subscript is the vertex/multiplicity bin
882 // where mixed is normalized such that the information about the number of pairs in same is kept
884 // returns a 2D histogram: deltaphi, deltaeta
887 // mixed: AliUEHist containing mixed event corresponding to this object
888 // <other parameters> : check documentation of AliUEHist::GetUEHist
890 // do not add this hists to the directory
891 Bool_t oldStatus = TH1::AddDirectoryStatus();
892 TH1::AddDirectory(kFALSE);
894 TH2* totalTracks = 0;
896 THnBase* trackSameAll = 0;
897 THnBase* trackMixedAll = 0;
898 TH2* eventSameAll = 0;
899 TH2* eventMixedAll = 0;
901 Int_t totalEvents = 0;
902 Int_t nCorrelationFunctions = 0;
904 GetHistsZVtxMult(step, region, ptLeadMin, ptLeadMax, &trackSameAll, &eventSameAll);
905 mixed->GetHistsZVtxMult(step, region, ptLeadMin, ptLeadMax, &trackMixedAll, &eventMixedAll);
907 // TH1* normParameters = new TH1F("normParameters", "", 100, 0, 2);
909 // trackSameAll->Dump();
911 TAxis* multAxis = trackSameAll->GetAxis(3);
913 if (multBinEnd < multBinBegin)
916 multBinEnd = multAxis->GetNbins();
919 for (Int_t multBin = TMath::Max(1, multBinBegin); multBin <= TMath::Min(multAxis->GetNbins(), multBinEnd); multBin++)
921 trackSameAll->GetAxis(3)->SetRange(multBin, multBin);
922 trackMixedAll->GetAxis(3)->SetRange(multBin, multBin);
924 // get mixed normalization correction factor: is independent of vertex bin if scaled with number of triggers
925 trackMixedAll->GetAxis(2)->SetRange(0, -1);
926 TH2* tracksMixed = trackMixedAll->Projection(1, 0, "E");
927 Float_t binWidthEta = tracksMixed->GetYaxis()->GetBinWidth(1);
929 // get mixed event normalization by assuming full acceptance at deta at 0 (integrate over dphi), excluding (0, 0)
930 Double_t mixedNormError = 0;
931 Double_t mixedNorm = tracksMixed->IntegralAndError(1, tracksMixed->GetYaxis()->FindBin(-0.01)-1, tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError);
932 Double_t mixedNormError2 = 0;
933 Double_t mixedNorm2 = tracksMixed->IntegralAndError(tracksMixed->GetYaxis()->FindBin(0.01)+1, tracksMixed->GetNbinsX(), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError2);
935 if (mixedNormError == 0 || mixedNormError2 == 0)
937 Printf("ERROR: Skipping multiplicity %d because mixed event is empty", multBin);
941 Int_t nBinsMixedNorm = (tracksMixed->GetYaxis()->FindBin(-0.01) - 1 - 1 + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);
942 mixedNorm /= nBinsMixedNorm;
943 mixedNormError /= nBinsMixedNorm;
945 Int_t nBinsMixedNorm2 = (tracksMixed->GetNbinsX() - tracksMixed->GetYaxis()->FindBin(0.01) - 1 + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);
946 mixedNorm2 /= nBinsMixedNorm2;
947 mixedNormError2 /= nBinsMixedNorm2;
949 mixedNorm = mixedNorm / mixedNormError / mixedNormError + mixedNorm2 / mixedNormError2 / mixedNormError2;
950 mixedNormError = TMath::Sqrt(1.0 / (1.0 / mixedNormError / mixedNormError + 1.0 / mixedNormError2 / mixedNormError2));
951 mixedNorm *= mixedNormError * mixedNormError;
953 /* Double_t mixedNorm = tracksMixed->IntegralAndError(tracksMixed->GetXaxis()->FindBin(-0.01), tracksMixed->GetXaxis()->FindBin(0.01), tracksMixed->GetYaxis()->FindBin(-0.01), tracksMixed->GetYaxis()->FindBin(0.01), mixedNormError);
954 Int_t nBinsMixedNorm = (tracksMixed->GetXaxis()->FindBin(0.01) - tracksMixed->GetXaxis()->FindBin(-0.01) + 1) * (tracksMixed->GetYaxis()->FindBin(0.01) - tracksMixed->GetYaxis()->FindBin(-0.01) + 1);*/
958 Float_t triggers = eventMixedAll->Integral(1, eventMixedAll->GetNbinsX(), multBin, multBin);
959 // Printf("%f +- %f | %f | %f", mixedNorm, mixedNormError, triggers, mixedNorm / triggers);
962 Printf("ERROR: Skipping multiplicity %d because mixed event is empty", multBin);
966 mixedNorm /= triggers;
967 mixedNormError /= triggers;
971 Printf("ERROR: Skipping multiplicity %d because mixed event is empty at (0,0)", multBin);
975 // finite bin correction
976 if (fTrackEtaCut > 0)
978 Double_t finiteBinCorrection = -1.0 / (2*fTrackEtaCut) * binWidthEta / 2 + 1;
979 Printf("Finite bin correction: %f", finiteBinCorrection);
980 mixedNorm *= finiteBinCorrection;
981 mixedNormError *= finiteBinCorrection;
985 Printf("ERROR: fTrackEtaCut not set. Finite bin correction cannot be applied. Continuing anyway...");
988 // Printf("Norm: %f +- %f", mixedNorm, mixedNormError);
990 // normParameters->Fill(mixedNorm);
992 TAxis* vertexAxis = trackSameAll->GetAxis(2);
993 for (Int_t vertexBin = 1; vertexBin <= vertexAxis->GetNbins(); vertexBin++)
994 // for (Int_t vertexBin = 3; vertexBin <= 8; vertexBin++)
996 trackSameAll->GetAxis(2)->SetRange(vertexBin, vertexBin);
997 trackMixedAll->GetAxis(2)->SetRange(vertexBin, vertexBin);
999 TH2* tracksSame = trackSameAll->Projection(1, 0, "E");
1000 tracksMixed = trackMixedAll->Projection(1, 0, "E");
1002 // asssume flat in dphi, gain in statistics
1003 // TH1* histMixedproj = mixedTwoD->ProjectionY();
1004 // histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX());
1006 // for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++)
1007 // for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++)
1008 // mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y));
1010 // delete histMixedproj;
1012 Float_t triggers2 = eventMixedAll->Integral(vertexBin, vertexBin, multBin, multBin);
1015 Printf("ERROR: Skipping multiplicity %d vertex bin %d because mixed event is empty", multBin, vertexBin);
1019 tracksMixed->Scale(1.0 / triggers2 / mixedNorm);
1020 // tracksSame->Scale(tracksMixed->Integral() / tracksSame->Integral());
1022 // new TCanvas; tracksSame->DrawClone("SURF1");
1023 // new TCanvas; tracksMixed->DrawClone("SURF1");
1025 // some code to judge the relative contribution of the different correlation functions to the overall uncertainty
1026 Double_t sums[] = { 0, 0, 0 };
1027 Double_t errors[] = { 0, 0, 0 };
1029 for (Int_t x=1; x<=tracksSame->GetNbinsX(); x++)
1030 for (Int_t y=1; y<=tracksSame->GetNbinsY(); y++)
1032 sums[0] += tracksSame->GetBinContent(x, y);
1033 errors[0] += tracksSame->GetBinError(x, y);
1034 sums[1] += tracksMixed->GetBinContent(x, y);
1035 errors[1] += tracksMixed->GetBinError(x, y);
1038 tracksSame->Divide(tracksMixed);
1040 for (Int_t x=1; x<=tracksSame->GetNbinsX(); x++)
1041 for (Int_t y=1; y<=tracksSame->GetNbinsY(); y++)
1043 sums[2] += tracksSame->GetBinContent(x, y);
1044 errors[2] += tracksSame->GetBinError(x, y);
1047 for (Int_t x=0; x<3; x++)
1049 errors[x] /= sums[x];
1051 Printf("The correlation function %d %d has uncertainties %f %f %f (Ratio S/M %f)", multBin, vertexBin, errors[0], errors[1], errors[2], (errors[1] > 0) ? errors[0] / errors[1] : -1);
1052 // code to draw contributions
1054 TH1* proj = tracksSame->ProjectionX("projx", tracksSame->GetYaxis()->FindBin(-1.59), tracksSame->GetYaxis()->FindBin(1.59));
1055 proj->SetTitle(Form("Bin %d", vertexBin));
1056 proj->SetLineColor(vertexBin);
1057 proj->DrawCopy((vertexBin > 1) ? "SAME" : "");
1061 totalTracks = (TH2*) tracksSame->Clone("totalTracks");
1063 totalTracks->Add(tracksSame);
1065 totalEvents += eventSameAll->GetBinContent(vertexBin, multBin);
1067 // new TCanvas; tracksMixed->DrawCopy("SURF1");
1073 nCorrelationFunctions++;
1078 Double_t sums[] = { 0, 0, 0 };
1079 Double_t errors[] = { 0, 0, 0 };
1081 for (Int_t x=1; x<=totalTracks->GetNbinsX(); x++)
1082 for (Int_t y=1; y<=totalTracks->GetNbinsY(); y++)
1084 sums[0] += totalTracks->GetBinContent(x, y);
1085 errors[0] += totalTracks->GetBinError(x, y);
1088 errors[0] /= sums[0];
1090 Printf("Dividing %f tracks by %d events (%d correlation function(s)) (error %f)", totalTracks->Integral(), totalEvents, nCorrelationFunctions, errors[0]);
1091 if (totalEvents > 0)
1092 totalTracks->Scale(1.0 / totalEvents);
1094 // normalizate to dphi width
1095 Float_t normalization = totalTracks->GetXaxis()->GetBinWidth(1);
1096 totalTracks->Scale(1.0 / normalization);
1099 delete trackSameAll;
1100 delete trackMixedAll;
1101 delete eventSameAll;
1102 delete eventMixedAll;
1104 // new TCanvas; normParameters->Draw();
1106 TH1::AddDirectory(oldStatus);
1112 TH2* AliUEHist::GetPtDistInPhiRegion(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Float_t phiBegin, Float_t phiEnd)
1114 // Returns pT,assoc distribution in the given pt,trig, multiplicity, phi region
1115 // Does not use sum of ratios for mixed event correction (TODO to be improved)
1116 // returns a 2D histogram: deltaphi, deltaeta
1119 // mixed: AliUEHist containing mixed event corresponding to this object
1122 //____________________________________________________________________
1123 TH2* AliUEHist::GetSumOfRatios(AliUEHist* mixed, AliUEHist::CFStep step, AliUEHist::Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Bool_t etaNorm, Bool_t useVertexBins)
1125 // Calls GetUEHist(...) for *each* multiplicity bin and performs a sum of ratios:
1126 // 1_N [ (same/mixed)_1 + (same/mixed)_2 + (same/mixed)_3 + ... ]
1127 // where N is the total number of events/trigger particles and the subscript is the multiplicity bin
1129 // Can only be used for the 2D histogram at present
1132 // mixed: AliUEHist containing mixed event corresponding to this object
1133 // <other parameters> : check documentation of AliUEHist::GetUEHist
1135 TH2* totalTracks = 0;
1136 Int_t totalEvents = 0;
1138 Int_t vertexBin = 1;
1139 TAxis* vertexAxis = fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(5);
1140 if (useVertexBins && !vertexAxis)
1142 Printf("Vertex axis requested but not available");
1151 SetZVtxRange(vertexAxis->GetBinLowEdge(vertexBin) + 0.01, vertexAxis->GetBinUpEdge(vertexBin) - 0.01);
1152 mixed->SetZVtxRange(vertexAxis->GetBinLowEdge(vertexBin) + 0.01, vertexAxis->GetBinUpEdge(vertexBin) - 0.01);
1156 // multiplicity loop
1157 Int_t multIter = multBinBegin;
1160 Int_t multBinBeginLocal = multBinBegin;
1161 Int_t multBinEndLocal = multBinEnd;
1163 if (multBinEnd >= multBinBegin)
1165 multBinBeginLocal = multIter;
1166 multBinEndLocal = multIter;
1171 TH2* tracks = (TH2*) GetUEHist(step, region, ptLeadMin, ptLeadMax, multBinBeginLocal, multBinEndLocal, 1, etaNorm, &nEvents);
1172 // undo normalization
1173 tracks->Scale(nEvents);
1174 totalEvents += nEvents;
1176 TH2* mixedTwoD = (TH2*) mixed->GetUEHist(step, region, ptLeadMin, ptLeadMax, multBinBeginLocal, multBinEndLocal, 1, etaNorm);
1178 // asssume flat in dphi, gain in statistics
1179 // TH1* histMixedproj = mixedTwoD->ProjectionY();
1180 // histMixedproj->Scale(1.0 / mixedTwoD->GetNbinsX());
1182 // for (Int_t x=1; x<=mixedTwoD->GetNbinsX(); x++)
1183 // for (Int_t y=1; y<=mixedTwoD->GetNbinsY(); y++)
1184 // mixedTwoD->SetBinContent(x, y, histMixedproj->GetBinContent(y));
1186 // delete histMixedproj;
1188 // get mixed event normalization by assuming full acceptance at deta of 0 (only works for flat dphi)
1189 /* Double_t mixedNorm = mixedTwoD->Integral(1, mixedTwoD->GetNbinsX(), mixedTwoD->GetYaxis()->FindBin(-0.01), mixedTwoD->GetYaxis()->FindBin(0.01));
1190 mixedNorm /= mixedTwoD->GetNbinsX() * (mixedTwoD->GetYaxis()->FindBin(0.01) - mixedTwoD->GetYaxis()->FindBin(-0.01) + 1);
1191 tracks->Scale(mixedNorm);*/
1193 tracks->Scale(mixedTwoD->Integral() / tracks->Integral());
1195 tracks->Divide(mixedTwoD);
1200 totalTracks = tracks;
1203 totalTracks->Add(tracks);
1207 if (multIter > multBinEnd)
1211 if (!useVertexBins || vertexBin > vertexAxis->GetNbins())
1216 totalEvents = vertexAxis->GetNbins();
1218 Printf("Dividing %f tracks by %d events", totalTracks->Integral(), totalEvents);
1219 if (totalEvents > 0)
1220 totalTracks->Scale(1.0 / totalEvents);
1225 //____________________________________________________________________
1226 TH1* AliUEHist::GetPtHist(CFStep step, Region region, Float_t ptLeadMin, Float_t ptLeadMax, Int_t multBinBegin, Int_t multBinEnd, Float_t phiMin, Float_t phiMax, Float_t etaMin, Float_t etaMax, Bool_t skipPhiNormalization)
1228 // returns a histogram projected to pT,assoc with the provided cuts
1231 ResetBinLimits(fTrackHist[region]->GetGrid(step));
1232 ResetBinLimits(fEventHist->GetGrid(step));
1236 // the efficiency to have find an event depends on leading pT and this is not corrected for because anyway per bin we calculate tracks over events
1237 // therefore the number density must be calculated here per leading pT bin and then summed
1239 if (multBinEnd > multBinBegin)
1240 Printf("Using multiplicity range %d --> %d", multBinBegin, multBinEnd);
1241 fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(3)->SetRange(multBinBegin, multBinEnd);
1242 fEventHist->GetGrid(step)->GetGrid()->GetAxis(1)->SetRange(multBinBegin, multBinEnd);
1244 fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(0)->SetRangeUser(etaMin, etaMax);
1245 fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->SetRangeUser(phiMin, phiMax);
1247 // get real phi cuts due to binning
1248 Float_t phiMinReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinLowEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetFirst());
1249 Float_t phiMaxReal = fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetBinUpEdge(fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(4)->GetLast());
1250 Printf("phi min = %.2f and max = %.2f requested; due to binning range will be from %.2f to %.2f, i.e. a %.0f%% larger window", phiMin, phiMax, phiMinReal, phiMaxReal, (phiMaxReal - phiMinReal) / (phiMax - phiMin) * 100 - 100);
1252 Int_t firstBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMin);
1253 Int_t lastBin = fTrackHist[region]->GetAxis(2, step)->FindBin(ptLeadMax);
1255 TH1D* events = fEventHist->ShowProjection(0, step);
1257 for (Int_t bin=firstBin; bin<=lastBin; bin++)
1259 fTrackHist[region]->GetGrid(step)->GetGrid()->GetAxis(2)->SetRange(bin, bin);
1261 // project to pT,assoc
1262 TH1D* tracksTmp = (TH1D*) fTrackHist[region]->GetGrid(step)->Project(1);
1264 Printf("Calculated histogram in bin %d --> %f tracks", bin, tracksTmp->Integral());
1265 fTrackHist[region]->GetGrid(step)->SetRangeUser(2, 0, -1);
1267 // normalize to get a density (deta dphi)
1268 Float_t normalization = 1;
1269 TAxis* axis = fTrackHist[region]->GetGrid(step)->GetAxis(0);
1270 if (strcmp(axis->GetTitle(), "#eta") == 0)
1272 Printf("Normalizing using eta-axis range");
1273 normalization *= axis->GetBinUpEdge(axis->GetLast()) - axis->GetBinLowEdge(axis->GetFirst());
1277 Printf("Normalizating assuming accepted range of |eta| < 0.8");
1278 normalization *= 0.8 * 2;
1282 if (!skipPhiNormalization)
1283 normalization *= phiMaxReal - phiMinReal;
1285 //Printf("Normalization: %f", normalization);
1286 tracksTmp->Scale(1.0 / normalization);
1288 // and now dpT (bins have different width)
1289 for (Int_t i=1; i<=tracksTmp->GetNbinsX(); i++)
1291 tracksTmp->SetBinContent(i, tracksTmp->GetBinContent(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
1292 tracksTmp->SetBinError (i, tracksTmp->GetBinError(i) / tracksTmp->GetXaxis()->GetBinWidth(i));
1295 Int_t nEvents = (Int_t) events->GetBinContent(bin);
1297 tracksTmp->Scale(1.0 / nEvents);
1298 Printf("Calculated histogram in bin %d --> %d events", bin, nEvents);
1304 tracks->Add(tracksTmp);
1311 ResetBinLimits(fTrackHist[region]->GetGrid(step));
1312 ResetBinLimits(fEventHist->GetGrid(step));
1317 void AliUEHist::MultiplyHistograms(THnSparse* grid, THnSparse* target, TH1* histogram, Int_t var1, Int_t var2)
1319 // multiplies <grid> with <histogram> and put the result in <target>
1320 // <grid> has usually more dimensions than histogram. The axis which are used to choose the value
1321 // from <histogram> are given in <var1> and <var2>
1323 // if <histogram> is 0, just copies content from step1 to step2
1325 // clear target histogram
1331 if (grid->GetAxis(var1)->GetNbins() != histogram->GetNbinsX())
1332 AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), histogram->GetNbinsX()));
1334 if (var2 >= 0 && grid->GetAxis(var2)->GetNbins() != histogram->GetNbinsY())
1335 AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), histogram->GetNbinsY()));
1338 if (grid->GetNdimensions() > 6)
1339 AliFatal("Too many dimensions in THnSparse");
1343 // optimized implementation
1344 for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
1346 Double_t value = grid->GetBinContent(binIdx, bins);
1347 Double_t error = grid->GetBinError(binIdx);
1353 value *= histogram->GetBinContent(bins[var1]);
1354 error *= histogram->GetBinContent(bins[var1]);
1358 value *= histogram->GetBinContent(bins[var1], bins[var2]);
1359 error *= histogram->GetBinContent(bins[var1], bins[var2]);
1363 target->SetBinContent(bins, value);
1364 target->SetBinError(bins, error);
1368 //____________________________________________________________________
1369 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, TH1* trackCorrection, Int_t var1, Int_t var2)
1371 // corrects from step1 to step2 by multiplying the tracks with trackCorrection
1372 // trackCorrection can be a function of eta (var1 == 0), pT (var1 == 1), leading pT (var1 == 2), multiplicity (var1 == 3), delta phi (var1 == 4)
1373 // if var2 >= 0 a two dimension correction is assumed in trackCorrection
1375 // if trackCorrection is 0, just copies content from step1 to step2
1377 for (UInt_t region=0; region<fkRegions; region++)
1378 CorrectTracks(step1, step2, region, trackCorrection, var1, var2);
1381 //____________________________________________________________________
1382 void AliUEHist::CorrectTracks(CFStep step1, CFStep step2, Int_t region, TH1* trackCorrection, Int_t var1, Int_t var2)
1385 // see documentation of CorrectTracks above
1388 if (!fTrackHist[region])
1391 THnSparse* grid = fTrackHist[region]->GetGrid(step1)->GetGrid();
1392 THnSparse* target = fTrackHist[region]->GetGrid(step2)->GetGrid();
1394 Float_t entriesBefore = grid->GetEntries();
1396 MultiplyHistograms(grid, target, trackCorrection, var1, var2);
1398 Printf("AliUEHist::CorrectTracks: Corrected from %f (step %d) to %f (step %d) entries. Correction histogram: %f entries (integral: %f)", entriesBefore, step1, target->GetEntries(), step2, (trackCorrection) ? trackCorrection->GetEntries() : -1.0, (trackCorrection) ? trackCorrection->Integral() : -1.0);
1401 //____________________________________________________________________
1402 void AliUEHist::CorrectEvents(CFStep step1, CFStep step2, TH1* eventCorrection, Int_t var1, Int_t var2)
1404 // corrects from step1 to step2 by multiplying the events with eventCorrection
1405 // eventCorrection is as function of leading pT (var == 0) or multiplicity (var == 1)
1407 // if eventCorrection is 0, just copies content from step1 to step2
1409 AliCFGridSparse* grid = fEventHist->GetGrid(step1);
1410 AliCFGridSparse* target = fEventHist->GetGrid(step2);
1412 Float_t entriesBefore = grid->GetEntries();
1414 MultiplyHistograms(grid->GetGrid(), target->GetGrid(), eventCorrection, var1, var2);
1416 Printf("AliUEHist::CorrectEvents: Corrected from %f (step %d) to %f (step %d) entries. Correction histogram: %f entries (integral: %f)", entriesBefore, step1, target->GetEntries(), step2, (eventCorrection) ? eventCorrection->GetEntries() : -1.0, (eventCorrection) ? eventCorrection->Integral() : -1.0);
1419 //____________________________________________________________________
1420 void AliUEHist::Correct(AliUEHist* corrections)
1422 // applies the given corrections to extract from the step kCFStepReconstructed all previous steps
1424 // in this object the data is expected in the step kCFStepReconstructed
1426 if (fHistogramType.Length() == 0)
1428 Printf("W-AliUEHist::Correct: fHistogramType not defined. Guessing histogram type...");
1429 if (fTrackHist[kToward]->GetNVar() < 5)
1431 if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}N_{ch}/d#varphid#eta") == 0)
1432 fHistogramType = "NumberDensitypT";
1433 else if (strcmp(fTrackHist[kToward]->GetTitle(), "d^{2}#Sigma p_{T}/d#varphid#eta") == 0)
1434 fHistogramType = "SumpT";
1436 else if (fTrackHist[kToward]->GetNVar() == 5)
1438 if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "multiplicity") == 0)
1439 fHistogramType = "NumberDensityPhi";
1440 else if (strcmp(fTrackHist[kToward]->GetGrid(0)->GetGrid()->GetAxis(3)->GetTitle(), "centrality") == 0)
1441 fHistogramType = "NumberDensityPhiCentrality";
1444 if (fHistogramType.Length() == 0)
1445 AliFatal("Cannot figure out histogram type. Try setting it manually...");
1448 Printf("AliUEHist::Correct: Correcting %s...", fHistogramType.Data());
1450 if (strcmp(fHistogramType, "NumberDensitypT") == 0 || strcmp(fHistogramType, "NumberDensityPhi") == 0 || strcmp(fHistogramType, "SumpT") == 0)
1454 // bias due to migration in leading pT (because the leading particle is not reconstructed, and the subleading is used)
1455 // extracted as function of leading pT
1456 Bool_t biasFromMC = kFALSE;
1457 const char* projAxis = "z";
1458 Int_t secondBin = -1;
1460 if (strcmp(fHistogramType, "NumberDensityPhi") == 0)
1466 for (UInt_t region = 0; region < fkRegions; region++)
1468 if (!fTrackHist[region])
1471 TH1* leadingBiasTracks = 0;
1474 leadingBiasTracks = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, region, projAxis, 0, -1, 1); // from MC
1475 Printf("WARNING: Using MC bias correction");
1478 leadingBiasTracks = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, region, projAxis, 0, -1, 1); // from data
1480 CorrectTracks(kCFStepReconstructed, kCFStepTracked, region, leadingBiasTracks, 2, secondBin);
1481 if (region == kMin && fCombineMinMax)
1483 CorrectTracks(kCFStepReconstructed, kCFStepTracked, kMax, leadingBiasTracks, 2, secondBin);
1484 delete leadingBiasTracks;
1487 delete leadingBiasTracks;
1490 TH1* leadingBiasEvents = 0;
1492 leadingBiasEvents = (TH1*) corrections->GetBias(kCFStepReconstructed, kCFStepTracked, kToward, projAxis, 0, -1, 2); // from MC
1494 leadingBiasEvents = (TH1*) GetBias(kCFStepBiasStudy, kCFStepReconstructed, kToward, projAxis, 0, -1, 2); // from data
1496 //new TCanvas; leadingBiasEvents->DrawCopy();
1497 CorrectEvents(kCFStepReconstructed, kCFStepTracked, leadingBiasEvents, 0);
1498 delete leadingBiasEvents;
1500 // --- contamination correction ---
1502 TH2D* contamination = corrections->GetTrackingContamination();
1503 if (corrections->fContaminationEnhancement)
1505 Printf("Applying contamination enhancement");
1507 for (Int_t x=1; x<=contamination->GetNbinsX(); x++)
1508 for (Int_t y=1; y<=contamination->GetNbinsX(); y++)
1509 contamination->SetBinContent(x, y, contamination->GetBinContent(x, y) * corrections->fContaminationEnhancement->GetBinContent(corrections->fContaminationEnhancement->GetXaxis()->FindBin(contamination->GetYaxis()->GetBinCenter(y))));
1511 CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 0, 1);
1512 CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
1513 delete contamination;
1515 // --- efficiency correction ---
1516 // correct single-particle efficiency for associated particles
1517 // in addition correct for efficiency on trigger particles (tracks AND events)
1519 TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrection();
1520 // use kCFStepVertex as a temporary step (will be overwritten later anyway)
1521 CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
1522 delete efficiencyCorrection;
1525 efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, -1, 2);
1526 CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 0);
1527 CorrectTracks(kCFStepVertex, kCFStepAnaTopology, efficiencyCorrection, 2);
1528 delete efficiencyCorrection;
1531 CorrectTracks(kCFStepAnaTopology, kCFStepVertex, 0, -1);
1532 CorrectEvents(kCFStepAnaTopology, kCFStepVertex, 0, 0);
1534 // vertex correction on the level of events as function of multiplicity, weighting tracks and events with the same factor
1535 // practically independent of low pT cut
1536 TH1D* vertexCorrection = (TH1D*) corrections->GetEventEfficiency(kCFStepVertex, kCFStepTriggered, 1);
1538 // convert stage from true multiplicity to observed multiplicity by simple conversion factor
1539 TH1D* vertexCorrectionObs = (TH1D*) vertexCorrection->Clone("vertexCorrection2");
1540 vertexCorrectionObs->Reset();
1542 TF1* func = new TF1("func", "[1]+[0]/(x-[2])");
1544 func->SetParameters(0.1, 1, -0.7);
1545 vertexCorrection->Fit(func, "0I", "", 0, 3);
1546 for (Int_t i=1; i<=vertexCorrectionObs->GetNbinsX(); i++)
1548 Float_t xPos = 1.0 / 0.77 * vertexCorrectionObs->GetXaxis()->GetBinCenter(i);
1550 vertexCorrectionObs->SetBinContent(i, func->Eval(xPos));
1552 vertexCorrectionObs->SetBinContent(i, vertexCorrection->Interpolate(xPos));
1557 vertexCorrection->DrawCopy();
1558 vertexCorrectionObs->SetLineColor(2);
1559 vertexCorrectionObs->DrawCopy("same");
1560 func->SetRange(0, 4);
1561 func->DrawClone("same");
1564 CorrectTracks(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 3);
1565 CorrectEvents(kCFStepVertex, kCFStepTriggered, vertexCorrectionObs, 1);
1566 delete vertexCorrectionObs;
1567 delete vertexCorrection;
1571 CorrectTracks(kCFStepTriggered, kCFStepAll, 0, -1);
1572 CorrectEvents(kCFStepTriggered, kCFStepAll, 0, 0);
1574 else if (strcmp(fHistogramType, "NumberDensityPhiCentrality") == 0)
1576 if (fTrackHist[0]->GetNVar() <= 5)
1578 // do corrections copying between steps
1579 // CFStep step = kCFStepReconstructed;
1580 CFStep step = kCFStepBiasStudy;
1583 CorrectTracks(step, kCFStepTracked, 0, -1);
1584 CorrectEvents(step, kCFStepTracked, 0, -1);
1586 // Dont use eta in the following, because it is a Delta-eta axis
1588 // contamination correction
1589 // correct single-particle contamination for associated particles
1591 TH1* contamination = corrections->GetTrackingContamination(1);
1595 Printf("Applying contamination enhancement");
1597 for (Int_t bin = 1; bin <= contamination->GetNbinsX(); bin++)
1599 printf("%f", contamination->GetBinContent(bin));
1600 if (contamination->GetBinContent(bin) > 0)
1601 contamination->SetBinContent(bin, 1.0 + 1.1 * (contamination->GetBinContent(bin) - 1.0));
1602 printf(" --> %f\n", contamination->GetBinContent(bin));
1606 CorrectTracks(kCFStepTracked, kCFStepTrackedOnlyPrim, contamination, 1);
1607 delete contamination;
1609 // correct for additional contamination due to trigger particle around phi ~ 0
1610 TH2* correlatedContamination = corrections->GetCorrelatedContamination();
1613 Printf("Applying contamination enhancement");
1615 for (Int_t bin = 1; bin <= correlatedContamination->GetNbinsX(); bin++)
1616 for (Int_t bin2 = 1; bin2 <= correlatedContamination->GetNbinsY(); bin2++)
1618 printf("%f", correlatedContamination->GetBinContent(bin, bin2));
1619 if (correlatedContamination->GetBinContent(bin, bin2) > 0)
1620 correlatedContamination->SetBinContent(bin, bin2, 1.0 + 1.1 * (correlatedContamination->GetBinContent(bin, bin2) - 1.0));
1621 printf(" --> %f\n", correlatedContamination->GetBinContent(bin, bin2));
1625 // new TCanvas; correlatedContamination->DrawCopy("COLZ");
1626 CorrectCorrelatedContamination(kCFStepTrackedOnlyPrim, 0, correlatedContamination);
1627 // Printf("\n\n\nWARNING ---> SKIPPING CorrectCorrelatedContamination\n\n\n");
1629 delete correlatedContamination;
1631 // TODO correct for contamination of trigger particles (for tracks AND events)
1632 CorrectEvents(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 0);
1634 // --- efficiency correction ---
1635 // correct single-particle efficiency for associated particles
1636 // in addition correct for efficiency on trigger particles (tracks AND events)
1638 // in bins of pT and centrality
1639 TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrectionCentrality();
1640 // new TCanvas; efficiencyCorrection->DrawCopy("COLZ");
1641 // use kCFStepAnaTopology as a temporary step
1642 CorrectTracks(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, efficiencyCorrection, 1, 3);
1643 delete efficiencyCorrection;
1645 // correct pT,T in bins of pT and centrality
1646 efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3, 2);
1647 CorrectEvents(kCFStepTrackedOnlyPrim, kCFStepVertex, efficiencyCorrection, 0, 1);
1648 CorrectTracks(kCFStepAnaTopology, kCFStepVertex, efficiencyCorrection, 2, 3);
1649 delete efficiencyCorrection;
1651 // no correction for vertex finding efficiency and trigger efficiency needed in PbPb
1653 CorrectTracks(kCFStepVertex, kCFStepAll, 0, -1);
1654 CorrectEvents(kCFStepVertex, kCFStepAll, 0, -1);
1658 // with 6 axes there is not enough memory, do the corrections in-place
1659 Printf(">>>>>>>> Applying corrections in place to reduce memory consumption");
1660 CFStep step = kCFStepBiasStudy;
1661 // CFStep step = kCFStepReconstructed;
1663 // Dont use eta in the following, because it is a Delta-eta axis
1665 // --- contamination correction ---
1666 // correct single-particle contamination for associated particles
1667 // correct contamination for trigger particles (tracks AND events)
1670 TH1* contamination = corrections->GetTrackingContamination(1);
1674 Printf("Applying contamination enhancement");
1676 for (Int_t bin = 1; bin <= contamination->GetNbinsX(); bin++)
1678 printf("%f", contamination->GetBinContent(bin));
1679 if (contamination->GetBinContent(bin) > 0)
1680 contamination->SetBinContent(bin, 1.0 + 1.1 * (contamination->GetBinContent(bin) - 1.0));
1681 printf(" --> %f\n", contamination->GetBinContent(bin));
1685 // correct pT,A in bins of pT
1686 CorrectTracks(step, step, contamination, 1);
1687 delete contamination;
1689 // correct pT,T in bins of pT (for tracks AND events)
1690 contamination = corrections->GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 1, -1, 2);
1691 new TCanvas; contamination->DrawCopy();
1692 CorrectEvents(step, step, contamination, 0);
1693 CorrectTracks(step, step, contamination, 2);
1694 delete contamination;
1696 // correct for additional contamination due to trigger particle around phi ~ 0
1699 TH2* correlatedContamination = corrections->GetCorrelatedContamination();
1702 Printf("Applying contamination enhancement");
1704 for (Int_t bin = 1; bin <= correlatedContamination->GetNbinsX(); bin++)
1705 for (Int_t bin2 = 1; bin2 <= correlatedContamination->GetNbinsY(); bin2++)
1707 printf("%f", correlatedContamination->GetBinContent(bin, bin2));
1708 if (correlatedContamination->GetBinContent(bin, bin2) > 0)
1709 correlatedContamination->SetBinContent(bin, bin2, 1.0 + 1.1 * (correlatedContamination->GetBinContent(bin, bin2) - 1.0));
1710 printf(" --> %f\n", correlatedContamination->GetBinContent(bin, bin2));
1714 // new TCanvas; correlatedContamination->DrawCopy("COLZ");
1715 CorrectCorrelatedContamination(step, 0, correlatedContamination);
1717 delete correlatedContamination;
1720 Printf("\n\n\nWARNING ---> SKIPPING CorrectCorrelatedContamination\n\n\n");
1722 // --- tracking efficiency correction ---
1723 // correct single-particle efficiency for associated particles
1724 // correct for efficiency on trigger particles (tracks AND events)
1726 // in bins of pT and centrality
1727 TH1* efficiencyCorrection = corrections->GetTrackingEfficiencyCorrectionCentrality();
1728 // new TCanvas; efficiencyCorrection->DrawCopy("COLZ");
1730 // correct pT,A in bins of pT and centrality
1731 CorrectTracks(step, step, efficiencyCorrection, 1, 3);
1732 delete efficiencyCorrection;
1734 // correct pT,T in bins of pT and centrality
1735 efficiencyCorrection = corrections->GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3, 2);
1736 CorrectEvents(step, step, efficiencyCorrection, 0, 1);
1737 CorrectTracks(step, step, efficiencyCorrection, 2, 3);
1739 delete efficiencyCorrection;
1743 AliFatal(Form("Unknown histogram for correction: %s", GetTitle()));
1746 //____________________________________________________________________
1747 THnBase* AliUEHist::GetTrackEfficiencyND(CFStep step1, CFStep step2)
1749 // creates a track-level efficiency by dividing step2 by step1
1750 // in all dimensions but the particle species one
1752 AliCFContainer* sourceContainer = fTrackHistEfficiency;
1753 // step offset because we start with kCFStepAnaTopology
1754 step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
1755 step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
1757 ResetBinLimits(sourceContainer->GetGrid(step1));
1758 ResetBinLimits(sourceContainer->GetGrid(step2));
1760 if (fEtaMax > fEtaMin)
1762 Printf("Restricted eta-range to %f %f", fEtaMin, fEtaMax);
1763 sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
1764 sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
1767 Int_t dimensions[] = { 0, 1, 3, 4 };
1768 THnBase* generated = sourceContainer->GetGrid(step1)->GetGrid()->ProjectionND(4, dimensions);
1769 THnBase* measured = sourceContainer->GetGrid(step2)->GetGrid()->ProjectionND(4, dimensions);
1771 // Printf("%d %d %f %f", step1, step2, generated->GetEntries(), measured->GetEntries());
1773 ResetBinLimits(sourceContainer->GetGrid(step1));
1774 ResetBinLimits(sourceContainer->GetGrid(step2));
1776 measured->Divide(measured, generated, 1, 1, "B");
1783 //____________________________________________________________________
1784 TH1* AliUEHist::GetTrackEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Int_t source, Int_t axis3)
1786 // creates a track-level efficiency by dividing step2 by step1
1787 // projected to axis1 and axis2 (optional if >= 0)
1789 // source: 0 = fTrackHist; 1 = fTrackHistEfficiency; 2 = fTrackHistEfficiency rebinned for pT,T / pT,lead binning
1791 // integrate over regions
1792 // cache it for efficiency (usually more than one efficiency is requested)
1794 AliCFContainer* sourceContainer = 0;
1800 fCache = (AliCFContainer*) fTrackHist[0]->Clone();
1801 for (UInt_t i = 1; i < fkRegions; i++)
1803 fCache->Add(fTrackHist[i]);
1805 sourceContainer = fCache;
1807 else if (source == 1 || source == 2)
1809 sourceContainer = fTrackHistEfficiency;
1810 // step offset because we start with kCFStepAnaTopology
1811 step1 = (CFStep) ((Int_t) step1 - (Int_t) kCFStepAnaTopology);
1812 step2 = (CFStep) ((Int_t) step2 - (Int_t) kCFStepAnaTopology);
1817 // reset all limits and set the right ones except those in axis1, axis2 and axis3
1818 ResetBinLimits(sourceContainer->GetGrid(step1));
1819 ResetBinLimits(sourceContainer->GetGrid(step2));
1820 if (fEtaMax > fEtaMin && axis1 != 0 && axis2 != 0 && axis3 != 0)
1822 Printf("Restricted eta-range to %f %f", fEtaMin, fEtaMax);
1823 sourceContainer->GetGrid(step1)->SetRangeUser(0, fEtaMin, fEtaMax);
1824 sourceContainer->GetGrid(step2)->SetRangeUser(0, fEtaMin, fEtaMax);
1826 if (fPtMax > fPtMin && axis1 != 1 && axis2 != 1 && axis3 != 1)
1828 Printf("Restricted pt-range to %f %f", fPtMin, fPtMax);
1829 sourceContainer->GetGrid(step1)->SetRangeUser(1, fPtMin, fPtMax);
1830 sourceContainer->GetGrid(step2)->SetRangeUser(1, fPtMin, fPtMax);
1832 if (fCentralityMax > fCentralityMin && axis1 != 3 && axis2 != 3 && axis3 != 3)
1834 Printf("Restricted centrality range to %f %f", fCentralityMin, fCentralityMax);
1835 sourceContainer->GetGrid(step1)->SetRangeUser(3, fCentralityMin, fCentralityMax);
1836 sourceContainer->GetGrid(step2)->SetRangeUser(3, fCentralityMin, fCentralityMax);
1838 if (fZVtxMax > fZVtxMin && axis1 != 4 && axis2 != 4 && axis3 != 4)
1840 Printf("Restricted z-vtx range to %f %f", fZVtxMin, fZVtxMax);
1841 sourceContainer->GetGrid(step1)->SetRangeUser(4, fZVtxMin, fZVtxMax);
1842 sourceContainer->GetGrid(step2)->SetRangeUser(4, fZVtxMin, fZVtxMax);
1850 generated = sourceContainer->Project(step1, axis1, axis2, axis3);
1851 measured = sourceContainer->Project(step2, axis1, axis2, axis3);
1853 else if (axis2 >= 0)
1855 generated = sourceContainer->Project(step1, axis1, axis2);
1856 measured = sourceContainer->Project(step2, axis1, axis2);
1860 generated = sourceContainer->Project(step1, axis1);
1861 measured = sourceContainer->Project(step2, axis1);
1864 // check for bins with less than 50 entries, print warning
1872 binEnd[0] = generated->GetNbinsX();
1873 binEnd[1] = generated->GetNbinsY();
1874 binEnd[2] = generated->GetNbinsZ();
1876 if (fEtaMax > fEtaMin)
1880 binBegin[0] = generated->GetXaxis()->FindBin(fEtaMin);
1881 binEnd[0] = generated->GetXaxis()->FindBin(fEtaMax);
1885 binBegin[1] = generated->GetYaxis()->FindBin(fEtaMin);
1886 binEnd[1] = generated->GetYaxis()->FindBin(fEtaMax);
1890 binBegin[2] = generated->GetZaxis()->FindBin(fEtaMin);
1891 binEnd[2] = generated->GetZaxis()->FindBin(fEtaMax);
1895 if (fPtMax > fPtMin)
1897 // TODO this is just checking up to 15 for now
1898 Float_t ptMax = TMath::Min((Float_t) 15., fPtMax);
1901 binBegin[0] = generated->GetXaxis()->FindBin(fPtMin);
1902 binEnd[0] = generated->GetXaxis()->FindBin(ptMax);
1906 binBegin[1] = generated->GetYaxis()->FindBin(fPtMin);
1907 binEnd[1] = generated->GetYaxis()->FindBin(ptMax);
1911 binBegin[2] = generated->GetZaxis()->FindBin(fPtMin);
1912 binEnd[2] = generated->GetZaxis()->FindBin(ptMax);
1920 for (Int_t i=0; i<3; i++)
1921 vars[i] = binBegin[i];
1923 const Int_t limit = 50;
1926 if (generated->GetDimension() == 1 && generated->GetBinContent(vars[0]) < limit)
1928 Printf("Empty bin at %s=%.2f (%.2f entries)", generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]), generated->GetBinContent(vars[0]));
1931 else if (generated->GetDimension() == 2 && generated->GetBinContent(vars[0], vars[1]) < limit)
1933 Printf("Empty bin at %s=%.2f %s=%.2f (%.2f entries)",
1934 generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
1935 generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
1936 generated->GetBinContent(vars[0], vars[1]));
1939 else if (generated->GetDimension() == 3 && generated->GetBinContent(vars[0], vars[1], vars[2]) < limit)
1941 Printf("Empty bin at %s=%.2f %s=%.2f %s=%.2f (%.2f entries)",
1942 generated->GetXaxis()->GetTitle(), generated->GetXaxis()->GetBinCenter(vars[0]),
1943 generated->GetYaxis()->GetTitle(), generated->GetYaxis()->GetBinCenter(vars[1]),
1944 generated->GetZaxis()->GetTitle(), generated->GetZaxis()->GetBinCenter(vars[2]),
1945 generated->GetBinContent(vars[0], vars[1], vars[2]));
1950 if (vars[2] == binEnd[2]+1)
1952 vars[2] = binBegin[2];
1956 if (vars[1] == binEnd[1]+1)
1958 vars[1] = binBegin[1];
1962 if (vars[0] == binEnd[0]+1)
1967 Printf("Correction has %d empty bins (out of %d bins)", count, total);
1969 // rebin if required
1972 TAxis* axis = fEventHist->GetGrid(0)->GetGrid()->GetAxis(0);
1974 if (axis->GetNbins() < measured->GetNbinsX())
1978 // 2d rebin with variable axis does not exist in root
1980 TH1* tmp = measured;
1981 measured = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
1982 for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
1983 for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
1985 ((TH2*) measured)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
1986 measured->SetBinError(x, y, 0); // cannot trust bin error, set to 0
1991 generated = new TH2D(Form("%s_rebinned", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetXbins()->GetArray(), tmp->GetNbinsY(), tmp->GetYaxis()->GetXbins()->GetArray());
1992 for (Int_t x=1; x<=tmp->GetNbinsX(); x++)
1993 for (Int_t y=1; y<=tmp->GetNbinsY(); y++)
1995 ((TH2*) generated)->Fill(tmp->GetXaxis()->GetBinCenter(x), tmp->GetYaxis()->GetBinCenter(y), tmp->GetBinContent(x, y));
1996 generated->SetBinError(x, y, 0); // cannot trust bin error, set to 0
2002 TH1* tmp = measured;
2003 measured = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
2007 generated = tmp->Rebin(axis->GetNbins(), Form("%s_rebinned", tmp->GetName()), axis->GetXbins()->GetArray());
2011 else if (axis->GetNbins() > measured->GetNbinsX())
2014 AliFatal("Rebinning only works for 1d at present");
2016 // this is an unfortunate case where the number of bins has to be increased in principle
2017 // there is a region where the binning is finner in one histogram and a region where it is the other way round
2018 // this cannot be resolved in principle, but as we only calculate the ratio the bin in the second region get the same entries
2019 // only a certain binning is supported here
2020 if (axis->GetNbins() != 100 || measured->GetNbinsX() != 39)
2021 AliFatal(Form("Invalid binning --> %d %d", axis->GetNbins(), measured->GetNbinsX()));
2023 Double_t newBins[] = {0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 100.0};
2025 // reduce binning below 5 GeV/c
2026 TH1* tmp = measured;
2027 measured = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
2030 // increase binning above 5 GeV/c
2032 measured = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
2033 for (Int_t bin=1; bin<=measured->GetNbinsX(); bin++)
2035 measured->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(measured->GetBinCenter(bin))));
2036 measured->SetBinError(bin, tmp->GetBinError(tmp->FindBin(measured->GetBinCenter(bin))));
2040 // reduce binning below 5 GeV/c
2042 generated = tmp->Rebin(27, Form("%s_rebinned", tmp->GetName()), newBins);
2045 // increase binning above 5 GeV/c
2047 generated = new TH1F(Form("%s_rebinned2", tmp->GetName()), tmp->GetTitle(), axis->GetNbins(), axis->GetBinLowEdge(1), axis->GetBinUpEdge(axis->GetNbins()));
2048 for (Int_t bin=1; bin<=generated->GetNbinsX(); bin++)
2050 generated->SetBinContent(bin, tmp->GetBinContent(tmp->FindBin(generated->GetBinCenter(bin))));
2051 generated->SetBinError(bin, tmp->GetBinError(tmp->FindBin(generated->GetBinCenter(bin))));
2057 measured->Divide(measured, generated, 1, 1, "B");
2061 ResetBinLimits(sourceContainer->GetGrid(step1));
2062 ResetBinLimits(sourceContainer->GetGrid(step2));
2067 //____________________________________________________________________
2068 TH1* AliUEHist::GetEventEfficiency(CFStep step1, CFStep step2, Int_t axis1, Int_t axis2, Float_t ptLeadMin, Float_t ptLeadMax)
2070 // creates a event-level efficiency by dividing step2 by step1
2071 // projected to axis1 and axis2 (optional if >= 0)
2073 if (ptLeadMax > ptLeadMin)
2075 fEventHist->GetGrid(step1)->SetRangeUser(0, ptLeadMin, ptLeadMax);
2076 fEventHist->GetGrid(step2)->SetRangeUser(0, ptLeadMin, ptLeadMax);
2084 generated = fEventHist->Project(step1, axis1, axis2);
2085 measured = fEventHist->Project(step2, axis1, axis2);
2089 generated = fEventHist->Project(step1, axis1);
2090 measured = fEventHist->Project(step2, axis1);
2093 measured->Divide(measured, generated, 1, 1, "B");
2097 if (ptLeadMax > ptLeadMin)
2099 fEventHist->GetGrid(step1)->SetRangeUser(0, 0, -1);
2100 fEventHist->GetGrid(step2)->SetRangeUser(0, 0, -1);
2106 //____________________________________________________________________
2107 void AliUEHist::WeightHistogram(TH3* hist1, TH1* hist2)
2109 // weights each entry of the 3d histogram hist1 with the 1d histogram hist2
2110 // where the matching is done of the z axis of hist1 with the x axis of hist2
2112 if (hist1->GetNbinsZ() != hist2->GetNbinsX())
2113 AliFatal(Form("Inconsistent binning %d %d", hist1->GetNbinsZ(), hist2->GetNbinsX()));
2115 for (Int_t x=1; x<=hist1->GetNbinsX(); x++)
2117 for (Int_t y=1; y<=hist1->GetNbinsY(); y++)
2119 for (Int_t z=1; z<=hist1->GetNbinsZ(); z++)
2121 if (hist2->GetBinContent(z) > 0)
2123 hist1->SetBinContent(x, y, z, hist1->GetBinContent(x, y, z) / hist2->GetBinContent(z));
2124 hist1->SetBinError(x, y, z, hist1->GetBinError(x, y, z) / hist2->GetBinContent(z));
2128 hist1->SetBinContent(x, y, z, 0);
2129 hist1->SetBinError(x, y, z, 0);
2136 //____________________________________________________________________
2137 TH1* AliUEHist::GetBias(CFStep step1, CFStep step2, Int_t region, const char* axis, Float_t leadPtMin, Float_t leadPtMax, Int_t weighting)
2139 // extracts the track-level bias (integrating out the multiplicity) between two steps (dividing step2 by step1)
2140 // in the given region (sum over all regions is calculated if region == -1)
2141 // done by weighting the track-level distribution with the number of events as function of leading pT
2142 // and then calculating the ratio between the distributions
2143 // projected to axis which is a TH3::Project3D string, e.g. "x", or "yx"
2144 // no projection is done if axis == 0
2145 // weighting: 0 = tracks weighted with events (as discussed above)
2146 // 1 = only track bias is returned
2147 // 2 = only event bias is returned
2149 AliCFContainer* tmp = 0;
2153 tmp = (AliCFContainer*) fTrackHist[0]->Clone();
2154 for (UInt_t i = 1; i < fkRegions; i++)
2156 tmp->Add(fTrackHist[i]);
2158 else if (region == kMin && fCombineMinMax)
2160 tmp = (AliCFContainer*) fTrackHist[kMin]->Clone();
2161 tmp->Add(fTrackHist[kMax]);
2164 tmp = fTrackHist[region];
2166 ResetBinLimits(tmp->GetGrid(step1));
2167 ResetBinLimits(fEventHist->GetGrid(step1));
2168 SetBinLimits(tmp->GetGrid(step1));
2170 ResetBinLimits(tmp->GetGrid(step2));
2171 ResetBinLimits(fEventHist->GetGrid(step2));
2172 SetBinLimits(tmp->GetGrid(step2));
2174 TH1D* events1 = (TH1D*)fEventHist->Project(step1, 0);
2175 TH3D* hist1 = (TH3D*)tmp->Project(step1, 0, tmp->GetNVar()-1, 2);
2177 WeightHistogram(hist1, events1);
2179 TH1D* events2 = (TH1D*)fEventHist->Project(step2, 0);
2180 TH3D* hist2 = (TH3D*)tmp->Project(step2, 0, tmp->GetNVar()-1, 2);
2182 WeightHistogram(hist2, events2);
2184 TH1* generated = hist1;
2185 TH1* measured = hist2;
2187 if (weighting == 0 || weighting == 1)
2191 if (leadPtMax > leadPtMin)
2193 hist1->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
2194 hist2->GetZaxis()->SetRangeUser(leadPtMin, leadPtMax);
2197 if (fEtaMax > fEtaMin && !TString(axis).Contains("x"))
2199 hist1->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
2200 hist2->GetXaxis()->SetRangeUser(fEtaMin, fEtaMax);
2203 generated = hist1->Project3D(axis);
2204 measured = hist2->Project3D(axis);
2206 // delete hists here if projection has been done
2213 else if (weighting == 2)
2217 generated = events1;
2221 measured->Divide(generated);
2225 ResetBinLimits(tmp->GetGrid(step1));
2226 ResetBinLimits(tmp->GetGrid(step2));
2228 if ((region == -1) || (region == kMin && fCombineMinMax))
2234 //____________________________________________________________________
2235 void AliUEHist::CorrectCorrelatedContamination(CFStep step, Int_t region, TH1* trackCorrection)
2237 // corrects for the given factor in a small delta-eta and delta-phi window as function of pT,A and pT,T
2239 if (!fTrackHist[region])
2242 THnSparse* grid = fTrackHist[region]->GetGrid(step)->GetGrid();
2247 if (grid->GetAxis(var1)->GetNbins() != trackCorrection->GetNbinsX())
2248 AliFatal(Form("Invalid binning (var1): %d %d", grid->GetAxis(var1)->GetNbins(), trackCorrection->GetNbinsX()));
2250 if (grid->GetAxis(var2)->GetNbins() != trackCorrection->GetNbinsY())
2251 AliFatal(Form("Invalid binning (var2): %d %d", grid->GetAxis(var2)->GetNbins(), trackCorrection->GetNbinsY()));
2253 // optimized implementation
2254 for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
2258 Double_t value = grid->GetBinContent(binIdx, bins);
2259 Double_t error = grid->GetBinError(binIdx);
2261 // check eta and phi axes
2262 if (TMath::Abs(grid->GetAxis(0)->GetBinCenter(bins[0])) > 0.1)
2264 if (TMath::Abs(grid->GetAxis(4)->GetBinCenter(bins[4])) > 0.1)
2267 value *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
2268 error *= trackCorrection->GetBinContent(bins[var1], bins[var2]);
2270 grid->SetBinContent(bins, value);
2271 grid->SetBinError(bins, error);
2274 Printf("AliUEHist::CorrectCorrelatedContamination: Corrected.");
2277 //____________________________________________________________________
2278 TH2* AliUEHist::GetCorrelatedContamination()
2280 // contamination correlated with the trigger particle is evaluated between step kCFStepTracked and kCFStepTrackedOnlyPrim in the region of delta eta and delta phi < 0.1 (smallest bin!)
2282 Int_t step1 = kCFStepTrackedOnlyPrim;
2283 Int_t step2 = kCFStepTracked;
2285 fTrackHist[0]->GetGrid(step1)->SetRangeUser(0, -0.01, 0.01); // delta eta
2286 fTrackHist[0]->GetGrid(step1)->SetRangeUser(4, -0.01, 0.01); // delta phi
2287 TH2* tracksStep1 = (TH2*) fTrackHist[0]->Project(step1, 1, 2);
2289 fTrackHist[0]->GetGrid(step2)->SetRangeUser(0, -0.01, 0.01); // delta eta
2290 fTrackHist[0]->GetGrid(step2)->SetRangeUser(4, -0.01, 0.01); // delta phi
2291 TH2* tracksStep2 = (TH2*) fTrackHist[0]->Project(step2, 1, 2);
2293 tracksStep1->Divide(tracksStep2);
2295 TH1* triggersStep1 = fEventHist->Project(step1, 0);
2296 TH1* triggersStep2 = fEventHist->Project(step2, 0);
2298 TH1* singleParticle = GetTrackingContamination(1);
2300 for (Int_t x=1; x<=tracksStep1->GetNbinsX(); x++)
2301 for (Int_t y=1; y<=tracksStep1->GetNbinsY(); y++)
2302 if (singleParticle->GetBinContent(x) > 0 && triggersStep1->GetBinContent(y) > 0)
2303 tracksStep1->SetBinContent(x, y, tracksStep1->GetBinContent(x, y) / triggersStep1->GetBinContent(y) * triggersStep2->GetBinContent(y) / singleParticle->GetBinContent(x));
2305 tracksStep1->SetBinContent(x, y, 0);
2307 delete singleParticle;
2309 delete triggersStep1;
2310 delete triggersStep2;
2315 //____________________________________________________________________
2316 TH2D* AliUEHist::GetTrackingEfficiency()
2318 // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
2319 // integrates over the regions and all other variables than pT and eta to increase the statistics
2321 // returned histogram has to be deleted by the user
2323 return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 0, 1));
2326 //____________________________________________________________________
2327 TH2D* AliUEHist::GetFakeRate()
2329 return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepReconstructed, 0, 1));
2332 //____________________________________________________________________
2333 TH2D* AliUEHist::GetTrackingEfficiencyCentrality()
2335 // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
2336 // integrates over the regions and all other variables than pT, centrality to increase the statistics
2338 // returned histogram has to be deleted by the user
2340 return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, 1, 3));
2343 //____________________________________________________________________
2344 TH1D* AliUEHist::GetTrackingEfficiency(Int_t axis)
2346 // extracts the tracking efficiency by calculating the efficiency from step kCFStepAnaTopology to kCFStepTrackedOnlyPrim
2347 // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2349 return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepAnaTopology, kCFStepTrackedOnlyPrim, axis));
2352 //____________________________________________________________________
2353 TH1D* AliUEHist::GetFakeRate(Int_t axis)
2355 return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepReconstructed, axis));
2357 //____________________________________________________________________
2358 TH2D* AliUEHist::GetTrackingCorrection()
2360 // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2361 // integrates over the regions and all other variables than pT and eta to increase the statistics
2363 // returned histogram has to be deleted by the user
2365 return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, 0, 1));
2368 //____________________________________________________________________
2369 TH1D* AliUEHist::GetTrackingCorrection(Int_t axis)
2371 // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2372 // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2374 return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepAnaTopology, axis));
2377 //____________________________________________________________________
2378 TH2D* AliUEHist::GetTrackingEfficiencyCorrection()
2380 // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2381 // integrates over the regions and all other variables than pT and eta to increase the statistics
2383 // returned histogram has to be deleted by the user
2385 return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 0, 1));
2388 //____________________________________________________________________
2389 TH2D* AliUEHist::GetTrackingEfficiencyCorrectionCentrality()
2391 // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2392 // integrates over the regions and all other variables than pT and centrality to increase the statistics
2394 // returned histogram has to be deleted by the user
2396 return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, 1, 3));
2399 //____________________________________________________________________
2400 TH1D* AliUEHist::GetTrackingEfficiencyCorrection(Int_t axis)
2402 // extracts the tracking correction by calculating the efficiency from step kCFStepAnaTopology to kCFStepTracked
2403 // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2405 return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTrackedOnlyPrim, kCFStepAnaTopology, axis));
2408 //____________________________________________________________________
2409 TH2D* AliUEHist::GetTrackingContamination()
2411 // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
2412 // integrates over the regions and all other variables than pT and eta to increase the statistics
2414 // returned histogram has to be deleted by the user
2416 return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 0, 1));
2419 //____________________________________________________________________
2420 TH2D* AliUEHist::GetTrackingContaminationCentrality()
2422 // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
2423 // integrates over the regions and all other variables than pT and centrality to increase the statistics
2425 // returned histogram has to be deleted by the user
2427 return dynamic_cast<TH2D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, 1, 3));
2430 //____________________________________________________________________
2431 TH1D* AliUEHist::GetTrackingContamination(Int_t axis)
2433 // extracts the tracking contamination by secondaries by calculating the efficiency from step kCFStepTrackedOnlyPrim to kCFStepTracked
2434 // integrates over the regions and all other variables than pT (axis == 0) and eta (axis == 1) to increase the statistics
2436 return dynamic_cast<TH1D*> (GetTrackEfficiency(kCFStepTracked, kCFStepTrackedOnlyPrim, axis));
2439 //____________________________________________________________________
2440 const char* AliUEHist::GetRegionTitle(Region region)
2442 // returns the name of the given region
2451 return (fCombineMinMax) ? "Transverse" : "Min";
2459 //____________________________________________________________________
2460 const char* AliUEHist::GetStepTitle(CFStep step)
2462 // returns the name of the given step
2467 return "All events";
2468 case kCFStepTriggered:
2471 return "Primary Vertex";
2472 case kCFStepAnaTopology:
2473 return "Required analysis topology";
2474 case kCFStepTrackedOnlyPrim:
2475 return "Tracked (matched MC, only primaries)";
2476 case kCFStepTracked:
2477 return "Tracked (matched MC, all)";
2478 case kCFStepReconstructed:
2479 return "Reconstructed";
2480 case kCFStepRealLeading:
2481 return "Correct leading particle identified";
2482 case kCFStepBiasStudy:
2483 return "Bias study applying tracking efficiency";
2484 case kCFStepBiasStudy2:
2485 return "Bias study applying tracking efficiency in two steps";
2486 case kCFStepCorrected:
2487 return "Corrected for efficiency on-the-fly";
2493 //____________________________________________________________________
2494 void AliUEHist::CopyReconstructedData(AliUEHist* from)
2496 // copies those histograms extracted from ESD to this object
2498 // TODO at present only the pointers are copied
2500 for (Int_t region=0; region<4; region++)
2502 if (!fTrackHist[region])
2505 fTrackHist[region]->SetGrid(AliUEHist::kCFStepReconstructed, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepReconstructed));
2506 //fTrackHist[region]->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
2507 fTrackHist[region]->SetGrid(AliUEHist::kCFStepBiasStudy, from->fTrackHist[region]->GetGrid(AliUEHist::kCFStepBiasStudy));
2510 fEventHist->SetGrid(AliUEHist::kCFStepReconstructed, from->fEventHist->GetGrid(AliUEHist::kCFStepReconstructed));
2511 //fEventHist->SetGrid(AliUEHist::kCFStepTrackedOnlyPrim, from->fEventHist->GetGrid(AliUEHist::kCFStepTrackedOnlyPrim));
2512 fEventHist->SetGrid(AliUEHist::kCFStepBiasStudy, from->fEventHist->GetGrid(AliUEHist::kCFStepBiasStudy));
2515 void AliUEHist::DeepCopy(AliUEHist* from)
2517 // copies the entries of this object's members from the object <from> to this object
2518 // fills using the fill function and thus allows that the objects have different binning
2520 for (Int_t region=0; region<4; region++)
2522 if (!fTrackHist[region] || !from->fTrackHist[region])
2525 for (Int_t step=0; step<fTrackHist[region]->GetNStep(); step++)
2527 Printf("Copying region %d step %d", region, step);
2528 THnSparse* target = fTrackHist[region]->GetGrid(step)->GetGrid();
2529 THnSparse* source = from->fTrackHist[region]->GetGrid(step)->GetGrid();
2532 target->RebinnedAdd(source);
2536 for (Int_t step=0; step<fEventHist->GetNStep(); step++)
2538 Printf("Ev: Copying step %d", step);
2539 THnSparse* target = fEventHist->GetGrid(step)->GetGrid();
2540 THnSparse* source = from->fEventHist->GetGrid(step)->GetGrid();
2543 target->RebinnedAdd(source);
2546 for (Int_t step=0; step<TMath::Min(fTrackHistEfficiency->GetNStep(), from->fTrackHistEfficiency->GetNStep()); step++)
2548 if (!from->fTrackHistEfficiency->GetGrid(step))
2551 Printf("Eff: Copying step %d", step);
2552 THnSparse* target = fTrackHistEfficiency->GetGrid(step)->GetGrid();
2553 THnSparse* source = from->fTrackHistEfficiency->GetGrid(step)->GetGrid();
2556 target->RebinnedAdd(source);
2560 //____________________________________________________________________
2561 void AliUEHist::ExtendTrackingEfficiency(Bool_t verbose)
2563 // fits the tracking efficiency at high pT with a constant and fills all bins with this tracking efficiency
2565 Float_t fitRangeBegin = 5.01;
2566 Float_t fitRangeEnd = 14.99;
2567 Float_t extendRangeBegin = 10.01;
2569 if (fTrackHistEfficiency->GetNVar() == 3)
2571 TH1* obj = GetTrackingEfficiency(1);
2579 obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
2581 Float_t trackingEff = obj->GetFunction("pol0")->GetParameter(0);
2583 obj = GetTrackingContamination(1);
2591 obj->Fit("pol0", (verbose) ? "+" : "0+", "SAME", fitRangeBegin, fitRangeEnd);
2593 Float_t trackingCont = obj->GetFunction("pol0")->GetParameter(0);
2595 Printf("AliUEHist::ExtendTrackingEfficiency: Fitted efficiency between %f and %f and got %f tracking efficiency and %f tracking contamination correction. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, trackingEff, trackingCont, extendRangeBegin, fEtaMin, fEtaMax);
2597 // extend for full pT range
2598 for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
2599 for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
2600 for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
2608 fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
2609 fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff);
2610 fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff / trackingCont);
2613 else if (fTrackHistEfficiency->GetNVar() == 4)
2615 // fit in centrality intervals of 20% for efficiency, one bin for contamination
2616 Float_t* trackingEff = 0;
2617 Float_t* trackingCont = 0;
2618 Float_t centralityBins[] = { 0, 10, 20, 40, 60, 100 };
2619 Int_t nCentralityBins = 5;
2621 Printf("AliUEHist::ExtendTrackingEfficiency: Fitting efficiencies between %f and %f. Extending from %f onwards (within %f < eta < %f)", fitRangeBegin, fitRangeEnd, extendRangeBegin, fEtaMin, fEtaMax);
2623 // 0 = eff; 1 = cont
2624 for (Int_t caseNo = 0; caseNo < 2; caseNo++)
2626 Float_t* target = 0;
2627 Int_t centralityBinsLocal = nCentralityBins;
2631 trackingEff = new Float_t[centralityBinsLocal];
2632 target = trackingEff;
2636 centralityBinsLocal = 1;
2637 trackingCont = new Float_t[centralityBinsLocal];
2638 target = trackingCont;
2641 for (Int_t i=0; i<centralityBinsLocal; i++)
2643 if (centralityBinsLocal == 1)
2644 SetCentralityRange(centralityBins[0] + 0.1, centralityBins[nCentralityBins] - 0.1);
2646 SetCentralityRange(centralityBins[i] + 0.1, centralityBins[i+1] - 0.1);
2647 TH1* proj = (caseNo == 0) ? GetTrackingEfficiency(1) : GetTrackingContamination(1);
2653 if ((Int_t) proj->Fit("pol0", (verbose) ? "+" : "Q0+", "SAME", fitRangeBegin, fitRangeEnd) == 0)
2654 target[i] = proj->GetFunction("pol0")->GetParameter(0);
2657 Printf("AliUEHist::ExtendTrackingEfficiency: case %d, centrality %d, eff %f", caseNo, i, target[i]);
2661 // extend for full pT range
2662 for (Int_t x = fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMin); x <= fTrackHistEfficiency->GetAxis(0, 0)->FindBin(fEtaMax); x++)
2663 for (Int_t y = fTrackHistEfficiency->GetAxis(1, 0)->FindBin(extendRangeBegin); y <= fTrackHistEfficiency->GetNBins(1); y++)
2664 for (Int_t z = 1; z <= fTrackHistEfficiency->GetNBins(2); z++) // particle type axis
2666 for (Int_t z2 = 1; z2 <= fTrackHistEfficiency->GetNBins(3); z2++) // centrality axis
2676 while (centralityBins[z2Bin+1] < fTrackHistEfficiency->GetAxis(3, 0)->GetBinCenter(z2))
2679 //Printf("%d %d", z2, z2Bin);
2681 fTrackHistEfficiency->GetGrid(0)->SetElement(bins, 100);
2682 fTrackHistEfficiency->GetGrid(1)->SetElement(bins, 100.0 * trackingEff[z2Bin]);
2683 if (trackingCont[0] > 0)
2684 fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 100.0 * trackingEff[z2Bin] / trackingCont[0]);
2686 fTrackHistEfficiency->GetGrid(2)->SetElement(bins, 0);
2690 delete[] trackingEff;
2691 delete[] trackingCont;
2694 SetCentralityRange(0, 0);
2697 void AliUEHist::AdditionalDPhiCorrection(Int_t step)
2699 // corrects the dphi distribution with an extra factor close to dphi ~ 0
2701 Printf("WARNING: In AliUEHist::AdditionalDPhiCorrection.");
2703 THnSparse* grid = fTrackHist[0]->GetGrid(step)->GetGrid();
2705 // optimized implementation
2706 for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
2709 Double_t value = grid->GetBinContent(binIdx, bins);
2710 Double_t error = grid->GetBinError(binIdx);
2712 Float_t binCenter = grid->GetAxis(4)->GetBinCenter(bins[4]);
2713 if (TMath::Abs(binCenter) < 0.2)
2718 else if (TMath::Abs(binCenter) < 0.3)
2724 grid->SetBinContent(bins, value);
2725 grid->SetBinError(bins, error);
2729 void AliUEHist::Scale(Double_t factor)
2731 // scales all contained histograms by the given factor
2733 for (Int_t i=0; i<4; i++)
2735 fTrackHist[i]->Scale(factor);
2737 fEventHist->Scale(factor);
2738 fTrackHistEfficiency->Scale(factor);
2741 void AliUEHist::Reset()
2743 // resets all contained histograms
2745 for (Int_t i=0; i<4; i++)
2747 for (Int_t step=0; step<fTrackHist[i]->GetNStep(); step++)
2748 fTrackHist[i]->GetGrid(step)->GetGrid()->Reset();
2750 for (Int_t step=0; step<fEventHist->GetNStep(); step++)
2751 fEventHist->GetGrid(step)->GetGrid()->Reset();
2753 for (Int_t step=0; step<fTrackHistEfficiency->GetNStep(); step++)
2754 fTrackHistEfficiency->GetGrid(step)->GetGrid()->Reset();
2757 THnBase* AliUEHist::ChangeToThn(THnBase* sparse)
2759 // change the object to THn for faster processing
2761 // convert to THn (SEGV's for some strange reason...)
2762 // x = THn::CreateHn("a", "a", sparse);
2764 // own implementation
2766 for (Int_t i=0; i<sparse->GetNdimensions(); i++)
2767 nBins[i] = sparse->GetAxis(i)->GetNbins();
2768 THn* tmpTHn = new THnF(Form("%s_thn", sparse->GetName()), sparse->GetTitle(), sparse->GetNdimensions(), nBins, 0, 0);
2769 for (Int_t i=0; i<sparse->GetNdimensions(); i++)
2771 tmpTHn->SetBinEdges(i, sparse->GetAxis(i)->GetXbins()->GetArray());
2772 tmpTHn->GetAxis(i)->SetTitle(sparse->GetAxis(i)->GetTitle());
2774 tmpTHn->RebinnedAdd(sparse);
2779 void AliUEHist::CondenseBin(THnSparse* grid, THnSparse* target, Int_t axis, Float_t targetValue, Float_t from, Float_t to)
2782 // loops through the histogram and moves all entries to a single point <targetValue> on the axis <axis>
2783 // if <from> and <to> are set, then moving only occurs if value on <axis> is betweem <from> and <to>
2786 if (grid->GetNdimensions() > 6)
2787 AliFatal("Too many dimensions in THnSparse");
2789 Int_t targetBin = grid->GetAxis(axis)->FindBin(targetValue);
2790 AliInfo(Form("Target bin on axis %d with value %f is %d", axis, targetValue, targetBin));
2793 Int_t toBin = grid->GetAxis(axis)->GetNbins();
2796 fromBin = grid->GetAxis(axis)->FindBin(from);
2797 toBin = grid->GetAxis(axis)->FindBin(to);
2798 AliInfo(Form("Only condensing between bin %d and %d", fromBin, toBin));
2802 for (Int_t binIdx = 0; binIdx < grid->GetNbins(); binIdx++)
2804 Double_t value = grid->GetBinContent(binIdx, bins);
2805 Double_t error = grid->GetBinError(binIdx);
2807 if (bins[axis] >= fromBin && bins[axis] <= toBin)
2808 bins[axis] = targetBin;
2810 value += target->GetBinContent(bins);
2811 error = TMath::Sqrt(error * error + target->GetBinError(bins) * target->GetBinError(bins));
2813 target->SetBinContent(bins, value);
2814 target->SetBinError(bins, error);
2818 void AliUEHist::CondenseBin(CFStep step, Int_t trackAxis, Int_t eventAxis, Float_t targetValue, Float_t from, Float_t to, CFStep tmpStep)
2820 // loops through the histogram at <step> and moves all entries to a single point <targetValue> on the axes
2821 // <trackAxis> and <eventAxis>. <tmpStep> is used to temporary store the data
2822 // This is useful e.g. to move bin content around for MC productions where the centrality selection did
2823 // not yield the desired result
2826 fEventHist->GetGrid(tmpStep)->GetGrid()->Reset();
2827 for (UInt_t i=0; i<fkRegions; i++)
2829 fTrackHist[i]->GetGrid(tmpStep)->GetGrid()->Reset();
2832 CorrectTracks(step, tmpStep, 0, -1);
2833 CorrectEvents(step, tmpStep, 0, -1);
2836 fEventHist->GetGrid(step)->GetGrid()->Reset();
2837 for (UInt_t i=0; i<fkRegions; i++)
2839 fTrackHist[i]->GetGrid(step)->GetGrid()->Reset();
2842 for (UInt_t i=0; i<fkRegions; i++)
2847 THnSparse* grid = fTrackHist[i]->GetGrid(tmpStep)->GetGrid();
2848 THnSparse* target = fTrackHist[i]->GetGrid(step)->GetGrid();
2850 CondenseBin(grid, target, trackAxis, targetValue, from, to);
2852 CondenseBin(fEventHist->GetGrid(tmpStep)->GetGrid(), fEventHist->GetGrid(step)->GetGrid(), eventAxis, targetValue, from, to);
2855 fEventHist->GetGrid(tmpStep)->GetGrid()->Reset();
2856 for (UInt_t i=0; i<fkRegions; i++)
2858 fTrackHist[i]->GetGrid(tmpStep)->GetGrid()->Reset();