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