Adding simple example to load default debug streamer
[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),
4c6d06dc 216 fIsZeroSuppressed(kFALSE),
467a4b25 217 fLastSector(-1),
467a4b25 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
4c6d06dc 255 fParam->Update();
467a4b25 256}
257//_____________________________________________________________________
258AliTPCCalibPulser::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),
4c6d06dc 271 fIsZeroSuppressed(sig.fIsZeroSuppressed),
467a4b25 272 fLastSector(-1),
467a4b25 273 fROC(AliTPCROC::Instance()),
04049c81 274 fMapping(NULL),
467a4b25 275 fParam(new AliTPCParam),
276 fPedestalTPC(0x0),
277 fPadNoiseTPC(0x0),
d8faeea7 278 fOutliers(0x0),
467a4b25 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),
3cd27a08 302// fEvent(-1),
467a4b25 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
4c6d06dc 342 fParam->Update();
467a4b25 343}
344//_____________________________________________________________________
345AliTPCCalibPulser& 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//_____________________________________________________________________
356AliTPCCalibPulser::~AliTPCCalibPulser()
357{
358 //
359 // destructor
360 //
361
d8faeea7 362 Reset();
363
364 if ( fDebugStreamer) delete fDebugStreamer;
365 delete fROC;
366 delete fParam;
367}
368void AliTPCCalibPulser::Reset()
369{
370 //
371 // Delete all information: Arrays, Histograms, CalRoc objects
372 //
467a4b25 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();
467a4b25 386}
387//_____________________________________________________________________
3cd27a08 388Int_t AliTPCCalibPulser::Update(const Int_t icsector,
467a4b25 389 const Int_t icRow,
390 const Int_t icPad,
391 const Int_t icTimeBin,
392 const Float_t csignal)
393{
394 //
3cd27a08 395 // Signal filling methode on the fly pedestal and time offset correction if necessary.
467a4b25 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 //
b401648b 399
400 if (icRow<0) return 0;
401 if (icPad<0) return 0;
402 if (icTimeBin<0) return 0;
467a4b25 403 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
404
d8faeea7 405 if ( icRow<0 || icPad<0 ){
406 AliWarning("Wrong Pad or Row number, skipping!");
407 return 0;
408 }
409
467a4b25 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//_____________________________________________________________________
436void 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 //
d8faeea7 493 // what if by chance histo[median] == 0 ?!?
467a4b25 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 }
4c6d06dc 517 fPadPedestal*=(Float_t)(!fIsZeroSuppressed);
467a4b25 518}
519//_____________________________________________________________________
520void 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;
d8faeea7 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
467a4b25 532 const Int_t kCemin = 2; // range for the analysis of the ce signal +- channels from the peak
533 const Int_t kCemax = 7;
d8faeea7 534 param[0] = ceQmax;
535 param[1] = ceTime;
536 param[2] = ceRMS;
537 qSum = ceQsum;
467a4b25 538
d8faeea7 539 if (cemaxpos>0){
540 ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
541 if ( ceQmax<ceMaxThreshold ) return;
467a4b25 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 }
d8faeea7 551 if (ceQsum>ceSumThreshold) {
467a4b25 552 ceTime/=ceQsum;
553 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
d8faeea7 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 }
467a4b25 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//_____________________________________________________________________
3cd27a08 581void AliTPCCalibPulser::ProcessPad()
467a4b25 582{
583 //
584 // Process data of current pad
585 //
586
587 FindPedestal();
588 TVectorD param(3);
3cd27a08 589 Float_t qSum;
590 FindPulserSignal(param, qSum);
467a4b25 591
592 Double_t meanT = param[1];
593 Double_t sigmaT = param[2];
594
d8faeea7 595
467a4b25 596 //Fill Event T0 counter
597 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
598
599 //Fill Q histogram
3cd27a08 600 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
467a4b25 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;
3cd27a08 610 (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
467a4b25 611 }
467a4b25 612 ResetPad();
613}
614//_____________________________________________________________________
3cd27a08 615void AliTPCCalibPulser::EndEvent()
467a4b25 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);
f0312c5d 627 if ( !vTimes || fVTime0OffsetCounter[iSec]==0 ) continue;
628 Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
467a4b25 629 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
3cd27a08 630 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
467a4b25 631
3cd27a08 632 GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
d8faeea7 633 //GetHistoT0(iSec,kTRUE)->Fill( time,iChannel );
467a4b25 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
3cd27a08 649 Float_t q = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
650 Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
467a4b25 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" <<
3cd27a08 669// "Event=" << fEvent <<
467a4b25 670 "Sector="<< sector <<
671 "Row=" << row<<
672 "Pad=" << pad <<
673 "PadC=" << padc <<
674 "PadSec="<< channel <<
3cd27a08 675 "Time0=" << time0 <<
676 "Time=" << time <<
677 "RMS=" << rms <<
678 "Sum=" << q <<
467a4b25 679 "hist.=" << h1 <<
680 "\n";
681
682 delete h1;
683 }
684 //Debug end
467a4b25 685 }
686 }
08205ed7 687}
688//_____________________________________________________________________
689Bool_t AliTPCCalibPulser::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
690{
691 //
692 // Event Processing loop - AliTPCRawStream
693 //
694 ResetEvent();
08205ed7 695 Bool_t withInput = kFALSE;
08205ed7 696 while ( rawStreamFast->NextDDL() ){
4c6d06dc 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 }
08205ed7 710 }
4c6d06dc 711 }
08205ed7 712 }
713 if (withInput){
4c6d06dc 714 EndEvent();
08205ed7 715 }
716 return withInput;
717}
718//_____________________________________________________________________
719Bool_t AliTPCCalibPulser::ProcessEventFast(AliRawReader *rawReader)
720{
721 //
722 // Event processing loop - AliRawReader
723 //
04049c81 724 AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
08205ed7 725 Bool_t res=ProcessEventFast(rawStreamFast);
726 delete rawStreamFast;
727 return res;
467a4b25 728}
729//_____________________________________________________________________
3cd27a08 730Bool_t AliTPCCalibPulser::ProcessEvent(AliTPCRawStream *rawStream)
467a4b25 731{
732 //
733 // Event Processing loop - AliTPCRawStream
734 //
735
467a4b25 736 ResetEvent();
737
738 Bool_t withInput = kFALSE;
739
740 while (rawStream->Next()) {
467a4b25 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 }
467a4b25 750 if (withInput){
751 EndEvent();
752 }
467a4b25 753 return withInput;
754}
755//_____________________________________________________________________
756Bool_t AliTPCCalibPulser::ProcessEvent(AliRawReader *rawReader)
757{
758 //
759 // Event processing loop - AliRawReader
760 //
761
762
04049c81 763 AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
467a4b25 764
765 rawReader->Select("TPC");
766
767 return ProcessEvent(&rawStream);
768}
769//_____________________________________________________________________
770Bool_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//_____________________________________________________________________
3cd27a08 782TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
467a4b25 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//_____________________________________________________________________
3cd27a08 808TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
467a4b25 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//_____________________________________________________________________
3cd27a08 818TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
467a4b25 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//_____________________________________________________________________
3cd27a08 828TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
467a4b25 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//_____________________________________________________________________
3cd27a08 838TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
467a4b25 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//_____________________________________________________________________
3cd27a08 852TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
467a4b25 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//_____________________________________________________________________
3cd27a08 862TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
467a4b25 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//_____________________________________________________________________
3cd27a08 874TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
467a4b25 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//_____________________________________________________________________
3cd27a08 885TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
467a4b25 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//_____________________________________________________________________
3cd27a08 896AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
467a4b25 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//_____________________________________________________________________
3cd27a08 913AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
467a4b25 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//_____________________________________________________________________
3cd27a08 923AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
467a4b25 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//_____________________________________________________________________
3cd27a08 933AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
467a4b25 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//_____________________________________________________________________
943AliTPCCalROC* 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//_____________________________________________________________________
3cd27a08 953void AliTPCCalibPulser::ResetEvent()
467a4b25 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//_____________________________________________________________________
3cd27a08 976void AliTPCCalibPulser::ResetPad()
467a4b25 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//_____________________________________________________________________
989void 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//_____________________________________________________________________
3cd27a08 1039void AliTPCCalibPulser::Analyse()
467a4b25 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
3cd27a08 1062 Short_t *arrayhQ = hQ->GetArray();
1063 Short_t *arrayhT0 = hT0->GetArray();
1064 Short_t *arrayhRMS = hRMS->GetArray();
467a4b25 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
467a4b25 1076 Float_t cogTime0 = -1000;
1077 Float_t cogQ = -1000;
1078 Float_t cogRMS = -1000;
1079 Float_t cogOut = 0;
1080
467a4b25 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;
467a4b25 1084/*
3cd27a08 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);
467a4b25 1088 cogQ = paramQ[1];
1089 cogTime0 = paramT0[1];
1090 cogRMS = paramRMS[1];
1091*/
3cd27a08 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);
467a4b25 1095
467a4b25 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
467a4b25 1135 }
467a4b25 1136 }
1137 delete fDebugStreamer;
1138 fDebugStreamer = 0x0;
1139}
1140//_____________________________________________________________________
1141void 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//_____________________________________________________________________
1170TObjArray* 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}