2 // This class calculates the exclusive charged particle density
3 // in each for the 5 FMD rings.
5 #include "AliFMDCorrector.h"
10 #include "AliForwardCorrectionManager.h"
11 #include "AliFMDCorrSecondaryMap.h"
12 #include "AliFMDCorrVertexBias.h"
13 #include "AliFMDCorrMergingEfficiency.h"
14 #include "AliFMDCorrAcceptance.h"
22 ClassImp(AliFMDCorrector)
27 //____________________________________________________________________
28 AliFMDCorrector::AliFMDCorrector()
31 fUseSecondaryMap(true),
32 fUseVertexBias(false),
33 fUseAcceptance(false),
34 fUseMergingEfficiency(false),
38 DGUARD(fDebug, 3, "Default CTOR of AliFMDCorrector");
41 //____________________________________________________________________
42 AliFMDCorrector::AliFMDCorrector(const char* title)
43 : TNamed("fmdCorrector", title),
45 fUseSecondaryMap(true),
46 fUseVertexBias(false),
47 fUseAcceptance(false),
48 fUseMergingEfficiency(false),
55 DGUARD(fDebug, 3, "Named CTOR of AliFMDCorrector: %s", title);
56 fRingHistos.SetName(GetName());
57 fRingHistos.Add(new RingHistos(1, 'I'));
58 fRingHistos.Add(new RingHistos(2, 'I'));
59 fRingHistos.Add(new RingHistos(2, 'O'));
60 fRingHistos.Add(new RingHistos(3, 'I'));
61 fRingHistos.Add(new RingHistos(3, 'O'));
64 //____________________________________________________________________
65 AliFMDCorrector::AliFMDCorrector(const AliFMDCorrector& o)
68 fUseSecondaryMap(o.fUseSecondaryMap),
69 fUseVertexBias(o.fUseVertexBias),
70 fUseAcceptance(o.fUseAcceptance),
71 fUseMergingEfficiency(o.fUseMergingEfficiency),
77 // o Object to copy from
78 DGUARD(fDebug, 3, "Copy CTOR of AliFMDCorrector");
79 TIter next(&o.fRingHistos);
81 while ((obj = next())) fRingHistos.Add(obj);
84 //____________________________________________________________________
85 AliFMDCorrector::~AliFMDCorrector()
90 DGUARD(fDebug, 3, "DTOR of AliFMDCorrector");
94 //____________________________________________________________________
96 AliFMDCorrector::operator=(const AliFMDCorrector& o)
98 // Assignment operator
101 // o Object to assign from
102 DGUARD(fDebug, 3, "Assignment of AliFMDCorrector");
103 if (&o == this) return *this;
104 TNamed::operator=(o);
107 fRingHistos.Delete();
108 fUseSecondaryMap = o.fUseSecondaryMap;
109 fUseVertexBias = o.fUseVertexBias;
110 fUseAcceptance = o.fUseAcceptance;
111 fUseMergingEfficiency = o.fUseMergingEfficiency;
112 TIter next(&o.fRingHistos);
114 while ((obj = next())) fRingHistos.Add(obj);
119 //____________________________________________________________________
121 AliFMDCorrector::SetupForData(const TAxis&)
124 // Initialize this object
127 // etaAxis Eta axis to use
129 DGUARD(fDebug, 1, "Initialization of AliFMDCorrector");
130 if (!fUseSecondaryMap)
131 AliWarning("Secondary maps not used - BE CAREFUL");
133 AliWarning("Vertex bias not used");
135 AliWarning("Acceptance from dead-channels not used");
138 //____________________________________________________________________
139 AliFMDCorrector::RingHistos*
140 AliFMDCorrector::GetRingHistos(UShort_t d, Char_t r) const
143 // Get the ring histogram container
149 // Ring histogram container
153 case 1: idx = 0; break;
154 case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
155 case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
157 if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
159 return static_cast<RingHistos*>(fRingHistos.At(idx));
162 //____________________________________________________________________
164 AliFMDCorrector::DivideMap(TH2* num, const TH2* denom,
165 Bool_t alsoUnderOver) const
168 // Implement TH1::Divide but
169 // - Assume compatible histograms
170 // - Unless third argument is true, do not divide over/under flow bins
172 if (!num || !denom) return;
174 Int_t first = (alsoUnderOver ? 0 : 1);
175 Int_t lastX = num->GetNbinsX() + (alsoUnderOver ? 1 : 0);
176 Int_t lastY = num->GetNbinsY() + (alsoUnderOver ? 1 : 0);
178 for (Int_t ix = first; ix <= lastX; ix++) {
179 for (Int_t iy = first; iy <= lastY; iy++) {
180 Int_t bin = num->GetBin(ix,iy);
181 Double_t c0 = num->GetBinContent(bin);
182 Double_t c1 = denom->GetBinContent(bin);
184 num->SetBinContent(bin,0);
185 num->SetBinError(bin, 0);
188 Double_t w = c0 / c1;
189 Double_t e0 = num->GetBinError(bin);
190 Double_t e1 = denom->GetBinError(bin);
191 Double_t c12 = c1*c1;
192 Double_t e2 = (e0*e0*c1*c1 + e1*e1*c0*c0)/(c12*c12);
194 num->SetBinContent(bin, w);
195 num->SetBinError(bin, TMath::Sqrt(e2));
199 //____________________________________________________________________
201 AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
205 // Do the calculations
207 // hists Cache of histograms
213 DGUARD(fDebug, 1, "Correct histograms of AliFMDCorrector");
214 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
216 UShort_t uvb = vtxbin;
217 for (UShort_t d=1; d<=3; d++) {
218 UShort_t nr = (d == 1 ? 1 : 2);
219 for (UShort_t q=0; q<nr; q++) {
220 Char_t r = (q == 0 ? 'I' : 'O');
221 TH2D* h = hists.Get(d,r);
222 RingHistos* rh = GetRingHistos(d,r);
224 if (fUseSecondaryMap) {
225 TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d, r, uvb);
227 AliWarning(Form("No secondary correction for FMDM%d%c "
228 "in vertex bin %d", d, r, uvb));
231 // Divide by primary/total ratio
232 DivideMap(h, bg, false);
234 if (fUseVertexBias) {
235 TH2D* ef = fcm.GetVertexBias()->GetCorrection(r, uvb);
237 AliWarning(Form("No event %s vertex bias correction in vertex bin %d",
238 (r == 'I' || r == 'i' ? "inner" : "outer"), uvb));
241 // Divide by the event selection efficiency
242 DivideMap(h, ef, false);
244 if (fUseAcceptance) {
245 TH2D* ac = fcm.GetAcceptance()->GetCorrection(d, r, uvb);
247 AliWarning(Form("No acceptance correction for FMD%d%c in "
248 "vertex bin %d", d, r, uvb));
251 // Fill overflow bin with ones
252 for (Int_t i = 1; i <= h->GetNbinsX(); i++)
253 h->SetBinContent(i, h->GetNbinsY()+1, 1);
255 // Divide by the acceptance correction -
256 DivideMap(h, ac, fcm.GetAcceptance()->HasOverflow());
259 if (fUseMergingEfficiency) {
260 if (!fcm.GetMergingEfficiency()) {
261 AliWarning("No merging efficiencies");
264 TH1D* sf = fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb);
265 if (!fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb)) {
266 AliWarning(Form("No merging efficiency for FMD%d%c at vertex bin %d",
272 for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
273 Float_t c = sf->GetBinContent(ieta);
274 Float_t ec = sf->GetBinError(ieta);
276 for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) {
278 h->SetBinContent(ieta,iphi,0);
279 h->SetBinError(ieta,iphi,0);
283 Double_t m = h->GetBinContent(ieta, iphi) / c;
284 Double_t em = h->GetBinError(ieta, iphi);
286 Double_t e = TMath::Sqrt(em * em + (m * ec) * (m * ec));
288 h->SetBinContent(ieta,iphi,m);
289 h->SetBinError(ieta,iphi,e);
294 rh->fDensity->Add(h);
301 //____________________________________________________________________
303 AliFMDCorrector::Terminate(const TList* dir, TList* output, Int_t nEvents)
306 // Scale the histograms to the total number of events
308 // dir Where the output is stored
309 // nEvents Number of events
311 DGUARD(fDebug, 1, "Scale histograms of AliFMDCorrector");
312 if (nEvents <= 0) return;
313 TList* d = static_cast<TList*>(dir->FindObject(GetName()));
316 TList* out = new TList;
317 out->SetName(d->GetName());
320 TIter next(&fRingHistos);
322 THStack* sums = new THStack("sums", "Sums of ring results");
323 while ((o = static_cast<RingHistos*>(next()))) {
324 o->Terminate(d, nEvents);
326 Warning("Terminate", "No density from %s", o->GetName());
329 TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
330 o->fDensity->GetNbinsY(),"e");
331 sum->Scale(1., "width");
332 sum->SetTitle(o->GetName());
333 sum->SetDirectory(0);
334 sum->SetYTitle("#sum N_{ch,primary}");
340 //____________________________________________________________________
342 AliFMDCorrector::CreateOutputObjects(TList* dir)
345 // Output diagnostic histograms to directory
348 // dir List to write in
350 DGUARD(fDebug, 1, "Define output of AliFMDCorrector");
351 TList* d = new TList;
352 d->SetName(GetName());
355 d->Add(AliForwardUtil::MakeParameter("secondary", fUseSecondaryMap));
356 d->Add(AliForwardUtil::MakeParameter("vertexBias", fUseVertexBias));
357 d->Add(AliForwardUtil::MakeParameter("acceptance", fUseAcceptance));
358 d->Add(AliForwardUtil::MakeParameter("merging", fUseMergingEfficiency));
361 TIter next(&fRingHistos);
363 while ((o = static_cast<RingHistos*>(next()))) {
364 o->CreateOutputObjects(d);
368 #define PF(N,V,...) \
369 AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
370 #define PFB(N,FLAG) \
372 AliForwardUtil::PrintName(N); \
373 std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
375 #define PFV(N,VALUE) \
377 AliForwardUtil::PrintName(N); \
378 std::cout << (VALUE) << std::endl; } while(false)
380 //____________________________________________________________________
382 AliFMDCorrector::Print(Option_t* /* option */) const
389 AliForwardUtil::PrintTask(*this);
390 gROOT->IncreaseDirLevel();
391 PFB("Use secondary maps", fUseSecondaryMap );
392 PFB("Use vertex bias", fUseVertexBias );
393 PFB("Use acceptance", fUseAcceptance );
394 PFB("Use merging efficiency", fUseMergingEfficiency);
395 gROOT->DecreaseDirLevel();
398 //====================================================================
399 AliFMDCorrector::RingHistos::RingHistos()
400 : AliForwardUtil::RingHistos(),
408 //____________________________________________________________________
409 AliFMDCorrector::RingHistos::RingHistos(UShort_t d, Char_t r)
410 : AliForwardUtil::RingHistos(d,r),
419 fDensity = new TH2D("primaryDensity",
420 "#sum N_{ch,primary}/(#Delta#eta#Delta#phi)",
421 200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
423 fDensity->SetDirectory(0);
424 fDensity->SetMarkerColor(Color());
426 fDensity->SetXTitle("#eta");
427 fDensity->SetYTitle("#phi [radians]");
428 fDensity->SetZTitle("Primary N_{ch} density");
430 //____________________________________________________________________
431 AliFMDCorrector::RingHistos::RingHistos(const RingHistos& o)
432 : AliForwardUtil::RingHistos(o),
438 // o Object to copy from
442 //____________________________________________________________________
443 AliFMDCorrector::RingHistos&
444 AliFMDCorrector::RingHistos::operator=(const RingHistos& o)
447 // Assignment operator
449 // o Object to assign from
454 if (&o == this) return *this;
455 AliForwardUtil::RingHistos::operator=(o);
457 if (fDensity) delete fDensity;
459 fDensity = static_cast<TH2D*>(o.fDensity->Clone());
463 //____________________________________________________________________
464 AliFMDCorrector::RingHistos::~RingHistos()
469 // if (fDensity) delete fDensity;
472 //____________________________________________________________________
474 AliFMDCorrector::RingHistos::CreateOutputObjects(TList* dir)
479 // dir Where to put it
481 TList* d = DefineOutputList(dir);
485 //____________________________________________________________________
487 AliFMDCorrector::RingHistos::Terminate(TList* dir, Int_t nEvents)
490 // Scale the histograms to the total number of events
492 // dir where the output is stored
493 // nEvents Number of events
495 TList* l = GetOutputList(dir);
498 TH2D* density = static_cast<TH2D*>(GetOutputHist(l,"primaryDensity"));
499 if (density) density->Scale(1./nEvents);
503 //____________________________________________________________________