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