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       array->ls();
 681       return kFALSE;
 682     }
 683     THnSparseD *hi = (THnSparseD*)tmph->Clone();
 684     fgHistPHQ->AddAt(hi, iter);
 685   }
 686
 687   return kTRUE;
 688 }
 689
 690 void AliTRDdEdxUtils::FillCalibHist(const Int_t ncls, const TVectorD *arrayQ, const TVectorD *arrayX, THnSparseD * hcalib, const Double_t scale)
 691 {
 692   //
 693   //fill calibration hist
 694   //
 695   if(!hcalib){printf("AliTRDdEdxUtils::FillCalibHist errro hcalib null!!\n"); exit(1);}
 696
 697   for(Int_t ii=0; ii<ncls; ii++){
 698     const Double_t dq = (*arrayQ)[ii];
 699     const Double_t xx = (*arrayX)[ii];
 700
 701     const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
 702     const Double_t xmin = hcalib->GetAxis(0)->GetXmin();
 703     const Double_t xmax = hcalib->GetAxis(0)->GetXmax();
 704
 705     if(xx>=xmax || xx<xmin){
 706       printf("AliTRDdEdxUtils::FillCalibHist error x overflow or underflow! %s %15f %15f %15f\n", hcalib->GetName(),  xx, xmin, xmax); exit(1);
 707     }
 708
 709     const Double_t var[]={xx, TMath::Min(dq, qmax)/scale};
 710     hcalib->Fill(var);
 711   }
 712 }
 713
 714 void AliTRDdEdxUtils::FillCalibHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale)
 715 {
 716   //
 717   //get cluster Q and fill calib hist, if kinvq = kTRUE, 1/Q is filled
 718   //
 719
 720   THnSparseD * hcalib = AliTRDdEdxUtils::GetHistPHQ(kinvq, mag, charge);
 721
 722   TVectorD arrayQ(200), arrayX(200);
 723   const Int_t ncls = AliTRDdEdxUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, trdv1);
 724   FillCalibHist(ncls, &arrayQ, &arrayX, hcalib, kinvq ? 1/scale : scale);
 725
 726   static Int_t kprint = 100;
 727   if(kprint<0){
 728     printf("\nAliTRDdEdxUtils::FillCalibHist summary: \n");
 729     printf("\nkinvq= %d;\n", kinvq);
 730     for(Int_t iq=0; iq<ncls; iq++){
 731       printf("arrayX[%3d] = %15.0f; arrayQ[%3d] =  %15f;\n", iq, arrayX[iq], iq, arrayQ[iq]);
 732     }
 733     printf("\n");
 734   }
 735   kprint++;
 736 }
 737
 738 Int_t AliTRDdEdxUtils::ApplyCalib(const Int_t nc0, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
 739 {
 740   //
 741   //apply calibration on arrayQ
 742   //
 743   if(!cobj){ printf("AliTRDdEdxUtils::ApplyCalib error gain array null!!\n"); exit(1);}
 744
 745   TVectorD tmpq(arrayQ->GetNrows());
 746   TVectorD tmpx(arrayX->GetNrows());
 747   Int_t ncls = 0;
 748
 749   const TVectorD * gain = (TVectorD*) cobj->At(0);
 750   for(Int_t ii=0; ii<nc0; ii++){
 751     const Double_t dq = (*arrayQ)[ii];
 752     const Int_t xx = (Int_t)(*arrayX)[ii];
 753     const Double_t gg = (*gain)[xx];
 754
 755     if(gg<EPSILON){
 756       continue;
 757     }
 758
 759     tmpq[ncls] = dq*gg;
 760     tmpx[ncls] = xx;
 761     ncls++;
 762   }
 763
 764   (*arrayQ)=tmpq;
 765   (*arrayX)=tmpx;
 766
 767   return ncls;
 768 }
 769
 770 void AliTRDdEdxUtils::GetPHCountMeanRMS(const TH1D *hnor, TH1D *&hmean)
 771 {
 772   //
 773   //calculate from the ph calib hist the (mean-3sigma) ph-count in the chamber, save in the TH1D output
 774   //
 775   const Int_t ndet = 540;
 776   TObjArray *obj=new TObjArray(ndet);
 777   obj->SetOwner();
 778   for(Int_t ii=0; ii<ndet; ii++){
 779     obj->Add(new TVectorD(AliTRDseedV1::kNtb));
 780   }
 781
 782   //ibin = binlowedge of bin(ibin+1) = the number fills this bin
 783   for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
 784     const Double_t stat = hnor->GetBinContent(ibin+1);
 785     if(stat<EPSILON){
 786       continue;
 787     }
 788
 789     const Int_t idet = ToDetector(ibin);
 790     const Int_t itb  = ToTimeBin(ibin);
 791     TVectorD *vec=(TVectorD *)obj->At(idet);
 792     (*vec)[itb] = stat;
 793   }
 794
 795   hmean = new TH1D(Form("%sdetmean", hnor->GetName()), "", hnor->GetNbinsX(), hnor->GetXaxis()->GetXmin(), hnor->GetXaxis()->GetXmax());
 796   for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
 797     const Int_t idet = ToDetector(ibin);
 798     const TVectorD *vec=(TVectorD *)obj->At(idet);
 799
 800     Int_t nonzero = 0;
 801     for(Int_t ii=0; ii<vec->GetNrows(); ii++){
 802       if((*vec)[ii]>EPSILON){
 803         nonzero++;
 804       }
 805     }
 806
 807     Double_t mean = 0;
 808     const Double_t lowfrac = 0.6;
 809     //if there are too many 0's, reject this chamber by settig mean=rms=0
 810     if(nonzero> (AliTRDseedV1::kNtb*(1-lowfrac)) ){
 811       //only highest (1-lowfrac)*31 timebins are used to estimate the mean and rms! important! otherwise the 0' will make rms very large!
 812       mean = TruncatedMean(AliTRDseedV1::kNtb, vec->GetMatrixArray(), lowfrac, 1);
 813     }
 814
 815     hmean->SetBinContent(ibin+1, mean);
 816   }
 817
 818   delete obj;
 819 }
 820
 821 void AliTRDdEdxUtils::CalibOutput(const TList *lin, Int_t run)
 822 {
 823   //
 824   //produce calibration objects
 825   //
 826
 827   TString objnames("TRDCalibHistGain TRDCalibHistT0 TRDCalibHistVd ");
 828   for(Int_t iter=0; iter<8; iter++){
 829     objnames+= GetPHQName(0, iter)+" ";
 830   }
 831
 832   TList * lout = new TList;
 833   lout->SetOwner();
 834
 835   TTreeSRedirector *calibStream = new TTreeSRedirector(Form("TRDCalibStream_%010d.root", run));
 836
 837   const Int_t nh=lin->GetEntries();
 838   for(Int_t ii=0; ii<nh; ii++){
 839     const THnSparseD *hh=(THnSparseD*)lin->At(ii);
 840     const TString hname = hh->GetName();
 841     if(!objnames.Contains(hname))
 842       continue;
 843
 844     TObjArray * cobj0 = GetCalibObj(hh, run, lout, calibStream);
 845     lout->Add(cobj0);
 846   }
 847
 848   //lout->ls();
 849
 850   //=============================================================
 851   //=============================================================
 852
 853   TFile *fout=new TFile(Form("TRDCalibObj_%010d.root", run),"recreate");
 854   fout->cd();
 855   const Int_t nout=lout->GetEntries();
 856   for(Int_t ii=0; ii<nout; ii++){
 857     const TString oname = lout->At(ii)->GetName();
 858     if(oname.Contains("Obj")){
 859       TObjArray * cobj = (TObjArray*) lout->At(ii);
 860       cobj->Write(oname, TObjArray::kSingleKey);
 861     }
 862   }
 863   fout->Save();
 864   fout->Close();
 865   delete fout;
 866
 867   fout=new TFile(Form("TRDCalibList_%010d.root", run),"recreate");
 868   fout->cd();
 869   lin->Write();
 870   lout->Write();
 871   fout->Save();
 872   fout->Close();
 873   delete fout;
 874
 875   delete calibStream;
 876
 877   /*
 878     http://root.cern.ch/root/html/TH1.html
 879     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:
 880
 881     h->SetDirectory(0);          for the current histogram h
 882     TH1::AddDirectory(kFALSE);   sets a global switch disabling the reference
 883
 884     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.
 885   */
 886   delete lout;
 887 }
 888
 889 TObjArray* AliTRDdEdxUtils::GetCalibObj(const THnSparseD *hh, Int_t run, TList *lout, TTreeSRedirector *calibStream)
 890 {
 891   //
 892   //produce calibration objects
 893   //
 894
 895   const TString hname = hh->GetName();
 896   const Bool_t kinvq = TString(hname(hname.First('Q')+1,1)).Atoi()&4;
 897
 898   //----------------------------------------
 899   //               Define nbin, tag, cobj0
 900   //----------------------------------------
 901   Int_t nbin =-999;
 902   if(hname.Contains("Gain") || hname.Contains("T0") || hname.Contains("Vd")){
 903     nbin = NTRDchamber();
 904   }
 905   else if(hname.Contains("PHQ")){
 906     nbin = NTRDtimebin();
 907   }
 908   else{
 909     printf("AliTRDdEdxUtils::GetCalibObj error wrong hname!! %s\n", hname.Data()); exit(1);
 910   }
 911
 912   TString tag(hname);
 913   tag.ReplaceAll("Hist","Obj");
 914
 915   TObjArray * cobj0 = new TObjArray(1);
 916   cobj0->SetOwner();
 917   cobj0->SetName(tag);
 918   cobj0->Add(new TVectorD(nbin));
 919
 920   //----------------------------------------
 921   //               Define lowFrac, highFrac
 922   //----------------------------------------
 923   Double_t lowFrac = -999, highFrac = -999;
 924   if(hname.Contains("Gain") || (hname.Contains("PHQ") && !kinvq) ){
 925     lowFrac = 0.01; highFrac = Q0Frac();
 926   }
 927   else if(hname.Contains("PHQ") && kinvq){
 928     lowFrac = Q1Frac(); highFrac = 0.99;
 929   }
 930   else{
 931     lowFrac = 0.01;
 932     highFrac = 0.99;
 933   }
 934
 935   //----------------------------------------
 936   //               Get analysis result
 937   //----------------------------------------
 938   TH1::AddDirectory(kFALSE);//important!
 939   TH1D *hnor=0x0, *hmpv=0x0, *hres=0x0, *hwid=0x0, *htrdphmean = 0x0;//if(!lout), these have to be deleted
 940   TH2D *hpj = hh->Projection(1,0);
 941   FitSlicesY(hpj, hnor, hmpv, hwid, hres, 0, lowFrac, highFrac);
 942   if(hname.Contains("PHQ")){
 943     GetPHCountMeanRMS(hnor, htrdphmean);
 944     if(lout){
 945       lout->Add(htrdphmean);
 946     }
 947   }
 948   delete hpj;
 949
 950   if(lout){
 951     lout->Add(hnor);
 952     lout->Add(hmpv);
 953     lout->Add(hwid);
 954     lout->Add(hres);
 955   }
 956
 957   //----------------------------------------
 958   //               Define Counter
 959   //----------------------------------------
 960   TVectorD *countDet=0x0;
 961   TObjArray *countSSL=0x0;
 962
 963   if(hname.Contains("PHQ") && !kinvq){
 964     countDet=new TVectorD(540);
 965     countSSL=new TObjArray(90);//SectorStackLayer
 966     countSSL->SetOwner();
 967     for(Int_t ii=0; ii<90; ii++){
 968       countSSL->Add(new TVectorD(6));
 969     }
 970   }
 971
 972   //----------------------------------------
 973   //               Fill cobj0
 974   //----------------------------------------
 975
 976   //ibin = binlowedge of bin(ibin+1) = the number fills this bin
 977   for(Int_t ibin=0; ibin<nbin; ibin++){
 978     Double_t gnor = hnor->GetBinContent(ibin+1);
 979     Double_t gmpv = hmpv->GetBinContent(ibin+1);
 980     Double_t gwid = hwid->GetBinContent(ibin+1);
 981     Double_t gres = hres->GetBinContent(ibin+1);
 982
 983     //--- set additional cut by kpass
 984     Bool_t kpass = kTRUE;
 985     Double_t gtrdphmean = -999;
 986     if(htrdphmean){
 987       gtrdphmean = htrdphmean->GetBinContent(ibin+1);
 988       //chamber no statistics (e.g. too many 0's), not usual, not seen in run 143237
 989       if(gtrdphmean<EPSILON){
 990         kpass = kFALSE;
 991       }
 992       if(gnor<TimeBinCountCut()*gtrdphmean){
 993         kpass = kFALSE;
 994       }
 995     }
 996
 997     //--- set calibration constant p0
 998     Double_t p0= 0;
 999
1000     //reason for gmpv=0:
1001     //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;
1002     //2)TruncatedMean(hist) out of range (only for Q0, not Q1).
1003
1004     if(gmpv>EPSILON && kpass){
1005       if(tag.Contains("T0")){
1006         p0 = gmpv;
1007       }
1008       else{
1009         p0 = 1/gmpv;
1010       }
1011       //printf("outcalibobj%s %d %15.6e\n", tag.Data(), ibin, p0);
1012     }
1013
1014     (*( (TVectorD*)cobj0->At(0) ))[ibin] = p0;
1015
1016     //--- save optional record
1017     if(p0>EPSILON && countDet && countSSL){
1018       const Int_t idet = ToDetector(ibin);
1019       (*countDet)[idet]=1;
1020
1021       const Int_t isector = ToSector(ibin);
1022       const Int_t istack = ToStack(ibin);
1023       const Int_t ilayer = ToLayer(ibin);
1024       TVectorD * vecsectorstack = (TVectorD*)countSSL->At(istack*18+isector);
1025       (*vecsectorstack)[ilayer]=1;
1026     }
1027
1028     if(calibStream){
1029       (*calibStream)<<tag<<
1030         "run="<<run<<
1031         "p0="<<p0<<
1032
1033         "nor="<<gnor<<
1034         "mpv="<<gmpv<<
1035         "wid="<<gwid<<
1036         "res="<<gres<<
1037         "gtrdphmean="<<gtrdphmean<<
1038
1039         "ibin="<<ibin<<
1040         "\n";
1041     }
1042   }
1043
1044   //----------------------------------------
1045   //               Status Report
1046   //----------------------------------------
1047   if(countDet && countSSL){
1048     TVectorD count2Dstack(90);
1049     for(Int_t ii=0; ii<90; ii++){
1050       TVectorD * vecsectorstack = (TVectorD*)countSSL->At(ii);
1051       const Int_t nlayer = (Int_t)vecsectorstack->Sum();
1052       if(nlayer==6){
1053         count2Dstack[ii]=1;
1054       }
1055     }
1056
1057     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());
1058   }
1059
1060   //----------------------------------------
1061   //               Clean Up
1062   //----------------------------------------
1063
1064   TH1D **hhs[]={&hnor, &hmpv, &hwid, &hres, &htrdphmean};
1065   const Int_t nhh=sizeof(hhs)/sizeof(TH1D**);
1066   for(Int_t ihh=0; ihh<nhh; ihh++){
1067     if(!lout){
1068       delete (*hhs[ihh]);
1069     }
1070   }
1071
1072   delete countDet;
1073   delete countSSL;
1074
1075   //----------------------------------------
1076
1077   return cobj0;
1078 }
1079
1080 //===================================================================================
1081 //                                   dEdx calculation
1082 //===================================================================================
1083 Double_t AliTRDdEdxUtils::GetSignal(const Int_t nch, const Int_t ncls, const Double_t qq)
1084 {
1085   //
1086   //signal = nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq)
1087   //
1088   if(ncls>1e3){
1089     printf("AliTRDdEdxUtils::GetSignal error ncls>1e3! %d\n", ncls); exit(1);
1090   }
1091
1092   return nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq);
1093 }
1094
1095 Int_t AliTRDdEdxUtils::GetNch(const Double_t signal)
1096 {
1097   //
1098   //signal = nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq)
1099   //
1100   return Int_t(signal/1e6);
1101
1102 }
1103
1104 Int_t AliTRDdEdxUtils::GetNcls(const Double_t signal)
1105 {
1106   //
1107   //signal = Nch*1e6 + Ncls*1e3 + (Q>=1e3? 999 : Q)
1108   //
1109
1110   return Int_t ( (signal-GetNch(signal)*1e6)/1e3 );
1111 }
1112
1113 Double_t AliTRDdEdxUtils::GetQ(const Double_t signal)
1114 {
1115   //
1116   //signal = nch*1e6 + ncls*1e3 + (qq>=1e3? 999 : qq)
1117   //
1118
1119   return signal-GetNch(signal)*1e6 - GetNcls(signal)*1e3;
1120 }
1121
1122 Double_t AliTRDdEdxUtils::ToyCook(const Bool_t kinvq, Int_t &ncluster, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
1123 {
1124   //
1125   //template for cookdedx
1126   //
1127   if(cobj){
1128     if(arrayQ && arrayX){
1129       ncluster = ApplyCalib(ncluster, arrayQ, arrayX, cobj);
1130     }
1131     else{
1132       printf("AliTRDdEdxUtils::ToyCook arrayQ arrayX null, applycalib can not be applied!\n"); exit(1);
1133     }
1134   }
1135
1136   Double_t lowFrac =-999, highFrac = -999;
1137   if(kinvq){
1138     lowFrac = Q1Frac(); highFrac = 0.99;
1139   }
1140   else{
1141     lowFrac = 0.01; highFrac = Q0Frac();
1142   }
1143
1144   Double_t meanQ = TruncatedMean(ncluster, arrayQ->GetMatrixArray(), lowFrac, highFrac);
1145   if(kinvq){
1146     if(meanQ>EPSILON){
1147       meanQ = 1/meanQ;
1148     }
1149   }
1150
1151   return meanQ;
1152 }
1153
1154 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)
1155 {
1156   //
1157   //combine tpc and trd dedx
1158   //
1159
1160   for(Int_t iq=0; iq<tpcncls; iq++){
1161     (*coarrayQ)[iq]=(*tpcarrayQ)[iq];
1162     if(tpcarrayX && trdarrayX && coarrayX){
1163       (*coarrayX)[iq]=(*tpcarrayX)[iq];
1164     }
1165   }
1166   for(Int_t iq=0; iq<trdncls; iq++){
1167     (*coarrayQ)[tpcncls+iq]=(*trdarrayQ)[iq];
1168     if(tpcarrayX && trdarrayX && coarrayX){
1169       (*coarrayX)[tpcncls+iq]=159+(*trdarrayX)[iq];
1170     }
1171   }
1172
1173   concls=trdncls+tpcncls;
1174
1175   const Double_t coQ = ToyCook(kinvq, concls, coarrayQ, coarrayX);
1176
1177   return coQ;
1178 }
1179
1180
1181 //===================================================================================
1182 //                                   dEdx Getter and Setter
1183 //===================================================================================
1184 Double_t AliTRDdEdxUtils::GetAngularCorrection(const AliTRDseedV1 *seed)
1185 {
1186   //
1187   //return angular normalization factor
1188   //
1189
1190   return TMath::Sqrt(1+seed->GetYref(1)*seed->GetYref(1)+seed->GetZref(1)*seed->GetZref(1));
1191 }
1192
1193 Double_t AliTRDdEdxUtils::GetClusterQ(const Bool_t kinvq, const AliTRDseedV1 * seed, const Int_t itb)
1194 {
1195   //
1196   //get cluster charge
1197   //
1198   Double_t dq = 0;
1199   const AliTRDcluster *cl = 0x0;
1200
1201   cl = seed->GetClusters(itb);                    if(cl) dq+=cl->GetRawQ();
1202   cl = seed->GetClusters(itb+AliTRDseedV1::kNtb); if(cl) dq+=cl->GetRawQ();
1203
1204   dq /= GetAngularCorrection(seed);
1205
1206   dq /= 45.;
1207
1208   if(kinvq){
1209     if(dq>EPSILON){
1210       dq = 1/dq;
1211     }
1212   }
1213
1214   return dq;
1215 }
1216
1217 Int_t AliTRDdEdxUtils::GetArrayClusterQ(const Bool_t kinvq, TVectorD *arrayQ, TVectorD *arrayX, const AliTRDtrackV1 *trdtrack, Int_t timeBin0, Int_t timeBin1, Int_t tstep)
1218 {
1219   //
1220   //return nclustter
1221   //(if kinvq, return 1/q array), size of array must be larger than 31*6
1222   //
1223   if(!arrayQ || arrayQ->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
1224     printf("AliTRDdEdxUtils::GetArrayClusterQ error arrayQ null or size too small! %d\n", arrayQ? arrayQ->GetNrows() : -999); exit(1);
1225   }
1226   if(!arrayX && arrayX->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
1227     printf("AliTRDdEdxUtils::GetArrayClusterQ error arrayX null or size too small! %d\n", arrayX? arrayX->GetNrows() : -999); exit(1);
1228   }
1229
1230   const Int_t mintb = 0;
1231   const Int_t maxtb = AliTRDseedV1::kNtb-1;
1232   if(timeBin0<mintb) timeBin0=mintb;
1233   if(timeBin1>maxtb) timeBin1=maxtb;
1234   if(tstep<=0) tstep=1;
1235
1236   //============
1237   Int_t tbN=0;
1238   Double_t tbQ[200];
1239   Int_t tbBin[200];
1240
1241   for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){
1242     const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
1243     if(!seed)
1244       continue;
1245
1246     const Int_t det = seed->GetDetector();
1247
1248     for(Int_t itb=timeBin0; itb<=timeBin1; itb+=tstep){
1249       const Double_t dq = GetClusterQ(kinvq, seed, itb);
1250       if(dq<EPSILON)
1251         continue;
1252
1253       const Int_t gtb = det * AliTRDseedV1::kNtb + itb;
1254
1255       tbQ[tbN]=dq;
1256       tbBin[tbN]=gtb;
1257       tbN++;
1258     }
1259   }
1260
1261   Int_t ncls = 0;
1262   for(Int_t iq=0; iq<tbN; iq++){
1263     if(tbQ[iq]<EPSILON)
1264       continue;
1265
1266     (*arrayQ)[ncls] = tbQ[iq];
1267     (*arrayX)[ncls] = tbBin[iq];
1268
1269     ncls++;
1270   }
1271
1272   Int_t kprint = 100;
1273   if(kprint<0){
1274     printf("\nAliTRDdEdxUtils::GetArrayClusterQ raw cluster-Q\n");
1275     for(Int_t iq=0; iq<ncls; iq++){
1276       const Int_t ichamber =  ToLayer((*arrayX)[iq]);
1277       const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
1278       if(!seed){
1279         printf("error seed null!!\n"); exit(1);
1280       }
1281       const Double_t rawq =  (*arrayQ)[iq] * 45. * GetAngularCorrection(seed);
1282       printf("esdid=%d; chamber=%d; timebin=%d; rawq= %.3f; myq[%d]= %e;\n", trdtrack->GetESDid(), ichamber, ToTimeBin((*arrayX)[iq]), rawq, iq, (*arrayQ)[iq]);
1283     }
1284     printf("\n");
1285   }
1286
1287   return ncls;
1288 }
1289
1290 Int_t AliTRDdEdxUtils::UpdateArrayX(const Int_t ncls, TVectorD* arrayX)
1291 {
1292   //
1293   //arrayX det*Ntb+itb -> itb
1294   //
1295
1296   TVectorD countChamber(6);
1297   for(Int_t ii = 0; ii<ncls; ii++){
1298     const Int_t xx = (Int_t)(*arrayX)[ii];
1299     const Int_t idet = ToDetector(xx);
1300
1301     const Double_t ich = AliTRDgeometry::GetLayer(idet);
1302     const Double_t itb = ToTimeBin(xx);
1303     (*arrayX)[ii] = ich*AliTRDseedV1::kNtb+itb;
1304
1305     countChamber[ich] = 1;
1306   }
1307
1308   const Double_t nch = countChamber.Sum();
1309   return (Int_t) nch;
1310 }
1311
1312 void AliTRDdEdxUtils::SetChamberQT(const AliTRDtrackV1 *trdtrack, const Int_t kcalib, THnSparseD * hgain, THnSparseD * ht0, THnSparseD * hvd)
1313 {
1314   //
1315   //CookdEdx at TRD track level, use chamber info, related calibrations: chamber-gain; T0, Vd based on raw PH distribution
1316   //
1317
1318   static Int_t kprint = 100;
1319
1320   fgNchamber = 0;
1321   for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){
1322     //initialize output, default values: 0, so that summation and weighting will automatically discard default quantities
1323     fgChamberQ[ichamber] = fgChamberTmean[ichamber] = 0;
1324
1325     const AliTRDseedV1 *seed = trdtrack->GetTracklet(ichamber);
1326     if (!seed)
1327       continue;
1328
1329     const Int_t idet = seed->GetDetector();
1330
1331     //-------------------------------------------------------------------------
1332
1333     Double_t qsum = 0, qtsum = 0, w2sum = 0;
1334     for(Int_t itb=0; itb<AliTRDseedV1::kNtb; itb++){
1335       const Double_t dq = GetClusterQ(0, seed, itb);
1336       if(dq<EPSILON)
1337         continue;
1338
1339       qsum += dq;
1340       qtsum += dq*itb;
1341       w2sum += dq*itb*itb;
1342     }
1343     if(qsum<EPSILON)
1344       continue;
1345
1346     //-------------------------------------------------------------------------
1347
1348     Double_t tbm, tbr = 0;
1349     tbm = GetMeanRMS(qsum, qtsum, w2sum, &tbr);
1350
1351     qsum /= 1.25e3/45.;
1352
1353     if(hgain){
1354       const Double_t var[]={idet, qsum};
1355       hgain->Fill(var);
1356     }
1357     if(ht0){
1358       const Double_t var[]={idet, tbm};
1359       ht0->Fill(var);
1360     }
1361     if(hvd){
1362       const Double_t var[]={idet, tbr};
1363       hvd->Fill(var);
1364     }
1365
1366     Double_t gain = 1, t0 = 0, vd = 1;
1367     if(kcalib){
1368       if(!fgObjGain) {printf("AliTRDdEdxUtils::SetChamberQT error Gain array null!!\n"); exit(1);}
1369       if(!  fgObjT0) {printf("AliTRDdEdxUtils::SetChamberQT error T0   array null!!\n"); exit(1);}
1370       if(!  fgObjVd) {printf("AliTRDdEdxUtils::SetChamberQT error Vd   array null!!\n"); exit(1);}
1371
1372       const TVectorD * gainvec = (TVectorD*) fgObjGain->At(0); gain = (*gainvec)[idet];
1373       const TVectorD *   t0vec = (TVectorD*)   fgObjT0->At(0);   t0 = (*  t0vec)[idet];
1374       const TVectorD *   vdvec = (TVectorD*)   fgObjVd->At(0);   vd = (*  vdvec)[idet];
1375     }
1376     if(kprint<0){
1377       printf("\nAliTRDdEdxUtils::CookdEdxV2\n");
1378       printf("idet = %d;\n", idet);
1379       printf("gain = %15f; t0 = %15f; vd = %15f;\n", gain, t0, vd);
1380       printf("\n");
1381     }
1382
1383     qsum *= gain;
1384     tbm = (tbm-t0)*vd;
1385
1386     if(qsum<EPSILON)
1387       continue;
1388
1389     //-------------------------------------------------------------------------
1390
1391     //should have non-zero value, initialized with default 0 (except for calibrated tbm, may be very close to 0)
1392     fgChamberQ[ichamber]  = qsum;
1393     fgChamberTmean[ichamber] = tbm;
1394     fgNchamber++;
1395   }
1396
1397   if(kprint<0){
1398     printf("\nAliTRDdEdxUtils::CookdEdxV2 summary:\n");
1399
1400     printf("\nfgNchamber = %d\n", fgNchamber);
1401     for(Int_t ich=0; ich<AliTRDtrackV1::kNplane; ich++){
1402       printf("fgChamberTmean[%d] = %15f; fgChamberQ[%d] = %15f;\n", ich, fgChamberTmean[ich], ich, fgChamberQ[ich]);
1403     }
1404   }
1405
1406   fgTrackTmean = -999;
1407   if(fgNchamber){
1408     fgTrackTmean = 0;
1409     for(Int_t ich=0; ich<AliTRDtrackV1::kNplane; ich++){
1410       fgTrackTmean += fgChamberTmean[ich];
1411     }
1412     fgTrackTmean /= fgNchamber;
1413   }
1414
1415   if(kprint<0){
1416     printf("\nAliTRDdEdxUtils::CookdEdxV2\n");
1417     printf("GetTrackTmean() %15f\n", GetTrackTmean());
1418     printf("\n");
1419   }
1420   kprint++;
1421
1422   return;
1423 }
1424
1425
1426 //===================================================================================
1427 //                                 dEdx Parameterization
1428 //===================================================================================
1429
1430 Double_t AliTRDdEdxUtils::Q0MeanTRDpp(const Double_t bg)
1431 {
1432   //
1433   //truncated Mean Q_{xx} in TRD
1434   //
1435
1436   Double_t par[8];
1437   //03132012161150
1438   //opt: ppQ0
1439 par[0]=   2.397001e-01;
1440 par[1]=   1.334697e+00;
1441 par[2]=   6.967470e+00;
1442 par[3]=   9.055289e-02;
1443 par[4]=   9.388760e+00;
1444 par[5]=   9.452742e-04;
1445 par[6]=  -1.866091e+00;
1446 par[7]=   1.403545e+00;
1447
1448   ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see2.log:hhtype2Q0b2c2 scale        0.428543 at ltbg        0.650000
1449   return   0.428543*MeandEdxTR(&bg, par);
1450 }
1451
1452 Double_t AliTRDdEdxUtils::Q0MeanTRDPbPb(const Double_t bg)
1453 {
1454   //
1455   //truncated Mean Q_{xx} in TRD
1456   //
1457
1458   Double_t par[8];
1459
1460   //03132012161259
1461   //opt: PbPbQ0
1462 par[0]=   1.844912e-01;
1463 par[1]=   2.509702e+00;
1464 par[2]=   6.744031e+00;
1465 par[3]=   7.355123e-02;
1466 par[4]=   1.166023e+01;
1467 par[5]=   1.736186e-04;
1468 par[6]=  -1.716063e+00;
1469 par[7]=   1.611366e+00;
1470
1471   ///u/xlu/.task/CommondEdx/myAnaData/Optimum/check11/Merged/LHC10e_plot/Separation/see4.log:hhtype4Q0b2c2 scale        0.460994 at ltbg        0.650000
1472   return   0.460994*MeandEdxTR(&bg, par);
1473 }
1474
1475 Double_t AliTRDdEdxUtils::Q1MeanTRDpp(const Double_t bg)
1476 {
1477   //
1478   //truncated Mean 1/(1/Q)_{xx} in TRD
1479   //
1480
1481   Double_t par[8];
1482
1483   //So 4. Mär 13:30:51 CET 2012
1484   //opt: trdppQ1
1485   par[0]=   2.434646e-01;
1486   par[1]=   1.400211e+00;
1487   par[2]=   6.937471e+00;
1488   par[3]=   7.758118e-02;
1489   par[4]=   1.097372e+01;
1490   par[5]=   4.297518e-04;
1491   par[6]=  -1.806266e+00;
1492   par[7]=   1.543811e+00;
1493
1494   //hhtype2Q1b2c2 scale        0.418629 at ltbg        0.650000
1495
1496   return  0.418629*MeandEdxTR(&bg, par);
1497 }
1498
1499 Double_t AliTRDdEdxUtils::Q1MeanTRDPbPb(const Double_t bg)
1500 {
1501   //
1502   //truncated Mean 1/(1/Q)_{xx} in TRD
1503   //
1504
1505   Double_t par[8];
1506
1507   //So 4. Mär 13:30:52 CET 2012
1508   //opt: trdPbPbQ1
1509   par[0]=   2.193660e-01;
1510   par[1]=   2.051864e+00;
1511   par[2]=   6.825112e+00;
1512   par[3]=   6.151693e-02;
1513   par[4]=   1.390343e+01;
1514   par[5]=   6.010032e-05;
1515   par[6]=  -1.676324e+00;
1516   par[7]=   1.838873e+00;
1517
1518   //hhtype4Q1b2c2 scale        0.457988 at ltbg        0.650000
1519
1520   return  0.457988*MeandEdxTR(&bg, par);
1521 }
1522
1523 Double_t AliTRDdEdxUtils::QMeanTPC(const Double_t bg)
1524 {
1525   //
1526   //bethe bloch in TPC
1527   //
1528
1529   Double_t par[5];
1530   //Mi 15. Feb 14:48:05 CET 2012
1531   //train_2012-02-13_1214.12001, tpcsig
1532   par[0]=       4.401269;
1533   par[1]=       9.725370;
1534   par[2]=       0.000178;
1535   par[3]=       1.904962;
1536   par[4]=       1.426576;
1537
1538   return MeandEdx(&bg, par);
1539 }
1540
1541 Double_t AliTRDdEdxUtils::MeandEdxTR(const Double_t * xx, const Double_t * pin)
1542 {
1543   //
1544   //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
1545   //npar = 8 = 3+5
1546   //
1547   Double_t ptr[4]={0};
1548   for(int ii=0; ii<3; ii++){
1549     ptr[ii+1]=pin[ii];
1550   }
1551   return MeanTR(xx,ptr) + MeandEdx(xx,&(pin[3]));
1552 }
1553
1554 Double_t AliTRDdEdxUtils::MeanTR(const Double_t * xx, const Double_t * par)
1555 {
1556   //
1557   //ALEPH+LOGISTIC parametrization for dEdx+TR, in unit of MIP
1558   //npar = 4
1559   //
1560
1561   const Double_t bg = xx[0];
1562   const Double_t gamma = sqrt(1+bg*bg);
1563
1564   const Double_t p0 = TMath::Abs(par[1]);
1565   const Double_t p1 = TMath::Abs(par[2]);
1566   const Double_t p2 = TMath::Abs(par[3]);
1567
1568   const Double_t zz = TMath::Log(gamma);
1569   const Double_t tryield = p0/( 1 + TMath::Exp(-p1*(zz-p2)) );
1570
1571   return par[0]+tryield;
1572 }
1573
1574 Double_t AliTRDdEdxUtils::MeandEdx(const Double_t * xx, const Double_t * par)
1575 {
1576   //
1577   //ALEPH parametrization for dEdx
1578   //npar = 5
1579   //
1580
1581   const Double_t bg = xx[0];
1582   const Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
1583
1584   const Double_t p0 = TMath::Abs(par[0]);
1585   const Double_t p1 = TMath::Abs(par[1]);
1586   const Double_t p2 = TMath::Abs(par[2]);
1587   const Double_t p3 = TMath::Abs(par[3]);
1588   const Double_t p4 = TMath::Abs(par[4]);
1589
1590   const Double_t aa = TMath::Power(beta, p3);
1591   const Double_t bb = TMath::Log( p2 + TMath::Power(1./bg, p4) );
1592
1593   //printf("test----- %f %f -- %f %f %f %f %f --- %f %f %f\n", bg, beta, p0, p1, p2, p3, p4, p0/aa, aa, bb);
1594
1595   return (p1-aa-bb)*p0/aa;
1596 }
1597
1598 Double_t AliTRDdEdxUtils::ToLogx(FFunc func, const Double_t * xx, const Double_t * par)
1599 {
1600   //
1601   //f(x)-> f(y) with y=log10(x)
1602   //
1603   const Double_t x2[]={TMath::Power(10, xx[0])};
1604   return func(x2, par);
1605 }
1606
1607 //===================================================================================
1608 //                                Detector, Data and Control Constant
1609 //===================================================================================
1610 Int_t  AliTRDdEdxUtils::ToDetector(const Int_t gtb)
1611 {
1612   //
1613   //gtb = det*Ntb+itb
1614   //
1615   return gtb/AliTRDseedV1::kNtb;
1616 }
1617
1618 Int_t  AliTRDdEdxUtils::ToTimeBin(const Int_t gtb)
1619 {
1620   //
1621   //gtb = det*Ntb+itb
1622   //
1623   return gtb%AliTRDseedV1::kNtb;
1624 }
1625
1626 Int_t  AliTRDdEdxUtils::ToSector(const Int_t gtb)
1627 {
1628   //
1629   //return sector
1630   //
1631   return AliTRDgeometry::GetSector(ToDetector(gtb));
1632 }
1633
1634 Int_t  AliTRDdEdxUtils::ToStack(const Int_t gtb)
1635 {
1636   //
1637   //return stack
1638   //
1639   return AliTRDgeometry::GetStack(ToDetector(gtb));
1640 }
1641
1642 Int_t  AliTRDdEdxUtils::ToLayer(const Int_t gtb)
1643 {
1644   //
1645   //return layer
1646   //
1647   return AliTRDgeometry::GetLayer(ToDetector(gtb));
1648 }
1649
1650 TString AliTRDdEdxUtils::GetRunType(const Int_t run)
1651 {
1652   //
1653   //return run type
1654   //
1655
1656   TString type;
1657   if(run>=121527 && run<= 126460)//LHC10d
1658     type="pp2010LHC10d";
1659   else if(run>=126461 && run<= 130930)//LHC10e
1660     type="pp2010LHC10e";
1661   else if(run>=136782 && run <= 139846)//LHC10h
1662     type="PbPb2010LHC10h";
1663   else if(run>= 143224 && run<= 143237)//2011Feb
1664     type="cosmic2011Feb";
1665   else if(run>= 150587 && run<= 154930){
1666     type="cosmic2011MayJun";
1667
1668     TString runstr(Form("%d",run));
1669     const TString listrun1kg("154601 154602 154629 154634 154636 154639 154643");
1670     if(listrun1kg.Contains(runstr)){
1671       type+="1kG";
1672     }
1673     else{
1674       type+="5kG";
1675     }
1676   }
1677   else{
1678     type="unknown";
1679   }
1680
1681   type.ToUpper();
1682   return type;
1683 }
1684
1685 void AliTRDdEdxUtils::PrintControl()
1686 {
1687   //
1688   //print out control variable
1689   //
1690   printf("\nAliTRDdEdxUtils::PrintControl Q0Frac %.1f, Q1Frac %.1f, TimeBinCountCut %.2f, CalibTPCnclsCut %d, IsExBOn %d\n\n", Q0Frac(), Q1Frac(), TimeBinCountCut(), CalibTPCnclsCut(), IsExBOn());
1691 }
1692 //===================================================================================
1693 //===================================================================================