AliPhysics  d565ceb (d565ceb)
TestELossDist.C
Go to the documentation of this file.
1 #include <iomanip>
2 #ifndef __CINT__
3 #include <iostream>
4 #include <TH1.h>
5 #include <TH2.h>
6 #include <TRandom.h>
7 #include <TCanvas.h>
8 #include <TArrayD.h>
9 #include <TStyle.h>
10 // #include <TFitResult.h>
11 #include <TGraphErrors.h>
12 #include <TF1.h>
13 #include <TMath.h>
14 #include <THStack.h>
15 #include <TList.h>
16 #include <TLatex.h>
17 #include <TProfile.h>
18 #include <TLegend.h>
19 #include <TLegendEntry.h>
20 #else
21 class TH1;
22 class TH2;
23 class TArrayD;
24 class TRandom;
25 class TCanvas;
26 class TF1;
27 #endif
28 
29 //___________________________________________________________________
30 static Double_t landauGaus1(Double_t* xp, Double_t* pp);
31 static Double_t landauGausN(Double_t* xp, Double_t* pp);
32 static Double_t landauGausI(Double_t* xp, Double_t* pp);
33 
34 //====================================================================
40 struct Function
41 {
46  enum {
48  kC,
52  kXi,
56  kN,
58  kA
59  };
61  static const Double_t fgkMPShift;
63  static const Double_t fgkInvRoot2Pi;
65  static const Double_t fgkConvNSigma;
67  static const Double_t fgkConvNSteps;
84  static Double_t Landau(Double_t x, Double_t delta, Double_t xi)
85  {
86  return TMath::Landau(x, delta - xi * fgkMPShift, xi);
87  }
119  Double_t xi, Double_t sigma)
120  {
121  Double_t deltap = delta - xi * fgkMPShift;
122  Double_t xlow = x - fgkConvNSigma * sigma;
123  Double_t xhigh = x + fgkConvNSigma * sigma;
124  Double_t step = (xhigh - xlow) / fgkConvNSteps;
125  Double_t sum = 0;
126 
127  for (Int_t i = 0; i < fgkConvNSteps / 2; i++) {
128  Double_t x1 = xlow + (i - .5) * step;
129  Double_t x2 = xhigh - (i - .5) * step;
130 
131  Double_t s1 =
132  TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma);
133  Double_t s2 =
134  TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma);
135  sum += s1 + s2;
136  }
137  Double_t ret = step * sum * fgkInvRoot2Pi / sigma;
138  return ret;
139  }
161  Double_t xi, Double_t sigma)
162  {
163  if (i <= 1) return LandauGaus(x, delta, xi, sigma);
164  Double_t di = i;
165  Double_t deltai = i * (delta + xi * TMath::Log(di));
166  Double_t xii = i * xi;
167  Double_t sigmai = TMath::Sqrt(di) * sigma;
168 
169  if (sigmai < 1e-10) return Landau(x, deltai, xii);
170  Double_t ret = LandauGaus(x, deltai, xii, sigmai);;
171  // Info("ILandauGaus", "Fi(%f;%d,%f,%f,%f)->%f",
172  // x, i, deltai, xii, sigmai, ret);
173  return ret;
174  }
207  UShort_t ipar, Double_t dPar,
208  Double_t delta, Double_t xi,
209  Double_t sigma)
210  {
211  if (dPar == 0) return 0;
212  Double_t dp = dPar;
213  Double_t d2 = dPar / 2;
214  Double_t deltaI = i * (delta + xi * TMath::Log(i));
215  Double_t xiI = i * xi;
216  Double_t si = TMath::Sqrt(Double_t(i));
217  Double_t sigmaI = si*sigma;
218  Double_t y1 = 0;
219  Double_t y2 = 0;
220  Double_t y3 = 0;
221  Double_t y4 = 0;
222  switch (ipar) {
223  case 0:
224  y1 = ILandauGaus(i, x, deltaI+i*dp, xiI, sigmaI);
225  y2 = ILandauGaus(i, x, deltaI+i*d2, xiI, sigmaI);
226  y3 = ILandauGaus(i, x, deltaI-i*d2, xiI, sigmaI);
227  y4 = ILandauGaus(i, x, deltaI-i*dp, xiI, sigmaI);
228  break;
229  case 1:
230  y1 = ILandauGaus(i, x, deltaI, xiI+i*dp, sigmaI);
231  y2 = ILandauGaus(i, x, deltaI, xiI+i*d2, sigmaI);
232  y3 = ILandauGaus(i, x, deltaI, xiI-i*d2, sigmaI);
233  y4 = ILandauGaus(i, x, deltaI, xiI-i*dp, sigmaI);
234  break;
235  case 2:
236  y1 = ILandauGaus(i, x, deltaI, xiI, sigmaI+si*dp);
237  y2 = ILandauGaus(i, x, deltaI, xiI, sigmaI+si*d2);
238  y3 = ILandauGaus(i, x, deltaI, xiI, sigmaI-si*d2);
239  y4 = ILandauGaus(i, x, deltaI, xiI, sigmaI-si*dp);
240  break;
241  default:
242  return 0;
243  }
244 
245  Double_t d0 = y1 - y4;
246  Double_t d1 = 2 * (y2 - y3);
247 
248  Double_t g = 1/(2*dp) * (4*d1 - d0) / 3;
249 
250  return g;
251  }
278  Double_t sigma, Int_t n, Double_t* a)
279  {
280  Double_t res = LandauGaus(x, delta, xi, sigma);
281  for (Int_t i = 2; i <= n; i++)
282  res += a[i-2] * ILandauGaus(i, x, delta, xi, sigma);
283  return res;
284  }
305  Double_t sigma, Int_t n, Double_t* a)
306  {
307  Double_t num = LandauGaus(x, delta, xi, sigma);
308  Double_t den = num;
309  for (Int_t i = 2; i <= n; i++) {
310  Double_t f = ILandauGaus(i, x, delta, xi, sigma);
311  num += a[i-2] * f * i;
312  den += a[i-2] * f;
313  }
314  if (den < 1e-4) return 0;
315  return num / den;
316  }
339  Double_t DEstimatePar(Double_t x, Int_t ipar, Double_t delta, Double_t xi,
341  Double_t dPar=0.001)
342  {
343  Double_t dp = dPar;
344  Double_t d2 = dPar / 2;
345  Double_t y1 = 0;
346  Double_t y2 = 0;
347  Double_t y3 = 0;
348  Double_t y4 = 0;
349  switch (ipar) {
350  case Function::kC:
351  case Function::kN:
352  return 0;
353  case Function::kDelta:
354  y1 = NEstimate(x, delta+dp, xi, sigma, n, a);
355  y2 = NEstimate(x, delta+d2, xi, sigma, n, a);
356  y3 = NEstimate(x, delta-d2, xi, sigma, n, a);
357  y4 = NEstimate(x, delta-dp, xi, sigma, n, a);
358  break;
359  case Function::kXi:
360  y1 = NEstimate(x, delta, xi+dp, sigma, n, a);
361  y2 = NEstimate(x, delta, xi+d2, sigma, n, a);
362  y3 = NEstimate(x, delta, xi-d2, sigma, n, a);
363  y4 = NEstimate(x, delta, xi-dp, sigma, n, a);
364  break;
365  case Function::kSigma:
366  y1 = NEstimate(x, delta, xi, sigma+dp, n, a);
367  y2 = NEstimate(x, delta, xi, sigma+d2, n, a);
368  y3 = NEstimate(x, delta, xi, sigma-d2, n, a);
369  y4 = NEstimate(x, delta, xi, sigma-dp, n, a);
370  break;
371  default: {
372  Int_t j = ipar-kA;
373  if (j+1 > n) return 0;
374  Double_t aa = a[j];
375  a[j] = aa+dp; y1 = NEstimate(x, delta, xi, sigma, n, a);
376  a[j] = aa+d2; y2 = NEstimate(x, delta, xi, sigma, n, a);
377  a[j] = aa-d2; y3 = NEstimate(x, delta, xi, sigma, n, a);
378  a[j] = aa-dp; y4 = NEstimate(x, delta, xi, sigma, n, a);
379  a[j] = aa;
380  }
381  break;
382  }
383  Double_t d0 = y1 - y4;
384  Double_t d1 = 2 * (y2 - y3);
385 
386  Double_t g = 1/(2*dp) * (4*d1 - d0) / 3;
387  return g;
388  }
403  static TF1* MakeFunc(Int_t n, Double_t c, Double_t delta,
404  Double_t xi, Double_t sigma, Double_t* a,
405  Double_t xmin, Double_t xmax)
406  {
407  Int_t nPar = kN + n;
408  TF1* f = new TF1(Form("landGaus%d",n), &landauGausN, xmin, xmax, nPar);
409  f->SetNpx(500);
410  f->SetParName(kC, "C");
411  f->SetParName(kDelta, "#Delta_{p}");
412  f->SetParName(kXi, "#xi");
413  f->SetParName(kSigma, "#sigma");
414  f->SetParName(kN, "N");
415  f->SetParameter(kC, c);
416  f->SetParameter(kDelta, delta);
417  f->SetParameter(kXi, xi);
418  f->SetParameter(kSigma, sigma);
419  f->FixParameter(kN, n);
420  f->SetParLimits(kDelta, xmin, 4);
421  f->SetParLimits(kXi, 0.00, 4);
422  f->SetParLimits(kSigma, 0.01, 0.11);
423 
424  for (Int_t i = 2; i <= n; i++) {
425  Int_t j = i-2;
426  f->SetParName(kA+j, Form("a_{%d}", i));
427  f->SetParameter(kA+j, a[j]);
428  f->SetParLimits(kA+j, 0, 1);
429  }
430  return f;
431  }
446  static TF1* MakeIFunc(Int_t i, Double_t c, Double_t delta,
447  Double_t xi, Double_t sigma,
448  Double_t xmin, Double_t xmax)
449  {
450  Int_t nPar = 5;
451  TF1* f = new TF1(Form("landGausI%d",i), &landauGausI, xmin, xmax, nPar);
452  f->SetNpx(100);
453  f->SetParName(kC, "C");
454  f->SetParName(kDelta, "#Delta_{p}");
455  f->SetParName(kXi, "#xi");
456  f->SetParName(kSigma, "#sigma");
457  f->SetParName(kN, "N");
458  f->SetParameter(kC, c);
459  f->SetParameter(kDelta, delta);
460  f->SetParameter(kXi, xi);
461  f->SetParameter(kSigma, sigma);
462  f->FixParameter(kN, i);
463  return f;
464  }
465 
466  // --- Object code -------------------------------------------------
481  Function(Int_t n, Double_t c, Double_t delta, Double_t xi, Double_t sigma,
482  Double_t* a, Double_t xmin, Double_t xmax)
483  : fF(MakeFunc(n, c, delta, xi, sigma, a, xmin, xmax))
484  {}
485  Function(TF1* f) : fF(f) {}
491  TF1* GetF1() { return fF; }
504  {
505  return fF->Eval(x);
506  }
530  {
531  Double_t delta = GetDelta();
532  Double_t xi = GetXi();
533  Double_t sigma = GetSigma();
534  Double_t dFdDelta2 = 0;
535  Double_t dFdXi2 = 0;
536  Double_t dFdSigma2 = 0;
537  e = 0;
538  for (Int_t i = 1; i <= GetN(); i++) {
539  Double_t a = GetA(i);
540  Double_t dDelta = a*IdLandauGausdPar(i, x, 1, 0.001, delta, xi, sigma);
541  Double_t dXi = a*IdLandauGausdPar(i, x, 2, 0.001, delta, xi, sigma);
542  Double_t dSigma = a*IdLandauGausdPar(i, x, 3, 0.001, delta, xi, sigma);
543  Double_t dFda = a*ILandauGaus(i, x, delta, xi, sigma);
544  dFdDelta2 += dDelta * dDelta;
545  dFdXi2 += dXi * dXi;
546  dFdSigma2 += dSigma * dSigma;
547  e += TMath::Power(dFda * GetEA(i),2);
548  }
549  Double_t edelta = GetEDelta();
550  Double_t exi = GetEXi();
551  Double_t esigma = GetESigma();
552  e += (dFdDelta2 * edelta * edelta +
553  dFdXi2 * exi * exi +
554  dFdSigma2 * esigma * esigma);
555  e = TMath::Sqrt(e);
556  return fF->Eval(x);
557  }
558 
577  {
578  UShort_t n = (maxN < 0 ? GetN() : TMath::Min(UShort_t(maxN), GetN()));
579  Double_t delta = GetDelta();
580  Double_t xi = GetXi();
581  Double_t sigma = GetSigma();
582  return NEstimate(x, delta, xi, sigma, n, GetAs());
583  }
616  {
617 
618  UShort_t n = (maxN < 0 ? GetN() : TMath::Min(UShort_t(maxN), GetN()));
619  Double_t delta = GetDelta();
620  Double_t xi = GetXi();
621  Double_t sigma = GetSigma();
622  Double_t* a = GetAs();
623  Double_t dDelta = (DEstimatePar(x,Function::kDelta,delta,xi,sigma,n,a,\
624  0.01 * GetEDelta()) * GetEDelta());
625  Double_t dXi = (DEstimatePar(x,Function::kXi, delta,xi,sigma,n,a,
626  0.01 * GetEXi()) * GetEXi());
627  Double_t dSigma = (DEstimatePar(x,Function::kSigma,delta,xi,sigma,n,a,
628  0.01 * GetESigma()) * GetESigma());
629  e = dDelta * dDelta + dXi * dXi + dSigma * dSigma;
630  for (Int_t i = 2; i <= n; i++) {
631  Double_t dAi = (DEstimatePar(x, Function::kA+i-2,delta,xi,sigma,n,a,
632  0.01 * GetEA(i)) * GetEA(i));
633  e += dAi * dAi;
634  }
635  e = TMath::Sqrt(e);
636  return NEstimate(x, GetDelta(), GetXi(), GetSigma(), n, GetAs());
637  }
655  {
656  UShort_t n = (maxN < 0 ? GetN() : TMath::Min(UShort_t(maxN), GetN()));
657  TArrayD p(n);
658  Double_t delta = GetDelta();
659  Double_t xi = GetXi();
660  Double_t sigma = GetSigma();
661 
662  if (Evaluate(x) < 1e-4) return 0;
663 
664  // std::cout << "RandomEstimateNParticles(" << x << "):";
665  for (Int_t i = 1; i <= n; i++) {
666  Double_t a = GetA(i);
667  Double_t f = ILandauGaus(i, x, delta, xi, sigma);
668  p[i-1] = a * f;
669  // std::cout << p[i-1] << ",";
670  }
671  Double_t r = gRandom->Uniform(p.GetSum());
672  Double_t l = 0;
673  // std::cout << "s=" << p.GetSum() << ",r=" << r << ",";
674  for (Int_t i = 1; i <= n; i++) {
675  if (r > l && r <= l + p[i-1]) {
676  // std::cout << "l=" << l << ",l+p=" << l+p[i-1] << "-> " << i
677  // << std::endl;
678  return i;
679  }
680  l += p[i-1];
681  }
682  return 0;
683  }
689  void Draw(Option_t* option="") { fF->Draw(option); }
691  Double_t GetC() const { return fF->GetParameter(kC); }
693  Double_t GetDelta() const { return fF->GetParameter(kDelta); }
695  Double_t GetXi() const { return fF->GetParameter(kXi); }
697  Double_t GetSigma() const { return fF->GetParameter(kSigma); }
699  UShort_t GetN() const { return fF->GetParameter(kN); }
701  Double_t GetA(Int_t i) const { return i == 1 ? 1 : fF->GetParameter(kA+i-2);}
703  Double_t* GetAs() const { return &(fF->GetParameters()[kA]);}
705  Double_t GetEC() const { return fF->GetParError(kC); }
707  Double_t GetEDelta() const { return fF->GetParError(kDelta); }
709  Double_t GetEXi() const { return fF->GetParError(kXi); }
711  Double_t GetESigma() const { return fF->GetParError(kSigma); }
713  UShort_t GetEN() const { return 0; }
715  Double_t GetEA(Int_t i) const { return i == 1 ? 0 : fF->GetParError(kA+i-2);}
717  Double_t* GetEAs() const { return &(fF->GetParErrors()[kA]);}
719  TF1* fF;
720 };
721 #ifndef __CINT__
722 const Double_t Function::fgkMPShift = -0.22278298;
723 const Double_t Function::fgkInvRoot2Pi = 1. / TMath::Sqrt(2*TMath::Pi());
726 #endif
727 
728 //====================================================================
735 struct Fitter
736 {
737  // --- Object code ------------------------------------------------
738  Fitter(Int_t nMinus)
739  : fNMinus(nMinus)
740  {
741  }
749  TF1* FitFirst(TH1* dist)
750  {
751  fFunctions.Clear();
752 
753  // Get upper edge
754  Int_t midBin = dist->GetMaximumBin();
755  Double_t midE = dist->GetBinLowEdge(midBin);
756 
757  // Get low edge
758  Int_t minBin = midBin - fNMinus;
759  Double_t minE = dist->GetBinCenter(minBin);
760  Double_t maxE = dist->GetBinCenter(midBin+2*fNMinus);
761  Double_t a[] = { 0 };
762  TF1* f = Function::MakeFunc(1, 0.5, midE, 0.07, 0.1, a, minE, maxE);
763 
764  // Do the fit
765  dist->Fit(f, "RNQS", "", minE, maxE);
766  fFunctions.AddAtAndExpand(f, 0);
767  return f;
768  }
777  Function* FitN(TH1* dist, Int_t n)
778  {
779  if (n == 1) return new Function(FitFirst(dist));
780  TF1* f1 = static_cast<TF1*>(fFunctions.At(0));
781  if (!f1) {
782  f1 = FitFirst(dist);
783  if (!f1) {
784  ::Warning("FitN", "First fit missing");
785  return 0;
786  }
787  }
788 
789  Double_t delta1 = f1->GetParameter(Function::kDelta);
790  Double_t xi1 = f1->GetParameter(Function::kXi);
791  Double_t maxEi = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
792  Double_t minE = f1->GetXmin();
793 
794  TArrayD a(n-1);
795  for (UShort_t i = 2; i <= n; i++) a.fArray[i-2] = (n==2? 0.05 : 0.000001);
796 
797  TF1* fn = Function::MakeFunc(n, f1->GetParameter(Function::kC),
798  delta1, xi1,
799  f1->GetParameter(Function::kSigma),
800  a.fArray, minE, maxEi);
801 
802  dist->Fit(fn, "RSNQ", "", minE, maxEi);
803  dist->Fit(fn, "RSNQM", "", minE, maxEi);
804  fFunctions.AddAtAndExpand(fn, n-1);
805 
806  return new Function(fn);
807  }
810 };
811 
812 //====================================================================
824 {
825  Double_t x = xp[0];
826  Double_t c = pp[Function::kC];
827  Double_t delta = pp[Function::kDelta];
828  Double_t xi = pp[Function::kXi];
830 
831  return c * Function::LandauGaus(x, delta, xi, sigma);
832 }
844 {
845  Double_t x = xp[0];
846  Double_t c = pp[Function::kC];
847  Double_t delta = pp[Function::kDelta];
848  Double_t xi = pp[Function::kXi];
850  Int_t n = Int_t(pp[Function::kN]);
851  Double_t* a = &pp[Function::kA];
852 
853  return c * Function::NLandauGaus(x, delta, xi, sigma, n, a);
854 }
866 {
867  Double_t x = xp[0];
868  Double_t c = pp[Function::kC];
869  Double_t delta = pp[Function::kDelta];
870  Double_t xi = pp[Function::kXi];
872  Int_t i = Int_t(pp[Function::kN]);
873 
874  return c * Function::ILandauGaus(i, x, delta, xi, sigma);
875 }
876 
877 //====================================================================
884 struct Generator
885 {
894  TH1* fN;
895  TArrayD fRP; // Raw probabilities
896  TArrayD fP; // Normalized probabilities
897  TF1* fF;
899 
906  Generator(Double_t delta=.55, Double_t xi=0.06, Double_t sigma=0.02,
907  Int_t maxN=5, const Double_t* a=0)
908  : fDelta(delta), fXi(xi), fSigma(sigma), fMaxN(maxN),
909  fSingle(0), fSummed(0), fSimple(0), fNSum(0), fRP(), fP(fMaxN), fF(0),
910  fNObs(0)
911  {
912  fMaxN = TMath::Min(fMaxN, 5);
913  const Double_t ps[] = { 1, .05, .02, .005, .002 };
914  fRP.Set(5, ps);
915  if (a) {
916  for (Int_t i = 0; i < maxN; i++) fRP[i] = a[i];
917  }
918 
919  Double_t sum = 0;
920  for (Int_t i = 0; i < fMaxN; i++) sum += fRP[i];
921  for (Int_t i = 0; i < fP.fN; i++) {
922  fP[i] = fRP[i] / sum;
923  if (i != 0) fP[i] += fP[i-1];
924  }
925 
926  fSingle = new TH1D("single", "Single signals", 500, 0, 10);
927  fSingle->SetXTitle("Signal");
928  fSingle->SetYTitle("Events");
929  fSingle->SetFillColor(kRed+1);
930  fSingle->SetMarkerColor(kRed+1);
931  fSingle->SetLineColor(kRed+1);
932  fSingle->SetFillStyle(3001);
933  fSingle->SetMarkerStyle(21);
934  fSingle->SetDirectory(0);
935  fSingle->Sumw2();
936 
937  fSummed = static_cast<TH1D*>(fSingle->Clone("summed"));
938  fSummed->SetTitle("Summed signals");
939  fSummed->SetFillColor(kBlue+1);
940  fSummed->SetMarkerColor(kBlue+1);
941  fSummed->SetLineColor(kBlue+1);
942  fSummed->SetDirectory(0);
943 
944  fSimple = static_cast<TH1D*>(fSingle->Clone("summed"));
945  fSimple->SetTitle("Summed signals");
946  fSimple->SetFillColor(kGreen+1);
947  fSimple->SetMarkerColor(kGreen+1);
948  fSimple->SetLineColor(kGreen+1);
949  fSimple->SetDirectory(0);
950 
951  fNSum = new TH2D("single", "Observation signals",
952  fMaxN, .5, fMaxN+.5, 100, 0, 10);
953  fNSum->SetMarkerColor(kMagenta+1);
954  fNSum->SetLineColor(kMagenta+1);
955  fNSum->SetFillColor(kMagenta+1);
956  fNSum->SetYTitle("Signal");
957  fNSum->SetXTitle("N");
958  fNSum->SetZTitle("Events");
959  fNSum->SetDirectory(0);
960  fNSum->Sumw2();
961 
962  fN = new TH1D("n", "Number of particles", fMaxN, .5, fMaxN+.5);
963  fN->SetXTitle("N");
964  fN->SetYTitle("Events");
965  fN->SetMarkerColor(kMagenta+1);
966  fN->SetLineColor(kMagenta+1);
967  fN->SetFillColor(kMagenta+1);
968  fN->SetMarkerStyle(22);
969  fN->SetFillStyle(3001);
970  fN->SetDirectory(0);
971 
972  std::cout << "Probabilities:" << std::setprecision(4)
973  << "\n ";
974  for (Int_t i = 0; i < fRP.fN; i++)
975  std::cout << std::setw(7) << fRP[i] << " ";
976  std::cout << " = " << fRP.GetSum() << "\n ";
977  for (Int_t i = 0; i < fP.fN; i++)
978  std::cout << std::setw(7) << fP[i] << " ";
979  std::cout << std::endl;
980 
981  Double_t aa[] = { 0 };
982  fF = Function::MakeFunc(1, 1, fDelta, fXi, fSigma, aa, 0, 10);
983  }
988  void Clear()
989  {
990  fSingle->Reset();
991  fSummed->Reset();
992  fNSum->Reset();
993  fN->Reset();
994  }
1002  {
1003  Double_t rndm = gRandom->Uniform();
1004  Int_t ret = 0;
1005  if ( rndm < fP[0]) ret = 1;
1006  else if (rndm >= fP[0] && rndm < fP[1]) ret = 2;
1007  else if (rndm >= fP[1] && rndm < fP[2]) ret = 3;
1008  else if (rndm >= fP[2] && rndm < fP[3]) ret = 4;
1009  else if (rndm >= fP[3] && rndm < fP[4]) ret = 5;
1010  // Printf("%f -> %d", rndm, ret);
1011  fN->Fill(ret);
1012  return ret;
1013  }
1018  {
1019  Double_t ret = 0;
1020  if (fF) ret = fF->GetRandom();
1021  else {
1022  Double_t rmpv = gRandom->Gaus(fDelta, fSigma);
1023  ret = gRandom->Landau(rmpv, fXi);
1024  }
1025  fSingle->Fill(ret);
1026  fSimple->Fill(gRandom->Landau(fDelta - fXi * Function::fgkMPShift, fXi));
1027  return ret;
1028  }
1035  {
1036  Double_t ret = 0;
1037  for (Int_t i = 0; i < n; i++) ret += Generate1Signal();
1038  fSummed->Fill(ret);
1039  fNSum->Fill(n, ret);
1040  return ret;
1041  }
1048  void Generate(Int_t nObs=1000, Bool_t reset=true)
1049  {
1050  if (reset) Clear();
1051  fNObs = nObs;
1052 
1053  for (Int_t i = 0; i < nObs; i++) {
1054  // if (((i+1) % (nObs/10)) == 0)
1055  // std::cout << "Event " << std::setw(6) << i << std::endl;
1056  Int_t n = GenerateN();
1057  GenerateNSignal(n);
1058  }
1059  Double_t m = fSingle->GetMaximum();
1060  fSingle->Scale(1. / m);
1061  m = fSummed->GetMaximum();
1062  fSummed->Scale(1. / m);
1063  m = fSimple->GetMaximum();
1064  fSimple->Scale(1. / m);
1065 
1066  std::cout << "Resulting probabilities:\n ";
1067  for (Int_t i = 1; i <= fN->GetNbinsX(); i++)
1068  std::cout << std::setprecision(4) << std::setw(7)
1069  << fN->GetBinContent(i) << " ";
1070  std::cout << std::endl;
1071  }
1079  void PrintParameter(Double_t input, TF1* f, Int_t i)
1080  {
1081  Double_t val = 0, err = 0;
1082  if (i < f->GetNpar()) {
1083  val = f->GetParameter(i);
1084  err = f->GetParError(i);
1085  }
1086  std::cout << " " << std::setw(16) << f->GetParName(i) << ": "
1087  << std::fixed
1088  << std::setprecision(4) << std::setw(6) << input << " -> "
1089  << std::setprecision(4) << std::setw(6) << val << " +/- "
1090  << std::setprecision(5) << std::setw(7) << err << " [delta:"
1091  << std::setw(3) << int(100 * (input-val) / input) << "%,err:"
1092  << std::setw(3) << int(100 * err / val) << "%]"
1093  << std::endl;
1094  }
1100  void PrintResult(TF1* f)
1101  {
1102  Double_t chi2 = f->GetChisquare();
1103  Int_t ndf = f->GetNDF();
1104  std::cout << "Chi^2/NDF = " << chi2 << "/" << ndf << "="
1105  << chi2/ndf << std::endl;
1106  PrintParameter(1, f, Function::kC);
1107  PrintParameter(fDelta, f, Function::kDelta);
1108  PrintParameter(fXi, f, Function::kXi);
1109  PrintParameter(fSigma, f, Function::kSigma);
1110  PrintParameter(fMaxN, f, Function::kN);
1111  for (Int_t i = 0; i < fMaxN-1; i++)
1112  PrintParameter(fRP[i+1], f, Function::kA+i);
1113  }
1123  void DrawParameter(Double_t x, Double_t y, Double_t input, TF1* f, Int_t i)
1124  {
1125  Double_t val = 0, err = 0;
1126  if (i < f->GetNpar()) {
1127  val = f->GetParameter(i);
1128  err = f->GetParError(i);
1129  }
1130 
1131  TLatex* ltx = new TLatex(x, y, f->GetParName(i));
1132  ltx->SetTextSize(.04);
1133  ltx->SetTextFont(132);
1134  ltx->SetTextAlign(11);
1135  ltx->SetNDC();
1136  ltx->Draw();
1137 
1138  ltx->DrawLatex(x+.05, y, Form("%5.3f", input));
1139  ltx->DrawLatex(x+.14, y, "#rightarrow");
1140  ltx->SetTextAlign(21);
1141  ltx->DrawLatex(x+.3, y, Form("%5.3f #pm %6.4f", val, err));
1142  ltx->SetTextAlign(11);
1143  ltx->DrawLatex(x+.48, y, Form("[#Delta=%3d%%, #delta=%3d%%]",
1144  int(100 * (input-val)/input),
1145  int(100*err/val)));
1146  }
1154  void DrawResult(Double_t x, Double_t y, TF1* f)
1155  {
1156  Double_t chi2 = f->GetChisquare();
1157  Int_t ndf = f->GetNDF();
1158  TLatex* ltx = new TLatex(x, y, Form("#chi^{2}/NDF=%5.2f/%3d=%6.3f",
1159  chi2, ndf, chi2/ndf));
1160  ltx->SetTextSize(.04);
1161  ltx->SetNDC();
1162  ltx->SetTextFont(132);
1163  ltx->Draw();
1164 
1165  x += .05;
1166  y -= .035; DrawParameter(x, y, 1, f, Function::kC);
1167  y -= .035; DrawParameter(x, y, fDelta, f, Function::kDelta);
1168  y -= .035; DrawParameter(x, y, fXi, f, Function::kXi);
1169  y -= .035; DrawParameter(x, y, fSigma, f, Function::kSigma);
1170  y -= .035; DrawParameter(x, y, fMaxN, f, Function::kN);
1171  for (Int_t i = 0; i < fMaxN-1; i++) {
1172  y -= .035; DrawParameter(x, y, fRP[i+1], f, Function::kA+i);
1173  }
1174  }
1179  void Draw(const char* type="png")
1180  {
1181  gStyle->SetOptTitle(0);
1182  gStyle->SetOptStat(0);
1183  gStyle->SetPalette(1);
1184 
1185  TCanvas* c = new TCanvas("c", "Generator distributions", 1200, 1000);
1186  c->SetFillColor(0);
1187  c->SetBorderSize(0);
1188  c->SetBorderMode(0);
1189  c->SetTopMargin(0.01);
1190  c->SetRightMargin(0.01);
1191 
1192  c->Divide(2, 2);
1193 
1194  // --- First pad: Data (multi, single, landau) + Fit -------------
1195  TVirtualPad* p = c->cd(1);
1196  p->SetLogy();
1197  p->SetBorderSize(0);
1198  p->SetBorderMode(0);
1199  p->SetTopMargin(0.01);
1200  p->SetRightMargin(0.01);
1201  fSummed->Draw();
1202  fSingle->Draw("same");
1203  fSimple->Draw("same");
1204 
1205  if (fF) {
1206  Double_t max = fF->GetMaximum();
1207  Double_t a[] = { 0 };
1208  TF1* fcopy = Function::MakeFunc(1,
1209  fF->GetParameter(Function::kC) /max,
1210  fF->GetParameter(Function::kDelta),
1211  fF->GetParameter(Function::kXi),
1212  fF->GetParameter(Function::kSigma),
1213  a, 0, 10);
1214  fcopy->SetLineColor(2);
1215  fcopy->SetLineStyle(2);
1216  fcopy->Draw("same");
1217  }
1218 
1219  Fitter* fitter = new Fitter(4);
1220  Function* fit = fitter->FitN(fSummed, 5);
1221  TF1* f = fit->GetF1();
1222  f->Draw("same");
1223  f->SetLineColor(kBlack);
1224  TF1* fst = static_cast<TF1*>(fitter->fFunctions.At(0));
1225  fst->SetLineWidth(3);
1226  fst->SetLineStyle(2);
1227  fst->SetLineColor(kMagenta);
1228  fst->Draw("same");
1229  PrintResult(f);
1230  DrawResult(.25, .95, f);
1231  TGraphErrors* gr = new TGraphErrors(fSummed->GetNbinsX());
1232  gr->SetName("fitError");
1233  gr->SetTitle("Fit w/errors");
1234  gr->SetFillColor(kYellow+1);
1235  gr->SetFillStyle(3001);
1236  for (Int_t j = 1; j <= fSummed->GetNbinsX(); j++) {
1237  Double_t x = fSummed->GetBinCenter(j);
1238  Double_t e = 0;
1239  Double_t y = fit->Evaluate(x, e);
1240  gr->SetPoint(j-1, x, y);
1241  gr->SetPointError(j-1, 0, e);
1242  }
1243  gr->Draw("same l e4");
1244 
1245  TLegend* l = new TLegend(.15, .11, .5, .4);
1246  l->SetBorderSize(0);
1247  l->SetFillColor(0);
1248  l->SetFillStyle(0);
1249  l->SetTextFont(132);
1250  l->AddEntry(fSimple, "L: Single-L", "lp");
1251  l->AddEntry(fSingle, "LG: Single-L#otimesG", "lp");
1252  l->AddEntry(fSummed,"NLG: Multi-L#otimesG", "lp");
1253  l->AddEntry(gr, "f_{N}: Fit to NLG");
1254  l->Draw();
1255 
1256  // --- Second pad: Data (1, 2, 3, ...) --------------------------
1257  p = c->cd(2);
1258  p->SetLogy();
1259  p->SetBorderSize(0);
1260  p->SetBorderMode(0);
1261  p->SetTopMargin(0.01);
1262  p->SetRightMargin(0.01);
1263  p->SetGridx();
1264  THStack* stack = new THStack(fNSum, "y");
1265  TIter next(stack->GetHists());
1266  TH1* h = 0;
1267  Int_t i = 2;
1268  Double_t max = 0;
1269  TLegend* l2 = new TLegend(.6, .5, .95, .95);
1270  l2->SetBorderSize(0);
1271  l2->SetFillStyle(0);
1272  l2->SetFillColor(0);
1273  l2->SetTextFont(132);
1274  while ((h = static_cast<TH1*>(next()))) {
1275  if (i == 2) max = h->GetMaximum();
1276  h->SetLineColor(i);
1277  h->SetFillColor(i);
1278  h->SetFillStyle(3001);
1279  h->Scale(1/max);
1280  l2->AddEntry(h, Form("%d particles", i - 1), "f");
1281  i++;
1282  }
1283  stack->Draw("hist nostack");
1284  stack->GetHistogram()->SetXTitle("Signal");
1285  stack->GetHistogram()->SetYTitle("Events");
1286  i = 2;
1287  for (Int_t j = 1; j <= 5; j++) {
1288  TF1* fi = Function::MakeIFunc(j, fRP[j-1], fDelta, fXi, fSigma, 0, 10);
1289  if (j == 1) max = fi->GetMaximum();
1290  fi->SetParameter(0, fRP[j-1]/max);
1291  fi->SetLineColor(i);
1292  fi->SetLineWidth(3);
1293  fi->SetLineStyle(2);
1294  fi->Draw("same");
1295 
1296  TF1* fj = Function::MakeIFunc(j, fit->GetC() * fit->GetA(j),
1297  fit->GetDelta(), fit->GetXi(),
1298  fit->GetSigma(), 0, 10);
1299  fj->SetLineColor(i);
1300  fj->SetLineWidth(2);
1301  fj->Draw("same");
1302 
1303  std::cout << "Component " << j << " scale=" << fi->GetParameter(0)
1304  << " (" << fRP[j-1] << ")" << " " << fj->GetParameter(0)
1305  << " (" << fit->GetC() << "*" << fit->GetA(j) << ")"
1306  << std::endl;
1307  i++;
1308  }
1309  TLegendEntry* ee = l2->AddEntry("", "Input f_{i}", "l");
1310  ee->SetLineStyle(2);
1311  ee->SetLineWidth(3);
1312  l2->AddEntry("", "Fit f_{i}", "l");
1313  l2->Draw();
1314  // fNSum->Draw("lego2");
1315 
1316  // --- Third pad: Mean multiplicity ------------------------------
1317  TH1* nhist1 = new TH1F("nhist1", "NHist", 5, 0.5, 5.5);
1318  nhist1->SetFillColor(kCyan+1);
1319  nhist1->SetMarkerColor(kCyan+1);
1320  nhist1->SetLineColor(kCyan+1);
1321  nhist1->SetFillStyle(3002);
1322  TH1* nhist2 = new TH1F("nhist2", "NHist", 5, 0.5, 5.5);
1323  nhist2->SetFillColor(kYellow+1);
1324  nhist2->SetMarkerColor(kYellow+1);
1325  nhist2->SetLineColor(kYellow+1);
1326  nhist2->SetFillStyle(3002);
1327 
1328  TH2* resp = new TH2F("resp", "Reponse", 100, 0, 10, 6, -.5, 5.5);
1329  resp->SetXTitle("Signal");
1330  resp->SetYTitle("# particles");
1331  for (Int_t j = 0; j < fNObs; j++) {
1332  if (((j+1) % (fNObs / 10)) == 0)
1333  std::cout << "Event # " << j+1 << std::endl;
1334  Double_t x = fSummed->GetRandom();
1335  Double_t y = fit->RandomEstimateNParticles(x);
1336  nhist1->Fill(y);
1337  resp->Fill(x, y);
1338  }
1339  TGraphErrors* graph = new TGraphErrors(fSummed->GetNbinsX());
1340  graph->SetName("evalWeighted");
1341  graph->SetTitle("Evaluated function");
1342  graph->SetLineColor(kBlack);
1343  graph->SetFillColor(kYellow+1);
1344  graph->SetFillStyle(3001);
1345  for (Int_t j = 1; j <= fSummed->GetNbinsX(); j++) {
1346  Double_t x = fSummed->GetBinCenter(j);
1347  Double_t e = 0;
1348  Double_t y = fit->EstimateNParticles(x, e);
1349  nhist2->Fill(y, fSummed->GetBinContent(j));
1350  graph->SetPoint(j-1, x, y);
1351  graph->SetPointError(j-1, 0, e);
1352  }
1353  nhist1->Scale(1. / nhist1->GetMaximum());
1354  nhist2->Scale(1. / nhist2->GetMaximum());
1355  fN->Scale(1. / fN->GetMaximum());
1356 
1357  p = c->cd(3);
1358  p->SetLogy();
1359  p->SetBorderSize(0);
1360  p->SetBorderMode(0);
1361  p->SetTopMargin(0.01);
1362  p->SetRightMargin(0.01);
1363  fN->Draw();
1364  nhist1->Draw("same");
1365  nhist2->Draw("same");
1366 
1367  TLegend* l3 = new TLegend(.3, .7, .95, .95);
1368  l3->SetBorderSize(0);
1369  l3->SetFillStyle(0);
1370  l3->SetFillColor(0);
1371  l3->SetTextFont(132);
1372  l3->AddEntry(fN,
1373  Form("Input n distribution: #bar{m}=%5.3f, s_{m}=%5.3f",
1374  fN->GetMean(), fN->GetRMS()), "f");
1375  l3->AddEntry(nhist1,
1376  Form("Random draw from fit: #bar{m}=%5.3f, s_{m}=%5.3f",
1377  nhist1->GetMean(), nhist1->GetRMS()), "f");
1378  l3->AddEntry(nhist2,
1379  Form("Weighted evaluation of fit: #bar{m}=%5.3f, s_{m}=%5.3f",
1380  nhist2->GetMean(), nhist2->GetRMS()), "f");
1381  l3->Draw();
1382 
1383  // --- Fourth pad: Reponse function ------------------------------
1384  p = c->cd(4);
1385  // p->SetLogy();
1386  p->SetBorderSize(0);
1387  p->SetBorderMode(0);
1388  p->SetTopMargin(0.01);
1389  p->SetRightMargin(0.01);
1390  p->SetGridx();
1391  TProfile* prof1 = fNSum->ProfileY();
1392  prof1->SetLineColor(kMagenta+1);
1393  prof1->Draw();
1394  TProfile* prof2 = resp->ProfileX();
1395  prof2->SetLineColor(kCyan+2);
1396  prof2->Draw("same");
1397  graph->SetLineWidth(2);
1398  graph->Draw("same LP E4");
1399 
1400  TLegend* l4 = new TLegend(.2, .11, .8, .4);
1401  l4->SetBorderSize(0);
1402  l4->SetFillStyle(0);
1403  l4->SetFillColor(0);
1404  l4->SetTextFont(132);
1405  l4->AddEntry(prof1, "Input distribution of N particles", "l");
1406  l4->AddEntry(prof2, "Random estimate of N particles", "l");
1407  l4->AddEntry(graph, "E: #sum_{i}^{N}i a_{i}f_{i}/"
1408  "#sum_{i}^{N}a_{i}f_{i} evaluated over NLG", "lf");
1409  l4->Draw();
1410 
1411  c->cd();
1412 
1413  if (type && type[0] != '\0')
1414  c->Print(Form("TestELossDist.%s", type));
1415  }
1416 };
1417 
1425 void
1426 TestELossDist(const char* type="png")
1427 {
1428  gRandom->SetSeed(12345);
1429  const Double_t a1[] = { 1, .05, .02, .005, .002 }; // 7TeV pp
1430  const Double_t a2[] = { 1, .1, .05, .01, .005 };
1431  const Double_t a3[] = { 1, .2, .10, .05, .01 };
1432 
1433  Generator* g = new Generator(0.55, 0.06, 0.06, 5, a1);
1434  g->Generate(100000);
1435  g->Draw(type);
1436 
1437 
1438 }
1439 
1440 
Double_t fXi
Double_t Evaluate(Double_t x) const
double Double_t
Definition: External.C:58
void PrintParameter(Double_t input, TF1 *f, Int_t i)
Int_t GenerateN()
static Double_t landauGausI(Double_t *xp, Double_t *pp)
static Double_t ILandauGaus(Int_t i, Double_t x, Double_t delta, Double_t xi, Double_t sigma)
Definition: External.C:236
Double_t fDelta
static const Double_t fgkConvNSteps
Number of steps in integration.
Definition: TestELossDist.C:67
void PrintResult(TF1 *f)
Double_t Generate1Signal()
Double_t GetEXi() const
Constant.
Definition: TestELossDist.C:48
TH1 * fSingle
static Double_t landauGausN(Double_t *xp, Double_t *pp)
void Draw(const char *type="png")
static const Double_t fgkInvRoot2Pi
1 / sqrt(2 * pi)
Definition: TestELossDist.C:63
Double_t fSigma
static Double_t Landau(Double_t x, Double_t delta, Double_t xi)
Definition: TestELossDist.C:84
TCanvas * c
Definition: TestFitELoss.C:172
Double_t GetESigma() const
TF1 * FitFirst(TH1 *dist)
Double_t GetC() const
AliStack * stack
TRandom * gRandom
TH1 * fSimple
static Double_t IdLandauGausdPar(Int_t i, Double_t x, UShort_t ipar, Double_t dPar, Double_t delta, Double_t xi, Double_t sigma)
UShort_t GetN() const
static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi, Double_t sigma, Int_t n, Double_t *a)
Number of particles.
Definition: TestELossDist.C:56
Double_t * sigma
Double_t EstimateNParticles(Double_t x, Short_t maxN=-1)
Double_t Evaluate(Double_t x, Double_t &e) const
Sigma of gaussian.
Definition: TestELossDist.C:54
int Int_t
Definition: External.C:63
Double_t GetEC() const
void DrawParameter(Double_t x, Double_t y, Double_t input, TF1 *f, Int_t i)
Double_t GetEDelta() const
Function * FitN(TH1 *dist, Int_t n)
Double_t GetXi() const
TArrayD fP
Definition: External.C:228
Definition: External.C:212
Generator(Double_t delta=.55, Double_t xi=0.06, Double_t sigma=0.02, Int_t maxN=5, const Double_t *a=0)
TList * fitter
Definition: DrawAnaELoss.C:26
Double_t RandomEstimateNParticles(Double_t x, Short_t maxN=-1)
void Clear()
TH1 * fSummed
Double_t * GetAs() const
Int_t fNMinus
Double_t GetDelta() const
short Short_t
Definition: External.C:23
Double_t GetSigma() const
Double_t GetA(Int_t i) const
Double_t NEstimate(Double_t x, Double_t delta, Double_t xi, Double_t sigma, Int_t n, Double_t *a)
Double_t EstimateNParticles(Double_t x, Double_t &e, Short_t maxN=-1)
Width of landau.
Definition: TestELossDist.C:52
TArrayD fRP
Double_t GetEA(Int_t i) const
Definition: External.C:220
static Double_t LandauGaus(Double_t x, Double_t delta, Double_t xi, Double_t sigma)
void Generate(Int_t nObs=1000, Bool_t reset=true)
UShort_t GetEN() const
unsigned short UShort_t
Definition: External.C:28
TObjArray fFunctions
void DrawResult(Double_t x, Double_t y, TF1 *f)
const char Option_t
Definition: External.C:48
Double_t GenerateNSignal(Int_t n)
bool Bool_t
Definition: External.C:53
static const Double_t fgkMPShift
MP shift.
Definition: TestELossDist.C:61
Double_t * GetEAs() const
MPV of landau.
Definition: TestELossDist.C:50
static Double_t landauGaus1(Double_t *xp, Double_t *pp)
Fitter(Int_t nMinus)
static const Double_t fgkConvNSigma
Number of sigmas to integrate over.
Definition: TestELossDist.C:65
Definition: External.C:196
void Draw(Option_t *option="")
Weights.
Definition: TestELossDist.C:58
void TestELossDist(const char *type="png")