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