StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StdEdxModel.cxx
1 #include <assert.h>
2 #include "Riostream.h"
3 #include <stdio.h>
4 #include "TROOT.h"
5 #include "TSystem.h"
6 #include "TAxis.h"
7 #include "TH1.h"
8 #include "TH2.h"
9 #include "TSystem.h"
10 #include "TMath.h"
11 #include "StdEdxModel.h"
12 #include "Bichsel.h"
13 #include "TCallf77.h"
14 #include "St_spline3C.h"
15 #ifndef WIN32
16 #define ggderiv ggderiv_
17 #else
18 #define ggderiv GGDERIV
19 #endif
20 extern "C" {
21  void type_of_call ggderiv(Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t *);
22 };
23 //#include "StMessMgr.h"
24 using namespace std;
25 ClassImp(StdEdxModel)
26 StdEdxModel *StdEdxModel::fgStdEdxModel = 0;
27 Int_t StdEdxModel::_debug = 1;
28 //________________________________________________________________________________
29 StdEdxModel* StdEdxModel::instance() {
30  if (! fgStdEdxModel) new StdEdxModel();
31  return fgStdEdxModel;
32 }
33 //________________________________________________________________________________
34 StdEdxModel::StdEdxModel() : fScale(1)
35  , fTmaxL10eV(5) // Tcut = 100 keV
36 {
37  // LOG_INFO << "StdEdxModel:: use StTpcRSMaker model for dE/dx calculations" << endm;
38  cout << "StdEdxModel:: use StTpcRSMaker model for dE/dx calculations" << endl;
39  memset(beg, 0, end-beg+1);
40  TDirectory *dir = gDirectory;
41  fgStdEdxModel = this;
42  Int_t NFiles = 2; // old
43  if ( ! Stspline3LndNdxL10::instance()) NFiles = 1;
44  if ( Stspline3LndNdxL10::instance()->IsValid()) NFiles = 1;
45 
46  const Char_t *path = ".:./StarDb/dEdxModel:$STAR/StarDb/dEdxModel";
47 #ifndef __HEED__
48  const Char_t *Files[2] = {"dNdE_Bichsel.root","dNdx_Bichsel.root"};
49 #else
50  const Char_t *Files[2] = {"dNdx_Heed.root","dNdx_Heed.root"};
51 #endif
52 #define __Warn_Hist__(__HIST__) {m ## __HIST__ = (TH1D *) pFile->Get(#__HIST__); \
53  if (m ## __HIST__) {Warning("StdEdxModel","Histogram %s/%s has been found",m ## __HIST__->GetName(),m ## __HIST__->GetTitle()); m ## __HIST__->SetDirectory(0);}}
54  for (Int_t i = 0; i < NFiles; i++) { // files
55  Char_t *file = gSystem->Which(path,Files[i],kReadPermission);
56  if (! file) Fatal("StdEdxModel","File %s has not been found in path %s",Files[i],path);
57  else Warning("StdEdxModel","File %s has been found as %s",Files[i],file);
58  TFile *pFile = new TFile(file);
59  if (i == 1) {
60  __Warn_Hist__(dNdx);
61  __Warn_Hist__(dNdxL10);
62  __Warn_Hist__(LndNdxL10);
63  __Warn_Hist__(LndNdxL10Smooth);
64  mLndNdxL10Spline5 = (TSpline5 *) pFile->Get("LndNdxL10Spline5"); if (mLndNdxL10Spline5) {Warning("StdEdxModel","TSpline5 %s has been found",mLndNdxL10Smooth->GetName());}
65  assert(mdNdx || mdNdxL10);
66  } else {
67  mdNdEL10 = (TH1D *) pFile->Get("dNdEL10"); assert(mdNdEL10); Warning("StdEdxModel","Histogram %s/%s has been found",mdNdEL10->GetName(),mdNdEL10->GetTitle()); mdNdEL10->SetDirectory(0);
68  }
69  delete pFile;
70  delete [] file;
71  }
72 #undef __Warn_Hist__
73  dir->cd();
74  fGGaus = new TF1("GGaus",ggaus, -1., 5., 5);
75  fGGaus->SetParNames("NormL","mu","sigma","alpha","k");
76  fGGaus->SetParameters(0,0,0.3,0.1,0);
77  fGGaus->FixParameter(4,0.0);
78 
79  fGausExp = new TF1("GausExp",gausexp, -5., 5., 5);
80  fGausExp->SetParNames("NormL","mu","sigma","k","l");
81  fGausExp->SetParameters(0,0,0.3,5.,0);
82  fGausExp->FixParameter(4,0.0);
83  InitPar();
84  // Set normalization point the same as for I70 (increase energy per conduction electron from 20 eB to 52 eV)
85  // Double_t dEdxMIPLog = TMath::Log(2.62463815285237434) + 0.0174824 + 3.57110e-03; ; //TMath::Log(2.39761562607903311); // [keV/cm] for dX = 2 cm
86  Double_t dEdxMIPLog = TMath::Log(2.62463815285237434); ; //TMath::Log(2.39761562607903311); // [keV/cm] for dX = 2 cm
87  Double_t MIPBetaGamma10 = TMath::Log10(4.);
88  // log2dx, charge mass
89  Double_t pars[3] = { 1.0, 1.0, 0.13956995};
90  Double_t dEdxLog = zMP(&MIPBetaGamma10, pars);
91  fLogkeVperElectron = dEdxMIPLog - dEdxLog;
92  cout << "StdEdxModel:: set scale = " << Form("%5.1f",1e3*keVperElectron()) << " eV/electron" << endl;
93 }
94 //________________________________________________________________________________
95 StdEdxModel::~StdEdxModel() {
96  fgStdEdxModel = 0;
97  SafeDelete(mdNdx);
98  SafeDelete(mdNdxL10);
99  SafeDelete(mLndNdxL10);
100  SafeDelete(mLndNdxL10Smooth);
101  SafeDelete(mLndNdxL10Spline5);
102  SafeDelete(fGGaus);
103  SafeDelete(fGausExp);
104  SafeDelete(fpol2F);
105  SafeDelete(fpol5F);
106  SafeDelete(fpol6F);
107  SafeDelete(fpol7F);
108 }
109 //________________________________________________________________________________
110 Double_t StdEdxModel::dNdx(Double_t poverm, Double_t charge) {
111  if (!fgStdEdxModel) instance();
112  fTmaxL10eV = tmaxL10eV(poverm);
113  Double_t Q_eff = TMath::Abs(charge);
114  if (Q_eff > 1) {
115  Double_t beta = poverm/TMath::Sqrt(1.0 + poverm*poverm);
116  // Effective charge from GEANT gthion.F
117  Double_t w1 = 1.034 - 0.1777*TMath::Exp(-0.08114*Q_eff);
118  Double_t w2 = beta*TMath::Power(Q_eff,-2./3.);
119  Double_t w3 = 121.4139*w2 + 0.0378*TMath::Sin(190.7165*w2);
120  Q_eff *= 1. -w1*TMath::Exp(-w3);
121  }
122  Double_t dNdx = 0;
123  if ( Stspline3LndNdxL10::instance() && Stspline3LndNdxL10::instance()->IsValid()) {
124  dNdx = TMath::Exp(Stspline3LndNdxL10::instance()->Func()->Eval(TMath::Log10(poverm)));
125  } else {// old
126  if (mLndNdxL10Spline5) {
127  dNdx = TMath::Exp(mLndNdxL10Spline5->Eval(TMath::Log10(poverm)));
128  } else if (mLndNdxL10Smooth) {
129  dNdx = TMath::Exp(mLndNdxL10Smooth->Interpolate(TMath::Log10(poverm)));
130  } else if (mLndNdxL10) {
131  dNdx = TMath::Exp(mLndNdxL10->Interpolate(TMath::Log10(poverm)));
132  } else if (mdNdxL10) {
133  dNdx = mdNdxL10->Interpolate(TMath::Log10(poverm));
134  } else if (mdNdx) {
135  dNdx = mdNdx->Interpolate(poverm);
136  }
137  }
138  return fScale*Q_eff*Q_eff*dNdx;
139 }
140 //________________________________________________________________________________
141 Double_t StdEdxModel::dNdxL10func(Double_t *x, Double_t *p) {
142  Double_t bgL10 = x[0];
143  Double_t bg = TMath::Power(10., bgL10);
144  return instance()->dNdx(bg, 1.);
145 }
146 //________________________________________________________________________________
147 TF1 *StdEdxModel::dNdxL10F() {
148  TF1 *f = new TF1("dNdxL10F",dNdxL10func,-2,5,0);
149  return f;
150 }
151 //________________________________________________________________________________
152 Double_t StdEdxModel::dNdE() {
153  static Double_t cLog10 = TMath::Log(10.);
154  Double_t dE = TMath::Exp(cLog10*mdNdEL10->GetRandom());
155  return dE;
156 }
157 //________________________________________________________________________________
158 Double_t StdEdxModel::gausw(Double_t *x, Double_t *p) {
159  // Skew normal distribution https://en.wikipedia.org/wiki/Skew_normal_distribution
160  Int_t k = p[4]; // switch between value and derivatives
161  Double_t X = x[0];
162  Double_t NormL = p[0];
163  Double_t ksi = p[1];
164  Double_t w = p[2];
165  Double_t alpha = p[3];
166  Double_t t = (X - ksi)/w;
167  Double_t v = t/TMath::Sqrt2();
168  Double_t G = TMath::Exp(NormL)*TMath::Gaus(t,0,1,kTRUE);
169  Double_t GA = TMath::Gaus(alpha*t,0,1,kTRUE);
170  Double_t E = (1. + TMath::Erf(alpha*v));
171  Double_t V = G/w*E;
172  if (k == 0) return V;
173  Double_t dVdNormL = V;
174  if (k == 1) return dVdNormL;
175  /*
176  dV/V = dG/G - dw/w + dE/E
177  dG/G = d(-t**2/2) = -t * dt;
178  dt = d((X - ksi)/w) = -dksi/w -dw *t/w
179  GA = 1/sqrt(2*pi) exp(-(alpha*t)**2/)
180  dE/E = GA/E*((alpha*t)*(dalpha*t + alpha*dt) = GA/E*((alpha*t)*(dalpha*t + alpha*(-dksi/w -dw *t/w)))
181  dV/V = -t * dt -dw/w + GA/E*((alpha*t)*(dalpha*t + alpha*(-dksi/w -dw *t/w)))
182 
183  dlog(V)/
184 
185  */
186  /*
187  V = G / w * E
188  dt = - 1/w * dksi - t/w *dw = - (dksi + dw)/w
189  dG = G * dt = - G * (dksi + dw)/w = - V * (dksi + dw)/E
190  dv = dt/sqrt(2) = - (dksi + dw)/w/sqrt(2)
191  GA = Gaus(alpha*v)
192  dE = GA*(dalpha * v + alpha*dv) = GA * ( dalpha *v - alpha * (dksi + dw)/w/sqrt(2))
193  dV = dG / w * E - G / w * E * dw/ w + G / w * dE
194  = - G * (dksi + dw)/w / w * E - V * dw / w + G / w * GA * ( dalpha *v - alpha * (dksi + dw)/w/sqrt(2))
195  = V ( - (dksi + dw)/w - dw / w + GA/E *( dalpha *v - alpha * (dksi + dw)/w/sqrt(2)))
196  = V ( - (dksi + 2* dw)/w + GA/E *( dalpha *v - alpha * (dksi + dw)/w/sqrt(2)))
197 
198  dV/dksi = V ( - dksi /w - GA/E *alpha * dksi /w/sqrt(2))
199  = V ( - 1 /w - GA/E *alpha /w/sqrt(2))
200  = - V/w*(1 + GA/E*alpha/sqrt(2))
201  dV/dw = V ( - 2* dw /w - GA/E alpha * + dw /w/sqrt(2))
202  = - V/w *(2 + GA/E*alpha/sqrt(2))
203  dV/dalpha = GA*v
204 
205  */
206 /* maxima
207 t(ksi,w) := (x - ksi)/w;
208 v(ksi,w) := t(ksi,w)/sqrt(2);
209 G(ksi,w) := gaussprob(t(ksi,w))/w;
210 E(ksi,w,alpha) := 1 + erf(alpha*v(ksi,w));
211 Val(ksi,w,alpha):= gaussprob(t(ksi,w))/w * ( 1 + erf(alpha*v(ksi,w)));
212 G : jacobian([Val(ksi,w,alpha)], [ksi,w,alpha]);
213 trigsimp(%);
214 fortran(%);
215 
216 */
217 //Double_t dVdksi = - V/w*(1 + GA/E*alpha/TMath::Sqrt2());
218  Double_t dVdksi = - (V + G/w*GA*alpha/TMath::Sqrt2())/w;
219  if (k == 2) return dVdksi;
220 //Double_t dVdw = - V/w*(2 + GA/E*alpha/TMath::Sqrt2());
221  Double_t dVdw = - (2*V + G/w*GA*alpha/TMath::Sqrt2())/w;
222  if (k == 3) return dVdw;
223  Double_t dVdalpha = GA*v;
224  return dVdalpha;
225 }
226 //_______________________________________________________________________________
227 Double_t StdEdxModel::ggaus(Double_t *x, Double_t *p) {
228  return ggausD(x,p,0);
229 }
230 //_______________________________________________________________________________
231  Double_t StdEdxModel::ggausD(Double_t *x, Double_t *p, Double_t *der) {
232  // Int_t k = p[4];
233  Double_t NormL = p[0];
234  Double_t mu = p[1]; // Mode = ksi + w *m_0(alpha); most propable value
235  Double_t sigma = p[2]; // Sqrt(Variance);
236  Double_t alpha = p[3];
237  Double_t ksi = mu;
238  Double_t w = sigma;
239  if (TMath::Abs(alpha) > 1e-7) {
240  Double_t delta =alpha/TMath::Sqrt(1 + alpha*alpha);
241  Double_t muz = delta/TMath::Sqrt(TMath::PiOver2());
242  Double_t sigmaz = TMath::Sqrt(1 - muz*muz);
243  Double_t gamma1 = (4 - TMath::Pi())/2 * TMath::Power(delta*TMath::Sqrt(2./TMath::Pi()), 3)
244  / TMath::Power(1 - 2*delta*delta/TMath::Pi(), 1.5);
245  Double_t m_0 = muz - gamma1*sigmaz/2 - TMath::Sign(1.,alpha)/2*TMath::Exp(-2*TMath::Pi()/TMath::Abs(alpha));
246  w = sigma/TMath::Sqrt(1 - 2* delta*delta/TMath::Pi());
247  ksi = mu - w*m_0;
248  // Double_t mean = ksi + w * muz;
249  }
250  Double_t par[5] = {NormL, ksi, w, alpha, 0};
251  Double_t V = gausw(x, par);
252  if (der) {
253 #if 0 /* Derivatives */
254  /* Maxima
255  load (f90) $
256  :lisp (setq *f90-output-line-length-max* 1000000000)
257  stardisp: true$
258  delta(alpha) := alpha/sqrt(1 + alpha*alpha);
259  muz(alpha) := delta(alpha)/sqrt(%pi/2);
260  sigmaz(alpha) := sqrt(1 - muz(alpha)*muz(alpha));
261  gamma1(alpha) := (4 - %pi)/2 * (delta(alpha)/sqrt(%pi/2))**3 / (1 - 2*delta(alpha)**2 /%pi)**1.5;
262  signn(x) := x/abs(x);
263  m_0(alpha) := muz(alpha) - gamma1(alpha)/sigmaz(alpha)/2 - signn(alpha)/2*exp(-2*%pi/abs(alpha));
264  w(sigma,alpha):= sigma/sqrt(1 - 2*delta(alpha)**2/%pi);
265  ksi(mu,sigma,alpha):= mu - w(sigma,alpha)*m_0(alpha);
266 
267  F : jacobian([ ksi(mu,sigma,alpha), w(sigma,alpha), alpha], [mu, sigma, alpha]);
268  trigsimp(%);
269  f90(%);%(1,1) = 1
270 %(1,2) = -(sqrt(%pi*alpha**2+%pi)*exp(-(2*%pi)/abs(alpha))*(sqrt(%pi)*sqrt(alpha**2+1)*((6*%pi**2-24*%pi+16)*alpha**5+(10*%pi**2-24*%pi)*alpha**3+4*%pi**2*alpha)*exp((2*%pi)/abs(alpha))*abs(alpha)+((-sqrt(2)*%pi**3)+2**(5.0d+0/2.0d+0)*%pi**2-2**(5.0d+0/2.0d+0)*%pi)*alpha**7+((-3*sqrt(2)*%pi**3)+2**(7.0d+0/2.0d+0)*%pi**2-2**(5.0d+0/2.0d+0)*%pi)*alpha**5+(2**(5.0d+0/2.0d+0)*%pi**2-3*sqrt(2)*%pi**3)*alpha**3-sqrt(2)*%pi**3*alpha))/(sqrt((%pi-2)*alpha**2+%pi)*((2**(3.0d+0/2.0d+0)*%pi**3-2**(7.0d+0/2.0d+0)*%pi**2+2**(7.0d+0/2.0d+0)*%pi)*alpha**6+(3*2**(3.0d+0/2.0d+0)*%pi**3-2**(9.0d+0/2.0d+0)*%pi**2+2**(7.0d+0/2.0d+0)*%pi)*alpha**4+(3*2**(3.0d+0/2.0d+0)*%pi**3-2**(7.0d+0/2.0d+0)*%pi**2)*alpha**2+2**(3.0d+0/2.0d+0)*%pi**3)*abs(alpha))
271 %(1,3) = (sqrt((%pi-2)*alpha**2+%pi)*sqrt(%pi*alpha**2+%pi)*exp(-(2*%pi)/abs(alpha))*(sqrt(%pi)*sqrt(alpha**2+1)*(((2*%pi**4-12*%pi**3+24*%pi**2-16*%pi)*alpha**8+(8*%pi**4-36*%pi**3+48*%pi**2-16*%pi)*alpha**6+(12*%pi**4-36*%pi**3+24*%pi**2)*alpha**4+(8*%pi**4-12*%pi**3)*alpha**2+2*%pi**4)*abs(alpha)+(2*%pi**2-8*%pi+8)*alpha**8+(4*%pi**2-8*%pi)*alpha**6+2*%pi**2*alpha**4)*sigma+(((-5*sqrt(2)*%pi**3)+2**(9.0d+0/2.0d+0)*%pi**2+2**(7.0d+0/2.0d+0)*%pi)*alpha**8+((-3*2**(5.0d+0/2.0d+0)*%pi**3)+9*2**(5.0d+0/2.0d+0)*%pi**2+2**(7.0d+0/2.0d+0)*%pi)*alpha**6+(5*2**(5.0d+0/2.0d+0)*%pi**2-9*sqrt(2)*%pi**3)*alpha**4-2**(3.0d+0/2.0d+0)*%pi**3*alpha**2)*exp((2*%pi)/abs(alpha))*abs(alpha)*sigma))/(sqrt(%pi)*sqrt(alpha**2+1)*((2*%pi**4-16*%pi**3+48*%pi**2-64*%pi+32)*alpha**12+(10*%pi**4-64*%pi**3+144*%pi**2-128*%pi+32)*alpha**10+(20*%pi**4-96*%pi**3+144*%pi**2-64*%pi)*alpha**8+(20*%pi**4-64*%pi**3+48*%pi**2)*alpha**6+(10*%pi**4-16*%pi**3)*alpha**4+2*%pi**4*alpha**2)*abs(alpha))
272 %(2,1) = 0
273 %(2,2) = sqrt(%pi*alpha**2+%pi)/sqrt((%pi-2)*alpha**2+%pi)
274 %(2,3) = (2*alpha*sqrt(%pi*alpha**2+%pi)*sigma)/(sqrt((%pi-2)*alpha**2+%pi)*((%pi-2)*alpha**4+(2*%pi-2)*alpha**2+%pi))
275 %(3,1) = 0
276 %(3,2) = 0
277 %(3,3) = 1
278 
279 
280  t(ksi,w) := (x - ksi)/w;
281  v(ksi,w) := t(ksi,w)/sqrt(2);
282  G(ksi,w) := 1./sqrt(2*%pi)*exp(-t(ksi,w)**2/2);
283  E(ksi,w,alpha) := 1 + erf(alpha*v(ksi,w));
284  Val(ksi,w,alpha):= G(ksi,w)/w *E(ksi,w,alpha);
285 
286  G : jacobian([Val(ksi,w,alpha)], [ksi,w,alpha]);
287  trigsimp(%);
288  f90(%);
289 
290  T: G . F;
291  trigsimp(%);
292  A : %;
293  f90(A);
294  with_stdout ("A.txt", f90(A));
295  */
296 #endif
297  der[0] = der[1] = der[2] = 0;
298  ggderiv(x[0], ksi, sigma, w, alpha, der);
299  }
300  return V;
301 }
302 //_______________________________________________________________________________
303 Double_t StdEdxModel::gausexp(Double_t *x, Double_t *p) {
304  return gausexpD(x,p);
305 }
306 //_______________________________________________________________________________
307 Double_t StdEdxModel::gausexpD(Double_t *x, Double_t *p, Double_t *der) {
308  // Souvik Das, "A simple alternative to the Crystal Ball function"
309  // https://arxiv.org/pdf/1603.08591.pdf
310  // Int_t l = p[4];//
311  Double_t normL = p[0];
312  Double_t mu = p[1];
313  Double_t sigma = p[2];
314  Double_t kh = p[3];
315  Double_t t = (x[0] - mu)/sigma;
316  Double_t V = 0.;
317  Double_t k = kh;
318  if (kh < 0) {
319  k = - kh;
320  t = - t;
321  }
322  Double_t dVdt = 0;
323  Double_t dVdk = 0;
324  if (t < k) {
325  V = TMath::Exp(-t*t/2); // dV/dt = - V * t => - k * V
326  dVdt = -V*t;
327  } else {
328  V = TMath::Exp(k*k/2 - k*t); // dV/dt = - k * V
329  dVdt = -V*k;
330  dVdk = V*(k - t);
331  }
332 // (%i5) integrate(exp(((-t)*t)/2),t,minf,k)
333 // k
334 // sqrt(%pi) erf(-------)
335 // sqrt(2) sqrt(%pi)
336 // (%o5) ---------------------- + ---------
337 // sqrt(2) sqrt(2)
338 // _
339 
340 // (%i6) integrate(exp((k*k)/2-k*t),t,k,inf)
341 // Is k positive, negative or zero?
342 
343 // positive
344 // ;
345 // 2
346 // k
347 // - --
348 // 2
349 // %e
350 // (%o6) ------
351 // k
352  static Double_t SQ2pi = TMath::Sqrt(TMath::Pi())/TMath::Sqrt2();
353  Double_t D = SQ2pi*(1 + TMath::Erf(k/TMath::Sqrt2()));
354  Double_t C = TMath::Exp(-k*k/2)/k;
355  Double_t N = (D + C)*sigma;
356  Double_t Prob = TMath::Exp(normL)*V/N;
357  if (der) {
358  Double_t dtdM = -1./sigma;
359  Double_t dtdS = -t /sigma;
360  Double_t dDdk = TMath::Exp(-k*k/2.);
361  Double_t dCdk = - C*(k + 1./(k*k));
362  Double_t dNdk = (dDdk + dCdk) * sigma;
363  Double_t dPdM = Prob * dVdt / V * dtdM; // over mu
364  Double_t dPdS = Prob * dVdt / V * dtdS; // over sigma
365  Double_t dPdk = Prob *(dVdk / V - dNdk / N);
366  der[0] = dPdM;
367  der[1] = dPdS;
368  der[2] = dPdk;
369  }
370  return Prob;
371 }
372 //________________________________________________________________________________
373 /* CERN-77-09 Sauli Yellow report
374  Gas dens(g/cm**3) Wi (ev) dE/dx (keV/cm) np (1/cm) nT(1/cm) nT/np
375  Ar 1.66e-3 26 2.44 29.4 94 3.2
376  CH4 6.70e-4 13.1 2.21 16 53 3.3
377  P10 1.561ee-3 25.44 2.43 28.82 92.24 3.2 TpcRS = <nT/np> = 2.47
378  P10: (dE/dx)/nT = 2.43/92.24 = 26.3 eV : W = 25.4 eV TpcRS W = 24.3
379  Ar: (dE/dx)/nT = 2.44/94 = 25.96 eV : W = 26 eV
380 */
381 //________________________________________________________________________________
382 void StdEdxModel::Parameters(Double_t Np, Double_t *parameters, Double_t *dPardNp) {
383  parameters[0] = parameters[1] = parameters[2] = 0;
384  if (Np <= 1.0) return;
385  for (Int_t l = 0; l < 3; l++) {
386  if (! dPardNp) {
387  parameters[l] = Parameter(Np, l);
388  } else {
389  parameters[l] = Parameter(Np, l, &dPardNp[l]);
390  }
391  }
392 }
393 //________________________________________________________________________________
394 Double_t StdEdxModel::Parameter(Double_t Np, Int_t l, Double_t *dPardNp) {
395  // parameters from dEdxFit::FitGG4
396  static Double_t parsA[2] = { 5.4634, -0.57598}; //alpha x
397  static Double_t parsS[3] = { 1.6924, -1.2912, 0.24698}; //sigma versus log(x)
398  static Double_t parsM[8] = { -4.3432, 4.6327, -1.9522, 0.4691, -0.066615, 0.0055111, -0.00024531, 4.5394e-06}; //mu pol7
399  Double_t x = TMath::Log(Np);
400  Double_t dxdNp = 1./Np;
401  if (l == 2) {
402  Double_t alpha = parsA[0] + x * parsA[1];
403  if (dPardNp) dPardNp[0] = parsA[1] * dxdNp;
404  return alpha;
405  } else if (l == 1) {
406  Double_t xx = (x > 0) ? TMath::Log(x) : 0;
407  Double_t sigma = parsS[0] + xx * ( parsS[1] + xx * parsS[2]);
408  if (dPardNp) {
409  Double_t dxxdx = 1./xx;
410  Double_t dxxdNp = dxxdx * dxdNp;
411  dPardNp[0] = (parsS[1] + 2 * xx * parsS[2]) * dxxdNp;
412  }
413  return sigma;
414  } else if (l == 0) {
415  Double_t mu = fpol7F->EvalPar(&x, parsM);
416  if (dPardNp) {
417  static Double_t parsMD[7] = {0};
418  if (!parsMD[0]) {
419  for (Int_t i = 0; i < 7; i++) {
420  parsMD[i] = (i+1)*parsM[i+1];
421  }
422  }
423  dPardNp[0] = fpol6F->EvalPar(&x, parsMD) * dxdNp;
424  }
425  return mu;
426  } else {
427  assert(0);
428  return 0;
429  }
430 }
431 //________________________________________________________________________________
432 Double_t StdEdxModel::MukeV(Double_t Np) {
433  return Parameter(Np, 0) + fLogkeVperElectron + TMath::Log(Np);
434 }
435 //________________________________________________________________________________
436 Double_t StdEdxModel::funParam(Double_t *x, Double_t *p) {
437  Int_t l = p[0];
438  if (l < 0 || l > 2) return 0;
439  Double_t Np = TMath::Exp(x[0]);
440  return StdEdxModel::instance()->Parameter(Np, l);
441 }
442 //________________________________________________________________________________
443 TF1 *StdEdxModel::FParam(Int_t l) {
444  const Char_t *fNames[3] = {"MuPar","sigmaPar","alphaPar"};
445  TF1 *f = 0;
446  if (l < 0 || l > 2) return f;
447  f = (TF1 *) gROOT->GetListOfFunctions()->FindObject(fNames[l]);
448  if (! f) {
449  f = new TF1(fNames[l],funParam,1.5,12,1);
450  f->SetParName(0,fNames[l]);
451  f->SetParameter(0,l);
452  cout << "Create FParam with name " << f->GetName() << endl;
453  }
454  return f;
455 
456 }
457 //________________________________________________________________________________
458 Double_t StdEdxModel::Prob(Double_t /* log(nE/Np) */ ee, Double_t Np, Double_t *der) { // GG: ggaus
459  Double_t params[5] = {0};
460  Double_t V = 0;
461  if (! der) {
462  Parameters(Np, &params[1]);
463  V = ggaus(&ee, params);
464  } else {
465  Double_t dPardNp[3] = {0};
466  Parameters(Np, &params[1], dPardNp);
467  Double_t dVdP[3] = {0};
468  V = ggausD(&ee, params, dVdP);
469  der[0] = 0;
470  for (Int_t l = 0; l < 3; l++) der[0] += dVdP[l]*dPardNp[l];
471  static Int_t _debug = 0;
472  if (_debug) {
473  Double_t xP = TMath::Log(Np);
474  Double_t D = instance()->FProbP()->Derivative(xP,&ee)/Np;
475  cout << "estimated derivative (xP = " << xP << ", ee = " << ee << ") = " << der[0] << " calculated derivative = " << D << endl;
476  }
477  }
478  return V;
479 }
480 //--------------------------------------------------------------------------------
481 Double_t StdEdxModel::funcProb(Double_t *x, Double_t *p) {
482  return instance()->Prob(x[0],p[0]);
483 }
484 //________________________________________________________________________________
485 TF1 *StdEdxModel::FProb() {
486  const Char_t *name = "GGProb";
487  TF1 *f = (TF1 *) gROOT->GetListOfFunctions()->FindObject(name);
488  if (! f) {
489  f = new TF1(name,funcProb,-1,4,1);
490  f->SetParName(0,"Np");
491  f->SetParameter(0,32);
492  cout << "Create FProb with name " << f->GetName() << endl;
493  }
494  return f;
495 }
496 //--------------------------------------------------------------------------------
497 Double_t StdEdxModel::funcProbP(Double_t *x, Double_t *p) {
498  if (p[0] < 1) return 0;
499  return instance()->Prob(TMath::Log(p[0]),x[0]);
500 }
501 //________________________________________________________________________________
502 TF1 *StdEdxModel::FProbP() {
503  const Char_t *name = "GGProbP";
504  TF1 *f = (TF1 *) gROOT->GetListOfFunctions()->FindObject(name);
505  if (! f) {
506  f = new TF1(name,funcProbP,2,12,1);
507  f->SetParName(0,"x");
508  f->SetParameter(0,1);
509  cout << "Create FProbP with name " << f->GetName() << endl;
510  }
511  return f;
512 }
513 //--------------------------------------------------------------------------------
514 Double_t StdEdxModel::funcProbDer(Double_t *x, Double_t *p) {
515  Double_t der;
516  instance()->Prob(x[0], p[0], &der);
517  return der;
518 }
519 //________________________________________________________________________________
520 TF1 *StdEdxModel::FProbDer() {
521  const Char_t *name = "GGProbDer";
522  TF1 *f = (TF1 *) gROOT->GetListOfFunctions()->FindObject(name);
523  if (! f) {
524  f = new TF1(name,funcProbDer,-1,4,1);
525  f->SetParName(0,"Np");
526  f->SetParameter(0,32);
527  cout << "Create FProb with name " << f->GetName() << endl;
528  }
529  return f;
530 }
531 #if 0
532 //________________________________________________________________________________
533 Double_t StdEdxModel::ProbEx(Double_t /* log(nE/Np) */ ee, Double_t Np, Double_t *der) { // GEX : gausexp
534  Double_t params[5] = {0};
535  Double_t V = 0;
536  if (! der) {
537  Parameters(Np, &params[1]);
538  V = gausexp(&ee, params);
539  } else {
540  Double_t dPardNp[3] = {0};
541  Parameters(Np, &params[1], dPardNp);
542  Double_t dVdP[3] = {0};
543  V = gausexpD(&ee, params, dVdP);
544  der[0] = 0;
545  for (Int_t l = 0; l < 3; l++) der[0] += dVdP[l]*dPardNp[l];
546  }
547  return V;
548 }
549 #endif
550 //________________________________________________________________________________
551 Double_t StdEdxModel::ProbdEGeVlog(Double_t dEGeVLog, Double_t Np, Double_t *der) {
552  Double_t ee = Logne(dEGeVLog) - TMath::Log(Np);
553  return Prob(ee, Np, der);
554 }
555 //________________________________________________________________________________
556 Double_t StdEdxModel::zMP(Double_t *x, Double_t *p) { // log(keV/cm)
557  Double_t bgL10 = x[0];
558  Double_t pOverMRC = TMath::Power(10., bgL10);
559  // Double_t bgL10 = TMath::Log10(pOverM);
560  Double_t log2dx = p[0];
561  Double_t charge = p[1];
562  Double_t mass = p[2];
563  Double_t dx = TMath::Power( 2., log2dx);
564  Double_t dNdx = StdEdxModel::instance()->dNdxEff(pOverMRC, charge, mass); // */dNdxVsBgC*.root [-1.5,5]
565  Double_t Np = dNdx*dx;
566  // Double_t NpLog = TMath::Log(Np);
567  // Double_t mu = instance()->Parameter(Np, 0);
568  // Double_t sigma = instance()->Parameter(Np, 1);
569  // Double_t alpha = instance()->Parameter(Np, 2);
570  // Double_t dEkeVLog = NpLog + mu -3.13746587897608142e+00 +1.78334647296254700e-01;// + 7.02725079814016507e+00;// - 3.13746587897608142e+00;// 43.4 eV/conducting electron
571  Double_t dEkeVLog = instance()->MukeV(Np); // Parameter(Np, 0);
572  Double_t dEdxLog = dEkeVLog - TMath::Log(dx);
573  Double_t dEdxCor = 0;
574 #if 0
575  StElectonsDEV_dEdx *EL = StElectonsDEV_dEdx::instance();
576  if ( EL && EL->IsValid() && EL->InRange(bgL10)) {
577  dEdxCor = StElectonsDEV_dEdx::instance()->Func()->Eval(bgL10);
578  static TF1 *elCor1 = 0;
579  if (! elCor1) {
580  Double_t pars1[4] = {-0.2100929, 0.1500455, -0.02743834, 0.001894849}; //electrons [2.1,3.3]
581  Double_t pars2[4] = {-0.04352509, 0.03437069, -0.002851349, -0.0003568138}; //electornsD
582  Double_t pars[4] = {0};
583  for (Int_t i = 0; i < 4; i++) pars[i] = pars1[i] + pars2[i];
584  elCor1 = new TF1("dEdxElCor1","pol3",2.0,3.5); elCor1->SetParameters(pars);
585  }
586  dEdxCor += elCor1->Eval(bgL10);
587  dEdxCor += -7.63891e-03 - 3.57110e-03 ;
588  } else { // pions and protons
589  static Double_t pionM = 0.13956995;
590  static Double_t protonM = 0.9382723;
591  static Double_t mPionL10 = TMath::Log10(pionM);
592  static Double_t mProtonL10 = TMath::Log10(protonM);
593  static Double_t dML10 = mProtonL10 - mPionL10;
594  Double_t dEdxCorPion = 0;
595  Double_t dEdxCorProton = 0;
596  StPionDEV_dEdx *PI = StPionDEV_dEdx::instance();
597  if ( PI && PI->IsValid() && PI->InRange(bgL10)) {
598  dEdxCorPion = PI->Func()->Eval(bgL10);
599  static TF1 *piCor1 = 0;
600  if (! piCor1) {
601  Double_t pars[4] = {0.008567411, 0.01817216, -0.004932309, -0.001434306}; //pionD
602  piCor1 = new TF1("dEdxPiCor1","pol3",-0.2,1.6); piCor1->SetParameters(pars);
603  }
604  dEdxCorPion += piCor1->Eval(bgL10);
605  }
606  StProtonDEV_dEdx *P = StProtonDEV_dEdx::instance();
607  if (P &&P->IsValid() &&P->InRange(bgL10)) {
608  dEdxCorProton = P->Func()->Eval(bgL10);
609  static TF1 *protonCor1 = 0;
610  if (! protonCor1) {
611  Double_t pars[6] = {0.01745018 + 0.005654033 - 3.57110e-03 , 0.005726225 -0.00228347, 0.004416636, -0.02814983, 0.1824491, -0.2114645}; //protonD
612  // Double_t pars2[2] = {0.005654033, -0.00228347}; //proton dEdxG
613  protonCor1 = new TF1("dEdxProtonCor1","pol5",-0.5,0.8); protonCor1->SetParameters(pars);
614  }
615  dEdxCorProton += protonCor1->Eval(bgL10);
616  }
617  Double_t mL10 = TMath::Log10(mass);
618  dEdxCor += dEdxCorPion + (dEdxCorProton - dEdxCorPion)*(mL10 - mPionL10)/dML10;
619  }
620  // dEdxCor = 1.66944e-02;
621  // dEdxCor += 0.0174824;
622 #endif
623  return dEdxLog + dEdxCor;
624 }
625 //________________________________________________________________________________
626 TF1 *StdEdxModel::ZMP(Double_t log2dx, Double_t charge, Double_t mass) {
627  TString fName(Form("New%i",(int)(2*(log2dx+2))));
628  TF1 *f = (TF1 *) gROOT->GetListOfFunctions()->FindObject(fName);
629  if (! f) {
630  f = new TF1(fName,zMP,-2,5,3);
631  f->SetLineStyle(4);
632  f->SetParNames("log2dx","charge","mass");
633  f->SetParameter(0,log2dx);
634  f->SetParameter(1, charge);
635  f->SetParameter(2, mass);
636  cout << "Create ZMPNew with name " << f->GetName() << " for log2dx = " << log2dx << " and mass = " << mass << endl;
637  }
638  return f;
639 }
640 //________________________________________________________________________________
641 void StdEdxModel::InitPar() {
642  fpol2F = (TF1 *) gROOT->GetListOfFunctions()->FindObject("pol2");
643  fpol5F = (TF1 *) gROOT->GetListOfFunctions()->FindObject("pol5");
644  fpol6F = (TF1 *) gROOT->GetListOfFunctions()->FindObject("pol6");
645  fpol7F = (TF1 *) gROOT->GetListOfFunctions()->FindObject("pol7");
646  if (! fpol2F || ! fpol5F || ! fpol6F || ! fpol7F) {
647  TF1::InitStandardFunctions();
648  fpol2F = (TF1 *) gROOT->GetListOfFunctions()->FindObject("pol2");
649  fpol5F = (TF1 *) gROOT->GetListOfFunctions()->FindObject("pol5");
650  fpol6F = (TF1 *) gROOT->GetListOfFunctions()->FindObject("pol6");
651  fpol7F = (TF1 *) gROOT->GetListOfFunctions()->FindObject("pol7");
652  }
653 }
654 //________________________________________________________________________________
655 Double_t StdEdxModel::tmaxL10eV(Double_t bg) {
656  static Double_t Tcut = 1e-4; // 100 keV maximum cluster size (~80 keV)
657  static Double_t m_e = 0.51099907e-3;
658  Double_t bg2 = bg*bg;
659  // Double_t gamma = TMath::Sqrt(bg2 + 1);
660  Double_t tMax = 2*m_e*bg2; // /(1. + mOverM*(2*gamma + mOverM));
661  return TMath::Log10(1e9*TMath::Min(Tcut, tMax));
662 }
663 //________________________________________________________________________________
664 Double_t StdEdxModel::saturationTanH(Double_t *x, Double_t *p) { // nP saturation versus beta*gamma from TpcRS (nP/dX - dN/dx_model)
665  // TF1 *s8 = new TF1("s8","[0]+[1]*TMath::TanH([2]+x*([3]+x*([4] + x*([5] +x*([6] + x*[7])))))",-2,5)
666  return p[0]+p[1]*TMath::TanH(p[2]+x[0]*(p[3]+x[0]*(p[4] + x[0]*(p[5] +x[0]*(p[6] + x[0]*p[7])))));
667 }
668 //________________________________________________________________________________
669 TF1 *StdEdxModel::SaturTanH() {
670  static TF1 *f = 0;
671  if (! f) {
672  f = new TF1("SaturTanHF",StdEdxModel::saturationTanH,-5,5,8);
673  Double_t pars[8] = { -0.060944, -0.014597, 1.6066, -3.4821, 3.4131, -1.3879, 0.26295, -0.019545}; //
674  f->SetParNames("offset","slope","a0","a1","a2","a3","a4","a5","a6");
675  f->SetParameters(pars);
676  }
677  return f;
678 }
679 //________________________________________________________________________________
680 Double_t StdEdxModel::saturationFunc(Double_t *x, Double_t *p) { // nP saturation versus beta*gamma from TpcRS (nP/dX - dN/dx_model)
681  return (p[0] + p[1]*TMath::TanH(p[2]*(x[0] - p[3])))*(1 + x[0]*(p[4] + x[0]*(p[5] + x[0]*p[6])));
682 }
683 //________________________________________________________________________________
684 TF1 *StdEdxModel::Saturation(Int_t particle) {
685  static TF1 *f = 0;
686  // "offset","slope","scale","shift","decay",decay2","decay3"
687  Double_t params[1][7] = {{ 0., 1., 1., 0., 0.0, 0.0, 0.0}};
688  /*
689 // /hlt/cephfs/fisyak/Fit Sat Sep 24 11:19:27 2022
690  Double_t pars[7] = { -0.077712, -0.0033412, -3.3287, 2.8284, -0.02002, 0, 0}; // electron-COL chisq = 400.355859 / NDF = 59
691  Double_t pars[7] = { -0.038331, -0.0039637, -2.1902, 2.7107, 0.56673, -0.08022, 0}; // electron-COL chisq = 172.642030 / NDF = 58
692  Double_t pars[7] = { -0.045392, -0.006256, -1.9007, 2.6929, 0.18645, 0.062376, -0.015684}; // electron-COL chisq = 168.361785 / NDF = 57
693  Double_t pars[7] = { -0.079957, -0.0026383, -3.3327, 2.9211, -0.026865, 0, 0}; // electron+COL chisq = 490.739132 / NDF = 59
694  Double_t pars[7] = { -0.043914, -0.0035152, -2.273, 2.7878, 0.43439, -0.064819, 0}; // electron+COL chisq = 192.651659 / NDF = 58
695  Double_t pars[7] = { -0.047536, -0.0041229, -2.1841, 2.7883, 0.29741, -0.019738, -0.004555}; // electron+COL chisq = 192.011784 / NDF = 57
696  Double_t pars[7] = { -0.078978, -0.0029429, -3.3389, 2.8729, -0.024103, 0, 0}; // electronCOL chisq = 843.389038 / NDF = 59
697  Double_t pars[7] = { -0.0415, -0.0037696, -2.2214, 2.749, 0.48601, -0.070646, 0}; // electronCOL chisq = 336.681579 / NDF = 58
698  Double_t pars[7] = { -0.047724, -0.0049973, -2.0704, 2.7477, 0.23655, 0.012412, -0.0084929}; // electronCOL chisq = 332.928617 / NDF = 57
699  Double_t pars[7] = { -0.079857, -0.0020054, -5.3142, 2.801, -0.024112, 0, 0}; // electron-FXT chisq = 453.043418 / NDF = 59
700  Double_t pars[7] = { -0.069282, -0.0025422, -4.043, 2.7721, 0.062833, -0.012376, 0}; // electron-FXT chisq = 390.341916 / NDF = 58
701  Double_t pars[7] = { -0.12051, -0.02907, -1.3553, 2.6704, -0.66218, 0.28614, -0.033912}; // electron-FXT chisq = 279.392893 / NDF = 57
702  Double_t pars[7] = { -0.082144, -0.0012493, -5.9997, 2.9024, -0.031554, 0, 0}; // electron+FXT chisq = 1093.876523 / NDF = 60
703  Double_t pars[7] = { -0.060189, -0.0022009, -3.7086, 2.8155, 0.17296, -0.029677, 0}; // electron+FXT chisq = 748.738627 / NDF = 59
704  Double_t pars[7] = { -0.13164, -0.057247, -1.0776, 2.4509, -0.88342, 0.41293, -0.049502}; // electron+FXT chisq = 453.312446 / NDF = 58
705  Double_t pars[7] = { -0.081057, 0.0015963, 5.5394, 2.8465, -0.028146, 0, 0}; // electronFXT chisq = 1477.459419 / NDF = 60
706  Double_t pars[7] = { -0.064431, 0.0023793, 3.7896, 2.7925, 0.11764, -0.020966, 0}; // electronFXT chisq = 1115.434973 / NDF = 59
707  Double_t pars[7] = { -0.12629, 0.045864, 1.1627, 2.5427, -0.80495, 0.36744, -0.043926}; // electronFXT chisq = 673.906498 / NDF = 58
708 
709  */
710  if (! f) {
711  f = new TF1("SaturationF",StdEdxModel::saturationFunc,-5,5,7);
712  f->SetParNames("offset","slope","scale","shift","decay","decay2","decay3");
713  }
714  f->SetParameters(params[particle]);
715  return f;
716 }
717 //________________________________________________________________________________
718 Double_t StdEdxModel::bgCorrected(Double_t bgRC) {
719  // Parameterization of correction from /hlt/cephfs/fisyak/Fit/*/dBGLGADCut23.root 09/25/2022
720  Double_t pars[3] = {-0.00020089, 0.0031976, 0.0062467}; //
721  static TF1 *pol2 = 0;
722  if (! pol2) {
723  pol2 = (TF1 *) gROOT->GetListOfFunctions()->FindObject("pol2");
724  if (! pol2) {
725  TF1::InitStandardFunctions();
726  pol2 = (TF1 *) gROOT->GetListOfFunctions()->FindObject("pol2");
727  }
728  assert(pol2);
729  }
730  Double_t bgRC2 = bgRC*bgRC;
731  Double_t beta2RC = bgRC2/(bgRC2 + 1);
732  Double_t betaRCL10 = 0.5*TMath::Log(beta2RC);
733  Double_t bgMC = bgRC;
734  if (betaRCL10 < -0.5) bgMC += pol2->EvalPar(&betaRCL10, pars);
735  return bgMC;
736 }
737 //________________________________________________________________________________
738 TH1D *StdEdxModel::protonEff() {
739 //========= Macro generated from object: Func/
740 //========= by ROOT version5.34/39
741 
742  TH1D *eff = new TH1D("protonEff","",100,-2.0,0.0);
743  Double_t corr[102] = { 0,
744  0.8300, 0.8335, 0.8367, 0.8396, 0.8422, 0.8446, 0.8467, 0.8485, 0.8501, 0.8516,
745  0.8575, 0.8635, 0.8691, 0.8741, 0.8784, 0.8822, 0.8853, 0.8879, 0.8901, 0.8920,
746  0.8935, 0.8948, 0.8959, 0.8968, 0.8975, 0.8981, 0.8987, 0.8991, 0.8995, 0.8999,
747  0.9002, 0.9004, 0.9006, 0.9008, 0.9010, 0.9012, 0.9013, 0.9014, 0.9016, 0.9018,
748  0.9020, 0.9021, 0.9022, 0.9023, 0.9024, 0.9025, 0.9025, 0.9026, 0.9026, 0.9027,
749  0.9027, 0.9028, 0.9028, 0.9028, 0.9029, 0.9029, 0.9029, 0.9029, 0.9029, 0.9030,
750  0.9030, 0.9030, 0.9030, 0.9030, 0.9030, 0.9030, 0.9030, 0.9030, 0.9030, 0.9031,
751  0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031,
752  0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031,
753  0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031, 0.9031,
754  0};
755  eff->Set(102, corr);
756  eff->SetDirectory(0);
757  eff->SetStats(0);
758  eff->SetFillColor(19);
759  eff->SetFillStyle(0);
760  eff->SetLineColor(9);
761  eff->SetLineWidth(3);
762  eff->SetMarkerStyle(20);
763  eff->GetXaxis()->SetTitleOffset(1.2);
764  // eff->Draw("");
765  return eff;
766 }
767 //________________________________________________________________________________
768 Double_t StdEdxModel::NpCorrection(Double_t betagamma) {
769  Double_t bgL10 = TMath::Log10(betagamma);
770  bgL10 = TMath::Max(-2.0, TMath::Min(-1e-3,bgL10));
771  static TH1D *eff = 0;
772  if (! eff) eff = protonEff();
773  // return 1.03*eff->Interpolate(bgL10);
774  return eff->Interpolate(bgL10);
775 }
776 //________________________________________________________________________________
777 Double_t StdEdxModel::dNdxEff(Double_t poverm, Double_t charge, Double_t mass) {
778  if (!fgStdEdxModel) instance();
779  Double_t bgMC = bgCorrected(poverm);
780  Double_t dNdxMC = dNdx(bgMC, charge);
781  Double_t dNdx = dNdxMC*NpCorrection(poverm);
782  return dNdx;
783 }
784 //________________________________________________________________________________
785 Double_t StdEdxModel::dNdxEffL10func(Double_t *x, Double_t *p) {
786  Double_t bg = TMath::Power(10., x[0]);
787  return instance()->dNdxEff(bg, 1., 0.13956995);
788 }
789 //________________________________________________________________________________
790 TF1 *StdEdxModel::dNdxEffL10F() {
791  TF1 *f = new TF1("dNdxEffL10F",dNdxEffL10func,-2,5,0);
792  return f;
793 }
794 //________________________________________________________________________________
795 Double_t StdEdxModel::extremevalueG(Double_t *x, Double_t *p) {
796  Double_t normL = p[0];
797  Double_t mu = p[1];
798  Double_t sigmaI = p[2];
799  Double_t phase = p[3];
800  Double_t sigmaG = p[4];
801  Double_t t = (mu - x[0])*sigmaI;
802  Double_t frac = TMath::Sin(phase);
803  frac *= frac;
804  return TMath::Exp(normL)*((1. - frac)*TMath::Abs(sigmaI)*TMath::Exp(t - TMath::Exp(t)) + frac*TMath::Gaus(t, 0., sigmaG, kTRUE));
805 }
806 //________________________________________________________________________________
807 TF1 *StdEdxModel::ExValG() {
808  TString fName("ExValG");
809  TF1 *f = (TF1 *) gROOT->GetListOfFunctions()->FindObject(fName);
810  if (! f) {
811  f = new TF1(fName, extremevalueG, -2, 5, 5);
812  f->SetParNames("normL","mu","sigmaI", "phase","sigmaG");
813  f->SetParLimits(2, 0.1, 10.0);
814  f->SetParLimits(3, 0., TMath::PiOver2());
815  // f->SetParLimits(4, 0.1, 1.0);
816  f->FixParameter(4, 1.0);
817  }
818  f->SetParameters(0., 0., 2.5, 0.75, 1.0);
819  return f;
820 }
StElectonsDEV_dEdx
Definition: St_spline3C.h:48
StPionDEV_dEdx
Definition: St_spline3C.h:58
StProtonDEV_dEdx
Definition: St_spline3C.h:68
Generated by   doxygen 1.8.5