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