next50 trigger mask in AliHLTGlobalEsdConverterComponent
[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-8
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){
75     //printf("AliTRDdEdxBaseUtils::BinLogX warning xmin < epsilon! nothing done, axis not set. %e\n", from);  exit(1);
76     return;
77   }
78
79   Double_t *new_bins = new Double_t[bins + 1];
80
81   new_bins[0] = from;
82   const Double_t factor = pow(to/from, 1./bins);
83
84   for (int i = 1; i <= bins; i++) {
85    new_bins[i] = factor * new_bins[i-1];
86   }
87   axis->Set(bins, new_bins);
88   delete [] new_bins;
89 }
90
91 void AliTRDdEdxBaseUtils::GetCDFCuts(const TH1D *hh, Int_t ncut, Double_t cuts[], const Double_t cdfs[], Double_t thres)
92 {
93   //
94   //counts of hh is sorted
95   //
96
97   for(Int_t ii=0; ii<ncut; ii++){
98     cuts[ii] = -999;
99   }
100
101   Int_t nsel = 0; 
102   const Int_t nbin = hh->GetNbinsX();
103   Double_t datas[nbin];
104   for(Int_t ii=1; ii<=nbin; ii++){
105     const Double_t res = hh->GetBinContent(ii);
106     if(res<thres){
107       continue;
108     }
109
110     datas[nsel] = res;
111     nsel++;
112   }
113   if(!nsel)
114     return;
115
116   Int_t id[nsel];
117   TMath::Sort(nsel, datas, id, kFALSE);
118
119   for(Int_t ii=0; ii<ncut; ii++){
120     const Double_t icdf = cdfs[ii];
121     if(icdf<0 || icdf>1){
122       printf("AliTRDdEdxBaseUtils::GetCDFCuts error cdfs[%d] %15f out of range!\n", ii, icdf); exit(1);
123     }
124     cuts[ii] = datas[id[Int_t(icdf*nsel)]];
125   }
126 }
127
128 Double_t AliTRDdEdxBaseUtils::GetMeanRMS(Double_t nn, Double_t sum, Double_t w2s, Double_t * grms, Double_t * gerr)
129 {
130   //
131   //calculate mean (with error) and rms from sum, w2s, nn
132   //if nn=0, mean, error, and rms all = 0
133   //
134
135   Double_t tmean = 0, trms = 0, terr = 0;
136
137   if(nn>EPSILON){
138     tmean = sum/nn;
139
140     const Double_t arg = w2s/nn-tmean*tmean;
141     if(TMath::Abs(arg)<EPSILON){
142       trms = 0;
143     }
144     else{
145       if( arg <0 ){
146         printf("AliTRDdEdxBaseUtils::GetMeanRMS error negative sqrt argument!! %e -- %e %e %f\n", arg, w2s, sum, nn); exit(1);
147       }
148       
149       trms = TMath::Sqrt(arg);
150     }
151
152     terr = trms/TMath::Sqrt(nn);
153   }
154
155   if(grms){
156     (*grms) = trms;
157   }
158
159   if(gerr){
160     (*gerr) = terr;
161   }
162
163   return tmean;
164 }
165
166 Double_t AliTRDdEdxBaseUtils::TruncatedMean(Int_t nx, const Double_t xdata[], Double_t lowfrac, Double_t highfrac, Double_t * grms, Double_t * gerr, Double_t *wws)
167 {
168   //
169   //calculate truncated mean
170   //return <x*w>_{low-high according to x}
171   //
172
173   /*
174   //test->
175   for(Int_t ii=0; ii<nx; ii++){
176     printf("test %d/%d %f\n", ii, nx, xdata[ii]);
177   }
178   //<--test
179   */
180
181   Int_t index[nx];
182   TMath::Sort(nx, xdata, index, kFALSE);
183
184   Int_t nused = 0;
185   Double_t sum = 0;
186   Double_t w2s = 0;
187   const Int_t istart = Int_t (nx*lowfrac);
188   const Int_t istop  = Int_t (nx*highfrac);
189
190   //=,< correct, because when low=0, high=1 it is correct
191   for(Int_t ii=istart; ii<istop; ii++){
192     Double_t weight = 1;
193     if(wws){
194       weight = wws[index[ii]];
195     }
196     const Double_t sx = xdata[index[ii]]*weight;
197
198     sum += sx;
199     w2s += sx*sx;
200
201     nused++;
202     //printf("test in loop %d/%d %f %f %f\n", ii, nused, sx, sum, w2s);
203     
204   }
205
206   return GetMeanRMS(nused, sum, w2s, grms, gerr);
207 }
208
209 Double_t AliTRDdEdxBaseUtils::TruncatedMean(const TH1 *hh, Double_t lowfrac, Double_t highfrac, Double_t * grms, Double_t * gerr)
210 {
211   //
212   //do truncation on histogram
213   //
214   //if hh is scaled, be sure Sumw2 is called before scaling!! then mean, rms and err will all be correct
215   
216   //with under-/over-flow
217   Double_t npreTrunc = 0;
218   for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
219     const Double_t be = hh->GetBinError(itmp);
220     const Double_t bc = hh->GetBinContent(itmp);
221     if(be<EPSILON){
222       if(bc>EPSILON){
223         printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
224       }
225       continue;
226     }
227     npreTrunc += bc*bc/be/be;
228   }
229
230   const Double_t nstart = npreTrunc*lowfrac;
231   const Double_t nstop = npreTrunc*highfrac;
232
233   //with Double_t this should also handle normalized hist
234   Double_t ntot = 0;
235   Double_t nused = 0;
236   Double_t sum = 0;
237   Double_t w2s = 0;
238   for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
239     const Double_t be = hh->GetBinError(itmp);
240     const Double_t bc = hh->GetBinContent(itmp);
241     if(be<EPSILON){
242       if(bc>EPSILON){
243         printf("AliTRDdEdxBaseUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
244       }
245       continue;
246     }
247     const Double_t weight = bc*bc/be/be;
248     ntot+=weight;
249     //<= correct, because when high=1, nstop has to be included
250     if(ntot>nstart && ntot<=nstop){
251
252       const Double_t val = hh->GetBinCenter(itmp);
253       sum += weight*val;
254       w2s += weight*val*val;
255     
256       nused += weight;
257
258       //printf("test %d %f %f --- %f %f -- %f %f\n", itmp, weight, val, sum, w2s, nused, nsample);
259     }
260     else if(ntot>nstop){
261       if(itmp>=hh->GetNbinsX()){
262         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);
263         nused = 0;
264         w2s = sum = 0;
265       }
266       break;
267     }
268   }
269
270   return GetMeanRMS(nused, sum, w2s, grms, gerr);
271 }
272
273 void AliTRDdEdxBaseUtils::FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, Double_t thres, Double_t lowfrac, Double_t highfrac)
274 {
275   //
276   //fit slices of hh using truncation
277   //
278
279   const Int_t x0 = hh->GetXaxis()->GetFirst();
280   const Int_t x1 = hh->GetXaxis()->GetLast();
281   const Int_t y0 = hh->GetYaxis()->GetFirst();
282   const Int_t y1 = hh->GetYaxis()->GetLast();
283
284   const Int_t nx = hh->GetNbinsX();
285   const Int_t ny = hh->GetNbinsY();
286   const Double_t xmin = hh->GetXaxis()->GetXmin();
287   const Double_t xmax = hh->GetXaxis()->GetXmax();
288   const Double_t ymin = hh->GetYaxis()->GetXmin();
289   const Double_t ymax = hh->GetYaxis()->GetXmax();
290
291   hnor = new TH1D(Form("%s_amp",hh->GetName()), "", nx, xmin, xmax);
292   hmpv = new TH1D(Form("%s_mpv",hh->GetName()), "", nx, xmin, xmax);
293   hwid = new TH1D(Form("%s_wid",hh->GetName()), "", nx, xmin, xmax);
294   hres = new TH1D(Form("%s_res",hh->GetName()), "", nx, xmin, xmax);
295
296   Double_t newbins[ny+1];
297   for(Int_t ii=1; ii<=ny+1; ii++){
298     const Double_t lowx= hh->GetYaxis()->GetBinLowEdge(ii);
299     newbins[ii-1]=lowx;
300   }
301
302   for(Int_t ix=x0; ix<=x1; ix++){
303     //to speed up
304     const Double_t rawcount = hh->Integral(ix,ix,0, ny+1);
305     if(rawcount<EPSILON){
306       continue;
307     }
308
309     TH1D *htmp = new TH1D(Form("FitSlicesY_%s_%d", hh->GetName(), ix),"",ny, ymin, ymax);
310     htmp->GetXaxis()->Set(ny, newbins);
311     Double_t ntot = 0;
312     for(Int_t iy=y0; iy<=y1; iy++){
313       const Double_t be = hh->GetBinError(ix,iy);
314       const Double_t bc = hh->GetBinContent(ix, iy);
315
316       if(be<EPSILON){
317         if(bc>EPSILON){
318           printf("AliTRDdEdxBaseUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1);
319         }
320         continue;
321       }
322
323       htmp->SetBinContent(iy, bc);
324       htmp->SetBinError(iy, be);
325
326       ntot += (bc/be)*(bc/be);
327
328       //if(be) printf("test %d %d : %f %f %f\n", ix, iy, bc, be, pow(bc/be,2));
329     }
330
331     hnor->SetBinContent(ix, ntot);
332     hnor->SetBinError(  ix, 0);
333     
334     if(ntot<thres || htmp->GetRMS()<EPSILON){
335       delete htmp;
336       continue;
337     }
338
339     //test htmp->Draw();
340     Double_t trms = -999, terr = -999;
341     const Double_t tmean = TruncatedMean(htmp, lowfrac, highfrac, &trms, &terr);
342
343     hmpv->SetBinContent(ix, tmean);
344     hmpv->SetBinError(  ix, terr);
345
346     hwid->SetBinContent(ix, trms);
347     hwid->SetBinError(  ix, 0);
348
349     hres->SetBinContent(ix, tmean>EPSILON ? trms/tmean:0);
350     hres->SetBinError(  ix, 0);
351
352     delete htmp;
353   }
354
355   TH1 *hhs[]={hnor, hmpv, hwid, hres};
356   const TString yt[]={"N", "MPV", "#sigma", "#sigma/MPV"};
357   const Int_t nh = sizeof(hhs)/sizeof(TH1*);
358   for(Int_t ii=0; ii<nh; ii++){
359     hhs[ii]->SetYTitle(Form("%s of %s", yt[ii].Data(), hh->GetYaxis()->GetTitle()));
360     hhs[ii]->SetXTitle(hh->GetXaxis()->GetTitle());
361     hhs[ii]->GetYaxis()->SetTitleOffset(hh->GetYaxis()->GetTitleOffset());
362     hhs[ii]->SetTitle(hh->GetTitle());
363   }
364 }
365
366 //===================================================================================
367 //                                TRD Analysis Fast Tool
368 //===================================================================================
369
370 Int_t AliTRDdEdxBaseUtils::GetNtracklet(const AliESDEvent *esd)
371 {
372   //
373   //number of trd tracklet in one esd event
374   //
375   const Int_t ntrk0 = esd->GetNumberOfTracks();
376   Int_t ntrdv1=0;
377   for(Int_t ii=0; ii<ntrk0; ii++){
378     ntrdv1 += esd->GetTrack(ii)->GetTRDntracklets();
379   }
380   return ntrdv1;
381 }
382
383 AliTRDtrackV1 * AliTRDdEdxBaseUtils::GetTRDtrackV1(const AliESDtrack * esdtrack)
384 {
385   //
386   //Get TRD friend track
387   //
388
389   AliESDfriendTrack *  friendtrk = (AliESDfriendTrack *)esdtrack->GetFriendTrack();
390   if(!friendtrk){
391     //printf("xlulog AliAnalysisTaskCosmicTRD::GetTRDtrack no friend!!\n"); exit(1);
392     return 0x0;
393   }
394
395   TObject *calibObject=0x0;
396   AliTRDtrackV1 * trdtrack=0x0;
397   for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) {
398     if( (trdtrack=dynamic_cast<AliTRDtrackV1*>(calibObject)) )
399       break;
400   }
401
402   return trdtrack;
403 }
404
405 Bool_t AliTRDdEdxBaseUtils::IsInSameStack(const AliTRDtrackV1 *trdtrack)
406 {
407   //
408   // to check if all tracklets are in the same stack, useful in cosmic
409   //
410
411   TVectorD secs(18), stks(5);
412
413   for(Int_t ilayer = 0; ilayer < 6; ilayer++){
414     AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
415     if(!tracklet)
416       continue;
417     
418     const Int_t det = tracklet->GetDetector();
419     const Int_t isector = AliTRDgeometry::GetSector(det);
420     const Int_t istack  = AliTRDgeometry::GetStack(det);
421
422     secs[isector] = 1;
423     stks[istack]  = 1;
424  }
425
426   if(secs.Sum()!=1 || stks.Sum()!=1){
427     return kFALSE;
428   }
429   else 
430     return kTRUE;
431 }
432
433 Double_t AliTRDdEdxBaseUtils::GetRedefinedPhi(Double_t phi)
434 {
435   //
436   //redefine phi in pi/2 - pi/2 + 2pi
437   //
438
439   if(phi>0 && phi < TMath::Pi()/2){
440     phi += 2*TMath::Pi();
441   }
442
443   return phi;
444 }
445
446 Double_t AliTRDdEdxBaseUtils::Getdydx(const AliTRDseedV1 *tracklet)
447 {
448   //
449   //get dydx
450   //
451   if(tracklet)
452     return tracklet->GetYref(1);
453   else
454     return -999;
455 }
456
457 Double_t AliTRDdEdxBaseUtils::Getdzdx(const AliTRDseedV1 *tracklet)
458 {
459   //
460   //get dzdx
461   //
462   if(tracklet)
463     return tracklet->GetZref(1);
464   else
465     return -999;
466 }
467
468 Double_t AliTRDdEdxBaseUtils::Getdldx(const AliTRDseedV1 *tracklet)
469 {
470   //
471   //return angular normalization factor = dldx
472   //
473   if(tracklet)
474     return TMath::Sqrt(1+tracklet->GetYref(1)*tracklet->GetYref(1)+tracklet->GetZref(1)*tracklet->GetZref(1));
475   else
476     return -999;
477 }
478
479 AliTRDseedV1 * AliTRDdEdxBaseUtils::GetFirstTracklet(const AliTRDtrackV1 *trdtrack)
480 {
481   //
482   //as function name
483   //
484
485   AliTRDseedV1 *tracklet = 0x0;
486   for(Int_t ilayer = 0; ilayer < 6; ilayer++){
487     tracklet=trdtrack->GetTracklet(ilayer);
488     if(!tracklet)
489       continue;
490     
491     break;
492   }
493
494   return tracklet;
495 }
496
497 AliTRDseedV1 * AliTRDdEdxBaseUtils::GetLastTracklet(const AliTRDtrackV1 *trdtrack)
498 {
499   //
500   //get last tracklet
501   //
502   AliTRDseedV1 *tracklet = 0x0;
503   
504   for(Int_t ilayer = 5; ilayer >= 0; ilayer--){
505     tracklet=trdtrack->GetTracklet(ilayer);
506     if(!tracklet)
507       continue;
508
509     break;
510   }
511
512   return tracklet;
513 }
514
515 void AliTRDdEdxBaseUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom)
516 {
517   //
518   //as function name
519   //
520   isec = istk = -999;
521   mom = -999;
522
523   AliTRDseedV1 *tracklet = GetFirstTracklet(trdtrack);
524   if(tracklet){
525     const Int_t det = tracklet->GetDetector();
526     isec = AliTRDgeometry::GetSector(det);
527     istk = AliTRDgeometry::GetStack(det);
528     
529     mom = tracklet->GetMomentum();
530   }
531
532   return;
533 }
534
535 //===================================================================================
536 //                                Detector and Data Constant 
537 //===================================================================================
538 Int_t  AliTRDdEdxBaseUtils::ToDetector(Int_t gtb)
539 {
540   //
541   //gtb = det*Ntb+itb
542   //
543   return gtb/AliTRDseedV1::kNtb;
544 }
545
546 Int_t  AliTRDdEdxBaseUtils::ToTimeBin(Int_t gtb)
547
548   //
549   //gtb = det*Ntb+itb
550   //
551   return gtb%AliTRDseedV1::kNtb;
552 }
553
554 Int_t  AliTRDdEdxBaseUtils::ToSector(Int_t gtb)
555 {
556   //
557   //return sector
558   //
559   return AliTRDgeometry::GetSector(ToDetector(gtb));
560 }
561
562 Int_t  AliTRDdEdxBaseUtils::ToStack(Int_t gtb)
563 {
564   //
565   //return stack
566   //
567   return AliTRDgeometry::GetStack(ToDetector(gtb));
568 }
569
570 Int_t  AliTRDdEdxBaseUtils::ToLayer(Int_t gtb)
571 {
572   //
573   //return layer
574   //
575   return AliTRDgeometry::GetLayer(ToDetector(gtb));
576 }
577
578 void AliTRDdEdxBaseUtils::CheckRunB(TString listrun1kg, Int_t run, TString & type)
579 {
580   //
581   //check run b field
582   //
583   if(listrun1kg.Contains(Form("%d",run))){
584     type+="1kG"; 
585   }
586   else {
587     type+="5kG";
588   }
589 }
590
591 TString AliTRDdEdxBaseUtils::GetRunType(Int_t run)
592 {
593   //
594   //return run type
595   //
596
597   TString type;
598   if(run>=121527 && run<= 126460)//LHC10d
599     type="2010ppLHC10d";
600   else if(run>=126461 && run<= 130930)//LHC10e
601     type="2010ppLHC10e";
602   else if(run>=136782 && run <= 139846)//LHC10h
603     type="2010PbPbLHC10h";
604   else if(run>= 143224 && run<= 143237)//2011Feb
605     type="2011cosmicFeb";
606   else if(run>= 150587 && run<= 154930){
607     type="2011cosmicMayJun";
608
609     CheckRunB("154601 154602 154629 154634 154636 154639 154643", run,  type);
610   }
611   else if(run>=173047 && run<=180164){
612     type="2012cosmic";
613
614     CheckRunB("173047 173049 173050 173065 173070 173092 173095 173113 173131 173160 174154 174156 174192 174194", run, type);
615   }
616   else{
617     type="unknown";
618   }
619
620   type.ToUpper();
621   return type;
622 }
623
624 void AliTRDdEdxBaseUtils::PrintControl()
625 {
626   //
627   //print out control variable
628   //
629   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());
630 }
631
632 //===================================================================================
633 //                                 dEdx Parameterization
634 //===================================================================================
635 void AliTRDdEdxBaseUtils::FastFitdEdxTR(TH1 * hh)
636 {
637   //
638   //fit dedx tr from mean
639   //
640
641   TF1 *ff=new TF1("ff", MeandEdxTRLogx, -0.5, 4.5, 8);
642   Double_t par[8];
643   par[0]=   2.397001e-01;
644   par[1]=   1.334697e+00;
645   par[2]=   6.967470e+00;
646   par[3]=   9.055289e-02;
647   par[4]=   9.388760e+00;
648   par[5]=   9.452742e-04;
649   par[6]=  -1.866091e+00;
650   par[7]=   1.403545e+00;
651
652   ff->SetParameters(par);
653   hh->Fit(ff,"N");
654
655   for(int ii=0; ii<8; ii++){
656     printf("par[%d]=%e;\n", ii, ff->GetParameter(ii));
657   }
658
659   TFile *fout=new TFile("fastfit.root","recreate");
660   hh->Write();
661   ff->Write();
662   fout->Save();
663   fout->Close();
664
665   delete fout;
666   delete ff;
667 }
668
669 Double_t AliTRDdEdxBaseUtils::Q0MeanTRDpp(Double_t bg)
670 {
671   //
672   //truncated Mean Q_{xx} in TRD
673   //
674  
675   Double_t par[8];
676   Double_t scale = 1;
677   
678   //2012-0605 /u/xlu/.task/CommondEdx/myAnaData/newTune_r56595/inuse/plot/fastFit
679   scale = 0.9257;
680 par[0]=8.316476e-01;
681 par[1]=1.600907e+00;
682 par[2]=7.728447e+00;
683 par[3]=6.235622e-02;
684 par[4]=1.136209e+01;
685 par[5]=-1.495764e-06;
686 par[6]=-2.536119e+00;
687 par[7]=3.416031e+00;
688
689   return   scale*MeandEdxTR(&bg, par);
690 }
691
692 Double_t AliTRDdEdxBaseUtils::Q0MeanTRDPbPb(Double_t bg)
693 {
694   //
695   //truncated Mean Q_{xx} in TRD
696   //
697  
698   Double_t par[8];
699
700   //03132012161259
701   //opt: PbPbQ0
702 par[0]=   1.844912e-01;
703 par[1]=   2.509702e+00;
704 par[2]=   6.744031e+00;
705 par[3]=   7.355123e-02;
706 par[4]=   1.166023e+01;
707 par[5]=   1.736186e-04;
708 par[6]=  -1.716063e+00;
709 par[7]=   1.611366e+00;
710
711   ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see4.log:hhtype4Q0b2c2 scale        0.460994 at ltbg        0.650000  
712   return   0.460994*MeandEdxTR(&bg, par);
713 }
714
715 Double_t AliTRDdEdxBaseUtils::Q1MeanTRDpp(Double_t bg)
716 {
717   //
718   //truncated Mean 1/(1/Q)_{xx} in TRD
719   //
720  
721   Double_t par[8];
722
723   //So 4. Mär 13:30:51 CET 2012
724   //opt: trdppQ1
725   par[0]=   2.434646e-01;
726   par[1]=   1.400211e+00;
727   par[2]=   6.937471e+00;
728   par[3]=   7.758118e-02;
729   par[4]=   1.097372e+01;
730   par[5]=   4.297518e-04;
731   par[6]=  -1.806266e+00;
732   par[7]=   1.543811e+00;
733
734   //hhtype2Q1b2c2 scale        0.418629 at ltbg        0.650000
735
736   return  0.418629*MeandEdxTR(&bg, par);
737 }
738
739 Double_t AliTRDdEdxBaseUtils::Q1MeanTRDPbPb(Double_t bg)
740 {
741   //
742   //truncated Mean 1/(1/Q)_{xx} in TRD
743   //
744  
745   Double_t par[8];
746
747   //So 4. Mär 13:30:52 CET 2012
748   //opt: trdPbPbQ1
749   par[0]=   2.193660e-01;
750   par[1]=   2.051864e+00;
751   par[2]=   6.825112e+00;
752   par[3]=   6.151693e-02;
753   par[4]=   1.390343e+01;
754   par[5]=   6.010032e-05;
755   par[6]=  -1.676324e+00;
756   par[7]=   1.838873e+00;
757
758   //hhtype4Q1b2c2 scale        0.457988 at ltbg        0.650000
759
760   return  0.457988*MeandEdxTR(&bg, par);
761 }
762
763 Double_t AliTRDdEdxBaseUtils::QMeanTPC(Double_t bg)
764 {
765   //
766   //bethe bloch in TPC
767   //
768
769   Double_t par[5];
770   //Mi 15. Feb 14:48:05 CET 2012
771   //train_2012-02-13_1214.12001, tpcsig
772   par[0]=       4.401269;
773   par[1]=       9.725370;
774   par[2]=       0.000178;
775   par[3]=       1.904962;
776   par[4]=       1.426576;
777
778   return MeandEdx(&bg, par);
779 }
780
781 Double_t AliTRDdEdxBaseUtils::MeandEdxTR(const Double_t * xx,  const Double_t * pin)
782 {
783   //
784   //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
785   //npar = 8 = 3+5
786   //
787   Double_t ptr[4]={0};
788   for(int ii=0; ii<3; ii++){
789     ptr[ii+1]=pin[ii];
790   }
791   return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
792 }
793
794 Double_t AliTRDdEdxBaseUtils::MeanTR(const Double_t * xx,  const Double_t * par)
795 {
796   //
797   //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
798   //npar = 4
799   //
800
801   const Double_t bg = xx[0];
802   const Double_t gamma = sqrt(1+bg*bg);
803
804   const Double_t p0 = TMath::Abs(par[1]);
805   const Double_t p1 = TMath::Abs(par[2]);
806   const Double_t p2 = TMath::Abs(par[3]);
807
808   const Double_t zz = TMath::Log(gamma);
809   const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
810
811   return par[0]+tryield;
812 }
813
814 Double_t AliTRDdEdxBaseUtils::MeandEdx(const Double_t * xx,  const Double_t * par)
815 {
816   return ALEPH(xx, par);
817 }
818
819 Double_t AliTRDdEdxBaseUtils::ALEPH(const Double_t * xx,  const Double_t * par)
820 {
821   //
822   //ALEPH parametrization for dEdx
823   //npar = 5
824   //
825
826   const Double_t bg = xx[0];
827   const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
828
829   const Double_t p0 = TMath::Abs(par[0]);
830   const Double_t p1 = TMath::Abs(par[1]);
831
832   //after redefining p2 (<0) -> ln P2
833   //const Double_t p2 = par[2];
834
835   //restore back
836   const Double_t p2 = TMath::Abs(par[2]);
837
838   const Double_t p3 = TMath::Abs(par[3]);
839   const Double_t p4 = TMath::Abs(par[4]);
840
841   const Double_t aa = TMath::Power(beta, p3);
842
843   //numerically very bad when p2~1e-15, bg~1e3, p4~5.
844   const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
845
846   //+1 to avoid the singularity when bg<1 and p4>>1
847   //very^inf important! with this hesse NOTPOS->OK and no nan or inf in migrad
848   //i.e. numerically very stable
849   //the reason is when bg<1 ln blows up very fast with p4>>1 and nan/inf appears in migrad search
850   //the fitting will adjust the parameters as if bg is not shifted, the fitted curve looks fine!!
851   //-- 2012 Aug 8
852   //----> fail for 10h, not used, restore back!
853   //const Double_t lbg = TMath::Log(bg);
854
855   //redefine p2(<0) -> ln P2
856   //corresponds to alephFit.C fAlephOPt = 11
857   //const Double_t bb = p2 + TMath::Log( 1 + TMath::Exp(-p2-p4*lbg) );
858   
859   //printf("test----- %f %f -- %f %f %f %f %f --- %f %f %f\n", bg, beta, p0, p1, p2, p3, p4, p0/aa, aa, bb);
860
861   return (p1-aa-bb)*p0/aa;
862 }
863
864 Double_t AliTRDdEdxBaseUtils::ToLogx(FFunc func, const Double_t * xx,  const Double_t * par)
865 {
866   //
867   //f(x)-> f(y) with y=log10(x)
868   //
869   const Double_t x2[]={TMath::Power(10, xx[0])};
870   return func(x2, par);
871 }
872