more secure string operations
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibSignal.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 <TMath.h>
40
41 #include <TDirectory.h>
42 #include <TSystem.h>
43 #include <TFile.h>
44
45 //AliRoot includes
46 #include "AliRawReader.h"
47 #include "AliRawReaderRoot.h"
48 #include "AliRawReaderDate.h"
49 #include "AliTPCRawStream.h"
50 #include "AliTPCCalROC.h"
51 #include "AliTPCCalPad.h"
52 #include "AliTPCROC.h"
53 #include "AliTPCParam.h"
54 #include "AliTPCCalibSignal.h"
55 #include "AliTPCcalibDB.h"
56 #include "AliMathBase.h"
57 #include "TTreeStream.h"
58
59 //date
60 #include "event.h"
61 ClassImp(AliTPCCalibSignal) /*FOLD00*/
62
63 AliTPCCalibSignal::AliTPCCalibSignal() : /*FOLD00*/
64     TObject(),
65     fFirstTimeBin(60),
66     fLastTimeBin(120),
67     fFirstTimeBinT0(-15),
68     fLastTimeBinT0(15),
69     fNbinsT0(200),
70     fXminT0(-2),
71     fXmaxT0(2),
72     fNbinsQ(200),
73     fXminQ(14),
74     fXmaxQ(55),
75     fNbinsRMS(100),
76     fXminRMS(0),
77     fXmaxRMS(5),
78     fLastSector(-1),
79     fROC(AliTPCROC::Instance()),
80     fParam(new AliTPCParam),
81     fPedestalTPC(0x0),
82     fBpedestal(kFALSE),
83     fCalRocArrayT0(72),
84     fCalRocArrayQ(72),
85     fCalRocArrayRMS(72),
86     fCalRocArrayOutliers(72),
87     fHistoQArray(72),
88     fHistoT0Array(72),
89     fHistoRMSArray(72),
90     fPadTimesArrayEvent(72),
91     fPadQArrayEvent(72),
92     fPadRMSArrayEvent(72),
93     fPadPedestalArrayEvent(72),
94     fCurrentChannel(-1),
95     fCurrentSector(-1),
96     fCurrentRow(-1),
97     fMaxPadSignal(-1),
98     fMaxTimeBin(-1),
99     fPadSignal(1024),
100     fVTime0Offset1(72),
101     fVTime0Offset1Counter(72),
102     fEvent(-1),
103     fDebugStreamer(0x0),
104     fDebugLevel(0)
105 {
106     //
107     // AliTPCSignal default constructor
108     //
109
110 }
111 //_____________________________________________________________________
112 AliTPCCalibSignal::AliTPCCalibSignal(const AliTPCCalibSignal &sig) :
113     TObject(sig),
114     fFirstTimeBin(sig.fFirstTimeBin),
115     fLastTimeBin(sig.fLastTimeBin),
116     fFirstTimeBinT0(sig.fFirstTimeBinT0),
117     fLastTimeBinT0(sig.fLastTimeBinT0),
118     fNbinsT0(sig.fNbinsT0),
119     fXminT0(sig.fXminT0),
120     fXmaxT0(sig.fXmaxT0),
121     fNbinsQ(sig.fNbinsQ),
122     fXminQ(sig.fXminQ),
123     fXmaxQ(sig.fXmaxQ),
124     fNbinsRMS(sig.fNbinsRMS),
125     fXminRMS(sig.fXminRMS),
126     fXmaxRMS(sig.fXmaxRMS),
127     fLastSector(-1),
128     fROC(AliTPCROC::Instance()),
129     fParam(new AliTPCParam),
130     fPedestalTPC(sig.fPedestalTPC),
131     fBpedestal(sig.fBpedestal),
132     fCalRocArrayT0(72),
133     fCalRocArrayQ(72),
134     fCalRocArrayRMS(72),
135     fCalRocArrayOutliers(72),
136     fHistoQArray(72),
137     fHistoT0Array(72),
138     fHistoRMSArray(72),
139     fPadTimesArrayEvent(72),
140     fPadQArrayEvent(72),
141     fPadRMSArrayEvent(72),
142     fPadPedestalArrayEvent(72),
143     fCurrentChannel(-1),
144     fCurrentSector(-1),
145     fCurrentRow(-1),
146     fMaxPadSignal(-1),
147     fMaxTimeBin(-1),
148     fPadSignal(1024),
149     fVTime0Offset1(72),
150     fVTime0Offset1Counter(72),
151     fEvent(-1),
152     fDebugStreamer(0x0),
153     fDebugLevel(sig.fDebugLevel)
154 {
155     //
156     // AliTPCSignal default constructor
157     //
158
159     for (Int_t iSec = 0; iSec < 72; iSec++){
160         const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
161         const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
162         const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
163         const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
164
165         const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
166         const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
167         const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
168
169         if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
170         if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
171         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
172         if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
173
174         if ( hQ != 0x0 ){
175             TH2S *hNew = new TH2S(*hQ);
176             hNew->SetDirectory(0);
177             fHistoQArray.AddAt(hNew,iSec);
178         }
179         if ( hT0 != 0x0 ){
180             TH2S *hNew = new TH2S(*hT0);
181             hNew->SetDirectory(0);
182             fHistoQArray.AddAt(hNew,iSec);
183         }
184         if ( hRMS != 0x0 ){
185             TH2S *hNew = new TH2S(*hRMS);
186             hNew->SetDirectory(0);
187             fHistoQArray.AddAt(hNew,iSec);
188         }
189     }
190
191 }
192 //_____________________________________________________________________
193 AliTPCCalibSignal& AliTPCCalibSignal::operator = (const  AliTPCCalibSignal &source)
194 {
195   //
196   // assignment operator
197   //
198   if (&source == this) return *this;
199   new (this) AliTPCCalibSignal(source);
200
201   return *this;
202 }
203 //_____________________________________________________________________
204 AliTPCCalibSignal::~AliTPCCalibSignal()
205 {
206     //
207     // destructor
208     //
209
210     fCalRocArrayT0.Delete();
211     fCalRocArrayQ.Delete();
212     fCalRocArrayRMS.Delete();
213
214     fHistoQArray.Delete();
215     fHistoT0Array.Delete();
216     fHistoRMSArray.Delete();
217
218     fPadTimesArrayEvent.Delete();
219     fPadQArrayEvent.Delete();
220     fPadRMSArrayEvent.Delete();
221     fPadPedestalArrayEvent.Delete();
222
223     if ( fDebugStreamer) delete fDebugStreamer;
224     delete fROC;
225     delete fParam;
226 }
227 //_____________________________________________________________________
228 Int_t AliTPCCalibSignal::Update(const Int_t icsector, /*FOLD00*/
229                                 const Int_t icRow,
230                                 const Int_t icPad,
231                                 const Int_t icTimeBin,
232                                 const Float_t csignal)
233 {
234     //
235     // Signal filling methode on the fly pedestal and Time offset correction if necessary.
236     // no extra analysis necessary. Assumes knowledge of the signal shape!
237     // assumes that it is looped over consecutive time bins of one pad
238     //
239     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
240
241     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
242
243     //init first pad and sector in this event
244     if ( fCurrentChannel == -1 ) {
245         fCurrentChannel = iChannel;
246         fCurrentSector  = icsector;
247         fCurrentRow     = icRow;
248     }
249
250     //process last pad if we change to a new one
251     if ( iChannel != fCurrentChannel ){
252         ProcessPad();
253         fCurrentChannel = iChannel;
254         fCurrentSector  = icsector;
255         fCurrentRow     = icRow;
256     }
257
258     //fill signals for current pad
259     fPadSignal[icTimeBin]=csignal;
260     if ( csignal > fMaxPadSignal ){
261         fMaxPadSignal = csignal;
262         fMaxTimeBin   = icTimeBin;
263     }
264     return 0;
265 }
266 //_____________________________________________________________________
267 void AliTPCCalibSignal::ProcessPad() /*FOLD00*/
268 {
269     //
270     //  Process data of current pad
271     //
272
273     Float_t pedestal = 0;
274
275     if ( fBpedestal ){
276         //!!!!!!! does not work like this
277         //use pedestal database
278         AliTPCCalROC *pedestalROC = 0x0;
279
280         //only load new pedestals if the sector has changed
281         if ( fCurrentSector!=fLastSector ){
282             pedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
283             fLastSector=fCurrentSector;
284         }
285
286         pedestal = pedestalROC->GetValue(fCurrentChannel);
287
288     } else {
289
290         //find pedestal for pad on the fly
291         //using a few timebins before the signal
292         Int_t pminus1 = 10, pminus2=5;
293         Float_t sumN=0;
294
295         for (Int_t i=fMaxTimeBin-pminus1; i<fMaxTimeBin-pminus2+1; i++){
296             if ( i>fFirstTimeBin && i<fLastTimeBin ){
297                 pedestal+=fPadSignal[i];
298                 sumN+=1.;
299             }
300         }
301
302         if ( sumN>0 ) pedestal/=sumN;
303     }
304
305
306
307     //!!!! check borders
308     //find signal mean and sigma
309     Int_t tminus = 2, tplus=7;
310     Double_t meanT=0, sigmaT=0, Qsum=0;
311
312
313     for (Int_t i=fMaxTimeBin-tminus; i<fMaxTimeBin+tplus; i++){
314         if ( i>=fFirstTimeBin && i<=fLastTimeBin ){
315             Double_t val=fPadSignal[i]-pedestal;
316             meanT+=val*(i+.5);      //+.5: center of the timebin
317             sigmaT+=val*(i+.5)*(i+.5);
318             Qsum+=val;
319         }
320     }
321
322
323
324     //!!!! What to do if Qsum == 0???
325     //!!!! Should there be some threshold for max - pedestal and/or Qsum???
326     //!!!! What if Qsum < 0
327     //!!!! only fill time0 offset if Qsum > 0???
328     if ( Qsum > 0 ){
329         meanT/=Qsum;
330         sigmaT/=Qsum;
331         sigmaT = TMath::Sqrt(TMath::Abs(meanT*meanT - sigmaT));
332
333         //fill Time0 offset data for this event
334         fVTime0Offset1[fCurrentSector]+=meanT;
335         fVTime0Offset1Counter[fCurrentSector]++;
336     } else {
337         Qsum=0;
338         meanT  = fLastTimeBinT0+1;               //put to overflow bin
339         sigmaT = fLastTimeBinT0-fFirstTimeBinT0; //put to overflow bin
340     }
341
342     //Fill Event T0 counter
343     (*GetPadTimesEvent(fCurrentSector,kTRUE))[fCurrentChannel] = meanT;
344
345
346     //Normalise Q to pad area of irocs
347     Float_t norm = fParam->GetPadPitchWidth(0)*fParam->GetPadPitchLength(0,0)/(
348         fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow));
349
350     //Fill Q histogram
351     GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(Qsum*norm), fCurrentChannel );
352
353     //Fill RMS histogram
354     GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
355
356
357     //Fill debugging info
358     if ( fDebugLevel>0 ){
359         (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=pedestal;
360         (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
361         (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=Qsum;
362     }
363
364     ResetPad();
365 }
366 //_____________________________________________________________________
367 void AliTPCCalibSignal::EndEvent() /*FOLD00*/
368 {
369     //
370     //  Process data of current pad
371     //
372     //check if last pad has allready been processed, if not do so
373     if ( fMaxTimeBin>-1 ) ProcessPad();
374
375     //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
376     for ( Int_t iSec = 0; iSec<72; iSec++ ){
377         TVectorF *vTimes = GetPadTimesEvent(iSec);
378         if ( !vTimes ) continue;
379
380         for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); iChannel++ ){
381             Float_t Time0 = fVTime0Offset1[iSec]/fVTime0Offset1Counter[iSec];
382             Float_t Time  = (*vTimes)[iChannel];
383
384             GetHistoT0(iSec,kTRUE)->Fill( Time-Time0,iChannel );
385
386
387             //Debug start
388             if ( fDebugLevel>0 ){
389                 if ( !fDebugStreamer ) {
390                         //debug stream
391                     TDirectory *backup = gDirectory;
392                     fDebugStreamer = new TTreeSRedirector("deb2.root");
393                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
394                 }
395
396                 Int_t row=0;
397                 Int_t pad=0;
398                 Int_t padc=0;
399
400                 Float_t Q   = (*GetPadQEvent(iSec))[iChannel];
401                 Float_t RMS = (*GetPadRMSEvent(iSec))[iChannel];
402
403                 UInt_t channel=iChannel;
404                 Int_t sector=iSec;
405
406                 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
407                 pad = channel-fROC->GetRowIndexes(sector)[row];
408                 padc = pad-(fROC->GetNPads(sector,row)/2);
409
410                 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
411                                     Form("hSignalD%d.%d.%d",sector,row,pad),
412                                     fLastTimeBin-fFirstTimeBin,
413                                     fFirstTimeBin,fLastTimeBin);
414                 h1->SetDirectory(0);
415
416                 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
417                     h1->Fill(i,fPadSignal(i));
418
419                 (*fDebugStreamer) << "DataPad" <<
420                     "Event=" << fEvent <<
421                     "Sector="<< sector <<
422                     "Row="   << row<<
423                     "Pad="   << pad <<
424                     "PadC="  << padc <<
425                     "PadSec="<< channel <<
426                     "Time0="  << Time0 <<
427                     "Time="  << Time <<
428                     "RMS="   << RMS <<
429                     "Sum="   << Q <<
430                     "hist.=" << h1 <<
431                     "\n";
432
433                 delete h1;
434             }
435             //Debug end
436
437         }
438     }
439
440 }
441 //_____________________________________________________________________
442 Bool_t AliTPCCalibSignal::ProcessEvent(AliTPCRawStream *rawStream) /*FOLD00*/
443 {
444   //
445   // Event Processing loop - AliTPCRawStream
446   //
447
448   ResetEvent();
449
450   Bool_t withInput = kFALSE;
451
452   while (rawStream->Next()) {
453
454       Int_t isector  = rawStream->GetSector();                       //  current sector
455       Int_t iRow     = rawStream->GetRow();                          //  current row
456       Int_t iPad     = rawStream->GetPad();                          //  current pad
457       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
458       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
459
460       Update(isector,iRow,iPad,iTimeBin,signal);
461       withInput = kTRUE;
462   }
463
464   if (withInput){
465       EndEvent();
466   }
467
468   return withInput;
469 }
470 //_____________________________________________________________________
471 Bool_t AliTPCCalibSignal::ProcessEvent(AliRawReader *rawReader)
472 {
473   //
474   //  Event processing loop - AliRawReader
475   //
476
477
478   AliTPCRawStream rawStream(rawReader);
479
480   rawReader->Select("TPC");
481
482   return ProcessEvent(&rawStream);
483 }
484 //_____________________________________________________________________
485 Bool_t AliTPCCalibSignal::ProcessEvent(eventHeaderStruct *event)
486 {
487   //
488   //  Event processing loop - date event
489   //
490     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
491     Bool_t result=ProcessEvent(rawReader);
492     delete rawReader;
493     return result;
494
495 }
496 //_____________________________________________________________________
497 TH2S* AliTPCCalibSignal::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
498                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
499                                   Char_t *type, Bool_t force)
500 {
501     //
502     // return pointer to Q histogram
503     // if force is true create a new histogram if it doesn't exist allready
504     //
505     if ( !force || arr->UncheckedAt(sector) )
506         return (TH2S*)arr->UncheckedAt(sector);
507
508     // if we are forced and histogram doesn't yes exist create it
509     Char_t name[255], title[255];
510
511     sprintf(name,"hCalib%s%.2d",type,sector);
512     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
513
514     // new histogram with Q calib information. One value for each pad!
515     TH2S* hist = new TH2S(name,title,
516                           nbinsY, ymin, ymax,
517                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
518     hist->SetDirectory(0);
519     arr->AddAt(hist,sector);
520     return hist;
521 }
522 //_____________________________________________________________________
523 TH2S* AliTPCCalibSignal::GetHistoT0(Int_t sector, Bool_t force) /*FOLD00*/
524 {
525     //
526     // return pointer to T0 histogram
527     // if force is true create a new histogram if it doesn't exist allready
528     //
529     TObjArray *arr = &fHistoT0Array;
530     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
531 }
532 //_____________________________________________________________________
533 TH2S* AliTPCCalibSignal::GetHistoQ(Int_t sector, Bool_t force) /*FOLD00*/
534 {
535     //
536     // return pointer to Q histogram
537     // if force is true create a new histogram if it doesn't exist allready
538     //
539     TObjArray *arr = &fHistoQArray;
540     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
541 }
542 //_____________________________________________________________________
543 TH2S* AliTPCCalibSignal::GetHistoRMS(Int_t sector, Bool_t force) /*FOLD00*/
544 {
545     //
546     // return pointer to Q histogram
547     // if force is true create a new histogram if it doesn't exist allready
548     //
549     TObjArray *arr = &fHistoRMSArray;
550     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
551 }
552 //_____________________________________________________________________
553 TVectorF* AliTPCCalibSignal::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force) /*FOLD00*/
554 {
555     //
556     // return pointer to Pad Info from 'arr' for the current event and sector
557     // if force is true create it if it doesn't exist allready
558     //
559     if ( !force || arr->UncheckedAt(sector) )
560         return (TVectorF*)arr->UncheckedAt(sector);
561
562     TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
563     arr->AddAt(vect,sector);
564     return vect;
565 }
566 //_____________________________________________________________________
567 TVectorF* AliTPCCalibSignal::GetPadTimesEvent(Int_t sector, Bool_t force) /*FOLD00*/
568 {
569     //
570     // return pointer to Pad Times Array for the current event and sector
571     // if force is true create it if it doesn't exist allready
572     //
573     TObjArray *arr = &fPadTimesArrayEvent;
574     return GetPadInfoEvent(sector,arr,force);
575 }
576 //_____________________________________________________________________
577 TVectorF* AliTPCCalibSignal::GetPadQEvent(Int_t sector, Bool_t force) /*FOLD00*/
578 {
579     //
580     // return pointer to Pad Q Array for the current event and sector
581     // if force is true create it if it doesn't exist allready
582     // for debugging purposes only
583     //
584
585     TObjArray *arr = &fPadQArrayEvent;
586     return GetPadInfoEvent(sector,arr,force);
587 }
588 //_____________________________________________________________________
589 TVectorF* AliTPCCalibSignal::GetPadRMSEvent(Int_t sector, Bool_t force) /*FOLD00*/
590 {
591     //
592     // return pointer to Pad RMS Array for the current event and sector
593     // if force is true create it if it doesn't exist allready
594     // for debugging purposes only
595     //
596     TObjArray *arr = &fPadRMSArrayEvent;
597     return GetPadInfoEvent(sector,arr,force);
598 }
599 //_____________________________________________________________________
600 TVectorF* AliTPCCalibSignal::GetPadPedestalEvent(Int_t sector, Bool_t force) /*FOLD00*/
601 {
602     //
603     // return pointer to Pad RMS Array for the current event and sector
604     // if force is true create it if it doesn't exist allready
605     // for debugging purposes only
606     //
607     TObjArray *arr = &fPadPedestalArrayEvent;
608     return GetPadInfoEvent(sector,arr,force);
609 }
610 //_____________________________________________________________________
611 AliTPCCalROC* AliTPCCalibSignal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
612 {
613     //
614     // return pointer to ROC Calibration
615     // if force is true create a new histogram if it doesn't exist allready
616     //
617     if ( !force || arr->UncheckedAt(sector) )
618         return (AliTPCCalROC*)arr->UncheckedAt(sector);
619
620     // if we are forced and histogram doesn't yes exist create it
621
622     // new AliTPCCalROC for T0 information. One value for each pad!
623     AliTPCCalROC *croc = new AliTPCCalROC(sector);
624     //init values
625     for ( UInt_t iChannel = 0; iChannel<croc->GetNchannels(); iChannel++){
626         croc->SetValue(iChannel, 0);
627     }
628     arr->AddAt(croc,sector);
629     return croc;
630 }
631 //_____________________________________________________________________
632 AliTPCCalROC* AliTPCCalibSignal::GetCalRocT0(Int_t sector, Bool_t force) /*FOLD00*/
633 {
634     //
635     // return pointer to Carge ROC Calibration
636     // if force is true create a new histogram if it doesn't exist allready
637     //
638     TObjArray *arr = &fCalRocArrayT0;
639     return GetCalRoc(sector, arr, force);
640 }
641 //_____________________________________________________________________
642 AliTPCCalROC* AliTPCCalibSignal::GetCalRocQ(Int_t sector, Bool_t force) /*FOLD00*/
643 {
644     //
645     // return pointer to T0 ROC Calibration
646     // if force is true create a new histogram if it doesn't exist allready
647     //
648     TObjArray *arr = &fCalRocArrayQ;
649     return GetCalRoc(sector, arr, force);
650 }
651 //_____________________________________________________________________
652 AliTPCCalROC* AliTPCCalibSignal::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
653 {
654     //
655     // return pointer to signal width ROC Calibration
656     // if force is true create a new histogram if it doesn't exist allready
657     //
658     TObjArray *arr = &fCalRocArrayRMS;
659     return GetCalRoc(sector, arr, force);
660 }
661 //_____________________________________________________________________
662 AliTPCCalROC* AliTPCCalibSignal::GetCalRocOutliers(Int_t sector, Bool_t force)
663 {
664     //
665     // return pointer to Outliers
666     // if force is true create a new histogram if it doesn't exist allready
667     //
668     TObjArray *arr = &fCalRocArrayOutliers;
669     return GetCalRoc(sector, arr, force);
670 }
671 //_____________________________________________________________________
672 void AliTPCCalibSignal::ResetEvent() /*FOLD00*/
673 {
674     //
675     //  Reset global counters  -- Should be called before each event is processed
676     //
677     fLastSector=-1;
678     fCurrentSector=-1;
679     fCurrentRow=-1;
680     fCurrentChannel=-1;
681
682     ResetPad();
683
684     fPadTimesArrayEvent.Delete();
685     fPadQArrayEvent.Delete();
686     fPadRMSArrayEvent.Delete();
687     fPadPedestalArrayEvent.Delete();
688
689     for ( Int_t i=0; i<72; i++ ){
690         fVTime0Offset1[i]=0;
691         fVTime0Offset1Counter[i]=0;
692     }
693 }
694 //_____________________________________________________________________
695 void AliTPCCalibSignal::ResetPad() /*FOLD00*/
696 {
697     //
698     //  Reset pad infos -- Should be called after a pad has been processed
699     //
700     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
701         fPadSignal[i] = 0;
702     fMaxTimeBin = -1;
703     fMaxPadSignal = -1;
704 }
705 //_____________________________________________________________________
706 void AliTPCCalibSignal::Analyse() /*FOLD00*/
707 {
708     //
709     //  Calculate calibration constants
710     //
711
712     TVectorD paramQ(3);
713     TVectorD paramT0(3);
714     TVectorD paramRMS(3);
715     TMatrixD dummy(3,3);
716
717     for (Int_t iSec=0; iSec<72; iSec++){
718         TH2S *hT0 = GetHistoT0(iSec);
719         if (!hT0 ) continue;
720
721         AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
722         AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
723         AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
724         AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
725
726         TH2S *hQ   = GetHistoQ(iSec);
727         TH2S *hRMS = GetHistoRMS(iSec);
728
729         Short_t *array_hQ   = hQ->GetArray();
730         Short_t *array_hT0  = hT0->GetArray();
731         Short_t *array_hRMS = hRMS->GetArray();
732
733         UInt_t nChannels = fROC->GetNChannels(iSec);
734
735         //debug
736         Int_t row=0;
737         Int_t pad=0;
738         Int_t padc=0;
739         //! debug
740
741         for (UInt_t iChannel=0; iChannel<nChannels; iChannel++){
742
743
744             Float_t cogTime0 = -1000;
745             Float_t cogQ     = -1000;
746             Float_t cogRMS   = -1000;
747             Float_t cogOut   = 0;
748
749
750             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
751             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
752             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
753
754 /*
755             AliMathBase::FitGaus(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
756             AliMathBase::FitGaus(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
757             AliMathBase::FitGaus(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
758             cogQ     = paramQ[1];
759             cogTime0 = paramT0[1];
760             cogRMS   = paramRMS[1];
761 */
762             cogQ     = AliMathBase::GetCOG(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
763             cogTime0 = AliMathBase::GetCOG(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
764             cogRMS   = AliMathBase::GetCOG(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
765
766
767
768             /*
769             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
770                 cogOut = 1;
771                 cogTime0 = 0;
772                 cogQ     = 0;
773                 cogRMS   = 0;
774             }
775 */
776             rocQ->SetValue(iChannel, cogQ*cogQ);
777             rocT0->SetValue(iChannel, cogTime0);
778             rocRMS->SetValue(iChannel, cogRMS);
779             rocOut->SetValue(iChannel, cogOut);
780
781
782             //debug
783             if ( fDebugLevel > 0 ){
784                 if ( !fDebugStreamer ) {
785                         //debug stream
786                     TDirectory *backup = gDirectory;
787                     fDebugStreamer = new TTreeSRedirector("deb2.root");
788                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
789                 }
790
791                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
792                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
793                 padc = pad-(fROC->GetNPads(iSec,row)/2);
794
795                 (*fDebugStreamer) << "DataEnd" <<
796                     "Sector="  << iSec      <<
797                     "Pad="     << pad       <<
798                     "PadC="    << padc      <<
799                     "Row="     << row       <<
800                     "PadSec="  << iChannel   <<
801                     "Q="       << cogQ      <<
802                     "T0="      << cogTime0  <<
803                     "RMS="     << cogRMS    <<
804                     "\n";
805             }
806             //! debug
807
808         }
809
810     }
811     delete fDebugStreamer;
812     fDebugStreamer = 0x0;
813 }
814 //_____________________________________________________________________
815 void AliTPCCalibSignal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
816 {
817     //
818     //  Write class to file
819     //
820
821     TString sDir(dir);
822     TString option;
823
824     if ( append )
825         option = "update";
826     else
827         option = "recreate";
828
829     TDirectory *backup = gDirectory;
830     TFile f(filename,option.Data());
831     f.cd();
832     if ( !sDir.IsNull() ){
833         f.mkdir(sDir.Data());
834         f.cd(sDir);
835     }
836     this->Write();
837     f.Close();
838
839     if ( backup ) backup->cd();
840 }