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