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