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