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