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