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 ) continue;
 625
 626         for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
 627             Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
 628             Float_t time  = (*vTimes).GetMatrixArray()[iChannel];
 629
 630             GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
 631             //GetHistoT0(iSec,kTRUE)->Fill( time,iChannel );
 632
 633
 634             //Debug start
 635             if ( fDebugLevel>0 ){
 636                 if ( !fDebugStreamer ) {
 637                         //debug stream
 638                     TDirectory *backup = gDirectory;
 639                     fDebugStreamer = new TTreeSRedirector("deb2.root");
 640                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 641                 }
 642
 643                 Int_t row=0;
 644                 Int_t pad=0;
 645                 Int_t padc=0;
 646
 647                 Float_t q   = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
 648                 Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
 649
 650                 UInt_t channel=iChannel;
 651                 Int_t sector=iSec;
 652
 653                 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
 654                 pad = channel-fROC->GetRowIndexes(sector)[row];
 655                 padc = pad-(fROC->GetNPads(sector,row)/2);
 656
 657                 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
 658                                     Form("hSignalD%d.%d.%d",sector,row,pad),
 659                                     fLastTimeBin-fFirstTimeBin,
 660                                     fFirstTimeBin,fLastTimeBin);
 661                 h1->SetDirectory(0);
 662
 663                 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
 664                     h1->Fill(i,fPadSignal(i));
 665
 666                 (*fDebugStreamer) << "DataPad" <<
 667 //                  "Event=" << fEvent <<
 668                     "Sector="<< sector <<
 669                     "Row="   << row<<
 670                     "Pad="   << pad <<
 671                     "PadC="  << padc <<
 672                     "PadSec="<< channel <<
 673                     "Time0="  << time0 <<
 674                     "Time="  << time <<
 675                     "RMS="   << rms <<
 676                     "Sum="   << q <<
 677                     "hist.=" << h1 <<
 678                     "\n";
 679
 680                 delete h1;
 681             }
 682             //Debug end
 683         }
 684     }
 685 }
 686 //_____________________________________________________________________
 687 Bool_t AliTPCCalibPulser::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
 688 {
 689   //
 690   // Event Processing loop - AliTPCRawStream
 691   //
 692   ResetEvent();
 693
 694   Bool_t withInput = kFALSE;
 695
 696   while ( rawStreamFast->NextDDL() ){
 697       while ( rawStreamFast->NextChannel() ){
 698           Int_t isector  = rawStreamFast->GetSector();                       //  current sector
 699           Int_t iRow     = rawStreamFast->GetRow();                          //  current row
 700           Int_t iPad     = rawStreamFast->GetPad();                          //  current pad
 701
 702           while ( rawStreamFast->NextBunch() ){
 703               Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
 704               Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
 705               for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
 706                   Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
 707                   Update(isector,iRow,iPad,iTimeBin+1,signal);
 708                   withInput = kTRUE;
 709               }
 710           }
 711       }
 712   }
 713   if (withInput){
 714       EndEvent();
 715   }
 716   return withInput;
 717 }
 718 //_____________________________________________________________________
 719 Bool_t AliTPCCalibPulser::ProcessEventFast(AliRawReader *rawReader)
 720 {
 721   //
 722   //  Event processing loop - AliRawReader
 723   //
 724   AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
 725   Bool_t res=ProcessEventFast(rawStreamFast);
 726   delete rawStreamFast;
 727   return res;
 728 }
 729 //_____________________________________________________________________
 730 Bool_t AliTPCCalibPulser::ProcessEvent(AliTPCRawStream *rawStream)
 731 {
 732   //
 733   // Event Processing loop - AliTPCRawStream
 734   //
 735
 736   rawStream->SetOldRCUFormat(fOldRCUformat);
 737
 738   ResetEvent();
 739
 740   Bool_t withInput = kFALSE;
 741
 742   while (rawStream->Next()) {
 743       Int_t isector  = rawStream->GetSector();                       //  current sector
 744       Int_t iRow     = rawStream->GetRow();                          //  current row
 745       Int_t iPad     = rawStream->GetPad();                          //  current pad
 746       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
 747       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
 748
 749       Update(isector,iRow,iPad,iTimeBin,signal);
 750       withInput = kTRUE;
 751   }
 752   if (withInput){
 753       EndEvent();
 754   }
 755   return withInput;
 756 }
 757 //_____________________________________________________________________
 758 Bool_t AliTPCCalibPulser::ProcessEvent(AliRawReader *rawReader)
 759 {
 760   //
 761   //  Event processing loop - AliRawReader
 762   //
 763
 764
 765   AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
 766
 767   rawReader->Select("TPC");
 768
 769   return ProcessEvent(&rawStream);
 770 }
 771 //_____________________________________________________________________
 772 Bool_t AliTPCCalibPulser::ProcessEvent(eventHeaderStruct *event)
 773 {
 774   //
 775   //  Event processing loop - date event
 776   //
 777     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 778     Bool_t result=ProcessEvent(rawReader);
 779     delete rawReader;
 780     return result;
 781
 782 }
 783 //_____________________________________________________________________
 784 TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
 785                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
 786                                   Char_t *type, Bool_t force)
 787 {
 788     //
 789     // return pointer to Q histogram
 790     // if force is true create a new histogram if it doesn't exist allready
 791     //
 792     if ( !force || arr->UncheckedAt(sector) )
 793         return (TH2S*)arr->UncheckedAt(sector);
 794
 795     // if we are forced and histogram doesn't yes exist create it
 796     Char_t name[255], title[255];
 797
 798     sprintf(name,"hCalib%s%.2d",type,sector);
 799     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
 800
 801     // new histogram with Q calib information. One value for each pad!
 802     TH2S* hist = new TH2S(name,title,
 803                           nbinsY, ymin, ymax,
 804                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
 805     hist->SetDirectory(0);
 806     arr->AddAt(hist,sector);
 807     return hist;
 808 }
 809 //_____________________________________________________________________
 810 TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
 811 {
 812     //
 813     // return pointer to T0 histogram
 814     // if force is true create a new histogram if it doesn't exist allready
 815     //
 816     TObjArray *arr = &fHistoT0Array;
 817     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
 818 }
 819 //_____________________________________________________________________
 820 TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
 821 {
 822     //
 823     // return pointer to Q histogram
 824     // if force is true create a new histogram if it doesn't exist allready
 825     //
 826     TObjArray *arr = &fHistoQArray;
 827     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
 828 }
 829 //_____________________________________________________________________
 830 TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
 831 {
 832     //
 833     // return pointer to Q histogram
 834     // if force is true create a new histogram if it doesn't exist allready
 835     //
 836     TObjArray *arr = &fHistoRMSArray;
 837     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
 838 }
 839 //_____________________________________________________________________
 840 TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
 841 {
 842     //
 843     // return pointer to Pad Info from 'arr' for the current event and sector
 844     // if force is true create it if it doesn't exist allready
 845     //
 846     if ( !force || arr->UncheckedAt(sector) )
 847         return (TVectorF*)arr->UncheckedAt(sector);
 848
 849     TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
 850     arr->AddAt(vect,sector);
 851     return vect;
 852 }
 853 //_____________________________________________________________________
 854 TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
 855 {
 856     //
 857     // return pointer to Pad Times Array for the current event and sector
 858     // if force is true create it if it doesn't exist allready
 859     //
 860     TObjArray *arr = &fPadTimesArrayEvent;
 861     return GetPadInfoEvent(sector,arr,force);
 862 }
 863 //_____________________________________________________________________
 864 TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
 865 {
 866     //
 867     // return pointer to Pad Q Array for the current event and sector
 868     // if force is true create it if it doesn't exist allready
 869     // for debugging purposes only
 870     //
 871
 872     TObjArray *arr = &fPadQArrayEvent;
 873     return GetPadInfoEvent(sector,arr,force);
 874 }
 875 //_____________________________________________________________________
 876 TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
 877 {
 878     //
 879     // return pointer to Pad RMS Array for the current event and sector
 880     // if force is true create it if it doesn't exist allready
 881     // for debugging purposes only
 882     //
 883     TObjArray *arr = &fPadRMSArrayEvent;
 884     return GetPadInfoEvent(sector,arr,force);
 885 }
 886 //_____________________________________________________________________
 887 TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
 888 {
 889     //
 890     // return pointer to Pad RMS Array for the current event and sector
 891     // if force is true create it if it doesn't exist allready
 892     // for debugging purposes only
 893     //
 894     TObjArray *arr = &fPadPedestalArrayEvent;
 895     return GetPadInfoEvent(sector,arr,force);
 896 }
 897 //_____________________________________________________________________
 898 AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
 899 {
 900     //
 901     // return pointer to ROC Calibration
 902     // if force is true create a new histogram if it doesn't exist allready
 903     //
 904     if ( !force || arr->UncheckedAt(sector) )
 905         return (AliTPCCalROC*)arr->UncheckedAt(sector);
 906
 907     // if we are forced and histogram doesn't yes exist create it
 908
 909     // new AliTPCCalROC for T0 information. One value for each pad!
 910     AliTPCCalROC *croc = new AliTPCCalROC(sector);
 911     arr->AddAt(croc,sector);
 912     return croc;
 913 }
 914 //_____________________________________________________________________
 915 AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
 916 {
 917     //
 918     // return pointer to Carge ROC Calibration
 919     // if force is true create a new histogram if it doesn't exist allready
 920     //
 921     TObjArray *arr = &fCalRocArrayT0;
 922     return GetCalRoc(sector, arr, force);
 923 }
 924 //_____________________________________________________________________
 925 AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
 926 {
 927     //
 928     // return pointer to T0 ROC Calibration
 929     // if force is true create a new histogram if it doesn't exist allready
 930     //
 931     TObjArray *arr = &fCalRocArrayQ;
 932     return GetCalRoc(sector, arr, force);
 933 }
 934 //_____________________________________________________________________
 935 AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
 936 {
 937     //
 938     // return pointer to signal width ROC Calibration
 939     // if force is true create a new histogram if it doesn't exist allready
 940     //
 941     TObjArray *arr = &fCalRocArrayRMS;
 942     return GetCalRoc(sector, arr, force);
 943 }
 944 //_____________________________________________________________________
 945 AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
 946 {
 947     //
 948     // return pointer to Outliers
 949     // if force is true create a new histogram if it doesn't exist allready
 950     //
 951     TObjArray *arr = &fCalRocArrayOutliers;
 952     return GetCalRoc(sector, arr, force);
 953 }
 954 //_____________________________________________________________________
 955 void AliTPCCalibPulser::ResetEvent()
 956 {
 957     //
 958     //  Reset global counters  -- Should be called before each event is processed
 959     //
 960     fLastSector=-1;
 961     fCurrentSector=-1;
 962     fCurrentRow=-1;
 963     fCurrentChannel=-1;
 964
 965     ResetPad();
 966
 967     fPadTimesArrayEvent.Delete();
 968     fPadQArrayEvent.Delete();
 969     fPadRMSArrayEvent.Delete();
 970     fPadPedestalArrayEvent.Delete();
 971
 972     for ( Int_t i=0; i<72; ++i ){
 973         fVTime0Offset[i]=0;
 974         fVTime0OffsetCounter[i]=0;
 975     }
 976 }
 977 //_____________________________________________________________________
 978 void AliTPCCalibPulser::ResetPad()
 979 {
 980     //
 981     //  Reset pad infos -- Should be called after a pad has been processed
 982     //
 983     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
 984         fPadSignal[i] = 0;
 985     fMaxTimeBin = -1;
 986     fMaxPadSignal = -1;
 987     fPadPedestal  = -1;
 988     fPadNoise     = -1;
 989 }
 990 //_____________________________________________________________________
 991 void AliTPCCalibPulser::Merge(AliTPCCalibPulser *sig)
 992 {
 993     //
 994     //  Merge reference histograms of sig to the current AliTPCCalibPulser
 995     //
 996
 997     //merge histograms
 998     for (Int_t iSec=0; iSec<72; ++iSec){
 999         TH2S *hRefQmerge   = sig->GetHistoQ(iSec);
1000         TH2S *hRefT0merge  = sig->GetHistoT0(iSec);
1001         TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
1002
1003
1004         if ( hRefQmerge ){
1005             TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1006             TH2S *hRefQ   = GetHistoQ(iSec);
1007             if ( hRefQ ) hRefQ->Add(hRefQmerge);
1008             else {
1009                 TH2S *hist = new TH2S(*hRefQmerge);
1010                 hist->SetDirectory(0);
1011                 fHistoQArray.AddAt(hist, iSec);
1012             }
1013             hRefQmerge->SetDirectory(dir);
1014         }
1015         if ( hRefT0merge ){
1016             TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1017             TH2S *hRefT0  = GetHistoT0(iSec);
1018             if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1019             else {
1020                 TH2S *hist = new TH2S(*hRefT0merge);
1021                 hist->SetDirectory(0);
1022                 fHistoT0Array.AddAt(hist, iSec);
1023             }
1024             hRefT0merge->SetDirectory(dir);
1025         }
1026         if ( hRefRMSmerge ){
1027             TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1028             TH2S *hRefRMS = GetHistoRMS(iSec);
1029             if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1030             else {
1031                 TH2S *hist = new TH2S(*hRefRMSmerge);
1032                 hist->SetDirectory(0);
1033                 fHistoRMSArray.AddAt(hist, iSec);
1034             }
1035             hRefRMSmerge->SetDirectory(dir);
1036         }
1037
1038     }
1039 }
1040 //_____________________________________________________________________
1041 void AliTPCCalibPulser::Analyse()
1042 {
1043     //
1044     //  Calculate calibration constants
1045     //
1046
1047     TVectorD paramQ(3);
1048     TVectorD paramT0(3);
1049     TVectorD paramRMS(3);
1050     TMatrixD dummy(3,3);
1051
1052     for (Int_t iSec=0; iSec<72; ++iSec){
1053         TH2S *hT0 = GetHistoT0(iSec);
1054         if (!hT0 ) continue;
1055
1056         AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
1057         AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
1058         AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1059         AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1060
1061         TH2S *hQ   = GetHistoQ(iSec);
1062         TH2S *hRMS = GetHistoRMS(iSec);
1063
1064         Short_t *arrayhQ   = hQ->GetArray();
1065         Short_t *arrayhT0  = hT0->GetArray();
1066         Short_t *arrayhRMS = hRMS->GetArray();
1067
1068         UInt_t nChannels = fROC->GetNChannels(iSec);
1069
1070         //debug
1071         Int_t row=0;
1072         Int_t pad=0;
1073         Int_t padc=0;
1074         //! debug
1075
1076         for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1077
1078             Float_t cogTime0 = -1000;
1079             Float_t cogQ     = -1000;
1080             Float_t cogRMS   = -1000;
1081             Float_t cogOut   = 0;
1082
1083             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1084             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1085             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1086 /*
1087             AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1088             AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1089             AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1090             cogQ     = paramQ[1];
1091             cogTime0 = paramT0[1];
1092             cogRMS   = paramRMS[1];
1093 */
1094             cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1095             cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1096             cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1097
1098             /*
1099             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1100                 cogOut = 1;
1101                 cogTime0 = 0;
1102                 cogQ     = 0;
1103                 cogRMS   = 0;
1104             }
1105 */
1106             rocQ->SetValue(iChannel, cogQ*cogQ);
1107             rocT0->SetValue(iChannel, cogTime0);
1108             rocRMS->SetValue(iChannel, cogRMS);
1109             rocOut->SetValue(iChannel, cogOut);
1110
1111
1112             //debug
1113             if ( fDebugLevel > 0 ){
1114                 if ( !fDebugStreamer ) {
1115                         //debug stream
1116                     TDirectory *backup = gDirectory;
1117                     fDebugStreamer = new TTreeSRedirector("deb2.root");
1118                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1119                 }
1120
1121                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1122                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1123                 padc = pad-(fROC->GetNPads(iSec,row)/2);
1124
1125                 (*fDebugStreamer) << "DataEnd" <<
1126                     "Sector="  << iSec      <<
1127                     "Pad="     << pad       <<
1128                     "PadC="    << padc      <<
1129                     "Row="     << row       <<
1130                     "PadSec="  << iChannel   <<
1131                     "Q="       << cogQ      <<
1132                     "T0="      << cogTime0  <<
1133                     "RMS="     << cogRMS    <<
1134                     "\n";
1135             }
1136             //! debug
1137         }
1138     }
1139     delete fDebugStreamer;
1140     fDebugStreamer = 0x0;
1141 }
1142 //_____________________________________________________________________
1143 void AliTPCCalibPulser::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1144 {
1145     //
1146     //  Write class to file
1147     //
1148
1149     TString sDir(dir);
1150     TString option;
1151
1152     if ( append )
1153         option = "update";
1154     else
1155         option = "recreate";
1156
1157     TDirectory *backup = gDirectory;
1158     TFile f(filename,option.Data());
1159     f.cd();
1160     if ( !sDir.IsNull() ){
1161         f.mkdir(sDir.Data());
1162         f.cd(sDir);
1163     }
1164     this->Write();
1165     f.Close();
1166
1167     if ( backup ) backup->cd();
1168 }
1169 //_____________________________________________________________________
1170 //_________________________  Test Functions ___________________________
1171 //_____________________________________________________________________
1172 TObjArray* AliTPCCalibPulser::TestBinning()
1173 {
1174     //
1175     //  Function to test the binning of the reference histograms
1176     //  type: T0, Q or RMS
1177     //  mode: 0 - number of filled bins per channel
1178     //        1 - number of empty bins between filled bins in one ROC
1179     //  returns TObjArray with the test histograms type*2+mode:
1180     //  position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1181
1182
1183     TObjArray *histArray = new TObjArray(6);
1184     const Char_t *type[] = {"T0","Q","RMS"};
1185     Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1186
1187     for (Int_t itype = 0; itype<3; ++itype){
1188         for (Int_t imode=0; imode<2; ++imode){
1189             Int_t icount = itype*2+imode;
1190             histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1191                                       Form("Test Binning of '%s', mode - %d",type[itype],imode),
1192                                       72,0,72),
1193                              icount);
1194         }
1195     }
1196
1197
1198     TH2S *hRef=0x0;
1199     Short_t *array=0x0;
1200     for (Int_t itype = 0; itype<3; ++itype){
1201         for (Int_t iSec=0; iSec<72; ++iSec){
1202             if ( itype == 0 ) hRef = GetHistoT0(iSec);
1203             if ( itype == 1 ) hRef = GetHistoQ(iSec);
1204             if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1205             if ( hRef == 0x0 ) continue;
1206             array = (hRef->GetArray());
1207             UInt_t nChannels = fROC->GetNChannels(iSec);
1208
1209             Int_t nempty=0;
1210             for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1211                 Int_t nfilled=0;
1212                 Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1213                 Int_t c1 = 0;
1214                 Int_t c2 = 0;
1215                 for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1216                     if ( array[offset+iBin]>0 ) {
1217                         nfilled++;
1218                         if ( c1 && c2 ) nempty++;
1219                         else c1 = 1;
1220                     }
1221                     else if ( c1 ) c2 = 1;
1222
1223                 }
1224                 ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1225             }
1226             ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);
1227         }
1228     }
1229     return histArray;
1230 }