Just something left from v2...now clean
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibPulser.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 /////////////////////////////////////////////////////////////////////////////////////////
19 //                                                                                     //
20 //                  Implementation of the TPC pulser calibration                       //
21 //                                                                                     //
22 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch   //
23 //                                                                                     //
24 /////////////////////////////////////////////////////////////////////////////////////////
25 /***************************************************************************
26  *                      Class Description                                  *
27  ***************************************************************************
28
29  The AliTPCCalibPulser class is used to get calibration data concerning the FEE using
30  runs performed with the calibration pulser.
31
32  The information retrieved is
33  - Time0 differences
34  - Signal width differences
35  - Amplification variations
36
37  the seen differences arise from the manufacturing tolerances of the PASAs and are very small within
38  one chip and somewhat large between different chips.
39
40  Histograms:
41    For each ROC three TH2S histos 'Reference Histograms'  (ROC channel vs. [Time0, signal width, Q sum]) is created when
42    it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are stored in the
43    TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.
44
45
46  Working principle:
47  ------------------
48  Raw calibration pulser data is processed by calling one of the ProcessEvent(...) functions
49  (see below). These in the end call the Update(...) function.
50
51  - the Update(...) function:
52    In this function the array fPadSignal is filled with the adc signals between the specified range
53    fFirstTimeBin and fLastTimeBin for the current pad.
54    before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
55    stored in fPadSignal.
56
57    - the ProcessPad() function:
58      Find Pedestal and Noise information
59      - use database information which has to be set by calling
60        SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)
61      - if no information from the pedestal data base
62        is available the informaion is calculated on the fly ( see FindPedestal() function )
63
64      Find the Pulser signal information
65      - calculate  mean = T0, RMS = signal width and Q sum in a range of -2+7 timebins around Q max
66        the Q sum is scaled by pad area
67        (see FindPulserSignal(...) function)
68
69      Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)
70      Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE)),
71
72  At the end of each event the EndEvent() function is called
73
74  - the EndEvent() function:
75    calculate the mean T0 for each ROC and fill the Time0 histogram with Time0-<Time0 for ROC>
76    This is done to overcome syncronisation problems between the trigger and the fec clock.
77
78  After accumulating the desired statistics the Analyse() function has to be called.
79  - the Analyse() function
80    Whithin this function the mean values of T0, RMS, Q are calculated for each pad, using
81    the AliMathBase::GetCOG(...) function, and the calibration
82    storage classes (AliTPCCalROC) are filled for each ROC.
83    The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and
84    fCalRocArrayQ;
85
86
87
88  User interface for filling data:
89  --------------------------------
90
91  To Fill information one of the following functions can be used:
92
93  Bool_t ProcessEvent(eventHeaderStruct *event);
94    - process Date event
95    - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
96
97  Bool_t ProcessEvent(AliRawReader *rawReader);
98    - process AliRawReader event
99    - use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)
100
101  Bool_t ProcessEvent(AliTPCRawStream *rawStream);
102    - process event from AliTPCRawStream
103    - call Update function for signal filling
104
105  Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
106               iPad,  const Int_t iTimeBin, const Float_t signal);
107    - directly  fill signal information (sector, row, pad, time bin, pad)
108      to the reference histograms
109
110  It is also possible to merge two independently taken calibrations using the function
111
112  void Merge(AliTPCCalibPulser *sig)
113    - copy histograms in 'sig' if the do not exist in this instance
114    - Add histograms in 'sig' to the histograms in this instance if the allready exist
115    - After merging call Analyse again!
116
117
118
119  -- example: filling data using root raw data:
120  void fillSignal(Char_t *filename)
121  {
122     rawReader = new AliRawReaderRoot(fileName);
123     if ( !rawReader ) return;
124     AliTPCCalibPulser *calib = new AliTPCCalibPulser;
125     while (rawReader->NextEvent()){
126       calib->ProcessEvent(rawReader);
127     }
128     calib->Analyse();
129     calib->DumpToFile("SignalData.root");
130     delete rawReader;
131     delete calib;
132  }
133
134
135  What kind of information is stored and how to retrieve them:
136  ------------------------------------------------------------
137
138  - Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):
139
140    TH2F *GetHistoT0(Int_t sector);
141    TH2F *GetHistoRMS(Int_t sector);
142    TH2F *GetHistoQ(Int_t sector);
143
144  - Accessing the calibration storage objects:
145
146    AliTPCCalROC *GetCalRocT0(Int_t sector);   // for the Time0 values
147    AliTPCCalROC *GetCalRocRMS(Int_t sector);  // for the signal width values
148    AliTPCCalROC *GetCalRocQ(Int_t sector);    // for the Q sum values
149
150    example for visualisation:
151    if the file "SignalData.root" was created using the above example one could do the following:
152
153    TFile fileSignal("SignalData.root")
154    AliTPCCalibPulser *sig = (AliTPCCalibPulser*)fileSignal->Get("AliTPCCalibPulser");
155    sig->GetCalRocT0(0)->Draw("colz");
156    sig->GetCalRocRMS(0)->Draw("colz");
157
158    or use the AliTPCCalPad functionality:
159    AliTPCCalPad padT0(ped->GetCalPadT0());
160    AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
161    padT0->MakeHisto2D()->Draw("colz");       //Draw A-Side Time0 Information
162    padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
163 */
164
165
166 //Root includes
167 #include <TObjArray.h>
168 #include <TH1F.h>
169 #include <TH2S.h>
170 #include <TString.h>
171 #include <TVectorF.h>
172 #include <TMath.h>
173
174 #include <TDirectory.h>
175 #include <TSystem.h>
176 #include <TROOT.h>
177 #include <TFile.h>
178
179 //AliRoot includes
180 #include "AliRawReader.h"
181 #include "AliRawReaderRoot.h"
182 #include "AliRawReaderDate.h"
183 #include "AliTPCRawStream.h"
184 #include "AliTPCRawStreamFast.h"
185 #include "AliTPCCalROC.h"
186 #include "AliTPCCalPad.h"
187 #include "AliTPCROC.h"
188 #include "AliTPCParam.h"
189 #include "AliTPCCalibPulser.h"
190 #include "AliTPCcalibDB.h"
191 #include "AliMathBase.h"
192 #include "AliLog.h"
193 #include "TTreeStream.h"
194
195 //date
196 #include "event.h"
197
198
199
200
201 ClassImp(AliTPCCalibPulser)
202
203 AliTPCCalibPulser::AliTPCCalibPulser() :
204     TObject(),
205     fFirstTimeBin(60),
206     fLastTimeBin(120),
207     fNbinsT0(200),
208     fXminT0(-2),
209     fXmaxT0(2),
210     fNbinsQ(200),
211     fXminQ(1),
212     fXmaxQ(40),
213     fNbinsRMS(100),
214     fXminRMS(0.1),
215     fXmaxRMS(5.1),
216     fIsZeroSuppressed(kFALSE),
217     fLastSector(-1),
218     fROC(AliTPCROC::Instance()),
219     fMapping(NULL),
220     fParam(new AliTPCParam),
221     fPedestalTPC(0x0),
222     fPadNoiseTPC(0x0),
223     fOutliers(0x0),
224     fPedestalROC(0x0),
225     fPadNoiseROC(0x0),
226     fCalRocArrayT0(72),
227     fCalRocArrayQ(72),
228     fCalRocArrayRMS(72),
229     fCalRocArrayOutliers(72),
230     fHistoQArray(72),
231     fHistoT0Array(72),
232     fHistoRMSArray(72),
233     fPadTimesArrayEvent(72),
234     fPadQArrayEvent(72),
235     fPadRMSArrayEvent(72),
236     fPadPedestalArrayEvent(72),
237     fCurrentChannel(-1),
238     fCurrentSector(-1),
239     fCurrentRow(-1),
240     fMaxPadSignal(-1),
241     fMaxTimeBin(-1),
242     fPadSignal(1024),
243     fPadPedestal(0),
244     fPadNoise(0),
245     fVTime0Offset(72),
246     fVTime0OffsetCounter(72),
247 //    fEvent(-1),
248     fDebugStreamer(0x0),
249     fDebugLevel(0)
250 {
251     //
252     // AliTPCSignal default constructor
253     //
254
255     fParam->Update();
256 }
257 //_____________________________________________________________________
258 AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
259     TObject(sig),
260     fFirstTimeBin(sig.fFirstTimeBin),
261     fLastTimeBin(sig.fLastTimeBin),
262     fNbinsT0(sig.fNbinsT0),
263     fXminT0(sig.fXminT0),
264     fXmaxT0(sig.fXmaxT0),
265     fNbinsQ(sig.fNbinsQ),
266     fXminQ(sig.fXminQ),
267     fXmaxQ(sig.fXmaxQ),
268     fNbinsRMS(sig.fNbinsRMS),
269     fXminRMS(sig.fXminRMS),
270     fXmaxRMS(sig.fXmaxRMS),
271     fIsZeroSuppressed(sig.fIsZeroSuppressed),
272     fLastSector(-1),
273     fROC(AliTPCROC::Instance()),
274     fMapping(NULL),
275     fParam(new AliTPCParam),
276     fPedestalTPC(0x0),
277     fPadNoiseTPC(0x0),
278     fOutliers(0x0),
279     fPedestalROC(0x0),
280     fPadNoiseROC(0x0),
281     fCalRocArrayT0(72),
282     fCalRocArrayQ(72),
283     fCalRocArrayRMS(72),
284     fCalRocArrayOutliers(72),
285     fHistoQArray(72),
286     fHistoT0Array(72),
287     fHistoRMSArray(72),
288     fPadTimesArrayEvent(72),
289     fPadQArrayEvent(72),
290     fPadRMSArrayEvent(72),
291     fPadPedestalArrayEvent(72),
292     fCurrentChannel(-1),
293     fCurrentSector(-1),
294     fCurrentRow(-1),
295     fMaxPadSignal(-1),
296     fMaxTimeBin(-1),
297     fPadSignal(1024),
298     fPadPedestal(0),
299     fPadNoise(0),
300     fVTime0Offset(72),
301     fVTime0OffsetCounter(72),
302 //    fEvent(-1),
303     fDebugStreamer(0x0),
304     fDebugLevel(sig.fDebugLevel)
305 {
306     //
307     // AliTPCSignal default constructor
308     //
309
310     for (Int_t iSec = 0; iSec < 72; ++iSec){
311         const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
312         const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
313         const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
314         const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
315
316         const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
317         const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
318         const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
319
320         if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
321         if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
322         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
323         if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
324
325         if ( hQ != 0x0 ){
326             TH2S *hNew = new TH2S(*hQ);
327             hNew->SetDirectory(0);
328             fHistoQArray.AddAt(hNew,iSec);
329         }
330         if ( hT0 != 0x0 ){
331             TH2S *hNew = new TH2S(*hT0);
332             hNew->SetDirectory(0);
333             fHistoQArray.AddAt(hNew,iSec);
334         }
335         if ( hRMS != 0x0 ){
336             TH2S *hNew = new TH2S(*hRMS);
337             hNew->SetDirectory(0);
338             fHistoQArray.AddAt(hNew,iSec);
339         }
340     }
341
342     fParam->Update();
343 }
344 //_____________________________________________________________________
345 AliTPCCalibPulser& AliTPCCalibPulser::operator = (const  AliTPCCalibPulser &source)
346 {
347   //
348   // assignment operator
349   //
350   if (&source == this) return *this;
351   new (this) AliTPCCalibPulser(source);
352
353   return *this;
354 }
355 //_____________________________________________________________________
356 AliTPCCalibPulser::~AliTPCCalibPulser()
357 {
358     //
359     // destructor
360     //
361
362     Reset();
363
364     if ( fDebugStreamer) delete fDebugStreamer;
365     delete fROC;
366     delete fParam;
367 }
368 void AliTPCCalibPulser::Reset()
369 {
370     //
371     // Delete all information: Arrays, Histograms, CalRoc objects
372     //
373     fCalRocArrayT0.Delete();
374     fCalRocArrayQ.Delete();
375     fCalRocArrayRMS.Delete();
376     fCalRocArrayOutliers.Delete();
377
378     fHistoQArray.Delete();
379     fHistoT0Array.Delete();
380     fHistoRMSArray.Delete();
381
382     fPadTimesArrayEvent.Delete();
383     fPadQArrayEvent.Delete();
384     fPadRMSArrayEvent.Delete();
385     fPadPedestalArrayEvent.Delete();
386 }
387 //_____________________________________________________________________
388 Int_t AliTPCCalibPulser::Update(const Int_t icsector,
389                                 const Int_t icRow,
390                                 const Int_t icPad,
391                                 const Int_t icTimeBin,
392                                 const Float_t csignal)
393 {
394     //
395     // Signal filling methode on the fly pedestal and time offset correction if necessary.
396     // no extra analysis necessary. Assumes knowledge of the signal shape!
397     // assumes that it is looped over consecutive time bins of one pad
398     //
399
400    if (icRow<0) return 0;
401    if (icPad<0) return 0;
402    if (icTimeBin<0) return 0;
403     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
404
405     if ( icRow<0 || icPad<0 ){
406         AliWarning("Wrong Pad or Row number, skipping!");
407         return 0;
408     }
409
410     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
411
412     //init first pad and sector in this event
413     if ( fCurrentChannel == -1 ) {
414         fCurrentChannel = iChannel;
415         fCurrentSector  = icsector;
416         fCurrentRow     = icRow;
417     }
418
419     //process last pad if we change to a new one
420     if ( iChannel != fCurrentChannel ){
421         ProcessPad();
422         fCurrentChannel = iChannel;
423         fCurrentSector  = icsector;
424         fCurrentRow     = icRow;
425     }
426
427     //fill signals for current pad
428     fPadSignal[icTimeBin]=csignal;
429     if ( csignal > fMaxPadSignal ){
430         fMaxPadSignal = csignal;
431         fMaxTimeBin   = icTimeBin;
432     }
433     return 0;
434 }
435 //_____________________________________________________________________
436 void AliTPCCalibPulser::FindPedestal(Float_t part)
437 {
438     //
439     // find pedestal and noise for the current pad. Use either database or
440     // truncated mean with part*100%
441     //
442     Bool_t noPedestal = kTRUE;;
443     if (fPedestalTPC&&fPadNoiseTPC){
444         //use pedestal database
445         //only load new pedestals if the sector has changed
446         if ( fCurrentSector!=fLastSector ){
447             fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
448             fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
449             fLastSector=fCurrentSector;
450         }
451
452         if ( fPedestalROC&&fPadNoiseROC ){
453             fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
454             fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
455             noPedestal   = kFALSE;
456         }
457
458     }
459
460     //if we are not running with pedestal database, or for the current sector there is no information
461     //available, calculate the pedestal and noise on the fly
462     if ( noPedestal ) {
463         const Int_t kPedMax = 100;  //maximum pedestal value
464         Float_t  max    =  0;
465         Float_t  maxPos =  0;
466         Int_t    median =  -1;
467         Int_t    count0 =  0;
468         Int_t    count1 =  0;
469         //
470         Float_t padSignal=0;
471         //
472         UShort_t histo[kPedMax];
473         memset(histo,0,kPedMax*sizeof(UShort_t));
474
475         for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
476             padSignal = fPadSignal.GetMatrixArray()[i];
477             if (padSignal<=0) continue;
478             if (padSignal>max && i>10) {
479                 max = padSignal;
480                 maxPos = i;
481             }
482             if (padSignal>kPedMax-1) continue;
483             histo[Int_t(padSignal+0.5)]++;
484             count0++;
485         }
486             //
487         for (Int_t i=1; i<kPedMax; ++i){
488             if (count1<count0*0.5) median=i;
489             count1+=histo[i];
490         }
491         // truncated mean
492         //
493         // what if by chance histo[median] == 0 ?!?
494         Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
495         //
496         for (Int_t idelta=1; idelta<10; ++idelta){
497             if (median-idelta<=0) continue;
498             if (median+idelta>kPedMax) continue;
499             if (count<part*count1){
500                 count+=histo[median-idelta];
501                 mean +=histo[median-idelta]*(median-idelta);
502                 rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
503                 count+=histo[median+idelta];
504                 mean +=histo[median+idelta]*(median+idelta);
505                 rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
506             }
507         }
508         fPadPedestal = 0;
509         fPadNoise    = 0;
510         if ( count > 0 ) {
511             mean/=count;
512             rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
513             fPadPedestal = mean;
514             fPadNoise    = rms;
515         } 
516     }
517     fPadPedestal*=(Float_t)(!fIsZeroSuppressed);
518 }
519 //_____________________________________________________________________
520 void AliTPCCalibPulser::FindPulserSignal(TVectorD &param, Float_t &qSum)
521 {
522     //
523     //  Find position, signal width and height of the CE signal (last signal)
524     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
525     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
526     //
527
528     Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
529     Int_t   cemaxpos       = fMaxTimeBin;
530     Float_t ceSumThreshold = 10.*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal sum
531     Float_t ceMaxThreshold = 5.*TMath::Max(fPadNoise,Float_t(1.));  // threshold for the signal max
532     const Int_t    kCemin  = 2;             // range for the analysis of the ce signal +- channels from the peak
533     const Int_t    kCemax  = 7;
534     param[0] = ceQmax;
535     param[1] = ceTime;
536     param[2] = ceRMS;
537     qSum     = ceQsum;
538
539     if (cemaxpos>0){
540         ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
541         if ( ceQmax<ceMaxThreshold ) return;
542         for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
543             Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
544             if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
545                 ceTime+=signal*(i+0.5);
546                 ceRMS +=signal*(i+0.5)*(i+0.5);
547                 ceQsum+=signal;
548             }
549         }
550     }
551     if (ceQsum>ceSumThreshold) {
552         ceTime/=ceQsum;
553         ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
554         //only fill the Time0Offset if pad was not marked as an outlier!
555         if ( !fOutliers ){
556             fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
557             fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
558         } else {
559             if ( !(fOutliers->GetCalROC(fCurrentSector)->GetValue(fCurrentChannel)) ){
560                 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
561                 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
562             }
563         }
564
565         //Normalise Q to the pad area
566         Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
567
568         ceQsum/=norm;
569     } else {
570         ceQmax=0;
571         ceTime=0;
572         ceRMS =0;
573         ceQsum=0;
574     }
575     param[0] = ceQmax;
576     param[1] = ceTime;
577     param[2] = ceRMS;
578     qSum     = ceQsum;
579 }
580 //_____________________________________________________________________
581 void AliTPCCalibPulser::ProcessPad()
582 {
583     //
584     //  Process data of current pad
585     //
586
587     FindPedestal();
588     TVectorD param(3);
589     Float_t  qSum;
590     FindPulserSignal(param, qSum);
591
592     Double_t meanT  = param[1];
593     Double_t sigmaT = param[2];
594
595
596     //Fill Event T0 counter
597     (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
598
599     //Fill Q histogram
600     GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
601
602     //Fill RMS histogram
603     GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
604
605
606     //Fill debugging info
607     if ( fDebugLevel>0 ){
608         (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
609         (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
610         (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
611     }
612     ResetPad();
613 }
614 //_____________________________________________________________________
615 void AliTPCCalibPulser::EndEvent()
616 {
617     //
618     //  Process data of current event
619     //
620
621     //check if last pad has allready been processed, if not do so
622     if ( fMaxTimeBin>-1 ) ProcessPad();
623
624     //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
625     for ( Int_t iSec = 0; iSec<72; ++iSec ){
626         TVectorF *vTimes = GetPadTimesEvent(iSec);
627         if ( !vTimes || fVTime0OffsetCounter[iSec]==0 ) continue;
628         Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
629         for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
630             Float_t time  = (*vTimes).GetMatrixArray()[iChannel];
631
632             GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
633             //GetHistoT0(iSec,kTRUE)->Fill( time,iChannel );
634
635
636             //Debug start
637             if ( fDebugLevel>0 ){
638                 if ( !fDebugStreamer ) {
639                         //debug stream
640                     TDirectory *backup = gDirectory;
641                     fDebugStreamer = new TTreeSRedirector("deb2.root");
642                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
643                 }
644
645                 Int_t row=0;
646                 Int_t pad=0;
647                 Int_t padc=0;
648
649                 Float_t q   = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
650                 Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
651
652                 UInt_t channel=iChannel;
653                 Int_t sector=iSec;
654
655                 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
656                 pad = channel-fROC->GetRowIndexes(sector)[row];
657                 padc = pad-(fROC->GetNPads(sector,row)/2);
658
659                 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
660                                     Form("hSignalD%d.%d.%d",sector,row,pad),
661                                     fLastTimeBin-fFirstTimeBin,
662                                     fFirstTimeBin,fLastTimeBin);
663                 h1->SetDirectory(0);
664
665                 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
666                     h1->Fill(i,fPadSignal(i));
667
668                 (*fDebugStreamer) << "DataPad" <<
669 //                  "Event=" << fEvent <<
670                     "Sector="<< sector <<
671                     "Row="   << row<<
672                     "Pad="   << pad <<
673                     "PadC="  << padc <<
674                     "PadSec="<< channel <<
675                     "Time0="  << time0 <<
676                     "Time="  << time <<
677                     "RMS="   << rms <<
678                     "Sum="   << q <<
679                     "hist.=" << h1 <<
680                     "\n";
681
682                 delete h1;
683             }
684             //Debug end
685         }
686     }
687 }
688 //_____________________________________________________________________
689 Bool_t AliTPCCalibPulser::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
690 {
691   //
692   // Event Processing loop - AliTPCRawStream
693   //
694   ResetEvent();
695   Bool_t withInput = kFALSE;
696   while ( rawStreamFast->NextDDL() ){
697     while ( rawStreamFast->NextChannel() ){
698       Int_t isector  = rawStreamFast->GetSector();                       //  current sector
699       Int_t iRow     = rawStreamFast->GetRow();                          //  current row
700       Int_t iPad     = rawStreamFast->GetPad();                          //  current pad
701
702       while ( rawStreamFast->NextBunch() ){
703         Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
704         Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
705         for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
706           Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
707           Update(isector,iRow,iPad,iTimeBin+1,signal);
708           withInput = kTRUE;
709         }
710       }
711     }
712   }
713   if (withInput){
714     EndEvent();
715   }
716   return withInput;
717 }
718 //_____________________________________________________________________
719 Bool_t AliTPCCalibPulser::ProcessEventFast(AliRawReader *rawReader)
720 {
721   //
722   //  Event processing loop - AliRawReader
723   //
724   AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
725   Bool_t res=ProcessEventFast(rawStreamFast);
726   delete rawStreamFast;
727   return res;
728 }
729 //_____________________________________________________________________
730 Bool_t AliTPCCalibPulser::ProcessEvent(AliTPCRawStream *rawStream)
731 {
732   //
733   // Event Processing loop - AliTPCRawStream
734   //
735
736   ResetEvent();
737
738   Bool_t withInput = kFALSE;
739
740   while (rawStream->Next()) {
741       Int_t isector  = rawStream->GetSector();                       //  current sector
742       Int_t iRow     = rawStream->GetRow();                          //  current row
743       Int_t iPad     = rawStream->GetPad();                          //  current pad
744       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
745       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
746
747       Update(isector,iRow,iPad,iTimeBin,signal);
748       withInput = kTRUE;
749   }
750   if (withInput){
751       EndEvent();
752   }
753   return withInput;
754 }
755 //_____________________________________________________________________
756 Bool_t AliTPCCalibPulser::ProcessEvent(AliRawReader *rawReader)
757 {
758   //
759   //  Event processing loop - AliRawReader
760   //
761
762
763   AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
764
765   rawReader->Select("TPC");
766
767   return ProcessEvent(&rawStream);
768 }
769 //_____________________________________________________________________
770 Bool_t AliTPCCalibPulser::ProcessEvent(eventHeaderStruct *event)
771 {
772   //
773   //  Event processing loop - date event
774   //
775     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
776     Bool_t result=ProcessEvent(rawReader);
777     delete rawReader;
778     return result;
779
780 }
781 //_____________________________________________________________________
782 TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
783                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
784                                   Char_t *type, Bool_t force)
785 {
786     //
787     // return pointer to Q histogram
788     // if force is true create a new histogram if it doesn't exist allready
789     //
790     if ( !force || arr->UncheckedAt(sector) )
791         return (TH2S*)arr->UncheckedAt(sector);
792
793     // if we are forced and histogram doesn't yes exist create it
794     Char_t name[255], title[255];
795
796     sprintf(name,"hCalib%s%.2d",type,sector);
797     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
798
799     // new histogram with Q calib information. One value for each pad!
800     TH2S* hist = new TH2S(name,title,
801                           nbinsY, ymin, ymax,
802                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
803     hist->SetDirectory(0);
804     arr->AddAt(hist,sector);
805     return hist;
806 }
807 //_____________________________________________________________________
808 TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
809 {
810     //
811     // return pointer to T0 histogram
812     // if force is true create a new histogram if it doesn't exist allready
813     //
814     TObjArray *arr = &fHistoT0Array;
815     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
816 }
817 //_____________________________________________________________________
818 TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
819 {
820     //
821     // return pointer to Q histogram
822     // if force is true create a new histogram if it doesn't exist allready
823     //
824     TObjArray *arr = &fHistoQArray;
825     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
826 }
827 //_____________________________________________________________________
828 TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
829 {
830     //
831     // return pointer to Q histogram
832     // if force is true create a new histogram if it doesn't exist allready
833     //
834     TObjArray *arr = &fHistoRMSArray;
835     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
836 }
837 //_____________________________________________________________________
838 TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
839 {
840     //
841     // return pointer to Pad Info from 'arr' for the current event and sector
842     // if force is true create it if it doesn't exist allready
843     //
844     if ( !force || arr->UncheckedAt(sector) )
845         return (TVectorF*)arr->UncheckedAt(sector);
846
847     TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
848     arr->AddAt(vect,sector);
849     return vect;
850 }
851 //_____________________________________________________________________
852 TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
853 {
854     //
855     // return pointer to Pad Times Array for the current event and sector
856     // if force is true create it if it doesn't exist allready
857     //
858     TObjArray *arr = &fPadTimesArrayEvent;
859     return GetPadInfoEvent(sector,arr,force);
860 }
861 //_____________________________________________________________________
862 TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
863 {
864     //
865     // return pointer to Pad Q Array for the current event and sector
866     // if force is true create it if it doesn't exist allready
867     // for debugging purposes only
868     //
869
870     TObjArray *arr = &fPadQArrayEvent;
871     return GetPadInfoEvent(sector,arr,force);
872 }
873 //_____________________________________________________________________
874 TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
875 {
876     //
877     // return pointer to Pad RMS Array for the current event and sector
878     // if force is true create it if it doesn't exist allready
879     // for debugging purposes only
880     //
881     TObjArray *arr = &fPadRMSArrayEvent;
882     return GetPadInfoEvent(sector,arr,force);
883 }
884 //_____________________________________________________________________
885 TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
886 {
887     //
888     // return pointer to Pad RMS Array for the current event and sector
889     // if force is true create it if it doesn't exist allready
890     // for debugging purposes only
891     //
892     TObjArray *arr = &fPadPedestalArrayEvent;
893     return GetPadInfoEvent(sector,arr,force);
894 }
895 //_____________________________________________________________________
896 AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
897 {
898     //
899     // return pointer to ROC Calibration
900     // if force is true create a new histogram if it doesn't exist allready
901     //
902     if ( !force || arr->UncheckedAt(sector) )
903         return (AliTPCCalROC*)arr->UncheckedAt(sector);
904
905     // if we are forced and histogram doesn't yes exist create it
906
907     // new AliTPCCalROC for T0 information. One value for each pad!
908     AliTPCCalROC *croc = new AliTPCCalROC(sector);
909     arr->AddAt(croc,sector);
910     return croc;
911 }
912 //_____________________________________________________________________
913 AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
914 {
915     //
916     // return pointer to Carge ROC Calibration
917     // if force is true create a new histogram if it doesn't exist allready
918     //
919     TObjArray *arr = &fCalRocArrayT0;
920     return GetCalRoc(sector, arr, force);
921 }
922 //_____________________________________________________________________
923 AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
924 {
925     //
926     // return pointer to T0 ROC Calibration
927     // if force is true create a new histogram if it doesn't exist allready
928     //
929     TObjArray *arr = &fCalRocArrayQ;
930     return GetCalRoc(sector, arr, force);
931 }
932 //_____________________________________________________________________
933 AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
934 {
935     //
936     // return pointer to signal width ROC Calibration
937     // if force is true create a new histogram if it doesn't exist allready
938     //
939     TObjArray *arr = &fCalRocArrayRMS;
940     return GetCalRoc(sector, arr, force);
941 }
942 //_____________________________________________________________________
943 AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
944 {
945     //
946     // return pointer to Outliers
947     // if force is true create a new histogram if it doesn't exist allready
948     //
949     TObjArray *arr = &fCalRocArrayOutliers;
950     return GetCalRoc(sector, arr, force);
951 }
952 //_____________________________________________________________________
953 void AliTPCCalibPulser::ResetEvent()
954 {
955     //
956     //  Reset global counters  -- Should be called before each event is processed
957     //
958     fLastSector=-1;
959     fCurrentSector=-1;
960     fCurrentRow=-1;
961     fCurrentChannel=-1;
962
963     ResetPad();
964
965     fPadTimesArrayEvent.Delete();
966     fPadQArrayEvent.Delete();
967     fPadRMSArrayEvent.Delete();
968     fPadPedestalArrayEvent.Delete();
969
970     for ( Int_t i=0; i<72; ++i ){
971         fVTime0Offset[i]=0;
972         fVTime0OffsetCounter[i]=0;
973     }
974 }
975 //_____________________________________________________________________
976 void AliTPCCalibPulser::ResetPad()
977 {
978     //
979     //  Reset pad infos -- Should be called after a pad has been processed
980     //
981     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
982         fPadSignal[i] = 0;
983     fMaxTimeBin = -1;
984     fMaxPadSignal = -1;
985     fPadPedestal  = -1;
986     fPadNoise     = -1;
987 }
988 //_____________________________________________________________________
989 void AliTPCCalibPulser::Merge(AliTPCCalibPulser *sig)
990 {
991     //
992     //  Merge reference histograms of sig to the current AliTPCCalibPulser
993     //
994
995     //merge histograms
996     for (Int_t iSec=0; iSec<72; ++iSec){
997         TH2S *hRefQmerge   = sig->GetHistoQ(iSec);
998         TH2S *hRefT0merge  = sig->GetHistoT0(iSec);
999         TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
1000
1001
1002         if ( hRefQmerge ){
1003             TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1004             TH2S *hRefQ   = GetHistoQ(iSec);
1005             if ( hRefQ ) hRefQ->Add(hRefQmerge);
1006             else {
1007                 TH2S *hist = new TH2S(*hRefQmerge);
1008                 hist->SetDirectory(0);
1009                 fHistoQArray.AddAt(hist, iSec);
1010             }
1011             hRefQmerge->SetDirectory(dir);
1012         }
1013         if ( hRefT0merge ){
1014             TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1015             TH2S *hRefT0  = GetHistoT0(iSec);
1016             if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1017             else {
1018                 TH2S *hist = new TH2S(*hRefT0merge);
1019                 hist->SetDirectory(0);
1020                 fHistoT0Array.AddAt(hist, iSec);
1021             }
1022             hRefT0merge->SetDirectory(dir);
1023         }
1024         if ( hRefRMSmerge ){
1025             TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1026             TH2S *hRefRMS = GetHistoRMS(iSec);
1027             if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1028             else {
1029                 TH2S *hist = new TH2S(*hRefRMSmerge);
1030                 hist->SetDirectory(0);
1031                 fHistoRMSArray.AddAt(hist, iSec);
1032             }
1033             hRefRMSmerge->SetDirectory(dir);
1034         }
1035
1036     }
1037 }
1038 //_____________________________________________________________________
1039 void AliTPCCalibPulser::Analyse()
1040 {
1041     //
1042     //  Calculate calibration constants
1043     //
1044
1045     TVectorD paramQ(3);
1046     TVectorD paramT0(3);
1047     TVectorD paramRMS(3);
1048     TMatrixD dummy(3,3);
1049
1050     for (Int_t iSec=0; iSec<72; ++iSec){
1051         TH2S *hT0 = GetHistoT0(iSec);
1052         if (!hT0 ) continue;
1053
1054         AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
1055         AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
1056         AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1057         AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1058
1059         TH2S *hQ   = GetHistoQ(iSec);
1060         TH2S *hRMS = GetHistoRMS(iSec);
1061
1062         Short_t *arrayhQ   = hQ->GetArray();
1063         Short_t *arrayhT0  = hT0->GetArray();
1064         Short_t *arrayhRMS = hRMS->GetArray();
1065
1066         UInt_t nChannels = fROC->GetNChannels(iSec);
1067
1068         //debug
1069         Int_t row=0;
1070         Int_t pad=0;
1071         Int_t padc=0;
1072         //! debug
1073
1074         for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1075
1076             Float_t cogTime0 = -1000;
1077             Float_t cogQ     = -1000;
1078             Float_t cogRMS   = -1000;
1079             Float_t cogOut   = 0;
1080
1081             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1082             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1083             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1084 /*
1085             AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1086             AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1087             AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1088             cogQ     = paramQ[1];
1089             cogTime0 = paramT0[1];
1090             cogRMS   = paramRMS[1];
1091 */
1092             cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1093             cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1094             cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1095
1096             /*
1097             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1098                 cogOut = 1;
1099                 cogTime0 = 0;
1100                 cogQ     = 0;
1101                 cogRMS   = 0;
1102             }
1103 */
1104             rocQ->SetValue(iChannel, cogQ*cogQ);
1105             rocT0->SetValue(iChannel, cogTime0);
1106             rocRMS->SetValue(iChannel, cogRMS);
1107             rocOut->SetValue(iChannel, cogOut);
1108
1109
1110             //debug
1111             if ( fDebugLevel > 0 ){
1112                 if ( !fDebugStreamer ) {
1113                         //debug stream
1114                     TDirectory *backup = gDirectory;
1115                     fDebugStreamer = new TTreeSRedirector("deb2.root");
1116                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1117                 }
1118
1119                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1120                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1121                 padc = pad-(fROC->GetNPads(iSec,row)/2);
1122
1123                 (*fDebugStreamer) << "DataEnd" <<
1124                     "Sector="  << iSec      <<
1125                     "Pad="     << pad       <<
1126                     "PadC="    << padc      <<
1127                     "Row="     << row       <<
1128                     "PadSec="  << iChannel   <<
1129                     "Q="       << cogQ      <<
1130                     "T0="      << cogTime0  <<
1131                     "RMS="     << cogRMS    <<
1132                     "\n";
1133             }
1134             //! debug
1135         }
1136     }
1137     delete fDebugStreamer;
1138     fDebugStreamer = 0x0;
1139 }
1140 //_____________________________________________________________________
1141 void AliTPCCalibPulser::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1142 {
1143     //
1144     //  Write class to file
1145     //
1146
1147     TString sDir(dir);
1148     TString option;
1149
1150     if ( append )
1151         option = "update";
1152     else
1153         option = "recreate";
1154
1155     TDirectory *backup = gDirectory;
1156     TFile f(filename,option.Data());
1157     f.cd();
1158     if ( !sDir.IsNull() ){
1159         f.mkdir(sDir.Data());
1160         f.cd(sDir);
1161     }
1162     this->Write();
1163     f.Close();
1164
1165     if ( backup ) backup->cd();
1166 }
1167 //_____________________________________________________________________
1168 //_________________________  Test Functions ___________________________
1169 //_____________________________________________________________________
1170 TObjArray* AliTPCCalibPulser::TestBinning()
1171 {
1172     //
1173     //  Function to test the binning of the reference histograms
1174     //  type: T0, Q or RMS
1175     //  mode: 0 - number of filled bins per channel
1176     //        1 - number of empty bins between filled bins in one ROC
1177     //  returns TObjArray with the test histograms type*2+mode:
1178     //  position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1179
1180
1181     TObjArray *histArray = new TObjArray(6);
1182     const Char_t *type[] = {"T0","Q","RMS"};
1183     Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1184
1185     for (Int_t itype = 0; itype<3; ++itype){
1186         for (Int_t imode=0; imode<2; ++imode){
1187             Int_t icount = itype*2+imode;
1188             histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1189                                       Form("Test Binning of '%s', mode - %d",type[itype],imode),
1190                                       72,0,72),
1191                              icount);
1192         }
1193     }
1194
1195
1196     TH2S *hRef=0x0;
1197     Short_t *array=0x0;
1198     for (Int_t itype = 0; itype<3; ++itype){
1199         for (Int_t iSec=0; iSec<72; ++iSec){
1200             if ( itype == 0 ) hRef = GetHistoT0(iSec);
1201             if ( itype == 1 ) hRef = GetHistoQ(iSec);
1202             if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1203             if ( hRef == 0x0 ) continue;
1204             array = (hRef->GetArray());
1205             UInt_t nChannels = fROC->GetNChannels(iSec);
1206
1207             Int_t nempty=0;
1208             for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1209                 Int_t nfilled=0;
1210                 Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1211                 Int_t c1 = 0;
1212                 Int_t c2 = 0;
1213                 for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1214                     if ( array[offset+iBin]>0 ) {
1215                         nfilled++;
1216                         if ( c1 && c2 ) nempty++;
1217                         else c1 = 1;
1218                     }
1219                     else if ( c1 ) c2 = 1;
1220
1221                 }
1222                 ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1223             }
1224             ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);
1225         }
1226     }
1227     return histArray;
1228 }