]> 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
Clean-up of code to conform with the style used
[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"
13#include <TH2D.h>
0bd4b00f 14#include <TProfile.h>
5bb5d1f6 15#include <THStack.h>
0bd4b00f 16#include <TROOT.h>
f7cfc454 17#include <TParameter.h>
0bd4b00f 18#include <iostream>
19#include <iomanip>
7e4038b5 20
21ClassImp(AliFMDDensityCalculator)
22#if 0
23; // For Emacs
24#endif
25
26//____________________________________________________________________
27AliFMDDensityCalculator::AliFMDDensityCalculator()
28 : TNamed(),
29 fRingHistos(),
dd497217 30 fSumOfWeights(0),
31 fWeightedSum(0),
ea3e5d95 32 fCorrections(0),
0bd4b00f 33 fMaxParticles(5),
9d05ffeb 34 fUsePoisson(false),
8411f7fe 35 fUsePhiAcceptance(kPhiCorrectNch),
0bd4b00f 36 fAccI(0),
37 fAccO(0),
1174780f 38 fFMD1iMax(0),
39 fFMD2iMax(0),
40 fFMD2oMax(0),
41 fFMD3iMax(0),
42 fFMD3oMax(0),
5bb5d1f6 43 fMaxWeights(0),
44 fLowCuts(0),
21d778b1 45 fEtaLumping(32),
46 fPhiLumping(4),
d2638bb7 47 fDebug(0),
6f4a5c0d 48 fCuts(),
49 fRecalculateEta(false)
7984e5f7 50{
51 //
52 // Constructor
53 //
6ab100ec 54 DGUARD(fDebug, 0, "Default CTOR of FMD density calculator");
7984e5f7 55}
7e4038b5 56
57//____________________________________________________________________
58AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
59 : TNamed("fmdDensityCalculator", title),
60 fRingHistos(),
dd497217 61 fSumOfWeights(0),
62 fWeightedSum(0),
ea3e5d95 63 fCorrections(0),
0bd4b00f 64 fMaxParticles(5),
9d05ffeb 65 fUsePoisson(false),
8411f7fe 66 fUsePhiAcceptance(kPhiCorrectNch),
0bd4b00f 67 fAccI(0),
68 fAccO(0),
1174780f 69 fFMD1iMax(0),
70 fFMD2iMax(0),
71 fFMD2oMax(0),
72 fFMD3iMax(0),
73 fFMD3oMax(0),
5bb5d1f6 74 fMaxWeights(0),
75 fLowCuts(0),
21d778b1 76 fEtaLumping(32),
77 fPhiLumping(4),
d2638bb7 78 fDebug(0),
6f4a5c0d 79 fCuts(),
80 fRecalculateEta(false)
7e4038b5 81{
7984e5f7 82 //
83 // Constructor
84 //
85 // Parameters:
86 // name Name of object
87 //
6ab100ec 88 DGUARD(fDebug, 0, "Named CTOR of FMD density calculator: %s", title);
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
7e4038b5 121}
122
123//____________________________________________________________________
124AliFMDDensityCalculator::AliFMDDensityCalculator(const
125 AliFMDDensityCalculator& o)
126 : TNamed(o),
127 fRingHistos(),
dd497217 128 fSumOfWeights(o.fSumOfWeights),
129 fWeightedSum(o.fWeightedSum),
ea3e5d95 130 fCorrections(o.fCorrections),
0bd4b00f 131 fMaxParticles(o.fMaxParticles),
9d05ffeb 132 fUsePoisson(o.fUsePoisson),
5bb5d1f6 133 fUsePhiAcceptance(o.fUsePhiAcceptance),
0bd4b00f 134 fAccI(o.fAccI),
135 fAccO(o.fAccO),
1174780f 136 fFMD1iMax(o.fFMD1iMax),
137 fFMD2iMax(o.fFMD2iMax),
138 fFMD2oMax(o.fFMD2oMax),
139 fFMD3iMax(o.fFMD3iMax),
140 fFMD3oMax(o.fFMD3oMax),
5bb5d1f6 141 fMaxWeights(o.fMaxWeights),
142 fLowCuts(o.fLowCuts),
b6b35c77 143 fEtaLumping(o.fEtaLumping),
144 fPhiLumping(o.fPhiLumping),
d2638bb7 145 fDebug(o.fDebug),
6f4a5c0d 146 fCuts(o.fCuts),
147 fRecalculateEta(o.fRecalculateEta)
7e4038b5 148{
7984e5f7 149 //
150 // Copy constructor
151 //
152 // Parameters:
153 // o Object to copy from
154 //
6ab100ec 155 DGUARD(fDebug, 0, "Copy CTOR of FMD density calculator");
7e4038b5 156 TIter next(&o.fRingHistos);
157 TObject* obj = 0;
158 while ((obj = next())) fRingHistos.Add(obj);
159}
160
161//____________________________________________________________________
162AliFMDDensityCalculator::~AliFMDDensityCalculator()
163{
7984e5f7 164 //
165 // Destructor
166 //
6ab100ec 167 DGUARD(fDebug, 3, "DTOR of FMD density calculator");
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 //
6ab100ec 184 DGUARD(fDebug, 3, "Assignment of FMD density calculator");
d015ecfe 185 if (&o == this) return *this;
ea3e5d95 186 TNamed::operator=(o);
7e4038b5 187
d2638bb7 188 fDebug = o.fDebug;
189 fMaxParticles = o.fMaxParticles;
190 fUsePoisson = o.fUsePoisson;
191 fUsePhiAcceptance = o.fUsePhiAcceptance;
192 fAccI = o.fAccI;
193 fAccO = o.fAccO;
194 fFMD1iMax = o.fFMD1iMax;
195 fFMD2iMax = o.fFMD2iMax;
196 fFMD2oMax = o.fFMD2oMax;
197 fFMD3iMax = o.fFMD3iMax;
198 fFMD3oMax = o.fFMD3oMax;
199 fMaxWeights = o.fMaxWeights;
200 fLowCuts = o.fLowCuts;
201 fEtaLumping = o.fEtaLumping;
202 fPhiLumping = o.fPhiLumping;
203 fCuts = o.fCuts;
6f4a5c0d 204 fRecalculateEta = o.fRecalculateEta;
d2638bb7 205
7e4038b5 206 fRingHistos.Delete();
207 TIter next(&o.fRingHistos);
208 TObject* obj = 0;
209 while ((obj = next())) fRingHistos.Add(obj);
210
211 return *this;
212}
213
1174780f 214//____________________________________________________________________
215void
9d05ffeb 216AliFMDDensityCalculator::Init(const TAxis& axis)
1174780f 217{
218 // Intialize this sub-algorithm
219 //
220 // Parameters:
221 // etaAxis Not used
6ab100ec 222 DGUARD(fDebug, 1, "Initialize FMD density calculator");
1174780f 223 CacheMaxWeights();
4077f3e8 224
9d05ffeb 225 TIter next(&fRingHistos);
226 RingHistos* o = 0;
d2638bb7 227 while ((o = static_cast<RingHistos*>(next()))) {
e18cb8bd 228 o->Init(axis);
229 // o->fMultCut = fCuts.GetFixedCut(o->fDet, o->fRing);
230 // o->fPoisson.Init(o->fDet,o->fRing,fEtaLumping, fPhiLumping);
d2638bb7 231 }
1174780f 232}
233
7e4038b5 234//____________________________________________________________________
235AliFMDDensityCalculator::RingHistos*
236AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
237{
7984e5f7 238 //
239 // Get the ring histogram container
240 //
241 // Parameters:
242 // d Detector
243 // r Ring
244 //
245 // Return:
246 // Ring histogram container
247 //
7e4038b5 248 Int_t idx = -1;
249 switch (d) {
250 case 1: idx = 0; break;
251 case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
252 case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
253 }
e308a636 254 if (idx < 0 || idx >= fRingHistos.GetEntries()) {
255 AliWarning(Form("Index %d of FMD%d%c out of range", idx, d, r));
256 return 0;
257 }
7e4038b5 258
259 return static_cast<RingHistos*>(fRingHistos.At(idx));
260}
5bb5d1f6 261
0bd4b00f 262//____________________________________________________________________
263Double_t
5bb5d1f6 264AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Int_t ieta,
265 Bool_t errors) const
0bd4b00f 266{
7984e5f7 267 //
268 // Get the multiplicity cut. If the user has set fMultCut (via
269 // SetMultCut) then that value is used. If not, then the lower
270 // value of the fit range for the enery loss fits is returned.
271 //
272 // Return:
273 // Lower cut on multiplicity
274 //
d2638bb7 275 return fCuts.GetMultCut(d,r,ieta,errors);
5bb5d1f6 276}
277
278//____________________________________________________________________
279Double_t
280AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Double_t eta,
8e400b14 281 Bool_t errors) const
5bb5d1f6 282{
283 //
284 // Get the multiplicity cut. If the user has set fMultCut (via
285 // SetMultCut) then that value is used. If not, then the lower
286 // value of the fit range for the enery loss fits is returned.
287 //
288 // Return:
289 // Lower cut on multiplicity
290 //
d2638bb7 291 return fCuts.GetMultCut(d,r,eta,errors);
0bd4b00f 292}
293
7e4038b5 294//____________________________________________________________________
295Bool_t
296AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd,
297 AliForwardUtil::Histos& hists,
0bd4b00f 298 UShort_t vtxbin,
21d778b1 299 Bool_t lowFlux,
b3d81bee 300 Double_t /*cent*/,
6f4a5c0d 301 Double_t zvtx)
7e4038b5 302{
7984e5f7 303 //
304 // Do the calculations
305 //
306 // Parameters:
307 // fmd AliESDFMD object (possibly) corrected for sharing
308 // hists Histogram cache
309 // vtxBin Vertex bin
310 // lowFlux Low flux flag.
311 //
312 // Return:
313 // true on successs
6ab100ec 314 DGUARD(fDebug, 1, "Calculate density in FMD density calculator");
21d778b1 315
7e4038b5 316 for (UShort_t d=1; d<=3; d++) {
317 UShort_t nr = (d == 1 ? 1 : 2);
318 for (UShort_t q=0; q<nr; q++) {
319 Char_t r = (q == 0 ? 'I' : 'O');
320 UShort_t ns= (q == 0 ? 20 : 40);
321 UShort_t nt= (q == 0 ? 512 : 256);
322 TH2D* h = hists.Get(d,r);
323 RingHistos* rh= GetRingHistos(d,r);
e308a636 324 if (!rh) {
325 AliError(Form("No ring histogram found for FMD%d%c", d, r));
326 fRingHistos.ls();
327 return false;
328 }
e18cb8bd 329 // rh->fPoisson.SetObject(d,r,vtxbin,cent);
d23503ee 330 rh->fPoisson.Reset(0);
821ffd28 331 // rh->ResetPoissonHistos(h, fEtaLumping, fPhiLumping);
9fde7142 332
7e4038b5 333 for (UShort_t s=0; s<ns; s++) {
334 for (UShort_t t=0; t<nt; t++) {
3e4a0875 335
5bb5d1f6 336 Float_t mult = fmd.Multiplicity(d,r,s,t);
337 Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
338 Float_t eta = fmd.Eta(d,r,s,t);
6f4a5c0d 339
340
341 if(fRecalculateEta)
342 eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,zvtx);
9fde7142 343
5bb5d1f6 344 if (mult == AliESDFMD::kInvalidMult || mult > 20) {
21d778b1 345 rh->fPoisson.Fill(t , s, false);
5bb5d1f6 346 rh->fEvsM->Fill(mult,0);
33b078ef 347 continue;
348 }
8411f7fe 349 if (fUsePhiAcceptance == kPhiCorrectELoss)
0082a8fc 350 mult *= AcceptanceCorrection(r,t);
5bb5d1f6 351
8e400b14 352 Double_t cut = 1024;
353 if (eta != AliESDFMD::kInvalidEta) cut = GetMultCut(d, r, eta,false);
354
5bb5d1f6 355 Double_t n = 0;
356 if (cut > 0 && mult > cut)
357 n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
9fde7142 358
f7cfc454 359 rh->fELoss->Fill(mult);
7e4038b5 360 rh->fEvsN->Fill(mult,n);
0bd4b00f 361 rh->fEtaVsN->Fill(eta, n);
9fde7142 362
5bb5d1f6 363 Double_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
dd497217 364 fCorrections->Fill(c);
7e4038b5 365 if (c > 0) n /= c;
366 rh->fEvsM->Fill(mult,n);
0bd4b00f 367 rh->fEtaVsM->Fill(eta, n);
368 rh->fCorr->Fill(eta, c);
f7cfc454 369
821ffd28 370 Bool_t hit = (n > 0.9 && c > 0);
371 if (hit) rh->fELossUsed->Fill(mult);
21d778b1 372 rh->fPoisson.Fill(t,s,hit,1./c);
7e4038b5 373 h->Fill(eta,phi,n);
5bb5d1f6 374 if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n);
7e4038b5 375 } // for t
376 } // for s
21d778b1 377
378 TH2D* hclone = static_cast<TH2D*>(h->Clone("hclone"));
379 if (!fUsePoisson) hclone->Reset();
380 if ( fUsePoisson) h->Reset();
381
821ffd28 382 TH2D* poisson = rh->fPoisson.Result();
d23503ee 383 for (Int_t t=0; t < poisson->GetNbinsX(); t++) {
384 for (Int_t s=0; s < poisson->GetNbinsY(); s++) {
21d778b1 385
386 Double_t poissonV = poisson->GetBinContent(t+1,s+1);
387 Double_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
388 Double_t eta = fmd.Eta(d,r,s,t);
6f4a5c0d 389 if(fRecalculateEta)
390 eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,zvtx);
21d778b1 391 if (fUsePoisson)
392 h->Fill(eta,phi,poissonV);
393 else
394 hclone->Fill(eta,phi,poissonV);
a7f31ad8 395 if (fUsePoisson) rh->fDensity->Fill(eta, phi, poissonV);
21d778b1 396 }
397 }
398
399 for (Int_t ieta=1; ieta <= h->GetNbinsX(); ieta++) {
400 for (Int_t iphi=1; iphi<= h->GetNbinsY(); iphi++) {
e83d0620 401
21d778b1 402 Double_t poissonV = 0; //h->GetBinContent(,s+1);
403 Double_t eLossV = 0;
404 if(fUsePoisson) {
405 poissonV = h->GetBinContent(ieta,iphi);
406 eLossV = hclone->GetBinContent(ieta,iphi);
407 }
408 else {
409 poissonV = hclone->GetBinContent(ieta,iphi);
410 eLossV = h->GetBinContent(ieta,iphi);
411 }
7b03929c 412
9d05ffeb 413 rh->fELossVsPoisson->Fill(eLossV, poissonV);
9d05ffeb 414 }
415 }
21d778b1 416 delete hclone;
417
7e4038b5 418 } // for q
419 } // for d
a7f31ad8 420
7e4038b5 421 return kTRUE;
422}
423
1174780f 424//_____________________________________________________________________
425Int_t
fb3430ac 426AliFMDDensityCalculator::FindMaxWeight(const AliFMDCorrELossFit* cor,
1174780f 427 UShort_t d, Char_t r, Int_t iEta) const
428{
fb3430ac 429 //
430 // Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
431 //
432 // Parameters:
433 // cor Correction
434 // d Detector
435 // r Ring
436 // iEta Eta bin
437 //
6ab100ec 438 DGUARD(fDebug, 3, "Find maximum weight in FMD density calculator");
f61b1745 439 if(!cor) return -1;
440
1174780f 441 AliFMDCorrELossFit::ELossFit* fit = cor->GetFit(d,r,iEta);
442 if (!fit) {
443 // AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
444 return -1;
445 }
1f480471 446 return TMath::Min(Int_t(fMaxParticles), fit->FindMaxWeight());
1174780f 447}
448
449//_____________________________________________________________________
450void
451AliFMDDensityCalculator::CacheMaxWeights()
452{
fb3430ac 453 //
454 // Find the max weights and cache them
455 //
6ab100ec 456 DGUARD(fDebug, 2, "Cache maximum weights in FMD density calculator");
1174780f 457 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
458 AliFMDCorrELossFit* cor = fcm.GetELossFit();
f61b1745 459 TAxis eta(240,-6,6);
460 if(cor)
461 eta.Set(cor->GetEtaAxis().GetNbins(),
462 cor->GetEtaAxis().GetXmin(),
463 cor->GetEtaAxis().GetXmax());
1174780f 464
465 Int_t nEta = eta.GetNbins();
466 fFMD1iMax.Set(nEta);
467 fFMD2iMax.Set(nEta);
468 fFMD2oMax.Set(nEta);
469 fFMD3iMax.Set(nEta);
470 fFMD3oMax.Set(nEta);
471
5bb5d1f6 472 fMaxWeights->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
473 fMaxWeights->GetYaxis()->SetBinLabel(1, "FMD1i");
474 fMaxWeights->GetYaxis()->SetBinLabel(2, "FMD2i");
475 fMaxWeights->GetYaxis()->SetBinLabel(3, "FMD2o");
476 fMaxWeights->GetYaxis()->SetBinLabel(4, "FMD3i");
477 fMaxWeights->GetYaxis()->SetBinLabel(5, "FMD3o");
478
479 AliInfo(Form("Get eta axis with %d bins from %f to %f",
480 nEta, eta.GetXmin(), eta.GetXmax()));
481 fLowCuts->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
482 fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
483 fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
484 fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
485 fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
486 fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
487
1174780f 488 for (Int_t i = 0; i < nEta; i++) {
5bb5d1f6 489 Double_t w[5];
490 w[0] = fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1);
491 w[1] = fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1);
492 w[2] = fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1);
493 w[3] = fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1);
494 w[4] = fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1);
495 Double_t l[5];
496 l[0] = GetMultCut(1, 'I', i+1, false);
497 l[1] = GetMultCut(2, 'I', i+1, false);
498 l[2] = GetMultCut(2, 'O', i+1, false);
499 l[3] = GetMultCut(3, 'I', i+1, false);
500 l[4] = GetMultCut(3, 'O', i+1, false);
501 for (Int_t j = 0; j < 5; j++) {
502 if (w[j] > 0) fMaxWeights->SetBinContent(i+1, j+1, w[j]);
503 if (l[j] > 0) fLowCuts->SetBinContent(i+1, j+1, l[j]);
504 }
1174780f 505 }
506}
507
508//_____________________________________________________________________
509Int_t
510AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const
511{
fb3430ac 512 //
513 // Find the (cached) maximum weight for FMD<i>dr</i> in
514 // @f$\eta@f$ bin @a iEta
515 //
516 // Parameters:
517 // d Detector
518 // r Ring
519 // iEta Eta bin
520 //
521 // Return:
522 // max weight or <= 0 in case of problems
523 //
1174780f 524 if (iEta < 0) return -1;
525
526 const TArrayI* max = 0;
527 switch (d) {
528 case 1: max = &fFMD1iMax; break;
529 case 2: max = (r == 'I' || r == 'i' ? &fFMD2iMax : &fFMD2oMax); break;
530 case 3: max = (r == 'I' || r == 'i' ? &fFMD3iMax : &fFMD3oMax); break;
531 }
532 if (!max) {
533 AliWarning(Form("No array for FMD%d%c", d, r));
534 return -1;
535 }
536
537 if (iEta >= max->fN) {
538 AliWarning(Form("Eta bin %3d out of bounds [0,%d]",
539 iEta, max->fN-1));
540 return -1;
541 }
542
543 AliDebug(30,Form("Max weight for FMD%d%c eta bin %3d: %d", d, r, iEta,
544 max->At(iEta)));
545 return max->At(iEta);
546}
547
548//_____________________________________________________________________
549Int_t
550AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const
551{
fb3430ac 552 //
553 // Find the (cached) maximum weight for FMD<i>dr</i> iat
554 // @f$\eta@f$
555 //
556 // Parameters:
557 // d Detector
558 // r Ring
559 // eta Eta bin
560 //
561 // Return:
562 // max weight or <= 0 in case of problems
563 //
1174780f 564 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
565 Int_t iEta = fcm.GetELossFit()->FindEtaBin(eta) -1;
566
567 return GetMaxWeight(d, r, iEta);
568}
569
7e4038b5 570//_____________________________________________________________________
571Float_t
572AliFMDDensityCalculator::NParticles(Float_t mult,
573 UShort_t d,
574 Char_t r,
575 UShort_t /*s*/,
576 UShort_t /*t*/,
0bd4b00f 577 UShort_t /*v*/,
7e4038b5 578 Float_t eta,
579 Bool_t lowFlux) const
580{
7984e5f7 581 //
582 // Get the number of particles corresponding to the signal mult
583 //
584 // Parameters:
585 // mult Signal
586 // d Detector
587 // r Ring
588 // s Sector
589 // t Strip (not used)
590 // v Vertex bin
591 // eta Pseudo-rapidity
592 // lowFlux Low-flux flag
593 //
594 // Return:
595 // The number of particles
596 //
5bb5d1f6 597 // if (mult <= GetMultCut()) return 0;
6ab100ec 598 DGUARD(fDebug, 3, "Calculate Nch in FMD density calculator");
7e4038b5 599 if (lowFlux) return 1;
0bd4b00f 600
0bd4b00f 601 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
602 AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
603 if (!fit) {
604 AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
605 return 0;
606 }
607
1174780f 608 Int_t m = GetMaxWeight(d,r,eta); // fit->FindMaxWeight();
0bd4b00f 609 if (m < 1) {
610 AliWarning(Form("No good fits for FMD%d%c at eta=%f", d, r, eta));
611 return 0;
612 }
79909b8b 613
0bd4b00f 614 UShort_t n = TMath::Min(fMaxParticles, UShort_t(m));
615 Double_t ret = fit->EvaluateWeighted(mult, n);
79909b8b 616
0bd4b00f 617 if (fDebug > 10) {
618 AliInfo(Form("FMD%d%c, eta=%7.4f, %8.5f -> %8.5f", d, r, eta, mult, ret));
619 }
79909b8b 620
0bd4b00f 621 fWeightedSum->Fill(ret);
622 fSumOfWeights->Fill(ret);
3e4a0875 623
0bd4b00f 624 return ret;
7e4038b5 625}
626
627//_____________________________________________________________________
628Float_t
0bd4b00f 629AliFMDDensityCalculator::Correction(UShort_t d,
630 Char_t r,
631 UShort_t /*s*/,
632 UShort_t t,
633 UShort_t /*v*/,
634 Float_t eta,
635 Bool_t lowFlux) const
7e4038b5 636{
7984e5f7 637 //
638 // Get the inverse correction factor. This consist of
639 //
640 // - acceptance correction (corners of sensors)
641 // - double hit correction (for low-flux events)
642 // - dead strip correction
643 //
644 // Parameters:
645 // d Detector
646 // r Ring
647 // s Sector
648 // t Strip (not used)
649 // v Vertex bin
650 // eta Pseudo-rapidity
651 // lowFlux Low-flux flag
652 //
653 // Return:
654 //
655 //
6ab100ec 656 DGUARD(fDebug, 3, "Apply correction in FMD density calculator");
0bd4b00f 657 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
7e4038b5 658
5bb5d1f6 659 Float_t correction = 1;
8411f7fe 660 if (fUsePhiAcceptance == kPhiCorrectNch)
661 correction = AcceptanceCorrection(r,t);
7e4038b5 662 if (lowFlux) {
0bd4b00f 663 TH1D* dblHitCor = 0;
664 if (fcm.GetDoubleHit())
665 dblHitCor = fcm.GetDoubleHit()->GetCorrection(d,r);
666
7e4038b5 667 if (dblHitCor) {
0bd4b00f 668 Double_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
669 if (dblC > 0) correction *= dblC;
7e4038b5 670 }
0bd4b00f 671 else {
7e4038b5 672 AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
0bd4b00f 673 }
7e4038b5 674 }
7e4038b5 675 return correction;
676}
677
0bd4b00f 678//_____________________________________________________________________
679TH1D*
680AliFMDDensityCalculator::GenerateAcceptanceCorrection(Char_t r) const
681{
7984e5f7 682 //
683 // Generate the acceptance corrections
684 //
685 // Parameters:
686 // r Ring to generate for
687 //
688 // Return:
689 // Newly allocated histogram of acceptance corrections
690 //
6ab100ec 691 DGUARD(fDebug, 3, "Make acceptance correction in FMD density calculator");
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 //
6ab100ec 786 DGUARD(fDebug, 1, "Scale histograms in FMD density calculator");
7e4038b5 787 if (nEvents <= 0) return;
9d99b0dd 788 TList* d = static_cast<TList*>(dir->FindObject(GetName()));
789 if (!d) return;
7e4038b5 790
791 TIter next(&fRingHistos);
792 RingHistos* o = 0;
5bb5d1f6 793 THStack* sums = new THStack("sums", "sums of ring signals");
794 while ((o = static_cast<RingHistos*>(next()))) {
9d99b0dd 795 o->ScaleHistograms(d, nEvents);
5bb5d1f6 796 TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
797 o->fDensity->GetNbinsY(),"e");
798 sum->Scale(1., "width");
799 sum->SetTitle(o->GetName());
800 sum->SetDirectory(0);
801 sum->SetYTitle("#sum N_{ch,incl}");
802 sums->Add(sum);
803 }
804 d->Add(sums);
7e4038b5 805}
806
807//____________________________________________________________________
808void
9d99b0dd 809AliFMDDensityCalculator::DefineOutput(TList* dir)
7e4038b5 810{
7984e5f7 811 //
812 // Output diagnostic histograms to directory
813 //
814 // Parameters:
815 // dir List to write in
816 //
6ab100ec 817 DGUARD(fDebug, 1, "Define output FMD density calculator");
7e4038b5 818 TList* d = new TList;
5bb5d1f6 819 d->SetOwner();
7e4038b5 820 d->SetName(GetName());
821 dir->Add(d);
dd497217 822 d->Add(fWeightedSum);
823 d->Add(fSumOfWeights);
824 d->Add(fCorrections);
0bd4b00f 825 d->Add(fAccI);
826 d->Add(fAccO);
5bb5d1f6 827 d->Add(fMaxWeights);
828 d->Add(fLowCuts);
0bd4b00f 829
d2638bb7 830 // TNamed* sigma = new TNamed("sigma",
831 // (fIncludeSigma ? "included" : "excluded"));
f7cfc454 832 TNamed* maxP = new TNamed("maxParticle", Form("%d", fMaxParticles));
833 TNamed* method = new TNamed("method",
834 (fUsePoisson ? "Poisson" : "Energy loss"));
835 TNamed* phiA = new TNamed("phiAcceptance",
0082a8fc 836 (fUsePhiAcceptance == 0 ? "disabled" :
837 fUsePhiAcceptance == 1 ? "particles" :
838 "energy loss"));
f7cfc454 839 TNamed* etaL = new TNamed("etaLumping", Form("%d", fEtaLumping));
840 TNamed* phiL = new TNamed("phiLumping", Form("%d", fPhiLumping));
d2638bb7 841 // TParameter<double>* nxi = new TParameter<double>("nXi", fNXi);
f7cfc454 842
d2638bb7 843 // d->Add(sigma);
f7cfc454 844 d->Add(maxP);
845 d->Add(method);
846 d->Add(phiA);
847 d->Add(etaL);
848 d->Add(phiL);
d2638bb7 849 // d->Add(nxi);
850 fCuts.Output(d,0);
f7cfc454 851
7e4038b5 852 TIter next(&fRingHistos);
853 RingHistos* o = 0;
854 while ((o = static_cast<RingHistos*>(next()))) {
e18cb8bd 855 o->fPoisson.SetLumping(fEtaLumping, fPhiLumping);
7e4038b5 856 o->Output(d);
857 }
858}
0bd4b00f 859//____________________________________________________________________
860void
1174780f 861AliFMDDensityCalculator::Print(Option_t* option) const
0bd4b00f 862{
7984e5f7 863 //
864 // Print information
865 //
866 // Parameters:
867 // option Not used
868 //
1174780f 869 char ind[gROOT->GetDirLevel()+3];
0bd4b00f 870 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
871 ind[gROOT->GetDirLevel()] = '\0';
1f480471 872 std::cout << ind << ClassName() << ": " << GetName() << '\n'
5bb5d1f6 873 << std::boolalpha
1f480471 874 << ind << " Max(particles): " << fMaxParticles << '\n'
5bb5d1f6 875 << ind << " Poisson method: " << fUsePoisson << '\n'
876 << ind << " Use phi acceptance: " << fUsePhiAcceptance << '\n'
1f480471 877 << ind << " Eta lumping: " << fEtaLumping << '\n'
878 << ind << " Phi lumping: " << fPhiLumping << '\n'
5bb5d1f6 879 << std::noboolalpha
1f480471 880 << std::flush;
d2638bb7 881 std::cout << ind << " Lower cut:" << std::endl;
882 fCuts.Print();
1174780f 883 TString opt(option);
884 opt.ToLower();
885 if (opt.Contains("nomax")) return;
886
887 std::cout << ind << " Max weights:\n";
888
889 for (UShort_t d=1; d<=3; d++) {
890 UShort_t nr = (d == 1 ? 1 : 2);
891 for (UShort_t q=0; q<nr; q++) {
892 ind[gROOT->GetDirLevel()] = ' ';
893 ind[gROOT->GetDirLevel()+1] = '\0';
894 Char_t r = (q == 0 ? 'I' : 'O');
895 std::cout << ind << " FMD" << d << r << ":";
896 ind[gROOT->GetDirLevel()+1] = ' ';
897 ind[gROOT->GetDirLevel()+2] = '\0';
898
899 const TArrayI& a = (d == 1 ? fFMD1iMax :
900 (d == 2 ? (r == 'I' ? fFMD2iMax : fFMD2oMax) :
901 (r == 'I' ? fFMD3iMax : fFMD3oMax)));
902 Int_t j = 0;
903 for (Int_t i = 0; i < a.fN; i++) {
904 if (a.fArray[i] < 1) continue;
905 if (j % 6 == 0) std::cout << "\n " << ind;
906 j++;
907 std::cout << " " << std::setw(3) << i << ": " << a.fArray[i];
908 }
909 std::cout << std::endl;
910 }
911 }
0bd4b00f 912}
7e4038b5 913
914//====================================================================
915AliFMDDensityCalculator::RingHistos::RingHistos()
9d99b0dd 916 : AliForwardUtil::RingHistos(),
7e4038b5 917 fEvsN(0),
918 fEvsM(0),
0bd4b00f 919 fEtaVsN(0),
920 fEtaVsM(0),
921 fCorr(0),
9d05ffeb 922 fDensity(0),
923 fELossVsPoisson(0),
821ffd28 924 fPoisson(),
f7cfc454 925 fELoss(0),
926 fELossUsed(0),
927 fMultCut(0)
7984e5f7 928{
929 //
930 // Default CTOR
931 //
932}
7e4038b5 933
934//____________________________________________________________________
935AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
9d99b0dd 936 : AliForwardUtil::RingHistos(d,r),
7e4038b5 937 fEvsN(0),
938 fEvsM(0),
0bd4b00f 939 fEtaVsN(0),
940 fEtaVsM(0),
941 fCorr(0),
9d05ffeb 942 fDensity(0),
943 fELossVsPoisson(0),
821ffd28 944 fPoisson("ignored"),
d2638bb7 945 fELoss(0),
f7cfc454 946 fELossUsed(0),
947 fMultCut(0)
7e4038b5 948{
7984e5f7 949 //
950 // Constructor
951 //
952 // Parameters:
953 // d detector
954 // r ring
955 //
5bb5d1f6 956 fEvsN = new TH2D("elossVsNnocorr",
957 "#Delta E/#Delta E_{mip} vs uncorrected inclusive N_{ch}",
958 250, -.5, 24.5, 250, -.5, 24.5);
7e4038b5 959 fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
960 fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
961 fEvsN->Sumw2();
962 fEvsN->SetDirectory(0);
5bb5d1f6 963
964 fEvsM = static_cast<TH2D*>(fEvsN->Clone("elossVsNcorr"));
965 fEvsM->SetTitle("#Delta E/#Delta E_{mip} vs corrected inclusive N_{ch}");
7e4038b5 966 fEvsM->SetDirectory(0);
967
5bb5d1f6 968 fEtaVsN = new TProfile("etaVsNnocorr",
969 "Average inclusive N_{ch} vs #eta (uncorrected)",
970 200, -4, 6);
0bd4b00f 971 fEtaVsN->SetXTitle("#eta");
972 fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
973 fEtaVsN->SetDirectory(0);
974 fEtaVsN->SetLineColor(Color());
975 fEtaVsN->SetFillColor(Color());
5bb5d1f6 976
977 fEtaVsM = static_cast<TProfile*>(fEtaVsN->Clone("etaVsNcorr"));
978 fEtaVsM->SetTitle("Average inclusive N_{ch} vs #eta (corrected)");
0bd4b00f 979 fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
980 fEtaVsM->SetDirectory(0);
0bd4b00f 981
982
5bb5d1f6 983 fCorr = new TProfile("corr", "Average correction", 200, -4, 6);
0bd4b00f 984 fCorr->SetXTitle("#eta");
985 fCorr->SetYTitle("#LT correction#GT");
986 fCorr->SetDirectory(0);
987 fCorr->SetLineColor(Color());
988 fCorr->SetFillColor(Color());
989
5bb5d1f6 990 fDensity = new TH2D("inclDensity", "Inclusive N_{ch} density",
7e4038b5 991 200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
992 0, 2*TMath::Pi());
993 fDensity->SetDirectory(0);
5bb5d1f6 994 fDensity->Sumw2();
995 fDensity->SetMarkerColor(Color());
7e4038b5 996 fDensity->SetXTitle("#eta");
997 fDensity->SetYTitle("#phi [radians]");
998 fDensity->SetZTitle("Inclusive N_{ch} density");
9d05ffeb 999
5bb5d1f6 1000 fELossVsPoisson = new TH2D("elossVsPoisson",
1001 "N_{ch} from energy loss vs from Poission",
3e4a0875 1002 500, 0, 100, 500, 0, 100);
9d05ffeb 1003 fELossVsPoisson->SetDirectory(0);
1004 fELossVsPoisson->SetXTitle("N_{ch} from #DeltaE");
1005 fELossVsPoisson->SetYTitle("N_{ch} from Poisson");
1006 fELossVsPoisson->SetZTitle("Correlation");
1007
f7cfc454 1008 fELoss = new TH1D("eloss", "#Delta/#Delta_{mip} in all strips",
1009 600, 0, 15);
1010 fELoss->SetXTitle("#Delta/#Delta_{mip} (selected)");
1011 fELoss->SetYTitle("P(#Delta/#Delta_{mip})");
1012 fELoss->SetFillColor(Color()-2);
1013 fELoss->SetFillStyle(3003);
1014 fELoss->SetLineColor(kBlack);
1015 fELoss->SetLineStyle(2);
1016 fELoss->SetLineWidth(2);
1017 fELoss->SetDirectory(0);
1018
1019 fELossUsed = static_cast<TH1D*>(fELoss->Clone("elossUsed"));
1020 fELossUsed->SetTitle("#Delta/#Delta_{mip} in used strips");
1021 fELossUsed->SetFillStyle(3002);
1022 fELossUsed->SetLineStyle(1);
1023 fELossUsed->SetDirectory(0);
1024
7e4038b5 1025}
1026//____________________________________________________________________
1027AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
9d99b0dd 1028 : AliForwardUtil::RingHistos(o),
7e4038b5 1029 fEvsN(o.fEvsN),
1030 fEvsM(o.fEvsM),
0bd4b00f 1031 fEtaVsN(o.fEtaVsN),
1032 fEtaVsM(o.fEtaVsM),
1033 fCorr(o.fCorr),
9d05ffeb 1034 fDensity(o.fDensity),
1035 fELossVsPoisson(o.fELossVsPoisson),
821ffd28 1036 fPoisson(o.fPoisson),
f7cfc454 1037 fELoss(o.fELoss),
1038 fELossUsed(o.fELossUsed),
1039 fMultCut(o.fMultCut)
7984e5f7 1040{
1041 //
1042 // Copy constructor
1043 //
1044 // Parameters:
1045 // o Object to copy from
1046 //
1047}
7e4038b5 1048
1049//____________________________________________________________________
1050AliFMDDensityCalculator::RingHistos&
1051AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
1052{
7984e5f7 1053 //
1054 // Assignment operator
1055 //
1056 // Parameters:
1057 // o Object to assign from
1058 //
1059 // Return:
1060 // Reference to this
1061 //
d015ecfe 1062 if (&o == this) return *this;
9d99b0dd 1063 AliForwardUtil::RingHistos::operator=(o);
7e4038b5 1064
9d05ffeb 1065 if (fEvsN) delete fEvsN;
1066 if (fEvsM) delete fEvsM;
1067 if (fEtaVsN) delete fEtaVsN;
1068 if (fEtaVsM) delete fEtaVsM;
1069 if (fCorr) delete fCorr;
1070 if (fDensity) delete fDensity;
1071 if (fELossVsPoisson) delete fELossVsPoisson;
7e4038b5 1072
9d05ffeb 1073 fEvsN = static_cast<TH2D*>(o.fEvsN->Clone());
1074 fEvsM = static_cast<TH2D*>(o.fEvsM->Clone());
1075 fEtaVsN = static_cast<TProfile*>(o.fEtaVsN->Clone());
1076 fEtaVsM = static_cast<TProfile*>(o.fEtaVsM->Clone());
1077 fCorr = static_cast<TProfile*>(o.fCorr->Clone());
1078 fDensity = static_cast<TH2D*>(o.fDensity->Clone());
f7cfc454 1079 fELossVsPoisson = static_cast<TH2D*>(o.fELossVsPoisson->Clone());
821ffd28 1080 fPoisson = o.fPoisson;
f7cfc454 1081 fELoss = static_cast<TH1D*>(o.fELoss->Clone());
1082 fELossUsed = static_cast<TH1D*>(o.fELossUsed->Clone());
7e4038b5 1083
1084 return *this;
1085}
1086//____________________________________________________________________
1087AliFMDDensityCalculator::RingHistos::~RingHistos()
1088{
7984e5f7 1089 //
1090 // Destructor
1091 //
9d05ffeb 1092}
1093
e308a636 1094
1095//____________________________________________________________________
9d05ffeb 1096void
9fde7142 1097AliFMDDensityCalculator::RingHistos::Init(const TAxis& /*eAxis*/)
9d05ffeb 1098{
e18cb8bd 1099 fPoisson.Init(-1,-1);
7e4038b5 1100}
1101
1102//____________________________________________________________________
1103void
1104AliFMDDensityCalculator::RingHistos::Output(TList* dir)
1105{
7984e5f7 1106 //
1107 // Make output
1108 //
1109 // Parameters:
1110 // dir Where to put it
1111 //
9d99b0dd 1112 TList* d = DefineOutputList(dir);
7e4038b5 1113 d->Add(fEvsN);
1114 d->Add(fEvsM);
0bd4b00f 1115 d->Add(fEtaVsN);
1116 d->Add(fEtaVsM);
1117 d->Add(fCorr);
7e4038b5 1118 d->Add(fDensity);
9d05ffeb 1119 d->Add(fELossVsPoisson);
821ffd28 1120 fPoisson.Output(d);
d23503ee 1121 fPoisson.GetOccupancy()->SetFillColor(Color());
1122 fPoisson.GetMean()->SetFillColor(Color());
1123 fPoisson.GetOccupancy()->SetFillColor(Color());
f7cfc454 1124 d->Add(fELoss);
1125 d->Add(fELossUsed);
d23503ee 1126
1127 Bool_t inner = (fRing == 'I' || fRing == 'i');
1128 Int_t nStr = inner ? 512 : 256;
1129 Int_t nSec = inner ? 20 : 40;
1130 TAxis x(nStr, -.5, nStr-.5);
1131 TAxis y(nSec, -.5, nSec-.5);
1132 x.SetTitle("strip");
1133 y.SetTitle("sector");
1134 fPoisson.Define(x, y);
1135
f7cfc454 1136 TParameter<double>* cut = new TParameter<double>("cut", fMultCut);
1137 d->Add(cut);
9d99b0dd 1138}
1139
1140//____________________________________________________________________
1141void
1142AliFMDDensityCalculator::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
1143{
7984e5f7 1144 //
1145 // Scale the histograms to the total number of events
1146 //
1147 // Parameters:
1148 // dir Where the output is
1149 // nEvents Number of events
1150 //
9d99b0dd 1151 TList* l = GetOutputList(dir);
1152 if (!l) return;
1153
5bb5d1f6 1154 TH1* density = GetOutputHist(l,"inclDensity");
9d99b0dd 1155 if (density) density->Scale(1./nEvents);
7e4038b5 1156}
1157
1158//____________________________________________________________________
1159//
1160// EOF
1161//
1162
1163
1164