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