]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/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
Changes for Root6: removing obsolete TH1 functions, corrected EINCLUDE, additional...
[u/mrichter/AliRoot.git] / PWGLF / 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>
5ca83fee 17#include <TVector3.h>
1f7aa5c7 18#include <TStopwatch.h>
bfab35d9 19#include <TParameter.h>
0bd4b00f 20#include <iostream>
21#include <iomanip>
1f7aa5c7 22#include <cstring>
7e4038b5 23
24ClassImp(AliFMDDensityCalculator)
25#if 0
26; // For Emacs
27#endif
28
5ca83fee 29//____________________________________________________________________
30const char* AliFMDDensityCalculator::fgkFolderName = "fmdDensityCalculator";
31
7e4038b5 32//____________________________________________________________________
33AliFMDDensityCalculator::AliFMDDensityCalculator()
34 : TNamed(),
35 fRingHistos(),
dd497217 36 fSumOfWeights(0),
37 fWeightedSum(0),
ea3e5d95 38 fCorrections(0),
0bd4b00f 39 fMaxParticles(5),
9d05ffeb 40 fUsePoisson(false),
8411f7fe 41 fUsePhiAcceptance(kPhiCorrectNch),
0bd4b00f 42 fAccI(0),
43 fAccO(0),
1174780f 44 fFMD1iMax(0),
45 fFMD2iMax(0),
46 fFMD2oMax(0),
47 fFMD3iMax(0),
48 fFMD3oMax(0),
5bb5d1f6 49 fMaxWeights(0),
50 fLowCuts(0),
21d778b1 51 fEtaLumping(32),
52 fPhiLumping(4),
d2638bb7 53 fDebug(0),
6f4a5c0d 54 fCuts(),
5ca83fee 55 fRecalculateEta(false),
8449e3e0 56 fRecalculatePhi(false),
1f7aa5c7 57 fMinQuality(10),
58 fCache(),
59 fDoTiming(false),
60 fHTiming(0)
7984e5f7 61{
62 //
63 // Constructor
64 //
5ca83fee 65 DGUARD(fDebug, 3, "Default CTOR of FMD density calculator");
7984e5f7 66}
7e4038b5 67
68//____________________________________________________________________
69AliFMDDensityCalculator::AliFMDDensityCalculator(const char* title)
5ca83fee 70 : TNamed(fgkFolderName, title),
7e4038b5 71 fRingHistos(),
dd497217 72 fSumOfWeights(0),
73 fWeightedSum(0),
ea3e5d95 74 fCorrections(0),
0bd4b00f 75 fMaxParticles(5),
9d05ffeb 76 fUsePoisson(false),
8411f7fe 77 fUsePhiAcceptance(kPhiCorrectNch),
0bd4b00f 78 fAccI(0),
79 fAccO(0),
1174780f 80 fFMD1iMax(0),
81 fFMD2iMax(0),
82 fFMD2oMax(0),
83 fFMD3iMax(0),
84 fFMD3oMax(0),
5bb5d1f6 85 fMaxWeights(0),
86 fLowCuts(0),
21d778b1 87 fEtaLumping(32),
88 fPhiLumping(4),
d2638bb7 89 fDebug(0),
6f4a5c0d 90 fCuts(),
5ca83fee 91 fRecalculateEta(false),
8449e3e0 92 fRecalculatePhi(false),
1f7aa5c7 93 fMinQuality(10),
94 fCache(),
95 fDoTiming(false),
96 fHTiming(0)
7e4038b5 97{
7984e5f7 98 //
99 // Constructor
100 //
101 // Parameters:
102 // name Name of object
103 //
5ca83fee 104 DGUARD(fDebug, 3, "Named CTOR of FMD density calculator: %s", title);
7e4038b5 105 fRingHistos.SetName(GetName());
e308a636 106 fRingHistos.SetOwner();
7e4038b5 107 fRingHistos.Add(new RingHistos(1, 'I'));
108 fRingHistos.Add(new RingHistos(2, 'I'));
109 fRingHistos.Add(new RingHistos(2, 'O'));
110 fRingHistos.Add(new RingHistos(3, 'I'));
111 fRingHistos.Add(new RingHistos(3, 'O'));
dd497217 112 fSumOfWeights = new TH1D("sumOfWeights", "Sum of Landau weights",
113 200, 0, 20);
0bd4b00f 114 fSumOfWeights->SetFillColor(kRed+1);
115 fSumOfWeights->SetXTitle("#sum_{i} a_{i} f_{i}(#Delta)");
dd497217 116 fWeightedSum = new TH1D("weightedSum", "Weighted sum of Landau propability",
117 200, 0, 20);
0bd4b00f 118 fWeightedSum->SetFillColor(kBlue+1);
119 fWeightedSum->SetXTitle("#sum_{i} i a_{i} f_{i}(#Delta)");
dd497217 120 fCorrections = new TH1D("corrections", "Distribution of corrections",
121 100, 0, 10);
0bd4b00f 122 fCorrections->SetFillColor(kBlue+1);
123 fCorrections->SetXTitle("correction");
dd497217 124
0bd4b00f 125 fAccI = GenerateAcceptanceCorrection('I');
126 fAccO = GenerateAcceptanceCorrection('O');
5bb5d1f6 127
128 fMaxWeights = new TH2D("maxWeights", "Maximum i of a_{i}'s to use",
129 1, 0, 1, 1, 0, 1);
130 fMaxWeights->SetXTitle("#eta");
131 fMaxWeights->SetDirectory(0);
132
133 fLowCuts = new TH2D("lowCuts", "Low cuts used", 1, 0, 1, 1, 0, 1);
134 fLowCuts->SetXTitle("#eta");
135 fLowCuts->SetDirectory(0);
136
7e4038b5 137}
138
139//____________________________________________________________________
140AliFMDDensityCalculator::AliFMDDensityCalculator(const
141 AliFMDDensityCalculator& o)
142 : TNamed(o),
143 fRingHistos(),
dd497217 144 fSumOfWeights(o.fSumOfWeights),
145 fWeightedSum(o.fWeightedSum),
ea3e5d95 146 fCorrections(o.fCorrections),
0bd4b00f 147 fMaxParticles(o.fMaxParticles),
9d05ffeb 148 fUsePoisson(o.fUsePoisson),
5bb5d1f6 149 fUsePhiAcceptance(o.fUsePhiAcceptance),
0bd4b00f 150 fAccI(o.fAccI),
151 fAccO(o.fAccO),
1174780f 152 fFMD1iMax(o.fFMD1iMax),
153 fFMD2iMax(o.fFMD2iMax),
154 fFMD2oMax(o.fFMD2oMax),
155 fFMD3iMax(o.fFMD3iMax),
156 fFMD3oMax(o.fFMD3oMax),
5bb5d1f6 157 fMaxWeights(o.fMaxWeights),
158 fLowCuts(o.fLowCuts),
b6b35c77 159 fEtaLumping(o.fEtaLumping),
160 fPhiLumping(o.fPhiLumping),
d2638bb7 161 fDebug(o.fDebug),
6f4a5c0d 162 fCuts(o.fCuts),
5ca83fee 163 fRecalculateEta(o.fRecalculateEta),
8449e3e0 164 fRecalculatePhi(o.fRecalculatePhi),
1f7aa5c7 165 fMinQuality(o.fMinQuality),
166 fCache(o.fCache),
167 fDoTiming(o.fDoTiming),
168 fHTiming(o.fHTiming)
7e4038b5 169{
7984e5f7 170 //
171 // Copy constructor
172 //
173 // Parameters:
174 // o Object to copy from
175 //
5ca83fee 176 DGUARD(fDebug, 3, "Copy CTOR of FMD density calculator");
7e4038b5 177 TIter next(&o.fRingHistos);
178 TObject* obj = 0;
179 while ((obj = next())) fRingHistos.Add(obj);
180}
181
182//____________________________________________________________________
183AliFMDDensityCalculator::~AliFMDDensityCalculator()
184{
7984e5f7 185 //
186 // Destructor
187 //
6ab100ec 188 DGUARD(fDebug, 3, "DTOR of FMD density calculator");
36ffcf83 189 // fRingHistos.Delete();
7e4038b5 190}
191
192//____________________________________________________________________
193AliFMDDensityCalculator&
194AliFMDDensityCalculator::operator=(const AliFMDDensityCalculator& o)
195{
7984e5f7 196 //
197 // Assignement operator
198 //
199 // Parameters:
200 // o Object to assign from
201 //
202 // Return:
203 // Reference to this object
204 //
6ab100ec 205 DGUARD(fDebug, 3, "Assignment of FMD density calculator");
d015ecfe 206 if (&o == this) return *this;
ea3e5d95 207 TNamed::operator=(o);
7e4038b5 208
d2638bb7 209 fDebug = o.fDebug;
210 fMaxParticles = o.fMaxParticles;
211 fUsePoisson = o.fUsePoisson;
212 fUsePhiAcceptance = o.fUsePhiAcceptance;
213 fAccI = o.fAccI;
214 fAccO = o.fAccO;
215 fFMD1iMax = o.fFMD1iMax;
216 fFMD2iMax = o.fFMD2iMax;
217 fFMD2oMax = o.fFMD2oMax;
218 fFMD3iMax = o.fFMD3iMax;
219 fFMD3oMax = o.fFMD3oMax;
220 fMaxWeights = o.fMaxWeights;
221 fLowCuts = o.fLowCuts;
222 fEtaLumping = o.fEtaLumping;
223 fPhiLumping = o.fPhiLumping;
224 fCuts = o.fCuts;
6f4a5c0d 225 fRecalculateEta = o.fRecalculateEta;
5ca83fee 226 fRecalculatePhi = o.fRecalculatePhi;
8449e3e0 227 fMinQuality = o.fMinQuality;
1f7aa5c7 228 fCache = o.fCache;
229 fDoTiming = o.fDoTiming;
230 fHTiming = o.fHTiming;
d2638bb7 231
7e4038b5 232 fRingHistos.Delete();
233 TIter next(&o.fRingHistos);
234 TObject* obj = 0;
235 while ((obj = next())) fRingHistos.Add(obj);
236
237 return *this;
238}
239
1174780f 240//____________________________________________________________________
241void
5934a3e3 242AliFMDDensityCalculator::SetupForData(const TAxis& axis)
1174780f 243{
244 // Intialize this sub-algorithm
245 //
246 // Parameters:
5ca83fee 247 // etaAxis Eta axis
6ab100ec 248 DGUARD(fDebug, 1, "Initialize FMD density calculator");
e2ebf8c4 249 CacheMaxWeights(axis);
4077f3e8 250
1f7aa5c7 251 fCache.Init(axis);
252
9d05ffeb 253 TIter next(&fRingHistos);
254 RingHistos* o = 0;
d2638bb7 255 while ((o = static_cast<RingHistos*>(next()))) {
5934a3e3 256 o->SetupForData(axis);
e18cb8bd 257 // o->fMultCut = fCuts.GetFixedCut(o->fDet, o->fRing);
258 // o->fPoisson.Init(o->fDet,o->fRing,fEtaLumping, fPhiLumping);
d2638bb7 259 }
1174780f 260}
261
7e4038b5 262//____________________________________________________________________
263AliFMDDensityCalculator::RingHistos*
264AliFMDDensityCalculator::GetRingHistos(UShort_t d, Char_t r) const
265{
7984e5f7 266 //
267 // Get the ring histogram container
268 //
269 // Parameters:
270 // d Detector
271 // r Ring
272 //
273 // Return:
274 // Ring histogram container
275 //
7e4038b5 276 Int_t idx = -1;
277 switch (d) {
278 case 1: idx = 0; break;
279 case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
280 case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
281 }
e308a636 282 if (idx < 0 || idx >= fRingHistos.GetEntries()) {
283 AliWarning(Form("Index %d of FMD%d%c out of range", idx, d, r));
284 return 0;
285 }
7e4038b5 286
287 return static_cast<RingHistos*>(fRingHistos.At(idx));
288}
5bb5d1f6 289
0bd4b00f 290//____________________________________________________________________
291Double_t
5bb5d1f6 292AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Int_t ieta,
293 Bool_t errors) const
0bd4b00f 294{
7984e5f7 295 //
296 // Get the multiplicity cut. If the user has set fMultCut (via
297 // SetMultCut) then that value is used. If not, then the lower
298 // value of the fit range for the enery loss fits is returned.
299 //
300 // Return:
301 // Lower cut on multiplicity
302 //
d2638bb7 303 return fCuts.GetMultCut(d,r,ieta,errors);
5bb5d1f6 304}
305
306//____________________________________________________________________
307Double_t
308AliFMDDensityCalculator::GetMultCut(UShort_t d, Char_t r, Double_t eta,
8e400b14 309 Bool_t errors) const
5bb5d1f6 310{
311 //
312 // Get the multiplicity cut. If the user has set fMultCut (via
313 // SetMultCut) then that value is used. If not, then the lower
314 // value of the fit range for the enery loss fits is returned.
315 //
316 // Return:
317 // Lower cut on multiplicity
318 //
d2638bb7 319 return fCuts.GetMultCut(d,r,eta,errors);
0bd4b00f 320}
321
7e4038b5 322//____________________________________________________________________
323Bool_t
324AliFMDDensityCalculator::Calculate(const AliESDFMD& fmd,
325 AliForwardUtil::Histos& hists,
21d778b1 326 Bool_t lowFlux,
b3d81bee 327 Double_t /*cent*/,
5ca83fee 328 const TVector3& ip)
7e4038b5 329{
7984e5f7 330 //
331 // Do the calculations
332 //
333 // Parameters:
334 // fmd AliESDFMD object (possibly) corrected for sharing
335 // hists Histogram cache
336 // vtxBin Vertex bin
337 // lowFlux Low flux flag.
338 //
339 // Return:
340 // true on successs
6ab100ec 341 DGUARD(fDebug, 1, "Calculate density in FMD density calculator");
5ca83fee 342
1f7aa5c7 343 TStopwatch timer;
344 TStopwatch totalT;
345
346 // First measurements of timing
347 // Re-calculation : fraction of sum 32.0% of total 18.1%
348 // N_{particle} : fraction of sum 15.2% of total 8.6%
349 // Correction : fraction of sum 26.4% of total 14.9%
350 // #phi acceptance : fraction of sum 0.2% of total 0.1%
351 // Copy to cache : fraction of sum 3.9% of total 2.2%
352 // Poisson calculation : fraction of sum 18.7% of total 10.6%
353 // Diagnostics : fraction of sum 3.7% of total 2.1%
354 Double_t reEtaTime = 0;
355 Double_t nPartTime = 0;
356 Double_t corrTime = 0;
357 Double_t rePhiTime = 0;
358 Double_t copyTime = 0;
359 Double_t poissonTime= 0;
360 Double_t diagTime = 0;
361 if (fDoTiming) totalT.Start(true);
362
363 Double_t etaCache[20*512]; // Same number of strips per ring
364 Double_t phiCache[20*512]; // whether it is inner our outer.
365 // We do not use TArrayD because we do not wont a bounds check
366 // TArrayD etaCache(20*512); // Same number of strips per ring
367 // TArrayD phiCache(20*512); // whether it is inner our outer.
368
5ca83fee 369 // --- Loop over detectors -----------------------------------------
7e4038b5 370 for (UShort_t d=1; d<=3; d++) {
371 UShort_t nr = (d == 1 ? 1 : 2);
372 for (UShort_t q=0; q<nr; q++) {
373 Char_t r = (q == 0 ? 'I' : 'O');
374 UShort_t ns= (q == 0 ? 20 : 40);
375 UShort_t nt= (q == 0 ? 512 : 256);
376 TH2D* h = hists.Get(d,r);
377 RingHistos* rh= GetRingHistos(d,r);
e308a636 378 if (!rh) {
379 AliError(Form("No ring histogram found for FMD%d%c", d, r));
380 fRingHistos.ls();
381 return false;
382 }
e18cb8bd 383 // rh->fPoisson.SetObject(d,r,vtxbin,cent);
d23503ee 384 rh->fPoisson.Reset(0);
5ca83fee 385 rh->fTotal->Reset();
386 rh->fGood->Reset();
821ffd28 387 // rh->ResetPoissonHistos(h, fEtaLumping, fPhiLumping);
1f7aa5c7 388
389 // Reset our eta cache
390 // for (Int_t i = 0; i < 20*512; i++)
391 // etaCache[i] = phiCache[i] = AliESDFMD::kInvalidEta;
392 memset(etaCache, 0, sizeof(Double_t)*20*512);
393 memset(phiCache, 0, sizeof(Double_t)*20*512);
394 // etaCache.Reset(AliESDFMD::kInvalidEta);
395 // phiCache.Reset(AliESDFMD::kInvalidEta);
396
5ca83fee 397 // --- Loop over sectors and strips ----------------------------
7e4038b5 398 for (UShort_t s=0; s<ns; s++) {
399 for (UShort_t t=0; t<nt; t++) {
3e4a0875 400
5ca83fee 401 Float_t mult = fmd.Multiplicity(d,r,s,t);
402 Float_t phi = fmd.Phi(d,r,s,t) * TMath::DegToRad();
403 Float_t eta = fmd.Eta(d,r,s,t);
404 Double_t oldPhi = phi;
405
406 // --- Re-calculate eta - needed for satelittes ------------
1f7aa5c7 407 if (fDoTiming) timer.Start(true);
408 if (eta == AliESDFMD::kInvalidEta || fRecalculateEta)
5ca83fee 409 eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,ip.Z());
1f7aa5c7 410 if (fDoTiming) reEtaTime += timer.CpuTime();
411 etaCache[s*nt+t] = eta;
5ca83fee 412
413 // --- Check this strip ------------------------------------
414 rh->fTotal->Fill(eta);
bfab35d9 415 if (mult == AliESDFMD::kInvalidMult) { // || mult > 20) {
8449e3e0 416 // Do not count invalid stuff
417 rh->fELoss->Fill(-1);
418 // rh->fEvsN->Fill(mult,-1);
419 // rh->fEvsM->Fill(mult,-1);
33b078ef 420 continue;
421 }
bfab35d9 422 if (mult > 20)
423 AliWarningF("Raw multiplicity of FMD%d%c[%02d,%03d] = %f > 20",
424 d, r, s, t, mult);
5ca83fee 425 // --- Automatic calculation of acceptance -----------------
426 rh->fGood->Fill(eta);
427
428 // --- If we asked to re-calculate phi for (x,y) IP --------
1f7aa5c7 429 if (fDoTiming) timer.Start(true);
5ca83fee 430 if (fRecalculatePhi) {
431 oldPhi = phi;
432 phi = AliForwardUtil::GetPhiFromStrip(r, t, phi, ip.X(), ip.Y());
433 }
1f7aa5c7 434 phiCache[s*nt+t] = phi;
435 if (fDoTiming) rePhiTime += timer.CpuTime();
5ca83fee 436
437 // --- Apply phi corner correction to eloss ----------------
8411f7fe 438 if (fUsePhiAcceptance == kPhiCorrectELoss)
0082a8fc 439 mult *= AcceptanceCorrection(r,t);
5bb5d1f6 440
5ca83fee 441 // --- Get the low multiplicity cut ------------------------
8e400b14 442 Double_t cut = 1024;
443 if (eta != AliESDFMD::kInvalidEta) cut = GetMultCut(d, r, eta,false);
bfab35d9 444 else AliWarningF("Eta for FMD%d%c[%02d,%03d] is invalid: %f",
445 d, r, s, t, eta);
8e400b14 446
5ca83fee 447 // --- Now caluculate Nch for this strip using fits --------
1f7aa5c7 448 if (fDoTiming) timer.Start(true);
5bb5d1f6 449 Double_t n = 0;
5ca83fee 450 if (cut > 0 && mult > cut) n = NParticles(mult,d,r,eta,lowFlux);
f7cfc454 451 rh->fELoss->Fill(mult);
8449e3e0 452 // rh->fEvsN->Fill(mult,n);
453 // rh->fEtaVsN->Fill(eta, n);
1f7aa5c7 454 if (fDoTiming) nPartTime += timer.CpuTime();
9fde7142 455
5ca83fee 456 // --- Calculate correction if needed ----------------------
1f7aa5c7 457 if (fDoTiming) timer.Start(true);
458 // Temporary stuff - remove Correction call
459 Double_t c = 1;
460 if (fUsePhiAcceptance == kPhiCorrectNch)
461 c = AcceptanceCorrection(r,t);
462 // Double_t c = Correction(d,r,t,eta,lowFlux);
463 if (fDoTiming) corrTime += timer.CpuTime();
dd497217 464 fCorrections->Fill(c);
7e4038b5 465 if (c > 0) n /= c;
8449e3e0 466 // rh->fEvsM->Fill(mult,n);
467 // rh->fEtaVsM->Fill(eta, n);
0bd4b00f 468 rh->fCorr->Fill(eta, c);
f7cfc454 469
5ca83fee 470 // --- Accumulate Poisson statistics -----------------------
821ffd28 471 Bool_t hit = (n > 0.9 && c > 0);
5ca83fee 472 if (hit) {
473 rh->fELossUsed->Fill(mult);
474 if (fRecalculatePhi) {
475 rh->fPhiBefore->Fill(oldPhi);
476 rh->fPhiAfter->Fill(phi);
477 }
478 }
21d778b1 479 rh->fPoisson.Fill(t,s,hit,1./c);
7e4038b5 480 h->Fill(eta,phi,n);
5ca83fee 481
482 // --- If we use ELoss fits, apply now ---------------------
5bb5d1f6 483 if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n);
7e4038b5 484 } // for t
485 } // for s
5ca83fee 486
487 // --- Automatic acceptance - Calculate as an efficiency -------
1f7aa5c7 488 // This is very fast, so we do not bother to time it
5ca83fee 489 rh->fGood->Divide(rh->fGood, rh->fTotal, 1, 1, "B");
490
491 // --- Make a copy and reset as needed -------------------------
1f7aa5c7 492 if (fDoTiming) timer.Start(true);
493 TH2D* hclone = fCache.Get(d,r);
494 // hclone->Reset();
495 // TH2D* hclone = static_cast<TH2D*>(h->Clone("hclone"));
21d778b1 496 if (!fUsePoisson) hclone->Reset();
1f7aa5c7 497 else {
498 for (Int_t i = 0; i <= h->GetNbinsX()+1; i++) {
499 for (Int_t j = 0; j <= h->GetNbinsY()+1; j++) {
500 hclone->SetBinContent(i,j,h->GetBinContent(i,j));
501 hclone->SetBinError(i,j,h->GetBinError(i,j));
502 }
503 }
504 // hclone->Add(h);
505 h->Reset();
506 }
507 if (fDoTiming) copyTime += timer.CpuTime();
21d778b1 508
5ca83fee 509 // --- Store Poisson result ------------------------------------
1f7aa5c7 510 if (fDoTiming) timer.Start(true);
821ffd28 511 TH2D* poisson = rh->fPoisson.Result();
d23503ee 512 for (Int_t t=0; t < poisson->GetNbinsX(); t++) {
513 for (Int_t s=0; s < poisson->GetNbinsY(); s++) {
21d778b1 514
515 Double_t poissonV = poisson->GetBinContent(t+1,s+1);
1f7aa5c7 516 // Use cached eta - since the calls to GetEtaFromStrip and
517 // GetPhiFromStrip are _very_ expensive
518 Double_t phi = phiCache[s*nt+t];
519 Double_t eta = etaCache[s*nt+t];
520 // Double_t phi = fmd.Phi(d,r,s,t) * TMath::DegToRad();
521 // Double_t eta = fmd.Eta(d,r,s,t);
522 // if (fRecalculateEta)
523 // eta = AliForwardUtil::GetEtaFromStrip(d,r,s,t,ip.Z());
524 // if (fRecalculatePhi)
525 // phi = AliForwardUtil::GetPhiFromStrip(r, t, phi, ip.X(), ip.Y());
526 if (fUsePoisson) {
21d778b1 527 h->Fill(eta,phi,poissonV);
1f7aa5c7 528 rh->fDensity->Fill(eta, phi, poissonV);
529 }
21d778b1 530 else
531 hclone->Fill(eta,phi,poissonV);
21d778b1 532 }
533 }
1f7aa5c7 534 if (fDoTiming) poissonTime += timer.CpuTime();
21d778b1 535
5ca83fee 536 // --- Make diagnostics - eloss vs poisson ---------------------
1f7aa5c7 537 if (fDoTiming) timer.Start(true);
5ca83fee 538 Int_t nY = h->GetNbinsY();
21d778b1 539 for (Int_t ieta=1; ieta <= h->GetNbinsX(); ieta++) {
5ca83fee 540 // Set the overflow bin to contain the phi acceptance
541 Double_t phiAcc = rh->fGood->GetBinContent(ieta);
542 Double_t phiAccE = rh->fGood->GetBinError(ieta);
543 h->SetBinContent(ieta, nY+1, phiAcc);
544 h->SetBinError(ieta, nY+1, phiAccE);
545 Double_t eta = h->GetXaxis()->GetBinCenter(ieta);
546 rh->fPhiAcc->Fill(eta, ip.Z(), phiAcc);
5ca83fee 547 for (Int_t iphi=1; iphi<= nY; iphi++) {
e83d0620 548
21d778b1 549 Double_t poissonV = 0; //h->GetBinContent(,s+1);
550 Double_t eLossV = 0;
551 if(fUsePoisson) {
552 poissonV = h->GetBinContent(ieta,iphi);
553 eLossV = hclone->GetBinContent(ieta,iphi);
554 }
555 else {
556 poissonV = hclone->GetBinContent(ieta,iphi);
557 eLossV = h->GetBinContent(ieta,iphi);
558 }
7b03929c 559
9d05ffeb 560 rh->fELossVsPoisson->Fill(eLossV, poissonV);
5ca83fee 561 rh->fDiffELossPoisson->Fill(poissonV < 1e-12 ? 0 :
562 (eLossV - poissonV) / poissonV);
563
9d05ffeb 564 }
565 }
1f7aa5c7 566 if (fDoTiming) diagTime += timer.CpuTime();
567 // delete hclone;
21d778b1 568
7e4038b5 569 } // for q
570 } // for d
a7f31ad8 571
1f7aa5c7 572 if (fDoTiming) {
573 fHTiming->Fill(1,reEtaTime);
574 fHTiming->Fill(2,nPartTime);
575 fHTiming->Fill(3,corrTime);
576 fHTiming->Fill(4,rePhiTime);
577 fHTiming->Fill(5,copyTime);
578 fHTiming->Fill(6,poissonTime);
579 fHTiming->Fill(7,diagTime);
580 fHTiming->Fill(8,totalT.CpuTime());
581 }
582
7e4038b5 583 return kTRUE;
584}
585
1174780f 586//_____________________________________________________________________
587Int_t
fb3430ac 588AliFMDDensityCalculator::FindMaxWeight(const AliFMDCorrELossFit* cor,
1174780f 589 UShort_t d, Char_t r, Int_t iEta) const
590{
fb3430ac 591 //
592 // Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
593 //
594 // Parameters:
595 // cor Correction
596 // d Detector
597 // r Ring
598 // iEta Eta bin
599 //
241cca4d 600 DGUARD(fDebug, 10, "Find maximum weight in FMD density calculator");
f61b1745 601 if(!cor) return -1;
602
8449e3e0 603 AliFMDCorrELossFit::ELossFit* fit = cor->FindFit(d,r,iEta, -1);
1174780f 604 if (!fit) {
605 // AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
606 return -1;
607 }
8449e3e0 608 return fit->FindMaxWeight(2*AliFMDCorrELossFit::ELossFit::fgMaxRelError,
609 AliFMDCorrELossFit::ELossFit::fgLeastWeight,
610 fMaxParticles);
1174780f 611}
612
613//_____________________________________________________________________
614void
e2ebf8c4 615AliFMDDensityCalculator::CacheMaxWeights(const TAxis& axis)
1174780f 616{
fb3430ac 617 //
618 // Find the max weights and cache them
619 //
6ab100ec 620 DGUARD(fDebug, 2, "Cache maximum weights in FMD density calculator");
1174780f 621 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
8449e3e0 622 const AliFMDCorrELossFit* cor = fcm.GetELossFit();
623 cor->CacheBins(fMinQuality);
a76fb27d 624
625 TAxis eta(axis.GetNbins(),
626 axis.GetXmin(),
627 axis.GetXmax());
628
f61b1745 629 if(cor)
630 eta.Set(cor->GetEtaAxis().GetNbins(),
631 cor->GetEtaAxis().GetXmin(),
632 cor->GetEtaAxis().GetXmax());
1174780f 633
634 Int_t nEta = eta.GetNbins();
635 fFMD1iMax.Set(nEta);
636 fFMD2iMax.Set(nEta);
637 fFMD2oMax.Set(nEta);
638 fFMD3iMax.Set(nEta);
639 fFMD3oMax.Set(nEta);
640
5bb5d1f6 641 fMaxWeights->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
642 fMaxWeights->GetYaxis()->SetBinLabel(1, "FMD1i");
643 fMaxWeights->GetYaxis()->SetBinLabel(2, "FMD2i");
644 fMaxWeights->GetYaxis()->SetBinLabel(3, "FMD2o");
645 fMaxWeights->GetYaxis()->SetBinLabel(4, "FMD3i");
646 fMaxWeights->GetYaxis()->SetBinLabel(5, "FMD3o");
647
648 AliInfo(Form("Get eta axis with %d bins from %f to %f",
649 nEta, eta.GetXmin(), eta.GetXmax()));
650 fLowCuts->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
651 fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
652 fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
653 fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
654 fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
655 fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
656
1174780f 657 for (Int_t i = 0; i < nEta; i++) {
5bb5d1f6 658 Double_t w[5];
659 w[0] = fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', i+1);
660 w[1] = fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', i+1);
661 w[2] = fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', i+1);
662 w[3] = fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', i+1);
663 w[4] = fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', i+1);
664 Double_t l[5];
665 l[0] = GetMultCut(1, 'I', i+1, false);
666 l[1] = GetMultCut(2, 'I', i+1, false);
667 l[2] = GetMultCut(2, 'O', i+1, false);
668 l[3] = GetMultCut(3, 'I', i+1, false);
669 l[4] = GetMultCut(3, 'O', i+1, false);
670 for (Int_t j = 0; j < 5; j++) {
671 if (w[j] > 0) fMaxWeights->SetBinContent(i+1, j+1, w[j]);
672 if (l[j] > 0) fLowCuts->SetBinContent(i+1, j+1, l[j]);
673 }
1174780f 674 }
675}
676
677//_____________________________________________________________________
678Int_t
679AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const
680{
fb3430ac 681 //
682 // Find the (cached) maximum weight for FMD<i>dr</i> in
683 // @f$\eta@f$ bin @a iEta
684 //
685 // Parameters:
686 // d Detector
687 // r Ring
688 // iEta Eta bin
689 //
690 // Return:
691 // max weight or <= 0 in case of problems
692 //
1174780f 693 if (iEta < 0) return -1;
694
695 const TArrayI* max = 0;
696 switch (d) {
697 case 1: max = &fFMD1iMax; break;
698 case 2: max = (r == 'I' || r == 'i' ? &fFMD2iMax : &fFMD2oMax); break;
699 case 3: max = (r == 'I' || r == 'i' ? &fFMD3iMax : &fFMD3oMax); break;
700 }
701 if (!max) {
702 AliWarning(Form("No array for FMD%d%c", d, r));
703 return -1;
704 }
705
706 if (iEta >= max->fN) {
707 AliWarning(Form("Eta bin %3d out of bounds [0,%d]",
708 iEta, max->fN-1));
709 return -1;
710 }
711
712 AliDebug(30,Form("Max weight for FMD%d%c eta bin %3d: %d", d, r, iEta,
713 max->At(iEta)));
714 return max->At(iEta);
715}
716
717//_____________________________________________________________________
718Int_t
719AliFMDDensityCalculator::GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const
720{
fb3430ac 721 //
722 // Find the (cached) maximum weight for FMD<i>dr</i> iat
723 // @f$\eta@f$
724 //
725 // Parameters:
726 // d Detector
727 // r Ring
728 // eta Eta bin
729 //
730 // Return:
731 // max weight or <= 0 in case of problems
732 //
1174780f 733 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
734 Int_t iEta = fcm.GetELossFit()->FindEtaBin(eta) -1;
735
736 return GetMaxWeight(d, r, iEta);
737}
738
7e4038b5 739//_____________________________________________________________________
740Float_t
741AliFMDDensityCalculator::NParticles(Float_t mult,
742 UShort_t d,
743 Char_t r,
7e4038b5 744 Float_t eta,
745 Bool_t lowFlux) const
746{
7984e5f7 747 //
748 // Get the number of particles corresponding to the signal mult
749 //
750 // Parameters:
751 // mult Signal
752 // d Detector
753 // r Ring
754 // s Sector
755 // t Strip (not used)
756 // v Vertex bin
757 // eta Pseudo-rapidity
758 // lowFlux Low-flux flag
759 //
760 // Return:
761 // The number of particles
762 //
5bb5d1f6 763 // if (mult <= GetMultCut()) return 0;
6ab100ec 764 DGUARD(fDebug, 3, "Calculate Nch in FMD density calculator");
7e4038b5 765 if (lowFlux) return 1;
0bd4b00f 766
0bd4b00f 767 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
8449e3e0 768 AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta, -1);
0bd4b00f 769 if (!fit) {
8449e3e0 770 AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f qual=%d",
771 d, r, eta, fMinQuality));
0bd4b00f 772 return 0;
773 }
774
1174780f 775 Int_t m = GetMaxWeight(d,r,eta); // fit->FindMaxWeight();
0bd4b00f 776 if (m < 1) {
777 AliWarning(Form("No good fits for FMD%d%c at eta=%f", d, r, eta));
778 return 0;
779 }
79909b8b 780
0bd4b00f 781 UShort_t n = TMath::Min(fMaxParticles, UShort_t(m));
782 Double_t ret = fit->EvaluateWeighted(mult, n);
79909b8b 783
0bd4b00f 784 if (fDebug > 10) {
785 AliInfo(Form("FMD%d%c, eta=%7.4f, %8.5f -> %8.5f", d, r, eta, mult, ret));
786 }
79909b8b 787
0bd4b00f 788 fWeightedSum->Fill(ret);
789 fSumOfWeights->Fill(ret);
3e4a0875 790
0bd4b00f 791 return ret;
7e4038b5 792}
793
794//_____________________________________________________________________
795Float_t
0bd4b00f 796AliFMDDensityCalculator::Correction(UShort_t d,
797 Char_t r,
0bd4b00f 798 UShort_t t,
0bd4b00f 799 Float_t eta,
800 Bool_t lowFlux) const
7e4038b5 801{
7984e5f7 802 //
803 // Get the inverse correction factor. This consist of
804 //
805 // - acceptance correction (corners of sensors)
806 // - double hit correction (for low-flux events)
807 // - dead strip correction
808 //
809 // Parameters:
810 // d Detector
811 // r Ring
812 // s Sector
813 // t Strip (not used)
814 // v Vertex bin
815 // eta Pseudo-rapidity
816 // lowFlux Low-flux flag
817 //
818 // Return:
819 //
820 //
241cca4d 821 DGUARD(fDebug, 10, "Apply correction in FMD density calculator");
5bb5d1f6 822 Float_t correction = 1;
8411f7fe 823 if (fUsePhiAcceptance == kPhiCorrectNch)
824 correction = AcceptanceCorrection(r,t);
7e4038b5 825 if (lowFlux) {
1f7aa5c7 826 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
827
0bd4b00f 828 TH1D* dblHitCor = 0;
829 if (fcm.GetDoubleHit())
830 dblHitCor = fcm.GetDoubleHit()->GetCorrection(d,r);
831
7e4038b5 832 if (dblHitCor) {
0bd4b00f 833 Double_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
834 if (dblC > 0) correction *= dblC;
7e4038b5 835 }
57522224 836 // else {
837 // AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
838 // }
7e4038b5 839 }
7e4038b5 840 return correction;
841}
842
0bd4b00f 843//_____________________________________________________________________
844TH1D*
845AliFMDDensityCalculator::GenerateAcceptanceCorrection(Char_t r) const
846{
7984e5f7 847 //
848 // Generate the acceptance corrections
849 //
850 // Parameters:
851 // r Ring to generate for
852 //
853 // Return:
854 // Newly allocated histogram of acceptance corrections
855 //
6ab100ec 856 DGUARD(fDebug, 3, "Make acceptance correction in FMD density calculator");
0bd4b00f 857 const Double_t ic1[] = { 4.9895, 15.3560 };
858 const Double_t ic2[] = { 1.8007, 17.2000 };
859 const Double_t oc1[] = { 4.2231, 26.6638 };
860 const Double_t oc2[] = { 1.8357, 27.9500 };
861 const Double_t* c1 = (r == 'I' || r == 'i' ? ic1 : oc1);
862 const Double_t* c2 = (r == 'I' || r == 'i' ? ic2 : oc2);
863 Double_t minR = (r == 'I' || r == 'i' ? 4.5213 : 15.4);
864 Double_t maxR = (r == 'I' || r == 'i' ? 17.2 : 28.0);
865 Int_t nStrips = (r == 'I' || r == 'i' ? 512 : 256);
866 Int_t nSec = (r == 'I' || r == 'i' ? 20 : 40);
867 Float_t basearc = 2 * TMath::Pi() / nSec;
868 Double_t rad = maxR - minR;
869 Float_t segment = rad / nStrips;
870 Float_t cr = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
871
872 // Numbers used to find end-point of strip.
873 // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
874 Float_t D = c1[0] * c2[1] - c1[1] * c2[0];
875 Float_t dx = c2[0] - c1[0];
876 Float_t dy = c2[1] - c1[1];
877 Float_t dr = TMath::Sqrt(dx*dx+dy*dy);
878
879 TH1D* ret = new TH1D(Form("acc%c", r),
880 Form("Acceptance correction for FMDx%c", r),
881 nStrips, -.5, nStrips-.5);
882 ret->SetXTitle("Strip");
883 ret->SetYTitle("#varphi acceptance");
884 ret->SetDirectory(0);
885 ret->SetFillColor(r == 'I' || r == 'i' ? kRed+1 : kBlue+1);
886 ret->SetFillStyle(3001);
887
888 for (Int_t t = 0; t < nStrips; t++) {
889 Float_t radius = minR + t * segment;
890
891 // If the radius of the strip is smaller than the radius corresponding
892 // to the first corner we have a full strip length
893 if (radius <= cr) {
894 ret->SetBinContent(t+1, 1);
895 continue;
896 }
897
898 // Next, we should find the end-point of the strip - that is,
899 // the coordinates where the line from c1 to c2 intersects a circle
900 // with radius given by the strip.
901 // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
902 // Calculate the determinant
903 Float_t det = radius * radius * dr * dr - D*D;
904
905 if (det <= 0) {
906 // <0 means No intersection
907 // =0 means Exactly tangent
908 ret->SetBinContent(t+1, 1);
909 continue;
910 }
911
912 // Calculate end-point and the corresponding opening angle
913 Float_t x = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
914 Float_t y = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
915 Float_t th = TMath::ATan2(x, y);
916
917 ret->SetBinContent(t+1, th / basearc);
918 }
919 return ret;
920}
7e4038b5 921
922//_____________________________________________________________________
923Float_t
924AliFMDDensityCalculator::AcceptanceCorrection(Char_t r, UShort_t t) const
925{
7984e5f7 926 //
927 // Get the acceptance correction for strip @a t in an ring of type @a r
928 //
929 // Parameters:
930 // r Ring type ('I' or 'O')
931 // t Strip number
932 //
933 // Return:
934 // Inverse acceptance correction
935 //
0bd4b00f 936 TH1D* acc = (r == 'I' || r == 'i' ? fAccI : fAccO);
937 return acc->GetBinContent(t+1);
7e4038b5 938}
939
940//____________________________________________________________________
941void
5934a3e3 942AliFMDDensityCalculator::Terminate(const TList* dir, TList* output,
943 Int_t nEvents)
7e4038b5 944{
7984e5f7 945 //
946 // Scale the histograms to the total number of events
947 //
948 // Parameters:
949 // dir where to put the output
950 // nEvents Number of events
951 //
6ab100ec 952 DGUARD(fDebug, 1, "Scale histograms in FMD density calculator");
7e4038b5 953 if (nEvents <= 0) return;
9d99b0dd 954 TList* d = static_cast<TList*>(dir->FindObject(GetName()));
955 if (!d) return;
7e4038b5 956
5934a3e3 957 TList* out = new TList;
958 out->SetName(d->GetName());
959 out->SetOwner();
960
7e4038b5 961 TIter next(&fRingHistos);
962 RingHistos* o = 0;
5bb5d1f6 963 THStack* sums = new THStack("sums", "sums of ring signals");
964 while ((o = static_cast<RingHistos*>(next()))) {
5934a3e3 965 o->Terminate(d, nEvents);
966 if (!o->fDensity) {
967 Warning("Terminate", "No density in %s", o->GetName());
968 continue;
969 }
5bb5d1f6 970 TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
971 o->fDensity->GetNbinsY(),"e");
972 sum->Scale(1., "width");
973 sum->SetTitle(o->GetName());
974 sum->SetDirectory(0);
975 sum->SetYTitle("#sum N_{ch,incl}");
976 sums->Add(sum);
977 }
5934a3e3 978 out->Add(sums);
979 output->Add(out);
7e4038b5 980}
981
982//____________________________________________________________________
983void
5934a3e3 984AliFMDDensityCalculator::CreateOutputObjects(TList* dir)
7e4038b5 985{
7984e5f7 986 //
987 // Output diagnostic histograms to directory
988 //
989 // Parameters:
990 // dir List to write in
991 //
6ab100ec 992 DGUARD(fDebug, 1, "Define output FMD density calculator");
7e4038b5 993 TList* d = new TList;
5bb5d1f6 994 d->SetOwner();
7e4038b5 995 d->SetName(GetName());
996 dir->Add(d);
dd497217 997 d->Add(fWeightedSum);
998 d->Add(fSumOfWeights);
999 d->Add(fCorrections);
0bd4b00f 1000 d->Add(fAccI);
1001 d->Add(fAccO);
5bb5d1f6 1002 d->Add(fMaxWeights);
1003 d->Add(fLowCuts);
0bd4b00f 1004
d2638bb7 1005 // TNamed* sigma = new TNamed("sigma",
1006 // (fIncludeSigma ? "included" : "excluded"));
241cca4d 1007 TObject* maxP = AliForwardUtil::MakeParameter("maxParticle", fMaxParticles);
1008 TObject* method = AliForwardUtil::MakeParameter("method", fUsePoisson);
1009 TObject* phiA = AliForwardUtil::MakeParameter("phiAcceptance",
1010 fUsePhiAcceptance);
1011 TObject* etaL = AliForwardUtil::MakeParameter("etaLumping", fEtaLumping);
1012 TObject* phiL = AliForwardUtil::MakeParameter("phiLumping", fPhiLumping);
5ca83fee 1013 TObject* reEt = AliForwardUtil::MakeParameter("recalcEta", fRecalculateEta);
1014 TObject* rePh = AliForwardUtil::MakeParameter("recalcPhi", fRecalculatePhi);
1015
bfab35d9 1016 TParameter<int>* nFiles = new TParameter<int>("nFiles", 1);
1017 nFiles->SetMergeMode('+');
1018
d2638bb7 1019 // d->Add(sigma);
f7cfc454 1020 d->Add(maxP);
1021 d->Add(method);
1022 d->Add(phiA);
1023 d->Add(etaL);
1024 d->Add(phiL);
5ca83fee 1025 d->Add(reEt);
1026 d->Add(rePh);
bfab35d9 1027 d->Add(nFiles);
d2638bb7 1028 // d->Add(nxi);
5ca83fee 1029 fCuts.Output(d,"lCuts");
f7cfc454 1030
7e4038b5 1031 TIter next(&fRingHistos);
1032 RingHistos* o = 0;
1033 while ((o = static_cast<RingHistos*>(next()))) {
e18cb8bd 1034 o->fPoisson.SetLumping(fEtaLumping, fPhiLumping);
5934a3e3 1035 o->CreateOutputObjects(d);
7e4038b5 1036 }
1f7aa5c7 1037
1038 if (!fDoTiming) return;
1039
1040 fHTiming = new TProfile("timing", "#LTt_{part}#GT", 8, .5, 8.5);
1041 fHTiming->SetDirectory(0);
1042 fHTiming->SetYTitle("#LTt_{part}#GT");
1043 fHTiming->SetXTitle("Part");
1044 fHTiming->SetFillColor(kRed+1);
1045 fHTiming->SetFillStyle(3001);
1046 fHTiming->SetMarkerStyle(20);
1047 fHTiming->SetMarkerColor(kBlack);
1048 fHTiming->SetLineColor(kBlack);
1049 fHTiming->SetStats(0);
1050 TAxis* xaxis = fHTiming->GetXaxis();
1051 xaxis->SetBinLabel(1, "Re-calculation of #eta");
1052 xaxis->SetBinLabel(2, "N_{particle}");
1053 xaxis->SetBinLabel(3, "Correction");
1054 xaxis->SetBinLabel(4, "Re-calculation of #phi");
1055 xaxis->SetBinLabel(5, "Copy to cache");
1056 xaxis->SetBinLabel(6, "Poisson calculation");
1057 xaxis->SetBinLabel(7, "Diagnostics");
1058 xaxis->SetBinLabel(8, "Total");
1059 d->Add(fHTiming);
7e4038b5 1060}
c8b1a7db 1061#define PF(N,V,...) \
1062 AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
1063#define PFB(N,FLAG) \
1064 do { \
1065 AliForwardUtil::PrintName(N); \
1066 std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
1067 } while(false)
1068#define PFV(N,VALUE) \
1069 do { \
1070 AliForwardUtil::PrintName(N); \
1071 std::cout << (VALUE) << std::endl; } while(false)
1072
0bd4b00f 1073//____________________________________________________________________
1074void
1174780f 1075AliFMDDensityCalculator::Print(Option_t* option) const
0bd4b00f 1076{
7984e5f7 1077 //
1078 // Print information
1079 //
1080 // Parameters:
1081 // option Not used
1082 //
c8b1a7db 1083 AliForwardUtil::PrintTask(*this);
1084 gROOT->IncreaseDirLevel();
1085
1086 TString phiM("none");
5ca83fee 1087 switch (fUsePhiAcceptance) {
c8b1a7db 1088 case kPhiNoCorrect: phiM = "none"; break;
1089 case kPhiCorrectNch: phiM = "correct Nch"; break;
1090 case kPhiCorrectELoss: phiM = "correct energy loss"; break;
5ca83fee 1091 }
c8b1a7db 1092
1093 PFV("Max(particles)", fMaxParticles );
1094 PFV("Poisson method", fUsePoisson );
1095 PFV("Eta lumping", fEtaLumping );
1096 PFV("Phi lumping", fPhiLumping );
1097 PFV("Recalculate eta", fRecalculateEta );
1098 PFV("Recalculate phi", fRecalculatePhi );
1099 PFV("Use phi acceptance", phiM);
1100 PFV("Lower cut", "");
d2638bb7 1101 fCuts.Print();
c8b1a7db 1102
1174780f 1103 TString opt(option);
1104 opt.ToLower();
c8b1a7db 1105 if (!opt.Contains("nomax")) {
1106 PFV("Max weights", "");
1107
1108 char ind[64];
1109 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
1110 ind[gROOT->GetDirLevel()] = '\0';
1111 for (UShort_t d=1; d<=3; d++) {
1112 UShort_t nr = (d == 1 ? 1 : 2);
1113 for (UShort_t q=0; q<nr; q++) {
1114 ind[gROOT->GetDirLevel()] = ' ';
1115 ind[gROOT->GetDirLevel()+1] = '\0';
1116 Char_t r = (q == 0 ? 'I' : 'O');
1117 std::cout << ind << " FMD" << d << r << ":";
1118 ind[gROOT->GetDirLevel()+1] = ' ';
1119 ind[gROOT->GetDirLevel()+2] = '\0';
1120
1121 const TArrayI& a = (d == 1 ? fFMD1iMax :
1122 (d == 2 ? (r == 'I' ? fFMD2iMax : fFMD2oMax) :
1123 (r == 'I' ? fFMD3iMax : fFMD3oMax)));
1124 Int_t j = 0;
1125 for (Int_t i = 0; i < a.fN; i++) {
1126 if (a.fArray[i] < 1) continue;
1127 if (j % 6 == 0) std::cout << "\n " << ind;
1128 j++;
1129 std::cout << " " << std::setw(3) << i << ": " << a.fArray[i];
1130 }
1131 std::cout << std::endl;
1174780f 1132 }
1174780f 1133 }
1134 }
c8b1a7db 1135 gROOT->DecreaseDirLevel();
0bd4b00f 1136}
7e4038b5 1137
1138//====================================================================
1139AliFMDDensityCalculator::RingHistos::RingHistos()
9d99b0dd 1140 : AliForwardUtil::RingHistos(),
5ca83fee 1141 fList(0),
8449e3e0 1142 // fEvsN(0),
1143 // fEvsM(0),
1144 // fEtaVsN(0),
1145 // fEtaVsM(0),
0bd4b00f 1146 fCorr(0),
9d05ffeb 1147 fDensity(0),
1148 fELossVsPoisson(0),
5ca83fee 1149 fDiffELossPoisson(0),
821ffd28 1150 fPoisson(),
f7cfc454 1151 fELoss(0),
1152 fELossUsed(0),
5ca83fee 1153 fMultCut(0),
1154 fTotal(0),
1155 fGood(0),
1156 fPhiAcc(0),
1157 fPhiBefore(0),
1158 fPhiAfter(0)
7984e5f7 1159{
1160 //
1161 // Default CTOR
1162 //
1163}
7e4038b5 1164
1165//____________________________________________________________________
1166AliFMDDensityCalculator::RingHistos::RingHistos(UShort_t d, Char_t r)
9d99b0dd 1167 : AliForwardUtil::RingHistos(d,r),
5ca83fee 1168 fList(0),
8449e3e0 1169 // fEvsN(0),
1170 // fEvsM(0),
1171 // fEtaVsN(0),
1172 // fEtaVsM(0),
0bd4b00f 1173 fCorr(0),
9d05ffeb 1174 fDensity(0),
1175 fELossVsPoisson(0),
5ca83fee 1176 fDiffELossPoisson(0),
821ffd28 1177 fPoisson("ignored"),
d2638bb7 1178 fELoss(0),
f7cfc454 1179 fELossUsed(0),
5ca83fee 1180 fMultCut(0),
1181 fTotal(0),
1182 fGood(0),
1183 fPhiAcc(0),
1184 fPhiBefore(0),
1185 fPhiAfter(0)
7e4038b5 1186{
7984e5f7 1187 //
1188 // Constructor
1189 //
1190 // Parameters:
1191 // d detector
1192 // r ring
1193 //
8449e3e0 1194#if 0
5bb5d1f6 1195 fEvsN = new TH2D("elossVsNnocorr",
1196 "#Delta E/#Delta E_{mip} vs uncorrected inclusive N_{ch}",
5ca83fee 1197 250, -.5, 24.5, 251, -1.5, 24.5);
7e4038b5 1198 fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
1199 fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
1200 fEvsN->Sumw2();
1201 fEvsN->SetDirectory(0);
5bb5d1f6 1202
1203 fEvsM = static_cast<TH2D*>(fEvsN->Clone("elossVsNcorr"));
1204 fEvsM->SetTitle("#Delta E/#Delta E_{mip} vs corrected inclusive N_{ch}");
7e4038b5 1205 fEvsM->SetDirectory(0);
1206
5bb5d1f6 1207 fEtaVsN = new TProfile("etaVsNnocorr",
1208 "Average inclusive N_{ch} vs #eta (uncorrected)",
1209 200, -4, 6);
0bd4b00f 1210 fEtaVsN->SetXTitle("#eta");
1211 fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
1212 fEtaVsN->SetDirectory(0);
1213 fEtaVsN->SetLineColor(Color());
1214 fEtaVsN->SetFillColor(Color());
5bb5d1f6 1215
1216 fEtaVsM = static_cast<TProfile*>(fEtaVsN->Clone("etaVsNcorr"));
1217 fEtaVsM->SetTitle("Average inclusive N_{ch} vs #eta (corrected)");
0bd4b00f 1218 fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
1219 fEtaVsM->SetDirectory(0);
8449e3e0 1220#endif
0bd4b00f 1221
5bb5d1f6 1222 fCorr = new TProfile("corr", "Average correction", 200, -4, 6);
0bd4b00f 1223 fCorr->SetXTitle("#eta");
1224 fCorr->SetYTitle("#LT correction#GT");
1225 fCorr->SetDirectory(0);
1226 fCorr->SetLineColor(Color());
1227 fCorr->SetFillColor(Color());
1228
5bb5d1f6 1229 fDensity = new TH2D("inclDensity", "Inclusive N_{ch} density",
7e4038b5 1230 200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
1231 0, 2*TMath::Pi());
1232 fDensity->SetDirectory(0);
5ca83fee 1233 fDensity->Sumw2(); fDensity->SetMarkerColor(Color());
7e4038b5 1234 fDensity->SetXTitle("#eta");
1235 fDensity->SetYTitle("#phi [radians]");
1236 fDensity->SetZTitle("Inclusive N_{ch} density");
9d05ffeb 1237
5bb5d1f6 1238 fELossVsPoisson = new TH2D("elossVsPoisson",
1239 "N_{ch} from energy loss vs from Poission",
3e4a0875 1240 500, 0, 100, 500, 0, 100);
9d05ffeb 1241 fELossVsPoisson->SetDirectory(0);
1242 fELossVsPoisson->SetXTitle("N_{ch} from #DeltaE");
1243 fELossVsPoisson->SetYTitle("N_{ch} from Poisson");
1244 fELossVsPoisson->SetZTitle("Correlation");
1245
5ca83fee 1246 fDiffELossPoisson = new TH1D("diffElossPoisson",
1247 "(N_{ch,#DeltaE}-N_{ch,Poisson})/N_{ch,Poisson}",
1248 100, -1, 1);
1249 fDiffELossPoisson->SetDirectory(0);
1250 fDiffELossPoisson->SetXTitle(fDiffELossPoisson->GetTitle());
1251 fDiffELossPoisson->SetYTitle("Frequency");
1252 fDiffELossPoisson->SetMarkerColor(Color());
1253 fDiffELossPoisson->SetFillColor(Color());
bfab35d9 1254 fDiffELossPoisson->SetFillStyle(3001);
5ca83fee 1255 fDiffELossPoisson->Sumw2();
1256
f7cfc454 1257 fELoss = new TH1D("eloss", "#Delta/#Delta_{mip} in all strips",
8449e3e0 1258 640, -1, 15);
f7cfc454 1259 fELoss->SetXTitle("#Delta/#Delta_{mip} (selected)");
1260 fELoss->SetYTitle("P(#Delta/#Delta_{mip})");
1261 fELoss->SetFillColor(Color()-2);
1262 fELoss->SetFillStyle(3003);
1263 fELoss->SetLineColor(kBlack);
1264 fELoss->SetLineStyle(2);
1265 fELoss->SetLineWidth(2);
1266 fELoss->SetDirectory(0);
1267
1268 fELossUsed = static_cast<TH1D*>(fELoss->Clone("elossUsed"));
1269 fELossUsed->SetTitle("#Delta/#Delta_{mip} in used strips");
1270 fELossUsed->SetFillStyle(3002);
1271 fELossUsed->SetLineStyle(1);
1272 fELossUsed->SetDirectory(0);
5ca83fee 1273
1274 fPhiBefore = new TH1D("phiBefore", "#phi distribution (before recalc)",
1275 (r == 'I' || r == 'i' ? 20 : 40), 0, 2*TMath::Pi());
1276 fPhiBefore->SetDirectory(0);
1277 fPhiBefore->SetXTitle("#phi");
1278 fPhiBefore->SetYTitle("Events");
1279 fPhiBefore->SetMarkerColor(Color());
1280 fPhiBefore->SetLineColor(Color());
1281 fPhiBefore->SetFillColor(Color());
bfab35d9 1282 fPhiBefore->SetFillStyle(3001);
5ca83fee 1283 fPhiBefore->SetMarkerStyle(20);
1284
1285 fPhiAfter = static_cast<TH1D*>(fPhiBefore->Clone("phiAfter"));
1286 fPhiAfter->SetTitle("#phi distribution (after re-calc)");
1287 fPhiAfter->SetDirectory(0);
f7cfc454 1288
7e4038b5 1289}
1290//____________________________________________________________________
1291AliFMDDensityCalculator::RingHistos::RingHistos(const RingHistos& o)
5ca83fee 1292 : AliForwardUtil::RingHistos(o),
1293 fList(o.fList),
8449e3e0 1294 // fEvsN(o.fEvsN),
1295 // fEvsM(o.fEvsM),
1296 // fEtaVsN(o.fEtaVsN),
1297 // fEtaVsM(o.fEtaVsM),
0bd4b00f 1298 fCorr(o.fCorr),
9d05ffeb 1299 fDensity(o.fDensity),
1300 fELossVsPoisson(o.fELossVsPoisson),
5ca83fee 1301 fDiffELossPoisson(o.fDiffELossPoisson),
821ffd28 1302 fPoisson(o.fPoisson),
f7cfc454 1303 fELoss(o.fELoss),
1304 fELossUsed(o.fELossUsed),
5ca83fee 1305 fMultCut(o.fMultCut),
1306 fTotal(o.fTotal),
1307 fGood(o.fGood),
1308 fPhiAcc(o.fPhiAcc),
1309 fPhiBefore(o.fPhiBefore),
1310 fPhiAfter(o.fPhiAfter)
7984e5f7 1311{
1312 //
1313 // Copy constructor
1314 //
1315 // Parameters:
1316 // o Object to copy from
1317 //
1318}
7e4038b5 1319
1320//____________________________________________________________________
1321AliFMDDensityCalculator::RingHistos&
1322AliFMDDensityCalculator::RingHistos::operator=(const RingHistos& o)
1323{
7984e5f7 1324 //
1325 // Assignment operator
1326 //
1327 // Parameters:
1328 // o Object to assign from
1329 //
1330 // Return:
1331 // Reference to this
1332 //
d015ecfe 1333 if (&o == this) return *this;
9d99b0dd 1334 AliForwardUtil::RingHistos::operator=(o);
7e4038b5 1335
8449e3e0 1336 // if (fEvsN) delete fEvsN;
1337 // if (fEvsM) delete fEvsM;
1338 // if (fEtaVsN) delete fEtaVsN;
1339 // if (fEtaVsM) delete fEtaVsM;
5ca83fee 1340 if (fCorr) delete fCorr;
1341 if (fDensity) delete fDensity;
1342 if (fELossVsPoisson) delete fELossVsPoisson;
1343 if (fDiffELossPoisson) delete fDiffELossPoisson;
1344 if (fTotal) delete fTotal;
1345 if (fGood) delete fGood;
1346 if (fPhiAcc) delete fPhiAcc;
1347 if (fPhiBefore) delete fPhiBefore;
1348 if (fPhiAfter) delete fPhiAfter;
1349
8449e3e0 1350 // fEvsN = static_cast<TH2D*>(o.fEvsN->Clone());
1351 // fEvsM = static_cast<TH2D*>(o.fEvsM->Clone());
1352 // fEtaVsN = static_cast<TProfile*>(o.fEtaVsN->Clone());
1353 // fEtaVsM = static_cast<TProfile*>(o.fEtaVsM->Clone());
5ca83fee 1354 fCorr = static_cast<TProfile*>(o.fCorr->Clone());
1355 fDensity = static_cast<TH2D*>(o.fDensity->Clone());
1356 fELossVsPoisson = static_cast<TH2D*>(o.fELossVsPoisson->Clone());
1357 fDiffELossPoisson = static_cast<TH1D*>(o.fDiffELossPoisson->Clone());
1358 fPoisson = o.fPoisson;
1359 fELoss = static_cast<TH1D*>(o.fELoss->Clone());
1360 fELossUsed = static_cast<TH1D*>(o.fELossUsed->Clone());
1361 fTotal = static_cast<TH1D*>(o.fTotal->Clone());
1362 fGood = static_cast<TH1D*>(o.fGood->Clone());
1363 fPhiAcc = static_cast<TH2D*>(o.fPhiAcc->Clone());
1364 fPhiBefore = static_cast<TH1D*>(o.fPhiBefore->Clone());
1365 fPhiAfter = static_cast<TH1D*>(o.fPhiAfter->Clone());
7e4038b5 1366 return *this;
1367}
1368//____________________________________________________________________
1369AliFMDDensityCalculator::RingHistos::~RingHistos()
1370{
7984e5f7 1371 //
1372 // Destructor
1373 //
9d05ffeb 1374}
1375
e308a636 1376
1377//____________________________________________________________________
9d05ffeb 1378void
5934a3e3 1379AliFMDDensityCalculator::RingHistos::SetupForData(const TAxis& eAxis)
9d05ffeb 1380{
5ca83fee 1381 // Initialize
1382 // This is called on first event
e18cb8bd 1383 fPoisson.Init(-1,-1);
5ca83fee 1384 fTotal = new TH1D("total", "Total # of strips per #eta",
1385 eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax());
1386 fTotal->SetDirectory(0);
1387 fTotal->SetXTitle("#eta");
1388 fTotal->SetYTitle("# of strips");
1389 fGood = static_cast<TH1D*>(fTotal->Clone("good"));
1390 fGood->SetTitle("# of good strips per #eta");
1391 fGood->SetDirectory(0);
1392
1393 fPhiAcc = new TH2D("phiAcc", "#phi acceptance vs Ip_{z}",
1394 eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
1395 10, -10, 10);
1396 fPhiAcc->SetXTitle("#eta");
1397 fPhiAcc->SetYTitle("v_{z} [cm]");
1398 fPhiAcc->SetZTitle("#phi acceptance");
1399 fPhiAcc->SetDirectory(0);
1400
1401 if (fList) fList->Add(fPhiAcc);
7e4038b5 1402}
1403
1404//____________________________________________________________________
1405void
5934a3e3 1406AliFMDDensityCalculator::RingHistos::CreateOutputObjects(TList* dir)
7e4038b5 1407{
7984e5f7 1408 //
5ca83fee 1409 // Make output. This is called as part of SlaveBegin
7984e5f7 1410 //
1411 // Parameters:
1412 // dir Where to put it
1413 //
9d99b0dd 1414 TList* d = DefineOutputList(dir);
8449e3e0 1415 // d->Add(fEvsN);
1416 // d->Add(fEvsM);
1417 // d->Add(fEtaVsN);
1418 // d->Add(fEtaVsM);
0bd4b00f 1419 d->Add(fCorr);
7e4038b5 1420 d->Add(fDensity);
9d05ffeb 1421 d->Add(fELossVsPoisson);
5ca83fee 1422 d->Add(fDiffELossPoisson);
821ffd28 1423 fPoisson.Output(d);
d23503ee 1424 fPoisson.GetOccupancy()->SetFillColor(Color());
1425 fPoisson.GetMean()->SetFillColor(Color());
1426 fPoisson.GetOccupancy()->SetFillColor(Color());
f7cfc454 1427 d->Add(fELoss);
1428 d->Add(fELossUsed);
5ca83fee 1429 d->Add(fPhiBefore);
1430 d->Add(fPhiAfter);
d23503ee 1431
1432 Bool_t inner = (fRing == 'I' || fRing == 'i');
1433 Int_t nStr = inner ? 512 : 256;
1434 Int_t nSec = inner ? 20 : 40;
1435 TAxis x(nStr, -.5, nStr-.5);
1436 TAxis y(nSec, -.5, nSec-.5);
1437 x.SetTitle("strip");
1438 y.SetTitle("sector");
1439 fPoisson.Define(x, y);
1440
241cca4d 1441 d->Add(AliForwardUtil::MakeParameter("cut", fMultCut));
5ca83fee 1442 fList = d;
9d99b0dd 1443}
1444
1445//____________________________________________________________________
1446void
5934a3e3 1447AliFMDDensityCalculator::RingHistos::Terminate(TList* dir, Int_t nEvents)
9d99b0dd 1448{
7984e5f7 1449 //
1450 // Scale the histograms to the total number of events
1451 //
1452 // Parameters:
1453 // dir Where the output is
1454 // nEvents Number of events
1455 //
9d99b0dd 1456 TList* l = GetOutputList(dir);
1457 if (!l) return;
1458
5934a3e3 1459 TH2D* density = static_cast<TH2D*>(GetOutputHist(l,"inclDensity"));
9d99b0dd 1460 if (density) density->Scale(1./nEvents);
5934a3e3 1461 fDensity = density;
7e4038b5 1462}
1463
1464//____________________________________________________________________
1465//
1466// EOF
1467//
1468
1469
1470