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