TRD/AliTRDdEdxUtils.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 // class to calculate TRD dEdx
  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
  30 #include "TF1.h"
  31 #include "TFile.h"
  32 #include "TH1D.h"
  33 #include "TH2D.h"
  34 #include "THnSparse.h"
  35 #include "TMath.h"
  36 #include "TMatrixD.h"
  37 #include "TMinuit.h"
  38 #include "TObjArray.h"
  39 #include "TRandom3.h"
  40 #include "TStopwatch.h"
  41 #include "TVectorD.h"
  42
  43 #include "TTreeStream.h"
  44
  45 #include "AliCDBId.h"
  46 #include "AliCDBMetaData.h"
  47 #include "AliCDBStorage.h"
  48 #include "AliESDEvent.h"
  49 #include "AliESDfriendTrack.h"
  50 #include "AliESDtrack.h"
  51 #include "AliTRDtrackV1.h"
  52
  53 #include "AliTRDdEdxUtils.h"
  54
  55 #define EPSILON 1e-12
  56
  57 THnSparseD * AliTRDdEdxUtils::fgHistGain=0x0;
  58 THnSparseD * AliTRDdEdxUtils::fgHistT0=0x0;
  59 THnSparseD * AliTRDdEdxUtils::fgHistVd=0x0;
  60 TObjArray * AliTRDdEdxUtils::fgHistPHQ=new TObjArray(8);
  61
  62 TString AliTRDdEdxUtils::fgCalibFile;
  63 TObjArray * AliTRDdEdxUtils::fgObjGain = 0x0;
  64 TObjArray * AliTRDdEdxUtils::fgObjT0 = 0x0;
  65 TObjArray * AliTRDdEdxUtils::fgObjVd = 0x0;
  66 TObjArray * AliTRDdEdxUtils::fgObjPHQ = 0x0;
  67
  68 Int_t AliTRDdEdxUtils::fgNchamber = -999;
  69 Double_t AliTRDdEdxUtils::fgChamberQ[6];
  70 Double_t AliTRDdEdxUtils::fgChamberTmean[6];
  71
  72 Double_t AliTRDdEdxUtils::fgTrackTmean = -999;
  73
  74 //===================================================================================
  75 //                                   Math and Histogram
  76 //===================================================================================
  77 void AliTRDdEdxUtils::GetCDFCuts(const TH1D *hh, const Int_t ncut, Double_t cuts[], const Double_t cdfs[], const Double_t thres)
  78 {
  79   //
  80   //counts of hh is sorted
  81   //
  82
  83   for(Int_t ii=0; ii<ncut; ii++){
  84     cuts[ii] = -999;
  85   }
  86
  87   Int_t nsel = 0;
  88   const Int_t nbin = hh->GetNbinsX();
  89   Double_t datas[nbin];
  90   for(Int_t ii=1; ii<=nbin; ii++){
  91     const Double_t res = hh->GetBinContent(ii);
  92     if(res<thres){
  93       continue;
  94     }
  95
  96     datas[nsel] = res;
  97     nsel++;
  98   }
  99   if(!nsel)
 100     return;
 101
 102   Int_t id[nsel];
 103   TMath::Sort(nsel, datas, id, kFALSE);
 104
 105   for(Int_t ii=0; ii<ncut; ii++){
 106     const Double_t icdf = cdfs[ii];
 107     if(icdf<0 || icdf>1){
 108       printf("AliTRDdEdxUtils::GetCDFCuts error cdfs[%d] %15f out of range!\n", ii, icdf); exit(1);
 109     }
 110     cuts[ii] = datas[id[Int_t(icdf*nsel)]];
 111   }
 112 }
 113
 114 Double_t AliTRDdEdxUtils::GetMeanRMS(const Double_t nn, const Double_t sum, const Double_t w2s, Double_t * grms, Double_t * gerr)
 115 {
 116   //
 117   //calculate mean (with error) and rms from sum, w2s, nn
 118   //if nn=0, mean, error, and rms all = 0
 119   //
 120
 121   Double_t tmean = 0, trms = 0, terr = 0;
 122
 123   if(nn>EPSILON){
 124     tmean = sum/nn;
 125
 126     const Double_t arg = w2s/nn-tmean*tmean;
 127     if(TMath::Abs(arg)<EPSILON){
 128       trms = 0;
 129     }
 130     else{
 131       if( arg <0 ){
 132         printf("AliTRDdEdxUtils::GetMeanRMS error negative sqrt argument!! %e -- %e %e %f\n", arg, w2s, sum, nn); exit(1);
 133       }
 134
 135       trms = TMath::Sqrt(arg);
 136     }
 137
 138     terr = trms/TMath::Sqrt(nn);
 139   }
 140
 141   if(grms){
 142     (*grms) = trms;
 143   }
 144
 145   if(gerr){
 146     (*gerr) = terr;
 147   }
 148
 149   return tmean;
 150 }
 151
 152 Double_t AliTRDdEdxUtils::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)
 153 {
 154   //
 155   //calculate truncated mean
 156   //return <x*w>_{low-high according to x}
 157   //
 158
 159   /*
 160   //test->
 161   for(Int_t ii=0; ii<nx; ii++){
 162     printf("test %d/%d %f\n", ii, nx, xdata[ii]);
 163   }
 164   //<--test
 165   */
 166
 167   Int_t index[nx];
 168   TMath::Sort(nx, xdata, index, kFALSE);
 169
 170   Int_t nused = 0;
 171   Double_t sum = 0;
 172   Double_t w2s = 0;
 173   const Int_t istart = Int_t (nx*lowfrac);
 174   const Int_t istop  = Int_t (nx*highfrac);
 175
 176   //=,< correct, because when low=0, high=1 it is correct
 177   for(Int_t ii=istart; ii<istop; ii++){
 178     Double_t weight = 1;
 179     if(wws){
 180       weight = wws[index[ii]];
 181     }
 182     const Double_t sx = xdata[index[ii]]*weight;
 183
 184     sum += sx;
 185     w2s += sx*sx;
 186
 187     nused++;
 188     //printf("test in loop %d/%d %f %f %f\n", ii, nused, sx, sum, w2s);
 189
 190   }
 191
 192   return GetMeanRMS(nused, sum, w2s, grms, gerr);
 193 }
 194
 195 Double_t AliTRDdEdxUtils::TruncatedMean(const TH1 *hh, const Double_t lowfrac, const Double_t highfrac, Double_t * grms, Double_t * gerr)
 196 {
 197   //
 198   //do truncation on histogram
 199   //
 200   //if hh is scaled, be sure Sumw2 is called before scaling!! then mean, rms and err will all be correct
 201
 202   //with under-/over-flow
 203   Double_t npreTrunc = 0;
 204   for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
 205     const Double_t be = hh->GetBinError(itmp);
 206     const Double_t bc = hh->GetBinContent(itmp);
 207     if(be<EPSILON){
 208       if(bc>EPSILON){
 209         printf("AliTRDdEdxUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
 210       }
 211       continue;
 212     }
 213     npreTrunc += bc*bc/be/be;
 214   }
 215
 216   const Double_t nstart = npreTrunc*lowfrac;
 217   const Double_t nstop = npreTrunc*highfrac;
 218
 219   //with Double_t this should also handle normalized hist
 220   Double_t ntot = 0;
 221   Double_t nused = 0;
 222   Double_t sum = 0;
 223   Double_t w2s = 0;
 224   for(Int_t itmp=0; itmp<=hh->GetNbinsX()+1; itmp++){
 225     const Double_t be = hh->GetBinError(itmp);
 226     const Double_t bc = hh->GetBinContent(itmp);
 227     if(be<EPSILON){
 228       if(bc>EPSILON){
 229         printf("AliTRDdEdxUtils::TruncatedMean (hist) error %e %e %d\n", bc, be, itmp); exit(1);
 230       }
 231       continue;
 232     }
 233     const Double_t weight = bc*bc/be/be;
 234     ntot+=weight;
 235     //<= correct, because when high=1, nstop has to be included
 236     if(ntot>nstart && ntot<=nstop){
 237
 238       const Double_t val = hh->GetBinCenter(itmp);
 239       sum += weight*val;
 240       w2s += weight*val*val;
 241
 242       nused += weight;
 243
 244       //printf("test %d %f %f --- %f %f -- %f %f\n", itmp, weight, val, sum, w2s, nused, nsample);
 245     }
 246     else if(ntot>nstop){
 247       if(itmp>=hh->GetNbinsX()){
 248         printf("AliTRDdEdxUtils::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);
 249         nused = 0;
 250         w2s = sum = 0;
 251       }
 252       break;
 253     }
 254   }
 255
 256   return GetMeanRMS(nused, sum, w2s, grms, gerr);
 257 }
 258
 259 void AliTRDdEdxUtils::FitSlicesY(const TH2D *hh, TH1D *&hnor, TH1D *&hmpv, TH1D *&hwid, TH1D *&hres, const Double_t thres, const Double_t lowfrac, const Double_t highfrac)
 260 {
 261   //
 262   //fit slices of hh using truncation
 263   //
 264
 265   const Int_t x0 = hh->GetXaxis()->GetFirst();
 266   const Int_t x1 = hh->GetXaxis()->GetLast();
 267   const Int_t y0 = hh->GetYaxis()->GetFirst();
 268   const Int_t y1 = hh->GetYaxis()->GetLast();
 269
 270   const Int_t nx = hh->GetNbinsX();
 271   const Int_t ny = hh->GetNbinsY();
 272   const Double_t xmin = hh->GetXaxis()->GetXmin();
 273   const Double_t xmax = hh->GetXaxis()->GetXmax();
 274   const Double_t ymin = hh->GetYaxis()->GetXmin();
 275   const Double_t ymax = hh->GetYaxis()->GetXmax();
 276
 277   hnor = new TH1D(Form("%s_amp",hh->GetName()), "", nx, xmin, xmax);
 278   hmpv = new TH1D(Form("%s_mpv",hh->GetName()), "", nx, xmin, xmax);
 279   hwid = new TH1D(Form("%s_wid",hh->GetName()), "", nx, xmin, xmax);
 280   hres = new TH1D(Form("%s_res",hh->GetName()), "", nx, xmin, xmax);
 281
 282   for(Int_t ix=x0; ix<=x1; ix++){
 283     //to speed up
 284     const Double_t rawcount = hh->Integral(ix,ix,0, ny+1);
 285     if(rawcount<EPSILON){
 286       continue;
 287     }
 288
 289     TH1D *htmp = new TH1D(Form("FitSlicesY_%s_%d", hh->GetName(), ix),"",ny, ymin, ymax);
 290     Double_t ntot = 0;
 291     for(Int_t iy=y0; iy<=y1; iy++){
 292       const Double_t be = hh->GetBinError(ix,iy);
 293       const Double_t bc = hh->GetBinContent(ix, iy);
 294
 295       if(be<EPSILON){
 296         if(bc>EPSILON){
 297           printf("AliTRDdEdxUtils::FitSlicesY error %d %d %e %e\n", ix, iy, be, bc); exit(1);
 298         }
 299         continue;
 300       }
 301
 302       htmp->SetBinContent(iy, bc);
 303       htmp->SetBinError(iy, be);
 304
 305       ntot += (bc/be)*(bc/be);
 306
 307       //if(be) printf("test %d %d : %f %f %f\n", ix, iy, bc, be, pow(bc/be,2));
 308     }
 309
 310     hnor->SetBinContent(ix, ntot);
 311     hnor->SetBinError(  ix, 0);
 312
 313     if(ntot<thres || htmp->GetRMS()<EPSILON){
 314       delete htmp;
 315       continue;
 316     }
 317
 318     //test htmp->Draw();
 319     Double_t trms = -999, terr = -999;
 320     const Double_t tmean = TruncatedMean(htmp, lowfrac, highfrac, &trms, &terr);
 321
 322     hmpv->SetBinContent(ix, tmean);
 323     hmpv->SetBinError(  ix, terr);
 324
 325     hwid->SetBinContent(ix, trms);
 326     hwid->SetBinError(  ix, 0);
 327
 328     hres->SetBinContent(ix, tmean>EPSILON ? trms/tmean:0);
 329     hres->SetBinError(  ix, 0);
 330
 331     delete htmp;
 332   }
 333
 334   TH1 *hhs[]={hnor, hmpv, hwid, hres};
 335   const TString yt[]={"N", "MPV", "#sigma", "#sigma/MPV"};
 336   const Int_t nh = sizeof(hhs)/sizeof(TH1*);
 337   for(Int_t ii=0; ii<nh; ii++){
 338     hhs[ii]->SetYTitle(Form("%s of %s", yt[ii].Data(), hh->GetYaxis()->GetTitle()));
 339     hhs[ii]->SetXTitle(hh->GetXaxis()->GetTitle());
 340     hhs[ii]->GetYaxis()->SetTitleOffset(hh->GetYaxis()->GetTitleOffset());
 341     hhs[ii]->SetTitle(hh->GetTitle());
 342   }
 343 }
 344
 345 //===================================================================================
 346 //                                TRD Analysis Fast Tool
 347 //===================================================================================
 348
 349 Int_t AliTRDdEdxUtils::GetNtracklet(const AliESDEvent *esd)
 350 {
 351   //
 352   //number of trd tracklet in one esd event
 353   //
 354   const Int_t ntrk0 = esd->GetNumberOfTracks();
 355   Int_t ntrdv1=0;
 356   for(Int_t ii=0; ii<ntrk0; ii++){
 357     ntrdv1 += esd->GetTrack(ii)->GetTRDntracklets();
 358   }
 359   return ntrdv1;
 360 }
 361
 362 AliTRDtrackV1 * AliTRDdEdxUtils::GetTRDtrackV1(const AliESDtrack * esdtrack)
 363 {
 364   //
 365   //Get TRD friend track
 366   //
 367
 368   AliESDfriendTrack *  friendtrk = (AliESDfriendTrack *)esdtrack->GetFriendTrack();
 369   if(!friendtrk){
 370     //printf("xlulog AliAnalysisTaskCosmicTRD::GetTRDtrack no friend!!\n"); exit(1);
 371     return 0x0;
 372   }
 373
 374   TObject *calibObject=0x0;
 375   AliTRDtrackV1 * trdtrack=0x0;
 376   for(Int_t l=0; (calibObject=friendtrk->GetCalibObject(l)); l++) {
 377     if( (trdtrack=dynamic_cast<AliTRDtrackV1*>(calibObject)) )
 378       break;
 379   }
 380
 381   return trdtrack;
 382 }
 383
 384 Bool_t AliTRDdEdxUtils::IsInSameStack(const AliTRDtrackV1 *trdtrack)
 385 {
 386   //
 387   // to check if all tracklets are in the same stack, useful in cosmic
 388   //
 389
 390   TVectorD secs(18), stks(5);
 391
 392   for(Int_t ilayer = 0; ilayer < 6; ilayer++){
 393     AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
 394     if(!tracklet)
 395       continue;
 396
 397     const Int_t det = tracklet->GetDetector();
 398     const Int_t isector = AliTRDgeometry::GetSector(det);
 399     const Int_t istack  = AliTRDgeometry::GetStack(det);
 400
 401     secs[isector] = 1;
 402     stks[istack]  = 1;
 403  }
 404
 405   if(secs.Sum()!=1 || stks.Sum()!=1){
 406     return kFALSE;
 407   }
 408   else
 409     return kTRUE;
 410 }
 411
 412 Bool_t AliTRDdEdxUtils::GetFirstSectorStackMomentum(const AliTRDtrackV1 *trdtrack, Int_t & isec, Int_t & istk, Double_t & mom)
 413 {
 414   //
 415   //as function name
 416   //
 417   isec = istk = -999;
 418   mom = -999;
 419
 420   for(Int_t ilayer = 0; ilayer < 6; ilayer++){
 421     AliTRDseedV1 *tracklet=trdtrack->GetTracklet(ilayer);
 422     if(!tracklet)
 423       continue;
 424
 425     const Int_t det = tracklet->GetDetector();
 426     isec = AliTRDgeometry::GetSector(det);
 427     istk = AliTRDgeometry::GetStack(det);
 428
 429     mom = tracklet->GetMomentum();
 430
 431     break;
 432   }
 433
 434   if(isec<0)
 435     return kFALSE;
 436   else
 437     return kTRUE;
 438 }
 439
 440 //===================================================================================
 441 //                                Calibration
 442 //===================================================================================
 443 Double_t AliTRDdEdxUtils::GetCalibTPCscale(const Int_t tpcncls, const Double_t tpcsig)
 444 {
 445   //
 446   //the scale used in calibration
 447   //
 448
 449   if(tpcncls < CalibTPCnclsCut())
 450     return -999;
 451
 452   if(tpcsig<EPSILON)
 453     return -999;
 454
 455   return tpcsig/120;
 456
 457 }
 458
 459 Int_t AliTRDdEdxUtils::GetPHQIterator(const Bool_t kinvq, const Double_t mag, const Int_t charge)
 460 {
 461   //
 462   //iterator for calib obj and hist
 463   //
 464   return kinvq*4 + (mag>0)*2 + (charge>0);
 465 }
 466
 467 TObjArray * AliTRDdEdxUtils::GetObjPHQ(const Bool_t kinvq, const Double_t mag, const Int_t charge)
 468 {
 469   //
 470   //return calib obj
 471   //
 472   if(!fgObjPHQ){
 473     printf("AliTRDdEdxUtils::GetObjPHQ error fgObjPHQ null!!\n"); exit(1);
 474   }
 475
 476   return (TObjArray*) fgObjPHQ->At(GetPHQIterator(kinvq, mag, charge));
 477 }
 478
 479 THnSparseD * AliTRDdEdxUtils::GetHistPHQ(const Bool_t kinvq, const Double_t mag, const Int_t charge)
 480 {
 481   //
 482   //return calib hist
 483   //
 484   return (THnSparseD*) fgHistPHQ->At(GetPHQIterator(kinvq, mag, charge));
 485 }
 486
 487 TString AliTRDdEdxUtils::GetPHQName(const Bool_t kobj, const Int_t iter)
 488 {
 489   //
 490   //get name of calib obj/hist of PHQ
 491   //
 492   return Form("TRDCalib%sPHQ%d", kobj?"Obj":"Hist", iter);
 493 }
 494
 495 void AliTRDdEdxUtils::DeleteCalibObj()
 496 {
 497   //
 498   //delete calib obj
 499   //
 500   delete fgObjGain;
 501   delete fgObjT0;
 502   delete fgObjVd;
 503
 504   fgObjGain = 0x0;
 505   fgObjT0 = 0x0;
 506   fgObjVd = 0x0;
 507
 508   if(fgObjPHQ){
 509     fgObjPHQ->SetOwner();
 510     delete fgObjPHQ;
 511     fgObjPHQ = 0x0;
 512   }
 513 }
 514
 515 Bool_t AliTRDdEdxUtils::GenerateDefaultPHQOCDB(const TString path)
 516 {
 517   //
 518   //generate default OCDB object PHQ, do like
 519   //AliTRDdEdxUtils::GenerateDefaultPHQOCDB("local://./")
 520   //
 521
 522   TObjArray * arr8 = new TObjArray(8);
 523   arr8->SetOwner();
 524
 525   for(Int_t ii=0; ii<8; ii++){
 526     TObjArray * arr1 = new TObjArray(1);
 527     arr1->SetOwner();
 528     arr1->SetName(GetPHQName(1, ii));
 529
 530     const Int_t nbins = NTRDtimebin();
 531     TVectorD * vec = new TVectorD(nbins);
 532     for(Int_t jj=0; jj<nbins; jj++){
 533       (*vec)[jj] = 1;
 534     }
 535     arr1->Add(vec);
 536     arr8->Add(arr1);
 537   }
 538
 539   AliCDBMetaData *metaData= new AliCDBMetaData();
 540   metaData->SetObjectClassName("TObjArray");
 541   metaData->SetResponsible("Raphaelle Bailhache and Xianguo Lu");
 542
 543   AliCDBId id1("TRD/Calib/PHQ", 0, 999999999);
 544   AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage(path);
 545   gStorage->Put(arr8, id1, metaData);
 546
 547   delete metaData;
 548   delete arr8;
 549
 550   return kTRUE;
 551 }
 552
 553 void AliTRDdEdxUtils::IniCalibObj()
 554 {
 555   //
 556   //set CalibObj from file, clone to static calib obj
 557   //
 558
 559   DeleteCalibObj();
 560
 561   TFile *cfile=new TFile(fgCalibFile);
 562   if(!cfile->IsOpen()){
 563     printf("AliTRDdEdxUtils::IniCalibObj error fgCalibFile not open! %s\n", fgCalibFile.Data());exit(1);
 564   }
 565
 566   printf("\nAliTRDdEdxUtils::IniCalibObj file: %s\n", fgCalibFile.Data());
 567
 568   //---
 569   const TString objnames[] ={"TRDCalibObjGain", "TRDCalibObjT0", "TRDCalibObjVd"};
 570   TObjArray ** gobjs[]={           &fgObjGain,        &fgObjT0,        &fgObjVd};
 571
 572   const Int_t nobj = sizeof(objnames)/sizeof(TString);
 573   for(Int_t iobj=0; iobj<nobj; iobj++){
 574     TObjArray *tmpo=0x0;
 575     cfile->GetObject(objnames[iobj], tmpo);
 576     if(!tmpo){
 577       printf("AliTRDdEdxUtils::IniCalibObj error obj %s not found!\n", objnames[iobj].Data()); exit(1);
 578     }
 579
 580     (*gobjs[iobj])=(TObjArray*)tmpo->Clone();
 581     (*gobjs[iobj])->SetOwner();
 582   }
 583
 584   fgObjPHQ = new TObjArray(8);
 585   for(Int_t iter=0; iter<8; iter++){
 586     const TString objn = GetPHQName(1, iter);
 587     TObjArray *tmpo=0x0;
 588     cfile->GetObject(objn, tmpo);
 589     if(!tmpo){
 590       printf("AliTRDdEdxUtils::IniCalibObj error obj %s not found!\n", objn.Data()); exit(1);
 591     }
 592
 593     TObjArray *obji=(TObjArray*) tmpo->Clone();
 594     obji->SetOwner();
 595     fgObjPHQ->AddAt(obji, iter);
 596   }
 597
 598   //---
 599
 600   cfile->Close();
 601   delete cfile;
 602 }
 603
 604 void AliTRDdEdxUtils::DeleteCalibHist()
 605 {
 606   //
 607   //delete calib hist
 608   //
 609   delete fgHistGain;
 610   delete fgHistT0;
 611   delete fgHistVd;
 612
 613   fgHistGain = 0x0;
 614   fgHistT0 = 0x0;
 615   fgHistVd = 0x0;
 616
 617   //fgHistPHQ owns the hists
 618   fgHistPHQ->SetOwner();
 619   fgHistPHQ->Clear();
 620 }
 621
 622 void AliTRDdEdxUtils::IniCalibHist(TList *list, const Bool_t kPHQonly)
 623 {
 624   //
 625   //initialize calib hist, list should not own the hist, or list->Clear/delete hist should not be called
 626   //
 627
 628   DeleteCalibHist();
 629
 630   Int_t nbin[2];
 631   const Double_t xmin[2]={0, 0};
 632   Double_t xmax[2];
 633
 634   nbin[0]=NTRDtimebin(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=20;
 635   for(Int_t iter=0; iter<8; iter++){
 636     const TString hn = GetPHQName(0, iter);
 637     THnSparseD *hi = new THnSparseD(hn, "", 2, nbin, xmin, xmax);
 638     //fgHistPHQ owns the hists
 639     fgHistPHQ->AddAt(hi, iter);
 640     list->Add(hi);
 641   }
 642
 643   if(kPHQonly)
 644     return;
 645
 646   nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=20;                 fgHistGain = new THnSparseD("TRDCalibHistGain", "", 2, nbin, xmin, xmax);
 647   nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=AliTRDseedV1::kNtb; fgHistT0   = new THnSparseD("TRDCalibHistT0",   "", 2, nbin, xmin, xmax);
 648   nbin[0]=NTRDchamber(); nbin[1]= 11250; xmax[0]=nbin[0]; xmax[1]=AliTRDseedV1::kNtb; fgHistVd   = new THnSparseD("TRDCalibHistVd",   "", 2, nbin, xmin, xmax);
 649
 650   list->Add(fgHistGain);
 651   list->Add(fgHistT0);
 652   list->Add(fgHistVd);
 653 }
 654
 655 Bool_t AliTRDdEdxUtils::ReadCalibHist(const TString filename, const TString listname)
 656 {
 657   //
 658   //used in AliTRDPreprocessorOffline
 659   //read in calib hist from file, only for PHQ
 660   //
 661   DeleteCalibHist();
 662
 663   //maybe already open by others... don't close
 664   TFile fcalib(filename);
 665
 666   TObjArray * array = (TObjArray*)fcalib.Get(listname);
 667
 668   for(Int_t iter=0; iter<8; iter++){
 669     const TString hn = GetPHQName(0, iter);
 670     THnSparseD * tmph=0x0;
 671     if(array){
 672       tmph = (THnSparseD *) array->FindObject(hn);
 673     }
 674     else{
 675       tmph = (THnSparseD *) fcalib.Get(hn);
 676     }
 677     if(!tmph){
 678       printf("AliTRDdEdxUtils::ReadCalibHist warning calib hist not found! %s %s\n", filename.Data(), listname.Data());
 679       fcalib.ls();
 680       if(array){
 681         array->ls();
 682       }
 683       return kFALSE;
 684     }
 685     THnSparseD *hi = (THnSparseD*)tmph->Clone();
 686     fgHistPHQ->AddAt(hi, iter);
 687   }
 688
 689   return kTRUE;
 690 }
 691
 692 void AliTRDdEdxUtils::FillCalibHist(const Int_t ncls, const TVectorD *arrayQ, const TVectorD *arrayX, THnSparseD * hcalib, const Double_t scale)
 693 {
 694   //
 695   //fill calibration hist
 696   //
 697   if(!hcalib){printf("AliTRDdEdxUtils::FillCalibHist errro hcalib null!!\n"); exit(1);}
 698
 699   for(Int_t ii=0; ii<ncls; ii++){
 700     const Double_t dq = (*arrayQ)[ii];
 701     const Double_t xx = (*arrayX)[ii];
 702
 703     const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
 704     const Double_t xmin = hcalib->GetAxis(0)->GetXmin();
 705     const Double_t xmax = hcalib->GetAxis(0)->GetXmax();
 706
 707     if(xx>=xmax || xx<xmin){
 708       printf("AliTRDdEdxUtils::FillCalibHist error x overflow or underflow! %s %15f %15f %15f\n", hcalib->GetName(),  xx, xmin, xmax); exit(1);
 709     }
 710
 711     const Double_t var[]={xx, TMath::Min(dq, qmax)/scale};
 712     hcalib->Fill(var);
 713   }
 714 }
 715
 716 void AliTRDdEdxUtils::FillCalibHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale)
 717 {
 718   //
 719   //get cluster Q and fill calib hist, if kinvq = kTRUE, 1/Q is filled
 720   //
 721
 722   THnSparseD * hcalib = AliTRDdEdxUtils::GetHistPHQ(kinvq, mag, charge);
 723
 724   TVectorD arrayQ(200), arrayX(200);
 725   const Int_t ncls = AliTRDdEdxUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, trdv1);
 726   FillCalibHist(ncls, &arrayQ, &arrayX, hcalib, kinvq ? 1/scale : scale);
 727
 728   static Int_t kprint = 100;
 729   if(kprint<0){
 730     printf("\nAliTRDdEdxUtils::FillCalibHist summary: \n");
 731     printf("\nkinvq= %d;\n", kinvq);
 732     for(Int_t iq=0; iq<ncls; iq++){
 733       printf("arrayX[%3d] = %15.0f; arrayQ[%3d] =  %15f;\n", iq, arrayX[iq], iq, arrayQ[iq]);
 734     }
 735     printf("\n");
 736   }
 737   kprint++;
 738 }
 739
 740 Int_t AliTRDdEdxUtils::ApplyCalib(const Int_t nc0, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
 741 {
 742   //
 743   //apply calibration on arrayQ
 744   //
 745   if(!cobj){ printf("AliTRDdEdxUtils::ApplyCalib error gain array null!!\n"); exit(1);}
 746
 747   TVectorD tmpq(arrayQ->GetNrows());
 748   TVectorD tmpx(arrayX->GetNrows());
 749   Int_t ncls = 0;
 750
 751   const TVectorD * gain = (TVectorD*) cobj->At(0);
 752   for(Int_t ii=0; ii<nc0; ii++){
 753     const Double_t dq = (*arrayQ)[ii];
 754     const Int_t xx = (Int_t)(*arrayX)[ii];
 755     const Double_t gg = (*gain)[xx];
 756
 757     if(gg<EPSILON){
 758       continue;
 759     }
 760
 761     tmpq[ncls] = dq*gg;
 762     tmpx[ncls] = xx;
 763     ncls++;
 764   }
 765
 766   (*arrayQ)=tmpq;
 767   (*arrayX)=tmpx;
 768
 769   return ncls;
 770 }
 771
 772 void AliTRDdEdxUtils::GetPHCountMeanRMS(const TH1D *hnor, TH1D *&hmean)
 773 {
 774   //
 775   //calculate from the ph calib hist the (mean-3sigma) ph-count in the chamber, save in the TH1D output
 776   //
 777   const Int_t ndet = 540;
 778   TObjArray *obj=new TObjArray(ndet);
 779   obj->SetOwner();
 780   for(Int_t ii=0; ii<ndet; ii++){
 781     obj->Add(new TVectorD(AliTRDseedV1::kNtb));
 782   }
 783
 784   //ibin = binlowedge of bin(ibin+1) = the number fills this bin
 785   for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
 786     const Double_t stat = hnor->GetBinContent(ibin+1);
 787     if(stat<EPSILON){
 788       continue;
 789     }
 790
 791     const Int_t idet = ToDetector(ibin);
 792     const Int_t itb  = ToTimeBin(ibin);
 793     TVectorD *vec=(TVectorD *)obj->At(idet);
 794     (*vec)[itb] = stat;
 795   }
 796
 797   hmean = new TH1D(Form("%sdetmean", hnor->GetName()), "", hnor->GetNbinsX(), hnor->GetXaxis()->GetXmin(), hnor->GetXaxis()->GetXmax());
 798   for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
 799     const Int_t idet = ToDetector(ibin);
 800     const TVectorD *vec=(TVectorD *)obj->At(idet);
 801
 802     Int_t nonzero = 0;
 803     for(Int_t ii=0; ii<vec->GetNrows(); ii++){
 804       if((*vec)[ii]>EPSILON){
 805         nonzero++;
 806       }
 807     }
 808
 809     Double_t mean = 0;
 810     const Double_t lowfrac = 0.6;
 811     //if there are too many 0's, reject this chamber by settig mean=rms=0
 812     if(nonzero> (AliTRDseedV1::kNtb*(1-lowfrac)) ){
 813       //only highest (1-lowfrac)*31 timebins are used to estimate the mean and rms! important! otherwise the 0' will make rms very large!
 814       mean = TruncatedMean(AliTRDseedV1::kNtb, vec->GetMatrixArray(), lowfrac, 1);
 815     }
 816
 817     hmean->SetBinContent(ibin+1, mean);
 818   }
 819
 820   delete obj;
 821 }
 822
 823 void AliTRDdEdxUtils::CalibOutput(const TList *lin, Int_t run)
 824 {
 825   //
 826   //produce calibration objects
 827   //
 828
 829   TString objnames("TRDCalibHistGain TRDCalibHistT0 TRDCalibHistVd ");
 830   for(Int_t iter=0; iter<8; iter++){
 831     objnames+= GetPHQName(0, iter)+" ";
 832   }
 833
 834   TList * lout = new TList;
 835   lout->SetOwner();
 836
 837   TTreeSRedirector *calibStream = new TTreeSRedirector(Form("TRDCalibStream_%010d.root", run));
 838
 839   const Int_t nh=lin->GetEntries();
 840   for(Int_t ii=0; ii<nh; ii++){
 841     const THnSparseD *hh=(THnSparseD*)lin->At(ii);
 842     const TString hname = hh->GetName();
 843     if(!objnames.Contains(hname))
 844       continue;
 845
 846     TObjArray * cobj0 = GetCalibObj(hh, run, lout, calibStream);
 847     lout->Add(cobj0);
 848   }
 849
 850   //lout->ls();
 851
 852   //=============================================================
 853   //=============================================================
 854
 855   TFile *fout=new TFile(Form("TRDCalibObj_%010d.root", run),"recreate");
 856   fout->cd();
 857   const Int_t nout=lout->GetEntries();
 858   for(Int_t ii=0; ii<nout; ii++){
 859     const TString oname = lout->At(ii)->GetName();
 860     if(oname.Contains("Obj")){
 861       TObjArray * cobj = (TObjArray*) lout->At(ii);
 862       cobj->Write(oname, TObjArray::kSingleKey);
 863     }
 864   }
 865   fout->Save();
 866   fout->Close();
 867   delete fout;
 868
 869   fout=new TFile(Form("TRDCalibList_%010d.root", run),"recreate");
 870   fout->cd();
 871   lin->Write();
 872   lout->Write();
 873   fout->Save();
 874   fout->Close();
 875   delete fout;
 876
 877   delete calibStream;
 878
 879   /*
 880     http://root.cern.ch/root/html/TH1.html
 881     When an histogram is created, a reference to it is automatically added to the list of in-memory objects for the current file or directory. This default behaviour can be changed by:
 882
 883     h->SetDirectory(0);          for the current histogram h
 884     TH1::AddDirectory(kFALSE);   sets a global switch disabling the reference
 885
 886     When the histogram is deleted, the reference to it is removed from the list of objects in memory. When a file is closed, all histograms in memory associated with this file are automatically deleted.
 887   */
 888   delete lout;
 889 }
 890
 891 TObjArray* AliTRDdEdxUtils::GetCalibObj(const THnSparseD *hh, Int_t run, TList *lout, TTreeSRedirector *calibStream)
 892 {
 893   //
 894   //produce calibration objects
 895   //
 896
 897   const TString hname = hh->GetName();
 898   const Bool_t kinvq = TString(hname(hname.First('Q')+1,1)).Atoi()&4;
 899
 900   //----------------------------------------
 901   //               Define nbin, tag, cobj0
 902   //----------------------------------------
 903   Int_t nbin =-999;
 904   if(hname.Contains("Gain") || hname.Contains("T0") || hname.Contains("Vd")){
 905     nbin = NTRDchamber();
 906   }
 907   else if(hname.Contains("PHQ")){
 908     nbin = NTRDtimebin();
 909   }
 910   else{
 911     printf("AliTRDdEdxUtils::GetCalibObj error wrong hname!! %s\n", hname.Data()); exit(1);
 912   }
 913
 914   TString tag(hname);
 915   tag.ReplaceAll("Hist","Obj");
 916
 917   TObjArray * cobj0 = new TObjArray(1);
 918   cobj0->SetOwner();
 919   cobj0->SetName(tag);
 920   cobj0->Add(new TVectorD(nbin));
 921
 922   //----------------------------------------
 923   //               Define lowFrac, highFrac
 924   //----------------------------------------
 925   Double_t lowFrac = -999, highFrac = -999;
 926   if(hname.Contains("Gain") || (hname.Contains("PHQ") && !kinvq) ){
 927     lowFrac = 0.01; highFrac = Q0Frac();
 928   }
 929   else if(hname.Contains("PHQ") && kinvq){
 930     lowFrac = Q1Frac(); highFrac = 0.99;
 931   }
 932   else{
 933     lowFrac = 0.01;
 934     highFrac = 0.99;
 935   }
 936
 937   //----------------------------------------
 938   //               Get analysis result
 939   //----------------------------------------
 940   TH1::AddDirectory(kFALSE);//important!
 941   TH1D *hnor=0x0, *hmpv=0x0, *hres=0x0, *hwid=0x0, *htrdphmean = 0x0;//if(!lout), these have to be deleted
 942   TH2D *hpj = hh->Projection(1,0);
 943   FitSlicesY(hpj, hnor, hmpv, hwid, hres, 0, lowFrac, highFrac);
 944   if(hname.Contains("PHQ")){
 945     GetPHCountMeanRMS(hnor, htrdphmean);
 946     if(lout){
 947       lout->Add(htrdphmean);
 948     }
 949   }
 950   delete hpj;
 951
 952   if(lout){
 953     lout->Add(hnor);
 954     lout->Add(hmpv);
 955     lout->Add(hwid);
 956     lout->Add(hres);
 957   }
 958
 959   //----------------------------------------
 960   //               Define Counter
 961   //----------------------------------------
 962   TVectorD *countDet=0x0;
 963   TObjArray *countSSL=0x0;
 964
 965   if(hname.Contains("PHQ") && !kinvq){
 966     countDet=new TVectorD(540);
 967     countSSL=new TObjArray(90);//SectorStackLayer
 968     countSSL->SetOwner();
 969     for(Int_t ii=0; ii<90; ii++){
 970       countSSL->Add(new TVectorD(6));
 971     }
 972   }
 973
 974   //----------------------------------------
 975   //               Fill cobj0
 976   //----------------------------------------
 977
 978   //ibin = binlowedge of bin(ibin+1) = the number fills this bin
 979   for(Int_t ibin=0; ibin<nbin; ibin++){
 980     Double_t gnor = hnor->GetBinContent(ibin+1);
 981     Double_t gmpv = hmpv->GetBinContent(ibin+1);
 982     Double_t gwid = hwid->GetBinContent(ibin+1);
 983     Double_t gres = hres->GetBinContent(ibin+1);
 984
 985     //--- set additional cut by kpass
 986     Bool_t kpass = kTRUE;
 987     Double_t gtrdphmean = -999;
 988     if(htrdphmean){
 989       gtrdphmean = htrdphmean->GetBinContent(ibin+1);
 990       //chamber no statistics (e.g. too many 0's), not usual, not seen in run 143237
 991       if(gtrdphmean<EPSILON){
 992         kpass = kFALSE;
 993       }
 994       if(gnor<TimeBinCountCut()*gtrdphmean){
 995         kpass = kFALSE;
 996       }
 997     }
 998
 999     //--- set calibration constant p0
1000     Double_t p0= 0;
1001
1002     //reason for gmpv=0:
1003     //1)gnor<=3; truncation in hist: (0, 0.6*ntot=1.8 with ntot=3]={1}, in hist entries can pile up so that ntot=2, or 3, and (ntot>nstart && ntot<=nstop) is skipped;
1004     //2)TruncatedMean(hist) out of range (only for Q0, not Q1).
1005
1006     if(gmpv>EPSILON && kpass){
1007       if(tag.Contains("T0")){
1008         p0 = gmpv;
1009       }
1010       else{
1011         p0 = 1/gmpv;
1012       }
1013       //printf("outcalibobj%s %d %15.6e\n", tag.Data(), ibin, p0);
1014     }
1015
1016     (*( (TVectorD*)cobj0->At(0) ))[ibin] = p0;
1017
1018     //--- save optional record
1019     if(p0>EPSILON && countDet && countSSL){
1020       const Int_t idet = ToDetector(ibin);
1021       (*countDet)[idet]=1;
1022
1023       const Int_t isector = ToSector(ibin);
1024       const Int_t istack = ToStack(ibin);
1025       const Int_t ilayer = ToLayer(ibin);
1026       TVectorD * vecsectorstack = (TVectorD*)countSSL->At(istack*18+isector);
1027       (*vecsectorstack)[ilayer]=1;
1028     }
1029
1030     if(calibStream){
1031       (*calibStream)<<tag<<
1032         "run="<<run<<
1033         "p0="<<p0<<
1034
1035         "nor="<<gnor<<
1036         "mpv="<<gmpv<<
1037         "wid="<<gwid<<
1038         "res="<<gres<<
1039         "gtrdphmean="<<gtrdphmean<<
1040
1041         "ibin="<<ibin<<
1042         "\n";
1043     }
1044   }
1045
1046   //----------------------------------------
1047   //               Status Report
1048   //----------------------------------------
1049   if(countDet && countSSL){
1050     TVectorD count2Dstack(90);
1051     for(Int_t ii=0; ii<90; ii++){
1052       TVectorD * vecsectorstack = (TVectorD*)countSSL->At(ii);
1053       const Int_t nlayer = (Int_t)vecsectorstack->Sum();
1054       if(nlayer==6){
1055         count2Dstack[ii]=1;
1056       }
1057     }
1058
1059     printf("\nAliTRDdEdxUtils::GetCalibObj Summary run: %d name: %s entries: %.0f ndetector: %03.0f n2dstack %02.0f\n\n", run, hname.Data(), hh->GetEntries(), countDet->Sum(), count2Dstack.Sum());
1060   }
1061
1062   //----------------------------------------
1063   //               Clean Up
1064   //----------------------------------------
1065
1066   TH1D **hhs[]={&hnor, &hmpv, &hwid, &hres, &htrdphmean};
1067   const Int_t nhh=sizeof(hhs)/sizeof(TH1D**);
1068   for(Int_t ihh=0; ihh<nhh; ihh++){
1069     if(!lout){
1070       delete (*hhs[ihh]);
1071     }
1072   }
1073
1074   delete countDet;
1075   delete countSSL;
1076
1077   //----------------------------------------
1078
1079   return cobj0;
1080 }
1081
1082 //===================================================================================
1083 //                                   dEdx calculation
1084 //===================================================================================
1085 Double_t AliTRDdEdxUtils::GetSignal(const Int_t nch, const Int_t ncls, const Double_t qq)
1086 {
1087   //
1088   //signal = nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq)
1089   //
1090   if(ncls>1e3){
1091     printf("AliTRDdEdxUtils::GetSignal error ncls>1e3! %d\n", ncls); exit(1);
1092   }
1093
1094   return nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq);
1095 }
1096
1097 Int_t AliTRDdEdxUtils::GetNch(const Double_t signal)
1098 {
1099   //
1100   //signal = nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq)
1101   //
1102   return Int_t(signal/1e6);
1103
1104 }
1105
1106 Int_t AliTRDdEdxUtils::GetNcls(const Double_t signal)
1107 {
1108   //
1109   //signal = Nch*1e6 + Ncls*1e3 + (Q>=1e3? 999 : Q)
1110   //
1111
1112   return Int_t ( (signal-GetNch(signal)*1e6)/1e3 );
1113 }
1114
1115 Double_t AliTRDdEdxUtils::GetQ(const Double_t signal)
1116 {
1117   //
1118   //signal = nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq)
1119   //
1120
1121   return signal-GetNch(signal)*1e6 - GetNcls(signal)*1e3;
1122 }
1123
1124 Double_t AliTRDdEdxUtils::ToyCook(const Bool_t kinvq, Int_t &ncluster, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
1125 {
1126   //
1127   //template for cookdedx
1128   //
1129   if(cobj){
1130     if(arrayQ && arrayX){
1131       ncluster = ApplyCalib(ncluster, arrayQ, arrayX, cobj);
1132     }
1133     else{
1134       printf("AliTRDdEdxUtils::ToyCook arrayQ arrayX null, applycalib can not be applied!\n"); exit(1);
1135     }
1136   }
1137
1138   Double_t lowFrac =-999, highFrac = -999;
1139   if(kinvq){
1140     lowFrac = Q1Frac(); highFrac = 0.99;
1141   }
1142   else{
1143     lowFrac = 0.01; highFrac = Q0Frac();
1144   }
1145
1146   Double_t meanQ = TruncatedMean(ncluster, arrayQ->GetMatrixArray(), lowFrac, highFrac);
1147   if(kinvq){
1148     if(meanQ>EPSILON){
1149       meanQ = 1/meanQ;
1150     }
1151   }
1152
1153   return meanQ;
1154 }
1155
1156 Double_t AliTRDdEdxUtils::CombineddEdx(const Bool_t kinvq, Int_t &concls, TVectorD *coarrayQ, TVectorD *coarrayX, const Int_t tpcncls, const TVectorD *tpcarrayQ, const TVectorD *tpcarrayX, const Int_t trdncls, const TVectorD *trdarrayQ, const TVectorD *trdarrayX)
1157 {
1158   //
1159   //combine tpc and trd dedx
1160   //
1161
1162   for(Int_t iq=0; iq<tpcncls; iq++){
1163     (*coarrayQ)[iq]=(*tpcarrayQ)[iq];
1164     if(tpcarrayX && trdarrayX && coarrayX){
1165       (*coarrayX)[iq]=(*tpcarrayX)[iq];
1166     }
1167   }
1168   for(Int_t iq=0; iq<trdncls; iq++){
1169     (*coarrayQ)[tpcncls+iq]=(*trdarrayQ)[iq];
1170     if(tpcarrayX && trdarrayX && coarrayX){
1171       (*coarrayX)[tpcncls+iq]=159+(*trdarrayX)[iq];
1172     }
1173   }
1174
1175   concls=trdncls+tpcncls;
1176
1177   const Double_t coQ = ToyCook(kinvq, concls, coarrayQ, coarrayX);
1178
1179   return coQ;
1180 }
1181
1182
1183 //===================================================================================
1184 //                                   dEdx Getter and Setter
1185 //===================================================================================
1186 Double_t AliTRDdEdxUtils::GetAngularCorrection(const AliTRDseedV1 *seed)
1187 {
1188   //
1189   //return angular normalization factor
1190   //
1191
1192   return TMath::Sqrt(1+seed->GetYref(1)*seed->GetYref(1)+seed->GetZref(1)*seed->GetZref(1));
1193 }
1194
1195 Double_t AliTRDdEdxUtils::GetClusterQ(const Bool_t kinvq, const AliTRDseedV1 * seed, const Int_t itb)
1196 {
1197   //
1198   //get cluster charge
1199   //
1200   Double_t dq = 0;
1201   const AliTRDcluster *cl = 0x0;
1202
1203   cl = seed->GetClusters(itb);                    if(cl) dq+=cl->GetRawQ();
1204   cl = seed->GetClusters(itb+AliTRDseedV1::kNtb); if(cl) dq+=cl->GetRawQ();
1205
1206   dq /= GetAngularCorrection(seed);
1207
1208   dq /= 45.;
1209
1210   if(kinvq){
1211     if(dq>EPSILON){
1212       dq = 1/dq;
1213     }
1214   }
1215
1216   return dq;
1217 }
1218
1219 Int_t AliTRDdEdxUtils::GetArrayClusterQ(const Bool_t kinvq, TVectorD *arrayQ, TVectorD *arrayX, const AliTRDtrackV1 *trdtrack, Int_t timeBin0, Int_t timeBin1, Int_t tstep)
1220 {
1221   //
1222   //return nclustter
1223   //(if kinvq, return 1/q array), size of array must be larger than 31*6
1224   //
1225   if(!arrayQ || arrayQ->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
1226     printf("AliTRDdEdxUtils::GetArrayClusterQ error arrayQ null or size too small! %d\n", arrayQ? arrayQ->GetNrows() : -999); exit(1);
1227   }
1228   if(!arrayX || arrayX->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
1229     printf("AliTRDdEdxUtils::GetArrayClusterQ error arrayX null or size too small! %d\n", arrayX? arrayX->GetNrows() : -999); exit(1);
1230   }
1231
1232   const Int_t mintb = 0;
1233   const Int_t maxtb = AliTRDseedV1::kNtb-1;
1234   if(timeBin0<mintb) timeBin0=mintb;
1235   if(timeBin1>maxtb) timeBin1=maxtb;
1236   if(tstep<=0) tstep=1;
1237
1238   //============
1239   Int_t tbN=0;
1240   Double_t tbQ[200];
1241   Int_t tbBin[200];
1242
1243   for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){
1244     const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
1245     if(!seed)
1246       continue;
1247
1248     const Int_t det = seed->GetDetector();
1249
1250     for(Int_t itb=timeBin0; itb<=timeBin1; itb+=tstep){
1251       const Double_t dq = GetClusterQ(kinvq, seed, itb);
1252       if(dq<EPSILON)
1253         continue;
1254
1255       const Int_t gtb = det * AliTRDseedV1::kNtb + itb;
1256
1257       tbQ[tbN]=dq;
1258       tbBin[tbN]=gtb;
1259       tbN++;
1260     }
1261   }
1262
1263   Int_t ncls = 0;
1264   for(Int_t iq=0; iq<tbN; iq++){
1265     if(tbQ[iq]<EPSILON)
1266       continue;
1267
1268     (*arrayQ)[ncls] = tbQ[iq];
1269     (*arrayX)[ncls] = tbBin[iq];
1270
1271     ncls++;
1272   }
1273
1274   static Int_t kprint = 100;
1275   if(kprint<0){
1276     printf("\nAliTRDdEdxUtils::GetArrayClusterQ raw cluster-Q\n");
1277     for(Int_t iq=0; iq<ncls; iq++){
1278       const Int_t ichamber =  ToLayer((*arrayX)[iq]);
1279       const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
1280       if(!seed){
1281         printf("error seed null!!\n"); exit(1);
1282       }
1283       const Double_t rawq =  (*arrayQ)[iq] * 45. * GetAngularCorrection(seed);
1284       printf("esdid=%d; chamber=%d; timebin=%d; rawq= %.3f; myq[%d]= %e;\n", trdtrack->GetESDid(), ichamber, ToTimeBin((*arrayX)[iq]), rawq, iq, (*arrayQ)[iq]);
1285     }
1286     printf("\n");
1287   }
1288   kprint++;
1289
1290   return ncls;
1291 }
1292
1293 Int_t AliTRDdEdxUtils::UpdateArrayX(const Int_t ncls, TVectorD* arrayX)
1294 {
1295   //
1296   //arrayX det*Ntb+itb -> itb
1297   //
1298
1299   TVectorD countChamber(6);
1300   for(Int_t ii = 0; ii<ncls; ii++){
1301     const Int_t xx = (Int_t)(*arrayX)[ii];
1302     const Int_t idet = ToDetector(xx);
1303
1304     const Double_t ich = AliTRDgeometry::GetLayer(idet);
1305     const Double_t itb = ToTimeBin(xx);
1306     (*arrayX)[ii] = ich*AliTRDseedV1::kNtb+itb;
1307
1308     countChamber[ich] = 1;
1309   }
1310
1311   const Double_t nch = countChamber.Sum();
1312   return (Int_t) nch;
1313 }
1314
1315 void AliTRDdEdxUtils::SetChamberQT(const AliTRDtrackV1 *trdtrack, const Int_t kcalib, THnSparseD * hgain, THnSparseD * ht0, THnSparseD * hvd)
1316 {
1317   //
1318   //CookdEdx at TRD track level, use chamber info, related calibrations: chamber-gain; T0, Vd based on raw PH distribution
1319   //
1320
1321   static Int_t kprint = 100;
1322
1323   fgNchamber = 0;
1324   for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){
1325     //initialize output, default values: 0, so that summation and weighting will automatically discard default quantities
1326     fgChamberQ[ichamber] = fgChamberTmean[ichamber] = 0;
1327
1328     const AliTRDseedV1 *seed = trdtrack->GetTracklet(ichamber);
1329     if (!seed)
1330       continue;
1331
1332     const Int_t idet = seed->GetDetector();
1333
1334     //-------------------------------------------------------------------------
1335
1336     Double_t qsum = 0, qtsum = 0, w2sum = 0;
1337     for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++){
1338       const Double_t dq = GetClusterQ(0, seed, itb);
1339       if(dq<EPSILON)
1340         continue;
1341
1342       qsum += dq;
1343       qtsum += dq*itb;
1344       w2sum += dq*itb*itb;
1345     }
1346     if(qsum<EPSILON)
1347       continue;
1348
1349     //-------------------------------------------------------------------------
1350
1351     Double_t tbm, tbr = 0;
1352     tbm = GetMeanRMS(qsum, qtsum, w2sum, &tbr);
1353
1354     qsum /= 1.25e3/45.;
1355
1356     if(hgain){
1357       const Double_t var[]={idet, qsum};
1358       hgain->Fill(var);
1359     }
1360     if(ht0){
1361       const Double_t var[]={idet, tbm};
1362       ht0->Fill(var);
1363     }
1364     if(hvd){
1365       const Double_t var[]={idet, tbr};
1366       hvd->Fill(var);
1367     }
1368
1369     Double_t gain = 1, t0 = 0, vd = 1;
1370     if(kcalib){
1371       if(!fgObjGain) {printf("AliTRDdEdxUtils::SetChamberQT error Gain array null!!\n"); exit(1);}
1372       if(!  fgObjT0) {printf("AliTRDdEdxUtils::SetChamberQT error T0   array null!!\n"); exit(1);}
1373       if(!  fgObjVd) {printf("AliTRDdEdxUtils::SetChamberQT error Vd   array null!!\n"); exit(1);}
1374
1375       const TVectorD * gainvec = (TVectorD*) fgObjGain->At(0); gain = (*gainvec)[idet];
1376       const TVectorD *   t0vec = (TVectorD*)   fgObjT0->At(0);   t0 = (*  t0vec)[idet];
1377       const TVectorD *   vdvec = (TVectorD*)   fgObjVd->At(0);   vd = (*  vdvec)[idet];
1378     }
1379     if(kprint<0){
1380       printf("\nAliTRDdEdxUtils::CookdEdxV2\n");
1381       printf("idet = %d;\n", idet);
1382       printf("gain = %15f; t0 = %15f; vd = %15f;\n", gain, t0, vd);
1383       printf("\n");
1384     }
1385
1386     qsum *= gain;
1387     tbm = (tbm-t0)*vd;
1388
1389     if(qsum<EPSILON)
1390       continue;
1391
1392     //-------------------------------------------------------------------------
1393
1394     //should have non-zero value, initialized with default 0 (except for calibrated tbm, may be very close to 0)
1395     fgChamberQ[ichamber]  = qsum;
1396     fgChamberTmean[ichamber] = tbm;
1397     fgNchamber++;
1398   }
1399
1400   if(kprint<0){
1401     printf("\nAliTRDdEdxUtils::CookdEdxV2 summary:\n");
1402
1403     printf("\nfgNchamber = %d\n", fgNchamber);
1404     for(Int_t ich=0; ich<AliTRDtrackV1::kNplane; ich++){
1405       printf("fgChamberTmean[%d] = %15f; fgChamberQ[%d] = %15f;\n", ich, fgChamberTmean[ich], ich, fgChamberQ[ich]);
1406     }
1407   }
1408
1409   fgTrackTmean = -999;
1410   if(fgNchamber){
1411     fgTrackTmean = 0;
1412     for(Int_t ich=0; ich<AliTRDtrackV1::kNplane; ich++){
1413       fgTrackTmean += fgChamberTmean[ich];
1414     }
1415     fgTrackTmean /= fgNchamber;
1416   }
1417
1418   if(kprint<0){
1419     printf("\nAliTRDdEdxUtils::CookdEdxV2\n");
1420     printf("GetTrackTmean() %15f\n", GetTrackTmean());
1421     printf("\n");
1422   }
1423   kprint++;
1424
1425   return;
1426 }
1427
1428
1429 //===================================================================================
1430 //                                 dEdx Parameterization
1431 //===================================================================================
1432
1433 Double_t AliTRDdEdxUtils::Q0MeanTRDpp(const Double_t bg)
1434 {
1435   //
1436   //truncated Mean Q_{xx} in TRD
1437   //
1438
1439   Double_t par[8];
1440   //03132012161150
1441   //opt: ppQ0
1442 par[0]=   2.397001e-01;
1443 par[1]=   1.334697e+00;
1444 par[2]=   6.967470e+00;
1445 par[3]=   9.055289e-02;
1446 par[4]=   9.388760e+00;
1447 par[5]=   9.452742e-04;
1448 par[6]=  -1.866091e+00;
1449 par[7]=   1.403545e+00;
1450
1451   ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see2.log:hhtype2Q0b2c2 scale        0.428543 at ltbg        0.650000
1452   return   0.428543*MeandEdxTR(&bg, par);
1453 }
1454
1455 Double_t AliTRDdEdxUtils::Q0MeanTRDPbPb(const Double_t bg)
1456 {
1457   //
1458   //truncated Mean Q_{xx} in TRD
1459   //
1460
1461   Double_t par[8];
1462
1463   //03132012161259
1464   //opt: PbPbQ0
1465 par[0]=   1.844912e-01;
1466 par[1]=   2.509702e+00;
1467 par[2]=   6.744031e+00;
1468 par[3]=   7.355123e-02;
1469 par[4]=   1.166023e+01;
1470 par[5]=   1.736186e-04;
1471 par[6]=  -1.716063e+00;
1472 par[7]=   1.611366e+00;
1473
1474   ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see4.log:hhtype4Q0b2c2 scale        0.460994 at ltbg        0.650000
1475   return   0.460994*MeandEdxTR(&bg, par);
1476 }
1477
1478 Double_t AliTRDdEdxUtils::Q1MeanTRDpp(const Double_t bg)
1479 {
1480   //
1481   //truncated Mean 1/(1/Q)_{xx} in TRD
1482   //
1483
1484   Double_t par[8];
1485
1486   //So 4. Mär 13:30:51 CET 2012
1487   //opt: trdppQ1
1488   par[0]=   2.434646e-01;
1489   par[1]=   1.400211e+00;
1490   par[2]=   6.937471e+00;
1491   par[3]=   7.758118e-02;
1492   par[4]=   1.097372e+01;
1493   par[5]=   4.297518e-04;
1494   par[6]=  -1.806266e+00;
1495   par[7]=   1.543811e+00;
1496
1497   //hhtype2Q1b2c2 scale        0.418629 at ltbg        0.650000
1498
1499   return  0.418629*MeandEdxTR(&bg, par);
1500 }
1501
1502 Double_t AliTRDdEdxUtils::Q1MeanTRDPbPb(const Double_t bg)
1503 {
1504   //
1505   //truncated Mean 1/(1/Q)_{xx} in TRD
1506   //
1507
1508   Double_t par[8];
1509
1510   //So 4. Mär 13:30:52 CET 2012
1511   //opt: trdPbPbQ1
1512   par[0]=   2.193660e-01;
1513   par[1]=   2.051864e+00;
1514   par[2]=   6.825112e+00;
1515   par[3]=   6.151693e-02;
1516   par[4]=   1.390343e+01;
1517   par[5]=   6.010032e-05;
1518   par[6]=  -1.676324e+00;
1519   par[7]=   1.838873e+00;
1520
1521   //hhtype4Q1b2c2 scale        0.457988 at ltbg        0.650000
1522
1523   return  0.457988*MeandEdxTR(&bg, par);
1524 }
1525
1526 Double_t AliTRDdEdxUtils::QMeanTPC(const Double_t bg)
1527 {
1528   //
1529   //bethe bloch in TPC
1530   //
1531
1532   Double_t par[5];
1533   //Mi 15. Feb 14:48:05 CET 2012
1534   //train_2012-02-13_1214.12001, tpcsig
1535   par[0]=       4.401269;
1536   par[1]=       9.725370;
1537   par[2]=       0.000178;
1538   par[3]=       1.904962;
1539   par[4]=       1.426576;
1540
1541   return MeandEdx(&bg, par);
1542 }
1543
1544 Double_t AliTRDdEdxUtils::MeandEdxTR(const Double_t * xx, const Double_t * pin)
1545 {
1546   //
1547   //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
1548   //npar = 8 = 3+5
1549   //
1550   Double_t ptr[4]={0};
1551   for(int ii=0; ii<3; ii++){
1552     ptr[ii+1]=pin[ii];
1553   }
1554   return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
1555 }
1556
1557 Double_t AliTRDdEdxUtils::MeanTR(const Double_t * xx, const Double_t * par)
1558 {
1559   //
1560   //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
1561   //npar = 4
1562   //
1563
1564   const Double_t bg = xx[0];
1565   const Double_t gamma = sqrt(1+bg*bg);
1566
1567   const Double_t p0 = TMath::Abs(par[1]);
1568   const Double_t p1 = TMath::Abs(par[2]);
1569   const Double_t p2 = TMath::Abs(par[3]);
1570
1571   const Double_t zz = TMath::Log(gamma);
1572   const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
1573
1574   return par[0]+tryield;
1575 }
1576
1577 Double_t AliTRDdEdxUtils::MeandEdx(const Double_t * xx, const Double_t * par)
1578 {
1579   //
1580   //ALEPH parametrization for dEdx
1581   //npar = 5
1582   //
1583
1584   const Double_t bg = xx[0];
1585   const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
1586
1587   const Double_t p0 = TMath::Abs(par[0]);
1588   const Double_t p1 = TMath::Abs(par[1]);
1589   const Double_t p2 = TMath::Abs(par[2]);
1590   const Double_t p3 = TMath::Abs(par[3]);
1591   const Double_t p4 = TMath::Abs(par[4]);
1592
1593   const Double_t aa = TMath::Power(beta, p3);
1594   const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
1595
1596   //printf("test----- %f %f -- %f %f %f %f %f --- %f %f %f\n", bg, beta, p0, p1, p2, p3, p4, p0/aa, aa, bb);
1597
1598   return (p1-aa-bb)*p0/aa;
1599 }
1600
1601 Double_t AliTRDdEdxUtils::ToLogx(FFunc func, const Double_t * xx, const Double_t * par)
1602 {
1603   //
1604   //f(x)-> f(y) with y=log10(x)
1605   //
1606   const Double_t x2[]={TMath::Power(10, xx[0])};
1607   return func(x2, par);
1608 }
1609
1610 //===================================================================================
1611 //                                Detector, Data and Control Constant
1612 //===================================================================================
1613 Int_t  AliTRDdEdxUtils::ToDetector(const Int_t gtb)
1614 {
1615   //
1616   //gtb = det*Ntb+itb
1617   //
1618   return gtb/AliTRDseedV1::kNtb;
1619 }
1620
1621 Int_t  AliTRDdEdxUtils::ToTimeBin(const Int_t gtb)
1622 {
1623   //
1624   //gtb = det*Ntb+itb
1625   //
1626   return gtb%AliTRDseedV1::kNtb;
1627 }
1628
1629 Int_t  AliTRDdEdxUtils::ToSector(const Int_t gtb)
1630 {
1631   //
1632   //return sector
1633   //
1634   return AliTRDgeometry::GetSector(ToDetector(gtb));
1635 }
1636
1637 Int_t  AliTRDdEdxUtils::ToStack(const Int_t gtb)
1638 {
1639   //
1640   //return stack
1641   //
1642   return AliTRDgeometry::GetStack(ToDetector(gtb));
1643 }
1644
1645 Int_t  AliTRDdEdxUtils::ToLayer(const Int_t gtb)
1646 {
1647   //
1648   //return layer
1649   //
1650   return AliTRDgeometry::GetLayer(ToDetector(gtb));
1651 }
1652
1653 TString AliTRDdEdxUtils::GetRunType(const Int_t run)
1654 {
1655   //
1656   //return run type
1657   //
1658
1659   TString type;
1660   if(run>=121527 && run<= 126460)//LHC10d
1661     type="pp2010LHC10d";
1662   else if(run>=126461 && run<= 130930)//LHC10e
1663     type="pp2010LHC10e";
1664   else if(run>=136782 && run <= 139846)//LHC10h
1665     type="PbPb2010LHC10h";
1666   else if(run>= 143224 && run<= 143237)//2011Feb
1667     type="cosmic2011Feb";
1668   else if(run>= 150587 && run<= 154930){
1669     type="cosmic2011MayJun";
1670
1671     TString runstr(Form("%d",run));
1672     const TString listrun1kg("154601 154602 154629 154634 154636 154639 154643");
1673     if(listrun1kg.Contains(runstr)){
1674       type+="1kG";
1675     }
1676     else{
1677       type+="5kG";
1678     }
1679   }
1680   else{
1681     type="unknown";
1682   }
1683
1684   type.ToUpper();
1685   return type;
1686 }
1687
1688 void AliTRDdEdxUtils::PrintControl()
1689 {
1690   //
1691   //print out control variable
1692   //
1693   printf("\nAliTRDdEdxUtils::PrintControl Q0Frac %.1f, Q1Frac %.1f, TimeBinCountCut %.2f, CalibTPCnclsCut %d, IsExBOn %d\n\n", Q0Frac(), Q1Frac(), TimeBinCountCut(), CalibTPCnclsCut(), IsExBOn());
1694 }
1695 //===================================================================================
1696 //===================================================================================