TPC/AliTPCCalibPulser.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 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////////////////////////////////
  19 //                                                                                     //
  20 //                  Implementation of the TPC pulser calibration                       //
  21 //                                                                                     //
  22 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch   //
  23 //                                                                                     //
  24 /////////////////////////////////////////////////////////////////////////////////////////
  25 /***************************************************************************
  26  *                      Class Description                                  *
  27  ***************************************************************************
  28
  29  The AliTPCCalibPulser class is used to get calibration data concerning the FEE using
  30  runs performed with the calibration pulser.
  31
  32  The information retrieved is
  33  - Time0 differences
  34  - Signal width differences
  35  - Amplification variations
  36
  37  the seen differences arise from the manufacturing tolerances of the PASAs and are very small within
  38  one chip and somewhat large between different chips.
  39
  40  Histograms:
  41    For each ROC three TH2S histos 'Reference Histograms'  (ROC channel vs. [Time0, signal width, Q sum]) is created when
  42    it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are stored in the
  43    TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.
  44
  45
  46  Working principle:
  47  ------------------
  48  Raw calibration pulser data is processed by calling one of the ProcessEvent(...) functions
  49  (see below). These in the end call the Update(...) function.
  50
  51  - the Update(...) function:
  52    In this function the array fPadSignal is filled with the adc signals between the specified range
  53    fFirstTimeBin and fLastTimeBin for the current pad.
  54    before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
  55    stored in fPadSignal.
  56
  57    - the ProcessPad() function:
  58      Find Pedestal and Noise information
  59      - use database information which has to be set by calling
  60        SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)
  61      - if no information from the pedestal data base
  62        is available the informaion is calculated on the fly ( see FindPedestal() function )
  63
  64      Find the Pulser signal information
  65      - calculate  mean = T0, RMS = signal width and Q sum in a range of -2+7 timebins around Q max
  66        the Q sum is scaled by pad area
  67        (see FindPulserSignal(...) function)
  68
  69      Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)
  70      Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE)),
  71
  72  At the end of each event the EndEvent() function is called
  73
  74  - the EndEvent() function:
  75    calculate the mean T0 for each ROC and fill the Time0 histogram with Time0-<Time0 for ROC>
  76    This is done to overcome syncronisation problems between the trigger and the fec clock.
  77
  78  After accumulating the desired statistics the Analyse() function has to be called.
  79  - the Analyse() function
  80    Whithin this function the mean values of T0, RMS, Q are calculated for each pad, using
  81    the AliMathBase::GetCOG(...) function, and the calibration
  82    storage classes (AliTPCCalROC) are filled for each ROC.
  83    The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and
  84    fCalRocArrayQ;
  85
  86
  87
  88  User interface for filling data:
  89  --------------------------------
  90
  91  To Fill information one of the following functions can be used:
  92
  93  Bool_t ProcessEvent(eventHeaderStruct *event);
  94    - process Date event
  95    - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
  96
  97  Bool_t ProcessEvent(AliRawReader *rawReader);
  98    - process AliRawReader event
  99    - use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)
 100
 101  Bool_t ProcessEvent(AliTPCRawStream *rawStream);
 102    - process event from AliTPCRawStream
 103    - call Update function for signal filling
 104
 105  Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
 106               iPad,  const Int_t iTimeBin, const Float_t signal);
 107    - directly  fill signal information (sector, row, pad, time bin, pad)
 108      to the reference histograms
 109
 110  It is also possible to merge two independently taken calibrations using the function
 111
 112  void Merge(AliTPCCalibPulser *sig)
 113    - copy histograms in 'sig' if the do not exist in this instance
 114    - Add histograms in 'sig' to the histograms in this instance if the allready exist
 115    - After merging call Analyse again!
 116
 117
 118
 119  -- example: filling data using root raw data:
 120  void fillSignal(Char_t *filename)
 121  {
 122     rawReader = new AliRawReaderRoot(fileName);
 123     if ( !rawReader ) return;
 124     AliTPCCalibPulser *calib = new AliTPCCalibPulser;
 125     while (rawReader->NextEvent()){
 126       calib->ProcessEvent(rawReader);
 127     }
 128     calib->Analyse();
 129     calib->DumpToFile("SignalData.root");
 130     delete rawReader;
 131     delete calib;
 132  }
 133
 134
 135  What kind of information is stored and how to retrieve them:
 136  ------------------------------------------------------------
 137
 138  - Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):
 139
 140    TH2F *GetHistoT0(Int_t sector);
 141    TH2F *GetHistoRMS(Int_t sector);
 142    TH2F *GetHistoQ(Int_t sector);
 143
 144  - Accessing the calibration storage objects:
 145
 146    AliTPCCalROC *GetCalRocT0(Int_t sector);   // for the Time0 values
 147    AliTPCCalROC *GetCalRocRMS(Int_t sector);  // for the signal width values
 148    AliTPCCalROC *GetCalRocQ(Int_t sector);    // for the Q sum values
 149
 150    example for visualisation:
 151    if the file "SignalData.root" was created using the above example one could do the following:
 152
 153    TFile fileSignal("SignalData.root")
 154    AliTPCCalibPulser *sig = (AliTPCCalibPulser*)fileSignal->Get("AliTPCCalibPulser");
 155    sig->GetCalRocT0(0)->Draw("colz");
 156    sig->GetCalRocRMS(0)->Draw("colz");
 157
 158    or use the AliTPCCalPad functionality:
 159    AliTPCCalPad padT0(ped->GetCalPadT0());
 160    AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
 161    padT0->MakeHisto2D()->Draw("colz");       //Draw A-Side Time0 Information
 162    padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
 163 */
 164
 165
 166 //Root includes
 167 #include <TObjArray.h>
 168 #include <TH1F.h>
 169 #include <TH2S.h>
 170 #include <TString.h>
 171 #include <TVectorF.h>
 172 #include <TMath.h>
 173
 174 #include <TDirectory.h>
 175 #include <TSystem.h>
 176 #include <TROOT.h>
 177 #include <TFile.h>
 178
 179 //AliRoot includes
 180 #include "AliRawReader.h"
 181 #include "AliRawReaderRoot.h"
 182 #include "AliRawReaderDate.h"
 183 #include "AliTPCRawStream.h"
 184 #include "AliTPCRawStreamFast.h"
 185 #include "AliTPCCalROC.h"
 186 #include "AliTPCCalPad.h"
 187 #include "AliTPCROC.h"
 188 #include "AliTPCParam.h"
 189 #include "AliTPCCalibPulser.h"
 190 #include "AliTPCcalibDB.h"
 191 #include "AliMathBase.h"
 192 #include "AliLog.h"
 193 #include "TTreeStream.h"
 194
 195 //date
 196 #include "event.h"
 197
 198
 199
 200
 201 ClassImp(AliTPCCalibPulser)
 202
 203 AliTPCCalibPulser::AliTPCCalibPulser() :
 204     TObject(),
 205     fFirstTimeBin(60),
 206     fLastTimeBin(120),
 207     fNbinsT0(200),
 208     fXminT0(-2),
 209     fXmaxT0(2),
 210     fNbinsQ(200),
 211     fXminQ(1),
 212     fXmaxQ(40),
 213     fNbinsRMS(100),
 214     fXminRMS(0.1),
 215     fXmaxRMS(5.1),
 216     fLastSector(-1),
 217     fOldRCUformat(kTRUE),
 218     fROC(AliTPCROC::Instance()),
 219     fMapping(NULL),
 220     fParam(new AliTPCParam),
 221     fPedestalTPC(0x0),
 222     fPadNoiseTPC(0x0),
 223     fOutliers(0x0),
 224     fPedestalROC(0x0),
 225     fPadNoiseROC(0x0),
 226     fCalRocArrayT0(72),
 227     fCalRocArrayQ(72),
 228     fCalRocArrayRMS(72),
 229     fCalRocArrayOutliers(72),
 230     fHistoQArray(72),
 231     fHistoT0Array(72),
 232     fHistoRMSArray(72),
 233     fPadTimesArrayEvent(72),
 234     fPadQArrayEvent(72),
 235     fPadRMSArrayEvent(72),
 236     fPadPedestalArrayEvent(72),
 237     fCurrentChannel(-1),
 238     fCurrentSector(-1),
 239     fCurrentRow(-1),
 240     fMaxPadSignal(-1),
 241     fMaxTimeBin(-1),
 242     fPadSignal(1024),
 243     fPadPedestal(0),
 244     fPadNoise(0),
 245     fVTime0Offset(72),
 246     fVTime0OffsetCounter(72),
 247 //    fEvent(-1),
 248     fDebugStreamer(0x0),
 249     fDebugLevel(0)
 250 {
 251     //
 252     // AliTPCSignal default constructor
 253     //
 254
 255 }
 256 //_____________________________________________________________________
 257 AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
 258     TObject(sig),
 259     fFirstTimeBin(sig.fFirstTimeBin),
 260     fLastTimeBin(sig.fLastTimeBin),
 261     fNbinsT0(sig.fNbinsT0),
 262     fXminT0(sig.fXminT0),
 263     fXmaxT0(sig.fXmaxT0),
 264     fNbinsQ(sig.fNbinsQ),
 265     fXminQ(sig.fXminQ),
 266     fXmaxQ(sig.fXmaxQ),
 267     fNbinsRMS(sig.fNbinsRMS),
 268     fXminRMS(sig.fXminRMS),
 269     fXmaxRMS(sig.fXmaxRMS),
 270     fLastSector(-1),
 271     fOldRCUformat(kTRUE),
 272     fROC(AliTPCROC::Instance()),
 273     fMapping(NULL),
 274     fParam(new AliTPCParam),
 275     fPedestalTPC(0x0),
 276     fPadNoiseTPC(0x0),
 277     fOutliers(0x0),
 278     fPedestalROC(0x0),
 279     fPadNoiseROC(0x0),
 280     fCalRocArrayT0(72),
 281     fCalRocArrayQ(72),
 282     fCalRocArrayRMS(72),
 283     fCalRocArrayOutliers(72),
 284     fHistoQArray(72),
 285     fHistoT0Array(72),
 286     fHistoRMSArray(72),
 287     fPadTimesArrayEvent(72),
 288     fPadQArrayEvent(72),
 289     fPadRMSArrayEvent(72),
 290     fPadPedestalArrayEvent(72),
 291     fCurrentChannel(-1),
 292     fCurrentSector(-1),
 293     fCurrentRow(-1),
 294     fMaxPadSignal(-1),
 295     fMaxTimeBin(-1),
 296     fPadSignal(1024),
 297     fPadPedestal(0),
 298     fPadNoise(0),
 299     fVTime0Offset(72),
 300     fVTime0OffsetCounter(72),
 301 //    fEvent(-1),
 302     fDebugStreamer(0x0),
 303     fDebugLevel(sig.fDebugLevel)
 304 {
 305     //
 306     // AliTPCSignal default constructor
 307     //
 308
 309     for (Int_t iSec = 0; iSec < 72; ++iSec){
 310         const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
 311         const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
 312         const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
 313         const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
 314
 315         const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
 316         const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
 317         const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
 318
 319         if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
 320         if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
 321         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
 322         if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
 323
 324         if ( hQ != 0x0 ){
 325             TH2S *hNew = new TH2S(*hQ);
 326             hNew->SetDirectory(0);
 327             fHistoQArray.AddAt(hNew,iSec);
 328         }
 329         if ( hT0 != 0x0 ){
 330             TH2S *hNew = new TH2S(*hT0);
 331             hNew->SetDirectory(0);
 332             fHistoQArray.AddAt(hNew,iSec);
 333         }
 334         if ( hRMS != 0x0 ){
 335             TH2S *hNew = new TH2S(*hRMS);
 336             hNew->SetDirectory(0);
 337             fHistoQArray.AddAt(hNew,iSec);
 338         }
 339     }
 340
 341 }
 342 //_____________________________________________________________________
 343 AliTPCCalibPulser& AliTPCCalibPulser::operator = (const  AliTPCCalibPulser &source)
 344 {
 345   //
 346   // assignment operator
 347   //
 348   if (&source == this) return *this;
 349   new (this) AliTPCCalibPulser(source);
 350
 351   return *this;
 352 }
 353 //_____________________________________________________________________
 354 AliTPCCalibPulser::~AliTPCCalibPulser()
 355 {
 356     //
 357     // destructor
 358     //
 359
 360     Reset();
 361
 362     if ( fDebugStreamer) delete fDebugStreamer;
 363     delete fROC;
 364     delete fParam;
 365 }
 366 void AliTPCCalibPulser::Reset()
 367 {
 368     //
 369     // Delete all information: Arrays, Histograms, CalRoc objects
 370     //
 371     fCalRocArrayT0.Delete();
 372     fCalRocArrayQ.Delete();
 373     fCalRocArrayRMS.Delete();
 374     fCalRocArrayOutliers.Delete();
 375
 376     fHistoQArray.Delete();
 377     fHistoT0Array.Delete();
 378     fHistoRMSArray.Delete();
 379
 380     fPadTimesArrayEvent.Delete();
 381     fPadQArrayEvent.Delete();
 382     fPadRMSArrayEvent.Delete();
 383     fPadPedestalArrayEvent.Delete();
 384 }
 385 //_____________________________________________________________________
 386 Int_t AliTPCCalibPulser::Update(const Int_t icsector,
 387                                 const Int_t icRow,
 388                                 const Int_t icPad,
 389                                 const Int_t icTimeBin,
 390                                 const Float_t csignal)
 391 {
 392     //
 393     // Signal filling methode on the fly pedestal and time offset correction if necessary.
 394     // no extra analysis necessary. Assumes knowledge of the signal shape!
 395     // assumes that it is looped over consecutive time bins of one pad
 396     //
 397
 398    if (icRow<0) return 0;
 399    if (icPad<0) return 0;
 400    if (icTimeBin<0) return 0;
 401     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
 402
 403     if ( icRow<0 || icPad<0 ){
 404         AliWarning("Wrong Pad or Row number, skipping!");
 405         return 0;
 406     }
 407
 408     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 409
 410     //init first pad and sector in this event
 411     if ( fCurrentChannel == -1 ) {
 412         fCurrentChannel = iChannel;
 413         fCurrentSector  = icsector;
 414         fCurrentRow     = icRow;
 415     }
 416
 417     //process last pad if we change to a new one
 418     if ( iChannel != fCurrentChannel ){
 419         ProcessPad();
 420         fCurrentChannel = iChannel;
 421         fCurrentSector  = icsector;
 422         fCurrentRow     = icRow;
 423     }
 424
 425     //fill signals for current pad
 426     fPadSignal[icTimeBin]=csignal;
 427     if ( csignal > fMaxPadSignal ){
 428         fMaxPadSignal = csignal;
 429         fMaxTimeBin   = icTimeBin;
 430     }
 431     return 0;
 432 }
 433 //_____________________________________________________________________
 434 void AliTPCCalibPulser::FindPedestal(Float_t part)
 435 {
 436     //
 437     // find pedestal and noise for the current pad. Use either database or
 438     // truncated mean with part*100%
 439     //
 440     Bool_t noPedestal = kTRUE;;
 441     if (fPedestalTPC&&fPadNoiseTPC){
 442         //use pedestal database
 443         //only load new pedestals if the sector has changed
 444         if ( fCurrentSector!=fLastSector ){
 445             fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
 446             fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
 447             fLastSector=fCurrentSector;
 448         }
 449
 450         if ( fPedestalROC&&fPadNoiseROC ){
 451             fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
 452             fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
 453             noPedestal   = kFALSE;
 454         }
 455
 456     }
 457
 458     //if we are not running with pedestal database, or for the current sector there is no information
 459     //available, calculate the pedestal and noise on the fly
 460     if ( noPedestal ) {
 461         const Int_t kPedMax = 100;  //maximum pedestal value
 462         Float_t  max    =  0;
 463         Float_t  maxPos =  0;
 464         Int_t    median =  -1;
 465         Int_t    count0 =  0;
 466         Int_t    count1 =  0;
 467         //
 468         Float_t padSignal=0;
 469         //
 470         UShort_t histo[kPedMax];
 471         memset(histo,0,kPedMax*sizeof(UShort_t));
 472
 473         for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
 474             padSignal = fPadSignal.GetMatrixArray()[i];
 475             if (padSignal<=0) continue;
 476             if (padSignal>max && i>10) {
 477                 max = padSignal;
 478                 maxPos = i;
 479             }
 480             if (padSignal>kPedMax-1) continue;
 481             histo[Int_t(padSignal+0.5)]++;
 482             count0++;
 483         }
 484             //
 485         for (Int_t i=1; i<kPedMax; ++i){
 486             if (count1<count0*0.5) median=i;
 487             count1+=histo[i];
 488         }
 489         // truncated mean
 490         //
 491         // what if by chance histo[median] == 0 ?!?
 492         Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
 493         //
 494         for (Int_t idelta=1; idelta<10; ++idelta){
 495             if (median-idelta<=0) continue;
 496             if (median+idelta>kPedMax) continue;
 497             if (count<part*count1){
 498                 count+=histo[median-idelta];
 499                 mean +=histo[median-idelta]*(median-idelta);
 500                 rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
 501                 count+=histo[median+idelta];
 502                 mean +=histo[median+idelta]*(median+idelta);
 503                 rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
 504             }
 505         }
 506         fPadPedestal = 0;
 507         fPadNoise    = 0;
 508         if ( count > 0 ) {
 509             mean/=count;
 510             rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
 511             fPadPedestal = mean;
 512             fPadNoise    = rms;
 513         }
 514     }
 515 }
 516 //_____________________________________________________________________
 517 void AliTPCCalibPulser::FindPulserSignal(TVectorD &param, Float_t &qSum)
 518 {
 519     //
 520     //  Find position, signal width and height of the CE signal (last signal)
 521     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
 522     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
 523     //
 524
 525     Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
 526     Int_t   cemaxpos       = fMaxTimeBin;
 527     Float_t ceSumThreshold = 10.*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal sum
 528     Float_t ceMaxThreshold = 5.*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal max
 529     const Int_t    kCemin  = 2;             // range for the analysis of the ce signal +- channels from the peak
 530     const Int_t    kCemax  = 7;
 531     param[0] = ceQmax;
 532     param[1] = ceTime;
 533     param[2] = ceRMS;
 534     qSum     = ceQsum;
 535
 536     if (cemaxpos>0){
 537         ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
 538         if ( ceQmax<ceMaxThreshold ) return;
 539         for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
 540             Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
 541             if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
 542                 ceTime+=signal*(i+0.5);
 543                 ceRMS +=signal*(i+0.5)*(i+0.5);
 544                 ceQsum+=signal;
 545             }
 546         }
 547     }
 548     if (ceQsum>ceSumThreshold) {
 549         ceTime/=ceQsum;
 550         ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
 551         //only fill the Time0Offset if pad was not marked as an outlier!
 552         if ( !fOutliers ){
 553             fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
 554             fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
 555         } else {
 556             if ( !(fOutliers->GetCalROC(fCurrentSector)->GetValue(fCurrentChannel)) ){
 557                 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
 558                 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
 559             }
 560         }
 561
 562         //Normalise Q to the pad area
 563         Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
 564
 565         ceQsum/=norm;
 566     } else {
 567         ceQmax=0;
 568         ceTime=0;
 569         ceRMS =0;
 570         ceQsum=0;
 571     }
 572     param[0] = ceQmax;
 573     param[1] = ceTime;
 574     param[2] = ceRMS;
 575     qSum     = ceQsum;
 576 }
 577 //_____________________________________________________________________
 578 void AliTPCCalibPulser::ProcessPad()
 579 {
 580     //
 581     //  Process data of current pad
 582     //
 583
 584     FindPedestal();
 585     TVectorD param(3);
 586     Float_t  qSum;
 587     FindPulserSignal(param, qSum);
 588
 589     Double_t meanT  = param[1];
 590     Double_t sigmaT = param[2];
 591
 592
 593     //Fill Event T0 counter
 594     (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
 595
 596     //Fill Q histogram
 597     GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
 598
 599     //Fill RMS histogram
 600     GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
 601
 602
 603     //Fill debugging info
 604     if ( fDebugLevel>0 ){
 605         (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
 606         (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
 607         (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
 608     }
 609     ResetPad();
 610 }
 611 //_____________________________________________________________________
 612 void AliTPCCalibPulser::EndEvent()
 613 {
 614     //
 615     //  Process data of current event
 616     //
 617
 618     //check if last pad has allready been processed, if not do so
 619     if ( fMaxTimeBin>-1 ) ProcessPad();
 620
 621     //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
 622     for ( Int_t iSec = 0; iSec<72; ++iSec ){
 623         TVectorF *vTimes = GetPadTimesEvent(iSec);
 624         if ( !vTimes || fVTime0OffsetCounter[iSec]==0 ) continue;
 625         Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
 626         for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
 627             Float_t time  = (*vTimes).GetMatrixArray()[iChannel];
 628
 629             GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
 630             //GetHistoT0(iSec,kTRUE)->Fill( time,iChannel );
 631
 632
 633             //Debug start
 634             if ( fDebugLevel>0 ){
 635                 if ( !fDebugStreamer ) {
 636                         //debug stream
 637                     TDirectory *backup = gDirectory;
 638                     fDebugStreamer = new TTreeSRedirector("deb2.root");
 639                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 640                 }
 641
 642                 Int_t row=0;
 643                 Int_t pad=0;
 644                 Int_t padc=0;
 645
 646                 Float_t q   = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
 647                 Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
 648
 649                 UInt_t channel=iChannel;
 650                 Int_t sector=iSec;
 651
 652                 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
 653                 pad = channel-fROC->GetRowIndexes(sector)[row];
 654                 padc = pad-(fROC->GetNPads(sector,row)/2);
 655
 656                 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
 657                                     Form("hSignalD%d.%d.%d",sector,row,pad),
 658                                     fLastTimeBin-fFirstTimeBin,
 659                                     fFirstTimeBin,fLastTimeBin);
 660                 h1->SetDirectory(0);
 661
 662                 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
 663                     h1->Fill(i,fPadSignal(i));
 664
 665                 (*fDebugStreamer) << "DataPad" <<
 666 //                  "Event=" << fEvent <<
 667                     "Sector="<< sector <<
 668                     "Row="   << row<<
 669                     "Pad="   << pad <<
 670                     "PadC="  << padc <<
 671                     "PadSec="<< channel <<
 672                     "Time0="  << time0 <<
 673                     "Time="  << time <<
 674                     "RMS="   << rms <<
 675                     "Sum="   << q <<
 676                     "hist.=" << h1 <<
 677                     "\n";
 678
 679                 delete h1;
 680             }
 681             //Debug end
 682         }
 683     }
 684 }
 685 //_____________________________________________________________________
 686 Bool_t AliTPCCalibPulser::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
 687 {
 688   //
 689   // Event Processing loop - AliTPCRawStream
 690   //
 691   ResetEvent();
 692
 693   Bool_t withInput = kFALSE;
 694
 695   while ( rawStreamFast->NextDDL() ){
 696       while ( rawStreamFast->NextChannel() ){
 697           Int_t isector  = rawStreamFast->GetSector();                       //  current sector
 698           Int_t iRow     = rawStreamFast->GetRow();                          //  current row
 699           Int_t iPad     = rawStreamFast->GetPad();                          //  current pad
 700
 701           while ( rawStreamFast->NextBunch() ){
 702               Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
 703               Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
 704               for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
 705                   Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
 706                   Update(isector,iRow,iPad,iTimeBin+1,signal);
 707                   withInput = kTRUE;
 708               }
 709           }
 710       }
 711   }
 712   if (withInput){
 713       EndEvent();
 714   }
 715   return withInput;
 716 }
 717 //_____________________________________________________________________
 718 Bool_t AliTPCCalibPulser::ProcessEventFast(AliRawReader *rawReader)
 719 {
 720   //
 721   //  Event processing loop - AliRawReader
 722   //
 723   AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
 724   Bool_t res=ProcessEventFast(rawStreamFast);
 725   delete rawStreamFast;
 726   return res;
 727 }
 728 //_____________________________________________________________________
 729 Bool_t AliTPCCalibPulser::ProcessEvent(AliTPCRawStream *rawStream)
 730 {
 731   //
 732   // Event Processing loop - AliTPCRawStream
 733   //
 734
 735   rawStream->SetOldRCUFormat(fOldRCUformat);
 736
 737   ResetEvent();
 738
 739   Bool_t withInput = kFALSE;
 740
 741   while (rawStream->Next()) {
 742       Int_t isector  = rawStream->GetSector();                       //  current sector
 743       Int_t iRow     = rawStream->GetRow();                          //  current row
 744       Int_t iPad     = rawStream->GetPad();                          //  current pad
 745       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
 746       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
 747
 748       Update(isector,iRow,iPad,iTimeBin,signal);
 749       withInput = kTRUE;
 750   }
 751   if (withInput){
 752       EndEvent();
 753   }
 754   return withInput;
 755 }
 756 //_____________________________________________________________________
 757 Bool_t AliTPCCalibPulser::ProcessEvent(AliRawReader *rawReader)
 758 {
 759   //
 760   //  Event processing loop - AliRawReader
 761   //
 762
 763
 764   AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
 765
 766   rawReader->Select("TPC");
 767
 768   return ProcessEvent(&rawStream);
 769 }
 770 //_____________________________________________________________________
 771 Bool_t AliTPCCalibPulser::ProcessEvent(eventHeaderStruct *event)
 772 {
 773   //
 774   //  Event processing loop - date event
 775   //
 776     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 777     Bool_t result=ProcessEvent(rawReader);
 778     delete rawReader;
 779     return result;
 780
 781 }
 782 //_____________________________________________________________________
 783 TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
 784                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
 785                                   Char_t *type, Bool_t force)
 786 {
 787     //
 788     // return pointer to Q histogram
 789     // if force is true create a new histogram if it doesn't exist allready
 790     //
 791     if ( !force || arr->UncheckedAt(sector) )
 792         return (TH2S*)arr->UncheckedAt(sector);
 793
 794     // if we are forced and histogram doesn't yes exist create it
 795     Char_t name[255], title[255];
 796
 797     sprintf(name,"hCalib%s%.2d",type,sector);
 798     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
 799
 800     // new histogram with Q calib information. One value for each pad!
 801     TH2S* hist = new TH2S(name,title,
 802                           nbinsY, ymin, ymax,
 803                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
 804     hist->SetDirectory(0);
 805     arr->AddAt(hist,sector);
 806     return hist;
 807 }
 808 //_____________________________________________________________________
 809 TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
 810 {
 811     //
 812     // return pointer to T0 histogram
 813     // if force is true create a new histogram if it doesn't exist allready
 814     //
 815     TObjArray *arr = &fHistoT0Array;
 816     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
 817 }
 818 //_____________________________________________________________________
 819 TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
 820 {
 821     //
 822     // return pointer to Q histogram
 823     // if force is true create a new histogram if it doesn't exist allready
 824     //
 825     TObjArray *arr = &fHistoQArray;
 826     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
 827 }
 828 //_____________________________________________________________________
 829 TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
 830 {
 831     //
 832     // return pointer to Q histogram
 833     // if force is true create a new histogram if it doesn't exist allready
 834     //
 835     TObjArray *arr = &fHistoRMSArray;
 836     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
 837 }
 838 //_____________________________________________________________________
 839 TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
 840 {
 841     //
 842     // return pointer to Pad Info from 'arr' for the current event and sector
 843     // if force is true create it if it doesn't exist allready
 844     //
 845     if ( !force || arr->UncheckedAt(sector) )
 846         return (TVectorF*)arr->UncheckedAt(sector);
 847
 848     TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
 849     arr->AddAt(vect,sector);
 850     return vect;
 851 }
 852 //_____________________________________________________________________
 853 TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
 854 {
 855     //
 856     // return pointer to Pad Times Array for the current event and sector
 857     // if force is true create it if it doesn't exist allready
 858     //
 859     TObjArray *arr = &fPadTimesArrayEvent;
 860     return GetPadInfoEvent(sector,arr,force);
 861 }
 862 //_____________________________________________________________________
 863 TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
 864 {
 865     //
 866     // return pointer to Pad Q Array for the current event and sector
 867     // if force is true create it if it doesn't exist allready
 868     // for debugging purposes only
 869     //
 870
 871     TObjArray *arr = &fPadQArrayEvent;
 872     return GetPadInfoEvent(sector,arr,force);
 873 }
 874 //_____________________________________________________________________
 875 TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
 876 {
 877     //
 878     // return pointer to Pad RMS Array for the current event and sector
 879     // if force is true create it if it doesn't exist allready
 880     // for debugging purposes only
 881     //
 882     TObjArray *arr = &fPadRMSArrayEvent;
 883     return GetPadInfoEvent(sector,arr,force);
 884 }
 885 //_____________________________________________________________________
 886 TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
 887 {
 888     //
 889     // return pointer to Pad RMS Array for the current event and sector
 890     // if force is true create it if it doesn't exist allready
 891     // for debugging purposes only
 892     //
 893     TObjArray *arr = &fPadPedestalArrayEvent;
 894     return GetPadInfoEvent(sector,arr,force);
 895 }
 896 //_____________________________________________________________________
 897 AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
 898 {
 899     //
 900     // return pointer to ROC Calibration
 901     // if force is true create a new histogram if it doesn't exist allready
 902     //
 903     if ( !force || arr->UncheckedAt(sector) )
 904         return (AliTPCCalROC*)arr->UncheckedAt(sector);
 905
 906     // if we are forced and histogram doesn't yes exist create it
 907
 908     // new AliTPCCalROC for T0 information. One value for each pad!
 909     AliTPCCalROC *croc = new AliTPCCalROC(sector);
 910     arr->AddAt(croc,sector);
 911     return croc;
 912 }
 913 //_____________________________________________________________________
 914 AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
 915 {
 916     //
 917     // return pointer to Carge ROC Calibration
 918     // if force is true create a new histogram if it doesn't exist allready
 919     //
 920     TObjArray *arr = &fCalRocArrayT0;
 921     return GetCalRoc(sector, arr, force);
 922 }
 923 //_____________________________________________________________________
 924 AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
 925 {
 926     //
 927     // return pointer to T0 ROC Calibration
 928     // if force is true create a new histogram if it doesn't exist allready
 929     //
 930     TObjArray *arr = &fCalRocArrayQ;
 931     return GetCalRoc(sector, arr, force);
 932 }
 933 //_____________________________________________________________________
 934 AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
 935 {
 936     //
 937     // return pointer to signal width ROC Calibration
 938     // if force is true create a new histogram if it doesn't exist allready
 939     //
 940     TObjArray *arr = &fCalRocArrayRMS;
 941     return GetCalRoc(sector, arr, force);
 942 }
 943 //_____________________________________________________________________
 944 AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
 945 {
 946     //
 947     // return pointer to Outliers
 948     // if force is true create a new histogram if it doesn't exist allready
 949     //
 950     TObjArray *arr = &fCalRocArrayOutliers;
 951     return GetCalRoc(sector, arr, force);
 952 }
 953 //_____________________________________________________________________
 954 void AliTPCCalibPulser::ResetEvent()
 955 {
 956     //
 957     //  Reset global counters  -- Should be called before each event is processed
 958     //
 959     fLastSector=-1;
 960     fCurrentSector=-1;
 961     fCurrentRow=-1;
 962     fCurrentChannel=-1;
 963
 964     ResetPad();
 965
 966     fPadTimesArrayEvent.Delete();
 967     fPadQArrayEvent.Delete();
 968     fPadRMSArrayEvent.Delete();
 969     fPadPedestalArrayEvent.Delete();
 970
 971     for ( Int_t i=0; i<72; ++i ){
 972         fVTime0Offset[i]=0;
 973         fVTime0OffsetCounter[i]=0;
 974     }
 975 }
 976 //_____________________________________________________________________
 977 void AliTPCCalibPulser::ResetPad()
 978 {
 979     //
 980     //  Reset pad infos -- Should be called after a pad has been processed
 981     //
 982     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
 983         fPadSignal[i] = 0;
 984     fMaxTimeBin = -1;
 985     fMaxPadSignal = -1;
 986     fPadPedestal  = -1;
 987     fPadNoise     = -1;
 988 }
 989 //_____________________________________________________________________
 990 void AliTPCCalibPulser::Merge(AliTPCCalibPulser *sig)
 991 {
 992     //
 993     //  Merge reference histograms of sig to the current AliTPCCalibPulser
 994     //
 995
 996     //merge histograms
 997     for (Int_t iSec=0; iSec<72; ++iSec){
 998         TH2S *hRefQmerge   = sig->GetHistoQ(iSec);
 999         TH2S *hRefT0merge  = sig->GetHistoT0(iSec);
1000         TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
1001
1002
1003         if ( hRefQmerge ){
1004             TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1005             TH2S *hRefQ   = GetHistoQ(iSec);
1006             if ( hRefQ ) hRefQ->Add(hRefQmerge);
1007             else {
1008                 TH2S *hist = new TH2S(*hRefQmerge);
1009                 hist->SetDirectory(0);
1010                 fHistoQArray.AddAt(hist, iSec);
1011             }
1012             hRefQmerge->SetDirectory(dir);
1013         }
1014         if ( hRefT0merge ){
1015             TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1016             TH2S *hRefT0  = GetHistoT0(iSec);
1017             if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1018             else {
1019                 TH2S *hist = new TH2S(*hRefT0merge);
1020                 hist->SetDirectory(0);
1021                 fHistoT0Array.AddAt(hist, iSec);
1022             }
1023             hRefT0merge->SetDirectory(dir);
1024         }
1025         if ( hRefRMSmerge ){
1026             TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1027             TH2S *hRefRMS = GetHistoRMS(iSec);
1028             if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1029             else {
1030                 TH2S *hist = new TH2S(*hRefRMSmerge);
1031                 hist->SetDirectory(0);
1032                 fHistoRMSArray.AddAt(hist, iSec);
1033             }
1034             hRefRMSmerge->SetDirectory(dir);
1035         }
1036
1037     }
1038 }
1039 //_____________________________________________________________________
1040 void AliTPCCalibPulser::Analyse()
1041 {
1042     //
1043     //  Calculate calibration constants
1044     //
1045
1046     TVectorD paramQ(3);
1047     TVectorD paramT0(3);
1048     TVectorD paramRMS(3);
1049     TMatrixD dummy(3,3);
1050
1051     for (Int_t iSec=0; iSec<72; ++iSec){
1052         TH2S *hT0 = GetHistoT0(iSec);
1053         if (!hT0 ) continue;
1054
1055         AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
1056         AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
1057         AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1058         AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1059
1060         TH2S *hQ   = GetHistoQ(iSec);
1061         TH2S *hRMS = GetHistoRMS(iSec);
1062
1063         Short_t *arrayhQ   = hQ->GetArray();
1064         Short_t *arrayhT0  = hT0->GetArray();
1065         Short_t *arrayhRMS = hRMS->GetArray();
1066
1067         UInt_t nChannels = fROC->GetNChannels(iSec);
1068
1069         //debug
1070         Int_t row=0;
1071         Int_t pad=0;
1072         Int_t padc=0;
1073         //! debug
1074
1075         for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1076
1077             Float_t cogTime0 = -1000;
1078             Float_t cogQ     = -1000;
1079             Float_t cogRMS   = -1000;
1080             Float_t cogOut   = 0;
1081
1082             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1083             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1084             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1085 /*
1086             AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1087             AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1088             AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1089             cogQ     = paramQ[1];
1090             cogTime0 = paramT0[1];
1091             cogRMS   = paramRMS[1];
1092 */
1093             cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1094             cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1095             cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1096
1097             /*
1098             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1099                 cogOut = 1;
1100                 cogTime0 = 0;
1101                 cogQ     = 0;
1102                 cogRMS   = 0;
1103             }
1104 */
1105             rocQ->SetValue(iChannel, cogQ*cogQ);
1106             rocT0->SetValue(iChannel, cogTime0);
1107             rocRMS->SetValue(iChannel, cogRMS);
1108             rocOut->SetValue(iChannel, cogOut);
1109
1110
1111             //debug
1112             if ( fDebugLevel > 0 ){
1113                 if ( !fDebugStreamer ) {
1114                         //debug stream
1115                     TDirectory *backup = gDirectory;
1116                     fDebugStreamer = new TTreeSRedirector("deb2.root");
1117                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1118                 }
1119
1120                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1121                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1122                 padc = pad-(fROC->GetNPads(iSec,row)/2);
1123
1124                 (*fDebugStreamer) << "DataEnd" <<
1125                     "Sector="  << iSec      <<
1126                     "Pad="     << pad       <<
1127                     "PadC="    << padc      <<
1128                     "Row="     << row       <<
1129                     "PadSec="  << iChannel   <<
1130                     "Q="       << cogQ      <<
1131                     "T0="      << cogTime0  <<
1132                     "RMS="     << cogRMS    <<
1133                     "\n";
1134             }
1135             //! debug
1136         }
1137     }
1138     delete fDebugStreamer;
1139     fDebugStreamer = 0x0;
1140 }
1141 //_____________________________________________________________________
1142 void AliTPCCalibPulser::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1143 {
1144     //
1145     //  Write class to file
1146     //
1147
1148     TString sDir(dir);
1149     TString option;
1150
1151     if ( append )
1152         option = "update";
1153     else
1154         option = "recreate";
1155
1156     TDirectory *backup = gDirectory;
1157     TFile f(filename,option.Data());
1158     f.cd();
1159     if ( !sDir.IsNull() ){
1160         f.mkdir(sDir.Data());
1161         f.cd(sDir);
1162     }
1163     this->Write();
1164     f.Close();
1165
1166     if ( backup ) backup->cd();
1167 }
1168 //_____________________________________________________________________
1169 //_________________________  Test Functions ___________________________
1170 //_____________________________________________________________________
1171 TObjArray* AliTPCCalibPulser::TestBinning()
1172 {
1173     //
1174     //  Function to test the binning of the reference histograms
1175     //  type: T0, Q or RMS
1176     //  mode: 0 - number of filled bins per channel
1177     //        1 - number of empty bins between filled bins in one ROC
1178     //  returns TObjArray with the test histograms type*2+mode:
1179     //  position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1180
1181
1182     TObjArray *histArray = new TObjArray(6);
1183     const Char_t *type[] = {"T0","Q","RMS"};
1184     Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1185
1186     for (Int_t itype = 0; itype<3; ++itype){
1187         for (Int_t imode=0; imode<2; ++imode){
1188             Int_t icount = itype*2+imode;
1189             histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1190                                       Form("Test Binning of '%s', mode - %d",type[itype],imode),
1191                                       72,0,72),
1192                              icount);
1193         }
1194     }
1195
1196
1197     TH2S *hRef=0x0;
1198     Short_t *array=0x0;
1199     for (Int_t itype = 0; itype<3; ++itype){
1200         for (Int_t iSec=0; iSec<72; ++iSec){
1201             if ( itype == 0 ) hRef = GetHistoT0(iSec);
1202             if ( itype == 1 ) hRef = GetHistoQ(iSec);
1203             if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1204             if ( hRef == 0x0 ) continue;
1205             array = (hRef->GetArray());
1206             UInt_t nChannels = fROC->GetNChannels(iSec);
1207
1208             Int_t nempty=0;
1209             for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1210                 Int_t nfilled=0;
1211                 Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1212                 Int_t c1 = 0;
1213                 Int_t c2 = 0;
1214                 for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1215                     if ( array[offset+iBin]>0 ) {
1216                         nfilled++;
1217                         if ( c1 && c2 ) nempty++;
1218                         else c1 = 1;
1219                     }
1220                     else if ( c1 ) c2 = 1;
1221
1222                 }
1223                 ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1224             }
1225             ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);
1226         }
1227     }
1228     return histArray;
1229 }