AliPhysics  fde8a9f (fde8a9f)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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),
80  fUsePhiAcceptance(kPhiCorrectNch),
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  100, 0, 10);
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(),
144  fSumOfWeights(o.fSumOfWeights),
145  fWeightedSum(o.fWeightedSum),
146  fCorrections(o.fCorrections),
147  fMaxParticles(o.fMaxParticles),
148  fUsePoisson(o.fUsePoisson),
149  fUsePhiAcceptance(o.fUsePhiAcceptance),
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),
157  fMaxWeights(o.fMaxWeights),
158  fLowCuts(o.fLowCuts),
159  fEtaLumping(o.fEtaLumping),
160  fPhiLumping(o.fPhiLumping),
161  fDebug(o.fDebug),
162  fCuts(o.fCuts),
163  fRecalculatePhi(o.fRecalculatePhi),
164  fMinQuality(o.fMinQuality),
165  fHitThreshold(o.fHitThreshold),
166  fCache(o.fCache),
167  fDoTiming(o.fDoTiming),
168  fHTiming(o.fHTiming),
169  fMaxOutliers(o.fMaxOutliers),
170  fOutlierCut(o.fOutlierCut)
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  }
530  rh->fPoisson.Fill(t,s,hit,1./c);
531  h->Fill(eta,phi,n);
532 
533  // --- If we use ELoss fits, apply now ---------------------
534  if (!fUsePoisson) rh->fDensity->Fill(eta,phi,n);
535  } // for t
536  } // for s
537 
538  // --- Automatic acceptance - Calculate as an efficiency -------
539  // This is very fast, so we do not bother to time it
540  rh->fGood->Divide(rh->fGood, rh->fTotal, 1, 1, "B");
541 
542  // --- Make a copy and reset as needed -------------------------
543  START_TIMER(timer);
544  TH2D* hclone = fCache.Get(d,r);
545  // hclone->Reset();
546  // TH2D* hclone = static_cast<TH2D*>(h->Clone("hclone"));
547  if (!fUsePoisson) hclone->Reset();
548  else {
549  for (Int_t i = 0; i <= h->GetNbinsX()+1; i++) {
550  for (Int_t j = 0; j <= h->GetNbinsY()+1; j++) {
551  hclone->SetBinContent(i,j,h->GetBinContent(i,j));
552  hclone->SetBinError(i,j,h->GetBinError(i,j));
553  }
554  }
555  // hclone->Add(h);
556  h->Reset();
557  }
558  ADD_TIMER(timer,copyTime);
559 
560  // --- Store Poisson result ------------------------------------
561  START_TIMER(timer);
562  TH2D* poisson = rh->fPoisson.Result();
563  for (Int_t t=0; t < poisson->GetNbinsX(); t++) {
564  for (Int_t s=0; s < poisson->GetNbinsY(); s++) {
565 
566  Double_t poissonV = poisson->GetBinContent(t+1,s+1);
567  // Use cached eta - since the calls to GetEtaFromStrip and
568  // GetPhiFromStrip are _very_ expensive
569  Double_t phi = phiCache[s*nt+t];
570  Double_t eta = etaCache[s*nt+t];
571  // Double_t phi = fmd.Phi(d,r,s,t) * TMath::DegToRad();
572  // Double_t eta = fmd.Eta(d,r,s,t);
573  if (fUsePoisson) {
574  h->Fill(eta,phi,poissonV);
575  rh->fDensity->Fill(eta, phi, poissonV);
576  }
577  else
578  hclone->Fill(eta,phi,poissonV);
579  }
580  }
581  ADD_TIMER(timer,poissonTime);
582 
583  // --- Make diagnostics - eloss vs poisson ---------------------
584  START_TIMER(timer);
585  Int_t nY = h->GetNbinsY();
586  Int_t nIn = 0; // Count non-outliers
587  Int_t nOut = 0; // Count outliers
588  for (Int_t ieta=1; ieta <= h->GetNbinsX(); ieta++) {
589  // Set the overflow bin to contain the phi acceptance
590  Double_t phiAcc = rh->fGood->GetBinContent(ieta);
591  Double_t phiAccE = rh->fGood->GetBinError(ieta);
592  h->SetBinContent(ieta, nY+1, phiAcc);
593  h->SetBinError(ieta, nY+1, phiAccE);
594  Double_t eta = h->GetXaxis()->GetBinCenter(ieta);
595  rh->fPhiAcc->Fill(eta, ip.Z(), phiAcc);
596  for (Int_t iphi=1; iphi<= nY; iphi++) {
597 
598  Double_t poissonV = 0; //h->GetBinContent(,s+1);
599  Double_t eLossV = 0;
600  if(fUsePoisson) {
601  poissonV = h->GetBinContent(ieta,iphi);
602  eLossV = hclone->GetBinContent(ieta,iphi);
603  }
604  else {
605  poissonV = hclone->GetBinContent(ieta,iphi);
606  eLossV = h->GetBinContent(ieta,iphi);
607  }
608 
609  if (poissonV < 1e-12 && eLossV < 1e-12)
610  // we do not care about trivially empty bins
611  continue;
612 
613  Bool_t outlier = CheckOutlier(eLossV, poissonV, fOutlierCut);
614  Double_t rel = eLossV < 1e-12 ? 0 : (poissonV - eLossV) / eLossV;
615  if (outlier) {
616  rh->fELossVsPoissonOut->Fill(eLossV, poissonV);
617  rh->fDiffELossPoissonOut->Fill(rel);
618  nOut++;
619  }
620  else {
621  rh->fELossVsPoisson->Fill(eLossV, poissonV);
622  rh->fDiffELossPoisson->Fill(rel);
623  nIn++;
624  } // if (outlier)
625  } // for (iphi)
626  } // for (ieta)
627  Int_t nTotal = (nIn+nOut);
628  Double_t outRatio = (nTotal > 0 ? Double_t(nOut) / nTotal : 0);
629  rh->fOutliers->Fill(outRatio);
630  if (outRatio < fMaxOutliers) rh->fPoisson.FillDiagnostics();
631  else h->SetBit(AliForwardUtil::kSkipRing);
632  ADD_TIMER(timer,diagTime);
633  // delete hclone;
634 
635  } // for q
636  } // for d
637 
638  if (fDoTiming) {
639  // fHTiming->Fill(1,reEtaTime);
640  fHTiming->Fill(2,nPartTime);
641  fHTiming->Fill(3,corrTime);
642  fHTiming->Fill(4,rePhiTime);
643  fHTiming->Fill(5,copyTime);
644  fHTiming->Fill(6,poissonTime);
645  fHTiming->Fill(7,diagTime);
646  fHTiming->Fill(8,totalT.CpuTime());
647  }
648 
649  return kTRUE;
650 }
651 
652 //_____________________________________________________________________
653 Bool_t
655  Double_t poisson,
656  Double_t cut) const
657 {
658  if (eloss < 1e-6) return true;
659  Double_t diff = TMath::Abs(poisson - eloss) / eloss;
660  return diff > cut;
661 }
662 //_____________________________________________________________________
663 Int_t
665  UShort_t d, Char_t r, Int_t iEta) const
666 {
667  //
668  // Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
669  //
670  // Parameters:
671  // cor Correction
672  // d Detector
673  // r Ring
674  // iEta Eta bin
675  //
676  DGUARD(fDebug, 10, "Find maximum weight in FMD density calculator");
677  if(!cor) return -1;
678 
679  AliFMDCorrELossFit::ELossFit* fit = cor->FindFit(d,r,iEta, -1);
680  if (!fit) {
681  // AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
682  return -1;
683  }
686  fMaxParticles);
687 }
688 //_____________________________________________________________________
689 Int_t
691  UShort_t d, Char_t r, Double_t eta) const
692 {
693  //
694  // Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
695  //
696  // Parameters:
697  // cor Correction
698  // d Detector
699  // r Ring
700  // eta Eta
701  //
702  DGUARD(fDebug, 10, "Find maximum weight in FMD density calculator");
703  if(!cor) return -1;
704 
705  AliFMDCorrELossFit::ELossFit* fit = cor->FindFit(d,r,eta, -1);
706  if (!fit) {
707  // AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f", d, r, eta));
708  return -1;
709  }
712  fMaxParticles);
713 }
714 
715 //_____________________________________________________________________
716 void
718 {
719  //
720  // Find the max weights and cache them
721  //
722  DGUARD(fDebug, 2, "Cache maximum weights in FMD density calculator");
724  const AliFMDCorrELossFit* cor = fcm.GetELossFit();
725  cor->CacheBins(fMinQuality);
726  cor->Print(fDebug > 5 ? "RCS" : "C");
727 
728  TAxis eta(axis.GetNbins(),
729  axis.GetXmin(),
730  axis.GetXmax());
731 
732  eta.Set(cor->GetEtaAxis().GetNbins(),
733  cor->GetEtaAxis().GetXmin(),
734  cor->GetEtaAxis().GetXmax());
735 
736  Int_t nEta = eta.GetNbins();
737  fFMD1iMax.Set(nEta);
738  fFMD2iMax.Set(nEta);
739  fFMD2oMax.Set(nEta);
740  fFMD3iMax.Set(nEta);
741  fFMD3oMax.Set(nEta);
742 
743  fMaxWeights->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
744  fMaxWeights->GetYaxis()->SetBinLabel(1, "FMD1i");
745  fMaxWeights->GetYaxis()->SetBinLabel(2, "FMD2i");
746  fMaxWeights->GetYaxis()->SetBinLabel(3, "FMD2o");
747  fMaxWeights->GetYaxis()->SetBinLabel(4, "FMD3i");
748  fMaxWeights->GetYaxis()->SetBinLabel(5, "FMD3o");
749 
750  AliInfo(Form("Get eta axis with %d bins from %f to %f",
751  nEta, eta.GetXmin(), eta.GetXmax()));
752  fLowCuts->SetBins(nEta, eta.GetXmin(), eta.GetXmax(), 5, .5, 5.5);
753  fLowCuts->GetYaxis()->SetBinLabel(1, "FMD1i");
754  fLowCuts->GetYaxis()->SetBinLabel(2, "FMD2i");
755  fLowCuts->GetYaxis()->SetBinLabel(3, "FMD2o");
756  fLowCuts->GetYaxis()->SetBinLabel(4, "FMD3i");
757  fLowCuts->GetYaxis()->SetBinLabel(5, "FMD3o");
758 
759  for (Int_t i = 0; i < nEta; i++) {
760  Double_t leta = fMaxWeights->GetXaxis()->GetBinCenter(i+1);
761  Double_t w[5];
762  w[0] = fFMD1iMax[i] = FindMaxWeight(cor, 1, 'I', leta);
763  w[1] = fFMD2iMax[i] = FindMaxWeight(cor, 2, 'I', leta);
764  w[2] = fFMD2oMax[i] = FindMaxWeight(cor, 2, 'O', leta);
765  w[3] = fFMD3iMax[i] = FindMaxWeight(cor, 3, 'I', leta);
766  w[4] = fFMD3oMax[i] = FindMaxWeight(cor, 3, 'O', leta);
767  for (Int_t j = 0; j < 5; j++)
768  if (w[j] > 0) fMaxWeights->SetBinContent(i+1, j+1, w[j]);
769  }
770 
771  // Cache cuts in histogram
773 }
774 
775 //_____________________________________________________________________
776 Int_t
778 {
779  //
780  // Find the (cached) maximum weight for FMD<i>dr</i> in
781  // @f$\eta@f$ bin @a iEta
782  //
783  // Parameters:
784  // d Detector
785  // r Ring
786  // iEta Eta bin
787  //
788  // Return:
789  // max weight or <= 0 in case of problems
790  //
791  if (iEta < 0) return -1;
792 
793  const TArrayI* max = 0;
794  switch (d) {
795  case 1: max = &fFMD1iMax; break;
796  case 2: max = (r == 'I' || r == 'i' ? &fFMD2iMax : &fFMD2oMax); break;
797  case 3: max = (r == 'I' || r == 'i' ? &fFMD3iMax : &fFMD3oMax); break;
798  }
799  if (!max) {
800  AliWarning(Form("No array for FMD%d%c", d, r));
801  return -1;
802  }
803 
804  if (iEta >= max->fN) {
805  AliWarning(Form("Eta bin %3d out of bounds [0,%d]",
806  iEta, max->fN-1));
807  return -1;
808  }
809 
810  AliDebug(30,Form("Max weight for FMD%d%c eta bin %3d: %d", d, r, iEta,
811  max->At(iEta)));
812  return max->At(iEta);
813 }
814 
815 //_____________________________________________________________________
816 Int_t
818 {
819  //
820  // Find the (cached) maximum weight for FMD<i>dr</i> iat
821  // @f$\eta@f$
822  //
823  // Parameters:
824  // d Detector
825  // r Ring
826  // eta Eta bin
827  //
828  // Return:
829  // max weight or <= 0 in case of problems
830  //
832  Int_t iEta = fcm.GetELossFit()->FindEtaBin(eta) -1;
833 
834  return GetMaxWeight(d, r, iEta);
835 }
836 
837 //_____________________________________________________________________
838 Float_t
840  UShort_t d,
841  Char_t r,
842  Float_t eta,
843  Bool_t lowFlux) const
844 {
845  //
846  // Get the number of particles corresponding to the signal mult
847  //
848  // Parameters:
849  // mult Signal
850  // d Detector
851  // r Ring
852  // s Sector
853  // t Strip (not used)
854  // v Vertex bin
855  // eta Pseudo-rapidity
856  // lowFlux Low-flux flag
857  //
858  // Return:
859  // The number of particles
860  //
861  // if (mult <= GetMultCut()) return 0;
862  DGUARD(fDebug, 3, "Calculate Nch in FMD density calculator");
863  if (lowFlux) return 1;
864 
866  AliFMDCorrELossFit::ELossFit* fit = fcm.GetELossFit()->FindFit(d,r,eta, -1);
867  if (!fit) {
868  AliWarning(Form("No energy loss fit for FMD%d%c at eta=%f qual=%d",
869  d, r, eta, fMinQuality));
870  return 0;
871  }
872 
873  Int_t m = GetMaxWeight(d,r,eta); // fit->FindMaxWeight();
874  if (m < 1) {
875  AliWarning(Form("No good fits for FMD%d%c at eta=%f", d, r, eta));
876  return 0;
877  }
878 
879  UShort_t n = TMath::Min(fMaxParticles, UShort_t(m));
880  Double_t ret = fit->EvaluateWeighted(mult, n);
881 
882  if (fDebug > 10) {
883  AliInfo(Form("FMD%d%c, eta=%7.4f, %8.5f -> %8.5f", d, r, eta, mult, ret));
884  }
885 
886  fWeightedSum->Fill(ret);
887  fSumOfWeights->Fill(ret);
888 
889  return ret;
890 }
891 
892 //_____________________________________________________________________
893 Float_t
895  Char_t r,
896  UShort_t t,
897  Float_t eta,
898  Bool_t lowFlux) const
899 {
900  //
901  // Get the inverse correction factor. This consist of
902  //
903  // - acceptance correction (corners of sensors)
904  // - double hit correction (for low-flux events)
905  // - dead strip correction
906  //
907  // Parameters:
908  // d Detector
909  // r Ring
910  // s Sector
911  // t Strip (not used)
912  // v Vertex bin
913  // eta Pseudo-rapidity
914  // lowFlux Low-flux flag
915  //
916  // Return:
917  //
918  //
919  DGUARD(fDebug, 10, "Apply correction in FMD density calculator");
920  Float_t correction = 1;
922  correction = AcceptanceCorrection(r,t);
923  if (lowFlux) {
925 
926  TH1D* dblHitCor = 0;
927  if (fcm.GetDoubleHit())
928  dblHitCor = fcm.GetDoubleHit()->GetCorrection(d,r);
929 
930  if (dblHitCor) {
931  Double_t dblC = dblHitCor->GetBinContent(dblHitCor->FindBin(eta));
932  if (dblC > 0) correction *= dblC;
933  }
934  // else {
935  // AliWarning(Form("Missing double hit correction for FMD%d%c",d,r));
936  // }
937  }
938  return correction;
939 }
940 
941 //_____________________________________________________________________
942 TH1D*
944 {
945  //
946  // Generate the acceptance corrections
947  //
948  // Parameters:
949  // r Ring to generate for
950  //
951  // Return:
952  // Newly allocated histogram of acceptance corrections
953  //
954  DGUARD(fDebug, 3, "Make acceptance correction in FMD density calculator");
955  const Double_t ic1[] = { 4.9895, 15.3560 };
956  const Double_t ic2[] = { 1.8007, 17.2000 };
957  const Double_t oc1[] = { 4.2231, 26.6638 };
958  const Double_t oc2[] = { 1.8357, 27.9500 };
959  const Double_t* c1 = (r == 'I' || r == 'i' ? ic1 : oc1);
960  const Double_t* c2 = (r == 'I' || r == 'i' ? ic2 : oc2);
961  Double_t minR = (r == 'I' || r == 'i' ? 4.5213 : 15.4);
962  Double_t maxR = (r == 'I' || r == 'i' ? 17.2 : 28.0);
963  Int_t nStrips = (r == 'I' || r == 'i' ? 512 : 256);
964  Int_t nSec = (r == 'I' || r == 'i' ? 20 : 40);
965  Float_t basearc = 2 * TMath::Pi() / nSec;
966  Double_t rad = maxR - minR;
967  Float_t segment = rad / nStrips;
968  Float_t cr = TMath::Sqrt(c1[0]*c1[0]+c1[1]*c1[1]);
969 
970  // Numbers used to find end-point of strip.
971  // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
972  Float_t D = c1[0] * c2[1] - c1[1] * c2[0];
973  Float_t dx = c2[0] - c1[0];
974  Float_t dy = c2[1] - c1[1];
975  Float_t dr = TMath::Sqrt(dx*dx+dy*dy);
976 
977  TH1D* ret = new TH1D(Form("acc%c", r),
978  Form("Acceptance correction for FMDx%c", r),
979  nStrips, -.5, nStrips-.5);
980  ret->SetXTitle("Strip");
981  ret->SetYTitle("#varphi acceptance");
982  ret->SetDirectory(0);
983  ret->SetFillColor(r == 'I' || r == 'i' ? kRed+1 : kBlue+1);
984  ret->SetFillStyle(3001);
985 
986  for (Int_t t = 0; t < nStrips; t++) {
987  Float_t radius = minR + t * segment;
988 
989  // If the radius of the strip is smaller than the radius corresponding
990  // to the first corner we have a full strip length
991  if (radius <= cr) {
992  ret->SetBinContent(t+1, 1);
993  continue;
994  }
995 
996  // Next, we should find the end-point of the strip - that is,
997  // the coordinates where the line from c1 to c2 intersects a circle
998  // with radius given by the strip.
999  // (See http://mathworld.wolfram.com/Circle-LineIntersection.html)
1000  // Calculate the determinant
1001  Float_t det = radius * radius * dr * dr - D*D;
1002 
1003  if (det <= 0) {
1004  // <0 means No intersection
1005  // =0 means Exactly tangent
1006  ret->SetBinContent(t+1, 1);
1007  continue;
1008  }
1009 
1010  // Calculate end-point and the corresponding opening angle
1011  Float_t x = (+D * dy + dx * TMath::Sqrt(det)) / dr / dr;
1012  Float_t y = (-D * dx + dy * TMath::Sqrt(det)) / dr / dr;
1013  Float_t th = TMath::ATan2(x, y);
1014 
1015  ret->SetBinContent(t+1, th / basearc);
1016  }
1017  return ret;
1018 }
1019 
1020 //_____________________________________________________________________
1021 Float_t
1023 {
1024  //
1025  // Get the acceptance correction for strip @a t in an ring of type @a r
1026  //
1027  // Parameters:
1028  // r Ring type ('I' or 'O')
1029  // t Strip number
1030  //
1031  // Return:
1032  // Inverse acceptance correction
1033  //
1034  TH1D* acc = (r == 'I' || r == 'i' ? fAccI : fAccO);
1035  return acc->GetBinContent(t+1);
1036 }
1037 
1038 //____________________________________________________________________
1039 void
1041  Int_t nEvents)
1042 {
1043  //
1044  // Scale the histograms to the total number of events
1045  //
1046  // Parameters:
1047  // dir where to put the output
1048  // nEvents Number of events
1049  //
1050  DGUARD(fDebug, 1, "Scale histograms in FMD density calculator");
1051  if (nEvents <= 0) return;
1052  TList* d = static_cast<TList*>(dir->FindObject(GetName()));
1053  if (!d) return;
1054 
1055  TList* out = new TList;
1056  out->SetName(d->GetName());
1057  out->SetOwner();
1058 
1059  fRingHistos.Add(new RingHistos(1, 'I'));
1060  fRingHistos.Add(new RingHistos(2, 'I'));
1061  fRingHistos.Add(new RingHistos(2, 'O'));
1062  fRingHistos.Add(new RingHistos(3, 'I'));
1063  fRingHistos.Add(new RingHistos(3, 'O'));
1064  TIter next(&fRingHistos);
1065  RingHistos* o = 0;
1066  THStack* sums = new THStack("sums", "sums of ring signals");
1067  while ((o = static_cast<RingHistos*>(next()))) {
1068  o->Terminate(d, nEvents);
1069  if (!o->fDensity) {
1070  Warning("Terminate", "No density in %s", o->GetName());
1071  continue;
1072  }
1073  TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
1074  o->fDensity->GetNbinsY(),"e");
1075  sum->Scale(1., "width");
1076  sum->SetTitle(o->GetName());
1077  sum->SetDirectory(0);
1078  sum->SetYTitle("#sum N_{ch,incl}");
1079  sums->Add(sum);
1080  }
1081  out->Add(sums);
1082  output->Add(out);
1083 }
1084 
1085 //____________________________________________________________________
1086 void
1088 {
1089  //
1090  // Output diagnostic histograms to directory
1091  //
1092  // Parameters:
1093  // dir List to write in
1094  //
1095  DGUARD(fDebug, 1, "Define output FMD density calculator");
1096  TList* d = new TList;
1097  d->SetOwner();
1098  d->SetName(GetName());
1099  dir->Add(d);
1100  d->Add(fWeightedSum);
1101  d->Add(fSumOfWeights);
1102  d->Add(fCorrections);
1103  d->Add(fAccI);
1104  d->Add(fAccO);
1105  d->Add(fMaxWeights);
1106  d->Add(fLowCuts);
1107 
1108  TParameter<int>* nFiles = new TParameter<int>("nFiles", 1);
1109  nFiles->SetMergeMode('+');
1110 
1111  // d->Add(sigma);
1112  d->Add(AliForwardUtil::MakeParameter("maxParticle", fMaxParticles));
1113  d->Add(AliForwardUtil::MakeParameter("minQuality", fMinQuality));
1114  d->Add(AliForwardUtil::MakeParameter("method", fUsePoisson));
1115  d->Add(AliForwardUtil::MakeParameter("phiAcceptance",fUsePhiAcceptance));
1116  d->Add(AliForwardUtil::MakeParameter("etaLumping", fEtaLumping));
1117  d->Add(AliForwardUtil::MakeParameter("phiLumping", fPhiLumping));
1118  d->Add(AliForwardUtil::MakeParameter("recalcPhi", fRecalculatePhi));
1119  d->Add(AliForwardUtil::MakeParameter("maxOutliers", fMaxOutliers));
1120  d->Add(AliForwardUtil::MakeParameter("outlierCut", fOutlierCut));
1121  d->Add(AliForwardUtil::MakeParameter("hitThreshold", fHitThreshold));
1122  d->Add(nFiles);
1123  // d->Add(nxi);
1124  fCuts.Output(d,"lCuts");
1125 
1126  fRingHistos.Add(new RingHistos(1, 'I'));
1127  fRingHistos.Add(new RingHistos(2, 'I'));
1128  fRingHistos.Add(new RingHistos(2, 'O'));
1129  fRingHistos.Add(new RingHistos(3, 'I'));
1130  fRingHistos.Add(new RingHistos(3, 'O'));
1131  TIter next(&fRingHistos);
1132  RingHistos* o = 0;
1133  while ((o = static_cast<RingHistos*>(next()))) {
1135  o->CreateOutputObjects(d);
1136  }
1137 
1138  if (!fDoTiming) return;
1139 
1140  fHTiming = new TProfile("timing", "#LTt_{part}#GT", 8, .5, 8.5);
1141  fHTiming->SetDirectory(0);
1142  fHTiming->SetYTitle("#LTt_{part}#GT");
1143  fHTiming->SetXTitle("Part");
1144  fHTiming->SetFillColor(kRed+1);
1145  fHTiming->SetFillStyle(3001);
1146  fHTiming->SetMarkerStyle(20);
1147  fHTiming->SetMarkerColor(kBlack);
1148  fHTiming->SetLineColor(kBlack);
1149  fHTiming->SetStats(0);
1150  TAxis* xaxis = fHTiming->GetXaxis();
1151  xaxis->SetBinLabel(1, "Re-calculation of #eta");
1152  xaxis->SetBinLabel(2, "N_{particle}");
1153  xaxis->SetBinLabel(3, "Correction");
1154  xaxis->SetBinLabel(4, "Re-calculation of #phi");
1155  xaxis->SetBinLabel(5, "Copy to cache");
1156  xaxis->SetBinLabel(6, "Poisson calculation");
1157  xaxis->SetBinLabel(7, "Diagnostics");
1158  xaxis->SetBinLabel(8, "Total");
1159  d->Add(fHTiming);
1160 }
1161 #define PF(N,V,...) \
1162  AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
1163 #define PFB(N,FLAG) \
1164  do { \
1165  AliForwardUtil::PrintName(N); \
1166  std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
1167  } while(false)
1168 #define PFV(N,VALUE) \
1169  do { \
1170  AliForwardUtil::PrintName(N); \
1171  std::cout << (VALUE) << std::endl; } while(false)
1172 
1173 //____________________________________________________________________
1174 void
1176 {
1177  //
1178  // Print information
1179  //
1180  // Parameters:
1181  // option Not used
1182  //
1184  gROOT->IncreaseDirLevel();
1185 
1186  TString phiM("none");
1187  switch (fUsePhiAcceptance) {
1188  case kPhiNoCorrect: phiM = "none"; break;
1189  case kPhiCorrectNch: phiM = "correct Nch"; break;
1190  case kPhiCorrectELoss: phiM = "correct energy loss"; break;
1191  }
1192 
1193  PFV("Max(particles)", fMaxParticles);
1194  PFB("Poisson method", fUsePoisson);
1195  PFV("Eta lumping", fEtaLumping);
1196  PFV("Phi lumping", fPhiLumping);
1197  PFB("Recalculate phi", fRecalculatePhi);
1198  PFB("Use phi acceptance", phiM);
1199  PFV("Min(quality)", fMinQuality);
1200  PFV("Threshold(hit)", fHitThreshold);
1201  PFV("Max(outliers)", fMaxOutliers);
1202  PFV("Cut(outlier)", fOutlierCut);
1203  PFV("Lower cut", "");
1204  fCuts.Print();
1205 
1206  TString opt(option);
1207  opt.ToLower();
1208  if (!opt.Contains("nomax")) {
1209  PFV("Max weights", "");
1210 
1211  char ind[64];
1212  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
1213  ind[gROOT->GetDirLevel()] = '\0';
1214  for (UShort_t d=1; d<=3; d++) {
1215  UShort_t nr = (d == 1 ? 1 : 2);
1216  for (UShort_t q=0; q<nr; q++) {
1217  ind[gROOT->GetDirLevel()] = ' ';
1218  ind[gROOT->GetDirLevel()+1] = '\0';
1219  Char_t r = (q == 0 ? 'I' : 'O');
1220  std::cout << ind << " FMD" << d << r << ":";
1221  ind[gROOT->GetDirLevel()+1] = ' ';
1222  ind[gROOT->GetDirLevel()+2] = '\0';
1223 
1224  const TArrayI& a = (d == 1 ? fFMD1iMax :
1225  (d == 2 ? (r == 'I' ? fFMD2iMax : fFMD2oMax) :
1226  (r == 'I' ? fFMD3iMax : fFMD3oMax)));
1227  Int_t j = 0;
1228  for (Int_t i = 0; i < a.fN; i++) {
1229  if (a.fArray[i] < 1) continue;
1230  if (j % 6 == 0) std::cout << "\n " << ind;
1231  j++;
1232  std::cout << " " << std::setw(3) << i << ": " << a.fArray[i];
1233  }
1234  std::cout << std::endl;
1235  }
1236  }
1237  }
1238  gROOT->DecreaseDirLevel();
1239 }
1240 
1241 //====================================================================
1244  fList(0),
1245  // fEvsN(0),
1246  // fEvsM(0),
1247  // fEtaVsN(0),
1248  // fEtaVsM(0),
1249  fCorr(0),
1250  fDensity(0),
1251  fELossVsPoisson(0),
1252  fDiffELossPoisson(0),
1253  fELossVsPoissonOut(0),
1254  fDiffELossPoissonOut(0),
1255  fOutliers(0),
1256  fPoisson(),
1257  fELoss(0),
1258  fELossUsed(0),
1259  fMultCut(0),
1260  fTotal(0),
1261  fGood(0),
1262  fPhiAcc(0),
1263  fPhiBefore(0),
1264  fPhiAfter(0),
1265  fEtaBefore(0),
1266  fEtaAfter(0)
1267 {
1268  //
1269  // Default CTOR
1270  //
1271 }
1272 
1273 //____________________________________________________________________
1275  : AliForwardUtil::RingHistos(d,r),
1276  fList(0),
1277  // fEvsN(0),
1278  // fEvsM(0),
1279  // fEtaVsN(0),
1280  // fEtaVsM(0),
1281  fCorr(0),
1282  fDensity(0),
1283  fELossVsPoisson(0),
1284  fDiffELossPoisson(0),
1285  fELossVsPoissonOut(0),
1286  fDiffELossPoissonOut(0),
1287  fOutliers(0),
1288  fPoisson("ignored"),
1289  fELoss(0),
1290  fELossUsed(0),
1291  fMultCut(0),
1292  fTotal(0),
1293  fGood(0),
1294  fPhiAcc(0),
1295  fPhiBefore(0),
1296  fPhiAfter(0),
1297  fEtaBefore(0),
1298  fEtaAfter(0)
1299 {
1300  //
1301  // Constructor
1302  //
1303  // Parameters:
1304  // d detector
1305  // r ring
1306  //
1307 #if 0
1308  fEvsN = new TH2D("elossVsNnocorr",
1309  "#Delta E/#Delta E_{mip} vs uncorrected inclusive N_{ch}",
1310  250, -.5, 24.5, 251, -1.5, 24.5);
1311  fEvsN->SetXTitle("#Delta E/#Delta E_{mip}");
1312  fEvsN->SetYTitle("Inclusive N_{ch} (uncorrected)");
1313  fEvsN->Sumw2();
1314  fEvsN->SetDirectory(0);
1315 
1316  fEvsM = static_cast<TH2D*>(fEvsN->Clone("elossVsNcorr"));
1317  fEvsM->SetTitle("#Delta E/#Delta E_{mip} vs corrected inclusive N_{ch}");
1318  fEvsM->SetDirectory(0);
1319 
1320  fEtaVsN = new TProfile("etaVsNnocorr",
1321  "Average inclusive N_{ch} vs #eta (uncorrected)",
1322  200, -4, 6);
1323  fEtaVsN->SetXTitle("#eta");
1324  fEtaVsN->SetYTitle("#LT N_{ch,incl}#GT (uncorrected)");
1325  fEtaVsN->SetDirectory(0);
1326  fEtaVsN->SetLineColor(Color());
1327  fEtaVsN->SetFillColor(Color());
1328 
1329  fEtaVsM = static_cast<TProfile*>(fEtaVsN->Clone("etaVsNcorr"));
1330  fEtaVsM->SetTitle("Average inclusive N_{ch} vs #eta (corrected)");
1331  fEtaVsM->SetYTitle("#LT N_{ch,incl}#GT (corrected)");
1332  fEtaVsM->SetDirectory(0);
1333 #endif
1334 
1335  fCorr = new TProfile("corr", "Average correction", 200, -4, 6);
1336  fCorr->SetXTitle("#eta");
1337  fCorr->SetYTitle("#LT correction#GT");
1338  fCorr->SetDirectory(0);
1339  fCorr->SetLineColor(Color());
1340  fCorr->SetFillColor(Color());
1341 
1342  fDensity = new TH2D("inclDensity", "Inclusive N_{ch} density",
1343  200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
1344  0, 2*TMath::Pi());
1345  fDensity->SetDirectory(0);
1346  fDensity->Sumw2(); fDensity->SetMarkerColor(Color());
1347  fDensity->SetXTitle("#eta");
1348  fDensity->SetYTitle("#phi [radians]");
1349  fDensity->SetZTitle("Inclusive N_{ch} density");
1350 
1351  // --- Create increasing sized bins --------------------------------
1352  TArrayD bins;
1353  // bins, lowest order, higest order, return array
1354  const char* nchP = "N_{ch}^{Poisson}";
1355  const char* nchE = "N_{ch}^{#Delta}";
1356  AliForwardUtil::MakeLogScale(300, 0, 2, bins);
1357  fELossVsPoisson = new TH2D("elossVsPoisson",
1358  "N_{ch} from energy loss vs from Poisson",
1359  bins.GetSize()-1, bins.GetArray(),
1360  bins.GetSize()-1, bins.GetArray());
1361  fELossVsPoisson->SetDirectory(0);
1362  fELossVsPoisson->SetXTitle(nchE);
1363  fELossVsPoisson->SetYTitle(nchP);
1364  fELossVsPoisson->SetZTitle("Correlation");
1366  static_cast<TH2D*>(fELossVsPoisson
1367  ->Clone(Form("%sOutlier",
1368  fELossVsPoisson->GetName())));
1369  fELossVsPoissonOut->SetDirectory(0);
1370  fELossVsPoissonOut->SetMarkerStyle(20);
1371  fELossVsPoissonOut->SetMarkerSize(0.3);
1372  fELossVsPoissonOut->SetMarkerColor(kBlack);
1373  fELossVsPoissonOut->SetTitle(Form("%s for outliers",
1374  fELossVsPoisson->GetTitle()));
1375 
1376  fDiffELossPoisson = new TH1D("diffElossPoisson",
1377  Form("(%s-%s)/%s", nchP, nchE, nchE),
1378  100, -1, 1);
1379  fDiffELossPoisson->SetDirectory(0);
1380  fDiffELossPoisson->SetXTitle(fDiffELossPoisson->GetTitle());
1381  fDiffELossPoisson->SetYTitle("Frequency");
1382  fDiffELossPoisson->SetMarkerColor(Color());
1383  fDiffELossPoisson->SetFillColor(Color());
1384  fDiffELossPoisson->SetFillStyle(3001);
1385  fDiffELossPoisson->Sumw2();
1386 
1388  static_cast<TH1D*>(fDiffELossPoisson
1389  ->Clone(Form("%sOutlier",fDiffELossPoisson->GetName())));
1390  fDiffELossPoissonOut->SetDirectory(0);
1391  fDiffELossPoissonOut->SetTitle(Form("%s for outliers",
1392  fDiffELossPoisson->GetTitle()));
1393  fDiffELossPoissonOut->SetMarkerColor(Color()-2);
1394  fDiffELossPoissonOut->SetFillColor(Color()-2);
1395  fDiffELossPoissonOut->SetFillStyle(3002);
1396 
1397  fOutliers = new TH1D("outliers", "Fraction of outliers", 100, 0, 1);
1398  fOutliers->SetDirectory(0);
1399  fOutliers->SetXTitle("N_{outlier}/(N_{outlier}+N_{inside})");
1400  fOutliers->SetYTitle("#sum_{events}#sum_{bins}");
1401  fOutliers->SetFillColor(Color());
1402  fOutliers->SetFillStyle(3001);
1403  fOutliers->SetLineColor(kBlack);
1404 
1405  fELoss = new TH1D("eloss", "#Delta/#Delta_{mip} in all strips",
1406  640, -1, 15);
1407  fELoss->SetXTitle("#Delta/#Delta_{mip} (selected)");
1408  fELoss->SetYTitle("P(#Delta/#Delta_{mip})");
1409  fELoss->SetFillColor(Color()-2);
1410  fELoss->SetFillStyle(3003);
1411  fELoss->SetLineColor(kBlack);
1412  fELoss->SetLineStyle(2);
1413  fELoss->SetLineWidth(1);
1414  fELoss->SetDirectory(0);
1415 
1416  fELossUsed = static_cast<TH1D*>(fELoss->Clone("elossUsed"));
1417  fELossUsed->SetTitle("#Delta/#Delta_{mip} in used strips");
1418  fELossUsed->SetFillStyle(3002);
1419  fELossUsed->SetLineStyle(1);
1420  fELossUsed->SetDirectory(0);
1421 
1422  fPhiBefore = new TH1D("phiBefore", "#phi distribution (before recalc)",
1423  (r == 'I' || r == 'i' ? 20 : 40), 0, 2*TMath::Pi());
1424  fPhiBefore->SetDirectory(0);
1425  fPhiBefore->SetXTitle("#phi");
1426  fPhiBefore->SetYTitle("Events");
1427  fPhiBefore->SetMarkerColor(Color());
1428  fPhiBefore->SetLineColor(Color());
1429  fPhiBefore->SetFillColor(Color());
1430  fPhiBefore->SetFillStyle(3001);
1431  fPhiBefore->SetMarkerStyle(20);
1432 
1433  fPhiAfter = static_cast<TH1D*>(fPhiBefore->Clone("phiAfter"));
1434  fPhiAfter->SetTitle("#phi distribution (after re-calc)");
1435  fPhiAfter->SetDirectory(0);
1436 
1437  fEtaBefore = new TH1D("etaBefore", "#eta distribution (before recalc)",
1438  200, -4, 6);
1439  fEtaBefore->SetDirectory(0);
1440  fEtaBefore->SetXTitle("#eta");
1441  fEtaBefore->SetYTitle("Events");
1442  fEtaBefore->SetMarkerColor(Color());
1443  fEtaBefore->SetLineColor(Color());
1444  fEtaBefore->SetFillColor(Color());
1445  fEtaBefore->SetFillStyle(3001);
1446  fEtaBefore->SetMarkerStyle(20);
1447 
1448  fEtaAfter = static_cast<TH1D*>(fEtaBefore->Clone("etaAfter"));
1449  fEtaAfter->SetTitle("#eta distribution (after re-calc)");
1450  fEtaAfter->SetDirectory(0);
1451 
1452 }
1453 //____________________________________________________________________
1455  : AliForwardUtil::RingHistos(o),
1456  fList(o.fList),
1457  // fEvsN(o.fEvsN),
1458  // fEvsM(o.fEvsM),
1459  // fEtaVsN(o.fEtaVsN),
1460  // fEtaVsM(o.fEtaVsM),
1461  fCorr(o.fCorr),
1462  fDensity(o.fDensity),
1463  fELossVsPoisson(o.fELossVsPoisson),
1464  fDiffELossPoisson(o.fDiffELossPoisson),
1465  fELossVsPoissonOut(o.fELossVsPoissonOut),
1466  fDiffELossPoissonOut(o.fDiffELossPoissonOut),
1467  fOutliers(o.fOutliers),
1468  fPoisson(o.fPoisson),
1469  fELoss(o.fELoss),
1470  fELossUsed(o.fELossUsed),
1471  fMultCut(o.fMultCut),
1472  fTotal(o.fTotal),
1473  fGood(o.fGood),
1474  fPhiAcc(o.fPhiAcc),
1475  fPhiBefore(o.fPhiBefore),
1476  fPhiAfter(o.fPhiAfter),
1477  fEtaBefore(o.fEtaBefore),
1478  fEtaAfter(o.fEtaAfter)
1479 {
1480  //
1481  // Copy constructor
1482  //
1483  // Parameters:
1484  // o Object to copy from
1485  //
1486 }
1487 
1488 //____________________________________________________________________
1491 {
1492  //
1493  // Assignment operator
1494  //
1495  // Parameters:
1496  // o Object to assign from
1497  //
1498  // Return:
1499  // Reference to this
1500  //
1501  if (&o == this) return *this;
1503 
1504  // if (fEvsN) delete fEvsN;
1505  // if (fEvsM) delete fEvsM;
1506  // if (fEtaVsN) delete fEtaVsN;
1507  // if (fEtaVsM) delete fEtaVsM;
1508  if (fCorr) delete fCorr;
1509  if (fDensity) delete fDensity;
1510  if (fELossVsPoisson) delete fELossVsPoisson;
1511  if (fDiffELossPoisson) delete fDiffELossPoisson;
1512  if (fTotal) delete fTotal;
1513  if (fGood) delete fGood;
1514  if (fPhiAcc) delete fPhiAcc;
1515  if (fPhiBefore) delete fPhiBefore;
1516  if (fPhiAfter) delete fPhiAfter;
1517  if (fEtaBefore) delete fEtaBefore;
1518  if (fEtaAfter) delete fEtaAfter;
1519 
1520  // fEvsN = static_cast<TH2D*>(o.fEvsN->Clone());
1521  // fEvsM = static_cast<TH2D*>(o.fEvsM->Clone());
1522  // fEtaVsN = static_cast<TProfile*>(o.fEtaVsN->Clone());
1523  // fEtaVsM = static_cast<TProfile*>(o.fEtaVsM->Clone());
1524  fCorr = static_cast<TProfile*>(o.fCorr->Clone());
1525  fDensity = static_cast<TH2D*>(o.fDensity->Clone());
1526  fELossVsPoisson = static_cast<TH2D*>(o.fELossVsPoisson->Clone());
1527  fDiffELossPoisson = static_cast<TH1D*>(o.fDiffELossPoisson->Clone());
1528  fELossVsPoissonOut = static_cast<TH2D*>(o.fELossVsPoisson->Clone());
1529  fDiffELossPoissonOut = static_cast<TH1D*>(o.fDiffELossPoisson->Clone());
1530  fOutliers = static_cast<TH1D*>(o.fOutliers->Clone());
1531  fPoisson = o.fPoisson;
1532  fELoss = static_cast<TH1D*>(o.fELoss->Clone());
1533  fELossUsed = static_cast<TH1D*>(o.fELossUsed->Clone());
1534  fTotal = static_cast<TH1D*>(o.fTotal->Clone());
1535  fGood = static_cast<TH1D*>(o.fGood->Clone());
1536  fPhiAcc = static_cast<TH2D*>(o.fPhiAcc->Clone());
1537  fPhiBefore = static_cast<TH1D*>(o.fPhiBefore->Clone());
1538  fPhiAfter = static_cast<TH1D*>(o.fPhiAfter->Clone());
1539  fEtaBefore = static_cast<TH1D*>(o.fEtaBefore->Clone());
1540  fEtaAfter = static_cast<TH1D*>(o.fEtaAfter->Clone());
1541  return *this;
1542 }
1543 //____________________________________________________________________
1545 {
1546  //
1547  // Destructor
1548  //
1549 }
1550 
1551 
1552 //____________________________________________________________________
1553 void
1555 {
1556  // Initialize
1557  // This is called on first event
1558  fPoisson.Init(-1,-1);
1559  fTotal = new TH1D("total", "Total # of strips per #eta",
1560  eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax());
1561  fTotal->SetDirectory(0);
1562  fTotal->SetXTitle("#eta");
1563  fTotal->SetYTitle("# of strips");
1564  fGood = static_cast<TH1D*>(fTotal->Clone("good"));
1565  fGood->SetTitle("# of good strips per #eta");
1566  fGood->SetDirectory(0);
1567 
1568  fPhiAcc = new TH2D("phiAcc", "#phi acceptance vs Ip_{z}",
1569  eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
1570  10, -10, 10);
1571  fPhiAcc->SetXTitle("#eta");
1572  fPhiAcc->SetYTitle("v_{z} [cm]");
1573  fPhiAcc->SetZTitle("#phi acceptance");
1574  fPhiAcc->SetDirectory(0);
1575 
1576  if (fList) fList->Add(fPhiAcc);
1577 }
1578 
1579 //____________________________________________________________________
1580 void
1582 {
1583  //
1584  // Make output. This is called as part of SlaveBegin
1585  //
1586  // Parameters:
1587  // dir Where to put it
1588  //
1589  TList* d = DefineOutputList(dir);
1590  // d->Add(fEvsN);
1591  // d->Add(fEvsM);
1592  // d->Add(fEtaVsN);
1593  // d->Add(fEtaVsM);
1594  d->Add(fCorr);
1595  d->Add(fDensity);
1596  d->Add(fELossVsPoisson);
1597  d->Add(fELossVsPoissonOut);
1598  d->Add(fDiffELossPoisson);
1599  d->Add(fDiffELossPoissonOut);
1600  d->Add(fOutliers);
1601  fPoisson.Output(d);
1602  fPoisson.GetOccupancy()->SetFillColor(Color());
1603  fPoisson.GetMean()->SetFillColor(Color());
1604  fPoisson.GetOccupancy()->SetFillColor(Color());
1605  d->Add(fELoss);
1606  d->Add(fELossUsed);
1607  d->Add(fPhiBefore);
1608  d->Add(fPhiAfter);
1609  d->Add(fEtaBefore);
1610  d->Add(fEtaAfter);
1611 
1612  TAxis x(NStrip(), -.5, NStrip()-.5);
1613  TAxis y(NSector(), -.5, NSector()-.5);
1614  x.SetTitle("strip");
1615  y.SetTitle("sector");
1616  fPoisson.Define(x, y);
1617 
1618  d->Add(AliForwardUtil::MakeParameter("cut", fMultCut));
1619  fList = d;
1620 }
1621 
1622 //____________________________________________________________________
1623 void
1625 {
1626  //
1627  // Scale the histograms to the total number of events
1628  //
1629  // Parameters:
1630  // dir Where the output is
1631  // nEvents Number of events
1632  //
1633  TList* l = GetOutputList(dir);
1634  if (!l) return;
1635 
1636  TH2D* density = static_cast<TH2D*>(GetOutputHist(l,"inclDensity"));
1637  if (density) density->Scale(1./nEvents);
1638  fDensity = density;
1639 }
1640 
1641 //____________________________________________________________________
1642 //
1643 // EOF
1644 //
1645 
1646 
1647 
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:26
void Reset(const TH2D *base)
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 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
RingHistos * GetRingHistos(UShort_t d, Char_t r) const
AliForwardUtil::Histos fCache
#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)
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
#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
void Terminate(TList *dir, Int_t nEvents)
Float_t nEvents[nProd]
Definition: External.C:220
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
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
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
static const char * fgkFolderName
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()