TPC/AliTPCCalibPedestal.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 /* $Id$ */
  18
  19
  20 //Root includes
  21 #include <TObjArray.h>
  22 #include <TH1F.h>
  23 #include <TH2F.h>
  24 #include <TString.h>
  25 #include <TMath.h>
  26 #include <TF1.h>
  27 #include <TRandom.h>
  28 #include <TDirectory.h>
  29 #include <TFile.h>
  30 //AliRoot includes
  31 #include "AliRawReader.h"
  32 #include "AliRawReaderRoot.h"
  33 #include "AliRawReaderDate.h"
  34 #include "AliTPCRawStream.h"
  35 #include "AliTPCCalROC.h"
  36 #include "AliTPCROC.h"
  37 #include "AliMathBase.h"
  38 #include "TTreeStream.h"
  39
  40 //date
  41 #include "event.h"
  42
  43 //header file
  44 #include "AliTPCCalibPedestal.h"
  45
  46
  47 ///////////////////////////////////////////////////////////////////////////////////////
  48 //          Implementation of the TPC pedestal and noise calibration
  49 //
  50 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
  51 //
  52 //
  53 // *************************************************************************************
  54 // *                                Class Description                                  *
  55 // *************************************************************************************
  56 //
  57 // Working principle:
  58 // ------------------
  59 // Raw pedestal data is processed by calling one of the ProcessEvent(...) functions
  60 // (see below). These in the end call the Update(...) function, where the data is filled
  61 // into histograms.
  62 //
  63 // For each ROC one TH2F histo (ROC channel vs. ADC channel) is created when
  64 // it is filled for the first time (GetHistoPedestal(ROC,kTRUE)). All histos are stored in the
  65 // TObjArray fHistoPedestalArray.
  66 //
  67 // For a fast filling of the histogram the corresponding bin number of the channel and ADC channel
  68 // is computed by hand and the histogram array is accessed directly via its pointer.
  69 // ATTENTION: Doing so the the entry counter of the histogram is not increased
  70 //            this means that e.g. the colz draw option gives an empty plot unless
  71 //          calling 'histo->SetEntries(1)' before drawing.
  72 //
  73 // After accumulating the desired statistics the Analyse() function has to be called.
  74 // Whithin this function the pedestal and noise values are calculated for each pad, using
  75 // the fast gaus fit function  AliMathBase::FitGaus(...), and the calibration
  76 // storage classes (AliTPCCalROC) are filled for each ROC.
  77 // The calibration information is stored in the TObjArrays fCalRocArrayPedestal and fCalRocArrayRMS;
  78 //
  79 //
  80 //
  81 // User interface for filling data:
  82 // --------------------------------
  83 //
  84 // To Fill information one of the following functions can be used:
  85 //
  86 // Bool_t ProcessEvent(eventHeaderStruct *event);
  87 //   - process Date event
  88 //   - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
  89 //
  90 // Bool_t ProcessEvent(AliRawReader *rawReader);
  91 //  - process AliRawReader event
  92 //   - use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)
  93 //
  94 // Bool_t ProcessEvent(AliTPCRawStream *rawStream);
  95 //   - process event from AliTPCRawStream
  96 //   - call Update function for signal filling
  97 //
  98 // Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
  99 //              iPad,  const Int_t iTimeBin, const Float_t signal);
 100 //   - directly  fill signal information (sector, row, pad, time bin, pad)
 101 //     to the reference histograms
 102 //
 103 // It is also possible to merge two independently taken calibrations using the function
 104 //
 105 // void Merge(AliTPCCalibPedestal *ped)
 106 //   - copy histograms in 'ped' if the do not exist in this instance
 107 //   - Add histograms in 'ped' to the histograms in this instance if the allready exist
 108 //   - After merging call Analyse again!
 109 //
 110 //
 111 //
 112 // -- example: filling data using root raw data:
 113 // void fillPedestal(Char_t *filename)
 114 // {
 115 //    rawReader = new AliRawReaderRoot(fileName);
 116 //    if ( !rawReader ) return;
 117 //    AliTPCCalibPedestal *calib = new AliTPCCalibPedestal;
 118 //    while (rawReader->NextEvent()){
 119 //      calib->ProcessEvent(rawReader);
 120 //    }
 121 //    calib->Analyse();
 122 //    calib->DumpToFile("PedestalData.root");
 123 //    delete rawReader;
 124 //    delete calib;
 125 // }
 126 //
 127 //
 128 // What kind of information is stored and how to retrieve them:
 129 // ------------------------------------------------------------
 130 //
 131 // - Accessing the 'Reference Histograms' (pedestal distribution histograms):
 132 //
 133 // TH2F *GetHistoPedestal(Int_t sector);
 134 //
 135 // - Accessing the calibration storage objects:
 136 //
 137 // AliTPCCalROC *GetCalRocPedestal(Int_t sector);  - for the pedestal values
 138 // AliTPCCalROC *GetCalRocNoise(Int_t sector);     - for the Noise values
 139 //
 140 // example for visualisation:
 141 // if the file "PedestalData.root" was created using the above example one could do the following:
 142 //
 143 // TFile filePedestal("PedestalData.root")
 144 // AliTPCCalibPedestal *ped = (AliTPCCalibPedestal*)filePedestal->Get("AliTPCCalibPedestal");
 145 // ped->GetCalRocPedestal(0)->Draw("colz");
 146 // ped->GetCalRocRMS(0)->Draw("colz");
 147 //
 148 // or use the AliTPCCalPad functionality:
 149 // AliTPCCalPad padPedestal(ped->GetCalPadPedestal());
 150 // AliTPCCalPad padNoise(ped->GetCalPadRMS());
 151 // padPedestal->MakeHisto2D()->Draw("colz");  //Draw A-Side Pedestal Information
 152 // padNoise->MakeHisto2D()->Draw("colz");  //Draw A-Side Noise Information
 153 //
 154 /*
 155  example: fill pedestal with gausschen noise
 156  AliTPCCalibPedestal ped;
 157  ped.TestEvent();
 158  ped.Analyse();
 159  //Draw output;
 160  TCanvas* c1 = new TCanvas;
 161  c1->Divide(1,2);
 162  c1->cd(1);
 163  ped.GetHistoPedestal(0)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
 164  ped.GetHistoPedestal(0)->Draw("colz");
 165  c1->cd(2);
 166  ped.GetHistoPedestal(36)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
 167  ped.GetHistoPedestal(36)->Draw("colz");
 168  TCanvas* c2 = new TCanvas;
 169  c2->Divide(2,2);
 170  c2->cd(1);
 171  ped.GetCalRocPedestal(0)->Draw("colz");
 172  c2->cd(2);
 173  ped.GetCalRocRMS(0)->Draw("colz");
 174  c2->cd(3);
 175  ped.GetCalRocPedestal(36)->Draw("colz");
 176  c2->cd(4);
 177  ped.GetCalRocRMS(36)->Draw("colz");
 178
 179
 180 */
 181
 182
 183
 184 ClassImp(AliTPCCalibPedestal)
 185
 186 AliTPCCalibPedestal::AliTPCCalibPedestal() : /*FOLD00*/
 187   TObject(),
 188   fFirstTimeBin(60),
 189   fLastTimeBin(1000),
 190   fAdcMin(1),
 191   fAdcMax(100),
 192   fOldRCUformat(kTRUE),
 193   fROC(AliTPCROC::Instance()),
 194   fCalRocArrayPedestal(72),
 195   fCalRocArrayRMS(72),
 196   fHistoPedestalArray(72)
 197 {
 198     //
 199     // default constructor
 200     //
 201 }
 202 //_____________________________________________________________________
 203 AliTPCCalibPedestal::AliTPCCalibPedestal(const AliTPCCalibPedestal &ped) : /*FOLD00*/
 204   TObject(ped),
 205   fFirstTimeBin(ped.GetFirstTimeBin()),
 206   fLastTimeBin(ped.GetLastTimeBin()),
 207   fAdcMin(ped.GetAdcMin()),
 208   fAdcMax(ped.GetAdcMax()),
 209   fOldRCUformat(kTRUE),
 210   fROC(AliTPCROC::Instance()),
 211   fCalRocArrayPedestal(72),
 212   fCalRocArrayRMS(72),
 213   fHistoPedestalArray(72)
 214 {
 215     //
 216     // copy constructor
 217     //
 218     for (Int_t iSec = 0; iSec < 72; ++iSec){
 219         const AliTPCCalROC *calPed = (AliTPCCalROC*)ped.fCalRocArrayPedestal.UncheckedAt(iSec);
 220         const AliTPCCalROC *calRMS = (AliTPCCalROC*)ped.fCalRocArrayRMS.UncheckedAt(iSec);
 221         const TH2F         *hPed   = (TH2F*)ped.fHistoPedestalArray.UncheckedAt(iSec);
 222
 223         if ( calPed != 0x0 ) fCalRocArrayPedestal.AddAt(new AliTPCCalROC(*calPed), iSec);
 224         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
 225
 226         if ( hPed != 0x0 ){
 227             TH2F *hNew = new TH2F(*hPed);
 228             hNew->SetDirectory(0);
 229             fHistoPedestalArray.AddAt(hNew,iSec);
 230         }
 231     }
 232 }
 233 //_____________________________________________________________________
 234 AliTPCCalibPedestal& AliTPCCalibPedestal::operator = (const  AliTPCCalibPedestal &source)
 235 {
 236   //
 237   // assignment operator
 238   //
 239   if (&source == this) return *this;
 240   new (this) AliTPCCalibPedestal(source);
 241
 242   return *this;
 243 }
 244 //_____________________________________________________________________
 245 AliTPCCalibPedestal::~AliTPCCalibPedestal() /*FOLD00*/
 246 {
 247   //
 248   // destructor
 249   //
 250
 251     fCalRocArrayPedestal.Delete();
 252     fCalRocArrayRMS.Delete();
 253     fHistoPedestalArray.Delete();
 254     delete fROC;
 255 }
 256 //_____________________________________________________________________
 257 Int_t AliTPCCalibPedestal::Update(const Int_t icsector, /*FOLD00*/
 258                                 const Int_t icRow,
 259                                 const Int_t icPad,
 260                                 const Int_t icTimeBin,
 261                                 const Float_t csignal)
 262 {
 263     //
 264     // Signal filling methode
 265     //
 266
 267     //return if we are out of the specified time bin or adc range
 268     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
 269     if ( ((Int_t)csignal>fAdcMax) || ((Int_t)csignal<fAdcMin)  ) return 0;
 270
 271     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 272
 273     // fast filling methode.
 274     // Attention: the entry counter of the histogram is not increased
 275     //            this means that e.g. the colz draw option gives an empty plot
 276     Int_t bin = (iChannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
 277
 278     GetHistoPedestal(icsector,kTRUE)->GetArray()[bin]++;
 279
 280     return 0;
 281 }
 282 //_____________________________________________________________________
 283 Bool_t AliTPCCalibPedestal::ProcessEvent(AliTPCRawStream *rawStream)
 284 {
 285   //
 286   // Event Processing loop - AliTPCRawStream
 287   //
 288
 289   rawStream->SetOldRCUFormat(fOldRCUformat);
 290
 291   Bool_t withInput = kFALSE;
 292
 293   while (rawStream->Next()) {
 294
 295       Int_t isector  = rawStream->GetSector();                       //  current sector
 296       Int_t iRow     = rawStream->GetRow();                          //  current row
 297       Int_t iPad     = rawStream->GetPad();                          //  current pad
 298       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
 299       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
 300
 301       Update(isector,iRow,iPad,iTimeBin,signal);
 302       withInput = kTRUE;
 303   }
 304
 305   return withInput;
 306 }
 307 //_____________________________________________________________________
 308 Bool_t AliTPCCalibPedestal::ProcessEvent(AliRawReader *rawReader)
 309 {
 310   //
 311   //  Event processing loop - AliRawReader
 312   //
 313
 314
 315   AliTPCRawStream rawStream(rawReader);
 316
 317   rawReader->Select("TPC");
 318
 319   return ProcessEvent(&rawStream);
 320 }
 321 //_____________________________________________________________________
 322 Bool_t AliTPCCalibPedestal::ProcessEvent(eventHeaderStruct *event)
 323 {
 324   //
 325   //  process date event
 326   //
 327     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 328     Bool_t result=ProcessEvent(rawReader);
 329     delete rawReader;
 330     return result;
 331 }
 332 //_____________________________________________________________________
 333 Bool_t AliTPCCalibPedestal::TestEvent() /*FOLD00*/
 334 {
 335   //
 336   //  Test event loop
 337   // fill one oroc and one iroc with random gaus
 338   //
 339
 340     gRandom->SetSeed(0);
 341
 342     for (UInt_t iSec=0; iSec<72; ++iSec){
 343         if (iSec%36>0) continue;
 344         for (UInt_t iRow=0; iRow < fROC->GetNRows(iSec); ++iRow){
 345             for (UInt_t iPad=0; iPad < fROC->GetNPads(iSec,iRow); ++iPad){
 346                 for (UInt_t iTimeBin=0; iTimeBin<1024; ++iTimeBin){
 347                     Float_t signal=(Int_t)(iRow+3+gRandom->Gaus(0,.7));
 348                     if ( signal>0 )Update(iSec,iRow,iPad,iTimeBin,signal);
 349                 }
 350             }
 351         }
 352     }
 353     return kTRUE;
 354 }
 355 //_____________________________________________________________________
 356 TH2F* AliTPCCalibPedestal::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
 357                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
 358                                   Char_t *type, Bool_t force)
 359 {
 360     //
 361     // return pointer to Q histogram
 362     // if force is true create a new histogram if it doesn't exist allready
 363     //
 364     if ( !force || arr->UncheckedAt(sector) )
 365         return (TH2F*)arr->UncheckedAt(sector);
 366
 367     // if we are forced and histogram doesn't yes exist create it
 368     Char_t name[255], title[255];
 369
 370     sprintf(name,"hCalib%s%.2d",type,sector);
 371     sprintf(title,"%s calibration histogram sector %.2d;ADC channel;Channel (pad)",type,sector);
 372
 373     // new histogram with Q calib information. One value for each pad!
 374     TH2F* hist = new TH2F(name,title,
 375                           nbinsY, ymin, ymax,
 376                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector)
 377                          );
 378     hist->SetDirectory(0);
 379     arr->AddAt(hist,sector);
 380     return hist;
 381 }
 382 //_____________________________________________________________________
 383 TH2F* AliTPCCalibPedestal::GetHistoPedestal(Int_t sector, Bool_t force) /*FOLD00*/
 384 {
 385     //
 386     // return pointer to T0 histogram
 387     // if force is true create a new histogram if it doesn't exist allready
 388     //
 389     TObjArray *arr = &fHistoPedestalArray;
 390     return GetHisto(sector, arr, fAdcMax-fAdcMin, fAdcMin, fAdcMax, "Pedestal", force);
 391 }
 392 //_____________________________________________________________________
 393 AliTPCCalROC* AliTPCCalibPedestal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
 394 {
 395     //
 396     // return pointer to ROC Calibration
 397     // if force is true create a new histogram if it doesn't exist allready
 398     //
 399     if ( !force || arr->UncheckedAt(sector) )
 400         return (AliTPCCalROC*)arr->UncheckedAt(sector);
 401
 402     // if we are forced and the histogram doesn't yet exist create it
 403
 404     // new AliTPCCalROC for T0 information. One value for each pad!
 405     AliTPCCalROC *croc = new AliTPCCalROC(sector);
 406     arr->AddAt(croc,sector);
 407     return croc;
 408 }
 409 //_____________________________________________________________________
 410 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocPedestal(Int_t sector, Bool_t force) /*FOLD00*/
 411 {
 412     //
 413     // return pointer to Carge ROC Calibration
 414     // if force is true create a new histogram if it doesn't exist allready
 415     //
 416     TObjArray *arr = &fCalRocArrayPedestal;
 417     return GetCalRoc(sector, arr, force);
 418 }
 419 //_____________________________________________________________________
 420 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
 421 {
 422     //
 423     // return pointer to signal width ROC Calibration
 424     // if force is true create a new histogram if it doesn't exist allready
 425     //
 426     TObjArray *arr = &fCalRocArrayRMS;
 427     return GetCalRoc(sector, arr, force);
 428 }
 429 //_____________________________________________________________________
 430 void AliTPCCalibPedestal::Merge(AliTPCCalibPedestal *ped)
 431 {
 432     //
 433     //  Merge reference histograms of sig to the current AliTPCCalibSignal
 434     //
 435
 436     //merge histograms
 437     for (Int_t iSec=0; iSec<72; ++iSec){
 438         TH2F *hRefPedMerge   = ped->GetHistoPedestal(iSec);
 439
 440
 441         if ( hRefPedMerge ){
 442             TDirectory *dir = hRefPedMerge->GetDirectory(); hRefPedMerge->SetDirectory(0);
 443             TH2F *hRefPed   = GetHistoPedestal(iSec);
 444             if ( hRefPed ) hRefPed->Add(hRefPedMerge);
 445             else {
 446                 TH2F *hist = new TH2F(*hRefPedMerge);
 447                 hist->SetDirectory(0);
 448                 fHistoPedestalArray.AddAt(hist, iSec);
 449             }
 450             hRefPedMerge->SetDirectory(dir);
 451         }
 452     }
 453 }
 454 //_____________________________________________________________________
 455 void AliTPCCalibPedestal::Analyse() /*FOLD00*/
 456 {
 457     //
 458     //  Calculate calibration constants
 459     //
 460
 461     Int_t nbinsAdc = fAdcMax-fAdcMin;
 462
 463     TVectorD param(3);
 464     TMatrixD dummy(3,3);
 465
 466     Float_t *array_hP=0;
 467
 468
 469     for (Int_t iSec=0; iSec<72; ++iSec){
 470         TH2F *hP = GetHistoPedestal(iSec);
 471         if ( !hP ) continue;
 472
 473         AliTPCCalROC *rocPedestal = GetCalRocPedestal(iSec,kTRUE);
 474         AliTPCCalROC *rocRMS      = GetCalRocRMS(iSec,kTRUE);
 475
 476         array_hP = hP->GetArray();
 477         UInt_t nChannels = fROC->GetNChannels(iSec);
 478
 479         for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
 480             Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
 481             Double_t ret = AliMathBase::FitGaus(array_hP+offset,nbinsAdc,fAdcMin,fAdcMax,&param,&dummy);
 482             // if the fitting failed set noise and pedestal to 0
 483             if ( ret == -4 ) {
 484                 param[1]=0;
 485                 param[2]=0;
 486             }
 487             rocPedestal->SetValue(iChannel,param[1]);
 488             rocRMS->SetValue(iChannel,param[2]);
 489         }
 490     }
 491 }
 492 //_____________________________________________________________________
 493 void AliTPCCalibPedestal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
 494 {
 495     //
 496     //  Write class to file
 497     //
 498
 499     TString sDir(dir);
 500     TString option;
 501
 502     if ( append )
 503         option = "update";
 504     else
 505         option = "recreate";
 506
 507     TDirectory *backup = gDirectory;
 508     TFile f(filename,option.Data());
 509     f.cd();
 510     if ( !sDir.IsNull() ){
 511         f.mkdir(sDir.Data());
 512         f.cd(sDir);
 513     }
 514     this->Write();
 515     f.Close();
 516
 517     if ( backup ) backup->cd();
 518 }