AliPhysics  6cf2591 (6cf2591)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Reference.C
Go to the documentation of this file.
1 #include <TH1.h>
2 #include <TF1.h>
3 #include <TMath.h>
4 #include <TGraph.h>
5 #include <TStyle.h>
6 #include <TLatex.h>
7 
8 struct Reference
9 {
10  static Double_t function(Double_t* x, Double_t* par)
11  {
12  const Double_t invSq2Pi = 1. / TMath::Sqrt(TMath::TwoPi());
13  const Double_t mpShift = -0.22278298;
14  const Int_t np = 100;
15  const Double_t sc = 5;
16 
17  Double_t sum = 0;
18  Double_t x0 = x[0];
19  Double_t xi = par[0];
20  Double_t mpc = par[1] - mpShift * xi;
21  Double_t area = par[2];
22  Double_t sigma = par[3];
23  Double_t xLow = x0 - sc * sigma;
24  Double_t xUpp = x0 + sc * sigma;
25  Double_t step = (xUpp - xLow) / np;
26 
27  for (Int_t i = 1; i <= np/2; i++) {
28  Double_t x1 = xLow + (i - .5) * step;
29  Double_t x2 = xUpp - (i - .5) * step;
30 
31  sum += TMath::Landau(x1, mpc, xi, true) * TMath::Gaus(x0,x1,sigma);
32  sum += TMath::Landau(x2, mpc, xi, true) * TMath::Gaus(x0,x2,sigma);
33  }
34 
35  return area * step * sum * invSq2Pi / sigma;
36  }
37 
38  static TF1* fit(TH1* hist,
39  Double_t* range,
40  Double_t* guess,
41  Double_t* lowLimits,
42  Double_t* uppLimits,
43  Double_t* params,
44  Double_t* errors,
45  Double_t& chi2,
46  Int_t& ndf)
47  {
48  TF1* func = new TF1(Form("func%s", hist->GetName()), &function,
49  range[0], range[1], 4);
50  func->SetParameters(guess);
51  func->SetParNames("#xi", "#Delta_{p}", "C", "#sigma");
52  for (Int_t i = 0; i < 4; i++)
53  func->SetParLimits(i, lowLimits[i], uppLimits[i]);
54 
55  hist->Fit(func, "RB0");
56 
57  for (Int_t i = 0; i < 4; i++) {
58  params[i] = func->GetParameter(i);
59  errors[i] = func->GetParError(i);
60  }
61  chi2 = func->GetChisquare();
62  ndf = func->GetNDF();
63 
64  return func;
65  }
66  static Bool_t search(TF1* f, Double_t start, Double_t step,
67  Bool_t peak, Double_t& res)
68  {
69  const Int_t maxCalls = 10000;
70  Double_t lOld = -2;
71  Double_t l = peak ? -1 : -1e300;
72  Int_t i = 0;
73  Double_t x = 0;
74  Double_t cur = start;
75 
76  while ((i < maxCalls) &&
77  TMath::Abs(l - lOld) > 1e-6 &&
78  TMath::Abs(step) > 1e-8) {
79  i++;
80  lOld = l;
81  x = cur + step;
82  l = f->Eval(x);
83  if (!peak) l = TMath::Abs(l-res);
84  // Printf("mode=%d x=%g step=%g l=%g lOld=%g", mode, x, step, l, lOld);
85 
86  if ((peak && l < lOld) || (!peak && l > lOld))
87  step = -step/10; // Go the other way in smaller steps
88 
89  cur += step;
90  }
91  if (i >= maxCalls) return false;
92  res = x;
93  return true;
94  }
95  static Int_t find(TF1* f, Int_t iXi, Int_t iMpc,
96  Double_t& maxX, Double_t& fwhm)
97  {
98 
99  Double_t xi = f->GetParameter(iXi); // params[0];
100  Double_t mpc = f->GetParameter(iMpc); // params[1];
101  if (!search(f, mpc-0.1*xi, 0.05 * xi, true, maxX)) return false;
102  Double_t left = f->Eval(maxX) / 2;
103  Double_t right = left;
104  if (!search(f, maxX-.5*xi, xi, false, left)) return false;
105  if (!search(f, maxX+xi, -xi, false, right)) return false;
106  fwhm = right - left;
107 
108  return true;
109  }
110  static TH1* histo()
111  {
112  Int_t data[100] = {0, 0, 0, 0, 0, 0, 2, 6, 11, 18,
113  18, 55, 90,141,255,323,454,563,681,737,
114  821,796,832,720,637,558,519,460,357,291,
115  279,241,212,153,164,139,106, 95, 91, 76,
116  80, 80, 59, 58, 51, 30, 49, 23, 35, 28,
117  23, 22, 27, 27, 24, 20, 16, 17, 14, 20,
118  12, 12, 13, 10, 17, 7, 6, 12, 6, 12,
119  4, 9, 9, 10, 3, 4, 5, 2, 4, 1,
120  5, 5, 1, 7, 1, 6, 3, 3, 3, 4,
121  5, 4, 4, 2, 2, 7, 2, 4};
122  TH1D *ret = new TH1D("snr","Signal-to-noise",100,0,100);
123  for (Int_t i = 0; i < 100; i++) ret->Fill(i, data[i]);
124 
125  ret->SetFillColor(kRed+1);
126  ret->SetFillStyle(3001);
127  ret->SetXTitle("#Delta");
128  ret->SetDirectory(0);
129 
130  return ret;
131  }
132  static void maxFwhm(TF1* f, Int_t iXi, Int_t iMpc)
133  {
134  printf("Find max and FWHM ...");
135  Double_t maxX = 0;
136  Double_t fwhm = 0;
137  if (!find(f, iXi, iMpc, maxX, fwhm))
138  Printf(" failed");
139  else {
140  Printf(" Max: %f FWHM: %f", maxX, fwhm);
141  Double_t y = f->Eval(maxX);
142  TGraph* maxG = new TGraph(1);
143  maxG->SetPoint(0, maxX, y);
144  maxG->SetMarkerStyle(20);
145  maxG->SetMarkerColor(kBlue+2);
146  maxG->Draw("p");
147 
148  TLatex* maxT = new TLatex(1.2*maxX, y, Form("Max: %f", maxX));
149  maxT->SetTextAlign(11);
150  maxT->Draw();
151 
152  TGraph* fwhmG = new TGraph(2);
153  fwhmG->SetPoint(0, maxX-fwhm/2, y/2);
154  fwhmG->SetPoint(1, maxX+fwhm/2, y/2);
155  fwhmG->SetMarkerStyle(21);
156  fwhmG->SetMarkerColor(kMagenta+2);
157  fwhmG->Draw("pl");
158 
159  TLatex* fwhmT = new TLatex(1.2*(maxX+fwhm/2), y/2,
160  Form("FWHM: %f", fwhm));
161  fwhmT->SetTextAlign(11);
162  fwhmT->Draw();
163  }
164 
165  }
166  static void test()
167  {
168  TH1* hist = histo();
169 
170  Double_t range[] = { .3 * hist->GetMean(), 3 * hist->GetMean() };
171  Double_t guess[] = { 1.8, 20., 5e4, 3.0 };
172  Double_t low[] = { 0.5, 5., 1., 0.4 };
173  Double_t upp[] = { 5.0, 50., 1e6, 5.0 };
174  Double_t chi2 = 0;
175  Int_t ndf = 0;
176  Double_t params[4];
177  Double_t errors[4];
178 
179  printf("Fitting ...");
180  TF1* f = fit(hist, range, guess, low, upp, params, errors, chi2, ndf);
181  Printf(" done");
182  f->Print();
183 
184  printf("Drawing ...");
185  gStyle->SetOptStat(1111);
186  gStyle->SetOptFit(111);
187  hist->Draw();
188  f->Draw("same");
189  Printf(" done");
190 
191  maxFwhm(f, 0, 1);
192  }
193 };
static Int_t find(TF1 *f, Int_t iXi, Int_t iMpc, Double_t &maxX, Double_t &fwhm)
Definition: Reference.C:95
double Double_t
Definition: External.C:58
static void test()
Definition: Reference.C:166
Double_t * sigma
static TF1 * fit(TH1 *hist, Double_t *range, Double_t *guess, Double_t *lowLimits, Double_t *uppLimits, Double_t *params, Double_t *errors, Double_t &chi2, Int_t &ndf)
Definition: Reference.C:38
int Int_t
Definition: External.C:63
Definition: External.C:212
static void maxFwhm(TF1 *f, Int_t iXi, Int_t iMpc)
Definition: Reference.C:132
static TH1 * histo()
Definition: Reference.C:110
bool Bool_t
Definition: External.C:53
Definition: External.C:196
static Bool_t search(TF1 *f, Double_t start, Double_t step, Bool_t peak, Double_t &res)
Definition: Reference.C:66