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