]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/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
Added task to do the analysis on track refs in the central
[u/mrichter/AliRoot.git] / PWG2 / 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"
13#include <TH2D.h>
0bd4b00f 14#include <TProfile.h>
5bb5d1f6 15#include <THStack.h>
0bd4b00f 16#include <TROOT.h>
17#include <iostream>
18#include <iomanip>
7e4038b5 19
20ClassImp(AliFMDDensityCalculator)
21#if 0
22; // For Emacs
23#endif
24
25//____________________________________________________________________
26AliFMDDensityCalculator::AliFMDDensityCalculator()
27 : TNamed(),
28 fRingHistos(),
0bd4b00f 29 fMultCut(0),
5bb5d1f6 30 fNXi(1),
31 fIncludeSigma(true),
dd497217 32 fSumOfWeights(0),
33 fWeightedSum(0),
ea3e5d95 34 fCorrections(0),
0bd4b00f 35 fMaxParticles(5),
9d05ffeb 36 fUsePoisson(false),
5bb5d1f6 37 fUsePhiAcceptance(false),
0bd4b00f 38 fAccI(0),
39 fAccO(0),
1174780f 40 fFMD1iMax(0),
41 fFMD2iMax(0),
42 fFMD2oMax(0),
43 fFMD3iMax(0),
44 fFMD3oMax(0),
5bb5d1f6 45 fMaxWeights(0),
46 fLowCuts(0),
b6b35c77 47 fEtaLumping(1),
48 fPhiLumping(1),
ea3e5d95 49 fDebug(0)
7984e5f7 50{
51 //
52 // Constructor
53 //
95985fc4 54
7984e5f7 55}
7e4038b5 56
57//____________________________________________________________________
58AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
59 : TNamed("fmdDensityCalculator", title),
60 fRingHistos(),
0bd4b00f 61 fMultCut(0),
5bb5d1f6 62 fNXi(1),
63 fIncludeSigma(true),
dd497217 64 fSumOfWeights(0),
65 fWeightedSum(0),
ea3e5d95 66 fCorrections(0),
0bd4b00f 67 fMaxParticles(5),
9d05ffeb 68 fUsePoisson(false),
5bb5d1f6 69 fUsePhiAcceptance(false),
0bd4b00f 70 fAccI(0),
71 fAccO(0),
1174780f 72 fFMD1iMax(0),
73 fFMD2iMax(0),
74 fFMD2oMax(0),
75 fFMD3iMax(0),
76 fFMD3oMax(0),
5bb5d1f6 77 fMaxWeights(0),
78 fLowCuts(0),
b6b35c77 79 fEtaLumping(5),
80 fPhiLumping(1),
ea3e5d95 81 fDebug(0)
7e4038b5 82{
7984e5f7 83 //
84 // Constructor
85 //
86 // Parameters:
87 // name Name of object
88 //
7e4038b5 89 fRingHistos.SetName(GetName());
e308a636 90 fRingHistos.SetOwner();
7e4038b5 91 fRingHistos.Add(new RingHistos(1, 'I'));
92 fRingHistos.Add(new RingHistos(2, 'I'));
93 fRingHistos.Add(new RingHistos(2, 'O'));
94 fRingHistos.Add(new RingHistos(3, 'I'));
95 fRingHistos.Add(new RingHistos(3, 'O'));
dd497217 96 fSumOfWeights = new TH1D("sumOfWeights", "Sum of Landau weights",
97 200, 0, 20);
0bd4b00f 98 fSumOfWeights->SetFillColor(kRed+1);
99 fSumOfWeights->SetXTitle("#sum_{i} a_{i} f_{i}(#Delta)");
dd497217 100 fWeightedSum = new TH1D("weightedSum", "Weighted sum of Landau propability",
101 200, 0, 20);
0bd4b00f 102 fWeightedSum->SetFillColor(kBlue+1);
103 fWeightedSum->SetXTitle("#sum_{i} i a_{i} f_{i}(#Delta)");
dd497217 104 fCorrections = new TH1D("corrections", "Distribution of corrections",
105 100, 0, 10);
0bd4b00f 106 fCorrections->SetFillColor(kBlue+1);
107 fCorrections->SetXTitle("correction");
dd497217 108
0bd4b00f 109 fAccI = GenerateAcceptanceCorrection('I');
110 fAccO = GenerateAcceptanceCorrection('O');
5bb5d1f6 111
112 fMaxWeights = new TH2D("maxWeights", "Maximum i of a_{i}'s to use",
113 1, 0, 1, 1, 0, 1);
114 fMaxWeights->SetXTitle("#eta");
115 fMaxWeights->SetDirectory(0);
116
117 fLowCuts = new TH2D("lowCuts", "Low cuts used", 1, 0, 1, 1, 0, 1);
118 fLowCuts->SetXTitle("#eta");
119 fLowCuts->SetDirectory(0);
120
121 for (Int_t i = 0; i < 5; i++) fMultCuts[i] = 0;
7e4038b5 122}
123
124//____________________________________________________________________
125AliFMDDensityCalculator::AliFMDDensityCalculator(const
126 AliFMDDensityCalculator& o)
127 : TNamed(o),
128 fRingHistos(),
dd497217 129 fMultCut(o.fMultCut),
5bb5d1f6 130 fNXi(o.fNXi),
131 fIncludeSigma(o.fIncludeSigma),
dd497217 132 fSumOfWeights(o.fSumOfWeights),
133 fWeightedSum(o.fWeightedSum),
ea3e5d95 134 fCorrections(o.fCorrections),
0bd4b00f 135 fMaxParticles(o.fMaxParticles),
9d05ffeb 136 fUsePoisson(o.fUsePoisson),
5bb5d1f6 137 fUsePhiAcceptance(o.fUsePhiAcceptance),
0bd4b00f 138 fAccI(o.fAccI),
139 fAccO(o.fAccO),
1174780f 140 fFMD1iMax(o.fFMD1iMax),
141 fFMD2iMax(o.fFMD2iMax),
142 fFMD2oMax(o.fFMD2oMax),
143 fFMD3iMax(o.fFMD3iMax),
144 fFMD3oMax(o.fFMD3oMax),
5bb5d1f6 145 fMaxWeights(o.fMaxWeights),
146 fLowCuts(o.fLowCuts),
b6b35c77 147 fEtaLumping(o.fEtaLumping),
148 fPhiLumping(o.fPhiLumping),
ea3e5d95 149 fDebug(o.fDebug)
7e4038b5 150{
7984e5f7 151 //
152 // Copy constructor
153 //
154 // Parameters:
155 // o Object to copy from
156 //
7e4038b5 157 TIter next(&o.fRingHistos);
158 TObject* obj = 0;
159 while ((obj = next())) fRingHistos.Add(obj);
160}
161
162//____________________________________________________________________
163AliFMDDensityCalculator::~AliFMDDensityCalculator()
164{
7984e5f7 165 //
166 // Destructor
167 //
7e4038b5 168 fRingHistos.Delete();
169}
170
171//____________________________________________________________________
172AliFMDDensityCalculator&
173AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
174{
7984e5f7 175 //
176 // Assignement operator
177 //
178 // Parameters:
179 // o Object to assign from
180 //
181 // Return:
182 // Reference to this object
183 //
ea3e5d95 184 TNamed::operator=(o);
7e4038b5 185
0bd4b00f 186 fMultCut = o.fMultCut;
5bb5d1f6 187 fNXi = o.fNXi;
188 fIncludeSigma = o.fIncludeSigma;
0bd4b00f 189 fDebug = o.fDebug;
190 fMaxParticles = o.fMaxParticles;
9d05ffeb 191 fUsePoisson = o.fUsePoisson;
5bb5d1f6 192 fUsePhiAcceptance = o.fUsePhiAcceptance;
0bd4b00f 193 fAccI = o.fAccI;
194 fAccO = o.fAccO;
1174780f 195 fFMD1iMax = o.fFMD1iMax;
196 fFMD2iMax = o.fFMD2iMax;
197 fFMD2oMax = o.fFMD2oMax;
198 fFMD3iMax = o.fFMD3iMax;
199 fFMD3oMax = o.fFMD3oMax;
5bb5d1f6 200 fMaxWeights = o.fMaxWeights;
201 fLowCuts = o.fLowCuts;
b6b35c77 202 fEtaLumping = o.fEtaLumping;
203 fPhiLumping = o.fPhiLumping;
7e4038b5 204 fRingHistos.Delete();
205 TIter next(&o.fRingHistos);
206 TObject* obj = 0;
207 while ((obj = next())) fRingHistos.Add(obj);
208
209 return *this;
210}
211
1174780f 212//____________________________________________________________________
213void
9d05ffeb 214AliFMDDensityCalculator::Init(const TAxis& axis)
1174780f 215{
216 // Intialize this sub-algorithm
217 //
218 // Parameters:
219 // etaAxis Not used
220 CacheMaxWeights();
9d05ffeb 221
222 TIter next(&fRingHistos);
223 RingHistos* o = 0;
224 while ((o = static_cast<RingHistos*>(next())))
225 o->Init(axis);
1174780f 226}
227
7e4038b5 228//____________________________________________________________________
229AliFMDDensityCalculator::RingHistos*
230AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
231{
7984e5f7 232 //
233 // Get the ring histogram container
234 //
235 // Parameters:
236 // d Detector
237 // r Ring
238 //
239 // Return:
240 // Ring histogram container
241 //
7e4038b5 242 Int_t idx = -1;
243 switch (d) {
244 case 1: idx = 0; break;
245 case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
246 case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
247 }
e308a636 248 if (idx < 0 || idx >= fRingHistos.GetEntries()) {
249 AliWarning(Form("Index %d of FMD%d%c out of range", idx, d, r));
250 return 0;
251 }
7e4038b5 252
253 return static_cast<RingHistos*>(fRingHistos.At(idx));
254}
5bb5d1f6 255
256//____________________________________________________________________
257void
258AliFMDDensityCalculator::SetMultCuts(Double_t fmd1i,
259 Double_t fmd2i,
260 Double_t fmd2o,
261 Double_t fmd3i,
262 Double_t fmd3o)
263{
264 fMultCuts[0] = fmd1i;
265 fMultCuts[1] = fmd2i;
266 fMultCuts[2] = fmd2o;
267 fMultCuts[3] = fmd3i;
268 fMultCuts[4] = fmd3o;
269 fMultCut = (fmd1i+fmd2i+fmd2o+fmd3i+fmd3o) / 5;
270}
271
0bd4b00f 272//____________________________________________________________________
273Double_t
5bb5d1f6 274AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Int_t ieta,
275 Bool_t errors) const
0bd4b00f 276{
7984e5f7 277 //
278 // Get the multiplicity cut. If the user has set fMultCut (via
279 // SetMultCut) then that value is used. If not, then the lower
280 // value of the fit range for the enery loss fits is returned.
281 //
282 // Return:
283 // Lower cut on multiplicity
284 //
5bb5d1f6 285 Int_t idx = (d == 1 ? 0 : 2*(d - 2) + 1 + ((r=='I' || r=='i') ? 0 : 1));
286 if (fMultCuts[idx] > 0) return fMultCuts[idx];
0bd4b00f 287 if (fMultCut > 0) return fMultCut;
288
289 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
290 AliFMDCorrELossFit* fits = fcm.GetELossFit();
5bb5d1f6 291 if (fNXi < 0) return fits->GetLowCut();
292
293 return fits->GetLowerBound(d, r, ieta, fNXi, errors, fIncludeSigma);
294}
295
296//____________________________________________________________________
297Double_t
298AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Double_t eta,
299 Bool_t /*errors*/) const
300{
301 //
302 // Get the multiplicity cut. If the user has set fMultCut (via
303 // SetMultCut) then that value is used. If not, then the lower
304 // value of the fit range for the enery loss fits is returned.
305 //
306 // Return:
307 // Lower cut on multiplicity
308 //
309 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
310 AliFMDCorrELossFit* fits = fcm.GetELossFit();
311 Int_t iEta = fits->FindEtaBin(eta);
312
313 return GetMultCut(d, r, iEta);
0bd4b00f 314}
315
7e4038b5 316//____________________________________________________________________
317Bool_t
318AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd,
319 AliForwardUtil::Histos& hists,
0bd4b00f 320 UShort_t vtxbin,
7e4038b5 321 Bool_t lowFlux)
322{
7984e5f7 323 //
324 // Do the calculations
325 //
326 // Parameters:
327 // fmd AliESDFMD object (possibly) corrected for sharing
328 // hists Histogram cache
329 // vtxBin Vertex bin
330 // lowFlux Low flux flag.
331 //
332 // Return:
333 // true on successs
334 //
7e4038b5 335 for (UShort_t d=1; d<=3; d++) {
336 UShort_t nr = (d == 1 ? 1 : 2);
337 for (UShort_t q=0; q<nr; q++) {
338 Char_t r = (q == 0 ? 'I' : 'O');
339 UShort_t ns= (q == 0 ? 20 : 40);
340 UShort_t nt= (q == 0 ? 512 : 256);
341 TH2D* h = hists.Get(d,r);
342 RingHistos* rh= GetRingHistos(d,r);
e308a636 343 if (!rh) {
344 AliError(Form("No ring histogram found for FMD%d%c", d, r));
345 fRingHistos.ls();
346 return false;
347 }
b6b35c77 348 rh->ResetPoissonHistos(h, fEtaLumping, fPhiLumping);
9fde7142 349
7e4038b5 350 for (UShort_t s=0; s<ns; s++) {
351 for (UShort_t t=0; t<nt; t++) {
5bb5d1f6 352 Float_t mult = fmd.Multiplicity(d,r,s,t);
353 Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
354 Float_t eta = fmd.Eta(d,r,s,t);
355 Double_t cut = 1024;
356 if (eta != AliESDFMD::kInvalidEta) cut = GetMultCut(d, r, eta,false);
9d05ffeb 357 rh->fTotalStrips->Fill(eta, phi);
9fde7142 358
5bb5d1f6 359 if (mult == AliESDFMD::kInvalidMult || mult > 20) {
33b078ef 360 rh->fEmptyStrips->Fill(eta,phi);
5bb5d1f6 361 rh->fEvsM->Fill(mult,0);
33b078ef 362 continue;
363 }
5bb5d1f6 364
365 Double_t n = 0;
366 if (cut > 0 && mult > cut)
367 n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
9fde7142 368
7e4038b5 369 rh->fEvsN->Fill(mult,n);
0bd4b00f 370 rh->fEtaVsN->Fill(eta, n);
9fde7142 371
5bb5d1f6 372 Double_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
dd497217 373 fCorrections->Fill(c);
7e4038b5 374 if (c > 0) n /= c;
375 rh->fEvsM->Fill(mult,n);
0bd4b00f 376 rh->fEtaVsM->Fill(eta, n);
377 rh->fCorr->Fill(eta, c);
9fde7142 378
b6b35c77 379 if (n > 0.9 && c > 0) rh->fBasicHits->Fill(eta,phi, 1./c);
380 else rh->fEmptyStrips->Fill(eta,phi);
9fde7142 381
7e4038b5 382 h->Fill(eta,phi,n);
5bb5d1f6 383 if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n);
7e4038b5 384 } // for t
385 } // for s
b6b35c77 386
387
388 // --- Loop over poisson histograms ----------------------------
9d05ffeb 389 for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
390 for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) {
391 Double_t eLossV = h->GetBinContent(ieta, iphi);
5bb5d1f6 392 Double_t eta = h->GetXaxis()->GetBinCenter(ieta);
393 Double_t phi = h->GetYaxis()->GetBinCenter(iphi);
b6b35c77 394 Int_t jeta = rh->fEmptyStrips->GetXaxis()->FindBin(eta);
395 Int_t jphi = rh->fEmptyStrips->GetYaxis()->FindBin(phi);
396 Double_t empty = rh->fEmptyStrips->GetBinContent(jeta, jphi);
397 Double_t total = rh->fTotalStrips->GetBinContent(jeta, jphi);
9fde7142 398 Double_t hits = rh->fBasicHits->GetBinContent(ieta,iphi);
5bb5d1f6 399
b6b35c77 400 // Mean in region of interest
401 Double_t poissonM = (total <= 0 || empty <= 0 ? 0 :
402 -TMath::Log(empty / total));
5bb5d1f6 403
404 // Note, that given filled=total-empty, and
405 //
406 // m = -log(empty/total)
407 // = -log(1 - filled/total)
408 //
409 // v = m / (1 - exp(-m))
410 // = -total/filled * (log(total-filled)-log(total))
411 // = -total / (total-empty) * log(empty/total)
412 // = total (log(total)-log(empty)) / (total-empty)
413 //
33b078ef 414 Double_t poissonV = hits;
b6b35c77 415 if(poissonM > 0)
416 // Correct for counting statistics and weight by counts
5bb5d1f6 417 poissonV *= poissonM / (1 - TMath::Exp(-1*poissonM));
33b078ef 418 Double_t poissonE = TMath::Sqrt(hits);
12fffad7 419 if(poissonV > 0) poissonE = TMath::Sqrt(poissonV);
7b03929c 420
9d05ffeb 421 rh->fELossVsPoisson->Fill(eLossV, poissonV);
422 rh->fEmptyVsTotal->Fill(total, empty);
423 if (fUsePoisson) {
424 h->SetBinContent(ieta,iphi,poissonV);
425 h->SetBinError(ieta,iphi,poissonE);
5bb5d1f6 426 rh->fDensity->Fill(eta, phi, poissonV);
9d05ffeb 427 }
428 }
429 }
430
7e4038b5 431 } // for q
432 } // for d
433
434 return kTRUE;
435}
436
1174780f 437//_____________________________________________________________________
438Int_t
fb3430ac 439AliFMDDensityCalculator::FindMaxWeight(const AliFMDCorrELossFit* cor,
1174780f 440 UShort_t d, Char_t r, Int_t iEta) const
441{
fb3430ac 442 //
443 // Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
444 //
445 // Parameters:
446 // cor Correction
447 // d Detector
448 // r Ring
449 // iEta Eta bin
450 //
1174780f 451 AliFMDCorrELossFit::ELossFit* fit = cor->GetFit(d,r,iEta);
452 if (!fit) {
453 // AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
454 return -1;
455 }
1f480471 456 return TMath::Min(Int_t(fMaxParticles), fit->FindMaxWeight());
1174780f 457}
458
459//_____________________________________________________________________
460void
461AliFMDDensityCalculator::CacheMaxWeights()
462{
fb3430ac 463 //
464 // Find the max weights and cache them
465 //
1174780f 466 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
467 AliFMDCorrELossFit* cor = fcm.GetELossFit();
468 const TAxis& eta = cor->GetEtaAxis();
469
470 Int_t nEta = eta.GetNbins();
471 fFMD1iMax.Set(nEta);
472 fFMD2iMax.Set(nEta);
473 fFMD2oMax.Set(nEta);
474 fFMD3iMax.Set(nEta);
475 fFMD3oMax.Set(nEta);
476
5bb5d1f6 477 fMaxWeights->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
478 fMaxWeights->GetYaxis()->SetBinLabel(1, "FMD1i");
479 fMaxWeights->GetYaxis()->SetBinLabel(2, "FMD2i");
480 fMaxWeights->GetYaxis()->SetBinLabel(3, "FMD2o");
481 fMaxWeights->GetYaxis()->SetBinLabel(4, "FMD3i");
482 fMaxWeights->GetYaxis()->SetBinLabel(5, "FMD3o");
483
484 AliInfo(Form("Get eta axis with %d bins from %f to %f",
485 nEta, eta.GetXmin(), eta.GetXmax()));
486 fLowCuts->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
487 fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
488 fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
489 fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
490 fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
491 fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
492
1174780f 493 for (Int_t i = 0; i < nEta; i++) {
5bb5d1f6 494 Double_t w[5];
495 w[0] = fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1);
496 w[1] = fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1);
497 w[2] = fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1);
498 w[3] = fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1);
499 w[4] = fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1);
500 Double_t l[5];
501 l[0] = GetMultCut(1, 'I', i+1, false);
502 l[1] = GetMultCut(2, 'I', i+1, false);
503 l[2] = GetMultCut(2, 'O', i+1, false);
504 l[3] = GetMultCut(3, 'I', i+1, false);
505 l[4] = GetMultCut(3, 'O', i+1, false);
506 for (Int_t j = 0; j < 5; j++) {
507 if (w[j] > 0) fMaxWeights->SetBinContent(i+1, j+1, w[j]);
508 if (l[j] > 0) fLowCuts->SetBinContent(i+1, j+1, l[j]);
509 }
1174780f 510 }
511}
512
513//_____________________________________________________________________
514Int_t
515AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const
516{
fb3430ac 517 //
518 // Find the (cached) maximum weight for FMD<i>dr</i> in
519 // @f$\eta@f$ bin @a iEta
520 //
521 // Parameters:
522 // d Detector
523 // r Ring
524 // iEta Eta bin
525 //
526 // Return:
527 // max weight or <= 0 in case of problems
528 //
1174780f 529 if (iEta < 0) return -1;
530
531 const TArrayI* max = 0;
532 switch (d) {
533 case 1: max = &fFMD1iMax; break;
534 case 2: max = (r == 'I' || r == 'i' ? &fFMD2iMax : &fFMD2oMax); break;
535 case 3: max = (r == 'I' || r == 'i' ? &fFMD3iMax : &fFMD3oMax); break;
536 }
537 if (!max) {
538 AliWarning(Form("No array for FMD%d%c", d, r));
539 return -1;
540 }
541
542 if (iEta >= max->fN) {
543 AliWarning(Form("Eta bin %3d out of bounds [0,%d]",
544 iEta, max->fN-1));
545 return -1;
546 }
547
548 AliDebug(30,Form("Max weight for FMD%d%c eta bin %3d: %d", d, r, iEta,
549 max->At(iEta)));
550 return max->At(iEta);
551}
552
553//_____________________________________________________________________
554Int_t
555AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const
556{
fb3430ac 557 //
558 // Find the (cached) maximum weight for FMD<i>dr</i> iat
559 // @f$\eta@f$
560 //
561 // Parameters:
562 // d Detector
563 // r Ring
564 // eta Eta bin
565 //
566 // Return:
567 // max weight or <= 0 in case of problems
568 //
1174780f 569 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
570 Int_t iEta = fcm.GetELossFit()->FindEtaBin(eta) -1;
571
572 return GetMaxWeight(d, r, iEta);
573}
574
7e4038b5 575//_____________________________________________________________________
576Float_t
577AliFMDDensityCalculator::NParticles(Float_t mult,
578 UShort_t d,
579 Char_t r,
580 UShort_t /*s*/,
581 UShort_t /*t*/,
0bd4b00f 582 UShort_t /*v*/,
7e4038b5 583 Float_t eta,
584 Bool_t lowFlux) const
585{
7984e5f7 586 //
587 // Get the number of particles corresponding to the signal mult
588 //
589 // Parameters:
590 // mult Signal
591 // d Detector
592 // r Ring
593 // s Sector
594 // t Strip (not used)
595 // v Vertex bin
596 // eta Pseudo-rapidity
597 // lowFlux Low-flux flag
598 //
599 // Return:
600 // The number of particles
601 //
5bb5d1f6 602 // if (mult <= GetMultCut()) return 0;
7e4038b5 603 if (lowFlux) return 1;
0bd4b00f 604
0bd4b00f 605 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
606 AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
607 if (!fit) {
608 AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
609 return 0;
610 }
611
1174780f 612 Int_t m = GetMaxWeight(d,r,eta); // fit->FindMaxWeight();
0bd4b00f 613 if (m < 1) {
614 AliWarning(Form("No good fits for FMD%d%c at eta=%f", d, r, eta));
615 return 0;
616 }
617 UShort_t n = TMath::Min(fMaxParticles, UShort_t(m));
618 Double_t ret = fit->EvaluateWeighted(mult, n);
619
620 if (fDebug > 10) {
621 AliInfo(Form("FMD%d%c, eta=%7.4f, %8.5f -> %8.5f", d, r, eta, mult, ret));
622 }
1174780f 623
0bd4b00f 624 fWeightedSum->Fill(ret);
625 fSumOfWeights->Fill(ret);
7e4038b5 626
0bd4b00f 627 return ret;
7e4038b5 628}
629
630//_____________________________________________________________________
631Float_t
0bd4b00f 632AliFMDDensityCalculator::Correction(UShort_t d,
633 Char_t r,
634 UShort_t /*s*/,
635 UShort_t t,
636 UShort_t /*v*/,
637 Float_t eta,
638 Bool_t lowFlux) const
7e4038b5 639{
7984e5f7 640 //
641 // Get the inverse correction factor. This consist of
642 //
643 // - acceptance correction (corners of sensors)
644 // - double hit correction (for low-flux events)
645 // - dead strip correction
646 //
647 // Parameters:
648 // d Detector
649 // r Ring
650 // s Sector
651 // t Strip (not used)
652 // v Vertex bin
653 // eta Pseudo-rapidity
654 // lowFlux Low-flux flag
655 //
656 // Return:
657 //
658 //
0bd4b00f 659 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
7e4038b5 660
5bb5d1f6 661 Float_t correction = 1;
662 if (fUsePhiAcceptance) correction = AcceptanceCorrection(r,t);
7e4038b5 663 if (lowFlux) {
0bd4b00f 664 TH1D* dblHitCor = 0;
665 if (fcm.GetDoubleHit())
666 dblHitCor = fcm.GetDoubleHit()->GetCorrection(d,r);
667
7e4038b5 668 if (dblHitCor) {
0bd4b00f 669 Double_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
670 if (dblC > 0) correction *= dblC;
7e4038b5 671 }
0bd4b00f 672 else {
7e4038b5 673 AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
0bd4b00f 674 }
7e4038b5 675 }
7e4038b5 676 return correction;
677}
678
0bd4b00f 679//_____________________________________________________________________
680TH1D*
681AliFMDDensityCalculator::GenerateAcceptanceCorrection(Char_t r) const
682{
7984e5f7 683 //
684 // Generate the acceptance corrections
685 //
686 // Parameters:
687 // r Ring to generate for
688 //
689 // Return:
690 // Newly allocated histogram of acceptance corrections
691 //
0bd4b00f 692 const Double_t ic1[] = { 4.9895, 15.3560 };
693 const Double_t ic2[] = { 1.8007, 17.2000 };
694 const Double_t oc1[] = { 4.2231, 26.6638 };
695 const Double_t oc2[] = { 1.8357, 27.9500 };
696 const Double_t* c1 = (r == 'I' || r == 'i' ? ic1 : oc1);
697 const Double_t* c2 = (r == 'I' || r == 'i' ? ic2 : oc2);
698 Double_t minR = (r == 'I' || r == 'i' ? 4.5213 : 15.4);
699 Double_t maxR = (r == 'I' || r == 'i' ? 17.2 : 28.0);
700 Int_t nStrips = (r == 'I' || r == 'i' ? 512 : 256);
701 Int_t nSec = (r == 'I' || r == 'i' ? 20 : 40);
702 Float_t basearc = 2 * TMath::Pi() / nSec;
703 Double_t rad = maxR - minR;
704 Float_t segment = rad / nStrips;
705 Float_t cr = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
706
707 // Numbers used to find end-point of strip.
708 // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
709 Float_t D = c1[0] * c2[1] - c1[1] * c2[0];
710 Float_t dx = c2[0] - c1[0];
711 Float_t dy = c2[1] - c1[1];
712 Float_t dr = TMath::Sqrt(dx*dx+dy*dy);
713
714 TH1D* ret = new TH1D(Form("acc%c", r),
715 Form("Acceptance correction for FMDx%c", r),
716 nStrips, -.5, nStrips-.5);
717 ret->SetXTitle("Strip");
718 ret->SetYTitle("#varphi acceptance");
719 ret->SetDirectory(0);
720 ret->SetFillColor(r == 'I' || r == 'i' ? kRed+1 : kBlue+1);
721 ret->SetFillStyle(3001);
722
723 for (Int_t t = 0; t < nStrips; t++) {
724 Float_t radius = minR + t * segment;
725
726 // If the radius of the strip is smaller than the radius corresponding
727 // to the first corner we have a full strip length
728 if (radius <= cr) {
729 ret->SetBinContent(t+1, 1);
730 continue;
731 }
732
733 // Next, we should find the end-point of the strip - that is,
734 // the coordinates where the line from c1 to c2 intersects a circle
735 // with radius given by the strip.
736 // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
737 // Calculate the determinant
738 Float_t det = radius * radius * dr * dr - D*D;
739
740 if (det <= 0) {
741 // <0 means No intersection
742 // =0 means Exactly tangent
743 ret->SetBinContent(t+1, 1);
744 continue;
745 }
746
747 // Calculate end-point and the corresponding opening angle
748 Float_t x = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
749 Float_t y = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
750 Float_t th = TMath::ATan2(x, y);
751
752 ret->SetBinContent(t+1, th / basearc);
753 }
754 return ret;
755}
7e4038b5 756
757//_____________________________________________________________________
758Float_t
759AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
760{
7984e5f7 761 //
762 // Get the acceptance correction for strip @a t in an ring of type @a r
763 //
764 // Parameters:
765 // r Ring type ('I' or 'O')
766 // t Strip number
767 //
768 // Return:
769 // Inverse acceptance correction
770 //
0bd4b00f 771 TH1D* acc = (r == 'I' || r == 'i' ? fAccI : fAccO);
772 return acc->GetBinContent(t+1);
7e4038b5 773}
774
775//____________________________________________________________________
776void
fb3430ac 777AliFMDDensityCalculator::ScaleHistograms(const TList* dir, Int_t nEvents)
7e4038b5 778{
7984e5f7 779 //
780 // Scale the histograms to the total number of events
781 //
782 // Parameters:
783 // dir where to put the output
784 // nEvents Number of events
785 //
7e4038b5 786 if (nEvents <= 0) return;
9d99b0dd 787 TList* d = static_cast<TList*>(dir->FindObject(GetName()));
788 if (!d) return;
7e4038b5 789
790 TIter next(&fRingHistos);
791 RingHistos* o = 0;
5bb5d1f6 792 THStack* sums = new THStack("sums", "sums of ring signals");
793 while ((o = static_cast<RingHistos*>(next()))) {
9d99b0dd 794 o->ScaleHistograms(d, nEvents);
5bb5d1f6 795 TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
796 o->fDensity->GetNbinsY(),"e");
797 sum->Scale(1., "width");
798 sum->SetTitle(o->GetName());
799 sum->SetDirectory(0);
800 sum->SetYTitle("#sum N_{ch,incl}");
801 sums->Add(sum);
802 }
803 d->Add(sums);
7e4038b5 804}
805
806//____________________________________________________________________
807void
9d99b0dd 808AliFMDDensityCalculator::DefineOutput(TList* dir)
7e4038b5 809{
7984e5f7 810 //
811 // Output diagnostic histograms to directory
812 //
813 // Parameters:
814 // dir List to write in
815 //
7e4038b5 816 TList* d = new TList;
5bb5d1f6 817 d->SetOwner();
7e4038b5 818 d->SetName(GetName());
819 dir->Add(d);
dd497217 820 d->Add(fWeightedSum);
821 d->Add(fSumOfWeights);
822 d->Add(fCorrections);
0bd4b00f 823 d->Add(fAccI);
824 d->Add(fAccO);
5bb5d1f6 825 d->Add(fMaxWeights);
826 d->Add(fLowCuts);
0bd4b00f 827
7e4038b5 828 TIter next(&fRingHistos);
829 RingHistos* o = 0;
830 while ((o = static_cast<RingHistos*>(next()))) {
831 o->Output(d);
832 }
833}
0bd4b00f 834//____________________________________________________________________
835void
1174780f 836AliFMDDensityCalculator::Print(Option_t* option) const
0bd4b00f 837{
7984e5f7 838 //
839 // Print information
840 //
841 // Parameters:
842 // option Not used
843 //
1174780f 844 char ind[gROOT->GetDirLevel()+3];
0bd4b00f 845 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
846 ind[gROOT->GetDirLevel()] = '\0';
1f480471 847 std::cout << ind << ClassName() << ": " << GetName() << '\n'
5bb5d1f6 848 << std::boolalpha
0bd4b00f 849 << ind << " Multiplicity cut: " << fMultCut << '\n'
5bb5d1f6 850 << ind << " # of (xi+sigma) factor: " << fNXi << '\n'
851 << ind << " Include sigma in cut: " << fIncludeSigma << '\n'
852 << ind << " Low cut method: "
853 << (fMultCut > 0 ? "fixed" :
854 (fNXi >= 0 ? "xi+sigma" : "fit range")) << '\n'
1f480471 855 << ind << " Max(particles): " << fMaxParticles << '\n'
5bb5d1f6 856 << ind << " Poisson method: " << fUsePoisson << '\n'
857 << ind << " Use phi acceptance: " << fUsePhiAcceptance << '\n'
1f480471 858 << ind << " Eta lumping: " << fEtaLumping << '\n'
859 << ind << " Phi lumping: " << fPhiLumping << '\n'
5bb5d1f6 860 << std::noboolalpha
1f480471 861 << std::flush;
1174780f 862 TString opt(option);
863 opt.ToLower();
864 if (opt.Contains("nomax")) return;
865
866 std::cout << ind << " Max weights:\n";
867
868 for (UShort_t d=1; d<=3; d++) {
869 UShort_t nr = (d == 1 ? 1 : 2);
870 for (UShort_t q=0; q<nr; q++) {
871 ind[gROOT->GetDirLevel()] = ' ';
872 ind[gROOT->GetDirLevel()+1] = '\0';
873 Char_t r = (q == 0 ? 'I' : 'O');
874 std::cout << ind << " FMD" << d << r << ":";
875 ind[gROOT->GetDirLevel()+1] = ' ';
876 ind[gROOT->GetDirLevel()+2] = '\0';
877
878 const TArrayI& a = (d == 1 ? fFMD1iMax :
879 (d == 2 ? (r == 'I' ? fFMD2iMax : fFMD2oMax) :
880 (r == 'I' ? fFMD3iMax : fFMD3oMax)));
881 Int_t j = 0;
882 for (Int_t i = 0; i < a.fN; i++) {
883 if (a.fArray[i] < 1) continue;
884 if (j % 6 == 0) std::cout << "\n " << ind;
885 j++;
886 std::cout << " " << std::setw(3) << i << ": " << a.fArray[i];
887 }
888 std::cout << std::endl;
889 }
890 }
0bd4b00f 891}
7e4038b5 892
893//====================================================================
894AliFMDDensityCalculator::RingHistos::RingHistos()
9d99b0dd 895 : AliForwardUtil::RingHistos(),
7e4038b5 896 fEvsN(0),
897 fEvsM(0),
0bd4b00f 898 fEtaVsN(0),
899 fEtaVsM(0),
900 fCorr(0),
9d05ffeb 901 fDensity(0),
902 fELossVsPoisson(0),
903 fTotalStrips(0),
904 fEmptyStrips(0),
9fde7142 905 fBasicHits(0),
9d05ffeb 906 fEmptyVsTotal(0)
7984e5f7 907{
908 //
909 // Default CTOR
910 //
911}
7e4038b5 912
913//____________________________________________________________________
914AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
9d99b0dd 915 : AliForwardUtil::RingHistos(d,r),
7e4038b5 916 fEvsN(0),
917 fEvsM(0),
0bd4b00f 918 fEtaVsN(0),
919 fEtaVsM(0),
920 fCorr(0),
9d05ffeb 921 fDensity(0),
922 fELossVsPoisson(0),
923 fTotalStrips(0),
924 fEmptyStrips(0),
9fde7142 925 fBasicHits(0),
9d05ffeb 926 fEmptyVsTotal(0)
7e4038b5 927{
7984e5f7 928 //
929 // Constructor
930 //
931 // Parameters:
932 // d detector
933 // r ring
934 //
5bb5d1f6 935 fEvsN = new TH2D("elossVsNnocorr",
936 "#Delta E/#Delta E_{mip} vs uncorrected inclusive N_{ch}",
937 250, -.5, 24.5, 250, -.5, 24.5);
7e4038b5 938 fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
939 fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
940 fEvsN->Sumw2();
941 fEvsN->SetDirectory(0);
5bb5d1f6 942
943 fEvsM = static_cast<TH2D*>(fEvsN->Clone("elossVsNcorr"));
944 fEvsM->SetTitle("#Delta E/#Delta E_{mip} vs corrected inclusive N_{ch}");
7e4038b5 945 fEvsM->SetDirectory(0);
946
5bb5d1f6 947 fEtaVsN = new TProfile("etaVsNnocorr",
948 "Average inclusive N_{ch} vs #eta (uncorrected)",
949 200, -4, 6);
0bd4b00f 950 fEtaVsN->SetXTitle("#eta");
951 fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
952 fEtaVsN->SetDirectory(0);
953 fEtaVsN->SetLineColor(Color());
954 fEtaVsN->SetFillColor(Color());
5bb5d1f6 955
956 fEtaVsM = static_cast<TProfile*>(fEtaVsN->Clone("etaVsNcorr"));
957 fEtaVsM->SetTitle("Average inclusive N_{ch} vs #eta (corrected)");
0bd4b00f 958 fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
959 fEtaVsM->SetDirectory(0);
0bd4b00f 960
961
5bb5d1f6 962 fCorr = new TProfile("corr", "Average correction", 200, -4, 6);
0bd4b00f 963 fCorr->SetXTitle("#eta");
964 fCorr->SetYTitle("#LT correction#GT");
965 fCorr->SetDirectory(0);
966 fCorr->SetLineColor(Color());
967 fCorr->SetFillColor(Color());
968
5bb5d1f6 969 fDensity = new TH2D("inclDensity", "Inclusive N_{ch} density",
7e4038b5 970 200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
971 0, 2*TMath::Pi());
972 fDensity->SetDirectory(0);
5bb5d1f6 973 fDensity->Sumw2();
974 fDensity->SetMarkerColor(Color());
7e4038b5 975 fDensity->SetXTitle("#eta");
976 fDensity->SetYTitle("#phi [radians]");
977 fDensity->SetZTitle("Inclusive N_{ch} density");
9d05ffeb 978
5bb5d1f6 979 fELossVsPoisson = new TH2D("elossVsPoisson",
980 "N_{ch} from energy loss vs from Poission",
981 100, 0, 20, 100, 0, 20);
9d05ffeb 982 fELossVsPoisson->SetDirectory(0);
983 fELossVsPoisson->SetXTitle("N_{ch} from #DeltaE");
984 fELossVsPoisson->SetYTitle("N_{ch} from Poisson");
985 fELossVsPoisson->SetZTitle("Correlation");
986
5bb5d1f6 987 fEmptyVsTotal = new TH2D("emptyVsTotal",
988 "# of empty strips vs. total # strips",
989 21, -.5, 20.5, 21, -0.5, 20.5);
9d05ffeb 990 fEmptyVsTotal->SetDirectory(0);
991 fEmptyVsTotal->SetXTitle("Total # strips");
992 fEmptyVsTotal->SetYTitle("Empty # strips");
e308a636 993 fEmptyVsTotal->SetZTitle("Correlation");
7e4038b5 994}
995//____________________________________________________________________
996AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
9d99b0dd 997 : AliForwardUtil::RingHistos(o),
7e4038b5 998 fEvsN(o.fEvsN),
999 fEvsM(o.fEvsM),
0bd4b00f 1000 fEtaVsN(o.fEtaVsN),
1001 fEtaVsM(o.fEtaVsM),
1002 fCorr(o.fCorr),
9d05ffeb 1003 fDensity(o.fDensity),
1004 fELossVsPoisson(o.fELossVsPoisson),
1005 fTotalStrips(o.fTotalStrips),
1006 fEmptyStrips(o.fEmptyStrips),
9fde7142 1007 fBasicHits(o.fBasicHits),
9d05ffeb 1008 fEmptyVsTotal(o.fEmptyVsTotal)
7984e5f7 1009{
1010 //
1011 // Copy constructor
1012 //
1013 // Parameters:
1014 // o Object to copy from
1015 //
1016}
7e4038b5 1017
1018//____________________________________________________________________
1019AliFMDDensityCalculator::RingHistos&
1020AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
1021{
7984e5f7 1022 //
1023 // Assignment operator
1024 //
1025 // Parameters:
1026 // o Object to assign from
1027 //
1028 // Return:
1029 // Reference to this
1030 //
9d99b0dd 1031 AliForwardUtil::RingHistos::operator=(o);
7e4038b5 1032
9d05ffeb 1033 if (fEvsN) delete fEvsN;
1034 if (fEvsM) delete fEvsM;
1035 if (fEtaVsN) delete fEtaVsN;
1036 if (fEtaVsM) delete fEtaVsM;
1037 if (fCorr) delete fCorr;
1038 if (fDensity) delete fDensity;
1039 if (fELossVsPoisson) delete fELossVsPoisson;
1040 if (fTotalStrips) delete fTotalStrips;
1041 if (fEmptyStrips) delete fEmptyStrips;
7e4038b5 1042
9d05ffeb 1043 fEvsN = static_cast<TH2D*>(o.fEvsN->Clone());
1044 fEvsM = static_cast<TH2D*>(o.fEvsM->Clone());
1045 fEtaVsN = static_cast<TProfile*>(o.fEtaVsN->Clone());
1046 fEtaVsM = static_cast<TProfile*>(o.fEtaVsM->Clone());
1047 fCorr = static_cast<TProfile*>(o.fCorr->Clone());
1048 fDensity = static_cast<TH2D*>(o.fDensity->Clone());
1049 fELossVsPoisson = static_cast<TH2D*>(o.fELossVsPoisson);
1050 fTotalStrips = static_cast<TH2D*>(o.fTotalStrips);
1051 fEmptyStrips = static_cast<TH2D*>(o.fEmptyStrips);
7e4038b5 1052
1053 return *this;
1054}
1055//____________________________________________________________________
1056AliFMDDensityCalculator::RingHistos::~RingHistos()
1057{
7984e5f7 1058 //
1059 // Destructor
1060 //
9d05ffeb 1061 if (fEvsN) delete fEvsN;
1062 if (fEvsM) delete fEvsM;
1063 if (fEtaVsN) delete fEtaVsN;
1064 if (fEtaVsM) delete fEtaVsM;
1065 if (fCorr) delete fCorr;
1066 if (fDensity) delete fDensity;
1067 if (fELossVsPoisson) delete fELossVsPoisson;
1068 if (fTotalStrips) delete fTotalStrips;
1069 if (fEmptyStrips) delete fEmptyStrips;
1070}
1071
e308a636 1072//____________________________________________________________________
1073void
b6b35c77 1074AliFMDDensityCalculator::RingHistos::ResetPoissonHistos(const TH2D* h,
1075 Int_t etaLumping,
1076 Int_t phiLumping)
e308a636 1077{
1078 if (fTotalStrips) {
1079 fTotalStrips->Reset();
1080 fEmptyStrips->Reset();
1081 fBasicHits->Reset();
1082 return;
1083 }
1084 //Inserted by HHD
1085
b6b35c77 1086 Int_t nEta = h->GetNbinsX() / etaLumping;
1087 Int_t nEtaF = h->GetNbinsX();
1088 Double_t etaMin = h->GetXaxis()->GetXmin();
1089 Double_t etaMax = h->GetXaxis()->GetXmax();
1090 Int_t nPhi = h->GetNbinsY() / phiLumping;
1091 Int_t nPhiF = h->GetNbinsY();
1092 Double_t phiMin = h->GetYaxis()->GetXmin();
1093 Double_t phiMax = h->GetYaxis()->GetXmax();
1094
e308a636 1095 fTotalStrips = new TH2D(Form("total%s", fName.Data()),
1096 Form("Total number of strips in %s", fName.Data()),
b6b35c77 1097 nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
e308a636 1098 fEmptyStrips = new TH2D(Form("empty%s", fName.Data()),
1099 Form("Empty number of strips in %s", fName.Data()),
b6b35c77 1100 nEta, etaMin, etaMax, nPhi, phiMin, phiMax);
e308a636 1101 fBasicHits = new TH2D(Form("basichits%s", fName.Data()),
1102 Form("Basic number of hits in %s", fName.Data()),
b6b35c77 1103 nEtaF, etaMin, etaMax, nPhiF, phiMin, phiMax);
e308a636 1104
1105 fTotalStrips->SetDirectory(0);
1106 fEmptyStrips->SetDirectory(0);
1107 fBasicHits->SetDirectory(0);
1108 fTotalStrips->SetXTitle("#eta");
1109 fEmptyStrips->SetXTitle("#eta");
1110 fBasicHits->SetXTitle("#eta");
1111 fTotalStrips->SetYTitle("#varphi [radians]");
1112 fEmptyStrips->SetYTitle("#varphi [radians]");
1113 fBasicHits->SetYTitle("#varphi [radians]");
1114 fTotalStrips->Sumw2();
1115 fEmptyStrips->Sumw2();
1116 fBasicHits->Sumw2();
1117}
1118
1119//____________________________________________________________________
9d05ffeb 1120void
9fde7142 1121AliFMDDensityCalculator::RingHistos::Init(const TAxis& /*eAxis*/)
9d05ffeb 1122{
7e4038b5 1123}
1124
1125//____________________________________________________________________
1126void
1127AliFMDDensityCalculator::RingHistos::Output(TList* dir)
1128{
7984e5f7 1129 //
1130 // Make output
1131 //
1132 // Parameters:
1133 // dir Where to put it
1134 //
9d99b0dd 1135 TList* d = DefineOutputList(dir);
7e4038b5 1136 d->Add(fEvsN);
1137 d->Add(fEvsM);
0bd4b00f 1138 d->Add(fEtaVsN);
1139 d->Add(fEtaVsM);
1140 d->Add(fCorr);
7e4038b5 1141 d->Add(fDensity);
9d05ffeb 1142 d->Add(fELossVsPoisson);
1143 d->Add(fEmptyVsTotal);
9fde7142 1144 d->Add(fTotalStrips);
1145 d->Add(fEmptyStrips);
1146 d->Add(fBasicHits);
1147
1148
9d99b0dd 1149}
1150
1151//____________________________________________________________________
1152void
1153AliFMDDensityCalculator::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
1154{
7984e5f7 1155 //
1156 // Scale the histograms to the total number of events
1157 //
1158 // Parameters:
1159 // dir Where the output is
1160 // nEvents Number of events
1161 //
9d99b0dd 1162 TList* l = GetOutputList(dir);
1163 if (!l) return;
1164
5bb5d1f6 1165 TH1* density = GetOutputHist(l,"inclDensity");
9d99b0dd 1166 if (density) density->Scale(1./nEvents);
7e4038b5 1167}
1168
1169//____________________________________________________________________
1170//
1171// EOF
1172//
1173
1174
1175