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