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