0b648e9e1acfa6658a007f50d5a8c5acde4dc36e
[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 //          Implementation of the TPC pulser calibration
22 //
23 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
24 // 
25 // 
26 //-------------------------------------------------------
27
28
29 /* $Id$ */
30
31
32
33 //Root includes
34 #include <TObjArray.h>
35 #include <TH1F.h>
36 #include <TH2S.h>
37 #include <TString.h>
38 #include <TVectorF.h>
39 #include <TVectorD.h>
40 #include <TMatrixD.h>
41 #include <TMath.h>
42 #include <TGraph.h>
43
44 #include <TDirectory.h>
45 #include <TSystem.h>
46 #include <TFile.h>
47
48 //AliRoot includes
49 #include "AliRawReader.h"
50 #include "AliRawReaderRoot.h"
51 #include "AliRawReaderDate.h"
52 #include "AliRawEventHeaderBase.h"
53 #include "AliTPCRawStream.h"
54 #include "AliTPCcalibDB.h"
55 #include "AliTPCCalROC.h"
56 #include "AliTPCCalPad.h"
57 #include "AliTPCROC.h"
58 #include "AliTPCParam.h"
59 #include "AliTPCCalibCE.h"
60 #include "AliMathBase.h"
61 #include "TTreeStream.h"
62
63 //date
64 #include "event.h"
65 ClassImp(AliTPCCalibCE) /*FOLD00*/
66
67 AliTPCCalibCE::AliTPCCalibCE() : /*FOLD00*/
68     TObject(),
69     fFirstTimeBin(650),
70     fLastTimeBin(1000),
71     fNbinsT0(200),
72     fXminT0(-5),
73     fXmaxT0(5),
74     fNbinsQ(200),
75     fXminQ(1),
76     fXmaxQ(40),
77     fNbinsRMS(100),
78     fXminRMS(0.1),
79     fXmaxRMS(5.1),
80     fLastSector(-1),
81     fOldRCUformat(kTRUE),
82     fROC(AliTPCROC::Instance()),
83     fParam(new AliTPCParam),
84     fPedestalTPC(0x0),
85     fPadNoiseTPC(0x0),
86     fPedestalROC(0x0),
87     fPadNoiseROC(0x0),
88     fCalRocArrayT0(72),
89     fCalRocArrayQ(72),
90     fCalRocArrayRMS(72),
91     fCalRocArrayOutliers(72),
92     fHistoQArray(72),
93     fHistoT0Array(72),
94     fHistoRMSArray(72),
95     fHistoTmean(72),
96     fParamArrayEvent(1000),
97     fParamArrayEventPol1(72),
98     fParamArrayEventPol2(72),
99     fTMeanArrayEvent(1000),
100     fQMeanArrayEvent(1000),
101     fVEventTime(1000),
102     fVEventNumber(1000),
103     fNevents(0),
104     fTimeStamp(0),
105     fRunNumber(-1),
106     fOldRunNumber(-1),
107     fPadTimesArrayEvent(72),
108     fPadQArrayEvent(72),
109     fPadRMSArrayEvent(72),
110     fPadPedestalArrayEvent(72),
111     fCurrentChannel(-1),
112     fCurrentSector(-1),
113     fCurrentRow(-1),
114     fMaxPadSignal(-1),
115     fMaxTimeBin(-1),
116     fPadSignal(1024),
117     fPadPedestal(0),
118     fPadNoise(0),
119     fVTime0Offset(72),
120     fVTime0OffsetCounter(72),
121     fVMeanQ(72),
122     fVMeanQCounter(72),
123 //    fHTime0(0x0),
124     fEvent(-1),
125     fDebugStreamer(0x0),
126     fDebugLevel(0)
127 {
128     //
129     // AliTPCSignal default constructor
130     //
131 //    fHTime0 = new TH1F("hTime0Event","hTime0Event",(fLastTimeBin-fFirstTimeBin)*10,fFirstTimeBin,fLastTimeBin);
132 }
133 //_____________________________________________________________________
134 AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
135     TObject(sig),
136     fFirstTimeBin(sig.fFirstTimeBin),
137     fLastTimeBin(sig.fLastTimeBin),
138     fNbinsT0(sig.fNbinsT0),
139     fXminT0(sig.fXminT0),
140     fXmaxT0(sig.fXmaxT0),
141     fNbinsQ(sig.fNbinsQ),
142     fXminQ(sig.fXminQ),
143     fXmaxQ(sig.fXmaxQ),
144     fNbinsRMS(sig.fNbinsRMS),
145     fXminRMS(sig.fXminRMS),
146     fXmaxRMS(sig.fXmaxRMS),
147     fLastSector(-1),
148     fOldRCUformat(kTRUE),
149     fROC(AliTPCROC::Instance()),
150     fParam(new AliTPCParam),
151     fPedestalTPC(0x0),
152     fPadNoiseTPC(0x0),
153     fPedestalROC(0x0),
154     fPadNoiseROC(0x0),
155     fCalRocArrayT0(72),
156     fCalRocArrayQ(72),
157     fCalRocArrayRMS(72),
158     fCalRocArrayOutliers(72),
159     fHistoQArray(72),
160     fHistoT0Array(72),
161     fHistoRMSArray(72),
162     fHistoTmean(72),
163     fParamArrayEvent(1000),
164     fParamArrayEventPol1(72),
165     fParamArrayEventPol2(72),
166     fTMeanArrayEvent(1000),
167     fQMeanArrayEvent(1000),
168     fVEventTime(1000),
169     fVEventNumber(1000),
170     fNevents(sig.fNevents),
171     fTimeStamp(0),
172     fRunNumber(-1),
173     fOldRunNumber(-1),
174     fPadTimesArrayEvent(72),
175     fPadQArrayEvent(72),
176     fPadRMSArrayEvent(72),
177     fPadPedestalArrayEvent(72),
178     fCurrentChannel(-1),
179     fCurrentSector(-1),
180     fCurrentRow(-1),
181     fMaxPadSignal(-1),
182     fMaxTimeBin(-1),
183     fPadSignal(1024),
184     fPadPedestal(0),
185     fPadNoise(0),
186     fVTime0Offset(72),
187     fVTime0OffsetCounter(72),
188     fVMeanQ(72),
189     fVMeanQCounter(72),
190 //    fHTime0(0x0),
191     fEvent(-1),
192     fDebugStreamer(0x0),
193     fDebugLevel(sig.fDebugLevel)
194 {
195     //
196     // AliTPCSignal default constructor
197     //
198
199     for (Int_t iSec = 0; iSec < 72; iSec++){
200         const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
201         const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
202         const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
203         const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
204
205         const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
206         const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
207         const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
208
209         if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
210         if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
211         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
212         if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
213
214         if ( hQ != 0x0 ){
215             TH2S *hNew = new TH2S(*hQ);
216             hNew->SetDirectory(0);
217             fHistoQArray.AddAt(hNew,iSec);
218         }
219         if ( hT0 != 0x0 ){
220             TH2S *hNew = new TH2S(*hT0);
221             hNew->SetDirectory(0);
222             fHistoQArray.AddAt(hNew,iSec);
223         }
224         if ( hRMS != 0x0 ){
225             TH2S *hNew = new TH2S(*hRMS);
226             hNew->SetDirectory(0);
227             fHistoQArray.AddAt(hNew,iSec);
228         }
229     }
230     for (Int_t iEvent=0; iEvent<sig.fParamArrayEvent.GetSize(); iEvent++)
231         fParamArrayEvent.AddAtAndExpand(sig.fParamArrayEvent.At(iEvent),iEvent);
232     Int_t nrows = sig.fVEventTime.GetNrows();
233     fVEventTime.ResizeTo(nrows);
234     for (Int_t iEvent=0; iEvent<nrows; iEvent++)
235         fVEventTime[iEvent] = sig.fVEventTime[iEvent];
236
237 //    fHTime0 = new TH1F("hTime0Event","hTime0Event",(fLastTimeBin-fFirstTimeBin)*10,fFirstTimeBin,fLastTimeBin);
238 }
239 //_____________________________________________________________________
240 AliTPCCalibCE& AliTPCCalibCE::operator = (const  AliTPCCalibCE &source)
241 {
242   //
243   // assignment operator
244   //
245   if (&source == this) return *this;
246   new (this) AliTPCCalibCE(source);
247
248   return *this;
249 }
250 //_____________________________________________________________________
251 AliTPCCalibCE::~AliTPCCalibCE()
252 {
253     //
254     // destructor
255     //
256
257     fCalRocArrayT0.Delete();
258     fCalRocArrayQ.Delete();
259     fCalRocArrayRMS.Delete();
260
261     fHistoQArray.Delete();
262     fHistoT0Array.Delete();
263     fHistoRMSArray.Delete();
264
265     fPadTimesArrayEvent.Delete();
266     fPadQArrayEvent.Delete();
267     fPadRMSArrayEvent.Delete();
268     fPadPedestalArrayEvent.Delete();
269
270     if ( fDebugStreamer) delete fDebugStreamer;
271 //    if ( fHTime0 ) delete fHTime0;
272     delete fROC;
273     delete fParam;
274 }
275 //_____________________________________________________________________
276 Int_t AliTPCCalibCE::Update(const Int_t icsector, /*FOLD00*/
277                                 const Int_t icRow,
278                                 const Int_t icPad,
279                                 const Int_t icTimeBin,
280                                 const Float_t csignal)
281 {
282     //
283     // Signal filling methode on the fly pedestal and Time offset correction if necessary.
284     // no extra analysis necessary. Assumes knowledge of the signal shape!
285     // assumes that it is looped over consecutive time bins of one pad
286     //
287     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
288
289     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
290
291     //init first pad and sector in this event
292     if ( fCurrentChannel == -1 ) {
293         fCurrentChannel = iChannel;
294         fCurrentSector  = icsector;
295         fCurrentRow     = icRow;
296     }
297
298     //process last pad if we change to a new one
299     if ( iChannel != fCurrentChannel ){
300         ProcessPad();
301         fCurrentChannel = iChannel;
302         fCurrentSector  = icsector;
303         fCurrentRow     = icRow;
304     }
305
306     //fill signals for current pad
307     fPadSignal.GetMatrixArray()[icTimeBin]=csignal;
308     if ( csignal > fMaxPadSignal ){
309         fMaxPadSignal = csignal;
310         fMaxTimeBin   = icTimeBin;
311     }
312     return 0;
313 }
314 //_____________________________________________________________________
315 void AliTPCCalibCE::FindPedestal(Float_t part)
316 {
317     //
318     // find pedestal and noise for the current pad. Use either database or
319     // truncated mean with part*100%
320     //
321     Bool_t noPedestal = kTRUE;;
322     if (fPedestalTPC&&fPadNoiseTPC){
323         //use pedestal database
324         //only load new pedestals if the sector has changed
325         if ( fCurrentSector!=fLastSector ){
326             fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
327             fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
328             fLastSector=fCurrentSector;
329         }
330
331         if ( fPedestalROC&&fPadNoiseROC ){
332             fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
333             fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
334             noPedestal   = kFALSE;
335         }
336
337     }
338
339     //if we are not running with pedestal database, or for the current sector there is no information
340     //available, calculate the pedestal and noise on the fly
341     if ( noPedestal ) {
342         const Int_t kPedMax = 100;  //maximum pedestal value
343         Float_t  max    =  0;
344         Float_t  maxPos =  0;
345         Int_t    median =  -1;
346         Int_t    count0 =  0;
347         Int_t    count1 =  0;
348         //
349         Float_t padSignal=0;
350         //
351         UShort_t histo[kPedMax];
352         memset(histo,0,kPedMax*sizeof(UShort_t));
353
354         for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; i++){
355             padSignal = fPadSignal.GetMatrixArray()[i];
356             if (padSignal<=0) continue;
357             if (padSignal>max && i>10) {
358                 max = padSignal;
359                 maxPos = i;
360             }
361             if (padSignal>kPedMax-1) continue;
362             histo[int(padSignal+0.5)]++;
363             count0++;
364         }
365             //
366         for (Int_t i=1; i<kPedMax; i++){
367             if (count1<count0*0.5) median=i;
368             count1+=histo[i];
369         }
370         // truncated mean
371         //
372         Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
373         //
374         for (Int_t idelta=1; idelta<10; idelta++){
375             if (median-idelta<=0) continue;
376             if (median+idelta>kPedMax) continue;
377             if (count<part*count1){
378                 count+=histo[median-idelta];
379                 mean +=histo[median-idelta]*(median-idelta);
380                 rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
381                 count+=histo[median+idelta];
382                 mean +=histo[median+idelta]*(median+idelta);
383                 rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
384             }
385         }
386         fPadPedestal = 0;
387         fPadNoise    = 0;
388         if ( count > 0 ) {
389             mean/=count;
390             rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
391             fPadPedestal = mean;
392             fPadNoise    = rms;
393         } 
394     }
395 }
396 //_____________________________________________________________________
397 void AliTPCCalibCE::FindCESignal(TVectorD &param, Float_t &qSum, const TVectorF maxima)
398 {
399     //
400     //  Find position, signal width and height of the CE signal (last signal)
401     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
402     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
403     //
404
405     Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
406     Int_t   cemaxpos       = 0;
407     Float_t ceSumThreshold = 8.*fPadNoise;  // threshold for the signal sum
408     const Int_t    kCemin  = 4;             // range for the analysis of the ce signal +- channels from the peak
409     const Int_t    kCemax  = 7;
410
411     Float_t minDist  = 25;  //initial minimum distance betweek roc mean ce signal and pad ce signal
412
413     // find maximum closest to the sector mean from the last event
414     for ( Int_t imax=0; imax<maxima.GetNrows(); imax++){
415         Float_t tmean = (*((TVectorF*)(fTMeanArrayEvent[fTMeanArrayEvent.GetLast()])))[fCurrentSector];
416             if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
417                 minDist  = tmean-maxima[imax];
418                 cemaxpos = (Int_t)maxima[imax];
419             }
420     }
421
422     if (cemaxpos!=0){
423         ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
424         for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; i++){
425             Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
426             if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
427                 ceTime+=signal*(i+0.5);
428                 ceRMS +=signal*(i+0.5)*(i+0.5);
429                 ceQsum+=signal;
430             }
431         }
432     }
433     if (ceQmax&&ceQsum>ceSumThreshold) {
434         ceTime/=ceQsum;
435         ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
436         fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
437         fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
438
439         //Normalise Q to pad area of irocs
440         Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
441
442         ceQsum/=norm;
443         fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
444         fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
445     } else {
446         ceQmax=0;
447         ceTime=0;
448         ceRMS =0;
449         ceQsum=0;
450     }
451     param[0] = ceQmax;
452     param[1] = ceTime;
453     param[2] = ceRMS;
454     qSum     = ceQsum;
455 }
456 //_____________________________________________________________________
457 Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus)
458 {
459     //
460     // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
461     // and 'tplus' timebins after 'pos'
462     //
463     for (Int_t iTime = pos; iTime>pos-tminus; iTime--)
464         if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
465     for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; iTime++,iTime2++){
466         if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
467             iTime2++;
468         if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
469     }
470     return kTRUE;
471 }
472 //_____________________________________________________________________
473 void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
474 {
475     //
476     // Find local maxima on the pad signal and Histogram them
477     //
478     Float_t ceThreshold = 5.*fPadNoise;  // threshold for the signal
479     Int_t   count       = 0;
480     Int_t   tminus      = 2;
481     Int_t   tplus       = 3;
482     for (Int_t i=fLastTimeBin-tplus-1; i>=fFirstTimeBin+tminus; i--){
483         if ( (fPadSignal[i]-fPadPedestal)>ceThreshold && IsPeak(i,tminus,tplus) ){
484             maxima.GetMatrixArray()[count++]=i;
485             GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
486         }
487     }
488 }
489 //_____________________________________________________________________
490 void AliTPCCalibCE::ProcessPad() /*FOLD00*/
491 {
492     //
493     //  Process data of current pad
494     //
495     FindPedestal();
496
497     TVectorF maxima(10);
498     FindLocalMaxima(maxima);
499     if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return;  // return because we don't have Time0 info for the CE yet
500
501
502
503     TVectorD param(3);
504     Float_t  Qsum;
505     FindCESignal(param, Qsum, maxima);
506
507     Double_t meanT  = param[1];
508     Double_t sigmaT = param[2];
509
510     //Fill Event T0 counter
511     (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
512
513     //Fill Q histogram
514     GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(Qsum), fCurrentChannel );
515
516     //Fill RMS histogram
517     GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
518
519
520     //Fill debugging info
521     if ( fDebugLevel>0 ){
522         (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
523         (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
524         (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=Qsum;
525     }
526
527     ResetPad();
528 }
529 //_____________________________________________________________________
530 void AliTPCCalibCE::EndEvent() /*FOLD00*/
531 {
532     //
533     //  Process data of current pad
534     //  The Functions 'SetTimeStamp' and 'SetRunNumber'  should be called
535     //  before the EndEvent function to set the event timestamp and number!!!
536     //  This is automatically done if the ProcessEvent(AliRawReader *rawReader)
537     //  function was called
538     //
539
540     //check if last pad has allready been processed, if not do so
541     if ( fMaxTimeBin>-1 ) ProcessPad();
542
543     TVectorD param(3);
544     TMatrixD dummy(3,3);
545     TVectorF vMeanTime(72);
546     TVectorF vMeanQ(72);
547     AliTPCCalROC calIroc(0);
548     AliTPCCalROC calOroc(36);
549
550     //find mean time0 offset for side A and C
551     Double_t time0Side[2];       //time0 for side A:0 and C:0
552     Double_t time0SideCount[2];  //time0 counter for side A:0 and C:0
553     time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
554     for ( Int_t iSec = 0; iSec<72; iSec++ ){
555         time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
556         time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
557     }
558     if ( time0SideCount[0] >0  )
559         time0Side[0]/=time0SideCount[0];
560     if ( time0SideCount[1] >0 )
561         time0Side[1]/=time0SideCount[1];
562     // end find time0 offset
563
564     //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
565     for ( Int_t iSec = 0; iSec<72; iSec++ ){
566
567       //find median and then calculate the mean around it
568         TH1S *hMeanT    = GetHistoTmean(iSec);
569         if ( !hMeanT ) continue;
570         Double_t entries = hMeanT->GetEntries();
571         Double_t sum     = 0;
572         Short_t *arr     = hMeanT->GetArray()+1;
573         Int_t ibin=0;
574         for ( ibin=0; ibin<hMeanT->GetNbinsX(); ibin++){
575             sum+=arr[ibin];
576             if ( sum>=(entries/2) ) break;
577         }
578         Int_t delta = 4;
579         Int_t firstBin = fFirstTimeBin+ibin-delta;
580         Int_t lastBin  = fFirstTimeBin+ibin+delta;
581         if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
582         if ( lastBin>fLastTimeBin   ) lastBin =fLastTimeBin;
583         Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
584         vMeanTime.GetMatrixArray()[iSec]=median;
585       // end find median
586
587         TVectorF *vTimes = GetPadTimesEvent(iSec);
588         if ( !vTimes ) continue;
589         AliTPCCalROC calIrocOutliers(0);
590         AliTPCCalROC calOrocOutliers(36);
591
592         // calculate mean Q of the sector
593         Float_t meanQ = 0;
594         if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
595         vMeanQ.GetMatrixArray()[iSec]=meanQ;
596
597         for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); iChannel++ ){
598             Float_t Time  = (*vTimes).GetMatrixArray()[iChannel];
599
600             //set values for temporary roc calibration class
601             if ( iSec < 36 ) {
602                 calIroc.SetValue(iChannel, Time);
603                 if ( Time == 0 ) calIrocOutliers.SetValue(iChannel,1);
604
605             } else {
606                 calOroc.SetValue(iChannel, Time);
607                 if ( Time == 0 ) calOrocOutliers.SetValue(iChannel,1);
608             }
609
610             if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
611                 GetHistoT0(iSec,kTRUE)->Fill( Time-time0Side[(iSec/18)%2],iChannel );
612
613
614
615             //-------------------------------  Debug start  ------------------------------
616             if ( fDebugLevel>0 ){
617                 if ( !fDebugStreamer ) {
618                         //debug stream
619                     TDirectory *backup = gDirectory;
620                     fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
621                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
622                 }
623
624                 Int_t row=0;
625                 Int_t pad=0;
626                 Int_t padc=0;
627
628                 Float_t Q   = (*GetPadQEvent(iSec))[iChannel];
629                 Float_t RMS = (*GetPadRMSEvent(iSec))[iChannel];
630
631                 UInt_t channel=iChannel;
632                 Int_t sector=iSec;
633
634                 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
635                 pad = channel-fROC->GetRowIndexes(sector)[row];
636                 padc = pad-(fROC->GetNPads(sector,row)/2);
637
638 //              TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
639 //                                  Form("hSignalD%d.%d.%d",sector,row,pad),
640 //                                  fLastTimeBin-fFirstTimeBin,
641 //                                  fFirstTimeBin,fLastTimeBin);
642 //              h1->SetDirectory(0);
643 //
644 //              for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
645 //                  h1->Fill(i,fPadSignal(i));
646
647                 Double_t T0Sec = 0;
648                 if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
649                     T0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
650                 Double_t T0Side = time0Side[(iSec/18)%2];
651                 (*fDebugStreamer) << "DataPad" <<
652                     "Event=" << fNevents <<
653                     "EventNumber=" << fRunNumber <<
654                     "TimeStamp="   << fTimeStamp <<
655                     "Sector="<< sector <<
656                     "Row="   << row<<
657                     "Pad="   << pad <<
658                     "PadC="  << padc <<
659                     "PadSec="<< channel <<
660                     "Time0Sec="  << T0Sec <<
661                     "Time0Side=" << T0Side <<
662                     "Time="  << Time <<
663                     "RMS="   << RMS <<
664                     "Sum="   << Q <<
665                     "MeanQ=" << meanQ <<
666                     //              "hist.=" << h1 <<
667                     "\n";
668
669                 //              delete h1;
670
671             }
672             //-----------------------------  Debug end  ------------------------------
673         }// end channel loop
674         hMeanT->Reset();
675
676         TVectorD paramPol1(3);
677         TVectorD paramPol2(6);
678         TMatrixD matPol1(3,3);
679         TMatrixD matPol2(6,6);
680         Float_t  chi2Pol1=0;
681         Float_t  chi2Pol2=0;
682
683         if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
684             if ( iSec < 36 ){
685                 calIroc.GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
686                 calIroc.GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
687             } else {
688                 calOroc.GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
689                 calOroc.GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
690             }
691
692             GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
693             GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
694         }
695 //      printf("events: %d -- size: %d\n",fNevents,GetParamArrayPol1(iSec)->GetSize());
696
697         //-------------------------------  Debug start  ------------------------------
698         if ( fDebugLevel>0 ){
699             if ( !fDebugStreamer ) {
700                 //debug stream
701                 TDirectory *backup = gDirectory;
702                 fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
703                 if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
704             }
705             (*fDebugStreamer) << "DataRoc" <<
706                 "Event=" << fEvent <<
707                 "EventNumber=" << fRunNumber <<
708                 "TimeStamp="   << fTimeStamp <<
709                 "Sector="<< iSec <<
710                 "hMeanT.=" << hMeanT <<
711                 "median=" << median <<
712                 "paramPol1.=" << &paramPol1 <<
713                 "paramPol2.=" << &paramPol2 <<
714                 "matPol1.="   << &matPol1 <<
715                 "matPol2.="   << &matPol2 <<
716                 "chi2Pol1="   << chi2Pol1 <<
717                 "chi2Pol2="   << chi2Pol2 <<
718                 "\n";
719         }
720         //-------------------------------  Debug end  ------------------------------
721     }// end sector loop
722
723 /*    AliMathBase::FitGaus(fHTime0->GetArray()+1,
724                          fHTime0->GetNbinsX(),
725                          fHTime0->GetXaxis()->GetXmin(),
726                          fHTime0->GetXaxis()->GetXmax(),
727                          &param, &dummy);*/
728 //    fHTime0->Reset();
729
730     //    fParamArrayEvent.AddAtAndExpand(new TVectorD(param),fNevents);
731     fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
732     fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
733     if ( fVEventTime.GetNrows() < fNevents ) {
734         fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+1000));
735         fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+1000));
736     }
737     fVEventTime[fNevents] = fTimeStamp;
738     fVEventNumber[fNevents] = fRunNumber;
739
740     fNevents++;
741     fOldRunNumber = fRunNumber;
742
743 }
744 //_____________________________________________________________________
745 Bool_t AliTPCCalibCE::ProcessEvent(AliTPCRawStream *rawStream) /*FOLD00*/
746 {
747   //
748   // Event Processing loop - AliTPCRawStream
749   // The Function 'SetTimeStamp' should be called for each event to set the event time stamp!!!
750   //
751
752   rawStream->SetOldRCUFormat(fOldRCUformat);
753
754   ResetEvent();
755
756   Bool_t withInput = kFALSE;
757
758   while (rawStream->Next()) {
759
760       Int_t isector  = rawStream->GetSector();                       //  current sector
761       Int_t iRow     = rawStream->GetRow();                          //  current row
762       Int_t iPad     = rawStream->GetPad();                          //  current pad
763       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
764       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
765
766       Update(isector,iRow,iPad,iTimeBin,signal);
767       withInput = kTRUE;
768   }
769
770   if (withInput){
771       EndEvent();
772   }
773
774   return withInput;
775 }
776 //_____________________________________________________________________
777 Bool_t AliTPCCalibCE::ProcessEvent(AliRawReader *rawReader)
778 {
779   //
780   //  Event processing loop - AliRawReader
781   //
782
783
784     AliTPCRawStream rawStream(rawReader);
785     AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
786     if (eventHeader){
787         fTimeStamp   = eventHeader->Get("Timestamp");
788         fRunNumber = eventHeader->Get("RunNb");
789     }
790
791
792   rawReader->Select("TPC");
793
794   return ProcessEvent(&rawStream);
795 }
796 //_____________________________________________________________________
797 Bool_t AliTPCCalibCE::ProcessEvent(eventHeaderStruct *event)
798 {
799   //
800   //  Event processing loop - date event
801   //
802     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
803     Bool_t result=ProcessEvent(rawReader);
804     delete rawReader;
805     return result;
806
807 }
808 //_____________________________________________________________________
809 TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
810                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
811                                   Char_t *type, Bool_t force)
812 {
813     //
814     // return pointer to TH2S histogram of 'type'
815     // if force is true create a new histogram if it doesn't exist allready
816     //
817     if ( !force || arr->UncheckedAt(sector) )
818         return (TH2S*)arr->UncheckedAt(sector);
819
820     // if we are forced and histogram doesn't yes exist create it
821     Char_t name[255], title[255];
822
823     sprintf(name,"hCalib%s%.2d",type,sector);
824     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
825
826     // new histogram with Q calib information. One value for each pad!
827     TH2S* hist = new TH2S(name,title,
828                           nbinsY, ymin, ymax,
829                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
830     hist->SetDirectory(0);
831     arr->AddAt(hist,sector);
832     return hist;
833 }
834 //_____________________________________________________________________
835 TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force) /*FOLD00*/
836 {
837     //
838     // return pointer to T0 histogram
839     // if force is true create a new histogram if it doesn't exist allready
840     //
841     TObjArray *arr = &fHistoT0Array;
842     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
843 }
844 //_____________________________________________________________________
845 TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force) /*FOLD00*/
846 {
847     //
848     // return pointer to Q histogram
849     // if force is true create a new histogram if it doesn't exist allready
850     //
851     TObjArray *arr = &fHistoQArray;
852     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
853 }
854 //_____________________________________________________________________
855 TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force) /*FOLD00*/
856 {
857     //
858     // return pointer to Q histogram
859     // if force is true create a new histogram if it doesn't exist allready
860     //
861     TObjArray *arr = &fHistoRMSArray;
862     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
863 }
864 //_____________________________________________________________________
865 TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
866                               Char_t *type, Bool_t force)
867 {
868     //
869     // return pointer to TH1S histogram
870     // if force is true create a new histogram if it doesn't exist allready
871     //
872     if ( !force || arr->UncheckedAt(sector) )
873         return (TH1S*)arr->UncheckedAt(sector);
874
875     // if we are forced and histogram doesn't yes exist create it
876     Char_t name[255], title[255];
877
878     sprintf(name,"hCalib%s%.2d",type,sector);
879     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
880
881     // new histogram with Q calib information. One value for each pad!
882     TH1S* hist = new TH1S(name,title,
883                           fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
884     hist->SetDirectory(0);
885     arr->AddAt(hist,sector);
886     return hist;
887 }
888 //_____________________________________________________________________
889 TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
890 {
891     //
892     // return pointer to Q histogram
893     // if force is true create a new histogram if it doesn't exist allready
894     //
895     TObjArray *arr = &fHistoTmean;
896     return GetHisto(sector, arr, "LastTmean", force);
897 }
898 //_____________________________________________________________________
899 TVectorF* AliTPCCalibCE::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force) /*FOLD00*/
900 {
901     //
902     // return pointer to Pad Info from 'arr' for the current event and sector
903     // if force is true create it if it doesn't exist allready
904     //
905     if ( !force || arr->UncheckedAt(sector) )
906         return (TVectorF*)arr->UncheckedAt(sector);
907
908     TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
909     arr->AddAt(vect,sector);
910     return vect;
911 }
912 //_____________________________________________________________________
913 TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force) /*FOLD00*/
914 {
915     //
916     // return pointer to Pad Times Array for the current event and sector
917     // if force is true create it if it doesn't exist allready
918     //
919     TObjArray *arr = &fPadTimesArrayEvent;
920     return GetPadInfoEvent(sector,arr,force);
921 }
922 //_____________________________________________________________________
923 TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force) /*FOLD00*/
924 {
925     //
926     // return pointer to Pad Q Array for the current event and sector
927     // if force is true create it if it doesn't exist allready
928     // for debugging purposes only
929     //
930
931     TObjArray *arr = &fPadQArrayEvent;
932     return GetPadInfoEvent(sector,arr,force);
933 }
934 //_____________________________________________________________________
935 TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force) /*FOLD00*/
936 {
937     //
938     // return pointer to Pad RMS Array for the current event and sector
939     // if force is true create it if it doesn't exist allready
940     // for debugging purposes only
941     //
942     TObjArray *arr = &fPadRMSArrayEvent;
943     return GetPadInfoEvent(sector,arr,force);
944 }
945 //_____________________________________________________________________
946 TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force) /*FOLD00*/
947 {
948     //
949     // return pointer to Pad RMS Array for the current event and sector
950     // if force is true create it if it doesn't exist allready
951     // for debugging purposes only
952     //
953     TObjArray *arr = &fPadPedestalArrayEvent;
954     return GetPadInfoEvent(sector,arr,force);
955 }
956 //_____________________________________________________________________
957 AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
958 {
959     //
960     // return pointer to ROC Calibration
961     // if force is true create a new histogram if it doesn't exist allready
962     //
963     if ( !force || arr->UncheckedAt(sector) )
964         return (AliTPCCalROC*)arr->UncheckedAt(sector);
965
966     // if we are forced and histogram doesn't yes exist create it
967
968     // new AliTPCCalROC for T0 information. One value for each pad!
969     AliTPCCalROC *croc = new AliTPCCalROC(sector);
970     //init values
971     for ( UInt_t iChannel = 0; iChannel<croc->GetNchannels(); iChannel++){
972         croc->SetValue(iChannel, 0);
973     }
974     arr->AddAt(croc,sector);
975     return croc;
976 }
977 //_____________________________________________________________________
978 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force) /*FOLD00*/
979 {
980     //
981     // return pointer to Carge ROC Calibration
982     // if force is true create a new histogram if it doesn't exist allready
983     //
984     TObjArray *arr = &fCalRocArrayT0;
985     return GetCalRoc(sector, arr, force);
986 }
987 //_____________________________________________________________________
988 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force) /*FOLD00*/
989 {
990     //
991     // return pointer to T0 ROC Calibration
992     // if force is true create a new histogram if it doesn't exist allready
993     //
994     TObjArray *arr = &fCalRocArrayQ;
995     return GetCalRoc(sector, arr, force);
996 }
997 //_____________________________________________________________________
998 AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
999 {
1000     //
1001     // return pointer to signal width ROC Calibration
1002     // if force is true create a new histogram if it doesn't exist allready
1003     //
1004     TObjArray *arr = &fCalRocArrayRMS;
1005     return GetCalRoc(sector, arr, force);
1006 }
1007 //_____________________________________________________________________
1008 AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1009 {
1010     //
1011     // return pointer to Outliers
1012     // if force is true create a new histogram if it doesn't exist allready
1013     //
1014     TObjArray *arr = &fCalRocArrayOutliers;
1015     return GetCalRoc(sector, arr, force);
1016 }
1017 //_____________________________________________________________________
1018 TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force)
1019 {
1020     //
1021     // return pointer to TObjArray of fit parameters
1022     // if force is true create a new histogram if it doesn't exist allready
1023     //
1024     if ( !force || arr->UncheckedAt(sector) )
1025         return (TObjArray*)arr->UncheckedAt(sector);
1026
1027     // if we are forced and array doesn't yes exist create it
1028
1029     // new TObjArray for parameters
1030     TObjArray *newArr = new TObjArray;
1031     arr->AddAt(newArr,sector);
1032     return newArr;
1033 }
1034 //_____________________________________________________________________
1035 TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1036 {
1037     //
1038     // return pointer to TObjArray of fit parameters from plane fit
1039     // if force is true create a new histogram if it doesn't exist allready
1040     //
1041     TObjArray *arr = &fParamArrayEventPol1;
1042     return GetParamArray(sector, arr, force);
1043 }
1044 //_____________________________________________________________________
1045 TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1046 {
1047     //
1048     // return pointer to TObjArray of fit parameters from parabola fit
1049     // if force is true create a new histogram if it doesn't exist allready
1050     //
1051     TObjArray *arr = &fParamArrayEventPol2;
1052     return GetParamArray(sector, arr, force);
1053 }
1054 //_____________________________________________________________________
1055 void AliTPCCalibCE::ResetEvent() /*FOLD00*/
1056 {
1057     //
1058     //  Reset global counters  -- Should be called before each event is processed
1059     //
1060     fLastSector=-1;
1061     fCurrentSector=-1;
1062     fCurrentRow=-1;
1063     fCurrentChannel=-1;
1064
1065     ResetPad();
1066
1067     fPadTimesArrayEvent.Delete();
1068     fPadQArrayEvent.Delete();
1069     fPadRMSArrayEvent.Delete();
1070     fPadPedestalArrayEvent.Delete();
1071
1072     for ( Int_t i=0; i<72; i++ ){
1073         fVTime0Offset.GetMatrixArray()[i]=0;
1074         fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1075         fVMeanQ.GetMatrixArray()[i]=0;
1076         fVMeanQCounter.GetMatrixArray()[i]=0;
1077     }
1078 }
1079 //_____________________________________________________________________
1080 void AliTPCCalibCE::ResetPad() /*FOLD00*/
1081 {
1082     //
1083     //  Reset pad infos -- Should be called after a pad has been processed
1084     //
1085     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
1086         fPadSignal.GetMatrixArray()[i] = 0;
1087     fMaxTimeBin   = -1;
1088     fMaxPadSignal = -1;
1089     fPadPedestal  = -1;
1090     fPadNoise     = -1;
1091 }
1092 //_____________________________________________________________________
1093 TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter) /*FOLD00*/
1094 {
1095     //
1096     // Make graph from fit parameters of pol1 or pol2 fit
1097     // xVariable:    0-run time, 1-run number, 2-internal event counter
1098     // fitType:      0-pol1 fit, 1-pol2 fit, 2-mean time, 2-mean Q
1099     // fitParameter: fit parameter ( 0-2 for pol1, 0-5 for pol2, 0 for mean time )
1100     //
1101
1102     Double_t *x = new Double_t[fNevents];
1103     Double_t *y = new Double_t[fNevents];
1104
1105     TVectorD *xVar = 0x0;
1106     TObjArray *aType = 0x0;
1107     Int_t npoints=0;
1108
1109     // sanity checks
1110     if ( (sector<0) || (sector>71) )      return 0x0;
1111     if ( (xVariable<0) || (xVariable>2) ) return 0x0;
1112     if ( (fitType<0) || (fitType>3) )     return 0x0;
1113     if ( fitType==0 ){
1114         if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1115         aType = &fParamArrayEventPol1;
1116         if ( aType->At(sector)==0x0 ) return 0x0;
1117     }
1118     else if ( fitType==1 ){
1119         if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1120         aType = &fParamArrayEventPol2;
1121         if ( aType->At(sector)==0x0 ) return 0x0;
1122     }
1123
1124
1125     if ( xVariable == 0 ) xVar = &fVEventTime;
1126     if ( xVariable == 1 ) xVar = &fVEventNumber;
1127     if ( xVariable == 2 ) {
1128         xVar = new TVectorD(fNevents);
1129         for ( Int_t i=0;i<fNevents; i++) (*xVar)[i]=i;
1130     }
1131
1132     for (Int_t ievent =0; ievent<fNevents; ievent++){
1133         if ( fitType<2 ){
1134             TObjArray *events = (TObjArray*)(aType->At(sector));
1135             if ( events->GetSize()<=ievent ) break;
1136             TVectorD *v = (TVectorD*)(events->At(ievent));
1137             if ( v!=0x0 ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1138         } else if (fitType == 2) {
1139             Double_t yValue=(*((TVectorF*)(fTMeanArrayEvent[ievent])))[sector];
1140             if ( yValue>0 ) { x[npoints]=(*xVar)[ievent]; y[npoints]=yValue;npoints++;}
1141         }else if (fitType == 3) {
1142             Double_t yValue=(*((TVectorF*)(fQMeanArrayEvent[ievent])))[sector];
1143             if ( yValue>0 ) { x[npoints]=(*xVar)[ievent]; y[npoints]=yValue;npoints++;}
1144         }
1145     }
1146
1147     TGraph *gr = new TGraph(npoints);
1148     //sort xVariable increasing
1149     Int_t    *sortIndex = new Int_t[npoints];
1150     TMath::Sort(npoints,x,sortIndex);
1151     for (Int_t i=0;i<npoints;i++){
1152         gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1153     }
1154
1155
1156     if ( xVariable == 2 ) delete xVar;
1157     delete x;
1158     delete y;
1159     delete sortIndex;
1160     return gr;
1161 }
1162 //_____________________________________________________________________
1163 void AliTPCCalibCE::Analyse()
1164 {
1165     //
1166     //  Calculate calibration constants
1167     //
1168
1169     TVectorD paramQ(3);
1170     TVectorD paramT0(3);
1171     TVectorD paramRMS(3);
1172     TMatrixD dummy(3,3);
1173
1174     for (Int_t iSec=0; iSec<72; iSec++){
1175         TH2S *hT0 = GetHistoT0(iSec);
1176         if (!hT0 ) continue;
1177
1178         AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
1179         AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
1180         AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1181         AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1182
1183         TH2S *hQ   = GetHistoQ(iSec);
1184         TH2S *hRMS = GetHistoRMS(iSec);
1185
1186         Short_t *array_hQ   = hQ->GetArray();
1187         Short_t *array_hT0  = hT0->GetArray();
1188         Short_t *array_hRMS = hRMS->GetArray();
1189
1190         UInt_t nChannels = fROC->GetNChannels(iSec);
1191
1192         //debug
1193         Int_t row=0;
1194         Int_t pad=0;
1195         Int_t padc=0;
1196         //! debug
1197
1198         for (UInt_t iChannel=0; iChannel<nChannels; iChannel++){
1199
1200
1201             Float_t cogTime0 = -1000;
1202             Float_t cogQ     = -1000;
1203             Float_t cogRMS   = -1000;
1204             Float_t cogOut   = 0;
1205
1206
1207             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1208             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1209             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1210
1211 /*
1212             AliMathBase::FitGaus(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1213             AliMathBase::FitGaus(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1214             AliMathBase::FitGaus(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1215             cogQ     = paramQ[1];
1216             cogTime0 = paramT0[1];
1217             cogRMS   = paramRMS[1];
1218 */
1219             cogQ     = AliMathBase::GetCOG(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1220             cogTime0 = AliMathBase::GetCOG(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1221             cogRMS   = AliMathBase::GetCOG(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1222
1223
1224
1225             /*
1226             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1227                 cogOut = 1;
1228                 cogTime0 = 0;
1229                 cogQ     = 0;
1230                 cogRMS   = 0;
1231             }
1232 */
1233             rocQ->SetValue(iChannel, cogQ*cogQ);
1234             rocT0->SetValue(iChannel, cogTime0);
1235             rocRMS->SetValue(iChannel, cogRMS);
1236             rocOut->SetValue(iChannel, cogOut);
1237
1238
1239             //debug
1240             if ( fDebugLevel > 0 ){
1241                 if ( !fDebugStreamer ) {
1242                         //debug stream
1243                     TDirectory *backup = gDirectory;
1244                     fDebugStreamer = new TTreeSRedirector("debugCalibCEAnalysis.root");
1245                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1246                 }
1247
1248                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1249                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1250                 padc = pad-(fROC->GetNPads(iSec,row)/2);
1251
1252                 (*fDebugStreamer) << "DataEnd" <<
1253                     "Sector="  << iSec      <<
1254                     "Pad="     << pad       <<
1255                     "PadC="    << padc      <<
1256                     "Row="     << row       <<
1257                     "PadSec="  << iChannel   <<
1258                     "Q="       << cogQ      <<
1259                     "T0="      << cogTime0  <<
1260                     "RMS="     << cogRMS    <<
1261                     "\n";
1262             }
1263             //! debug
1264
1265         }
1266
1267     }
1268     fDebugStreamer->GetFile()->Write();
1269 //    delete fDebugStreamer;
1270 //    fDebugStreamer = 0x0;
1271 }
1272 //_____________________________________________________________________
1273 void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1274 {
1275     //
1276     //  Write class to file
1277     //
1278
1279     TString sDir(dir);
1280     TString option;
1281
1282     if ( append )
1283         option = "update";
1284     else
1285         option = "recreate";
1286
1287     TDirectory *backup = gDirectory;
1288     TFile f(filename,option.Data());
1289     f.cd();
1290     if ( !sDir.IsNull() ){
1291         f.mkdir(sDir.Data());
1292         f.cd(sDir);
1293     }
1294     this->Write();
1295     f.Close();
1296
1297     if ( backup ) backup->cd();
1298 }