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