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