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