1 // This class calculates the inclusive charged particle density
2 // in each for the 5 FMD rings.
4 #include "AliFMDDensityCalculator.h"
9 #include "AliForwardCorrectionManager.h"
10 #include "AliFMDCorrDoubleHit.h"
11 #include "AliFMDCorrELossFit.h"
17 #include <TParameter.h>
21 ClassImp(AliFMDDensityCalculator)
26 //____________________________________________________________________
27 AliFMDDensityCalculator::AliFMDDensityCalculator()
35 fUsePhiAcceptance(kPhiCorrectNch),
49 fRecalculateEta(false)
54 DGUARD(fDebug, 0, "Default CTOR of FMD density calculator");
57 //____________________________________________________________________
58 AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
59 : TNamed("fmdDensityCalculator", title),
66 fUsePhiAcceptance(kPhiCorrectNch),
80 fRecalculateEta(false)
86 // name Name of object
88 DGUARD(fDebug, 0, "Named CTOR of FMD density calculator: %s", title);
89 fRingHistos.SetName(GetName());
90 fRingHistos.SetOwner();
91 fRingHistos.Add(new RingHistos(1, 'I'));
92 fRingHistos.Add(new RingHistos(2, 'I'));
93 fRingHistos.Add(new RingHistos(2, 'O'));
94 fRingHistos.Add(new RingHistos(3, 'I'));
95 fRingHistos.Add(new RingHistos(3, 'O'));
96 fSumOfWeights = new TH1D("sumOfWeights", "Sum of Landau weights",
98 fSumOfWeights->SetFillColor(kRed+1);
99 fSumOfWeights->SetXTitle("#sum_{i} a_{i} f_{i}(#Delta)");
100 fWeightedSum = new TH1D("weightedSum", "Weighted sum of Landau propability",
102 fWeightedSum->SetFillColor(kBlue+1);
103 fWeightedSum->SetXTitle("#sum_{i} i a_{i} f_{i}(#Delta)");
104 fCorrections = new TH1D("corrections", "Distribution of corrections",
106 fCorrections->SetFillColor(kBlue+1);
107 fCorrections->SetXTitle("correction");
109 fAccI = GenerateAcceptanceCorrection('I');
110 fAccO = GenerateAcceptanceCorrection('O');
112 fMaxWeights = new TH2D("maxWeights", "Maximum i of a_{i}'s to use",
114 fMaxWeights->SetXTitle("#eta");
115 fMaxWeights->SetDirectory(0);
117 fLowCuts = new TH2D("lowCuts", "Low cuts used", 1, 0, 1, 1, 0, 1);
118 fLowCuts->SetXTitle("#eta");
119 fLowCuts->SetDirectory(0);
123 //____________________________________________________________________
124 AliFMDDensityCalculator::AliFMDDensityCalculator(const
125 AliFMDDensityCalculator& o)
128 fSumOfWeights(o.fSumOfWeights),
129 fWeightedSum(o.fWeightedSum),
130 fCorrections(o.fCorrections),
131 fMaxParticles(o.fMaxParticles),
132 fUsePoisson(o.fUsePoisson),
133 fUsePhiAcceptance(o.fUsePhiAcceptance),
136 fFMD1iMax(o.fFMD1iMax),
137 fFMD2iMax(o.fFMD2iMax),
138 fFMD2oMax(o.fFMD2oMax),
139 fFMD3iMax(o.fFMD3iMax),
140 fFMD3oMax(o.fFMD3oMax),
141 fMaxWeights(o.fMaxWeights),
142 fLowCuts(o.fLowCuts),
143 fEtaLumping(o.fEtaLumping),
144 fPhiLumping(o.fPhiLumping),
147 fRecalculateEta(o.fRecalculateEta)
153 // o Object to copy from
155 DGUARD(fDebug, 0, "Copy CTOR of FMD density calculator");
156 TIter next(&o.fRingHistos);
158 while ((obj = next())) fRingHistos.Add(obj);
161 //____________________________________________________________________
162 AliFMDDensityCalculator::~AliFMDDensityCalculator()
167 DGUARD(fDebug, 3, "DTOR of FMD density calculator");
168 fRingHistos.Delete();
171 //____________________________________________________________________
172 AliFMDDensityCalculator&
173 AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
176 // Assignement operator
179 // o Object to assign from
182 // Reference to this object
184 DGUARD(fDebug, 3, "Assignment of FMD density calculator");
185 if (&o == this) return *this;
186 TNamed::operator=(o);
189 fMaxParticles = o.fMaxParticles;
190 fUsePoisson = o.fUsePoisson;
191 fUsePhiAcceptance = o.fUsePhiAcceptance;
194 fFMD1iMax = o.fFMD1iMax;
195 fFMD2iMax = o.fFMD2iMax;
196 fFMD2oMax = o.fFMD2oMax;
197 fFMD3iMax = o.fFMD3iMax;
198 fFMD3oMax = o.fFMD3oMax;
199 fMaxWeights = o.fMaxWeights;
200 fLowCuts = o.fLowCuts;
201 fEtaLumping = o.fEtaLumping;
202 fPhiLumping = o.fPhiLumping;
204 fRecalculateEta = o.fRecalculateEta;
206 fRingHistos.Delete();
207 TIter next(&o.fRingHistos);
209 while ((obj = next())) fRingHistos.Add(obj);
214 //____________________________________________________________________
216 AliFMDDensityCalculator::Init(const TAxis& axis)
218 // Intialize this sub-algorithm
222 DGUARD(fDebug, 1, "Initialize FMD density calculator");
223 CacheMaxWeights(axis);
225 TIter next(&fRingHistos);
227 while ((o = static_cast<RingHistos*>(next()))) {
229 // o->fMultCut = fCuts.GetFixedCut(o->fDet, o->fRing);
230 // o->fPoisson.Init(o->fDet,o->fRing,fEtaLumping, fPhiLumping);
234 //____________________________________________________________________
235 AliFMDDensityCalculator::RingHistos*
236 AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
239 // Get the ring histogram container
246 // Ring histogram container
250 case 1: idx = 0; break;
251 case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
252 case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
254 if (idx < 0 || idx >= fRingHistos.GetEntries()) {
255 AliWarning(Form("Index %d of FMD%d%c out of range", idx, d, r));
259 return static_cast<RingHistos*>(fRingHistos.At(idx));
262 //____________________________________________________________________
264 AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Int_t ieta,
268 // Get the multiplicity cut. If the user has set fMultCut (via
269 // SetMultCut) then that value is used. If not, then the lower
270 // value of the fit range for the enery loss fits is returned.
273 // Lower cut on multiplicity
275 return fCuts.GetMultCut(d,r,ieta,errors);
278 //____________________________________________________________________
280 AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Double_t eta,
284 // Get the multiplicity cut. If the user has set fMultCut (via
285 // SetMultCut) then that value is used. If not, then the lower
286 // value of the fit range for the enery loss fits is returned.
289 // Lower cut on multiplicity
291 return fCuts.GetMultCut(d,r,eta,errors);
294 //____________________________________________________________________
296 AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd,
297 AliForwardUtil::Histos& hists,
304 // Do the calculations
307 // fmd AliESDFMD object (possibly) corrected for sharing
308 // hists Histogram cache
310 // lowFlux Low flux flag.
314 DGUARD(fDebug, 1, "Calculate density in FMD density calculator");
316 for (UShort_t d=1; d<=3; d++) {
317 UShort_t nr = (d == 1 ? 1 : 2);
318 for (UShort_t q=0; q<nr; q++) {
319 Char_t r = (q == 0 ? 'I' : 'O');
320 UShort_t ns= (q == 0 ? 20 : 40);
321 UShort_t nt= (q == 0 ? 512 : 256);
322 TH2D* h = hists.Get(d,r);
323 RingHistos* rh= GetRingHistos(d,r);
325 AliError(Form("No ring histogram found for FMD%d%c", d, r));
329 // rh->fPoisson.SetObject(d,r,vtxbin,cent);
330 rh->fPoisson.Reset(0);
331 // rh->ResetPoissonHistos(h, fEtaLumping, fPhiLumping);
333 for (UShort_t s=0; s<ns; s++) {
334 for (UShort_t t=0; t<nt; t++) {
336 Float_t mult = fmd.Multiplicity(d,r,s,t);
337 Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
338 Float_t eta = fmd.Eta(d,r,s,t);
342 eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,zvtx);
344 if (mult == AliESDFMD::kInvalidMult || mult > 20) {
345 rh->fPoisson.Fill(t , s, false);
346 rh->fEvsM->Fill(mult,0);
349 if (fUsePhiAcceptance == kPhiCorrectELoss)
350 mult *= AcceptanceCorrection(r,t);
353 if (eta != AliESDFMD::kInvalidEta) cut = GetMultCut(d, r, eta,false);
356 if (cut > 0 && mult > cut)
357 n = NParticles(mult,d,r,s,t,vtxbin,eta,lowFlux);
359 rh->fELoss->Fill(mult);
360 rh->fEvsN->Fill(mult,n);
361 rh->fEtaVsN->Fill(eta, n);
363 Double_t c = Correction(d,r,s,t,vtxbin,eta,lowFlux);
364 fCorrections->Fill(c);
366 rh->fEvsM->Fill(mult,n);
367 rh->fEtaVsM->Fill(eta, n);
368 rh->fCorr->Fill(eta, c);
370 Bool_t hit = (n > 0.9 && c > 0);
371 if (hit) rh->fELossUsed->Fill(mult);
372 rh->fPoisson.Fill(t,s,hit,1./c);
374 if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n);
378 TH2D* hclone = static_cast<TH2D*>(h->Clone("hclone"));
379 if (!fUsePoisson) hclone->Reset();
380 if ( fUsePoisson) h->Reset();
382 TH2D* poisson = rh->fPoisson.Result();
383 for (Int_t t=0; t < poisson->GetNbinsX(); t++) {
384 for (Int_t s=0; s < poisson->GetNbinsY(); s++) {
386 Double_t poissonV = poisson->GetBinContent(t+1,s+1);
387 Double_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
388 Double_t eta = fmd.Eta(d,r,s,t);
390 eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,zvtx);
392 h->Fill(eta,phi,poissonV);
394 hclone->Fill(eta,phi,poissonV);
395 if (fUsePoisson) rh->fDensity->Fill(eta, phi, poissonV);
399 for (Int_t ieta=1; ieta <= h->GetNbinsX(); ieta++) {
400 for (Int_t iphi=1; iphi<= h->GetNbinsY(); iphi++) {
402 Double_t poissonV = 0; //h->GetBinContent(,s+1);
405 poissonV = h->GetBinContent(ieta,iphi);
406 eLossV = hclone->GetBinContent(ieta,iphi);
409 poissonV = hclone->GetBinContent(ieta,iphi);
410 eLossV = h->GetBinContent(ieta,iphi);
413 rh->fELossVsPoisson->Fill(eLossV, poissonV);
424 //_____________________________________________________________________
426 AliFMDDensityCalculator::FindMaxWeight(const AliFMDCorrELossFit* cor,
427 UShort_t d, Char_t r, Int_t iEta) const
430 // Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
438 DGUARD(fDebug, 3, "Find maximum weight in FMD density calculator");
441 AliFMDCorrELossFit::ELossFit* fit = cor->GetFit(d,r,iEta);
443 // AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
446 return TMath::Min(Int_t(fMaxParticles), fit->FindMaxWeight());
449 //_____________________________________________________________________
451 AliFMDDensityCalculator::CacheMaxWeights(const TAxis& axis)
454 // Find the max weights and cache them
456 DGUARD(fDebug, 2, "Cache maximum weights in FMD density calculator");
457 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
458 AliFMDCorrELossFit* cor = fcm.GetELossFit();
460 TAxis eta(axis.GetNbins(),
465 eta.Set(cor->GetEtaAxis().GetNbins(),
466 cor->GetEtaAxis().GetXmin(),
467 cor->GetEtaAxis().GetXmax());
469 Int_t nEta = eta.GetNbins();
476 fMaxWeights->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
477 fMaxWeights->GetYaxis()->SetBinLabel(1, "FMD1i");
478 fMaxWeights->GetYaxis()->SetBinLabel(2, "FMD2i");
479 fMaxWeights->GetYaxis()->SetBinLabel(3, "FMD2o");
480 fMaxWeights->GetYaxis()->SetBinLabel(4, "FMD3i");
481 fMaxWeights->GetYaxis()->SetBinLabel(5, "FMD3o");
483 AliInfo(Form("Get eta axis with %d bins from %f to %f",
484 nEta, eta.GetXmin(), eta.GetXmax()));
485 fLowCuts->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
486 fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
487 fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
488 fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
489 fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
490 fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
492 for (Int_t i = 0; i < nEta; i++) {
494 w[0] = fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1);
495 w[1] = fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1);
496 w[2] = fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1);
497 w[3] = fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1);
498 w[4] = fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1);
500 l[0] = GetMultCut(1, 'I', i+1, false);
501 l[1] = GetMultCut(2, 'I', i+1, false);
502 l[2] = GetMultCut(2, 'O', i+1, false);
503 l[3] = GetMultCut(3, 'I', i+1, false);
504 l[4] = GetMultCut(3, 'O', i+1, false);
505 for (Int_t j = 0; j < 5; j++) {
506 if (w[j] > 0) fMaxWeights->SetBinContent(i+1, j+1, w[j]);
507 if (l[j] > 0) fLowCuts->SetBinContent(i+1, j+1, l[j]);
512 //_____________________________________________________________________
514 AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const
517 // Find the (cached) maximum weight for FMD<i>dr</i> in
518 // @f$\eta@f$ bin @a iEta
526 // max weight or <= 0 in case of problems
528 if (iEta < 0) return -1;
530 const TArrayI* max = 0;
532 case 1: max = &fFMD1iMax; break;
533 case 2: max = (r == 'I' || r == 'i' ? &fFMD2iMax : &fFMD2oMax); break;
534 case 3: max = (r == 'I' || r == 'i' ? &fFMD3iMax : &fFMD3oMax); break;
537 AliWarning(Form("No array for FMD%d%c", d, r));
541 if (iEta >= max->fN) {
542 AliWarning(Form("Eta bin %3d out of bounds [0,%d]",
547 AliDebug(30,Form("Max weight for FMD%d%c eta bin %3d: %d", d, r, iEta,
549 return max->At(iEta);
552 //_____________________________________________________________________
554 AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const
557 // Find the (cached) maximum weight for FMD<i>dr</i> iat
566 // max weight or <= 0 in case of problems
568 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
569 Int_t iEta = fcm.GetELossFit()->FindEtaBin(eta) -1;
571 return GetMaxWeight(d, r, iEta);
574 //_____________________________________________________________________
576 AliFMDDensityCalculator::NParticles(Float_t mult,
583 Bool_t lowFlux) const
586 // Get the number of particles corresponding to the signal mult
593 // t Strip (not used)
595 // eta Pseudo-rapidity
596 // lowFlux Low-flux flag
599 // The number of particles
601 // if (mult <= GetMultCut()) return 0;
602 DGUARD(fDebug, 3, "Calculate Nch in FMD density calculator");
603 if (lowFlux) return 1;
605 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
606 AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta);
608 AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
612 Int_t m = GetMaxWeight(d,r,eta); // fit->FindMaxWeight();
614 AliWarning(Form("No good fits for FMD%d%c at eta=%f", d, r, eta));
618 UShort_t n = TMath::Min(fMaxParticles, UShort_t(m));
619 Double_t ret = fit->EvaluateWeighted(mult, n);
622 AliInfo(Form("FMD%d%c, eta=%7.4f, %8.5f -> %8.5f", d, r, eta, mult, ret));
625 fWeightedSum->Fill(ret);
626 fSumOfWeights->Fill(ret);
631 //_____________________________________________________________________
633 AliFMDDensityCalculator::Correction(UShort_t d,
639 Bool_t lowFlux) const
642 // Get the inverse correction factor. This consist of
644 // - acceptance correction (corners of sensors)
645 // - double hit correction (for low-flux events)
646 // - dead strip correction
652 // t Strip (not used)
654 // eta Pseudo-rapidity
655 // lowFlux Low-flux flag
660 DGUARD(fDebug, 3, "Apply correction in FMD density calculator");
661 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
663 Float_t correction = 1;
664 if (fUsePhiAcceptance == kPhiCorrectNch)
665 correction = AcceptanceCorrection(r,t);
668 if (fcm.GetDoubleHit())
669 dblHitCor = fcm.GetDoubleHit()->GetCorrection(d,r);
672 Double_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
673 if (dblC > 0) correction *= dblC;
676 AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
682 //_____________________________________________________________________
684 AliFMDDensityCalculator::GenerateAcceptanceCorrection(Char_t r) const
687 // Generate the acceptance corrections
690 // r Ring to generate for
693 // Newly allocated histogram of acceptance corrections
695 DGUARD(fDebug, 3, "Make acceptance correction in FMD density calculator");
696 const Double_t ic1[] = { 4.9895, 15.3560 };
697 const Double_t ic2[] = { 1.8007, 17.2000 };
698 const Double_t oc1[] = { 4.2231, 26.6638 };
699 const Double_t oc2[] = { 1.8357, 27.9500 };
700 const Double_t* c1 = (r == 'I' || r == 'i' ? ic1 : oc1);
701 const Double_t* c2 = (r == 'I' || r == 'i' ? ic2 : oc2);
702 Double_t minR = (r == 'I' || r == 'i' ? 4.5213 : 15.4);
703 Double_t maxR = (r == 'I' || r == 'i' ? 17.2 : 28.0);
704 Int_t nStrips = (r == 'I' || r == 'i' ? 512 : 256);
705 Int_t nSec = (r == 'I' || r == 'i' ? 20 : 40);
706 Float_t basearc = 2 * TMath::Pi() / nSec;
707 Double_t rad = maxR - minR;
708 Float_t segment = rad / nStrips;
709 Float_t cr = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
711 // Numbers used to find end-point of strip.
712 // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
713 Float_t D = c1[0] * c2[1] - c1[1] * c2[0];
714 Float_t dx = c2[0] - c1[0];
715 Float_t dy = c2[1] - c1[1];
716 Float_t dr = TMath::Sqrt(dx*dx+dy*dy);
718 TH1D* ret = new TH1D(Form("acc%c", r),
719 Form("Acceptance correction for FMDx%c", r),
720 nStrips, -.5, nStrips-.5);
721 ret->SetXTitle("Strip");
722 ret->SetYTitle("#varphi acceptance");
723 ret->SetDirectory(0);
724 ret->SetFillColor(r == 'I' || r == 'i' ? kRed+1 : kBlue+1);
725 ret->SetFillStyle(3001);
727 for (Int_t t = 0; t < nStrips; t++) {
728 Float_t radius = minR + t * segment;
730 // If the radius of the strip is smaller than the radius corresponding
731 // to the first corner we have a full strip length
733 ret->SetBinContent(t+1, 1);
737 // Next, we should find the end-point of the strip - that is,
738 // the coordinates where the line from c1 to c2 intersects a circle
739 // with radius given by the strip.
740 // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
741 // Calculate the determinant
742 Float_t det = radius * radius * dr * dr - D*D;
745 // <0 means No intersection
746 // =0 means Exactly tangent
747 ret->SetBinContent(t+1, 1);
751 // Calculate end-point and the corresponding opening angle
752 Float_t x = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
753 Float_t y = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
754 Float_t th = TMath::ATan2(x, y);
756 ret->SetBinContent(t+1, th / basearc);
761 //_____________________________________________________________________
763 AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
766 // Get the acceptance correction for strip @a t in an ring of type @a r
769 // r Ring type ('I' or 'O')
773 // Inverse acceptance correction
775 TH1D* acc = (r == 'I' || r == 'i' ? fAccI : fAccO);
776 return acc->GetBinContent(t+1);
779 //____________________________________________________________________
781 AliFMDDensityCalculator::ScaleHistograms(const TList* dir, Int_t nEvents)
784 // Scale the histograms to the total number of events
787 // dir where to put the output
788 // nEvents Number of events
790 DGUARD(fDebug, 1, "Scale histograms in FMD density calculator");
791 if (nEvents <= 0) return;
792 TList* d = static_cast<TList*>(dir->FindObject(GetName()));
795 TIter next(&fRingHistos);
797 THStack* sums = new THStack("sums", "sums of ring signals");
798 while ((o = static_cast<RingHistos*>(next()))) {
799 o->ScaleHistograms(d, nEvents);
800 TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
801 o->fDensity->GetNbinsY(),"e");
802 sum->Scale(1., "width");
803 sum->SetTitle(o->GetName());
804 sum->SetDirectory(0);
805 sum->SetYTitle("#sum N_{ch,incl}");
811 //____________________________________________________________________
813 AliFMDDensityCalculator::DefineOutput(TList* dir)
816 // Output diagnostic histograms to directory
819 // dir List to write in
821 DGUARD(fDebug, 1, "Define output FMD density calculator");
822 TList* d = new TList;
824 d->SetName(GetName());
826 d->Add(fWeightedSum);
827 d->Add(fSumOfWeights);
828 d->Add(fCorrections);
834 // TNamed* sigma = new TNamed("sigma",
835 // (fIncludeSigma ? "included" : "excluded"));
836 TNamed* maxP = new TNamed("maxParticle", Form("%d", fMaxParticles));
837 TNamed* method = new TNamed("method",
838 (fUsePoisson ? "Poisson" : "Energy loss"));
839 TNamed* phiA = new TNamed("phiAcceptance",
840 (fUsePhiAcceptance == 0 ? "disabled" :
841 fUsePhiAcceptance == 1 ? "particles" :
843 TNamed* etaL = new TNamed("etaLumping", Form("%d", fEtaLumping));
844 TNamed* phiL = new TNamed("phiLumping", Form("%d", fPhiLumping));
845 // TParameter<double>* nxi = new TParameter<double>("nXi", fNXi);
856 TIter next(&fRingHistos);
858 while ((o = static_cast<RingHistos*>(next()))) {
859 o->fPoisson.SetLumping(fEtaLumping, fPhiLumping);
863 //____________________________________________________________________
865 AliFMDDensityCalculator::Print(Option_t* option) const
873 char ind[gROOT->GetDirLevel()+3];
874 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
875 ind[gROOT->GetDirLevel()] = '\0';
876 std::cout << ind << ClassName() << ": " << GetName() << '\n'
878 << ind << " Max(particles): " << fMaxParticles << '\n'
879 << ind << " Poisson method: " << fUsePoisson << '\n'
880 << ind << " Use phi acceptance: " << fUsePhiAcceptance << '\n'
881 << ind << " Eta lumping: " << fEtaLumping << '\n'
882 << ind << " Phi lumping: " << fPhiLumping << '\n'
885 std::cout << ind << " Lower cut:" << std::endl;
889 if (opt.Contains("nomax")) return;
891 std::cout << ind << " Max weights:\n";
893 for (UShort_t d=1; d<=3; d++) {
894 UShort_t nr = (d == 1 ? 1 : 2);
895 for (UShort_t q=0; q<nr; q++) {
896 ind[gROOT->GetDirLevel()] = ' ';
897 ind[gROOT->GetDirLevel()+1] = '\0';
898 Char_t r = (q == 0 ? 'I' : 'O');
899 std::cout << ind << " FMD" << d << r << ":";
900 ind[gROOT->GetDirLevel()+1] = ' ';
901 ind[gROOT->GetDirLevel()+2] = '\0';
903 const TArrayI& a = (d == 1 ? fFMD1iMax :
904 (d == 2 ? (r == 'I' ? fFMD2iMax : fFMD2oMax) :
905 (r == 'I' ? fFMD3iMax : fFMD3oMax)));
907 for (Int_t i = 0; i < a.fN; i++) {
908 if (a.fArray[i] < 1) continue;
909 if (j % 6 == 0) std::cout << "\n " << ind;
911 std::cout << " " << std::setw(3) << i << ": " << a.fArray[i];
913 std::cout << std::endl;
918 //====================================================================
919 AliFMDDensityCalculator::RingHistos::RingHistos()
920 : AliForwardUtil::RingHistos(),
938 //____________________________________________________________________
939 AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
940 : AliForwardUtil::RingHistos(d,r),
960 fEvsN = new TH2D("elossVsNnocorr",
961 "#Delta E/#Delta E_{mip} vs uncorrected inclusive N_{ch}",
962 250, -.5, 24.5, 250, -.5, 24.5);
963 fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
964 fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
966 fEvsN->SetDirectory(0);
968 fEvsM = static_cast<TH2D*>(fEvsN->Clone("elossVsNcorr"));
969 fEvsM->SetTitle("#Delta E/#Delta E_{mip} vs corrected inclusive N_{ch}");
970 fEvsM->SetDirectory(0);
972 fEtaVsN = new TProfile("etaVsNnocorr",
973 "Average inclusive N_{ch} vs #eta (uncorrected)",
975 fEtaVsN->SetXTitle("#eta");
976 fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
977 fEtaVsN->SetDirectory(0);
978 fEtaVsN->SetLineColor(Color());
979 fEtaVsN->SetFillColor(Color());
981 fEtaVsM = static_cast<TProfile*>(fEtaVsN->Clone("etaVsNcorr"));
982 fEtaVsM->SetTitle("Average inclusive N_{ch} vs #eta (corrected)");
983 fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
984 fEtaVsM->SetDirectory(0);
987 fCorr = new TProfile("corr", "Average correction", 200, -4, 6);
988 fCorr->SetXTitle("#eta");
989 fCorr->SetYTitle("#LT correction#GT");
990 fCorr->SetDirectory(0);
991 fCorr->SetLineColor(Color());
992 fCorr->SetFillColor(Color());
994 fDensity = new TH2D("inclDensity", "Inclusive N_{ch} density",
995 200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
997 fDensity->SetDirectory(0);
999 fDensity->SetMarkerColor(Color());
1000 fDensity->SetXTitle("#eta");
1001 fDensity->SetYTitle("#phi [radians]");
1002 fDensity->SetZTitle("Inclusive N_{ch} density");
1004 fELossVsPoisson = new TH2D("elossVsPoisson",
1005 "N_{ch} from energy loss vs from Poission",
1006 500, 0, 100, 500, 0, 100);
1007 fELossVsPoisson->SetDirectory(0);
1008 fELossVsPoisson->SetXTitle("N_{ch} from #DeltaE");
1009 fELossVsPoisson->SetYTitle("N_{ch} from Poisson");
1010 fELossVsPoisson->SetZTitle("Correlation");
1012 fELoss = new TH1D("eloss", "#Delta/#Delta_{mip} in all strips",
1014 fELoss->SetXTitle("#Delta/#Delta_{mip} (selected)");
1015 fELoss->SetYTitle("P(#Delta/#Delta_{mip})");
1016 fELoss->SetFillColor(Color()-2);
1017 fELoss->SetFillStyle(3003);
1018 fELoss->SetLineColor(kBlack);
1019 fELoss->SetLineStyle(2);
1020 fELoss->SetLineWidth(2);
1021 fELoss->SetDirectory(0);
1023 fELossUsed = static_cast<TH1D*>(fELoss->Clone("elossUsed"));
1024 fELossUsed->SetTitle("#Delta/#Delta_{mip} in used strips");
1025 fELossUsed->SetFillStyle(3002);
1026 fELossUsed->SetLineStyle(1);
1027 fELossUsed->SetDirectory(0);
1030 //____________________________________________________________________
1031 AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
1032 : AliForwardUtil::RingHistos(o),
1038 fDensity(o.fDensity),
1039 fELossVsPoisson(o.fELossVsPoisson),
1040 fPoisson(o.fPoisson),
1042 fELossUsed(o.fELossUsed),
1043 fMultCut(o.fMultCut)
1049 // o Object to copy from
1053 //____________________________________________________________________
1054 AliFMDDensityCalculator::RingHistos&
1055 AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
1058 // Assignment operator
1061 // o Object to assign from
1064 // Reference to this
1066 if (&o == this) return *this;
1067 AliForwardUtil::RingHistos::operator=(o);
1069 if (fEvsN) delete fEvsN;
1070 if (fEvsM) delete fEvsM;
1071 if (fEtaVsN) delete fEtaVsN;
1072 if (fEtaVsM) delete fEtaVsM;
1073 if (fCorr) delete fCorr;
1074 if (fDensity) delete fDensity;
1075 if (fELossVsPoisson) delete fELossVsPoisson;
1077 fEvsN = static_cast<TH2D*>(o.fEvsN->Clone());
1078 fEvsM = static_cast<TH2D*>(o.fEvsM->Clone());
1079 fEtaVsN = static_cast<TProfile*>(o.fEtaVsN->Clone());
1080 fEtaVsM = static_cast<TProfile*>(o.fEtaVsM->Clone());
1081 fCorr = static_cast<TProfile*>(o.fCorr->Clone());
1082 fDensity = static_cast<TH2D*>(o.fDensity->Clone());
1083 fELossVsPoisson = static_cast<TH2D*>(o.fELossVsPoisson->Clone());
1084 fPoisson = o.fPoisson;
1085 fELoss = static_cast<TH1D*>(o.fELoss->Clone());
1086 fELossUsed = static_cast<TH1D*>(o.fELossUsed->Clone());
1090 //____________________________________________________________________
1091 AliFMDDensityCalculator::RingHistos::~RingHistos()
1099 //____________________________________________________________________
1101 AliFMDDensityCalculator::RingHistos::Init(const TAxis& /*eAxis*/)
1103 fPoisson.Init(-1,-1);
1106 //____________________________________________________________________
1108 AliFMDDensityCalculator::RingHistos::Output(TList* dir)
1114 // dir Where to put it
1116 TList* d = DefineOutputList(dir);
1123 d->Add(fELossVsPoisson);
1125 fPoisson.GetOccupancy()->SetFillColor(Color());
1126 fPoisson.GetMean()->SetFillColor(Color());
1127 fPoisson.GetOccupancy()->SetFillColor(Color());
1131 Bool_t inner = (fRing == 'I' || fRing == 'i');
1132 Int_t nStr = inner ? 512 : 256;
1133 Int_t nSec = inner ? 20 : 40;
1134 TAxis x(nStr, -.5, nStr-.5);
1135 TAxis y(nSec, -.5, nSec-.5);
1136 x.SetTitle("strip");
1137 y.SetTitle("sector");
1138 fPoisson.Define(x, y);
1140 TParameter<double>* cut = new TParameter<double>("cut", fMultCut);
1144 //____________________________________________________________________
1146 AliFMDDensityCalculator::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
1149 // Scale the histograms to the total number of events
1152 // dir Where the output is
1153 // nEvents Number of events
1155 TList* l = GetOutputList(dir);
1158 TH1* density = GetOutputHist(l,"inclDensity");
1159 if (density) density->Scale(1./nEvents);
1162 //____________________________________________________________________