changes from fzhou
[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 Double_t AliTRDdEdxBaseUtils::GetDeltaPhi(const AliTRDseedV1 *tracklet)
430 {
431   //
432   //get phi difference to normal incidence
433   //
434   if(tracklet)
435     return TMath::ATan(tracklet->GetYref(1));
436   else
437     return -999;
438 }
439
440 AliTRDseedV1 * AliTRDdEdxBaseUtils::GetLastTracklet(const AliTRDtrackV1 *trdtrack)
441 {
442   //
443   //get last tracklet
444   //
445   AliTRDseedV1 *tracklet = 0x0;
446   
447   for(Int_t ilayer = 5; ilayer >= 0; ilayer--){
448     tracklet=trdtrack->GetTracklet(ilayer);
449     if(!tracklet)
450       continue;
451
452     break;
453   }
454
455   return tracklet;
456 }
457
458 AliTRDseedV1 * AliTRDdEdxBaseUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom)
459 {
460   //
461   //as function name
462   //
463   isec = istk = -999;
464   mom = -999;
465
466   AliTRDseedV1 *tracklet = 0x0;
467   for(Int_t ilayer = 0; ilayer < 6; ilayer++){
468     tracklet=trdtrack->GetTracklet(ilayer);
469     if(!tracklet)
470       continue;
471     
472     const Int_t det = tracklet->GetDetector();
473     isec = AliTRDgeometry::GetSector(det);
474     istk = AliTRDgeometry::GetStack(det);
475
476     mom = tracklet->GetMomentum();
477
478     break;
479   }
480
481   return tracklet;
482 }
483
484 //===================================================================================
485 //                                Detector and Data Constant 
486 //===================================================================================
487 Int_t  AliTRDdEdxBaseUtils::ToDetector(const Int_t gtb)
488 {
489   //
490   //gtb = det*Ntb+itb
491   //
492   return gtb/AliTRDseedV1::kNtb;
493 }
494
495 Int_t  AliTRDdEdxBaseUtils::ToTimeBin(const Int_t gtb)
496
497   //
498   //gtb = det*Ntb+itb
499   //
500   return gtb%AliTRDseedV1::kNtb;
501 }
502
503 Int_t  AliTRDdEdxBaseUtils::ToSector(const Int_t gtb)
504 {
505   //
506   //return sector
507   //
508   return AliTRDgeometry::GetSector(ToDetector(gtb));
509 }
510
511 Int_t  AliTRDdEdxBaseUtils::ToStack(const Int_t gtb)
512 {
513   //
514   //return stack
515   //
516   return AliTRDgeometry::GetStack(ToDetector(gtb));
517 }
518
519 Int_t  AliTRDdEdxBaseUtils::ToLayer(const Int_t gtb)
520 {
521   //
522   //return layer
523   //
524   return AliTRDgeometry::GetLayer(ToDetector(gtb));
525 }
526
527 TString AliTRDdEdxBaseUtils::GetRunType(const Int_t run)
528 {
529   //
530   //return run type
531   //
532
533   TString type;
534   if(run>=121527 && run<= 126460)//LHC10d
535     type="pp2010LHC10d";
536   else if(run>=126461 && run<= 130930)//LHC10e
537     type="pp2010LHC10e";
538   else if(run>=136782 && run <= 139846)//LHC10h
539     type="PbPb2010LHC10h";
540   else if(run>= 143224 && run<= 143237)//2011Feb
541     type="cosmic2011Feb";
542   else if(run>= 150587 && run<= 154930){
543     type="cosmic2011MayJun";
544
545     TString runstr(Form("%d",run));
546     const TString listrun1kg("154601 154602 154629 154634 154636 154639 154643");
547     if(listrun1kg.Contains(runstr)){
548       type+="1kG";
549     }
550     else{
551       type+="5kG";
552     }      
553   }
554   else{
555     type="unknown";
556   }
557
558   type.ToUpper();
559   return type;
560 }
561
562 void AliTRDdEdxBaseUtils::PrintControl()
563 {
564   //
565   //print out control variable
566   //
567   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());
568 }
569
570 //===================================================================================
571 //                                 dEdx Parameterization
572 //===================================================================================
573 void AliTRDdEdxBaseUtils::FastFitdEdxTR(TH1 * hh)
574 {
575   //
576   //fit dedx tr from mean
577   //
578
579   TF1 *ff=new TF1("ff", MeandEdxTRLogx, -0.5, 4.5, 8);
580   Double_t par[8];
581   par[0]=   2.397001e-01;
582   par[1]=   1.334697e+00;
583   par[2]=   6.967470e+00;
584   par[3]=   9.055289e-02;
585   par[4]=   9.388760e+00;
586   par[5]=   9.452742e-04;
587   par[6]=  -1.866091e+00;
588   par[7]=   1.403545e+00;
589
590   ff->SetParameters(par);
591   hh->Fit(ff,"N");
592
593   for(int ii=0; ii<8; ii++){
594     printf("par[%d]=%e;\n", ii, ff->GetParameter(ii));
595   }
596
597   TFile *fout=new TFile("fastfit.root","recreate");
598   hh->Write();
599   ff->Write();
600   fout->Save();
601   fout->Close();
602
603   delete fout;
604   delete ff;
605 }
606
607 Double_t AliTRDdEdxBaseUtils::Q0MeanTRDpp(const Double_t bg)
608 {
609   //
610   //truncated Mean Q_{xx} in TRD
611   //
612  
613   Double_t par[8];
614   Double_t scale = 1;
615   
616   //2012-0605 /u/xlu/.task/CommondEdx/myAnaData/newTune_r56595/inuse/plot/fastFit
617   scale = 0.9257;
618 par[0]=8.316476e-01;
619 par[1]=1.600907e+00;
620 par[2]=7.728447e+00;
621 par[3]=6.235622e-02;
622 par[4]=1.136209e+01;
623 par[5]=-1.495764e-06;
624 par[6]=-2.536119e+00;
625 par[7]=3.416031e+00;
626
627   return   scale*MeandEdxTR(&bg, par);
628 }
629
630 Double_t AliTRDdEdxBaseUtils::Q0MeanTRDPbPb(const Double_t bg)
631 {
632   //
633   //truncated Mean Q_{xx} in TRD
634   //
635  
636   Double_t par[8];
637
638   //03132012161259
639   //opt: PbPbQ0
640 par[0]=   1.844912e-01;
641 par[1]=   2.509702e+00;
642 par[2]=   6.744031e+00;
643 par[3]=   7.355123e-02;
644 par[4]=   1.166023e+01;
645 par[5]=   1.736186e-04;
646 par[6]=  -1.716063e+00;
647 par[7]=   1.611366e+00;
648
649   ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see4.log:hhtype4Q0b2c2 scale        0.460994 at ltbg        0.650000  
650   return   0.460994*MeandEdxTR(&bg, par);
651 }
652
653 Double_t AliTRDdEdxBaseUtils::Q1MeanTRDpp(const Double_t bg)
654 {
655   //
656   //truncated Mean 1/(1/Q)_{xx} in TRD
657   //
658  
659   Double_t par[8];
660
661   //So 4. Mär 13:30:51 CET 2012
662   //opt: trdppQ1
663   par[0]=   2.434646e-01;
664   par[1]=   1.400211e+00;
665   par[2]=   6.937471e+00;
666   par[3]=   7.758118e-02;
667   par[4]=   1.097372e+01;
668   par[5]=   4.297518e-04;
669   par[6]=  -1.806266e+00;
670   par[7]=   1.543811e+00;
671
672   //hhtype2Q1b2c2 scale        0.418629 at ltbg        0.650000
673
674   return  0.418629*MeandEdxTR(&bg, par);
675 }
676
677 Double_t AliTRDdEdxBaseUtils::Q1MeanTRDPbPb(const Double_t bg)
678 {
679   //
680   //truncated Mean 1/(1/Q)_{xx} in TRD
681   //
682  
683   Double_t par[8];
684
685   //So 4. Mär 13:30:52 CET 2012
686   //opt: trdPbPbQ1
687   par[0]=   2.193660e-01;
688   par[1]=   2.051864e+00;
689   par[2]=   6.825112e+00;
690   par[3]=   6.151693e-02;
691   par[4]=   1.390343e+01;
692   par[5]=   6.010032e-05;
693   par[6]=  -1.676324e+00;
694   par[7]=   1.838873e+00;
695
696   //hhtype4Q1b2c2 scale        0.457988 at ltbg        0.650000
697
698   return  0.457988*MeandEdxTR(&bg, par);
699 }
700
701 Double_t AliTRDdEdxBaseUtils::QMeanTPC(const Double_t bg)
702 {
703   //
704   //bethe bloch in TPC
705   //
706
707   Double_t par[5];
708   //Mi 15. Feb 14:48:05 CET 2012
709   //train_2012-02-13_1214.12001, tpcsig
710   par[0]=       4.401269;
711   par[1]=       9.725370;
712   par[2]=       0.000178;
713   par[3]=       1.904962;
714   par[4]=       1.426576;
715
716   return MeandEdx(&bg, par);
717 }
718
719 Double_t AliTRDdEdxBaseUtils::MeandEdxTR(const Double_t * xx, const Double_t * pin)
720 {
721   //
722   //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
723   //npar = 8 = 3+5
724   //
725   Double_t ptr[4]={0};
726   for(int ii=0; ii<3; ii++){
727     ptr[ii+1]=pin[ii];
728   }
729   return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
730 }
731
732 Double_t AliTRDdEdxBaseUtils::MeanTR(const Double_t * xx, const Double_t * par)
733 {
734   //
735   //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
736   //npar = 4
737   //
738
739   const Double_t bg = xx[0];
740   const Double_t gamma = sqrt(1+bg*bg);
741
742   const Double_t p0 = TMath::Abs(par[1]);
743   const Double_t p1 = TMath::Abs(par[2]);
744   const Double_t p2 = TMath::Abs(par[3]);
745
746   const Double_t zz = TMath::Log(gamma);
747   const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
748
749   return par[0]+tryield;
750 }
751
752 Double_t AliTRDdEdxBaseUtils::MeandEdx(const Double_t * xx, const Double_t * par)
753 {
754   //
755   //ALEPH parametrization for dEdx
756   //npar = 5
757   //
758
759   const Double_t bg = xx[0];
760   const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
761
762   const Double_t p0 = TMath::Abs(par[0]);
763   const Double_t p1 = TMath::Abs(par[1]);
764   const Double_t p2 = TMath::Abs(par[2]);
765   const Double_t p3 = TMath::Abs(par[3]);
766   const Double_t p4 = TMath::Abs(par[4]);
767
768   const Double_t aa = TMath::Power(beta, p3);
769   const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
770
771   //printf("test----- %f %f -- %f %f %f %f %f --- %f %f %f\n", bg, beta, p0, p1, p2, p3, p4, p0/aa, aa, bb);
772
773   return (p1-aa-bb)*p0/aa;
774 }
775
776 Double_t AliTRDdEdxBaseUtils::ToLogx(FFunc func, const Double_t * xx, const Double_t * par)
777 {
778   //
779   //f(x)-> f(y) with y=log10(x)
780   //
781   const Double_t x2[]={TMath::Power(10, xx[0])};
782   return func(x2, par);
783 }
784