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