* The calib hist binning is changed to log binning
[u/mrichter/AliRoot.git] / TRD / AliTRDdEdxBaseUtils.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //
17 // TRD dEdx base utils
18 // xx
19 // xx
20 // xx
21 // xx
22 //
23 //  Xianguo Lu 
24 //  lu@physi.uni-heidelberg.de
25 //  Xianguo.Lu@cern.ch
26 //  
27 //
28
29 #include "TF1.h"
30 #include "TFile.h"
31 #include "TH1D.h"
32 #include "TH2D.h"
33 #include "THnSparse.h"
34 #include "TMath.h"
35 #include "TMatrixD.h"
36 #include "TMinuit.h"
37 #include "TObjArray.h"
38 #include "TRandom3.h"
39 #include "TStopwatch.h"
40 #include "TVectorD.h"
41
42 #include "TTreeStream.h"
43
44 #include "AliCDBId.h"
45 #include "AliCDBMetaData.h"
46 #include "AliCDBStorage.h"
47 #include "AliESDEvent.h"
48 #include "AliESDfriendTrack.h"
49 #include "AliESDtrack.h"
50 #include "AliTRDcalibDB.h"
51 #include "AliTRDCalROC.h"
52 #include "AliTRDtrackV1.h"
53
54 #include "AliTRDdEdxBaseUtils.h"
55
56 #define EPSILON 1e-12
57
58 Double_t AliTRDdEdxBaseUtils::fgQ0Frac = 0.3;
59 Double_t AliTRDdEdxBaseUtils::fgQ1Frac = 0.5;
60 Double_t AliTRDdEdxBaseUtils::fgTimeBinCountCut = 0.0; 
61 Int_t    AliTRDdEdxBaseUtils::fgCalibTPCnclsCut = 70;
62 Bool_t   AliTRDdEdxBaseUtils::fgExBOn = kTRUE; 
63 Bool_t   AliTRDdEdxBaseUtils::fgPadGainOn = kTRUE;
64 Double_t AliTRDdEdxBaseUtils::fgQScale = 45;
65
66 //===================================================================================
67 //                                   Math and Histogram
68 //===================================================================================
69 void AliTRDdEdxBaseUtils::BinLogX(TAxis *axis) 
70 {
71   //
72   // Method for the correct logarithmic binning of histograms
73   // copied and modified from AliTPCcalibBase
74
75   const Int_t bins = axis->GetNbins();
76
77   const Double_t from = axis->GetXmin();
78   const Double_t to = axis->GetXmax();
79   if (from<EPSILON) return;
80   Double_t *new_bins = new Double_t[bins + 1];
81
82   new_bins[0] = from;
83   const Double_t factor = pow(to/from, 1./bins);
84
85   for (int i = 1; i <= bins; i++) {
86    new_bins[i] = factor * new_bins[i-1];
87   }
88   axis->Set(bins, new_bins);
89   delete [] new_bins;
90 }
91
92 void AliTRDdEdxBaseUtils::GetCDFCuts(const TH1D *hh, const Int_t ncut, Double_t cuts[], const Double_t cdfs[], const Double_t thres)
93 {
94   //
95   //counts of hh is sorted
96   //
97
98   for(Int_t ii=0; ii<ncut; ii++){
99     cuts[ii] = -999;
100   }
101
102   Int_t nsel = 0; 
103   const Int_t nbin = hh->GetNbinsX();
104   Double_t datas[nbin];
105   for(Int_t ii=1; ii<=nbin; ii++){
106     const Double_t res = hh->GetBinContent(ii);
107     if(res<thres){
108       continue;
109     }
110
111     datas[nsel] = res;
112     nsel++;
113   }
114   if(!nsel)
115     return;
116
117   Int_t id[nsel];
118   TMath::Sort(nsel, datas, id, kFALSE);
119
120   for(Int_t ii=0; ii<ncut; ii++){
121     const Double_t icdf = cdfs[ii];
122     if(icdf<0 || icdf>1){
123       printf("AliTRDdEdxBaseUtils::GetCDFCuts error cdfs[%d] %15f out of range!\n", ii, icdf); exit(1);
124     }
125     cuts[ii] = datas[id[Int_t(icdf*nsel)]];
126   }
127 }
128
129 Double_t AliTRDdEdxBaseUtils::GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr)
130 {
131   //
132   //calculate mean (with error) and rms from sum, w2s, nn
133   //if nn=0, mean, error, and rms all = 0
134   //
135
136   Double_t tmean = 0, trms = 0, terr = 0;
137
138   if(nn>EPSILON){
139     tmean = sum/nn;
140
141     const Double_t arg = w2s/nn-tmean*tmean;
142     if(TMath::Abs(arg)<EPSILON){
143       trms = 0;
144     }
145     else{
146       if( arg <0 ){
147         printf("AliTRDdEdxBaseUtils::GetMeanRMS error negative sqrt argument!! %e -- %e %e %f\n", arg, w2s, sum, nn); exit(1);
148       }
149       
150       trms = TMath::Sqrt(arg);
151     }
152
153     terr = trms/TMath::Sqrt(nn);
154   }
155
156   if(grms){
157     (*grms) = trms;
158   }
159
160   if(gerr){
161     (*gerr) = terr;
162   }
163
164   return tmean;
165 }
166
167 Double_t AliTRDdEdxBaseUtils::TruncatedMean(const Int_t nx, const Double_t xdata[], const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr, Double_t *wws)
168 {
169   //
170   //calculate truncated mean
171   //return <x*w>_{low-high according to x}
172   //
173
174   /*
175   //test->
176   for(Int_t ii=0; ii<nx; ii++){
177     printf("test %d/%d %f\n", ii, nx, xdata[ii]);
178   }
179   //<--test
180   */
181
182   Int_t index[nx];
183   TMath::Sort(nx, xdata, index, kFALSE);
184
185   Int_t nused = 0;
186   Double_t sum = 0;
187   Double_t w2s = 0;
188   const Int_t istart = Int_t (nx*lowfrac);
189   const Int_t istop  = Int_t (nx*highfrac);
190
191   //=,< correct, because when low=0, high=1 it is correct
192   for(Int_t ii=istart; ii<istop; ii++){
193     Double_t weight = 1;
194     if(wws){
195       weight = wws[index[ii]];
196     }
197     const Double_t sx = xdata[index[ii]]*weight;
198
199     sum += sx;
200     w2s += sx*sx;
201
202     nused++;
203     //printf("test in loop %d/%d %f %f %f\n", ii, nused, sx, sum, w2s);
204     
205   }
206
207   return GetMeanRMS(nused, sum, w2s, grms, gerr);
208 }
209
210 Double_t AliTRDdEdxBaseUtils::TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr)
211 {
212   //
213   //do truncation on histogram
214   //
215   //if hh is scaled, be sure Sumw2 is called before scaling!! then mean, rms and err will all be correct
216   
217   //with under-/over-flow
218   Double_t npreTrunc = 0;
219   for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
220     const Double_t be = hh->GetBinError(itmp);
221     const Double_t bc = hh->GetBinContent(itmp);
222     if(be<EPSILON){
223       if(bc>EPSILON){
224         printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
225       }
226       continue;
227     }
228     npreTrunc += bc*bc/be/be;
229   }
230
231   const Double_t nstart = npreTrunc*lowfrac;
232   const Double_t nstop = npreTrunc*highfrac;
233
234   //with Double_t this should also handle normalized hist
235   Double_t ntot = 0;
236   Double_t nused = 0;
237   Double_t sum = 0;
238   Double_t w2s = 0;
239   for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
240     const Double_t be = hh->GetBinError(itmp);
241     const Double_t bc = hh->GetBinContent(itmp);
242     if(be<EPSILON){
243       if(bc>EPSILON){
244         printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
245       }
246       continue;
247     }
248     const Double_t weight = bc*bc/be/be;
249     ntot+=weight;
250     //<= correct, because when high=1, nstop has to be included
251     if(ntot>nstart && ntot<=nstop){
252
253       const Double_t val = hh->GetBinCenter(itmp);
254       sum += weight*val;
255       w2s += weight*val*val;
256     
257       nused += weight;
258
259       //printf("test %d %f %f --- %f %f -- %f %f\n", itmp, weight, val, sum, w2s, nused, nsample);
260     }
261     else if(ntot>nstop){
262       if(itmp>=hh->GetNbinsX()){
263         printf("AliTRDdEdxBaseUtils::TruncatedMean warning hist range too small %s %f %f %d %d, %15f %15f %15f; nused w2s sum set to 0\n", hh->GetName(), hh->GetBinLowEdge(1), hh->GetBinLowEdge(itmp), itmp, hh->GetNbinsX(), hh->GetBinContent(hh->GetNbinsX())/hh->Integral(0,hh->GetNbinsX()+1), hh->GetBinContent(hh->GetNbinsX()), hh->Integral(0,hh->GetNbinsX()+1)); //exit(1);
264         nused = 0;
265         w2s = sum = 0;
266       }
267       break;
268     }
269   }
270
271   return GetMeanRMS(nused, sum, w2s, grms, gerr);
272 }
273
274 void AliTRDdEdxBaseUtils::FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, const Double_t thres, const Double_t lowfrac, const Double_t highfrac)
275 {
276   //
277   //fit slices of hh using truncation
278   //
279
280   const Int_t x0 = hh->GetXaxis()->GetFirst();
281   const Int_t x1 = hh->GetXaxis()->GetLast();
282   const Int_t y0 = hh->GetYaxis()->GetFirst();
283   const Int_t y1 = hh->GetYaxis()->GetLast();
284
285   const Int_t nx = hh->GetNbinsX();
286   const Int_t ny = hh->GetNbinsY();
287   const Double_t xmin = hh->GetXaxis()->GetXmin();
288   const Double_t xmax = hh->GetXaxis()->GetXmax();
289   const Double_t ymin = hh->GetYaxis()->GetXmin();
290   const Double_t ymax = hh->GetYaxis()->GetXmax();
291
292   hnor = new TH1D(Form("%s_amp",hh->GetName()), "", nx, xmin, xmax);
293   hmpv = new TH1D(Form("%s_mpv",hh->GetName()), "", nx, xmin, xmax);
294   hwid = new TH1D(Form("%s_wid",hh->GetName()), "", nx, xmin, xmax);
295   hres = new TH1D(Form("%s_res",hh->GetName()), "", nx, xmin, xmax);
296
297   for(Int_t ix=x0; ix<=x1; ix++){
298     //to speed up
299     const Double_t rawcount = hh->Integral(ix,ix,0, ny+1);
300     if(rawcount<EPSILON){
301       continue;
302     }
303
304     TH1D *htmp = new TH1D(Form("FitSlicesY_%s_%d", hh->GetName(), ix),"",ny, ymin, ymax);
305     Double_t ntot = 0;
306     for(Int_t iy=y0; iy<=y1; iy++){
307       const Double_t be = hh->GetBinError(ix,iy);
308       const Double_t bc = hh->GetBinContent(ix, iy);
309
310       if(be<EPSILON){
311         if(bc>EPSILON){
312           printf("AliTRDdEdxBaseUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1);
313         }
314         continue;
315       }
316
317       htmp->SetBinContent(iy, bc);
318       htmp->SetBinError(iy, be);
319
320       ntot += (bc/be)*(bc/be);
321
322       //if(be) printf("test %d %d : %f %f %f\n", ix, iy, bc, be, pow(bc/be,2));
323     }
324
325     hnor->SetBinContent(ix, ntot);
326     hnor->SetBinError(  ix, 0);
327     
328     if(ntot<thres || htmp->GetRMS()<EPSILON){
329       delete htmp;
330       continue;
331     }
332
333     //test htmp->Draw();
334     Double_t trms = -999, terr = -999;
335     const Double_t tmean = TruncatedMean(htmp, lowfrac, highfrac, &trms, &terr);
336
337     hmpv->SetBinContent(ix, tmean);
338     hmpv->SetBinError(  ix, terr);
339
340     hwid->SetBinContent(ix, trms);
341     hwid->SetBinError(  ix, 0);
342
343     hres->SetBinContent(ix, tmean>EPSILON ? trms/tmean:0);
344     hres->SetBinError(  ix, 0);
345
346     delete htmp;
347   }
348
349   TH1 *hhs[]={hnor, hmpv, hwid, hres};
350   const TString yt[]={"N", "MPV", "#sigma", "#sigma/MPV"};
351   const Int_t nh = sizeof(hhs)/sizeof(TH1*);
352   for(Int_t ii=0; ii<nh; ii++){
353     hhs[ii]->SetYTitle(Form("%s of %s", yt[ii].Data(), hh->GetYaxis()->GetTitle()));
354     hhs[ii]->SetXTitle(hh->GetXaxis()->GetTitle());
355     hhs[ii]->GetYaxis()->SetTitleOffset(hh->GetYaxis()->GetTitleOffset());
356     hhs[ii]->SetTitle(hh->GetTitle());
357   }
358 }
359
360 //===================================================================================
361 //                                TRD Analysis Fast Tool
362 //===================================================================================
363
364 Int_t AliTRDdEdxBaseUtils::GetNtracklet(const AliESDEvent *esd)
365 {
366   //
367   //number of trd tracklet in one esd event
368   //
369   const Int_t ntrk0 = esd->GetNumberOfTracks();
370   Int_t ntrdv1=0;
371   for(Int_t ii=0; ii<ntrk0; ii++){
372     ntrdv1 += esd->GetTrack(ii)->GetTRDntracklets();
373   }
374   return ntrdv1;
375 }
376
377 AliTRDtrackV1 * AliTRDdEdxBaseUtils::GetTRDtrackV1(const AliESDtrack * esdtrack)
378 {
379   //
380   //Get TRD friend track
381   //
382
383   AliESDfriendTrack *  friendtrk = (AliESDfriendTrack *)esdtrack->GetFriendTrack();
384   if(!friendtrk){
385     //printf("xlulog AliAnalysisTaskCosmicTRD::GetTRDtrack no friend!!\n"); exit(1);
386     return 0x0;
387   }
388
389   TObject *calibObject=0x0;
390   AliTRDtrackV1 * trdtrack=0x0;
391   for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) {
392     if( (trdtrack=dynamic_cast<AliTRDtrackV1*>(calibObject)) )
393       break;
394   }
395
396   return trdtrack;
397 }
398
399 Bool_t AliTRDdEdxBaseUtils::IsInSameStack(const AliTRDtrackV1 *trdtrack)
400 {
401   //
402   // to check if all tracklets are in the same stack, useful in cosmic
403   //
404
405   TVectorD secs(18), stks(5);
406
407   for(Int_t ilayer = 0; ilayer < 6; ilayer++){
408     AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
409     if(!tracklet)
410       continue;
411     
412     const Int_t det = tracklet->GetDetector();
413     const Int_t isector = AliTRDgeometry::GetSector(det);
414     const Int_t istack  = AliTRDgeometry::GetStack(det);
415
416     secs[isector] = 1;
417     stks[istack]  = 1;
418  }
419
420   if(secs.Sum()!=1 || stks.Sum()!=1){
421     return kFALSE;
422   }
423   else 
424     return kTRUE;
425 }
426
427 Bool_t AliTRDdEdxBaseUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom)
428 {
429   //
430   //as function name
431   //
432   isec = istk = -999;
433   mom = -999;
434
435   for(Int_t ilayer = 0; ilayer < 6; ilayer++){
436     AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
437     if(!tracklet)
438       continue;
439     
440     const Int_t det = tracklet->GetDetector();
441     isec = AliTRDgeometry::GetSector(det);
442     istk = AliTRDgeometry::GetStack(det);
443
444     mom = tracklet->GetMomentum();
445
446     break;
447   }
448
449   if(isec<0)
450     return kFALSE;
451   else 
452     return kTRUE;
453 }
454
455 //===================================================================================
456 //                                Detector and Data Constant 
457 //===================================================================================
458 Int_t  AliTRDdEdxBaseUtils::ToDetector(const Int_t gtb)
459 {
460   //
461   //gtb = det*Ntb+itb
462   //
463   return gtb/AliTRDseedV1::kNtb;
464 }
465
466 Int_t  AliTRDdEdxBaseUtils::ToTimeBin(const Int_t gtb)
467
468   //
469   //gtb = det*Ntb+itb
470   //
471   return gtb%AliTRDseedV1::kNtb;
472 }
473
474 Int_t  AliTRDdEdxBaseUtils::ToSector(const Int_t gtb)
475 {
476   //
477   //return sector
478   //
479   return AliTRDgeometry::GetSector(ToDetector(gtb));
480 }
481
482 Int_t  AliTRDdEdxBaseUtils::ToStack(const Int_t gtb)
483 {
484   //
485   //return stack
486   //
487   return AliTRDgeometry::GetStack(ToDetector(gtb));
488 }
489
490 Int_t  AliTRDdEdxBaseUtils::ToLayer(const Int_t gtb)
491 {
492   //
493   //return layer
494   //
495   return AliTRDgeometry::GetLayer(ToDetector(gtb));
496 }
497
498 TString AliTRDdEdxBaseUtils::GetRunType(const Int_t run)
499 {
500   //
501   //return run type
502   //
503
504   TString type;
505   if(run>=121527 && run<= 126460)//LHC10d
506     type="pp2010LHC10d";
507   else if(run>=126461 && run<= 130930)//LHC10e
508     type="pp2010LHC10e";
509   else if(run>=136782 && run <= 139846)//LHC10h
510     type="PbPb2010LHC10h";
511   else if(run>= 143224 && run<= 143237)//2011Feb
512     type="cosmic2011Feb";
513   else if(run>= 150587 && run<= 154930){
514     type="cosmic2011MayJun";
515
516     TString runstr(Form("%d",run));
517     const TString listrun1kg("154601 154602 154629 154634 154636 154639 154643");
518     if(listrun1kg.Contains(runstr)){
519       type+="1kG";
520     }
521     else{
522       type+="5kG";
523     }      
524   }
525   else{
526     type="unknown";
527   }
528
529   type.ToUpper();
530   return type;
531 }
532
533 void AliTRDdEdxBaseUtils::PrintControl()
534 {
535   //
536   //print out control variable
537   //
538   printf("\nAliTRDdEdxBaseUtils::PrintControl Q0Frac %.1f, Q1Frac %.1f, TimeBinCountCut %.2f, CalibTPCnclsCut %d, IsExBOn %d, IsPadGainOn %d, QScale %.2f\n\n", Q0Frac(), Q1Frac(), TimeBinCountCut(), CalibTPCnclsCut(), IsExBOn(), IsPadGainOn(), QScale());
539 }
540
541 //===================================================================================
542 //                                 dEdx Parameterization
543 //===================================================================================
544
545 Double_t AliTRDdEdxBaseUtils::Q0MeanTRDpp(const Double_t bg)
546 {
547   //
548   //truncated Mean Q_{xx} in TRD
549   //
550  
551   Double_t par[8];
552   //03132012161150
553   //opt: ppQ0
554 par[0]=   2.397001e-01;
555 par[1]=   1.334697e+00;
556 par[2]=   6.967470e+00;
557 par[3]=   9.055289e-02;
558 par[4]=   9.388760e+00;
559 par[5]=   9.452742e-04;
560 par[6]=  -1.866091e+00;
561 par[7]=   1.403545e+00;
562
563   ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see2.log:hhtype2Q0b2c2 scale        0.428543 at ltbg        0.650000
564   return   0.428543*MeandEdxTR(&bg, par);
565 }
566
567 Double_t AliTRDdEdxBaseUtils::Q0MeanTRDPbPb(const Double_t bg)
568 {
569   //
570   //truncated Mean Q_{xx} in TRD
571   //
572  
573   Double_t par[8];
574
575   //03132012161259
576   //opt: PbPbQ0
577 par[0]=   1.844912e-01;
578 par[1]=   2.509702e+00;
579 par[2]=   6.744031e+00;
580 par[3]=   7.355123e-02;
581 par[4]=   1.166023e+01;
582 par[5]=   1.736186e-04;
583 par[6]=  -1.716063e+00;
584 par[7]=   1.611366e+00;
585
586   ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see4.log:hhtype4Q0b2c2 scale        0.460994 at ltbg        0.650000  
587   return   0.460994*MeandEdxTR(&bg, par);
588 }
589
590 Double_t AliTRDdEdxBaseUtils::Q1MeanTRDpp(const Double_t bg)
591 {
592   //
593   //truncated Mean 1/(1/Q)_{xx} in TRD
594   //
595  
596   Double_t par[8];
597
598   //So 4. Mär 13:30:51 CET 2012
599   //opt: trdppQ1
600   par[0]=   2.434646e-01;
601   par[1]=   1.400211e+00;
602   par[2]=   6.937471e+00;
603   par[3]=   7.758118e-02;
604   par[4]=   1.097372e+01;
605   par[5]=   4.297518e-04;
606   par[6]=  -1.806266e+00;
607   par[7]=   1.543811e+00;
608
609   //hhtype2Q1b2c2 scale        0.418629 at ltbg        0.650000
610
611   return  0.418629*MeandEdxTR(&bg, par);
612 }
613
614 Double_t AliTRDdEdxBaseUtils::Q1MeanTRDPbPb(const Double_t bg)
615 {
616   //
617   //truncated Mean 1/(1/Q)_{xx} in TRD
618   //
619  
620   Double_t par[8];
621
622   //So 4. Mär 13:30:52 CET 2012
623   //opt: trdPbPbQ1
624   par[0]=   2.193660e-01;
625   par[1]=   2.051864e+00;
626   par[2]=   6.825112e+00;
627   par[3]=   6.151693e-02;
628   par[4]=   1.390343e+01;
629   par[5]=   6.010032e-05;
630   par[6]=  -1.676324e+00;
631   par[7]=   1.838873e+00;
632
633   //hhtype4Q1b2c2 scale        0.457988 at ltbg        0.650000
634
635   return  0.457988*MeandEdxTR(&bg, par);
636 }
637
638 Double_t AliTRDdEdxBaseUtils::QMeanTPC(const Double_t bg)
639 {
640   //
641   //bethe bloch in TPC
642   //
643
644   Double_t par[5];
645   //Mi 15. Feb 14:48:05 CET 2012
646   //train_2012-02-13_1214.12001, tpcsig
647   par[0]=       4.401269;
648   par[1]=       9.725370;
649   par[2]=       0.000178;
650   par[3]=       1.904962;
651   par[4]=       1.426576;
652
653   return MeandEdx(&bg, par);
654 }
655
656 Double_t AliTRDdEdxBaseUtils::MeandEdxTR(const Double_t * xx, const Double_t * pin)
657 {
658   //
659   //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
660   //npar = 8 = 3+5
661   //
662   Double_t ptr[4]={0};
663   for(int ii=0; ii<3; ii++){
664     ptr[ii+1]=pin[ii];
665   }
666   return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
667 }
668
669 Double_t AliTRDdEdxBaseUtils::MeanTR(const Double_t * xx, const Double_t * par)
670 {
671   //
672   //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
673   //npar = 4
674   //
675
676   const Double_t bg = xx[0];
677   const Double_t gamma = sqrt(1+bg*bg);
678
679   const Double_t p0 = TMath::Abs(par[1]);
680   const Double_t p1 = TMath::Abs(par[2]);
681   const Double_t p2 = TMath::Abs(par[3]);
682
683   const Double_t zz = TMath::Log(gamma);
684   const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
685
686   return par[0]+tryield;
687 }
688
689 Double_t AliTRDdEdxBaseUtils::MeandEdx(const Double_t * xx, const Double_t * par)
690 {
691   //
692   //ALEPH parametrization for dEdx
693   //npar = 5
694   //
695
696   const Double_t bg = xx[0];
697   const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
698
699   const Double_t p0 = TMath::Abs(par[0]);
700   const Double_t p1 = TMath::Abs(par[1]);
701   const Double_t p2 = TMath::Abs(par[2]);
702   const Double_t p3 = TMath::Abs(par[3]);
703   const Double_t p4 = TMath::Abs(par[4]);
704
705   const Double_t aa = TMath::Power(beta, p3);
706   const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
707
708   //printf("test----- %f %f -- %f %f %f %f %f --- %f %f %f\n", bg, beta, p0, p1, p2, p3, p4, p0/aa, aa, bb);
709
710   return (p1-aa-bb)*p0/aa;
711 }
712
713 Double_t AliTRDdEdxBaseUtils::ToLogx(FFunc func, const Double_t * xx, const Double_t * par)
714 {
715   //
716   //f(x)-> f(y) with y=log10(x)
717   //
718   const Double_t x2[]={TMath::Power(10, xx[0])};
719   return func(x2, par);
720 }
721