]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/FORWARD/analysis2/AliFMDDensityCalculator.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliFMDDensityCalculator.cxx
CommitLineData
7984e5f7 1// This class calculates the inclusive charged particle density
2// in each for the 5 FMD rings.
3//
7e4038b5 4#include "AliFMDDensityCalculator.h"
5#include <AliESDFMD.h>
6#include <TAxis.h>
7#include <TList.h>
8#include <TMath.h>
0bd4b00f 9#include "AliForwardCorrectionManager.h"
fb3430ac 10#include "AliFMDCorrDoubleHit.h"
11#include "AliFMDCorrELossFit.h"
7e4038b5 12#include "AliLog.h"
77f97e3f 13#include "AliForwardUtil.h"
7e4038b5 14#include <TH2D.h>
0bd4b00f 15#include <TProfile.h>
5bb5d1f6 16#include <THStack.h>
0bd4b00f 17#include <TROOT.h>
5ca83fee 18#include <TVector3.h>
1f7aa5c7 19#include <TStopwatch.h>
bfab35d9 20#include <TParameter.h>
0bd4b00f 21#include <iostream>
22#include <iomanip>
1f7aa5c7 23#include <cstring>
7e4038b5 24
25ClassImp(AliFMDDensityCalculator)
26#if 0
27; // For Emacs
28#endif
29
5ca83fee 30//____________________________________________________________________
31const char* AliFMDDensityCalculator::fgkFolderName = "fmdDensityCalculator";
32
7e4038b5 33//____________________________________________________________________
34AliFMDDensityCalculator::AliFMDDensityCalculator()
35 : TNamed(),
36 fRingHistos(),
dd497217 37 fSumOfWeights(0),
38 fWeightedSum(0),
ea3e5d95 39 fCorrections(0),
0bd4b00f 40 fMaxParticles(5),
9d05ffeb 41 fUsePoisson(false),
8411f7fe 42 fUsePhiAcceptance(kPhiCorrectNch),
0bd4b00f 43 fAccI(0),
44 fAccO(0),
1174780f 45 fFMD1iMax(0),
46 fFMD2iMax(0),
47 fFMD2oMax(0),
48 fFMD3iMax(0),
49 fFMD3oMax(0),
5bb5d1f6 50 fMaxWeights(0),
51 fLowCuts(0),
21d778b1 52 fEtaLumping(32),
53 fPhiLumping(4),
d2638bb7 54 fDebug(0),
6f4a5c0d 55 fCuts(),
5ca83fee 56 fRecalculateEta(false),
8449e3e0 57 fRecalculatePhi(false),
1f7aa5c7 58 fMinQuality(10),
59 fCache(),
60 fDoTiming(false),
77f97e3f
CHC		 61 fHTiming(0),
62 fMaxOutliers(0.05),
63 fOutlierCut(0.50)
7984e5f7 64{
65 //
66 // Constructor
67 //
5ca83fee 68 DGUARD(fDebug, 3, "Default CTOR of FMD density calculator");
7984e5f7 69}
7e4038b5 70
71//____________________________________________________________________
72AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
5ca83fee 73 : TNamed(fgkFolderName, title),
7e4038b5 74 fRingHistos(),
dd497217 75 fSumOfWeights(0),
76 fWeightedSum(0),
ea3e5d95 77 fCorrections(0),
0bd4b00f 78 fMaxParticles(5),
9d05ffeb 79 fUsePoisson(false),
8411f7fe 80 fUsePhiAcceptance(kPhiCorrectNch),
0bd4b00f 81 fAccI(0),
82 fAccO(0),
1174780f 83 fFMD1iMax(0),
84 fFMD2iMax(0),
85 fFMD2oMax(0),
86 fFMD3iMax(0),
87 fFMD3oMax(0),
5bb5d1f6 88 fMaxWeights(0),
89 fLowCuts(0),
21d778b1 90 fEtaLumping(32),
91 fPhiLumping(4),
d2638bb7 92 fDebug(0),
6f4a5c0d 93 fCuts(),
5ca83fee 94 fRecalculateEta(false),
8449e3e0 95 fRecalculatePhi(false),
1f7aa5c7 96 fMinQuality(10),
97 fCache(),
98 fDoTiming(false),
77f97e3f
CHC		 99 fHTiming(0),
100 fMaxOutliers(0.05),
101 fOutlierCut(0.50)
7e4038b5 102{
7984e5f7 103 //
104 // Constructor
105 //
106 // Parameters:
107 // name Name of object
108 //
5ca83fee 109 DGUARD(fDebug, 3, "Named CTOR of FMD density calculator: %s", title);
7e4038b5 110 fRingHistos.SetName(GetName());
e308a636 111 fRingHistos.SetOwner();
7e4038b5 112 fRingHistos.Add(new RingHistos(1, 'I'));
113 fRingHistos.Add(new RingHistos(2, 'I'));
114 fRingHistos.Add(new RingHistos(2, 'O'));
115 fRingHistos.Add(new RingHistos(3, 'I'));
116 fRingHistos.Add(new RingHistos(3, 'O'));
dd497217 117 fSumOfWeights = new TH1D("sumOfWeights", "Sum of Landau weights",
118 200, 0, 20);
0bd4b00f 119 fSumOfWeights->SetFillColor(kRed+1);
120 fSumOfWeights->SetXTitle("#sum_{i} a_{i} f_{i}(#Delta)");
dd497217 121 fWeightedSum = new TH1D("weightedSum", "Weighted sum of Landau propability",
122 200, 0, 20);
0bd4b00f 123 fWeightedSum->SetFillColor(kBlue+1);
124 fWeightedSum->SetXTitle("#sum_{i} i a_{i} f_{i}(#Delta)");
dd497217 125 fCorrections = new TH1D("corrections", "Distribution of corrections",
126 100, 0, 10);
0bd4b00f 127 fCorrections->SetFillColor(kBlue+1);
128 fCorrections->SetXTitle("correction");
dd497217 129
0bd4b00f 130 fAccI = GenerateAcceptanceCorrection('I');
131 fAccO = GenerateAcceptanceCorrection('O');
5bb5d1f6 132
133 fMaxWeights = new TH2D("maxWeights", "Maximum i of a_{i}'s to use",
134 1, 0, 1, 1, 0, 1);
135 fMaxWeights->SetXTitle("#eta");
136 fMaxWeights->SetDirectory(0);
137
138 fLowCuts = new TH2D("lowCuts", "Low cuts used", 1, 0, 1, 1, 0, 1);
139 fLowCuts->SetXTitle("#eta");
140 fLowCuts->SetDirectory(0);
141
7e4038b5 142}
143
144//____________________________________________________________________
145AliFMDDensityCalculator::AliFMDDensityCalculator(const
146 AliFMDDensityCalculator& o)
147 : TNamed(o),
148 fRingHistos(),
dd497217 149 fSumOfWeights(o.fSumOfWeights),
150 fWeightedSum(o.fWeightedSum),
ea3e5d95 151 fCorrections(o.fCorrections),
0bd4b00f 152 fMaxParticles(o.fMaxParticles),
9d05ffeb 153 fUsePoisson(o.fUsePoisson),
5bb5d1f6 154 fUsePhiAcceptance(o.fUsePhiAcceptance),
0bd4b00f 155 fAccI(o.fAccI),
156 fAccO(o.fAccO),
1174780f 157 fFMD1iMax(o.fFMD1iMax),
158 fFMD2iMax(o.fFMD2iMax),
159 fFMD2oMax(o.fFMD2oMax),
160 fFMD3iMax(o.fFMD3iMax),
161 fFMD3oMax(o.fFMD3oMax),
5bb5d1f6 162 fMaxWeights(o.fMaxWeights),
163 fLowCuts(o.fLowCuts),
b6b35c77 164 fEtaLumping(o.fEtaLumping),
165 fPhiLumping(o.fPhiLumping),
d2638bb7 166 fDebug(o.fDebug),
6f4a5c0d 167 fCuts(o.fCuts),
5ca83fee 168 fRecalculateEta(o.fRecalculateEta),
8449e3e0 169 fRecalculatePhi(o.fRecalculatePhi),
1f7aa5c7 170 fMinQuality(o.fMinQuality),
171 fCache(o.fCache),
172 fDoTiming(o.fDoTiming),
77f97e3f
CHC		 173 fHTiming(o.fHTiming),
174 fMaxOutliers(o.fMaxOutliers),
175 fOutlierCut(o.fOutlierCut)
7e4038b5 176{
7984e5f7 177 //
178 // Copy constructor
179 //
180 // Parameters:
181 // o Object to copy from
182 //
5ca83fee 183 DGUARD(fDebug, 3, "Copy CTOR of FMD density calculator");
7e4038b5 184 TIter next(&o.fRingHistos);
185 TObject* obj = 0;
186 while ((obj = next())) fRingHistos.Add(obj);
187}
188
189//____________________________________________________________________
190AliFMDDensityCalculator::~AliFMDDensityCalculator()
191{
7984e5f7 192 //
193 // Destructor
194 //
6ab100ec 195 DGUARD(fDebug, 3, "DTOR of FMD density calculator");
36ffcf83 196 // fRingHistos.Delete();
7e4038b5 197}
198
199//____________________________________________________________________
200AliFMDDensityCalculator&
201AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
202{
7984e5f7 203 //
204 // Assignement operator
205 //
206 // Parameters:
207 // o Object to assign from
208 //
209 // Return:
210 // Reference to this object
211 //
6ab100ec 212 DGUARD(fDebug, 3, "Assignment of FMD density calculator");
d015ecfe 213 if (&o == this) return *this;
ea3e5d95 214 TNamed::operator=(o);
7e4038b5 215
d2638bb7 216 fDebug = o.fDebug;
217 fMaxParticles = o.fMaxParticles;
218 fUsePoisson = o.fUsePoisson;
219 fUsePhiAcceptance = o.fUsePhiAcceptance;
220 fAccI = o.fAccI;
221 fAccO = o.fAccO;
222 fFMD1iMax = o.fFMD1iMax;
223 fFMD2iMax = o.fFMD2iMax;
224 fFMD2oMax = o.fFMD2oMax;
225 fFMD3iMax = o.fFMD3iMax;
226 fFMD3oMax = o.fFMD3oMax;
227 fMaxWeights = o.fMaxWeights;
228 fLowCuts = o.fLowCuts;
229 fEtaLumping = o.fEtaLumping;
230 fPhiLumping = o.fPhiLumping;
231 fCuts = o.fCuts;
6f4a5c0d 232 fRecalculateEta = o.fRecalculateEta;
5ca83fee 233 fRecalculatePhi = o.fRecalculatePhi;
8449e3e0 234 fMinQuality = o.fMinQuality;
1f7aa5c7 235 fCache = o.fCache;
236 fDoTiming = o.fDoTiming;
237 fHTiming = o.fHTiming;
77f97e3f
CHC		 238 fMaxOutliers = o.fMaxOutliers;
239 fOutlierCut = o.fOutlierCut;
d2638bb7 240
7e4038b5 241 fRingHistos.Delete();
242 TIter next(&o.fRingHistos);
243 TObject* obj = 0;
244 while ((obj = next())) fRingHistos.Add(obj);
245
246 return *this;
247}
248
1174780f 249//____________________________________________________________________
250void
5934a3e3 251AliFMDDensityCalculator::SetupForData(const TAxis& axis)
1174780f 252{
253 // Intialize this sub-algorithm
254 //
255 // Parameters:
5ca83fee 256 // etaAxis Eta axis
6ab100ec 257 DGUARD(fDebug, 1, "Initialize FMD density calculator");
e2ebf8c4 258 CacheMaxWeights(axis);
4077f3e8 259
1f7aa5c7 260 fCache.Init(axis);
261
9d05ffeb 262 TIter next(&fRingHistos);
263 RingHistos* o = 0;
d2638bb7 264 while ((o = static_cast<RingHistos*>(next()))) {
5934a3e3 265 o->SetupForData(axis);
e18cb8bd 266 // o->fMultCut = fCuts.GetFixedCut(o->fDet, o->fRing);
267 // o->fPoisson.Init(o->fDet,o->fRing,fEtaLumping, fPhiLumping);
d2638bb7 268 }
1174780f 269}
270
7e4038b5 271//____________________________________________________________________
272AliFMDDensityCalculator::RingHistos*
273AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
274{
7984e5f7 275 //
276 // Get the ring histogram container
277 //
278 // Parameters:
279 // d Detector
280 // r Ring
281 //
282 // Return:
283 // Ring histogram container
284 //
7e4038b5 285 Int_t idx = -1;
286 switch (d) {
287 case 1: idx = 0; break;
288 case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
289 case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
290 }
e308a636 291 if (idx < 0 || idx >= fRingHistos.GetEntries()) {
292 AliWarning(Form("Index %d of FMD%d%c out of range", idx, d, r));
293 return 0;
294 }
7e4038b5 295
296 return static_cast<RingHistos*>(fRingHistos.At(idx));
297}
5bb5d1f6 298
7095962e
CHC		 299namespace {
300 Double_t Rng2Cut(UShort_t d, Char_t r, Int_t xbin, TH2* h)
301 {
302 Double_t ret = 1024;
303 if (xbin < 1 && xbin > h->GetXaxis()->GetNbins()) return ret;
304 Int_t ybin = 0;
305 switch(d) {
306 case 1: ybin = 1; break;
307 case 2: ybin = (r=='i' || r=='I') ? 2 : 3; break;
308 case 3: ybin = (r=='i' || r=='I') ? 4 : 5; break;
309 default: return ret;
310 }
311 ret = h->GetBinContent(xbin,ybin);
312 return ret;
313 }
314}
315
0bd4b00f 316//____________________________________________________________________
317Double_t
5bb5d1f6 318AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Int_t ieta,
7095962e 319 Bool_t /*errors*/) const
0bd4b00f 320{
7984e5f7 321 //
322 // Get the multiplicity cut. If the user has set fMultCut (via
323 // SetMultCut) then that value is used. If not, then the lower
324 // value of the fit range for the enery loss fits is returned.
325 //
326 // Return:
327 // Lower cut on multiplicity
328 //
7095962e
CHC		 329 return Rng2Cut(d, r, ieta, fLowCuts);
330 // return fCuts.GetMultCut(d,r,ieta,errors);
5bb5d1f6 331}
332
333//____________________________________________________________________
334Double_t
335AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Double_t eta,
7095962e 336 Bool_t /*errors*/) const
5bb5d1f6 337{
338 //
339 // Get the multiplicity cut. If the user has set fMultCut (via
340 // SetMultCut) then that value is used. If not, then the lower
341 // value of the fit range for the enery loss fits is returned.
342 //
343 // Return:
344 // Lower cut on multiplicity
345 //
7095962e
CHC		 346 Int_t ieta = fLowCuts->GetXaxis()->FindBin(eta);
347 return Rng2Cut(d, r, ieta, fLowCuts);
348 // return fCuts.GetMultCut(d,r,eta,errors);
0bd4b00f 349}
350
7e4038b5 351//____________________________________________________________________
352Bool_t
353AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd,
354 AliForwardUtil::Histos& hists,
21d778b1 355 Bool_t lowFlux,
b3d81bee 356 Double_t /*cent*/,
5ca83fee 357 const TVector3& ip)
7e4038b5 358{
7984e5f7 359 //
360 // Do the calculations
361 //
362 // Parameters:
363 // fmd AliESDFMD object (possibly) corrected for sharing
364 // hists Histogram cache
365 // vtxBin Vertex bin
366 // lowFlux Low flux flag.
367 //
368 // Return:
369 // true on successs
6ab100ec 370 DGUARD(fDebug, 1, "Calculate density in FMD density calculator");
5ca83fee 371
1f7aa5c7 372 TStopwatch timer;
373 TStopwatch totalT;
374
375 // First measurements of timing
376 // Re-calculation : fraction of sum 32.0% of total 18.1%
377 // N_{particle} : fraction of sum 15.2% of total 8.6%
378 // Correction : fraction of sum 26.4% of total 14.9%
379 // #phi acceptance : fraction of sum 0.2% of total 0.1%
380 // Copy to cache : fraction of sum 3.9% of total 2.2%
381 // Poisson calculation : fraction of sum 18.7% of total 10.6%
382 // Diagnostics : fraction of sum 3.7% of total 2.1%
383 Double_t reEtaTime = 0;
384 Double_t nPartTime = 0;
385 Double_t corrTime = 0;
386 Double_t rePhiTime = 0;
387 Double_t copyTime = 0;
388 Double_t poissonTime= 0;
389 Double_t diagTime = 0;
390 if (fDoTiming) totalT.Start(true);
391
392 Double_t etaCache[20*512]; // Same number of strips per ring
393 Double_t phiCache[20*512]; // whether it is inner our outer.
394 // We do not use TArrayD because we do not wont a bounds check
395 // TArrayD etaCache(20*512); // Same number of strips per ring
396 // TArrayD phiCache(20*512); // whether it is inner our outer.
397
5ca83fee 398 // --- Loop over detectors -----------------------------------------
7e4038b5 399 for (UShort_t d=1; d<=3; d++) {
400 UShort_t nr = (d == 1 ? 1 : 2);
401 for (UShort_t q=0; q<nr; q++) {
402 Char_t r = (q == 0 ? 'I' : 'O');
403 UShort_t ns= (q == 0 ? 20 : 40);
404 UShort_t nt= (q == 0 ? 512 : 256);
405 TH2D* h = hists.Get(d,r);
406 RingHistos* rh= GetRingHistos(d,r);
e308a636 407 if (!rh) {
408 AliError(Form("No ring histogram found for FMD%d%c", d, r));
409 fRingHistos.ls();
410 return false;
411 }
e18cb8bd 412 // rh->fPoisson.SetObject(d,r,vtxbin,cent);
d23503ee 413 rh->fPoisson.Reset(0);
5ca83fee 414 rh->fTotal->Reset();
415 rh->fGood->Reset();
821ffd28 416 // rh->ResetPoissonHistos(h, fEtaLumping, fPhiLumping);
1f7aa5c7 417
418 // Reset our eta cache
419 // for (Int_t i = 0; i < 20*512; i++)
420 // etaCache[i] = phiCache[i] = AliESDFMD::kInvalidEta;
421 memset(etaCache, 0, sizeof(Double_t)*20*512);
422 memset(phiCache, 0, sizeof(Double_t)*20*512);
423 // etaCache.Reset(AliESDFMD::kInvalidEta);
424 // phiCache.Reset(AliESDFMD::kInvalidEta);
425
5ca83fee 426 // --- Loop over sectors and strips ----------------------------
7e4038b5 427 for (UShort_t s=0; s<ns; s++) {
428 for (UShort_t t=0; t<nt; t++) {
3e4a0875 429
5ca83fee 430 Float_t mult = fmd.Multiplicity(d,r,s,t);
431 Float_t phi = fmd.Phi(d,r,s,t) * TMath::DegToRad();
432 Float_t eta = fmd.Eta(d,r,s,t);
433 Double_t oldPhi = phi;
434
435 // --- Re-calculate eta - needed for satelittes ------------
1f7aa5c7 436 if (fDoTiming) timer.Start(true);
437 if (eta == AliESDFMD::kInvalidEta || fRecalculateEta)
5ca83fee 438 eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,ip.Z());
1f7aa5c7 439 if (fDoTiming) reEtaTime += timer.CpuTime();
440 etaCache[s*nt+t] = eta;
5ca83fee 441
442 // --- Check this strip ------------------------------------
443 rh->fTotal->Fill(eta);
bfab35d9 444 if (mult == AliESDFMD::kInvalidMult) { // || mult > 20) {
8449e3e0 445 // Do not count invalid stuff
446 rh->fELoss->Fill(-1);
447 // rh->fEvsN->Fill(mult,-1);
448 // rh->fEvsM->Fill(mult,-1);
33b078ef 449 continue;
450 }
bfab35d9 451 if (mult > 20)
452 AliWarningF("Raw multiplicity of FMD%d%c[%02d,%03d] = %f > 20",
453 d, r, s, t, mult);
5ca83fee 454 // --- Automatic calculation of acceptance -----------------
455 rh->fGood->Fill(eta);
456
457 // --- If we asked to re-calculate phi for (x,y) IP --------
1f7aa5c7 458 if (fDoTiming) timer.Start(true);
5ca83fee 459 if (fRecalculatePhi) {
460 oldPhi = phi;
461 phi = AliForwardUtil::GetPhiFromStrip(r, t, phi, ip.X(), ip.Y());
462 }
1f7aa5c7 463 phiCache[s*nt+t] = phi;
464 if (fDoTiming) rePhiTime += timer.CpuTime();
5ca83fee 465
466 // --- Apply phi corner correction to eloss ----------------
8411f7fe 467 if (fUsePhiAcceptance == kPhiCorrectELoss)
0082a8fc 468 mult *= AcceptanceCorrection(r,t);
5bb5d1f6 469
5ca83fee 470 // --- Get the low multiplicity cut ------------------------
8e400b14 471 Double_t cut = 1024;
472 if (eta != AliESDFMD::kInvalidEta) cut = GetMultCut(d, r, eta,false);
bfab35d9 473 else AliWarningF("Eta for FMD%d%c[%02d,%03d] is invalid: %f",
474 d, r, s, t, eta);
8e400b14 475
5ca83fee 476 // --- Now caluculate Nch for this strip using fits --------
1f7aa5c7 477 if (fDoTiming) timer.Start(true);
5bb5d1f6 478 Double_t n = 0;
5ca83fee 479 if (cut > 0 && mult > cut) n = NParticles(mult,d,r,eta,lowFlux);
f7cfc454 480 rh->fELoss->Fill(mult);
8449e3e0 481 // rh->fEvsN->Fill(mult,n);
482 // rh->fEtaVsN->Fill(eta, n);
1f7aa5c7 483 if (fDoTiming) nPartTime += timer.CpuTime();
9fde7142 484
5ca83fee 485 // --- Calculate correction if needed ----------------------
1f7aa5c7 486 if (fDoTiming) timer.Start(true);
487 // Temporary stuff - remove Correction call
488 Double_t c = 1;
489 if (fUsePhiAcceptance == kPhiCorrectNch)
490 c = AcceptanceCorrection(r,t);
491 // Double_t c = Correction(d,r,t,eta,lowFlux);
492 if (fDoTiming) corrTime += timer.CpuTime();
dd497217 493 fCorrections->Fill(c);
7e4038b5 494 if (c > 0) n /= c;
8449e3e0 495 // rh->fEvsM->Fill(mult,n);
496 // rh->fEtaVsM->Fill(eta, n);
0bd4b00f 497 rh->fCorr->Fill(eta, c);
f7cfc454 498
5ca83fee 499 // --- Accumulate Poisson statistics -----------------------
821ffd28 500 Bool_t hit = (n > 0.9 && c > 0);
5ca83fee 501 if (hit) {
502 rh->fELossUsed->Fill(mult);
503 if (fRecalculatePhi) {
504 rh->fPhiBefore->Fill(oldPhi);
505 rh->fPhiAfter->Fill(phi);
506 }
507 }
21d778b1 508 rh->fPoisson.Fill(t,s,hit,1./c);
7e4038b5 509 h->Fill(eta,phi,n);
5ca83fee 510
511 // --- If we use ELoss fits, apply now ---------------------
5bb5d1f6 512 if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n);
7e4038b5 513 } // for t
514 } // for s
5ca83fee 515
516 // --- Automatic acceptance - Calculate as an efficiency -------
1f7aa5c7 517 // This is very fast, so we do not bother to time it
5ca83fee 518 rh->fGood->Divide(rh->fGood, rh->fTotal, 1, 1, "B");
519
520 // --- Make a copy and reset as needed -------------------------
1f7aa5c7 521 if (fDoTiming) timer.Start(true);
522 TH2D* hclone = fCache.Get(d,r);
523 // hclone->Reset();
524 // TH2D* hclone = static_cast<TH2D*>(h->Clone("hclone"));
21d778b1 525 if (!fUsePoisson) hclone->Reset();
1f7aa5c7 526 else {
527 for (Int_t i = 0; i <= h->GetNbinsX()+1; i++) {
528 for (Int_t j = 0; j <= h->GetNbinsY()+1; j++) {
529 hclone->SetBinContent(i,j,h->GetBinContent(i,j));
530 hclone->SetBinError(i,j,h->GetBinError(i,j));
531 }
532 }
533 // hclone->Add(h);
534 h->Reset();
535 }
536 if (fDoTiming) copyTime += timer.CpuTime();
21d778b1 537
5ca83fee 538 // --- Store Poisson result ------------------------------------
1f7aa5c7 539 if (fDoTiming) timer.Start(true);
821ffd28 540 TH2D* poisson = rh->fPoisson.Result();
d23503ee 541 for (Int_t t=0; t < poisson->GetNbinsX(); t++) {
542 for (Int_t s=0; s < poisson->GetNbinsY(); s++) {
21d778b1 543
544 Double_t poissonV = poisson->GetBinContent(t+1,s+1);
1f7aa5c7 545 // Use cached eta - since the calls to GetEtaFromStrip and
546 // GetPhiFromStrip are _very_ expensive
547 Double_t phi = phiCache[s*nt+t];
548 Double_t eta = etaCache[s*nt+t];
549 // Double_t phi = fmd.Phi(d,r,s,t) * TMath::DegToRad();
550 // Double_t eta = fmd.Eta(d,r,s,t);
551 // if (fRecalculateEta)
552 // eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,ip.Z());
553 // if (fRecalculatePhi)
554 // phi = AliForwardUtil::GetPhiFromStrip(r, t, phi, ip.X(), ip.Y());
555 if (fUsePoisson) {
21d778b1 556 h->Fill(eta,phi,poissonV);
1f7aa5c7 557 rh->fDensity->Fill(eta, phi, poissonV);
558 }
21d778b1 559 else
560 hclone->Fill(eta,phi,poissonV);
21d778b1 561 }
562 }
1f7aa5c7 563 if (fDoTiming) poissonTime += timer.CpuTime();
21d778b1 564
5ca83fee 565 // --- Make diagnostics - eloss vs poisson ---------------------
1f7aa5c7 566 if (fDoTiming) timer.Start(true);
5ca83fee 567 Int_t nY = h->GetNbinsY();
77f97e3f
CHC		 568 Int_t nIn = 0; // Count non-outliers
569 Int_t nOut = 0; // Count outliers
21d778b1 570 for (Int_t ieta=1; ieta <= h->GetNbinsX(); ieta++) {
5ca83fee 571 // Set the overflow bin to contain the phi acceptance
572 Double_t phiAcc = rh->fGood->GetBinContent(ieta);
573 Double_t phiAccE = rh->fGood->GetBinError(ieta);
574 h->SetBinContent(ieta, nY+1, phiAcc);
575 h->SetBinError(ieta, nY+1, phiAccE);
576 Double_t eta = h->GetXaxis()->GetBinCenter(ieta);
577 rh->fPhiAcc->Fill(eta, ip.Z(), phiAcc);
5ca83fee 578 for (Int_t iphi=1; iphi<= nY; iphi++) {
e83d0620 579
21d778b1 580 Double_t poissonV = 0; //h->GetBinContent(,s+1);
581 Double_t eLossV = 0;
582 if(fUsePoisson) {
583 poissonV = h->GetBinContent(ieta,iphi);
584 eLossV = hclone->GetBinContent(ieta,iphi);
585 }
586 else {
587 poissonV = hclone->GetBinContent(ieta,iphi);
588 eLossV = h->GetBinContent(ieta,iphi);
589 }
7b03929c 590
77f97e3f
CHC		 591 if (poissonV < 1e-12 && eLossV < 1e-12)
592 // we do not care about trivially empty bins
593 continue;
5ca83fee 594
77f97e3f
CHC		 595 Bool_t outlier = CheckOutlier(eLossV, poissonV, fOutlierCut);
596 Double_t rel = eLossV < 1e-12 ? 0 : (poissonV - eLossV) / eLossV;
597 if (outlier) {
598 rh->fELossVsPoissonOut->Fill(eLossV, poissonV);
599 rh->fDiffELossPoissonOut->Fill(rel);
600 nOut++;
9d05ffeb 601 }
77f97e3f
CHC		 602 else {
603 rh->fELossVsPoisson->Fill(eLossV, poissonV);
604 rh->fDiffELossPoisson->Fill(rel);
605 nIn++;
606 } // if (outlier)
607 } // for (iphi)
608 } // for (ieta)
67a4bb96
CHC		 609 Int_t nTotal = (nIn+nOut);
610 Double_t outRatio = (nTotal > 0 ? Double_t(nOut) / nTotal : 0);
77f97e3f
CHC		 611 rh->fOutliers->Fill(outRatio);
612 if (outRatio < fMaxOutliers) rh->fPoisson.FillDiagnostics();
613 else h->SetBit(AliForwardUtil::kSkipRing);
1f7aa5c7 614 if (fDoTiming) diagTime += timer.CpuTime();
615 // delete hclone;
21d778b1 616
7e4038b5 617 } // for q
618 } // for d
a7f31ad8 619
1f7aa5c7 620 if (fDoTiming) {
621 fHTiming->Fill(1,reEtaTime);
622 fHTiming->Fill(2,nPartTime);
623 fHTiming->Fill(3,corrTime);
624 fHTiming->Fill(4,rePhiTime);
625 fHTiming->Fill(5,copyTime);
626 fHTiming->Fill(6,poissonTime);
627 fHTiming->Fill(7,diagTime);
628 fHTiming->Fill(8,totalT.CpuTime());
629 }
630
7e4038b5 631 return kTRUE;
632}
633
77f97e3f
CHC		 634//_____________________________________________________________________
635Bool_t
636AliFMDDensityCalculator::CheckOutlier(Double_t eloss,
637 Double_t poisson,
638 Double_t cut) const
639{
640 if (eloss < 1e-6) return true;
641 Double_t diff = TMath::Abs(poisson - eloss) / eloss;
642 return diff > cut;
643}
1174780f 644//_____________________________________________________________________
645Int_t
fb3430ac 646AliFMDDensityCalculator::FindMaxWeight(const AliFMDCorrELossFit* cor,
1174780f 647 UShort_t d, Char_t r, Int_t iEta) const
648{
fb3430ac 649 //
650 // Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
651 //
652 // Parameters:
653 // cor Correction
654 // d Detector
655 // r Ring
656 // iEta Eta bin
657 //
241cca4d 658 DGUARD(fDebug, 10, "Find maximum weight in FMD density calculator");
f61b1745 659 if(!cor) return -1;
660
8449e3e0 661 AliFMDCorrELossFit::ELossFit* fit = cor->FindFit(d,r,iEta, -1);
1174780f 662 if (!fit) {
663 // AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
664 return -1;
665 }
8449e3e0 666 return fit->FindMaxWeight(2*AliFMDCorrELossFit::ELossFit::fgMaxRelError,
667 AliFMDCorrELossFit::ELossFit::fgLeastWeight,
668 fMaxParticles);
1174780f 669}
7095962e
CHC		 670//_____________________________________________________________________
671Int_t
672AliFMDDensityCalculator::FindMaxWeight(const AliFMDCorrELossFit* cor,
673 UShort_t d, Char_t r, Double_t eta) const
674{
675 //
676 // Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
677 //
678 // Parameters:
679 // cor Correction
680 // d Detector
681 // r Ring
682 // eta Eta
683 //
684 DGUARD(fDebug, 10, "Find maximum weight in FMD density calculator");
685 if(!cor) return -1;
686
687 AliFMDCorrELossFit::ELossFit* fit = cor->FindFit(d,r,eta, -1);
688 if (!fit) {
689 // AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
690 return -1;
691 }
692 return fit->FindMaxWeight(2*AliFMDCorrELossFit::ELossFit::fgMaxRelError,
693 AliFMDCorrELossFit::ELossFit::fgLeastWeight,
694 fMaxParticles);
695}
1174780f 696
697//_____________________________________________________________________
698void
e2ebf8c4 699AliFMDDensityCalculator::CacheMaxWeights(const TAxis& axis)
1174780f 700{
fb3430ac 701 //
702 // Find the max weights and cache them
703 //
6ab100ec 704 DGUARD(fDebug, 2, "Cache maximum weights in FMD density calculator");
1174780f 705 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
8449e3e0 706 const AliFMDCorrELossFit* cor = fcm.GetELossFit();
707 cor->CacheBins(fMinQuality);
7095962e 708 cor->Print(fDebug > 5 ? "RCS" : "C");
a76fb27d 709
710 TAxis eta(axis.GetNbins(),
711 axis.GetXmin(),
712 axis.GetXmax());
713
f61b1745 714 if(cor)
715 eta.Set(cor->GetEtaAxis().GetNbins(),
716 cor->GetEtaAxis().GetXmin(),
717 cor->GetEtaAxis().GetXmax());
1174780f 718
719 Int_t nEta = eta.GetNbins();
720 fFMD1iMax.Set(nEta);
721 fFMD2iMax.Set(nEta);
722 fFMD2oMax.Set(nEta);
723 fFMD3iMax.Set(nEta);
724 fFMD3oMax.Set(nEta);
725
5bb5d1f6 726 fMaxWeights->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
727 fMaxWeights->GetYaxis()->SetBinLabel(1, "FMD1i");
728 fMaxWeights->GetYaxis()->SetBinLabel(2, "FMD2i");
729 fMaxWeights->GetYaxis()->SetBinLabel(3, "FMD2o");
730 fMaxWeights->GetYaxis()->SetBinLabel(4, "FMD3i");
731 fMaxWeights->GetYaxis()->SetBinLabel(5, "FMD3o");
732
733 AliInfo(Form("Get eta axis with %d bins from %f to %f",
734 nEta, eta.GetXmin(), eta.GetXmax()));
735 fLowCuts->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
736 fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
737 fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
738 fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
739 fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
740 fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
7095962e 741
1174780f 742 for (Int_t i = 0; i < nEta; i++) {
7095962e 743 Double_t leta = fMaxWeights->GetXaxis()->GetBinCenter(i+1);
5bb5d1f6 744 Double_t w[5];
7095962e
CHC		 745 w[0] = fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', leta);
746 w[1] = fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', leta);
747 w[2] = fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', leta);
748 w[3] = fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', leta);
749 w[4] = fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', leta);
750 for (Int_t j = 0; j < 5; j++)
5bb5d1f6 751 if (w[j] > 0) fMaxWeights->SetBinContent(i+1, j+1, w[j]);
1174780f 752 }
7095962e
CHC		 753
754 // Cache cuts in histogram
755 fCuts.FillHistogram(fLowCuts);
1174780f 756}
757
758//_____________________________________________________________________
759Int_t
760AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const
761{
fb3430ac 762 //
763 // Find the (cached) maximum weight for FMD<i>dr</i> in
764 // @f$\eta@f$ bin @a iEta
765 //
766 // Parameters:
767 // d Detector
768 // r Ring
769 // iEta Eta bin
770 //
771 // Return:
772 // max weight or <= 0 in case of problems
773 //
1174780f 774 if (iEta < 0) return -1;
775
776 const TArrayI* max = 0;
777 switch (d) {
778 case 1: max = &fFMD1iMax; break;
779 case 2: max = (r == 'I' || r == 'i' ? &fFMD2iMax : &fFMD2oMax); break;
780 case 3: max = (r == 'I' || r == 'i' ? &fFMD3iMax : &fFMD3oMax); break;
781 }
782 if (!max) {
783 AliWarning(Form("No array for FMD%d%c", d, r));
784 return -1;
785 }
786
787 if (iEta >= max->fN) {
788 AliWarning(Form("Eta bin %3d out of bounds [0,%d]",
789 iEta, max->fN-1));
790 return -1;
791 }
792
793 AliDebug(30,Form("Max weight for FMD%d%c eta bin %3d: %d", d, r, iEta,
794 max->At(iEta)));
795 return max->At(iEta);
796}
797
798//_____________________________________________________________________
799Int_t
800AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const
801{
fb3430ac 802 //
803 // Find the (cached) maximum weight for FMD<i>dr</i> iat
804 // @f$\eta@f$
805 //
806 // Parameters:
807 // d Detector
808 // r Ring
809 // eta Eta bin
810 //
811 // Return:
812 // max weight or <= 0 in case of problems
813 //
1174780f 814 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
815 Int_t iEta = fcm.GetELossFit()->FindEtaBin(eta) -1;
816
817 return GetMaxWeight(d, r, iEta);
818}
819
7e4038b5 820//_____________________________________________________________________
821Float_t
822AliFMDDensityCalculator::NParticles(Float_t mult,
823 UShort_t d,
824 Char_t r,
7e4038b5 825 Float_t eta,
826 Bool_t lowFlux) const
827{
7984e5f7 828 //
829 // Get the number of particles corresponding to the signal mult
830 //
831 // Parameters:
832 // mult Signal
833 // d Detector
834 // r Ring
835 // s Sector
836 // t Strip (not used)
837 // v Vertex bin
838 // eta Pseudo-rapidity
839 // lowFlux Low-flux flag
840 //
841 // Return:
842 // The number of particles
843 //
5bb5d1f6 844 // if (mult <= GetMultCut()) return 0;
6ab100ec 845 DGUARD(fDebug, 3, "Calculate Nch in FMD density calculator");
7e4038b5 846 if (lowFlux) return 1;
0bd4b00f 847
0bd4b00f 848 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
8449e3e0 849 AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta, -1);
0bd4b00f 850 if (!fit) {
8449e3e0 851 AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f qual=%d",
852 d, r, eta, fMinQuality));
0bd4b00f 853 return 0;
854 }
855
1174780f 856 Int_t m = GetMaxWeight(d,r,eta); // fit->FindMaxWeight();
0bd4b00f 857 if (m < 1) {
858 AliWarning(Form("No good fits for FMD%d%c at eta=%f", d, r, eta));
859 return 0;
860 }
79909b8b 861
0bd4b00f 862 UShort_t n = TMath::Min(fMaxParticles, UShort_t(m));
863 Double_t ret = fit->EvaluateWeighted(mult, n);
79909b8b 864
0bd4b00f 865 if (fDebug > 10) {
866 AliInfo(Form("FMD%d%c, eta=%7.4f, %8.5f -> %8.5f", d, r, eta, mult, ret));
867 }
79909b8b 868
0bd4b00f 869 fWeightedSum->Fill(ret);
870 fSumOfWeights->Fill(ret);
3e4a0875 871
0bd4b00f 872 return ret;
7e4038b5 873}
874
875//_____________________________________________________________________
876Float_t
0bd4b00f 877AliFMDDensityCalculator::Correction(UShort_t d,
878 Char_t r,
0bd4b00f 879 UShort_t t,
0bd4b00f 880 Float_t eta,
881 Bool_t lowFlux) const
7e4038b5 882{
7984e5f7 883 //
884 // Get the inverse correction factor. This consist of
885 //
886 // - acceptance correction (corners of sensors)
887 // - double hit correction (for low-flux events)
888 // - dead strip correction
889 //
890 // Parameters:
891 // d Detector
892 // r Ring
893 // s Sector
894 // t Strip (not used)
895 // v Vertex bin
896 // eta Pseudo-rapidity
897 // lowFlux Low-flux flag
898 //
899 // Return:
900 //
901 //
241cca4d 902 DGUARD(fDebug, 10, "Apply correction in FMD density calculator");
5bb5d1f6 903 Float_t correction = 1;
8411f7fe 904 if (fUsePhiAcceptance == kPhiCorrectNch)
905 correction = AcceptanceCorrection(r,t);
7e4038b5 906 if (lowFlux) {
1f7aa5c7 907 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
908
0bd4b00f 909 TH1D* dblHitCor = 0;
910 if (fcm.GetDoubleHit())
911 dblHitCor = fcm.GetDoubleHit()->GetCorrection(d,r);
912
7e4038b5 913 if (dblHitCor) {
0bd4b00f 914 Double_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
915 if (dblC > 0) correction *= dblC;
7e4038b5 916 }
57522224 917 // else {
918 // AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
919 // }
7e4038b5 920 }
7e4038b5 921 return correction;
922}
923
0bd4b00f 924//_____________________________________________________________________
925TH1D*
926AliFMDDensityCalculator::GenerateAcceptanceCorrection(Char_t r) const
927{
7984e5f7 928 //
929 // Generate the acceptance corrections
930 //
931 // Parameters:
932 // r Ring to generate for
933 //
934 // Return:
935 // Newly allocated histogram of acceptance corrections
936 //
6ab100ec 937 DGUARD(fDebug, 3, "Make acceptance correction in FMD density calculator");
0bd4b00f 938 const Double_t ic1[] = { 4.9895, 15.3560 };
939 const Double_t ic2[] = { 1.8007, 17.2000 };
940 const Double_t oc1[] = { 4.2231, 26.6638 };
941 const Double_t oc2[] = { 1.8357, 27.9500 };
942 const Double_t* c1 = (r == 'I' || r == 'i' ? ic1 : oc1);
943 const Double_t* c2 = (r == 'I' || r == 'i' ? ic2 : oc2);
944 Double_t minR = (r == 'I' || r == 'i' ? 4.5213 : 15.4);
945 Double_t maxR = (r == 'I' || r == 'i' ? 17.2 : 28.0);
946 Int_t nStrips = (r == 'I' || r == 'i' ? 512 : 256);
947 Int_t nSec = (r == 'I' || r == 'i' ? 20 : 40);
948 Float_t basearc = 2 * TMath::Pi() / nSec;
949 Double_t rad = maxR - minR;
950 Float_t segment = rad / nStrips;
951 Float_t cr = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
952
953 // Numbers used to find end-point of strip.
954 // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
955 Float_t D = c1[0] * c2[1] - c1[1] * c2[0];
956 Float_t dx = c2[0] - c1[0];
957 Float_t dy = c2[1] - c1[1];
958 Float_t dr = TMath::Sqrt(dx*dx+dy*dy);
959
960 TH1D* ret = new TH1D(Form("acc%c", r),
961 Form("Acceptance correction for FMDx%c", r),
962 nStrips, -.5, nStrips-.5);
963 ret->SetXTitle("Strip");
964 ret->SetYTitle("#varphi acceptance");
965 ret->SetDirectory(0);
966 ret->SetFillColor(r == 'I' || r == 'i' ? kRed+1 : kBlue+1);
967 ret->SetFillStyle(3001);
968
969 for (Int_t t = 0; t < nStrips; t++) {
970 Float_t radius = minR + t * segment;
971
972 // If the radius of the strip is smaller than the radius corresponding
973 // to the first corner we have a full strip length
974 if (radius <= cr) {
975 ret->SetBinContent(t+1, 1);
976 continue;
977 }
978
979 // Next, we should find the end-point of the strip - that is,
980 // the coordinates where the line from c1 to c2 intersects a circle
981 // with radius given by the strip.
982 // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
983 // Calculate the determinant
984 Float_t det = radius * radius * dr * dr - D*D;
985
986 if (det <= 0) {
987 // <0 means No intersection
988 // =0 means Exactly tangent
989 ret->SetBinContent(t+1, 1);
990 continue;
991 }
992
993 // Calculate end-point and the corresponding opening angle
994 Float_t x = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
995 Float_t y = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
996 Float_t th = TMath::ATan2(x, y);
997
998 ret->SetBinContent(t+1, th / basearc);
999 }
1000 return ret;
1001}
7e4038b5 1002
1003//_____________________________________________________________________
1004Float_t
1005AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
1006{
7984e5f7 1007 //
1008 // Get the acceptance correction for strip @a t in an ring of type @a r
1009 //
1010 // Parameters:
1011 // r Ring type ('I' or 'O')
1012 // t Strip number
1013 //
1014 // Return:
1015 // Inverse acceptance correction
1016 //
0bd4b00f 1017 TH1D* acc = (r == 'I' || r == 'i' ? fAccI : fAccO);
1018 return acc->GetBinContent(t+1);
7e4038b5 1019}
1020
1021//____________________________________________________________________
1022void
5934a3e3 1023AliFMDDensityCalculator::Terminate(const TList* dir, TList* output,
1024 Int_t nEvents)
7e4038b5 1025{
7984e5f7 1026 //
1027 // Scale the histograms to the total number of events
1028 //
1029 // Parameters:
1030 // dir where to put the output
1031 // nEvents Number of events
1032 //
6ab100ec 1033 DGUARD(fDebug, 1, "Scale histograms in FMD density calculator");
7e4038b5 1034 if (nEvents <= 0) return;
9d99b0dd 1035 TList* d = static_cast<TList*>(dir->FindObject(GetName()));
1036 if (!d) return;
7e4038b5 1037
5934a3e3 1038 TList* out = new TList;
1039 out->SetName(d->GetName());
1040 out->SetOwner();
1041
7e4038b5 1042 TIter next(&fRingHistos);
1043 RingHistos* o = 0;
5bb5d1f6 1044 THStack* sums = new THStack("sums", "sums of ring signals");
1045 while ((o = static_cast<RingHistos*>(next()))) {
5934a3e3 1046 o->Terminate(d, nEvents);
1047 if (!o->fDensity) {
1048 Warning("Terminate", "No density in %s", o->GetName());
1049 continue;
1050 }
5bb5d1f6 1051 TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
1052 o->fDensity->GetNbinsY(),"e");
1053 sum->Scale(1., "width");
1054 sum->SetTitle(o->GetName());
1055 sum->SetDirectory(0);
1056 sum->SetYTitle("#sum N_{ch,incl}");
1057 sums->Add(sum);
1058 }
5934a3e3 1059 out->Add(sums);
1060 output->Add(out);
7e4038b5 1061}
1062
1063//____________________________________________________________________
1064void
5934a3e3 1065AliFMDDensityCalculator::CreateOutputObjects(TList* dir)
7e4038b5 1066{
7984e5f7 1067 //
1068 // Output diagnostic histograms to directory
1069 //
1070 // Parameters:
1071 // dir List to write in
1072 //
6ab100ec 1073 DGUARD(fDebug, 1, "Define output FMD density calculator");
7e4038b5 1074 TList* d = new TList;
5bb5d1f6 1075 d->SetOwner();
7e4038b5 1076 d->SetName(GetName());
1077 dir->Add(d);
dd497217 1078 d->Add(fWeightedSum);
1079 d->Add(fSumOfWeights);
1080 d->Add(fCorrections);
0bd4b00f 1081 d->Add(fAccI);
1082 d->Add(fAccO);
5bb5d1f6 1083 d->Add(fMaxWeights);
1084 d->Add(fLowCuts);
0bd4b00f 1085
bfab35d9 1086 TParameter<int>* nFiles = new TParameter<int>("nFiles", 1);
1087 nFiles->SetMergeMode('+');
1088
d2638bb7 1089 // d->Add(sigma);
77f97e3f 1090 d->Add(AliForwardUtil::MakeParameter("maxParticle", fMaxParticles));
7095962e 1091 d->Add(AliForwardUtil::MakeParameter("minQuality", fMinQuality));
77f97e3f
CHC		 1092 d->Add(AliForwardUtil::MakeParameter("method", fUsePoisson));
1093 d->Add(AliForwardUtil::MakeParameter("phiAcceptance",fUsePhiAcceptance));
1094 d->Add(AliForwardUtil::MakeParameter("etaLumping", fEtaLumping));
1095 d->Add(AliForwardUtil::MakeParameter("phiLumping", fPhiLumping));
1096 d->Add(AliForwardUtil::MakeParameter("recalcEta", fRecalculateEta));
1097 d->Add(AliForwardUtil::MakeParameter("recalcPhi", fRecalculatePhi));
1098 d->Add(AliForwardUtil::MakeParameter("maxOutliers", fMaxOutliers));
1099 d->Add(AliForwardUtil::MakeParameter("outlierCut", fOutlierCut));
bfab35d9 1100 d->Add(nFiles);
d2638bb7 1101 // d->Add(nxi);
5ca83fee 1102 fCuts.Output(d,"lCuts");
f7cfc454 1103
7e4038b5 1104 TIter next(&fRingHistos);
1105 RingHistos* o = 0;
1106 while ((o = static_cast<RingHistos*>(next()))) {
e18cb8bd 1107 o->fPoisson.SetLumping(fEtaLumping, fPhiLumping);
5934a3e3 1108 o->CreateOutputObjects(d);
7e4038b5 1109 }
1f7aa5c7 1110
1111 if (!fDoTiming) return;
1112
1113 fHTiming = new TProfile("timing", "#LTt_{part}#GT", 8, .5, 8.5);
1114 fHTiming->SetDirectory(0);
1115 fHTiming->SetYTitle("#LTt_{part}#GT");
1116 fHTiming->SetXTitle("Part");
1117 fHTiming->SetFillColor(kRed+1);
1118 fHTiming->SetFillStyle(3001);
1119 fHTiming->SetMarkerStyle(20);
1120 fHTiming->SetMarkerColor(kBlack);
1121 fHTiming->SetLineColor(kBlack);
1122 fHTiming->SetStats(0);
1123 TAxis* xaxis = fHTiming->GetXaxis();
1124 xaxis->SetBinLabel(1, "Re-calculation of #eta");
1125 xaxis->SetBinLabel(2, "N_{particle}");
1126 xaxis->SetBinLabel(3, "Correction");
1127 xaxis->SetBinLabel(4, "Re-calculation of #phi");
1128 xaxis->SetBinLabel(5, "Copy to cache");
1129 xaxis->SetBinLabel(6, "Poisson calculation");
1130 xaxis->SetBinLabel(7, "Diagnostics");
1131 xaxis->SetBinLabel(8, "Total");
1132 d->Add(fHTiming);
7e4038b5 1133}
c8b1a7db 1134#define PF(N,V,...) \
1135 AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
1136#define PFB(N,FLAG) \
1137 do { \
1138 AliForwardUtil::PrintName(N); \
1139 std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
1140 } while(false)
1141#define PFV(N,VALUE) \
1142 do { \
1143 AliForwardUtil::PrintName(N); \
1144 std::cout << (VALUE) << std::endl; } while(false)
1145
0bd4b00f 1146//____________________________________________________________________
1147void
1174780f 1148AliFMDDensityCalculator::Print(Option_t* option) const
0bd4b00f 1149{
7984e5f7 1150 //
1151 // Print information
1152 //
1153 // Parameters:
1154 // option Not used
1155 //
c8b1a7db 1156 AliForwardUtil::PrintTask(*this);
1157 gROOT->IncreaseDirLevel();
1158
1159 TString phiM("none");
5ca83fee 1160 switch (fUsePhiAcceptance) {
c8b1a7db 1161 case kPhiNoCorrect: phiM = "none"; break;
1162 case kPhiCorrectNch: phiM = "correct Nch"; break;
1163 case kPhiCorrectELoss: phiM = "correct energy loss"; break;
5ca83fee 1164 }
c8b1a7db 1165
1166 PFV("Max(particles)", fMaxParticles );
7095962e 1167 PFV("Min(quality)", fMinQuality );
c8b1a7db 1168 PFV("Poisson method", fUsePoisson );
1169 PFV("Eta lumping", fEtaLumping );
1170 PFV("Phi lumping", fPhiLumping );
1171 PFV("Recalculate eta", fRecalculateEta );
1172 PFV("Recalculate phi", fRecalculatePhi );
1173 PFV("Use phi acceptance", phiM);
1174 PFV("Lower cut", "");
d2638bb7 1175 fCuts.Print();
c8b1a7db 1176
1174780f 1177 TString opt(option);
1178 opt.ToLower();
c8b1a7db 1179 if (!opt.Contains("nomax")) {
1180 PFV("Max weights", "");
1181
1182 char ind[64];
1183 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
1184 ind[gROOT->GetDirLevel()] = '\0';
1185 for (UShort_t d=1; d<=3; d++) {
1186 UShort_t nr = (d == 1 ? 1 : 2);
1187 for (UShort_t q=0; q<nr; q++) {
1188 ind[gROOT->GetDirLevel()] = ' ';
1189 ind[gROOT->GetDirLevel()+1] = '\0';
1190 Char_t r = (q == 0 ? 'I' : 'O');
1191 std::cout << ind << " FMD" << d << r << ":";
1192 ind[gROOT->GetDirLevel()+1] = ' ';
1193 ind[gROOT->GetDirLevel()+2] = '\0';
1194
1195 const TArrayI& a = (d == 1 ? fFMD1iMax :
1196 (d == 2 ? (r == 'I' ? fFMD2iMax : fFMD2oMax) :
1197 (r == 'I' ? fFMD3iMax : fFMD3oMax)));
1198 Int_t j = 0;
1199 for (Int_t i = 0; i < a.fN; i++) {
1200 if (a.fArray[i] < 1) continue;
1201 if (j % 6 == 0) std::cout << "\n " << ind;
1202 j++;
1203 std::cout << " " << std::setw(3) << i << ": " << a.fArray[i];
1204 }
1205 std::cout << std::endl;
1174780f 1206 }
1174780f 1207 }
1208 }
c8b1a7db 1209 gROOT->DecreaseDirLevel();
0bd4b00f 1210}
7e4038b5 1211
1212//====================================================================
1213AliFMDDensityCalculator::RingHistos::RingHistos()
9d99b0dd 1214 : AliForwardUtil::RingHistos(),
5ca83fee 1215 fList(0),
8449e3e0 1216 // fEvsN(0),
1217 // fEvsM(0),
1218 // fEtaVsN(0),
1219 // fEtaVsM(0),
0bd4b00f 1220 fCorr(0),
9d05ffeb 1221 fDensity(0),
1222 fELossVsPoisson(0),
5ca83fee 1223 fDiffELossPoisson(0),
77f97e3f
CHC		 1224 fELossVsPoissonOut(0),
1225 fDiffELossPoissonOut(0),
1226 fOutliers(0),
821ffd28 1227 fPoisson(),
f7cfc454 1228 fELoss(0),
1229 fELossUsed(0),
5ca83fee 1230 fMultCut(0),
1231 fTotal(0),
1232 fGood(0),
1233 fPhiAcc(0),
1234 fPhiBefore(0),
1235 fPhiAfter(0)
7984e5f7 1236{
1237 //
1238 // Default CTOR
1239 //
1240}
7e4038b5 1241
1242//____________________________________________________________________
1243AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
9d99b0dd 1244 : AliForwardUtil::RingHistos(d,r),
5ca83fee 1245 fList(0),
8449e3e0 1246 // fEvsN(0),
1247 // fEvsM(0),
1248 // fEtaVsN(0),
1249 // fEtaVsM(0),
0bd4b00f 1250 fCorr(0),
9d05ffeb 1251 fDensity(0),
1252 fELossVsPoisson(0),
5ca83fee 1253 fDiffELossPoisson(0),
77f97e3f
CHC		 1254 fELossVsPoissonOut(0),
1255 fDiffELossPoissonOut(0),
1256 fOutliers(0),
821ffd28 1257 fPoisson("ignored"),
d2638bb7 1258 fELoss(0),
f7cfc454 1259 fELossUsed(0),
5ca83fee 1260 fMultCut(0),
1261 fTotal(0),
1262 fGood(0),
1263 fPhiAcc(0),
1264 fPhiBefore(0),
1265 fPhiAfter(0)
7e4038b5 1266{
7984e5f7 1267 //
1268 // Constructor
1269 //
1270 // Parameters:
1271 // d detector
1272 // r ring
1273 //
8449e3e0 1274#if 0
5bb5d1f6 1275 fEvsN = new TH2D("elossVsNnocorr",
1276 "#Delta E/#Delta E_{mip} vs uncorrected inclusive N_{ch}",
5ca83fee 1277 250, -.5, 24.5, 251, -1.5, 24.5);
7e4038b5 1278 fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
1279 fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
1280 fEvsN->Sumw2();
1281 fEvsN->SetDirectory(0);
5bb5d1f6 1282
1283 fEvsM = static_cast<TH2D*>(fEvsN->Clone("elossVsNcorr"));
1284 fEvsM->SetTitle("#Delta E/#Delta E_{mip} vs corrected inclusive N_{ch}");
7e4038b5 1285 fEvsM->SetDirectory(0);
1286
5bb5d1f6 1287 fEtaVsN = new TProfile("etaVsNnocorr",
1288 "Average inclusive N_{ch} vs #eta (uncorrected)",
1289 200, -4, 6);
0bd4b00f 1290 fEtaVsN->SetXTitle("#eta");
1291 fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
1292 fEtaVsN->SetDirectory(0);
1293 fEtaVsN->SetLineColor(Color());
1294 fEtaVsN->SetFillColor(Color());
5bb5d1f6 1295
1296 fEtaVsM = static_cast<TProfile*>(fEtaVsN->Clone("etaVsNcorr"));
1297 fEtaVsM->SetTitle("Average inclusive N_{ch} vs #eta (corrected)");
0bd4b00f 1298 fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
1299 fEtaVsM->SetDirectory(0);
8449e3e0 1300#endif
0bd4b00f 1301
5bb5d1f6 1302 fCorr = new TProfile("corr", "Average correction", 200, -4, 6);
0bd4b00f 1303 fCorr->SetXTitle("#eta");
1304 fCorr->SetYTitle("#LT correction#GT");
1305 fCorr->SetDirectory(0);
1306 fCorr->SetLineColor(Color());
1307 fCorr->SetFillColor(Color());
1308
5bb5d1f6 1309 fDensity = new TH2D("inclDensity", "Inclusive N_{ch} density",
7e4038b5 1310 200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
1311 0, 2*TMath::Pi());
1312 fDensity->SetDirectory(0);
5ca83fee 1313 fDensity->Sumw2(); fDensity->SetMarkerColor(Color());
7e4038b5 1314 fDensity->SetXTitle("#eta");
1315 fDensity->SetYTitle("#phi [radians]");
1316 fDensity->SetZTitle("Inclusive N_{ch} density");
9d05ffeb 1317
77f97e3f
CHC		 1318 // --- Create increasing sized bins --------------------------------
1319 TArrayD bins;
1320 // bins, lowest order, higest order, return array
1321 const char* nchP = "N_{ch}^{Poisson}";
1322 const char* nchE = "N_{ch}^{#Delta}";
1323 AliForwardUtil::MakeLogScale(300, 0, 2, bins);
5bb5d1f6 1324 fELossVsPoisson = new TH2D("elossVsPoisson",
77f97e3f
CHC		 1325 "N_{ch} from energy loss vs from Poisson",
1326 bins.GetSize()-1, bins.GetArray(),
1327 bins.GetSize()-1, bins.GetArray());
9d05ffeb 1328 fELossVsPoisson->SetDirectory(0);
77f97e3f
CHC		 1329 fELossVsPoisson->SetXTitle(nchE);
1330 fELossVsPoisson->SetYTitle(nchP);
9d05ffeb 1331 fELossVsPoisson->SetZTitle("Correlation");
77f97e3f
CHC		 1332 fELossVsPoissonOut =
1333 static_cast<TH2D*>(fELossVsPoisson
1334 ->Clone(Form("%sOutlier",
1335 fELossVsPoisson->GetName())));
1336 fELossVsPoissonOut->SetDirectory(0);
1337 fELossVsPoissonOut->SetMarkerStyle(20);
1338 fELossVsPoissonOut->SetMarkerSize(0.3);
1339 fELossVsPoissonOut->SetMarkerColor(kBlack);
1340 fELossVsPoissonOut->SetTitle(Form("%s for outliers",
1341 fELossVsPoisson->GetTitle()));
9d05ffeb 1342
5ca83fee 1343 fDiffELossPoisson = new TH1D("diffElossPoisson",
77f97e3f 1344 Form("(%s-%s)/%s", nchP, nchE, nchE),
5ca83fee 1345 100, -1, 1);
1346 fDiffELossPoisson->SetDirectory(0);
1347 fDiffELossPoisson->SetXTitle(fDiffELossPoisson->GetTitle());
1348 fDiffELossPoisson->SetYTitle("Frequency");
1349 fDiffELossPoisson->SetMarkerColor(Color());
1350 fDiffELossPoisson->SetFillColor(Color());
bfab35d9 1351 fDiffELossPoisson->SetFillStyle(3001);
5ca83fee 1352 fDiffELossPoisson->Sumw2();
1353
77f97e3f
CHC		 1354 fDiffELossPoissonOut =
1355 static_cast<TH1D*>(fDiffELossPoisson
1356 ->Clone(Form("%sOutlier",fDiffELossPoisson->GetName())));
1357 fDiffELossPoissonOut->SetDirectory(0);
1358 fDiffELossPoissonOut->SetTitle(Form("%s for outliers",
1359 fDiffELossPoisson->GetTitle()));
1360 fDiffELossPoissonOut->SetMarkerColor(Color()-2);
1361 fDiffELossPoissonOut->SetFillColor(Color()-2);
1362 fDiffELossPoissonOut->SetFillStyle(3002);
1363
1364 fOutliers = new TH1D("outliers", "Fraction of outliers", 100, 0, 1);
1365 fOutliers->SetDirectory(0);
1366 fOutliers->SetXTitle("N_{outlier}/(N_{outlier}+N_{inside})");
1367 fOutliers->SetYTitle("#sum_{events}#sum_{bins}");
1368 fOutliers->SetFillColor(Color());
1369 fOutliers->SetFillStyle(3001);
1370 fOutliers->SetLineColor(kBlack);
1371
f7cfc454 1372 fELoss = new TH1D("eloss", "#Delta/#Delta_{mip} in all strips",
8449e3e0 1373 640, -1, 15);
f7cfc454 1374 fELoss->SetXTitle("#Delta/#Delta_{mip} (selected)");
1375 fELoss->SetYTitle("P(#Delta/#Delta_{mip})");
1376 fELoss->SetFillColor(Color()-2);
1377 fELoss->SetFillStyle(3003);
1378 fELoss->SetLineColor(kBlack);
1379 fELoss->SetLineStyle(2);
77f97e3f 1380 fELoss->SetLineWidth(1);
f7cfc454 1381 fELoss->SetDirectory(0);
1382
1383 fELossUsed = static_cast<TH1D*>(fELoss->Clone("elossUsed"));
1384 fELossUsed->SetTitle("#Delta/#Delta_{mip} in used strips");
1385 fELossUsed->SetFillStyle(3002);
1386 fELossUsed->SetLineStyle(1);
1387 fELossUsed->SetDirectory(0);
5ca83fee 1388
1389 fPhiBefore = new TH1D("phiBefore", "#phi distribution (before recalc)",
1390 (r == 'I' || r == 'i' ? 20 : 40), 0, 2*TMath::Pi());
1391 fPhiBefore->SetDirectory(0);
1392 fPhiBefore->SetXTitle("#phi");
1393 fPhiBefore->SetYTitle("Events");
1394 fPhiBefore->SetMarkerColor(Color());
1395 fPhiBefore->SetLineColor(Color());
1396 fPhiBefore->SetFillColor(Color());
bfab35d9 1397 fPhiBefore->SetFillStyle(3001);
5ca83fee 1398 fPhiBefore->SetMarkerStyle(20);
1399
1400 fPhiAfter = static_cast<TH1D*>(fPhiBefore->Clone("phiAfter"));
1401 fPhiAfter->SetTitle("#phi distribution (after re-calc)");
1402 fPhiAfter->SetDirectory(0);
f7cfc454 1403
7e4038b5 1404}
1405//____________________________________________________________________
1406AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
5ca83fee 1407 : AliForwardUtil::RingHistos(o),
1408 fList(o.fList),
8449e3e0 1409 // fEvsN(o.fEvsN),
1410 // fEvsM(o.fEvsM),
1411 // fEtaVsN(o.fEtaVsN),
1412 // fEtaVsM(o.fEtaVsM),
0bd4b00f 1413 fCorr(o.fCorr),
9d05ffeb 1414 fDensity(o.fDensity),
1415 fELossVsPoisson(o.fELossVsPoisson),
5ca83fee 1416 fDiffELossPoisson(o.fDiffELossPoisson),
77f97e3f
CHC		 1417 fELossVsPoissonOut(o.fELossVsPoissonOut),
1418 fDiffELossPoissonOut(o.fDiffELossPoissonOut),
1419 fOutliers(o.fOutliers),
821ffd28 1420 fPoisson(o.fPoisson),
f7cfc454 1421 fELoss(o.fELoss),
1422 fELossUsed(o.fELossUsed),
5ca83fee 1423 fMultCut(o.fMultCut),
1424 fTotal(o.fTotal),
1425 fGood(o.fGood),
1426 fPhiAcc(o.fPhiAcc),
1427 fPhiBefore(o.fPhiBefore),
1428 fPhiAfter(o.fPhiAfter)
7984e5f7 1429{
1430 //
1431 // Copy constructor
1432 //
1433 // Parameters:
1434 // o Object to copy from
1435 //
1436}
7e4038b5 1437
1438//____________________________________________________________________
1439AliFMDDensityCalculator::RingHistos&
1440AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
1441{
7984e5f7 1442 //
1443 // Assignment operator
1444 //
1445 // Parameters:
1446 // o Object to assign from
1447 //
1448 // Return:
1449 // Reference to this
1450 //
d015ecfe 1451 if (&o == this) return *this;
9d99b0dd 1452 AliForwardUtil::RingHistos::operator=(o);
7e4038b5 1453
8449e3e0 1454 // if (fEvsN) delete fEvsN;
1455 // if (fEvsM) delete fEvsM;
1456 // if (fEtaVsN) delete fEtaVsN;
1457 // if (fEtaVsM) delete fEtaVsM;
5ca83fee 1458 if (fCorr) delete fCorr;
1459 if (fDensity) delete fDensity;
1460 if (fELossVsPoisson) delete fELossVsPoisson;
1461 if (fDiffELossPoisson) delete fDiffELossPoisson;
1462 if (fTotal) delete fTotal;
1463 if (fGood) delete fGood;
1464 if (fPhiAcc) delete fPhiAcc;
1465 if (fPhiBefore) delete fPhiBefore;
1466 if (fPhiAfter) delete fPhiAfter;
1467
8449e3e0 1468 // fEvsN = static_cast<TH2D*>(o.fEvsN->Clone());
1469 // fEvsM = static_cast<TH2D*>(o.fEvsM->Clone());
1470 // fEtaVsN = static_cast<TProfile*>(o.fEtaVsN->Clone());
1471 // fEtaVsM = static_cast<TProfile*>(o.fEtaVsM->Clone());
5ca83fee 1472 fCorr = static_cast<TProfile*>(o.fCorr->Clone());
1473 fDensity = static_cast<TH2D*>(o.fDensity->Clone());
1474 fELossVsPoisson = static_cast<TH2D*>(o.fELossVsPoisson->Clone());
1475 fDiffELossPoisson = static_cast<TH1D*>(o.fDiffELossPoisson->Clone());
77f97e3f
CHC		 1476 fELossVsPoissonOut = static_cast<TH2D*>(o.fELossVsPoisson->Clone());
1477 fDiffELossPoissonOut = static_cast<TH1D*>(o.fDiffELossPoisson->Clone());
1478 fOutliers = static_cast<TH1D*>(o.fOutliers->Clone());
5ca83fee 1479 fPoisson = o.fPoisson;
1480 fELoss = static_cast<TH1D*>(o.fELoss->Clone());
1481 fELossUsed = static_cast<TH1D*>(o.fELossUsed->Clone());
1482 fTotal = static_cast<TH1D*>(o.fTotal->Clone());
1483 fGood = static_cast<TH1D*>(o.fGood->Clone());
1484 fPhiAcc = static_cast<TH2D*>(o.fPhiAcc->Clone());
1485 fPhiBefore = static_cast<TH1D*>(o.fPhiBefore->Clone());
1486 fPhiAfter = static_cast<TH1D*>(o.fPhiAfter->Clone());
7e4038b5 1487 return *this;
1488}
1489//____________________________________________________________________
1490AliFMDDensityCalculator::RingHistos::~RingHistos()
1491{
7984e5f7 1492 //
1493 // Destructor
1494 //
9d05ffeb 1495}
1496
e308a636 1497
1498//____________________________________________________________________
9d05ffeb 1499void
5934a3e3 1500AliFMDDensityCalculator::RingHistos::SetupForData(const TAxis& eAxis)
9d05ffeb 1501{
5ca83fee 1502 // Initialize
1503 // This is called on first event
e18cb8bd 1504 fPoisson.Init(-1,-1);
5ca83fee 1505 fTotal = new TH1D("total", "Total # of strips per #eta",
1506 eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax());
1507 fTotal->SetDirectory(0);
1508 fTotal->SetXTitle("#eta");
1509 fTotal->SetYTitle("# of strips");
1510 fGood = static_cast<TH1D*>(fTotal->Clone("good"));
1511 fGood->SetTitle("# of good strips per #eta");
1512 fGood->SetDirectory(0);
1513
1514 fPhiAcc = new TH2D("phiAcc", "#phi acceptance vs Ip_{z}",
1515 eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
1516 10, -10, 10);
1517 fPhiAcc->SetXTitle("#eta");
1518 fPhiAcc->SetYTitle("v_{z} [cm]");
1519 fPhiAcc->SetZTitle("#phi acceptance");
1520 fPhiAcc->SetDirectory(0);
1521
1522 if (fList) fList->Add(fPhiAcc);
7e4038b5 1523}
1524
1525//____________________________________________________________________
1526void
5934a3e3 1527AliFMDDensityCalculator::RingHistos::CreateOutputObjects(TList* dir)
7e4038b5 1528{
7984e5f7 1529 //
5ca83fee 1530 // Make output. This is called as part of SlaveBegin
7984e5f7 1531 //
1532 // Parameters:
1533 // dir Where to put it
1534 //
9d99b0dd 1535 TList* d = DefineOutputList(dir);
8449e3e0 1536 // d->Add(fEvsN);
1537 // d->Add(fEvsM);
1538 // d->Add(fEtaVsN);
1539 // d->Add(fEtaVsM);
0bd4b00f 1540 d->Add(fCorr);
7e4038b5 1541 d->Add(fDensity);
9d05ffeb 1542 d->Add(fELossVsPoisson);
77f97e3f 1543 d->Add(fELossVsPoissonOut);
5ca83fee 1544 d->Add(fDiffELossPoisson);
77f97e3f
CHC		 1545 d->Add(fDiffELossPoissonOut);
1546 d->Add(fOutliers);
821ffd28 1547 fPoisson.Output(d);
d23503ee 1548 fPoisson.GetOccupancy()->SetFillColor(Color());
1549 fPoisson.GetMean()->SetFillColor(Color());
1550 fPoisson.GetOccupancy()->SetFillColor(Color());
f7cfc454 1551 d->Add(fELoss);
1552 d->Add(fELossUsed);
5ca83fee 1553 d->Add(fPhiBefore);
1554 d->Add(fPhiAfter);
d23503ee 1555
77f97e3f
CHC		 1556 TAxis x(NStrip(), -.5, NStrip()-.5);
1557 TAxis y(NSector(), -.5, NSector()-.5);
d23503ee 1558 x.SetTitle("strip");
1559 y.SetTitle("sector");
1560 fPoisson.Define(x, y);
1561
241cca4d 1562 d->Add(AliForwardUtil::MakeParameter("cut", fMultCut));
5ca83fee 1563 fList = d;
9d99b0dd 1564}
1565
1566//____________________________________________________________________
1567void
5934a3e3 1568AliFMDDensityCalculator::RingHistos::Terminate(TList* dir, Int_t nEvents)
9d99b0dd 1569{
7984e5f7 1570 //
1571 // Scale the histograms to the total number of events
1572 //
1573 // Parameters:
1574 // dir Where the output is
1575 // nEvents Number of events
1576 //
9d99b0dd 1577 TList* l = GetOutputList(dir);
1578 if (!l) return;
1579
5934a3e3 1580 TH2D* density = static_cast<TH2D*>(GetOutputHist(l,"inclDensity"));
9d99b0dd 1581 if (density) density->Scale(1./nEvents);
5934a3e3 1582 fDensity = density;
7e4038b5 1583}
1584
1585//____________________________________________________________________
1586//
1587// EOF
1588//
1589
1590
1591