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 <TH1F.h>
  22 #include <TH2F.h>
  23 #include <TString.h>
  24 #include <TMath.h>
  25 #include <TF1.h>
  26 #include <TRandom.h>
  27 #include <TDirectory.h>
  28 #include <TFile.h>
  29 #include <TMap.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 #include "AliTPCRawStreamFast.h"
  40
  41 //date
  42 #include "event.h"
  43
  44 //header file
  45 #include "AliTPCCalibPedestal.h"
  46
  47
  48 ///////////////////////////////////////////////////////////////////////////////////////
  49 //          Implementation of the TPC pedestal and noise calibration
  50 //
  51 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
  52 //
  53 //
  54 // *************************************************************************************
  55 // *                                Class Description                                  *
  56 // *************************************************************************************
  57 //
  58 // Working principle:
  59 // ------------------
  60 // Raw pedestal data is processed by calling one of the ProcessEvent(...) functions
  61 // (see below). These in the end call the Update(...) function, where the data is filled
  62 // into histograms.
  63 //
  64 // For each ROC one TH2F histo (ROC channel vs. ADC channel) is created when
  65 // it is filled for the first time (GetHistoPedestal(ROC,kTRUE)). All histos are stored in the
  66 // TObjArray fHistoPedestalArray.
  67 //
  68 // For a fast filling of the histogram the corresponding bin number of the channel and ADC channel
  69 // is computed by hand and the histogram array is accessed directly via its pointer.
  70 // ATTENTION: Doing so the the entry counter of the histogram is not increased
  71 //            this means that e.g. the colz draw option gives an empty plot unless
  72 //          calling 'histo->SetEntries(1)' before drawing.
  73 //
  74 // After accumulating the desired statistics the Analyse() function has to be called.
  75 // Whithin this function the pedestal and noise values are calculated for each pad, using
  76 // the fast gaus fit function  AliMathBase::FitGaus(...), and the calibration
  77 // storage classes (AliTPCCalROC) are filled for each ROC.
  78 // The calibration information is stored in the TObjArrays fCalRocArrayPedestal and fCalRocArrayRMS;
  79 //
  80 //
  81 //
  82 // User interface for filling data:
  83 // --------------------------------
  84 //
  85 // To Fill information one of the following functions can be used:
  86 //
  87 // Bool_t ProcessEvent(eventHeaderStruct *event);
  88 //   - process Date event
  89 //   - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
  90 //
  91 // Bool_t ProcessEvent(AliRawReader *rawReader);
  92 //  - process AliRawReader event
  93 //   - use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)
  94 //
  95 // Bool_t ProcessEvent(AliTPCRawStream *rawStream);
  96 //   - process event from AliTPCRawStream
  97 //   - call Update function for signal filling
  98 //
  99 // Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
 100 //              iPad,  const Int_t iTimeBin, const Float_t signal);
 101 //   - directly  fill signal information (sector, row, pad, time bin, pad)
 102 //     to the reference histograms
 103 //
 104 // It is also possible to merge two independently taken calibrations using the function
 105 //
 106 // void Merge(AliTPCCalibPedestal *ped)
 107 //   - copy histograms in 'ped' if the do not exist in this instance
 108 //   - Add histograms in 'ped' to the histograms in this instance if the allready exist
 109 //   - After merging call Analyse again!
 110 //
 111 //
 112 //
 113 // -- example: filling data using root raw data:
 114 // void fillPedestal(Char_t *filename)
 115 // {
 116 //    rawReader = new AliRawReaderRoot(fileName);
 117 //    if ( !rawReader ) return;
 118 //    AliTPCCalibPedestal *calib = new AliTPCCalibPedestal;
 119 //    while (rawReader->NextEvent()){
 120 //      calib->ProcessEvent(rawReader);
 121 //    }
 122 //    calib->Analyse();
 123 //    calib->DumpToFile("PedestalData.root");
 124 //    delete rawReader;
 125 //    delete calib;
 126 // }
 127 //
 128 //
 129 // What kind of information is stored and how to retrieve them:
 130 // ------------------------------------------------------------
 131 //
 132 // - Accessing the 'Reference Histograms' (pedestal distribution histograms):
 133 //
 134 // TH2F *GetHistoPedestal(Int_t sector);
 135 //
 136 // - Accessing the calibration storage objects:
 137 //
 138 // AliTPCCalROC *GetCalRocPedestal(Int_t sector);  - for the pedestal values, mean from gaus fit
 139 // AliTPCCalROC *GetCalRocSigma(Int_t sector);     - for the Noise values, sigma from guas fit
 140 // AliTPCCalROC *GetCalRocMean(Int_t sector);  - for the pedestal values, truncated mean
 141 // AliTPCCalROC *GetCalRocRMS(Int_t sector);     - for the Noise values, rms from truncated mean
 142 //
 143 // example for visualisation:
 144 // if the file "PedestalData.root" was created using the above example one could do the following:
 145 //
 146 // TFile filePedestal("PedestalData.root")
 147 // AliTPCCalibPedestal *ped = (AliTPCCalibPedestal*)filePedestal->Get("AliTPCCalibPedestal");
 148 // ped->GetCalRocPedestal(0)->Draw("colz");
 149 // ped->GetCalRocRMS(0)->Draw("colz");
 150 //
 151 // or use the AliTPCCalPad functionality:
 152 // AliTPCCalPad padPedestal(ped->GetCalPadPedestal());
 153 // AliTPCCalPad padNoise(ped->GetCalPadRMS());
 154 // padPedestal->MakeHisto2D()->Draw("colz");  //Draw A-Side Pedestal Information
 155 // padNoise->MakeHisto2D()->Draw("colz");  //Draw A-Side Noise Information
 156 //
 157 /*
 158  example: fill pedestal with gausschen noise
 159  AliTPCCalibPedestal ped;
 160  ped.TestEvent();
 161  ped.Analyse();
 162  //Draw output;
 163  TCanvas* c1 = new TCanvas;
 164  c1->Divide(1,2);
 165  c1->cd(1);
 166  ped.GetHistoPedestal(0)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
 167  ped.GetHistoPedestal(0)->Draw("colz");
 168  c1->cd(2);
 169  ped.GetHistoPedestal(36)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
 170  ped.GetHistoPedestal(36)->Draw("colz");
 171  TCanvas* c2 = new TCanvas;
 172  c2->Divide(2,2);
 173  c2->cd(1);
 174  ped.GetCalRocPedestal(0)->Draw("colz");
 175  c2->cd(2);
 176  ped.GetCalRocRMS(0)->Draw("colz");
 177  c2->cd(3);
 178  ped.GetCalRocPedestal(36)->Draw("colz");
 179  c2->cd(4);
 180  ped.GetCalRocRMS(36)->Draw("colz");
 181 */
 182 //
 183 // Time dependent pedestals:
 184 //
 185 // If wished there is the possibility to calculate for each channel and time bin
 186 // the mean pedestal [pedestals(t)]. This is done by
 187 //
 188 // 1) setting SetTimeAnalysis(kTRUE),
 189 // 2) processing the data by looping over the events using ProcessEvent(..)
 190 // 3) calling the Analyse() and AnalyseTime(nevents) functions (providing nevents)
 191 // 4) getting the pedestals(t) using   TArrayF **timePed = calibPedestal.GetTimePedestals();
 192 // 5) looking at values using   timePed[row][pad].At(timebin)
 193 //
 194 // This functionality is intended to be used on an LDC bu the detector algorithm
 195 // (TPCPEDESTALda) to generate a data set used for configuration of the pattern
 196 // memory for baseline subtraction in the ALTROs. Later the information should also
 197 // be stored as reference data.
 198 //
 199
 200
 201 ClassImp(AliTPCCalibPedestal)
 202
 203 AliTPCCalibPedestal::AliTPCCalibPedestal() :
 204   TObject(),
 205   fFirstTimeBin(60),
 206   fLastTimeBin(1000),
 207   fAdcMin(1),
 208   fAdcMax(100),
 209   fAnaMeanDown(0.),
 210   fAnaMeanUp(1.),
 211   fTimeAnalysis(kFALSE),
 212   fROC(AliTPCROC::Instance()),
 213   fMapping(NULL),
 214   fCalRocArrayPedestal(72),
 215   fCalRocArraySigma(72),
 216   fHistoPedestalArray(72),
 217   fTimeSignal(NULL),
 218   fCalRocArrayMean(72),
 219   fCalRocArrayRMS(72)
 220 {
 221   //
 222   // default constructor
 223   //
 224 }
 225
 226
 227 //_____________________________________________________________________
 228 AliTPCCalibPedestal::AliTPCCalibPedestal(const AliTPCCalibPedestal &ped) :
 229   TObject(ped),
 230   fFirstTimeBin(ped.GetFirstTimeBin()),
 231   fLastTimeBin(ped.GetLastTimeBin()),
 232   fAdcMin(ped.GetAdcMin()),
 233   fAdcMax(ped.GetAdcMax()),
 234   fAnaMeanDown(ped.fAnaMeanDown),
 235   fAnaMeanUp(ped.fAnaMeanUp),
 236   fTimeAnalysis(ped.fTimeAnalysis),
 237   fROC(AliTPCROC::Instance()),
 238   fMapping(NULL),
 239   fCalRocArrayPedestal(72),
 240   fCalRocArraySigma(72),
 241   fHistoPedestalArray(72),
 242   fTimeSignal(ped.fTimeSignal),
 243   fCalRocArrayMean(72),
 244   fCalRocArrayRMS(72)
 245 {
 246   //
 247   // copy constructor
 248   //
 249   for (Int_t iSec = 0; iSec < 72; ++iSec){
 250     const AliTPCCalROC *calPed = (AliTPCCalROC*)ped.fCalRocArrayPedestal.UncheckedAt(iSec);
 251     const AliTPCCalROC *calRMS = (AliTPCCalROC*)ped.fCalRocArrayRMS.UncheckedAt(iSec);
 252     const TH2F         *hPed   = (TH2F*)ped.fHistoPedestalArray.UncheckedAt(iSec);
 253
 254     if ( calPed != 0x0 ) fCalRocArrayPedestal.AddAt(new AliTPCCalROC(*calPed), iSec);
 255     if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
 256
 257     if ( hPed != 0x0 ){
 258       TH2F *hNew = new TH2F(*hPed);
 259       hNew->SetDirectory(0);
 260       fHistoPedestalArray.AddAt(hNew,iSec);
 261     }
 262   }
 263 }
 264 AliTPCCalibPedestal::AliTPCCalibPedestal(const TMap *config):
 265   TObject(),
 266   fFirstTimeBin(60),
 267   fLastTimeBin(1000),
 268   fAdcMin(1),
 269   fAdcMax(100),
 270   fAnaMeanDown(0.),
 271   fAnaMeanUp(1.),
 272   fTimeAnalysis(kFALSE),
 273   fROC(AliTPCROC::Instance()),
 274   fMapping(NULL),
 275   fCalRocArrayPedestal(72),
 276   fCalRocArraySigma(72),
 277   fHistoPedestalArray(72),
 278   fTimeSignal(NULL),
 279   fCalRocArrayMean(72),
 280   fCalRocArrayRMS(72)
 281 {
 282  //
 283  // This constructor uses a TMap for setting some parametes
 284  //
 285   if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
 286   if (config->GetValue("LastTimeBin"))  fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
 287   if (config->GetValue("AdcMin"))       fAdcMin = ((TObjString*)config->GetValue("AdcMin"))->GetString().Atoi();
 288   if (config->GetValue("AdcMax"))       fAdcMax = ((TObjString*)config->GetValue("AdcMax"))->GetString().Atoi();
 289   if (config->GetValue("TimeAnalysis")) SetTimeAnalysis(((TObjString*)config->GetValue("TimeAnalysis"))->GetString().Atoi());
 290
 291
 292 }
 293
 294
 295 //_____________________________________________________________________
 296 AliTPCCalibPedestal& AliTPCCalibPedestal::operator = (const  AliTPCCalibPedestal &source)
 297 {
 298   //
 299   // assignment operator
 300   //
 301   if (&source == this) return *this;
 302   new (this) AliTPCCalibPedestal(source);
 303
 304   return *this;
 305 }
 306
 307
 308 //_____________________________________________________________________
 309 AliTPCCalibPedestal::~AliTPCCalibPedestal()
 310 {
 311   //
 312   // destructor
 313   //
 314
 315   fCalRocArrayPedestal.Delete();
 316   fCalRocArrayRMS.Delete();
 317   fCalRocArraySigma.Delete();
 318   fHistoPedestalArray.Delete();
 319
 320   if ( fTimeSignal ) {
 321     for (Int_t i = 0; i < 159; i++) {
 322       delete [] fTimeSignal[i];
 323       fTimeSignal[i] = 0;
 324     }
 325     delete [] fTimeSignal;
 326     fTimeSignal = 0;
 327   }
 328
 329   // do not delete fMapping, because we do not own it.
 330
 331 }
 332
 333
 334 //_____________________________________________________________________
 335 void AliTPCCalibPedestal::SetTimeAnalysis(Bool_t time)
 336 {
 337   //
 338   // Use time dependent analysis: Pedestals are analysed as a function
 339   // of the drift time. There is one mean value generated for each time
 340   // bin and each channel. It can be used as reference data and for
 341   // configuration of the ALTRO pattern memory for baseline subtraction.
 342   //
 343   // ATTENTION: Use only on LDC in TPCPEDESTALda! On a LDC we get data
 344   // only from one sector. For the full TPC we would need a lot of
 345   // memory (36*159*140*1024*4bytes = 3.3GB)!
 346   //
 347
 348   fTimeAnalysis = time;
 349
 350   if ( !fTimeAnalysis ) return;
 351
 352   // prepare array for one sector (159*140*1024*4bytes = 92MB):
 353   fTimeSignal = new TArrayF*[159];
 354   for (Int_t i = 0; i < 159; i++) {  // padrows
 355     fTimeSignal[i] = new TArrayF[140];
 356     for (Int_t j = 0; j < 140; j++) {  // pads per row
 357       fTimeSignal[i][j].Set(1024);
 358       for (Int_t k = 0; k < 1024; k++) {  // time bins per pad
 359         fTimeSignal[i][j].AddAt(0., k);
 360       }
 361     }
 362   }
 363 }
 364
 365
 366 //_____________________________________________________________________
 367 Int_t AliTPCCalibPedestal::Update(const Int_t icsector,
 368                                   const Int_t icRow,
 369                                   const Int_t icPad,
 370                                   const Int_t icTimeBin,
 371                                   const Float_t csignal)
 372 {
 373   //
 374   // Signal filling method
 375   //
 376   if (icRow<0) return 0;
 377   if (icPad<0) return 0;
 378   if (icTimeBin<0) return 0;
 379
 380   // Time dependent pedestals
 381   if ( fTimeAnalysis ) {
 382     if ( icsector < 36 ) // IROC
 383       fTimeSignal[icRow][icPad].AddAt(fTimeSignal[icRow][icPad].At(icTimeBin)+csignal, icTimeBin);
 384     else
 385       fTimeSignal[icRow+63][icPad].AddAt(fTimeSignal[icRow+63][icPad].At(icTimeBin)+csignal, icTimeBin);
 386   }
 387   //return if we are out of the specified time bin or adc range
 388   if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
 389   if ( ((Int_t)csignal>fAdcMax) || ((Int_t)csignal<fAdcMin)  ) return 0;
 390
 391   Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 392
 393   // fast filling method
 394   // Attention: the entry counter of the histogram is not increased
 395   //            this means that e.g. the colz draw option gives an empty plot
 396   Int_t bin = (iChannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
 397
 398   GetHistoPedestal(icsector,kTRUE)->GetArray()[bin]++;
 399
 400   return 0;
 401 }
 402
 403
 404 //_____________________________________________________________________
 405 Bool_t AliTPCCalibPedestal::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
 406 {
 407   //
 408   // Event Processing loop - AliTPCRawStream
 409   //
 410   Bool_t withInput = kFALSE;
 411
 412   while ( rawStreamFast->NextDDL() ){
 413       while ( rawStreamFast->NextChannel() ){
 414           Int_t isector  = rawStreamFast->GetSector();                       //  current sector
 415           Int_t iRow     = rawStreamFast->GetRow();                          //  current row
 416           Int_t iPad     = rawStreamFast->GetPad();                          //  current pad
 417
 418           while ( rawStreamFast->NextBunch() ){
 419             Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
 420             Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
 421             for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
 422                   Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
 423                   Update(isector,iRow,iPad,iTimeBin+1,signal);
 424                   withInput = kTRUE;
 425               }
 426           }
 427       }
 428   }
 429
 430   return withInput;
 431 }
 432 //_____________________________________________________________________
 433 Bool_t AliTPCCalibPedestal::ProcessEventFast(AliRawReader *rawReader)
 434 {
 435   //
 436   //  Event processing loop - AliRawReader
 437   //
 438   AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
 439   Bool_t res=ProcessEventFast(rawStreamFast);
 440   delete rawStreamFast;
 441   return res;
 442 }
 443
 444 //_____________________________________________________________________
 445 Bool_t AliTPCCalibPedestal::ProcessEvent(AliTPCRawStream *rawStream)
 446 {
 447   //
 448   // Event Processing loop - AliTPCRawStream
 449   //
 450
 451   Bool_t withInput = kFALSE;
 452
 453   while (rawStream->Next()) {
 454
 455     Int_t iSector  = rawStream->GetSector();      //  current ROC
 456     Int_t iRow     = rawStream->GetRow();         //  current row
 457     Int_t iPad     = rawStream->GetPad();         //  current pad
 458     Int_t iTimeBin = rawStream->GetTime();        //  current time bin
 459     Float_t signal = rawStream->GetSignal();      //  current ADC signal
 460
 461     Update(iSector,iRow,iPad,iTimeBin,signal);
 462     withInput = kTRUE;
 463   }
 464
 465   return withInput;
 466 }
 467
 468
 469 //_____________________________________________________________________
 470 Bool_t AliTPCCalibPedestal::ProcessEvent(AliRawReader *rawReader)
 471 {
 472   //
 473   //  Event processing loop - AliRawReader
 474   //
 475
 476   // if fMapping is NULL the rawstream will crate its own mapping
 477   AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
 478   rawReader->Select("TPC");
 479   return ProcessEvent(&rawStream);
 480 }
 481
 482
 483 //_____________________________________________________________________
 484 Bool_t AliTPCCalibPedestal::ProcessEvent(eventHeaderStruct *event)
 485 {
 486   //
 487   //  process date event
 488   //
 489
 490   AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 491   Bool_t result=ProcessEvent(rawReader);
 492   delete rawReader;
 493   return result;
 494 }
 495
 496
 497 //_____________________________________________________________________
 498 Bool_t AliTPCCalibPedestal::TestEvent()
 499 {
 500   //
 501   //  Test event loop
 502   // fill one oroc and one iroc with random gaus
 503   //
 504
 505     gRandom->SetSeed(0);
 506
 507     for (UInt_t iSec=0; iSec<72; ++iSec){
 508         if (iSec%36>0) continue;
 509         for (UInt_t iRow=0; iRow < fROC->GetNRows(iSec); ++iRow){
 510             for (UInt_t iPad=0; iPad < fROC->GetNPads(iSec,iRow); ++iPad){
 511                 for (UInt_t iTimeBin=0; iTimeBin<1024; ++iTimeBin){
 512                     Float_t signal=(Int_t)(iRow+3+gRandom->Gaus(0,.7));
 513                     if ( signal>0 )Update(iSec,iRow,iPad,iTimeBin,signal);
 514                 }
 515             }
 516         }
 517     }
 518     return kTRUE;
 519 }
 520
 521
 522 //_____________________________________________________________________
 523 TH2F* AliTPCCalibPedestal::GetHisto(Int_t sector, TObjArray *arr,
 524                                     Int_t nbinsY, Float_t ymin, Float_t ymax,
 525                                     Char_t *type, Bool_t force)
 526 {
 527     //
 528     // return pointer to Q histogram
 529     // if force is true create a new histogram if it doesn't exist allready
 530     //
 531     if ( !force || arr->UncheckedAt(sector) )
 532         return (TH2F*)arr->UncheckedAt(sector);
 533
 534     // if we are forced and histogram doesn't yes exist create it
 535     Char_t name[255], title[255];
 536
 537     sprintf(name,"hCalib%s%.2d",type,sector);
 538     sprintf(title,"%s calibration histogram sector %.2d;ADC channel;Channel (pad)",type,sector);
 539
 540     // new histogram with Q calib information. One value for each pad!
 541     TH2F* hist = new TH2F(name,title,
 542                           nbinsY, ymin, ymax,
 543                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector)
 544                          );
 545     hist->SetDirectory(0);
 546     arr->AddAt(hist,sector);
 547     return hist;
 548 }
 549
 550
 551 //_____________________________________________________________________
 552 TH2F* AliTPCCalibPedestal::GetHistoPedestal(Int_t sector, Bool_t force)
 553 {
 554     //
 555     // return pointer to T0 histogram
 556     // if force is true create a new histogram if it doesn't exist allready
 557     //
 558     TObjArray *arr = &fHistoPedestalArray;
 559     return GetHisto(sector, arr, fAdcMax-fAdcMin, fAdcMin, fAdcMax, "Pedestal", force);
 560 }
 561
 562
 563 //_____________________________________________________________________
 564 AliTPCCalROC* AliTPCCalibPedestal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force)
 565 {
 566     //
 567     // return pointer to ROC Calibration
 568     // if force is true create a new histogram if it doesn't exist allready
 569     //
 570     if ( !force || arr->UncheckedAt(sector) )
 571         return (AliTPCCalROC*)arr->UncheckedAt(sector);
 572
 573     // if we are forced and the histogram doesn't yet exist create it
 574
 575     // new AliTPCCalROC for T0 information. One value for each pad!
 576     AliTPCCalROC *croc = new AliTPCCalROC(sector);
 577     arr->AddAt(croc,sector);
 578     return croc;
 579 }
 580
 581
 582 //_____________________________________________________________________
 583 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocPedestal(Int_t sector, Bool_t force)
 584 {
 585     //
 586     // return pointer to ROC with Pedestal data
 587     // if force is true create a new histogram if it doesn't exist allready
 588     //
 589     TObjArray *arr = &fCalRocArrayPedestal;
 590     return GetCalRoc(sector, arr, force);
 591 }
 592
 593
 594 //_____________________________________________________________________
 595 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocSigma(Int_t sector, Bool_t force)
 596 {
 597     //
 598     // return pointer to  ROC with signal witdth in sigma
 599     // if force is true create a new histogram if it doesn't exist allready
 600     //
 601     TObjArray *arr = &fCalRocArraySigma;
 602     return GetCalRoc(sector, arr, force);
 603 }
 604 //_____________________________________________________________________
 605 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocMean(Int_t sector, Bool_t force)
 606 {
 607   //
 608     // return pointer to ROC with signal mean information
 609     // if force is true create a new histogram if it doesn't exist allready
 610   //
 611   TObjArray *arr = &fCalRocArrayMean;
 612   return GetCalRoc(sector, arr, force);
 613 }
 614
 615 //_____________________________________________________________________
 616 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocRMS(Int_t sector, Bool_t force)
 617 {
 618   //
 619     // return pointer to signal width ROC Calibration
 620     // if force is true create a new histogram if it doesn't exist allready
 621   //
 622   TObjArray *arr = &fCalRocArrayRMS;
 623   return GetCalRoc(sector, arr, force);
 624 }
 625
 626
 627 //_____________________________________________________________________
 628 void AliTPCCalibPedestal::Merge(AliTPCCalibPedestal *ped)
 629 {
 630   //
 631   //  Merge reference histograms of sig to the current AliTPCCalibSignal
 632   //
 633
 634   // merge histograms
 635   for (Int_t iSec=0; iSec<72; ++iSec){
 636     TH2F *hRefPedMerge   = ped->GetHistoPedestal(iSec);
 637
 638     if ( hRefPedMerge ){
 639       TDirectory *dir = hRefPedMerge->GetDirectory(); hRefPedMerge->SetDirectory(0);
 640       TH2F *hRefPed   = GetHistoPedestal(iSec);
 641       if ( hRefPed ) hRefPed->Add(hRefPedMerge);
 642       else {
 643         TH2F *hist = new TH2F(*hRefPedMerge);
 644         hist->SetDirectory(0);
 645         fHistoPedestalArray.AddAt(hist, iSec);
 646       }
 647       hRefPedMerge->SetDirectory(dir);
 648     }
 649   }
 650
 651   // merge array
 652   // ...
 653
 654 }
 655
 656
 657 //_____________________________________________________________________
 658 void AliTPCCalibPedestal::Analyse()
 659 {
 660   //
 661   //  Calculate calibration constants
 662   //
 663
 664   Int_t nbinsAdc = fAdcMax-fAdcMin;
 665
 666   TVectorD param(4);
 667   TMatrixD dummy(3,3);
 668
 669   TH1F *hChannel=new TH1F("hChannel","hChannel",nbinsAdc,fAdcMin,fAdcMax);
 670
 671   Float_t *array_hP=0;
 672
 673   for (Int_t iSec=0; iSec<72; ++iSec){
 674     TH2F *hP = GetHistoPedestal(iSec);
 675     if ( !hP ) continue;
 676
 677     AliTPCCalROC *rocPedestal = GetCalRocPedestal(iSec,kTRUE);
 678     AliTPCCalROC *rocSigma    = GetCalRocSigma(iSec,kTRUE);
 679     AliTPCCalROC *rocMean     = GetCalRocMean(iSec,kTRUE);
 680     AliTPCCalROC *rocRMS      = GetCalRocRMS(iSec,kTRUE);
 681
 682     array_hP = hP->GetArray();
 683     UInt_t nChannels = fROC->GetNChannels(iSec);
 684
 685     for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
 686       Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
 687       //calculate mean and sigma using a gaus fit
 688       //Double_t ret =
 689       AliMathBase::FitGaus(array_hP+offset,nbinsAdc,fAdcMin,fAdcMax,&param,&dummy);
 690       // if the fitting failed set noise and pedestal to 0
 691       // is now done in AliMathBase::FitGaus !
 692 //       if ( ret == -4 ) {
 693 //      param[1]=0;
 694 //      param[2]=0;
 695 //       }
 696       if ( param[1]<fAdcMin || param[1]>fAdcMax ){
 697         param[1]=0;
 698         param[2]=0;
 699       }
 700       rocPedestal->SetValue(iChannel,param[1]);
 701       rocSigma->SetValue(iChannel,param[2]);
 702       //calculate mean and RMS using a truncated means
 703       hChannel->Set(nbinsAdc+2,array_hP+offset-1);
 704       hChannel->SetEntries(param[3]);
 705       param[1]=0;
 706       param[2]=0;
 707       if ( param[3]>0 ) AliMathBase::TruncatedMean(hChannel,&param,fAnaMeanDown,fAnaMeanUp);
 708       rocMean->SetValue(iChannel,param[1]);
 709       rocRMS->SetValue(iChannel,param[2]);
 710     }
 711   }
 712   delete hChannel;
 713 }
 714
 715
 716 //_____________________________________________________________________
 717 void AliTPCCalibPedestal::AnalyseTime(Int_t nevents)
 718 {
 719   //
 720   // Calculate for each channel and time bin the mean pedestal. This
 721   // is used on LDC by TPCPEDESTALda to generate data used for configuration
 722   // of the pattern memory for baseline subtraction in the ALTROs.
 723   //
 724
 725   if ( nevents <= 0 ) return;
 726   if ( fTimeAnalysis ) {
 727     for (Int_t i = 0; i < 159; i++) {  // padrows
 728       for (Int_t j = 0; j < 140; j++) {  // pads per row
 729         for (Int_t k = 0; k < 1024; k++) {  // time bins per pad
 730           fTimeSignal[i][j].AddAt(fTimeSignal[i][j].At(k)/(Float_t)nevents, k);
 731         }
 732       }
 733     }
 734   }
 735 }
 736
 737
 738 //_____________________________________________________________________
 739 void AliTPCCalibPedestal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
 740 {
 741   //
 742   //  Write class to file
 743   //
 744
 745   TString sDir(dir);
 746   TString option;
 747
 748   if ( append )
 749     option = "update";
 750   else
 751     option = "recreate";
 752
 753   TDirectory *backup = gDirectory;
 754   TFile f(filename,option.Data());
 755   f.cd();
 756   if ( !sDir.IsNull() ){
 757     f.mkdir(sDir.Data());
 758     f.cd(sDir);
 759   }
 760   this->Write();
 761   f.Close();
 762
 763   if ( backup ) backup->cd();
 764 }