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