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     Double_t dq = (*arrayQ)[ii]/scale;
 222     const Double_t qmin = hcalib->GetAxis(1)->GetXmin() +0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
 223     const Double_t qmax = hcalib->GetAxis(1)->GetXmax() -0.5 * hcalib->GetAxis(1)->GetBinWidth(1);
 224
 225     if(dq<qmin)
 226       dq = qmin;
 227     if(dq>qmax)
 228       dq = qmax;
 229
 230     const Double_t xx = (*arrayX)[ii];
 231     const Double_t xmin = hcalib->GetAxis(0)->GetXmin();
 232     const Double_t xmax = hcalib->GetAxis(0)->GetXmax();
 233
 234     if(xx>=xmax || xx<xmin){
 235       printf("AliTRDdEdxCalibUtils::FillCalibHist error x overflow or underflow! %s %15f %15f %15f\n", hcalib->GetName(),  xx, xmin, xmax); exit(1);
 236     }
 237
 238     const Double_t var[]={xx, dq};
 239     hcalib->Fill(var);
 240   }
 241 }
 242
 243 void AliTRDdEdxCalibUtils::FillHist(const AliTRDtrackV1 *trdv1, const Bool_t kinvq, const Double_t mag, const Int_t charge, const Double_t scale)
 244 {
 245   //
 246   //get cluster Q and fill calib hist, if kinvq = kTRUE, 1/Q is filled
 247   //
 248   if(!fgHistArray){
 249     printf("AliTRDdEdxCalibUtils::FillHist fgHistArray not initialized!!\n"); exit(1);
 250   }
 251
 252   THnBase * hcalib = fgHistArray->GetHist(kinvq, mag, charge);
 253
 254   TVectorD arrayQ(200), arrayX(200);
 255   const Int_t ncls = AliTRDdEdxReconUtils::GetArrayClusterQ(kinvq, &arrayQ, &arrayX, trdv1);
 256   FillHist(ncls, &arrayQ, &arrayX, hcalib, kinvq ? 1/scale : scale);
 257
 258   static Int_t kprint = 100;
 259   if(kprint<0){
 260     printf("\nAliTRDdEdxCalibUtils::FillHist summary: \n");
 261     printf("\nkinvq= %d;\n", kinvq);
 262     for(Int_t iq=0; iq<ncls; iq++){
 263       printf("arrayX[%3d] = %15.0f; arrayQ[%3d] =  %15f;\n", iq, arrayX[iq], iq, arrayQ[iq]);
 264     }
 265     printf("\n");
 266   }
 267   kprint++;
 268 }
 269
 270 //============================================================
 271 //
 272 //============================================================
 273
 274 Double_t AliTRDdEdxCalibUtils::GetCalibTPCscale(const Int_t tpcncls, const Double_t tpcsig)
 275 {
 276   //
 277   //the scale used in calibration
 278   //
 279
 280   if(tpcncls < AliTRDdEdxBaseUtils::CalibTPCnclsCut())
 281     return -999;
 282
 283   if(tpcsig<EPSILON)
 284     return -999;
 285
 286   return tpcsig/120;
 287
 288 }
 289
 290 void AliTRDdEdxCalibUtils::GetPHCountMeanRMS(const TH1D *hnor, TH1D *&hmean)
 291 {
 292   //
 293   //calculate from the ph calib hist the (mean-3sigma) ph-count in the chamber, save in the TH1D output
 294   //
 295   const Int_t ndet = 540;
 296   TObjArray *obj=new TObjArray(ndet);
 297   obj->SetOwner();
 298   for(Int_t ii=0; ii<ndet; ii++){
 299     obj->Add(new TVectorD(AliTRDseedV1::kNtb));
 300   }
 301
 302   //ibin = binlowedge of bin(ibin+1) = the number fills this bin
 303   for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
 304     const Double_t stat = hnor->GetBinContent(ibin+1);
 305     if(stat<EPSILON){
 306       continue;
 307     }
 308
 309     const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
 310     const Int_t itb  = AliTRDdEdxBaseUtils::ToTimeBin(ibin);
 311     TVectorD *vec=(TVectorD *)obj->At(idet);
 312     (*vec)[itb] = stat;
 313   }
 314
 315   hmean = new TH1D(Form("%sdetmean", hnor->GetName()), "", hnor->GetNbinsX(), hnor->GetXaxis()->GetXmin(), hnor->GetXaxis()->GetXmax());
 316   for(Int_t ibin=0; ibin<hnor->GetNbinsX(); ibin++){
 317     const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
 318     const TVectorD *vec=(TVectorD *)obj->At(idet);
 319
 320     Int_t nonzero = 0;
 321     for(Int_t ii=0; ii<vec->GetNrows(); ii++){
 322       if((*vec)[ii]>EPSILON){
 323         nonzero++;
 324       }
 325     }
 326
 327     Double_t mean = 0;
 328     const Double_t lowfrac = 0.6;
 329     //if there are too many 0's, reject this chamber by settig mean=rms=0
 330     if(nonzero> (AliTRDseedV1::kNtb*(1-lowfrac)) ){
 331       //only highest (1-lowfrac)*31 timebins are used to estimate the mean and rms! important! otherwise the 0' will make rms very large!
 332       mean = AliTRDdEdxBaseUtils::TruncatedMean(AliTRDseedV1::kNtb, vec->GetMatrixArray(), lowfrac, 1);
 333     }
 334
 335     hmean->SetBinContent(ibin+1, mean);
 336   }
 337
 338   delete obj;
 339 }
 340
 341 void AliTRDdEdxCalibUtils::Output(const TList *lin, Int_t run)
 342 {
 343   //
 344   //produce calibration objects
 345   //
 346
 347   TString objnames;
 348   for(Int_t iter=0; iter<8; iter++){
 349     objnames+= AliTRDdEdxCalibHistArray::GetNameAt(iter)+" ";
 350   }
 351
 352   TList * lout = new TList;
 353   lout->SetOwner();
 354
 355   TTreeSRedirector *calibStream = new TTreeSRedirector(Form("TRDCalibStream_%010d.root", run));
 356
 357   const Int_t nh=lin->GetEntries();
 358   for(Int_t ii=0; ii<nh; ii++){
 359     const THnBase *hh=(THnBase*)lin->At(ii);
 360     const TString hname = hh->GetName();
 361     if(!objnames.Contains(hname))
 362       continue;
 363
 364     TObjArray * cobj0 = HistToObj(hh, run, lout, calibStream);
 365     lout->Add(cobj0);
 366   }
 367
 368   //lout->ls();
 369
 370   //=============================================================
 371   //=============================================================
 372
 373   TFile *fout=new TFile(Form("TRDCalibObj_%010d.root", run),"recreate");
 374   fout->cd();
 375   const Int_t nout=lout->GetEntries();
 376   for(Int_t ii=0; ii<nout; ii++){
 377     const TString oname = lout->At(ii)->GetName();
 378     if(oname.Contains("Obj")){
 379       TObjArray * cobj = (TObjArray*) lout->At(ii);
 380       cobj->Write(oname, TObjArray::kSingleKey);
 381     }
 382   }
 383   fout->Save();
 384   fout->Close();
 385   delete fout;
 386
 387   fout=new TFile(Form("TRDCalibList_%010d.root", run),"recreate");
 388   fout->cd();
 389   lin->Write();
 390   lout->Write();
 391   fout->Save();
 392   fout->Close();
 393   delete fout;
 394
 395   delete calibStream;
 396
 397   /*
 398     http://root.cern.ch/root/html/TH1.html
 399     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:
 400
 401     h->SetDirectory(0);          for the current histogram h
 402     TH1::AddDirectory(kFALSE);   sets a global switch disabling the reference
 403
 404     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.
 405   */
 406   delete lout;
 407 }
 408
 409 TObjArray* AliTRDdEdxCalibUtils::HistToObj(const THnBase *hh, Int_t run, TList *lout, TTreeSRedirector *calibStream)
 410 {
 411   //
 412   //produce calibration objects
 413   //
 414
 415   const TString hname = hh->GetName();
 416   const Bool_t kinvq = TString(hname(hname.First('Q')+1,1)).Atoi()&4;
 417
 418   //----------------------------------------
 419   //               Define nbin, tag, cobj0
 420   //----------------------------------------
 421   const Int_t nbin = AliTRDdEdxBaseUtils::NTRDtimebin();
 422
 423
 424   TString tag(hname);
 425   tag.ReplaceAll("Hist","Obj");
 426
 427   TObjArray * cobj0 = new TObjArray(1);
 428   cobj0->SetOwner();
 429   cobj0->SetName(tag);
 430   cobj0->Add(new TVectorD(nbin));
 431
 432   //----------------------------------------
 433   //               Define lowFrac, highFrac
 434   //----------------------------------------
 435   Double_t lowFrac = -999, highFrac = -999;
 436   if(!kinvq) {
 437     lowFrac = 0.01; highFrac = AliTRDdEdxBaseUtils::Q0Frac();
 438   }
 439   else{
 440     lowFrac = AliTRDdEdxBaseUtils::Q1Frac(); highFrac = 0.99;
 441   }
 442
 443   //----------------------------------------
 444   //               Get analysis result
 445   //----------------------------------------
 446   TH1::AddDirectory(kFALSE);//important!
 447   TH1D *hnor=0x0, *hmpv=0x0, *hres=0x0, *hwid=0x0, *htrdphmean = 0x0;//if(!lout), these have to be deleted
 448   TH2D *hpj = hh->Projection(1,0);
 449   AliTRDdEdxBaseUtils::FitSlicesY(hpj, hnor, hmpv, hwid, hres, 0, lowFrac, highFrac);
 450   GetPHCountMeanRMS(hnor, htrdphmean);
 451   if(lout){
 452     lout->Add(htrdphmean);
 453   }
 454   delete hpj;
 455
 456   if(lout){
 457     lout->Add(hnor);
 458     lout->Add(hmpv);
 459     lout->Add(hwid);
 460     lout->Add(hres);
 461   }
 462
 463   //----------------------------------------
 464   //               Define Counter
 465   //----------------------------------------
 466   TVectorD *countDet=0x0;
 467   TObjArray *countSSL=0x0;
 468
 469   if(!kinvq){
 470     countDet=new TVectorD(540);
 471     countSSL=new TObjArray(90);//SectorStackLayer
 472     countSSL->SetOwner();
 473     for(Int_t ii=0; ii<90; ii++){
 474       countSSL->Add(new TVectorD(6));
 475     }
 476   }
 477
 478   //----------------------------------------
 479   //               Fill cobj0
 480   //----------------------------------------
 481
 482   //ibin = binlowedge of bin(ibin+1) = the number fills this bin
 483   for(Int_t ibin=0; ibin<nbin; ibin++){
 484     Double_t gnor = hnor->GetBinContent(ibin+1);
 485     Double_t gmpv = hmpv->GetBinContent(ibin+1);
 486     Double_t gwid = hwid->GetBinContent(ibin+1);
 487     Double_t gres = hres->GetBinContent(ibin+1);
 488
 489     //--- set additional cut by kpass
 490     Bool_t kpass = kTRUE;
 491     Double_t gtrdphmean = -999;
 492     if(htrdphmean){
 493       gtrdphmean = htrdphmean->GetBinContent(ibin+1);
 494       //chamber no statistics (e.g. too many 0's), not usual, not seen in run 143237
 495       if(gtrdphmean<EPSILON){
 496         kpass = kFALSE;
 497       }
 498       if(gnor<AliTRDdEdxBaseUtils::TimeBinCountCut()*gtrdphmean){
 499         kpass = kFALSE;
 500       }
 501     }
 502
 503     //--- set calibration constant p0
 504     Double_t p0= 0;
 505
 506     //reason for gmpv=0:
 507     //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;
 508     //2)TruncatedMean(hist) out of range (only for Q0, not Q1).
 509
 510     if(gmpv>EPSILON && kpass){
 511       if(tag.Contains("T0")){
 512         p0 = gmpv;
 513       }
 514       else{
 515         p0 = 1/gmpv;
 516       }
 517       //printf("outcalibobj%s %d %15.6e\n", tag.Data(), ibin, p0);
 518     }
 519
 520     (*( (TVectorD*)cobj0->At(0) ))[ibin] = p0;
 521
 522     //--- save optional record
 523     if(p0>EPSILON && countDet && countSSL){
 524       const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(ibin);
 525       (*countDet)[idet]=1;
 526
 527       const Int_t isector = AliTRDdEdxBaseUtils::ToSector(ibin);
 528       const Int_t istack = AliTRDdEdxBaseUtils::ToStack(ibin);
 529       const Int_t ilayer = AliTRDdEdxBaseUtils::ToLayer(ibin);
 530       TVectorD * vecsectorstack = (TVectorD*)countSSL->At(istack*18+isector);
 531       (*vecsectorstack)[ilayer]=1;
 532     }
 533
 534     if(calibStream){
 535       (*calibStream)<<tag<<
 536         "run="<<run<<
 537         "p0="<<p0<<
 538
 539         "nor="<<gnor<<
 540         "mpv="<<gmpv<<
 541         "wid="<<gwid<<
 542         "res="<<gres<<
 543         "gtrdphmean="<<gtrdphmean<<
 544
 545         "ibin="<<ibin<<
 546         "\n";
 547     }
 548   }
 549
 550   //----------------------------------------
 551   //               Status Report
 552   //----------------------------------------
 553   if(countDet && countSSL){
 554     TVectorD count2Dstack(90);
 555     for(Int_t ii=0; ii<90; ii++){
 556       TVectorD * vecsectorstack = (TVectorD*)countSSL->At(ii);
 557       const Int_t nlayer = (Int_t)vecsectorstack->Sum();
 558       if(nlayer==6){
 559         count2Dstack[ii]=1;
 560       }
 561     }
 562
 563     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());
 564   }
 565
 566   //----------------------------------------
 567   //               Clean Up
 568   //----------------------------------------
 569
 570   TH1D **hhs[]={&hnor, &hmpv, &hwid, &hres, &htrdphmean};
 571   const Int_t nhh=sizeof(hhs)/sizeof(TH1D**);
 572   for(Int_t ihh=0; ihh<nhh; ihh++){
 573     if(!lout){
 574       delete (*hhs[ihh]);
 575     }
 576   }
 577
 578   delete countDet;
 579   delete countSSL;
 580
 581   //----------------------------------------
 582
 583   return cobj0;
 584 }
 585