TRD/AliTRDdEdxCalibUtils.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 "AliTRDcalibDB.h"
  52 #include "AliTRDCalROC.h"
  53 #include "AliTRDtrackV1.h"
  54
  55 #include "AliTRDdEdxBaseUtils.h"
  56 #include "AliTRDdEdxReconUtils.h"
  57 #include "AliTRDdEdxCalibHistArray.h"
  58 #include "AliTRDdEdxCalibUtils.h"
  59
  60 #define EPSILON 1e-12
  61
  62 AliTRDdEdxCalibHistArray * AliTRDdEdxCalibUtils::fgHistArray = 0x0;
  63 TObjArray * AliTRDdEdxCalibUtils::fgObjArray = 0x0;
  64
  65 //============================================================
  66 //                           CalibObj
  67 //============================================================
  68 TObjArray * AliTRDdEdxCalibUtils::GetObjArray()
  69 {
  70   //
  71   //return fgObjArray, initialized if null
  72   //
  73
  74   if(!fgObjArray){
  75     fgObjArray = new TObjArray(8);
  76   }
  77
  78   return fgObjArray;
  79 }
  80
  81 TObjArray * AliTRDdEdxCalibUtils::GetObj(const Bool_t kinvq, const Double_t mag, const Int_t charge)
  82 {
  83   //
  84   //return calib obj
  85   //
  86   if(!fgObjArray){
  87     printf("AliTRDdEdxCalibUtils::GetObjArray(kinvq, mag, charge) error fgObjArray null!!\n"); exit(1);
  88   }
  89
  90   return (TObjArray*) fgObjArray->At(AliTRDdEdxCalibHistArray::GetIterator(kinvq, mag, charge));
  91 }
  92
  93 void AliTRDdEdxCalibUtils::DeleteObjArray()
  94 {
  95   //
  96   //delete calib obj
  97   //
  98   if(fgObjArray){
  99     fgObjArray->SetOwner();
 100     delete fgObjArray;
 101     fgObjArray = 0x0;
 102   }
 103 }
 104
 105 Bool_t AliTRDdEdxCalibUtils::GenerateDefaultOCDB(const TString path)
 106 {
 107   //
 108   //generate default OCDB object PHQ, do like
 109   //AliTRDdEdxCalibUtils::GenerateDefaultPHQOCDB("local://./")
 110   //
 111
 112   TObjArray * arr8 = new TObjArray(8);
 113   arr8->SetOwner();
 114
 115   for(Int_t ii=0; ii<8; ii++){
 116     TObjArray * arr1 = new TObjArray(1);
 117     arr1->SetOwner();
 118     TString objn(AliTRDdEdxCalibHistArray::GetNameAt(ii));
 119     objn.ReplaceAll("Hist","Obj");
 120     arr1->SetName(objn);
 121
 122     const Int_t nbins = AliTRDdEdxBaseUtils::NTRDtimebin();
 123     TVectorD * vec = new TVectorD(nbins);
 124     for(Int_t jj=0; jj<nbins; jj++){
 125       (*vec)[jj] = 1;
 126     }
 127     arr1->Add(vec);
 128     arr8->Add(arr1);
 129   }
 130
 131   AliCDBMetaData *metaData= new AliCDBMetaData();
 132   metaData->SetObjectClassName("TObjArray");
 133   metaData->SetResponsible("Raphaelle Bailhache and Xianguo Lu");
 134
 135   AliCDBId id1("TRD/Calib/PHQ", 0, 999999999);
 136   AliCDBStorage * gStorage = AliCDBManager::Instance()->GetStorage(path);
 137   gStorage->Put(arr8, id1, metaData);
 138
 139   delete metaData;
 140   delete arr8;
 141
 142   return kTRUE;
 143 }
 144
 145 //============================================================
 146 //                           CalibHist
 147 //============================================================
 148 void AliTRDdEdxCalibUtils::DeleteHistArray()
 149 {
 150   //
 151   //delete calib hist
 152   //
 153   delete fgHistArray;
 154   fgHistArray = 0x0;
 155 }
 156
 157 THnBase * AliTRDdEdxCalibUtils::GetHistAt(const Int_t iter)
 158 {
 159   //
 160   //
 161   //
 162   if(iter<fgHistArray->GetSize())
 163     return (THnBase*)fgHistArray->At(iter);
 164   else
 165     return 0x0;
 166 }
 167
 168 void AliTRDdEdxCalibUtils::IniHistArray(TList *list, const Bool_t kNoInv)
 169 {
 170   //
 171   //initialize calib hist, list should not own the hist, or list->Clear/delete hist should not be called
 172   //
 173
 174   delete fgHistArray;
 175   fgHistArray = new AliTRDdEdxCalibHistArray(kNoInv);
 176   list->Add(fgHistArray);
 177 }
 178
 179 Bool_t AliTRDdEdxCalibUtils::ReadHistArray(const TString filename, const TString listname)
 180 {
 181   //
 182   //used in AliTRDPreprocessorOffline
 183   //read in calib hist from file, only for PHQ
 184   //
 185
 186   //maybe already open by others... don't close
 187   TFile fcalib(filename);
 188   TObjArray * array = (TObjArray*)fcalib.Get(listname);
 189   const TString histname = AliTRDdEdxCalibHistArray::GetArrayName();
 190
 191   AliTRDdEdxCalibHistArray * tmph=0x0;
 192   if(array){
 193     tmph = (AliTRDdEdxCalibHistArray *) array->FindObject(histname);
 194   }
 195   else{
 196     tmph = (AliTRDdEdxCalibHistArray *) fcalib.Get(histname);
 197   }
 198   if(!tmph){
 199     printf("AliTRDdEdxCalibUtils::ReadCalibHist warning calib hist not found! %s %s %s\n", filename.Data(), listname.Data(), histname.Data());
 200     fcalib.ls();
 201     if(array){
 202       array->ls();
 203     }
 204     return kFALSE;
 205   }
 206
 207   delete fgHistArray;
 208   fgHistArray = new AliTRDdEdxCalibHistArray(*tmph);
 209
 210   return kTRUE;
 211 }
 212
 213 void AliTRDdEdxCalibUtils::FillHist(const Int_t ncls, const TVectorD *arrayQ, const TVectorD *arrayX, THnBase * hcalib, const Double_t scale)
 214 {
 215   //
 216   //fill calibration hist
 217   //
 218   if(!hcalib){printf("AliTRDdEdxCalibUtils::FillCalibHist errro hcalib null!!\n"); exit(1);}
 219
 220   for(Int_t ii=0; ii<ncls; ii++){
 221     const Double_t dq = (*arrayQ)[ii];
 222     const Double_t xx = (*arrayX)[ii];
 223
 224     const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
 225     const Double_t xmin = hcalib->GetAxis(0)->GetXmin();
 226     const Double_t xmax = hcalib->GetAxis(0)->GetXmax();
 227
 228     if(xx>=xmax || xx<xmin){
 229       printf("AliTRDdEdxCalibUtils::FillCalibHist error x overflow or underflow! %s %15f %15f %15f\n", hcalib->GetName(),  xx, xmin, xmax); exit(1);
 230     }
 231
 232     const Double_t var[]={xx, TMath::Min(dq, qmax)/scale};
 233     hcalib->Fill(var);
 234   }
 235 }
 236
 237 void AliTRDdEdxCalibUtils::FillHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale)
 238 {
 239   //
 240   //get cluster Q and fill calib hist, if kinvq = kTRUE, 1/Q is filled
 241   //
 242   if(!fgHistArray){
 243     printf("AliTRDdEdxCalibUtils::FillHist fgHistArray not initialized!!\n"); exit(1);
 244   }
 245
 246   THnBase * hcalib = fgHistArray->GetHist(kinvq, mag, charge);
 247
 248   TVectorD arrayQ(200), arrayX(200);
 249   const Int_t ncls = AliTRDdEdxReconUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, trdv1);
 250   FillHist(ncls, &arrayQ, &arrayX, hcalib, kinvq ? 1/scale : scale);
 251
 252   static Int_t kprint = 100;
 253   if(kprint<0){
 254     printf("\nAliTRDdEdxCalibUtils::FillHist summary: \n");
 255     printf("\nkinvq= %d;\n", kinvq);
 256     for(Int_t iq=0; iq<ncls; iq++){
 257       printf("arrayX[%3d] = %15.0f; arrayQ[%3d] =  %15f;\n", iq, arrayX[iq], iq, arrayQ[iq]);
 258     }
 259     printf("\n");
 260   }
 261   kprint++;
 262 }
 263
 264 //============================================================
 265 //
 266 //============================================================
 267
 268 Double_t AliTRDdEdxCalibUtils::GetCalibTPCscale(const Int_t tpcncls, const Double_t tpcsig)
 269 {
 270   //
 271   //the scale used in calibration
 272   //
 273
 274   if(tpcncls < AliTRDdEdxBaseUtils::CalibTPCnclsCut())
 275     return -999;
 276
 277   if(tpcsig<EPSILON)
 278     return -999;
 279
 280   return tpcsig/120;
 281
 282 }
 283
 284 void AliTRDdEdxCalibUtils::GetPHCountMeanRMS(const TH1D *hnor, TH1D *&hmean)
 285 {
 286   //
 287   //calculate from the ph calib hist the (mean-3sigma) ph-count in the chamber, save in the TH1D output
 288   //
 289   const Int_t ndet = 540;
 290   TObjArray *obj=new TObjArray(ndet);
 291   obj->SetOwner();
 292   for(Int_t ii=0; ii<ndet; ii++){
 293     obj->Add(new TVectorD(AliTRDseedV1::kNtb));
 294   }
 295
 296   //ibin = binlowedge of bin(ibin+1) = the number fills this bin
 297   for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
 298     const Double_t stat = hnor->GetBinContent(ibin+1);
 299     if(stat<EPSILON){
 300       continue;
 301     }
 302
 303     const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
 304     const Int_t itb  = AliTRDdEdxBaseUtils::ToTimeBin(ibin);
 305     TVectorD *vec=(TVectorD *)obj->At(idet);
 306     (*vec)[itb] = stat;
 307   }
 308
 309   hmean = new TH1D(Form("%sdetmean", hnor->GetName()), "", hnor->GetNbinsX(), hnor->GetXaxis()->GetXmin(), hnor->GetXaxis()->GetXmax());
 310   for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
 311     const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
 312     const TVectorD *vec=(TVectorD *)obj->At(idet);
 313
 314     Int_t nonzero = 0;
 315     for(Int_t ii=0; ii<vec->GetNrows(); ii++){
 316       if((*vec)[ii]>EPSILON){
 317         nonzero++;
 318       }
 319     }
 320
 321     Double_t mean = 0;
 322     const Double_t lowfrac = 0.6;
 323     //if there are too many 0's, reject this chamber by settig mean=rms=0
 324     if(nonzero> (AliTRDseedV1::kNtb*(1-lowfrac)) ){
 325       //only highest (1-lowfrac)*31 timebins are used to estimate the mean and rms! important! otherwise the 0' will make rms very large!
 326       mean = AliTRDdEdxBaseUtils::TruncatedMean(AliTRDseedV1::kNtb, vec->GetMatrixArray(), lowfrac, 1);
 327     }
 328
 329     hmean->SetBinContent(ibin+1, mean);
 330   }
 331
 332   delete obj;
 333 }
 334
 335 void AliTRDdEdxCalibUtils::Output(const TList *lin, Int_t run)
 336 {
 337   //
 338   //produce calibration objects
 339   //
 340
 341   TString objnames;
 342   for(Int_t iter=0; iter<8; iter++){
 343     objnames+= AliTRDdEdxCalibHistArray::GetNameAt(iter)+" ";
 344   }
 345
 346   TList * lout = new TList;
 347   lout->SetOwner();
 348
 349   TTreeSRedirector *calibStream = new TTreeSRedirector(Form("TRDCalibStream_%010d.root", run));
 350
 351   const Int_t nh=lin->GetEntries();
 352   for(Int_t ii=0; ii<nh; ii++){
 353     const THnBase *hh=(THnBase*)lin->At(ii);
 354     const TString hname = hh->GetName();
 355     if(!objnames.Contains(hname))
 356       continue;
 357
 358     TObjArray * cobj0 = HistToObj(hh, run, lout, calibStream);
 359     lout->Add(cobj0);
 360   }
 361
 362   //lout->ls();
 363
 364   //=============================================================
 365   //=============================================================
 366
 367   TFile *fout=new TFile(Form("TRDCalibObj_%010d.root", run),"recreate");
 368   fout->cd();
 369   const Int_t nout=lout->GetEntries();
 370   for(Int_t ii=0; ii<nout; ii++){
 371     const TString oname = lout->At(ii)->GetName();
 372     if(oname.Contains("Obj")){
 373       TObjArray * cobj = (TObjArray*) lout->At(ii);
 374       cobj->Write(oname, TObjArray::kSingleKey);
 375     }
 376   }
 377   fout->Save();
 378   fout->Close();
 379   delete fout;
 380
 381   fout=new TFile(Form("TRDCalibList_%010d.root", run),"recreate");
 382   fout->cd();
 383   lin->Write();
 384   lout->Write();
 385   fout->Save();
 386   fout->Close();
 387   delete fout;
 388
 389   delete calibStream;
 390
 391   /*
 392     http://root.cern.ch/root/html/TH1.html
 393     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:
 394
 395     h->SetDirectory(0);          for the current histogram h
 396     TH1::AddDirectory(kFALSE);   sets a global switch disabling the reference
 397
 398     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.
 399   */
 400   delete lout;
 401 }
 402
 403 TObjArray* AliTRDdEdxCalibUtils::HistToObj(const THnBase *hh, Int_t run, TList *lout, TTreeSRedirector *calibStream)
 404 {
 405   //
 406   //produce calibration objects
 407   //
 408
 409   const TString hname = hh->GetName();
 410   const Bool_t kinvq = TString(hname(hname.First('Q')+1,1)).Atoi()&4;
 411
 412   //----------------------------------------
 413   //               Define nbin, tag, cobj0
 414   //----------------------------------------
 415   const Int_t nbin = AliTRDdEdxBaseUtils::NTRDtimebin();
 416
 417
 418   TString tag(hname);
 419   tag.ReplaceAll("Hist","Obj");
 420
 421   TObjArray * cobj0 = new TObjArray(1);
 422   cobj0->SetOwner();
 423   cobj0->SetName(tag);
 424   cobj0->Add(new TVectorD(nbin));
 425
 426   //----------------------------------------
 427   //               Define lowFrac, highFrac
 428   //----------------------------------------
 429   Double_t lowFrac = -999, highFrac = -999;
 430   if(!kinvq) {
 431     lowFrac = 0.01; highFrac = AliTRDdEdxBaseUtils::Q0Frac();
 432   }
 433   else{
 434     lowFrac = AliTRDdEdxBaseUtils::Q1Frac(); highFrac = 0.99;
 435   }
 436
 437   //----------------------------------------
 438   //               Get analysis result
 439   //----------------------------------------
 440   TH1::AddDirectory(kFALSE);//important!
 441   TH1D *hnor=0x0, *hmpv=0x0, *hres=0x0, *hwid=0x0, *htrdphmean = 0x0;//if(!lout), these have to be deleted
 442   TH2D *hpj = hh->Projection(1,0);
 443   AliTRDdEdxBaseUtils::FitSlicesY(hpj, hnor, hmpv, hwid, hres, 0, lowFrac, highFrac);
 444   GetPHCountMeanRMS(hnor, htrdphmean);
 445   if(lout){
 446     lout->Add(htrdphmean);
 447   }
 448   delete hpj;
 449
 450   if(lout){
 451     lout->Add(hnor);
 452     lout->Add(hmpv);
 453     lout->Add(hwid);
 454     lout->Add(hres);
 455   }
 456
 457   //----------------------------------------
 458   //               Define Counter
 459   //----------------------------------------
 460   TVectorD *countDet=0x0;
 461   TObjArray *countSSL=0x0;
 462
 463   if(!kinvq){
 464     countDet=new TVectorD(540);
 465     countSSL=new TObjArray(90);//SectorStackLayer
 466     countSSL->SetOwner();
 467     for(Int_t ii=0; ii<90; ii++){
 468       countSSL->Add(new TVectorD(6));
 469     }
 470   }
 471
 472   //----------------------------------------
 473   //               Fill cobj0
 474   //----------------------------------------
 475
 476   //ibin = binlowedge of bin(ibin+1) = the number fills this bin
 477   for(Int_t ibin=0; ibin<nbin; ibin++){
 478     Double_t gnor = hnor->GetBinContent(ibin+1);
 479     Double_t gmpv = hmpv->GetBinContent(ibin+1);
 480     Double_t gwid = hwid->GetBinContent(ibin+1);
 481     Double_t gres = hres->GetBinContent(ibin+1);
 482
 483     //--- set additional cut by kpass
 484     Bool_t kpass = kTRUE;
 485     Double_t gtrdphmean = -999;
 486     if(htrdphmean){
 487       gtrdphmean = htrdphmean->GetBinContent(ibin+1);
 488       //chamber no statistics (e.g. too many 0's), not usual, not seen in run 143237
 489       if(gtrdphmean<EPSILON){
 490         kpass = kFALSE;
 491       }
 492       if(gnor<AliTRDdEdxBaseUtils::TimeBinCountCut()*gtrdphmean){
 493         kpass = kFALSE;
 494       }
 495     }
 496
 497     //--- set calibration constant p0
 498     Double_t p0= 0;
 499
 500     //reason for gmpv=0:
 501     //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;
 502     //2)TruncatedMean(hist) out of range (only for Q0, not Q1).
 503
 504     if(gmpv>EPSILON && kpass){
 505       if(tag.Contains("T0")){
 506         p0 = gmpv;
 507       }
 508       else{
 509         p0 = 1/gmpv;
 510       }
 511       //printf("outcalibobj%s %d %15.6e\n", tag.Data(), ibin, p0);
 512     }
 513
 514     (*( (TVectorD*)cobj0->At(0) ))[ibin] = p0;
 515
 516     //--- save optional record
 517     if(p0>EPSILON && countDet && countSSL){
 518       const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
 519       (*countDet)[idet]=1;
 520
 521       const Int_t isector = AliTRDdEdxBaseUtils::ToSector(ibin);
 522       const Int_t istack = AliTRDdEdxBaseUtils::ToStack(ibin);
 523       const Int_t ilayer = AliTRDdEdxBaseUtils::ToLayer(ibin);
 524       TVectorD * vecsectorstack = (TVectorD*)countSSL->At(istack*18+isector);
 525       (*vecsectorstack)[ilayer]=1;
 526     }
 527
 528     if(calibStream){
 529       (*calibStream)<<tag<<
 530         "run="<<run<<
 531         "p0="<<p0<<
 532
 533         "nor="<<gnor<<
 534         "mpv="<<gmpv<<
 535         "wid="<<gwid<<
 536         "res="<<gres<<
 537         "gtrdphmean="<<gtrdphmean<<
 538
 539         "ibin="<<ibin<<
 540         "\n";
 541     }
 542   }
 543
 544   //----------------------------------------
 545   //               Status Report
 546   //----------------------------------------
 547   if(countDet && countSSL){
 548     TVectorD count2Dstack(90);
 549     for(Int_t ii=0; ii<90; ii++){
 550       TVectorD * vecsectorstack = (TVectorD*)countSSL->At(ii);
 551       const Int_t nlayer = (Int_t)vecsectorstack->Sum();
 552       if(nlayer==6){
 553         count2Dstack[ii]=1;
 554       }
 555     }
 556
 557     printf("\nAliTRDdEdxCalibUtils::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());
 558   }
 559
 560   //----------------------------------------
 561   //               Clean Up
 562   //----------------------------------------
 563
 564   TH1D **hhs[]={&hnor, &hmpv, &hwid, &hres, &htrdphmean};
 565   const Int_t nhh=sizeof(hhs)/sizeof(TH1D**);
 566   for(Int_t ihh=0; ihh<nhh; ihh++){
 567     if(!lout){
 568       delete (*hhs[ihh]);
 569     }
 570   }
 571
 572   delete countDet;
 573   delete countSSL;
 574
 575   //----------------------------------------
 576
 577   return cobj0;
 578 }
 579