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