TRD/AliTRDdEdxBaseUtils.cxx

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