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