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