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