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