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