4842f327eb96f7a862ba77a5e762a11944c60b21
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibCE.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
17
18
19
20
21
22 /* $Id$ */
23
24
25
26 //Root includes
27 #include <TObjArray.h>
28 #include <TH1F.h>
29 #include <TH2S.h>
30 #include <TString.h>
31 #include <TVectorF.h>
32 #include <TVectorD.h>
33 #include <TMatrixD.h>
34 #include <TMath.h>
35 #include <TGraph.h>
36 #include <TString.h>
37
38 #include <TDirectory.h>
39 #include <TSystem.h>
40 #include <TFile.h>
41
42 //AliRoot includes
43 #include "AliRawReader.h"
44 #include "AliRawReaderRoot.h"
45 #include "AliRawReaderDate.h"
46 #include "AliRawEventHeaderBase.h"
47 #include "AliTPCRawStream.h"
48 #include "AliTPCcalibDB.h"
49 #include "AliTPCCalROC.h"
50 #include "AliTPCCalPad.h"
51 #include "AliTPCROC.h"
52 #include "AliTPCParam.h"
53 #include "AliTPCCalibCE.h"
54 #include "AliMathBase.h"
55 #include "TTreeStream.h"
56
57 //date
58 #include "event.h"
59 ClassImp(AliTPCCalibCE)
60
61 //////////////////////////////////////////////////////////////////////////////////////
62 //          Implementation of the TPC Central Electrode calibration
63 //
64 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
65 // 
66 //
67 //
68 // *************************************************************************************
69 // *                                Class Description                                  *
70 // *************************************************************************************
71 //
72 /* BEGIN_HTML
73  <h4>The AliTPCCalibCE class is used to get calibration data from the Central Electrode
74  using laser runs.</h4>
75
76  The information retrieved is
77  <ul style="list-style-type: square;">
78    <li>Time arrival from the CE</li>
79    <li>Signal width</li>
80    <li>Signal sum</li>
81  </ul>
82
83 <h4>Overview:</h4>
84  <ol style="list-style-type: upper-roman;">
85    <li><a href="#working">Working principle</a></li>
86    <li><a href="#user">User interface for filling data</a></li>
87    <li><a href="#info">Stored information</a></li>
88  </ol>
89
90  <h3><a name="working">I. Working principle</a></h3>
91
92  <h4>Raw laser data is processed by calling one of the ProcessEvent(...) functions
93  (see below). These in the end call the Update(...) function.</h4>
94
95  <ul style="list-style-type: square;">
96    <li>the Update(...) function:<br />
97        In this function the array fPadSignal is filled with the adc signals between the specified range
98        fFirstTimeBin and fLastTimeBin for the current pad.
99        before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
100        stored in fPadSignal.
101    </li>
102    <ul style="list-style-type: square;">
103    <li>the ProcessPad() function:</li>
104    <ol style="list-style-type: decimal;">
105      <li>Find Pedestal and Noise information</li>
106      <ul style="list-style-type: square;">
107        <li>use database information which has to be set by calling<br />
108            SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)</li>
109        <li>if no information from the pedestal data base
110            is available the informaion is calculated on the fly
111            ( see FindPedestal() function )</li>
112      </ul>
113      <li>Find local maxima of the pad signal</li>
114      <ul style="list-style-type: square;">
115        <li>maxima arise from the laser tracks, the CE and also periodic postpeaks after the CE signal have
116        have been observed ( see FindLocalMaxima(...) )</li>
117      </ul>
118      <li>Find the CE signal information</li>
119      <ul style="list-style-type: square;">
120        <li>to find the position of the CE signal the Tmean information from the previos event is used
121            as the CE signal the local maximum closest to this Tmean is identified</li>
122        <li>calculate  mean = T0, RMS = signal width and Q sum in a range of -4+7 timebins around Q max position
123            the Q sum is scaled by pad area (see FindPulserSignal(...) function)</li>
124      </ul>
125      <li>Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)</li>
126      <li>Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE))</li>
127      </ol>
128    </ul>
129  </ul>
130
131  <h4>At the end of each event the EndEvent() function is called</h4>
132
133  <ul style="list-style-type: square;">
134    <li>the EndEvent() function:</li>
135    <ul style="list-style-type: square;">
136      <li>calculate the mean T0 for side A and side C. Fill T0 histogram with Time0-<Time0 for side[A,C]>
137          This is done to overcome syncronisation problems between the trigger and the fec clock.</li>
138      <li>calculate Mean T for each ROC using the COG aroud the median of the LocalMaxima distribution in one sector</li>
139      <li>calculate Mean Q</li>
140      <li>calculate Global fit parameters for Pol1 and Pol2 fits</li>
141    </ul>
142  </ul>
143
144  <h4>After accumulating the desired statistics the Analyse() function has to be called.</h4>
145   <ul style="list-style-type: square;">
146   <li>the Analyse() function:</li>
147     <ul style="list-style-type: square;">
148       <li>calculate the mean values of T0, RMS, Q for each pad, using
149           the AliMathBase::GetCOG(...) function</li>
150       <li>fill the calibration storage classes (AliTPCCalROC) for each ROC</li>
151          (The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and fCalRocArrayQ</li>
152     </ul>
153   </ul>
154
155  <h3><a name="user">II. User interface for filling data</a></h3>
156
157  <h4>To Fill information one of the following functions can be used:</h4>
158
159  <ul style="list-style-type: none;">
160   <li> Bool_t ProcessEvent(eventHeaderStruct *event);</li>
161     <ul style="list-style-type: square;">
162       <li>process Date event</li>
163       <li>use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)</li>
164     </ul>
165     <br />
166
167   <li> Bool_t ProcessEvent(AliRawReader *rawReader);</li>
168     <ul style="list-style-type: square;">
169       <li>process AliRawReader event</li>
170       <li>use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)</li>
171     </ul>
172     <br />
173
174   <li> Bool_t ProcessEvent(AliTPCRawStream *rawStream);</li>
175     <ul style="list-style-type: square;">
176       <li>process event from AliTPCRawStream</li>
177       <li>call Update function for signal filling</li>
178     </ul>
179     <br />
180
181   <li> Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
182               iPad,  const Int_t iTimeBin, const Float_t signal);</li>
183     <ul style="list-style-type: square;">
184       <li>directly  fill signal information (sector, row, pad, time bin, pad)
185           to the reference histograms</li>
186     </ul>
187  </ul>
188
189  <h4>It is also possible to merge two independently taken calibrations using the function</h4>
190
191  <ul style="list-style-type: none;">
192   <li> void Merge(AliTPCCalibSignal *sig)</li>
193     <ul style="list-style-type: square;">
194       <li>copy histograms in 'sig' if they do not exist in this instance</li>
195       <li>Add histograms in 'sig' to the histograms in this instance if the allready exist</li>
196       <li>After merging call Analyse again!</li>
197     </ul>
198  </ul>
199
200
201  <h4>example: filling data using root raw data:</h4>
202  <pre> 
203  void fillCE(Char_t *filename)
204  {
205     rawReader = new AliRawReaderRoot(fileName);
206     if ( !rawReader ) return;
207     AliTPCCalibCE *calib = new AliTPCCalibCE;
208     while (rawReader->NextEvent()){
209       calib->ProcessEvent(rawReader);
210     }
211     calib->Analyse();
212     calib->DumpToFile("CEData.root");
213     delete rawReader;
214     delete calib;
215  }
216  </pre>
217
218  <h3><a name="info">III. What kind of information is stored and how to retrieve it</a></h4>
219
220  <h4><a name="info:stored">III.1 Stored information</a></h4>
221  <ul style="list-style-type: none;">
222   <li>Histograms:</li>
223   <ul style="list-style-type: none;">
224     <li>For each ROC three TH2S histos 'Reference Histograms'  (ROC channel vs. [Time0, signal width, Q sum])
225         is created when it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are
226         stored in the TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.</li>
227   </ul>
228   <br />
229
230  <li>Calibration Data:</li>
231  <ul style="list-style-type: none;">
232       <li>For each ROC three types of calibration data (AliTPCCalROC) is stored: for the mean arrival Time,
233           the signal width and the signal Sum. The AliTPCCalROC objects are stored in the TObjArrays
234           fCalRocArrayT0, fCalRocArrayRMS , fCalRocArrayQ. The object for each roc is created the first time it
235           is accessed (GetCalRoc[T0,RMS,Q](ROC,kTRUE));</li>
236  </ul>
237  <br />
238
239  <li>For each event the following information is stored:</li>
240    
241  <ul style="list-style-type: square;">
242    <li>event time ( TVectorD  fVEventTime )</li>
243    <li>event id   ( TVectorD  fVEventNumber )</li>
244    <br />
245    <li>mean arrival time for each ROC                ( TObjArray fTMeanArrayEvent )</li>
246    <li>mean Q for each ROC                           ( TObjArray fQMeanArrayEvent )</li>
247    <li>parameters of a plane fit for each ROC        ( TObjArray fParamArrayEventPol1 )</li>
248    <li>parameters of a 2D parabola fit for each ROC  ( TObjArray fParamArrayEventPol2 )</li>
249   </ul>
250  </ul>
251
252  <h4><a name="info:retrieve">III.2 Retrieving information</a></h4>
253  <ul style="list-style-type: none;">
254   <li>Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):</li>
255     <ul style="list-style-type: square;">
256       <li>TH2F *GetHistoT0(Int_t sector);</li>
257       <li>TH2F *GetHistoRMS(Int_t sector);</li>
258       <li>TH2F *GetHistoQ(Int_t sector);</li>
259     </ul>
260     <br />
261     
262   <li>Accessing the calibration storage objects:</li>
263     <ul style="list-style-type: square;">
264       <li>AliTPCCalROC *GetCalRocT0(Int_t sector);   // for the Time0 values</li>
265       <li>AliTPCCalROC *GetCalRocRMS(Int_t sector);  // for the signal width values</li>
266       <li>AliTPCCalROC *GetCalRocQ(Int_t sector);    // for the Q sum values</li>
267     </ul>
268     <br />
269
270   <li>Accessin the event by event information:</li>
271     <ul style="list-style-type: square;">
272       <li>The event by event information can be displayed using the</li>
273       <li>MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)</li>
274       <li>which creates a graph from the specified variables</li>
275     </ul>
276   </ul>
277   
278   <h4>example for visualisation:</h4>
279   <pre>
280   //if the file "CEData.root" was created using the above example one could do the following:
281   TFile fileCE("CEData.root")
282   AliTPCCalibCE *ce = (AliTPCCalibCE*)fileCE->Get("AliTPCCalibCE");
283   ce->GetCalRocT0(0)->Draw("colz");
284   ce->GetCalRocRMS(0)->Draw("colz");
285
286   //or use the AliTPCCalPad functionality:
287   AliTPCCalPad padT0(ped->GetCalPadT0());
288   AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
289   padT0->MakeHisto2D()->Draw("colz");       //Draw A-Side Time0 Information
290   padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
291
292   //display event by event information:
293   //Draw mean arrival time as a function of the event time for oroc sector A00
294   ce->MakeGraphTimeCE(36, 0, 2)->Draw("alp");
295   //Draw first derivative in local x from a plane fit as a function of the event time for oroc sector A00
296   ce->MakeGraphTimeCE(36, 0, 0, 1)->Draw("alp");  
297   </pre>
298 END_HTML */
299 //////////////////////////////////////////////////////////////////////////////////////
300
301
302 AliTPCCalibCE::AliTPCCalibCE() :
303     TObject(),
304     fFirstTimeBin(650),
305     fLastTimeBin(1000),
306     fNbinsT0(200),
307     fXminT0(-5),
308     fXmaxT0(5),
309     fNbinsQ(200),
310     fXminQ(1),
311     fXmaxQ(40),
312     fNbinsRMS(100),
313     fXminRMS(0.1),
314     fXmaxRMS(5.1),
315     fLastSector(-1),
316     fOldRCUformat(kTRUE),
317     fROC(AliTPCROC::Instance()),
318     fParam(new AliTPCParam),
319     fPedestalTPC(0x0),
320     fPadNoiseTPC(0x0),
321     fPedestalROC(0x0),
322     fPadNoiseROC(0x0),
323     fCalRocArrayT0(72),
324     fCalRocArrayQ(72),
325     fCalRocArrayRMS(72),
326     fCalRocArrayOutliers(72),
327     fHistoQArray(72),
328     fHistoT0Array(72),
329     fHistoRMSArray(72),
330     fHistoTmean(72),
331     fParamArrayEventPol1(72),
332     fParamArrayEventPol2(72),
333     fTMeanArrayEvent(72),
334     fQMeanArrayEvent(72),
335     fVEventTime(10),
336     fVEventNumber(10),
337     fNevents(0),
338     fTimeStamp(0),
339     fEventId(-1),
340     fRunNumber(-1),
341     fOldRunNumber(-1),
342     fPadTimesArrayEvent(72),
343     fPadQArrayEvent(72),
344     fPadRMSArrayEvent(72),
345     fPadPedestalArrayEvent(72),
346     fCurrentChannel(-1),
347     fCurrentSector(-1),
348     fCurrentRow(-1),
349     fMaxPadSignal(-1),
350     fMaxTimeBin(-1),
351     fPadSignal(1024),
352     fPadPedestal(0),
353     fPadNoise(0),
354     fVTime0Offset(72),
355     fVTime0OffsetCounter(72),
356     fVMeanQ(72),
357     fVMeanQCounter(72),
358     fEvent(-1),
359     fDebugStreamer(0x0),
360     fDebugLevel(0)
361 {
362     //
363     // AliTPCSignal default constructor
364     //
365 //    fHTime0 = new TH1F("hTime0Event","hTime0Event",(fLastTimeBin-fFirstTimeBin)*10,fFirstTimeBin,fLastTimeBin);
366 }
367 //_____________________________________________________________________
368 AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
369     TObject(sig),
370     fFirstTimeBin(sig.fFirstTimeBin),
371     fLastTimeBin(sig.fLastTimeBin),
372     fNbinsT0(sig.fNbinsT0),
373     fXminT0(sig.fXminT0),
374     fXmaxT0(sig.fXmaxT0),
375     fNbinsQ(sig.fNbinsQ),
376     fXminQ(sig.fXminQ),
377     fXmaxQ(sig.fXmaxQ),
378     fNbinsRMS(sig.fNbinsRMS),
379     fXminRMS(sig.fXminRMS),
380     fXmaxRMS(sig.fXmaxRMS),
381     fLastSector(-1),
382     fOldRCUformat(kTRUE),
383     fROC(AliTPCROC::Instance()),
384     fParam(new AliTPCParam),
385     fPedestalTPC(0x0),
386     fPadNoiseTPC(0x0),
387     fPedestalROC(0x0),
388     fPadNoiseROC(0x0),
389     fCalRocArrayT0(72),
390     fCalRocArrayQ(72),
391     fCalRocArrayRMS(72),
392     fCalRocArrayOutliers(72),
393     fHistoQArray(72),
394     fHistoT0Array(72),
395     fHistoRMSArray(72),
396     fHistoTmean(72),
397     fParamArrayEventPol1(72),
398     fParamArrayEventPol2(72),
399     fTMeanArrayEvent(72),
400     fQMeanArrayEvent(72),
401     fVEventTime(1000),
402     fVEventNumber(1000),
403     fNevents(sig.fNevents),
404     fTimeStamp(0),
405     fEventId(-1),
406     fRunNumber(-1),
407     fOldRunNumber(-1),
408     fPadTimesArrayEvent(72),
409     fPadQArrayEvent(72),
410     fPadRMSArrayEvent(72),
411     fPadPedestalArrayEvent(72),
412     fCurrentChannel(-1),
413     fCurrentSector(-1),
414     fCurrentRow(-1),
415     fMaxPadSignal(-1),
416     fMaxTimeBin(-1),
417     fPadSignal(1024),
418     fPadPedestal(0),
419     fPadNoise(0),
420     fVTime0Offset(72),
421     fVTime0OffsetCounter(72),
422     fVMeanQ(72),
423     fVMeanQCounter(72),
424     fEvent(-1),
425     fDebugStreamer(0x0),
426     fDebugLevel(sig.fDebugLevel)
427 {
428     //
429     // AliTPCSignal copy constructor
430     //
431
432     for (Int_t iSec = 0; iSec < 72; ++iSec){
433         const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
434         const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
435         const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
436         const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
437
438         const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
439         const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
440         const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
441
442         if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
443         if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
444         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
445         if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
446
447         if ( hQ != 0x0 ){
448 //          TDirectory *dir = hQ->GetDirectory();
449 //          hQ->SetDirectory(0);
450             TH2S *hNew = new TH2S(*hQ);
451             hNew->SetDirectory(0);
452             fHistoQArray.AddAt(hNew,iSec);
453 //            hQ->SetDirectory(dir);
454         }
455         if ( hT0 != 0x0 ){
456 //          TDirectory *dir = hT0->GetDirectory();
457 //          hT0->SetDirectory(0);
458             TH2S *hNew = new TH2S(*hT0);
459             hNew->SetDirectory(0);
460             fHistoT0Array.AddAt(hNew,iSec);
461 //            hT0->SetDirectory(dir);
462         }
463         if ( hRMS != 0x0 ){
464 //          TDirectory *dir = hRMS->GetDirectory();
465 //          hRMS->SetDirectory(0);
466             TH2S *hNew = new TH2S(*hRMS);
467             hNew->SetDirectory(0);
468             fHistoRMSArray.AddAt(hNew,iSec);
469 //            hRMS->SetDirectory(dir);
470         }
471     }
472
473     //copy fit parameters event by event
474     TObjArray *arr=0x0;
475     for (Int_t iSec=0; iSec<72; ++iSec){
476         arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
477         if ( arr ){
478             TObjArray *arrEvents = new TObjArray(arr->GetSize());
479             fParamArrayEventPol1.AddAt(arrEvents, iSec);
480             for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
481                 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
482                     arrEvents->AddAt(new TVectorD(*vec),iEvent);
483         }
484
485         arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
486         if ( arr ){
487             TObjArray *arrEvents = new TObjArray(arr->GetSize());
488             fParamArrayEventPol2.AddAt(arrEvents, iSec);
489             for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
490                 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
491                     arrEvents->AddAt(new TVectorD(*vec),iEvent);
492         }
493
494         TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
495         TVectorF *vMeanQ    = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
496         if ( vMeanTime )
497             fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
498         if ( vMeanQ )
499             fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
500     }
501
502
503     fVEventTime.ResizeTo(sig.fVEventTime);
504     fVEventNumber.ResizeTo(sig.fVEventNumber);
505     fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
506     fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
507
508 }
509 //_____________________________________________________________________
510 AliTPCCalibCE& AliTPCCalibCE::operator = (const  AliTPCCalibCE &source)
511 {
512   //
513   // assignment operator
514   //
515   if (&source == this) return *this;
516   new (this) AliTPCCalibCE(source);
517
518   return *this;
519 }
520 //_____________________________________________________________________
521 AliTPCCalibCE::~AliTPCCalibCE()
522 {
523     //
524     // destructor
525     //
526
527     fCalRocArrayT0.Delete();
528     fCalRocArrayQ.Delete();
529     fCalRocArrayRMS.Delete();
530     fCalRocArrayOutliers.Delete();
531
532     fHistoQArray.Delete();
533     fHistoT0Array.Delete();
534     fHistoRMSArray.Delete();
535
536     fHistoTmean.Delete();
537
538     fParamArrayEventPol1.Delete();
539     fParamArrayEventPol2.Delete();
540     fTMeanArrayEvent.Delete();
541     fQMeanArrayEvent.Delete();
542
543     fPadTimesArrayEvent.Delete();
544     fPadQArrayEvent.Delete();
545     fPadRMSArrayEvent.Delete();
546     fPadPedestalArrayEvent.Delete();
547
548     if ( fDebugStreamer) delete fDebugStreamer;
549 //    if ( fHTime0 ) delete fHTime0;
550     delete fParam;
551 }
552 //_____________________________________________________________________
553 Int_t AliTPCCalibCE::Update(const Int_t icsector,
554                                 const Int_t icRow,
555                                 const Int_t icPad,
556                                 const Int_t icTimeBin,
557                                 const Float_t csignal)
558 {
559     //
560     // Signal filling methode on the fly pedestal and Time offset correction if necessary.
561     // no extra analysis necessary. Assumes knowledge of the signal shape!
562     // assumes that it is looped over consecutive time bins of one pad
563     //
564     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
565
566     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
567
568     //init first pad and sector in this event
569     if ( fCurrentChannel == -1 ) {
570         fCurrentChannel = iChannel;
571         fCurrentSector  = icsector;
572         fCurrentRow     = icRow;
573     }
574
575     //process last pad if we change to a new one
576     if ( iChannel != fCurrentChannel ){
577         ProcessPad();
578         fCurrentChannel = iChannel;
579         fCurrentSector  = icsector;
580         fCurrentRow     = icRow;
581     }
582
583     //fill signals for current pad
584     fPadSignal.GetMatrixArray()[icTimeBin]=csignal;
585     if ( csignal > fMaxPadSignal ){
586         fMaxPadSignal = csignal;
587         fMaxTimeBin   = icTimeBin;
588     }
589     return 0;
590 }
591 //_____________________________________________________________________
592 void AliTPCCalibCE::FindPedestal(Float_t part)
593 {
594     //
595     // find pedestal and noise for the current pad. Use either database or
596     // truncated mean with part*100%
597     //
598     Bool_t noPedestal = kTRUE;
599
600     //use pedestal database if set
601     if (fPedestalTPC&&fPadNoiseTPC){
602         //only load new pedestals if the sector has changed
603         if ( fCurrentSector!=fLastSector ){
604             fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
605             fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
606             fLastSector=fCurrentSector;
607         }
608
609         if ( fPedestalROC&&fPadNoiseROC ){
610             fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
611             fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
612             noPedestal   = kFALSE;
613         }
614
615     }
616
617     //if we are not running with pedestal database, or for the current sector there is no information
618     //available, calculate the pedestal and noise on the fly
619     if ( noPedestal ) {
620         const Int_t kPedMax = 100;  //maximum pedestal value
621         Float_t  max    =  0;
622         Float_t  maxPos =  0;
623         Int_t    median =  -1;
624         Int_t    count0 =  0;
625         Int_t    count1 =  0;
626         //
627         Float_t padSignal=0;
628         //
629         UShort_t histo[kPedMax];
630         memset(histo,0,kPedMax*sizeof(UShort_t));
631
632         //fill pedestal histogram
633         for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
634             padSignal = fPadSignal.GetMatrixArray()[i];
635             if (padSignal<=0) continue;
636             if (padSignal>max && i>10) {
637                 max = padSignal;
638                 maxPos = i;
639             }
640             if (padSignal>kPedMax-1) continue;
641             histo[int(padSignal+0.5)]++;
642             count0++;
643         }
644         //find median
645         for (Int_t i=1; i<kPedMax; ++i){
646             if (count1<count0*0.5) median=i;
647             count1+=histo[i];
648         }
649         // truncated mean
650         //
651         Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
652         //
653         for (Int_t idelta=1; idelta<10; ++idelta){
654             if (median-idelta<=0) continue;
655             if (median+idelta>kPedMax) continue;
656             if (count<part*count1){
657                 count+=histo[median-idelta];
658                 mean +=histo[median-idelta]*(median-idelta);
659                 rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
660                 count+=histo[median+idelta];
661                 mean +=histo[median+idelta]*(median+idelta);
662                 rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
663             }
664         }
665         fPadPedestal = 0;
666         fPadNoise    = 0;
667         if ( count > 0 ) {
668             mean/=count;
669             rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
670             fPadPedestal = mean;
671             fPadNoise    = rms;
672         } 
673     }
674 }
675 //_____________________________________________________________________
676 void AliTPCCalibCE::FindCESignal(TVectorD &param, Float_t &qSum, const TVectorF maxima)
677 {
678     //
679     //  Find position, signal width and height of the CE signal (last signal)
680     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
681     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
682     //
683
684     Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
685     Int_t   cemaxpos       = 0;
686     Float_t ceSumThreshold = 8.*fPadNoise;  // threshold for the signal sum
687     const Int_t    kCemin  = 4;             // range for the analysis of the ce signal +- channels from the peak
688     const Int_t    kCemax  = 7;
689
690     Float_t minDist  = 25;  //initial minimum distance betweek roc mean ce signal and pad ce signal
691
692     // find maximum closest to the sector mean from the last event
693     for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
694         // get sector mean of last event
695         Float_t tmean = (*GetTMeanEvents(fCurrentSector))[fNevents-1];
696             if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
697                 minDist  = tmean-maxima[imax];
698                 cemaxpos = (Int_t)maxima[imax];
699             }
700     }
701
702     if (cemaxpos!=0){
703         ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
704         for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
705             if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
706                 Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
707                 if (signal>0) {
708                     ceTime+=signal*(i+0.5);
709                     ceRMS +=signal*(i+0.5)*(i+0.5);
710                     ceQsum+=signal;
711                 }
712             }
713         }
714     }
715     if (ceQmax&&ceQsum>ceSumThreshold) {
716         ceTime/=ceQsum;
717         ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
718         fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
719         fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
720
721         //Normalise Q to pad area of irocs
722         Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
723
724         ceQsum/=norm;
725         fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
726         fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
727     } else {
728         ceQmax=0;
729         ceTime=0;
730         ceRMS =0;
731         ceQsum=0;
732     }
733     param[0] = ceQmax;
734     param[1] = ceTime;
735     param[2] = ceRMS;
736     qSum     = ceQsum;
737 }
738 //_____________________________________________________________________
739 Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus)
740 {
741     //
742     // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
743     // and 'tplus' timebins after 'pos'
744     //
745     if ( (pos-tminus)<fFirstTimeBin || (pos+tplus)>fLastTimeBin ) return kFALSE;
746     for (Int_t iTime = pos; iTime>pos-tminus; --iTime)
747         if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
748     for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; ++iTime, ++iTime2){
749         if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
750             ++iTime2;
751         if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
752     }
753     return kTRUE;
754 }
755 //_____________________________________________________________________
756 void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
757 {
758     //
759     // Find local maxima on the pad signal and Histogram them
760     //
761     Float_t ceThreshold = 5.*fPadNoise;  // threshold for the signal
762     Int_t   count       = 0;
763     Int_t   tminus      = 2;
764     Int_t   tplus       = 3;
765     for (Int_t i=fLastTimeBin-tplus-1; i>=fFirstTimeBin+tminus; --i){
766         if ( (fPadSignal[i]-fPadPedestal)>ceThreshold && IsPeak(i,tminus,tplus) ){
767           if (count<maxima.GetNrows()){
768             maxima.GetMatrixArray()[count++]=i;
769             GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
770           }
771         }
772     }
773 }
774 //_____________________________________________________________________
775 void AliTPCCalibCE::ProcessPad()
776 {
777     //
778     //  Process data of current pad
779     //
780     FindPedestal();
781
782     TVectorF maxima(15);     // the expected maximum number of maxima in the complete TPC should be 8 laser beam layers
783                              // + central electrode and possibly post peaks from the CE signal
784                              // however if we are on a high noise pad a lot more peaks due to the noise might occur
785     FindLocalMaxima(maxima);
786     if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return;  // return because we don't have Time0 info for the CE yet
787
788
789
790     TVectorD param(3);
791     Float_t  Qsum;
792     FindCESignal(param, Qsum, maxima);
793
794     Double_t meanT  = param[1];
795     Double_t sigmaT = param[2];
796
797     //Fill Event T0 counter
798     (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
799
800     //Fill Q histogram
801     GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(Qsum), fCurrentChannel );
802
803     //Fill RMS histogram
804     GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
805
806
807     //Fill debugging info
808     if ( fDebugLevel>0 ){
809         (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
810         (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
811         (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=Qsum;
812     }
813
814     ResetPad();
815 }
816 //_____________________________________________________________________
817 void AliTPCCalibCE::EndEvent()
818 {
819     //
820     //  Process data of current pad
821     //  The Functions 'SetTimeStamp' and 'SetRunNumber'  should be called
822     //  before the EndEvent function to set the event timestamp and number!!!
823     //  This is automatically done if the ProcessEvent(AliRawReader *rawReader)
824     //  function was called
825     //
826
827     //check if last pad has allready been processed, if not do so
828     if ( fMaxTimeBin>-1 ) ProcessPad();
829
830     TVectorD param(3);
831     TMatrixD dummy(3,3);
832 //    TVectorF vMeanTime(72);
833 //    TVectorF vMeanQ(72);
834     AliTPCCalROC *calIroc=new AliTPCCalROC(0);
835     AliTPCCalROC *calOroc=new AliTPCCalROC(36);
836
837     //find mean time0 offset for side A and C
838     Double_t time0Side[2];       //time0 for side A:0 and C:0
839     Double_t time0SideCount[2];  //time0 counter for side A:0 and C:0
840     time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
841     for ( Int_t iSec = 0; iSec<72; ++iSec ){
842         time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
843         time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
844     }
845     if ( time0SideCount[0] >0  )
846         time0Side[0]/=time0SideCount[0];
847     if ( time0SideCount[1] >0 )
848         time0Side[1]/=time0SideCount[1];
849     // end find time0 offset
850
851     //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
852     for ( Int_t iSec = 0; iSec<72; ++iSec ){
853       //find median and then calculate the mean around it
854         TH1S *hMeanT    = GetHistoTmean(iSec); //histogram with local maxima position information
855         if ( !hMeanT ) continue;
856
857         Double_t entries = hMeanT->GetEntries();
858         Double_t sum     = 0;
859         Short_t *arr     = hMeanT->GetArray()+1;
860         Int_t ibin=0;
861         for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
862             sum+=arr[ibin];
863             if ( sum>=(entries/2) ) break;
864         }
865         Int_t delta = 4;
866         Int_t firstBin = fFirstTimeBin+ibin-delta;
867         Int_t lastBin  = fFirstTimeBin+ibin+delta;
868         if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
869         if ( lastBin>fLastTimeBin   ) lastBin =fLastTimeBin;
870         Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
871
872         // check boundaries for ebye info of mean time
873         TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
874         Int_t vSize=vMeanTime->GetNrows();
875         if ( vSize < fNevents+1 )
876             vMeanTime->ResizeTo(vSize+100);
877
878         vMeanTime->GetMatrixArray()[fNevents]=median;
879       // end find median
880
881         TVectorF *vTimes = GetPadTimesEvent(iSec);
882         if ( !vTimes ) continue;                     //continue if no time information for this sector is available
883
884
885         AliTPCCalROC calIrocOutliers(0);
886         AliTPCCalROC calOrocOutliers(36);
887
888         // calculate mean Q of the sector
889         Float_t meanQ = 0;
890         if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
891         TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
892         if ( vSize < fNevents+1 )           // vSize is the same as for vMeanTime!
893             vMeanQ->ResizeTo(vSize+100);
894
895         vMeanQ->GetMatrixArray()[fNevents]=meanQ;
896
897         for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
898             Float_t Time  = (*vTimes).GetMatrixArray()[iChannel];
899
900             //set values for temporary roc calibration class
901             if ( iSec < 36 ) {
902                 calIroc->SetValue(iChannel, Time);
903                 if ( Time == 0 ) calIrocOutliers.SetValue(iChannel,1);
904
905             } else {
906                 calOroc->SetValue(iChannel, Time);
907                 if ( Time == 0 ) calOrocOutliers.SetValue(iChannel,1);
908             }
909
910             if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
911                 GetHistoT0(iSec,kTRUE)->Fill( Time-time0Side[(iSec/18)%2],iChannel );
912
913
914
915             //-------------------------------  Debug start  ------------------------------
916             if ( fDebugLevel>0 ){
917                 if ( !fDebugStreamer ) {
918                         //debug stream
919                     TDirectory *backup = gDirectory;
920                     fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
921                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
922                 }
923
924                 Int_t row=0;
925                 Int_t pad=0;
926                 Int_t padc=0;
927
928                 Float_t Q   = (*GetPadQEvent(iSec))[iChannel];
929                 Float_t RMS = (*GetPadRMSEvent(iSec))[iChannel];
930
931                 UInt_t channel=iChannel;
932                 Int_t sector=iSec;
933
934                 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
935                 pad = channel-fROC->GetRowIndexes(sector)[row];
936                 padc = pad-(fROC->GetNPads(sector,row)/2);
937
938 //              TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
939 //                                  Form("hSignalD%d.%d.%d",sector,row,pad),
940 //                                  fLastTimeBin-fFirstTimeBin,
941 //                                  fFirstTimeBin,fLastTimeBin);
942 //              h1->SetDirectory(0);
943 //
944 //              for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
945 //                  h1->Fill(i,fPadSignal(i));
946
947                 Double_t T0Sec = 0;
948                 if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
949                     T0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
950                 Double_t T0Side = time0Side[(iSec/18)%2];
951                 (*fDebugStreamer) << "DataPad" <<
952                     "Event=" << fNevents <<
953                     "RunNumber=" << fRunNumber <<
954                     "TimeStamp="   << fTimeStamp <<
955                     "Sector="<< sector <<
956                     "Row="   << row<<
957                     "Pad="   << pad <<
958                     "PadC="  << padc <<
959                     "PadSec="<< channel <<
960                     "Time0Sec="  << T0Sec <<
961                     "Time0Side=" << T0Side <<
962                     "Time="  << Time <<
963                     "RMS="   << RMS <<
964                     "Sum="   << Q <<
965                     "MeanQ=" << meanQ <<
966                     //              "hist.=" << h1 <<
967                     "\n";
968
969                 //              delete h1;
970
971             }
972             //-----------------------------  Debug end  ------------------------------
973         }// end channel loop
974         hMeanT->Reset();
975
976         TVectorD paramPol1(3);
977         TVectorD paramPol2(6);
978         TMatrixD matPol1(3,3);
979         TMatrixD matPol2(6,6);
980         Float_t  chi2Pol1=0;
981         Float_t  chi2Pol2=0;
982
983         if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
984             if ( iSec < 36 ){
985                 calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
986                 calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
987             } else {
988                 calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
989                 calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
990             }
991
992             GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
993             GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
994         }
995 //      printf("events: %d -- size: %d\n",fNevents,GetParamArrayPol1(iSec)->GetSize());
996
997         //-------------------------------  Debug start  ------------------------------
998         if ( fDebugLevel>0 ){
999             if ( !fDebugStreamer ) {
1000                 //debug stream
1001                 TDirectory *backup = gDirectory;
1002                 fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
1003                 if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1004             }
1005             (*fDebugStreamer) << "DataRoc" <<
1006                 "Event=" << fEvent <<
1007                 "RunNumber=" << fRunNumber <<
1008                 "TimeStamp="   << fTimeStamp <<
1009                 "Sector="<< iSec <<
1010                 "hMeanT.=" << hMeanT <<
1011                 "median=" << median <<
1012                 "paramPol1.=" << &paramPol1 <<
1013                 "paramPol2.=" << &paramPol2 <<
1014                 "matPol1.="   << &matPol1 <<
1015                 "matPol2.="   << &matPol2 <<
1016                 "chi2Pol1="   << chi2Pol1 <<
1017                 "chi2Pol2="   << chi2Pol2 <<
1018                 "\n";
1019         }
1020         //-------------------------------  Debug end  ------------------------------
1021     }// end sector loop
1022
1023 //    fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1024 //    fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
1025     if ( fVEventTime.GetNrows() < fNevents+1 ) {
1026         fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1027         fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
1028     }
1029     fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1030     fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
1031
1032     fNevents++;
1033     fOldRunNumber = fRunNumber;
1034
1035     delete calIroc;
1036     delete calOroc;
1037 }
1038 //_____________________________________________________________________
1039 Bool_t AliTPCCalibCE::ProcessEvent(AliTPCRawStream *rawStream)
1040 {
1041   //
1042   // Event Processing loop - AliTPCRawStream
1043   // The Function 'SetTimeStamp' should be called for each event to set the event time stamp!!!
1044   //
1045
1046   rawStream->SetOldRCUFormat(fOldRCUformat);
1047
1048   ResetEvent();
1049
1050   Bool_t withInput = kFALSE;
1051
1052   while (rawStream->Next()) {
1053
1054       Int_t isector  = rawStream->GetSector();                       //  current sector
1055       Int_t iRow     = rawStream->GetRow();                          //  current row
1056       Int_t iPad     = rawStream->GetPad();                          //  current pad
1057       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
1058       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
1059
1060       Update(isector,iRow,iPad,iTimeBin,signal);
1061       withInput = kTRUE;
1062   }
1063
1064   if (withInput){
1065       EndEvent();
1066   }
1067
1068   return withInput;
1069 }
1070 //_____________________________________________________________________
1071 Bool_t AliTPCCalibCE::ProcessEvent(AliRawReader *rawReader)
1072 {
1073   //
1074   //  Event processing loop - AliRawReader
1075   //
1076
1077
1078     AliTPCRawStream rawStream(rawReader);
1079     AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1080     if (eventHeader){
1081         fTimeStamp   = eventHeader->Get("Timestamp");
1082         fRunNumber = eventHeader->Get("RunNb");
1083     }
1084     fEventId = *rawReader->GetEventId();
1085
1086
1087     rawReader->Select("TPC");
1088
1089   return ProcessEvent(&rawStream);
1090 }
1091 //_____________________________________________________________________
1092 Bool_t AliTPCCalibCE::ProcessEvent(eventHeaderStruct *event)
1093 {
1094   //
1095   //  Event processing loop - date event
1096   //
1097     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
1098     Bool_t result=ProcessEvent(rawReader);
1099     delete rawReader;
1100     return result;
1101
1102 }
1103 //_____________________________________________________________________
1104 TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1105                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
1106                                   Char_t *type, Bool_t force)
1107 {
1108     //
1109     // return pointer to TH2S histogram of 'type'
1110     // if force is true create a new histogram if it doesn't exist allready
1111     //
1112     if ( !force || arr->UncheckedAt(sector) )
1113         return (TH2S*)arr->UncheckedAt(sector);
1114
1115     // if we are forced and histogram doesn't exist yet create it
1116     Char_t name[255], title[255];
1117
1118     sprintf(name,"hCalib%s%.2d",type,sector);
1119     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1120
1121     // new histogram with Q calib information. One value for each pad!
1122     TH2S* hist = new TH2S(name,title,
1123                           nbinsY, ymin, ymax,
1124                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1125     hist->SetDirectory(0);
1126     arr->AddAt(hist,sector);
1127     return hist;
1128 }
1129 //_____________________________________________________________________
1130 TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
1131 {
1132     //
1133     // return pointer to T0 histogram
1134     // if force is true create a new histogram if it doesn't exist allready
1135     //
1136     TObjArray *arr = &fHistoT0Array;
1137     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1138 }
1139 //_____________________________________________________________________
1140 TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
1141 {
1142     //
1143     // return pointer to Q histogram
1144     // if force is true create a new histogram if it doesn't exist allready
1145     //
1146     TObjArray *arr = &fHistoQArray;
1147     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1148 }
1149 //_____________________________________________________________________
1150 TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
1151 {
1152     //
1153     // return pointer to Q histogram
1154     // if force is true create a new histogram if it doesn't exist allready
1155     //
1156     TObjArray *arr = &fHistoRMSArray;
1157     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1158 }
1159 //_____________________________________________________________________
1160 TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1161                               Char_t *type, Bool_t force)
1162 {
1163     //
1164     // return pointer to TH1S histogram
1165     // if force is true create a new histogram if it doesn't exist allready
1166     //
1167     if ( !force || arr->UncheckedAt(sector) )
1168         return (TH1S*)arr->UncheckedAt(sector);
1169
1170     // if we are forced and histogram doesn't yes exist create it
1171     Char_t name[255], title[255];
1172
1173     sprintf(name,"hCalib%s%.2d",type,sector);
1174     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1175
1176     // new histogram with Q calib information. One value for each pad!
1177     TH1S* hist = new TH1S(name,title,
1178                           fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1179     hist->SetDirectory(0);
1180     arr->AddAt(hist,sector);
1181     return hist;
1182 }
1183 //_____________________________________________________________________
1184 TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1185 {
1186     //
1187     // return pointer to Q histogram
1188     // if force is true create a new histogram if it doesn't exist allready
1189     //
1190     TObjArray *arr = &fHistoTmean;
1191     return GetHisto(sector, arr, "LastTmean", force);
1192 }
1193 //_____________________________________________________________________
1194 TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force)
1195 {
1196     //
1197     // return pointer to Pad Info from 'arr' for the current event and sector
1198     // if force is true create it if it doesn't exist allready
1199     //
1200     if ( !force || arr->UncheckedAt(sector) )
1201         return (TVectorF*)arr->UncheckedAt(sector);
1202
1203     TVectorF *vect = new TVectorF(size);
1204     arr->AddAt(vect,sector);
1205     return vect;
1206 }
1207 //_____________________________________________________________________
1208 TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
1209 {
1210     //
1211     // return pointer to Pad Times Array for the current event and sector
1212     // if force is true create it if it doesn't exist allready
1213     //
1214     TObjArray *arr = &fPadTimesArrayEvent;
1215     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1216 }
1217 //_____________________________________________________________________
1218 TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
1219 {
1220     //
1221     // return pointer to Pad Q Array for the current event and sector
1222     // if force is true create it if it doesn't exist allready
1223     // for debugging purposes only
1224     //
1225
1226     TObjArray *arr = &fPadQArrayEvent;
1227     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1228 }
1229 //_____________________________________________________________________
1230 TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
1231 {
1232     //
1233     // return pointer to Pad RMS Array for the current event and sector
1234     // if force is true create it if it doesn't exist allready
1235     // for debugging purposes only
1236     //
1237     TObjArray *arr = &fPadRMSArrayEvent;
1238     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1239 }
1240 //_____________________________________________________________________
1241 TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
1242 {
1243     //
1244     // return pointer to Pad RMS Array for the current event and sector
1245     // if force is true create it if it doesn't exist allready
1246     // for debugging purposes only
1247     //
1248     TObjArray *arr = &fPadPedestalArrayEvent;
1249     return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1250 }
1251 //_____________________________________________________________________
1252 TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1253 {
1254     //
1255     // return pointer to the EbyE info of the mean arrival time for 'sector'
1256     // if force is true create it if it doesn't exist allready
1257     //
1258     TObjArray *arr = &fTMeanArrayEvent;
1259     return GetVectSector(sector,arr,100,force);
1260 }
1261 //_____________________________________________________________________
1262 TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1263 {
1264     //
1265     // return pointer to the EbyE info of the mean arrival time for 'sector'
1266     // if force is true create it if it doesn't exist allready
1267     //
1268     TObjArray *arr = &fQMeanArrayEvent;
1269     return GetVectSector(sector,arr,100,force);
1270 }
1271 //_____________________________________________________________________
1272 AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force)
1273 {
1274     //
1275     // return pointer to ROC Calibration
1276     // if force is true create a new histogram if it doesn't exist allready
1277     //
1278     if ( !force || arr->UncheckedAt(sector) )
1279         return (AliTPCCalROC*)arr->UncheckedAt(sector);
1280
1281     // if we are forced and histogram doesn't yes exist create it
1282
1283     // new AliTPCCalROC for T0 information. One value for each pad!
1284     AliTPCCalROC *croc = new AliTPCCalROC(sector);
1285     arr->AddAt(croc,sector);
1286     return croc;
1287 }
1288 //_____________________________________________________________________
1289 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
1290 {
1291     //
1292     // return pointer to Carge ROC Calibration
1293     // if force is true create a new histogram if it doesn't exist allready
1294     //
1295     TObjArray *arr = &fCalRocArrayT0;
1296     return GetCalRoc(sector, arr, force);
1297 }
1298 //_____________________________________________________________________
1299 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
1300 {
1301     //
1302     // return pointer to T0 ROC Calibration
1303     // if force is true create a new histogram if it doesn't exist allready
1304     //
1305     TObjArray *arr = &fCalRocArrayQ;
1306     return GetCalRoc(sector, arr, force);
1307 }
1308 //_____________________________________________________________________
1309 AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
1310 {
1311     //
1312     // return pointer to signal width ROC Calibration
1313     // if force is true create a new histogram if it doesn't exist allready
1314     //
1315     TObjArray *arr = &fCalRocArrayRMS;
1316     return GetCalRoc(sector, arr, force);
1317 }
1318 //_____________________________________________________________________
1319 AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1320 {
1321     //
1322     // return pointer to Outliers
1323     // if force is true create a new histogram if it doesn't exist allready
1324     //
1325     TObjArray *arr = &fCalRocArrayOutliers;
1326     return GetCalRoc(sector, arr, force);
1327 }
1328 //_____________________________________________________________________
1329 TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force)
1330 {
1331     //
1332     // return pointer to TObjArray of fit parameters
1333     // if force is true create a new histogram if it doesn't exist allready
1334     //
1335     if ( !force || arr->UncheckedAt(sector) )
1336         return (TObjArray*)arr->UncheckedAt(sector);
1337
1338     // if we are forced and array doesn't yes exist create it
1339
1340     // new TObjArray for parameters
1341     TObjArray *newArr = new TObjArray;
1342     arr->AddAt(newArr,sector);
1343     return newArr;
1344 }
1345 //_____________________________________________________________________
1346 TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1347 {
1348     //
1349     // return pointer to TObjArray of fit parameters from plane fit
1350     // if force is true create a new histogram if it doesn't exist allready
1351     //
1352     TObjArray *arr = &fParamArrayEventPol1;
1353     return GetParamArray(sector, arr, force);
1354 }
1355 //_____________________________________________________________________
1356 TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1357 {
1358     //
1359     // return pointer to TObjArray of fit parameters from parabola fit
1360     // if force is true create a new histogram if it doesn't exist allready
1361     //
1362     TObjArray *arr = &fParamArrayEventPol2;
1363     return GetParamArray(sector, arr, force);
1364 }
1365 //_____________________________________________________________________
1366 void AliTPCCalibCE::ResetEvent()
1367 {
1368     //
1369     //  Reset global counters  -- Should be called before each event is processed
1370     //
1371     fLastSector=-1;
1372     fCurrentSector=-1;
1373     fCurrentRow=-1;
1374     fCurrentChannel=-1;
1375
1376     ResetPad();
1377
1378     fPadTimesArrayEvent.Delete();
1379     fPadQArrayEvent.Delete();
1380     fPadRMSArrayEvent.Delete();
1381     fPadPedestalArrayEvent.Delete();
1382
1383     for ( Int_t i=0; i<72; ++i ){
1384         fVTime0Offset.GetMatrixArray()[i]=0;
1385         fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1386         fVMeanQ.GetMatrixArray()[i]=0;
1387         fVMeanQCounter.GetMatrixArray()[i]=0;
1388     }
1389 }
1390 //_____________________________________________________________________
1391 void AliTPCCalibCE::ResetPad()
1392 {
1393     //
1394     //  Reset pad infos -- Should be called after a pad has been processed
1395     //
1396     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1397         fPadSignal.GetMatrixArray()[i] = 0;
1398     fMaxTimeBin   = -1;
1399     fMaxPadSignal = -1;
1400     fPadPedestal  = -1;
1401     fPadNoise     = -1;
1402 }
1403 //_____________________________________________________________________
1404 void AliTPCCalibCE::Merge(AliTPCCalibCE *ce)
1405 {
1406     //
1407     //  Merge ce to the current AliTPCCalibCE
1408     //
1409
1410     //merge histograms
1411     for (Int_t iSec=0; iSec<72; ++iSec){
1412         TH2S *hRefQmerge   = ce->GetHistoQ(iSec);
1413         TH2S *hRefT0merge  = ce->GetHistoT0(iSec);
1414         TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1415
1416
1417         if ( hRefQmerge ){
1418             TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1419             TH2S *hRefQ   = GetHistoQ(iSec);
1420             if ( hRefQ ) hRefQ->Add(hRefQmerge);
1421             else {
1422                 TH2S *hist = new TH2S(*hRefQmerge);
1423                 hist->SetDirectory(0);
1424                 fHistoQArray.AddAt(hist, iSec);
1425             }
1426             hRefQmerge->SetDirectory(dir);
1427         }
1428         if ( hRefT0merge ){
1429             TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1430             TH2S *hRefT0  = GetHistoT0(iSec);
1431             if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1432             else {
1433                 TH2S *hist = new TH2S(*hRefT0merge);
1434                 hist->SetDirectory(0);
1435                 fHistoT0Array.AddAt(hist, iSec);
1436             }
1437             hRefT0merge->SetDirectory(dir);
1438         }
1439         if ( hRefRMSmerge ){
1440             TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1441             TH2S *hRefRMS = GetHistoRMS(iSec);
1442             if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1443             else {
1444                 TH2S *hist = new TH2S(*hRefRMSmerge);
1445                 hist->SetDirectory(0);
1446                 fHistoRMSArray.AddAt(hist, iSec);
1447             }
1448             hRefRMSmerge->SetDirectory(dir);
1449         }
1450
1451     }
1452
1453     // merge time information
1454
1455
1456     Int_t nCEevents = ce->GetNeventsProcessed();
1457     for (Int_t iSec=0; iSec<72; ++iSec){
1458         TObjArray *arrPol1CE  = ce->GetParamArrayPol1(iSec);
1459         TObjArray *arrPol2CE  = ce->GetParamArrayPol2(iSec);
1460         TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1461         TVectorF *vMeanQCE    = ce->GetQMeanEvents(iSec);
1462
1463         TObjArray *arrPol1  = 0x0;
1464         TObjArray *arrPol2  = 0x0;
1465         TVectorF *vMeanTime = 0x0;
1466         TVectorF *vMeanQ    = 0x0;
1467
1468         //resize arrays
1469         if ( arrPol1CE && arrPol2CE ){
1470             arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1471             arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1472             arrPol1->Expand(fNevents+nCEevents);
1473             arrPol2->Expand(fNevents+nCEevents);
1474         }
1475         if ( vMeanTimeCE && vMeanQCE ){
1476             vMeanTime = GetTMeanEvents(iSec);
1477             vMeanQCE  = GetQMeanEvents(iSec);
1478             vMeanTime->ResizeTo(fNevents+nCEevents);
1479             vMeanQ->ResizeTo(fNevents+nCEevents);
1480         }
1481
1482
1483         for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1484             if ( arrPol1CE && arrPol2CE ){
1485                 TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1486                 TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1487                 if ( paramPol1 && paramPol2 ){
1488                     GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1489                     GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1490                 }
1491             }
1492             if ( vMeanTimeCE && vMeanQCE ){
1493                 vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1494                 vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1495             }
1496         }
1497     }
1498
1499
1500
1501     TVectorD*  eventTimes  = ce->GetEventTimes();
1502     TVectorD*  eventIds  = ce->GetEventIds();
1503     fVEventTime.ResizeTo(fNevents+nCEevents);
1504     fVEventNumber.ResizeTo(fNevents+nCEevents);
1505
1506     for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1507         Double_t evTime     = eventTimes->GetMatrixArray()[iEvent];
1508         Double_t evId       = eventIds->GetMatrixArray()[iEvent];
1509         fVEventTime.GetMatrixArray()[fNevents+iEvent] = evTime;
1510         fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1511     }
1512     fNevents+=nCEevents; //increase event counter
1513
1514 }
1515 //_____________________________________________________________________
1516 TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
1517 {
1518     //
1519     // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
1520     // xVariable:    0-event time, 1-event id, 2-internal event counter
1521     // fitType:      0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1522     // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1523     //                               0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1524     //                               not used for mean time and mean Q )
1525     // for an example see class description at the beginning
1526     //
1527
1528     Double_t *x = new Double_t[fNevents];
1529     Double_t *y = new Double_t[fNevents];
1530
1531     TVectorD *xVar = 0x0;
1532     TObjArray *aType = 0x0;
1533     Int_t npoints=0;
1534
1535     // sanity checks
1536     if ( (sector<0) || (sector>71) )      return 0x0;
1537     if ( (xVariable<0) || (xVariable>2) ) return 0x0;
1538     if ( (fitType<0) || (fitType>3) )     return 0x0;
1539     if ( fitType==0 ){
1540         if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1541         aType = &fParamArrayEventPol1;
1542         if ( aType->At(sector)==0x0 ) return 0x0;
1543     }
1544     else if ( fitType==1 ){
1545         if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1546         aType = &fParamArrayEventPol2;
1547         if ( aType->At(sector)==0x0 ) return 0x0;
1548     }
1549
1550
1551     if ( xVariable == 0 ) xVar = &fVEventTime;
1552     if ( xVariable == 1 ) xVar = &fVEventNumber;
1553     if ( xVariable == 2 ) {
1554         xVar = new TVectorD(fNevents);
1555         for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1556     }
1557
1558     for (Int_t ievent =0; ievent<fNevents; ++ievent){
1559         if ( fitType<2 ){
1560             TObjArray *events = (TObjArray*)(aType->At(sector));
1561             if ( events->GetSize()<=ievent ) break;
1562             TVectorD *v = (TVectorD*)(events->At(ievent));
1563             if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1564         } else if (fitType == 2) {
1565             Double_t xValue=(*xVar)[ievent];
1566             Double_t yValue=(*GetTMeanEvents(sector))[ievent];
1567             if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1568         }else if (fitType == 3) {
1569             Double_t xValue=(*xVar)[ievent];
1570             Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1571             if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1572         }
1573     }
1574
1575     TGraph *gr = new TGraph(npoints);
1576     //sort xVariable increasing
1577     Int_t    *sortIndex = new Int_t[npoints];
1578     TMath::Sort(npoints,x,sortIndex);
1579     for (Int_t i=0;i<npoints;++i){
1580         gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1581     }
1582
1583
1584     if ( xVariable == 2 ) delete xVar;
1585     delete x;
1586     delete y;
1587     delete sortIndex;
1588     return gr;
1589 }
1590 //_____________________________________________________________________
1591 void AliTPCCalibCE::Analyse()
1592 {
1593     //
1594     //  Calculate calibration constants
1595     //
1596
1597     TVectorD paramQ(3);
1598     TVectorD paramT0(3);
1599     TVectorD paramRMS(3);
1600     TMatrixD dummy(3,3);
1601
1602     for (Int_t iSec=0; iSec<72; ++iSec){
1603         TH2S *hT0 = GetHistoT0(iSec);
1604         if (!hT0 ) continue;
1605
1606         AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
1607         AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
1608         AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1609         AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1610
1611         TH2S *hQ   = GetHistoQ(iSec);
1612         TH2S *hRMS = GetHistoRMS(iSec);
1613
1614         Short_t *array_hQ   = hQ->GetArray();
1615         Short_t *array_hT0  = hT0->GetArray();
1616         Short_t *array_hRMS = hRMS->GetArray();
1617
1618         UInt_t nChannels = fROC->GetNChannels(iSec);
1619
1620         //debug
1621         Int_t row=0;
1622         Int_t pad=0;
1623         Int_t padc=0;
1624         //! debug
1625
1626         for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1627
1628
1629             Float_t cogTime0 = -1000;
1630             Float_t cogQ     = -1000;
1631             Float_t cogRMS   = -1000;
1632             Float_t cogOut   = 0;
1633
1634
1635             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1636             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1637             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1638
1639             cogQ     = AliMathBase::GetCOG(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1640             cogTime0 = AliMathBase::GetCOG(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1641             cogRMS   = AliMathBase::GetCOG(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1642
1643
1644
1645             /*
1646              //outlier specifications
1647             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1648                 cogOut = 1;
1649                 cogTime0 = 0;
1650                 cogQ     = 0;
1651                 cogRMS   = 0;
1652             }
1653 */
1654             rocQ->SetValue(iChannel, cogQ*cogQ);
1655             rocT0->SetValue(iChannel, cogTime0);
1656             rocRMS->SetValue(iChannel, cogRMS);
1657             rocOut->SetValue(iChannel, cogOut);
1658
1659
1660             //debug
1661             if ( fDebugLevel > 0 ){
1662                 if ( !fDebugStreamer ) {
1663                         //debug stream
1664                     TDirectory *backup = gDirectory;
1665                     fDebugStreamer = new TTreeSRedirector("debugCalibCEAnalysis.root");
1666                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1667                 }
1668
1669                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1670                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1671                 padc = pad-(fROC->GetNPads(iSec,row)/2);
1672
1673                 (*fDebugStreamer) << "DataEnd" <<
1674                     "Sector="  << iSec      <<
1675                     "Pad="     << pad       <<
1676                     "PadC="    << padc      <<
1677                     "Row="     << row       <<
1678                     "PadSec="  << iChannel   <<
1679                     "Q="       << cogQ      <<
1680                     "T0="      << cogTime0  <<
1681                     "RMS="     << cogRMS    <<
1682                     "\n";
1683             }
1684             //! debug
1685
1686         }
1687
1688     }
1689     if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
1690 //    delete fDebugStreamer;
1691 //    fDebugStreamer = 0x0;
1692 }
1693 //_____________________________________________________________________
1694 void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1695 {
1696     //
1697     //  Write class to file
1698     //
1699
1700     TString sDir(dir);
1701     TString option;
1702
1703     if ( append )
1704         option = "update";
1705     else
1706         option = "recreate";
1707
1708     TDirectory *backup = gDirectory;
1709     TFile f(filename,option.Data());
1710     f.cd();
1711     if ( !sDir.IsNull() ){
1712         f.mkdir(sDir.Data());
1713         f.cd(sDir);
1714     }
1715     this->Write();
1716     f.Close();
1717
1718     if ( backup ) backup->cd();
1719 }