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