2faef43afe65fea73c0ae22064450c77f84e13cf
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibPulser.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //Root includes
19 #include <TH1F.h>
20 #include <TH2S.h>
21 #include <TString.h>
22 #include <TVectorF.h>
23 #include <TMath.h>
24
25 #include <TDirectory.h>
26 #include <TSystem.h>
27 #include <TFile.h>
28
29 //AliRoot includes
30 #include "AliRawReader.h"
31 #include "AliRawReaderRoot.h"
32 #include "AliRawReaderDate.h"
33 #include "AliTPCRawStream.h"
34 #include "AliTPCCalROC.h"
35 #include "AliTPCCalPad.h"
36 #include "AliTPCROC.h"
37 #include "AliTPCParam.h"
38 #include "AliTPCCalibPulser.h"
39 #include "AliTPCcalibDB.h"
40 #include "AliMathBase.h"
41 #include "TTreeStream.h"
42
43 //date
44 #include "event.h"
45
46
47 ///////////////////////////////////////////////////////////////////////////////////////
48 //          Implementation of the TPC pulser calibration
49 //
50 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
51 // 
52 // 
53 /***************************************************************************
54  *                      Class Description                                  *
55  ***************************************************************************
56
57  The AliTPCCalibPulser class is used to get calibration data concerning the FEE using
58  runs performed with the calibration pulser.
59
60  The information retrieved is
61  - Time0 differences
62  - Signal width differences
63  - Amplification variations
64
65  the seen differences arise from the manufacturing tolerances of the PASAs and are very small within
66  one chip and somewhat large between different chips.
67
68  Histograms:
69    For each ROC three TH2S histos 'Reference Histograms'  (ROC channel vs. [Time0, signal width, Q sum]) is created when
70    it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are stored in the
71    TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.
72
73
74  Working principle:
75  ------------------
76  Raw calibration pulser data is processed by calling one of the ProcessEvent(...) functions
77  (see below). These in the end call the Update(...) function.
78
79  - the Update(...) function:
80    In this function the array fPadSignal is filled with the adc signals between the specified range
81    fFirstTimeBin and fLastTimeBin for the current pad.
82    before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
83    stored in fPadSignal.
84
85    - the ProcessPad() function:
86      Find Pedestal and Noise information
87      - use database information which has to be set by calling
88        SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)
89      - if no information from the pedestal data base
90        is available the informaion is calculated on the fly ( see FindPedestal() function )
91
92      Find the Pulser signal information
93      - calculate  mean = T0, RMS = signal width and Q sum in a range of -2+7 timebins around Q max
94        the Q sum is scaled by pad area
95        (see FindPulserSignal(...) function)
96
97      Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)
98      Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE)),
99
100  At the end of each event the EndEvent() function is called
101
102  - the EndEvent() function:
103    calculate the mean T0 for each ROC and fill the Time0 histogram with Time0-<Time0 for ROC>
104    This is done to overcome syncronisation problems between the trigger and the fec clock.
105
106  After accumulating the desired statistics the Analyse() function has to be called.
107  - the Analyse() function
108    Whithin this function the mean values of T0, RMS, Q are calculated for each pad, using
109    the AliMathBase::GetCOG(...) function, and the calibration
110    storage classes (AliTPCCalROC) are filled for each ROC.
111    The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and
112    fCalRocArrayQ;
113
114
115
116  User interface for filling data:
117  --------------------------------
118
119  To Fill information one of the following functions can be used:
120
121  Bool_t ProcessEvent(eventHeaderStruct *event);
122    - process Date event
123    - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
124
125  Bool_t ProcessEvent(AliRawReader *rawReader);
126    - process AliRawReader event
127    - use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)
128
129  Bool_t ProcessEvent(AliTPCRawStream *rawStream);
130    - process event from AliTPCRawStream
131    - call Update function for signal filling
132
133  Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
134               iPad,  const Int_t iTimeBin, const Float_t signal);
135    - directly  fill signal information (sector, row, pad, time bin, pad)
136      to the reference histograms
137
138  It is also possible to merge two independently taken calibrations using the function
139
140  void Merge(AliTPCCalibPulser *sig)
141    - copy histograms in 'sig' if the do not exist in this instance
142    - Add histograms in 'sig' to the histograms in this instance if the allready exist
143    - After merging call Analyse again!
144
145
146
147  -- example: filling data using root raw data:
148  void fillSignal(Char_t *filename)
149  {
150     rawReader = new AliRawReaderRoot(fileName);
151     if ( !rawReader ) return;
152     AliTPCCalibPulser *calib = new AliTPCCalibPulser;
153     while (rawReader->NextEvent()){
154       calib->ProcessEvent(rawReader);
155     }
156     calib->Analyse();
157     calib->DumpToFile("SignalData.root");
158     delete rawReader;
159     delete calib;
160  }
161
162
163  What kind of information is stored and how to retrieve them:
164  ------------------------------------------------------------
165
166  - Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):
167
168    TH2F *GetHistoT0(Int_t sector);
169    TH2F *GetHistoRMS(Int_t sector);
170    TH2F *GetHistoQ(Int_t sector);
171
172  - Accessing the calibration storage objects:
173
174    AliTPCCalROC *GetCalRocT0(Int_t sector);   // for the Time0 values
175    AliTPCCalROC *GetCalRocRMS(Int_t sector);  // for the signal width values
176    AliTPCCalROC *GetCalRocQ(Int_t sector);    // for the Q sum values
177
178    example for visualisation:
179    if the file "SignalData.root" was created using the above example one could do the following:
180
181    TFile fileSignal("SignalData.root")
182    AliTPCCalibPulser *sig = (AliTPCCalibPulser*)fileSignal->Get("AliTPCCalibPulser");
183    sig->GetCalRocT0(0)->Draw("colz");
184    sig->GetCalRocRMS(0)->Draw("colz");
185
186    or use the AliTPCCalPad functionality:
187    AliTPCCalPad padT0(ped->GetCalPadT0());
188    AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
189    padT0->MakeHisto2D()->Draw("colz");       //Draw A-Side Time0 Information
190    padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
191 */
192
193
194 ClassImp(AliTPCCalibPulser) /*FOLD00*/
195
196 AliTPCCalibPulser::AliTPCCalibPulser() : /*FOLD00*/
197     TObject(),
198     fFirstTimeBin(60),
199     fLastTimeBin(120),
200     fNbinsT0(200),
201     fXminT0(-2),
202     fXmaxT0(2),
203     fNbinsQ(200),
204     fXminQ(1),
205     fXmaxQ(40),
206     fNbinsRMS(100),
207     fXminRMS(0.1),
208     fXmaxRMS(5.1),
209     fLastSector(-1),
210     fOldRCUformat(kTRUE),
211     fROC(AliTPCROC::Instance()),
212     fParam(new AliTPCParam),
213     fPedestalTPC(0x0),
214     fPadNoiseTPC(0x0),
215     fPedestalROC(0x0),
216     fPadNoiseROC(0x0),
217     fCalRocArrayT0(72),
218     fCalRocArrayQ(72),
219     fCalRocArrayRMS(72),
220     fCalRocArrayOutliers(72),
221     fHistoQArray(72),
222     fHistoT0Array(72),
223     fHistoRMSArray(72),
224     fPadTimesArrayEvent(72),
225     fPadQArrayEvent(72),
226     fPadRMSArrayEvent(72),
227     fPadPedestalArrayEvent(72),
228     fCurrentChannel(-1),
229     fCurrentSector(-1),
230     fCurrentRow(-1),
231     fMaxPadSignal(-1),
232     fMaxTimeBin(-1),
233     fPadSignal(1024),
234     fPadPedestal(0),
235     fPadNoise(0),
236     fVTime0Offset(72),
237     fVTime0OffsetCounter(72),
238     fEvent(-1),
239     fDebugStreamer(0x0),
240     fDebugLevel(0)
241 {
242     //
243     // AliTPCSignal default constructor
244     //
245
246 }
247 //_____________________________________________________________________
248 AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
249     TObject(sig),
250     fFirstTimeBin(sig.fFirstTimeBin),
251     fLastTimeBin(sig.fLastTimeBin),
252     fNbinsT0(sig.fNbinsT0),
253     fXminT0(sig.fXminT0),
254     fXmaxT0(sig.fXmaxT0),
255     fNbinsQ(sig.fNbinsQ),
256     fXminQ(sig.fXminQ),
257     fXmaxQ(sig.fXmaxQ),
258     fNbinsRMS(sig.fNbinsRMS),
259     fXminRMS(sig.fXminRMS),
260     fXmaxRMS(sig.fXmaxRMS),
261     fLastSector(-1),
262     fOldRCUformat(kTRUE),
263     fROC(AliTPCROC::Instance()),
264     fParam(new AliTPCParam),
265     fPedestalTPC(0x0),
266     fPadNoiseTPC(0x0),
267     fPedestalROC(0x0),
268     fPadNoiseROC(0x0),
269     fCalRocArrayT0(72),
270     fCalRocArrayQ(72),
271     fCalRocArrayRMS(72),
272     fCalRocArrayOutliers(72),
273     fHistoQArray(72),
274     fHistoT0Array(72),
275     fHistoRMSArray(72),
276     fPadTimesArrayEvent(72),
277     fPadQArrayEvent(72),
278     fPadRMSArrayEvent(72),
279     fPadPedestalArrayEvent(72),
280     fCurrentChannel(-1),
281     fCurrentSector(-1),
282     fCurrentRow(-1),
283     fMaxPadSignal(-1),
284     fMaxTimeBin(-1),
285     fPadSignal(1024),
286     fPadPedestal(0),
287     fPadNoise(0),
288     fVTime0Offset(72),
289     fVTime0OffsetCounter(72),
290     fEvent(-1),
291     fDebugStreamer(0x0),
292     fDebugLevel(sig.fDebugLevel)
293 {
294     //
295     // AliTPCSignal default constructor
296     //
297
298     for (Int_t iSec = 0; iSec < 72; ++iSec){
299         const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
300         const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
301         const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
302         const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
303
304         const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
305         const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
306         const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
307
308         if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
309         if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
310         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
311         if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
312
313         if ( hQ != 0x0 ){
314             TH2S *hNew = new TH2S(*hQ);
315             hNew->SetDirectory(0);
316             fHistoQArray.AddAt(hNew,iSec);
317         }
318         if ( hT0 != 0x0 ){
319             TH2S *hNew = new TH2S(*hT0);
320             hNew->SetDirectory(0);
321             fHistoQArray.AddAt(hNew,iSec);
322         }
323         if ( hRMS != 0x0 ){
324             TH2S *hNew = new TH2S(*hRMS);
325             hNew->SetDirectory(0);
326             fHistoQArray.AddAt(hNew,iSec);
327         }
328     }
329
330 }
331 //_____________________________________________________________________
332 AliTPCCalibPulser& AliTPCCalibPulser::operator = (const  AliTPCCalibPulser &source)
333 {
334   //
335   // assignment operator
336   //
337   if (&source == this) return *this;
338   new (this) AliTPCCalibPulser(source);
339
340   return *this;
341 }
342 //_____________________________________________________________________
343 AliTPCCalibPulser::~AliTPCCalibPulser()
344 {
345     //
346     // destructor
347     //
348
349     fCalRocArrayT0.Delete();
350     fCalRocArrayQ.Delete();
351     fCalRocArrayRMS.Delete();
352     fCalRocArrayOutliers.Delete();
353
354     fHistoQArray.Delete();
355     fHistoT0Array.Delete();
356     fHistoRMSArray.Delete();
357
358     fPadTimesArrayEvent.Delete();
359     fPadQArrayEvent.Delete();
360     fPadRMSArrayEvent.Delete();
361     fPadPedestalArrayEvent.Delete();
362
363     if ( fDebugStreamer) delete fDebugStreamer;
364     delete fParam;
365 }
366 //_____________________________________________________________________
367 Int_t AliTPCCalibPulser::Update(const Int_t icsector, /*FOLD00*/
368                                 const Int_t icRow,
369                                 const Int_t icPad,
370                                 const Int_t icTimeBin,
371                                 const Float_t csignal)
372 {
373     //
374     // Signal filling methode on the fly pedestal and Time offset correction if necessary.
375     // no extra analysis necessary. Assumes knowledge of the signal shape!
376     // assumes that it is looped over consecutive time bins of one pad
377     //
378     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
379
380     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
381
382     //init first pad and sector in this event
383     if ( fCurrentChannel == -1 ) {
384         fCurrentChannel = iChannel;
385         fCurrentSector  = icsector;
386         fCurrentRow     = icRow;
387     }
388
389     //process last pad if we change to a new one
390     if ( iChannel != fCurrentChannel ){
391         ProcessPad();
392         fCurrentChannel = iChannel;
393         fCurrentSector  = icsector;
394         fCurrentRow     = icRow;
395     }
396
397     //fill signals for current pad
398     fPadSignal[icTimeBin]=csignal;
399     if ( csignal > fMaxPadSignal ){
400         fMaxPadSignal = csignal;
401         fMaxTimeBin   = icTimeBin;
402     }
403     return 0;
404 }
405 //_____________________________________________________________________
406 void AliTPCCalibPulser::FindPedestal(Float_t part)
407 {
408     //
409     // find pedestal and noise for the current pad. Use either database or
410     // truncated mean with part*100%
411     //
412     Bool_t noPedestal = kTRUE;;
413     if (fPedestalTPC&&fPadNoiseTPC){
414         //use pedestal database
415         //only load new pedestals if the sector has changed
416         if ( fCurrentSector!=fLastSector ){
417             fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
418             fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
419             fLastSector=fCurrentSector;
420         }
421
422         if ( fPedestalROC&&fPadNoiseROC ){
423             fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
424             fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
425             noPedestal   = kFALSE;
426         }
427
428     }
429
430     //if we are not running with pedestal database, or for the current sector there is no information
431     //available, calculate the pedestal and noise on the fly
432     if ( noPedestal ) {
433         const Int_t kPedMax = 100;  //maximum pedestal value
434         Float_t  max    =  0;
435         Float_t  maxPos =  0;
436         Int_t    median =  -1;
437         Int_t    count0 =  0;
438         Int_t    count1 =  0;
439         //
440         Float_t padSignal=0;
441         //
442         UShort_t histo[kPedMax];
443         memset(histo,0,kPedMax*sizeof(UShort_t));
444
445         for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
446             padSignal = fPadSignal.GetMatrixArray()[i];
447             if (padSignal<=0) continue;
448             if (padSignal>max && i>10) {
449                 max = padSignal;
450                 maxPos = i;
451             }
452             if (padSignal>kPedMax-1) continue;
453             histo[Int_t(padSignal+0.5)]++;
454             count0++;
455         }
456             //
457         for (Int_t i=1; i<kPedMax; ++i){
458             if (count1<count0*0.5) median=i;
459             count1+=histo[i];
460         }
461         // truncated mean
462         //
463         Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
464         //
465         for (Int_t idelta=1; idelta<10; ++idelta){
466             if (median-idelta<=0) continue;
467             if (median+idelta>kPedMax) continue;
468             if (count<part*count1){
469                 count+=histo[median-idelta];
470                 mean +=histo[median-idelta]*(median-idelta);
471                 rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
472                 count+=histo[median+idelta];
473                 mean +=histo[median+idelta]*(median+idelta);
474                 rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
475             }
476         }
477         fPadPedestal = 0;
478         fPadNoise    = 0;
479         if ( count > 0 ) {
480             mean/=count;
481             rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
482             fPadPedestal = mean;
483             fPadNoise    = rms;
484         } 
485     }
486 }
487 //_____________________________________________________________________
488 void AliTPCCalibPulser::FindPulserSignal(TVectorD &param, Float_t &qSum)
489 {
490     //
491     //  Find position, signal width and height of the CE signal (last signal)
492     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
493     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
494     //
495
496     Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
497     Int_t   cemaxpos       = fMaxTimeBin;
498     Float_t ceSumThreshold = 8.*fPadNoise;  // threshold for the signal sum
499     const Int_t    kCemin  = 2;             // range for the analysis of the ce signal +- channels from the peak
500     const Int_t    kCemax  = 7;
501
502     if (cemaxpos!=0){
503         ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
504         for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
505             Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
506             if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
507                 ceTime+=signal*(i+0.5);
508                 ceRMS +=signal*(i+0.5)*(i+0.5);
509                 ceQsum+=signal;
510             }
511         }
512     }
513     if (ceQmax&&ceQsum>ceSumThreshold) {
514         ceTime/=ceQsum;
515         ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
516         fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
517         fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
518
519         //Normalise Q to the pad area
520         Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
521
522         ceQsum/=norm;
523     } else {
524         ceQmax=0;
525         ceTime=0;
526         ceRMS =0;
527         ceQsum=0;
528     }
529     param[0] = ceQmax;
530     param[1] = ceTime;
531     param[2] = ceRMS;
532     qSum     = ceQsum;
533 }
534 //_____________________________________________________________________
535 void AliTPCCalibPulser::ProcessPad() /*FOLD00*/
536 {
537     //
538     //  Process data of current pad
539     //
540
541     FindPedestal();
542     TVectorD param(3);
543     Float_t  Qsum;
544     FindPulserSignal(param, Qsum);
545
546     Double_t meanT  = param[1];
547     Double_t sigmaT = param[2];
548
549     //Fill Event T0 counter
550     (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
551
552     //Fill Q histogram
553     GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(Qsum), fCurrentChannel );
554
555     //Fill RMS histogram
556     GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
557
558
559     //Fill debugging info
560     if ( fDebugLevel>0 ){
561         (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
562         (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
563         (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=Qsum;
564     }
565
566     ResetPad();
567 }
568 //_____________________________________________________________________
569 void AliTPCCalibPulser::EndEvent() /*FOLD00*/
570 {
571     //
572     //  Process data of current event
573     //
574
575     //check if last pad has allready been processed, if not do so
576     if ( fMaxTimeBin>-1 ) ProcessPad();
577
578     //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
579     for ( Int_t iSec = 0; iSec<72; ++iSec ){
580         TVectorF *vTimes = GetPadTimesEvent(iSec);
581         if ( !vTimes ) continue;
582
583         for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
584             Float_t Time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
585             Float_t Time  = (*vTimes).GetMatrixArray()[iChannel];
586
587             GetHistoT0(iSec,kTRUE)->Fill( Time-Time0,iChannel );
588
589
590             //Debug start
591             if ( fDebugLevel>0 ){
592                 if ( !fDebugStreamer ) {
593                         //debug stream
594                     TDirectory *backup = gDirectory;
595                     fDebugStreamer = new TTreeSRedirector("deb2.root");
596                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
597                 }
598
599                 Int_t row=0;
600                 Int_t pad=0;
601                 Int_t padc=0;
602
603                 Float_t Q   = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
604                 Float_t RMS = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
605
606                 UInt_t channel=iChannel;
607                 Int_t sector=iSec;
608
609                 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
610                 pad = channel-fROC->GetRowIndexes(sector)[row];
611                 padc = pad-(fROC->GetNPads(sector,row)/2);
612
613                 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
614                                     Form("hSignalD%d.%d.%d",sector,row,pad),
615                                     fLastTimeBin-fFirstTimeBin,
616                                     fFirstTimeBin,fLastTimeBin);
617                 h1->SetDirectory(0);
618
619                 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
620                     h1->Fill(i,fPadSignal(i));
621
622                 (*fDebugStreamer) << "DataPad" <<
623                     "Event=" << fEvent <<
624                     "Sector="<< sector <<
625                     "Row="   << row<<
626                     "Pad="   << pad <<
627                     "PadC="  << padc <<
628                     "PadSec="<< channel <<
629                     "Time0="  << Time0 <<
630                     "Time="  << Time <<
631                     "RMS="   << RMS <<
632                     "Sum="   << Q <<
633                     "hist.=" << h1 <<
634                     "\n";
635
636                 delete h1;
637             }
638             //Debug end
639
640         }
641     }
642
643 }
644 //_____________________________________________________________________
645 Bool_t AliTPCCalibPulser::ProcessEvent(AliTPCRawStream *rawStream) /*FOLD00*/
646 {
647   //
648   // Event Processing loop - AliTPCRawStream
649   //
650
651   rawStream->SetOldRCUFormat(fOldRCUformat);
652
653   ResetEvent();
654
655   Bool_t withInput = kFALSE;
656
657   while (rawStream->Next()) {
658
659       Int_t isector  = rawStream->GetSector();                       //  current sector
660       Int_t iRow     = rawStream->GetRow();                          //  current row
661       Int_t iPad     = rawStream->GetPad();                          //  current pad
662       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
663       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
664
665       Update(isector,iRow,iPad,iTimeBin,signal);
666       withInput = kTRUE;
667   }
668
669   if (withInput){
670       EndEvent();
671   }
672
673   return withInput;
674 }
675 //_____________________________________________________________________
676 Bool_t AliTPCCalibPulser::ProcessEvent(AliRawReader *rawReader)
677 {
678   //
679   //  Event processing loop - AliRawReader
680   //
681
682
683   AliTPCRawStream rawStream(rawReader);
684
685   rawReader->Select("TPC");
686
687   return ProcessEvent(&rawStream);
688 }
689 //_____________________________________________________________________
690 Bool_t AliTPCCalibPulser::ProcessEvent(eventHeaderStruct *event)
691 {
692   //
693   //  Event processing loop - date event
694   //
695     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
696     Bool_t result=ProcessEvent(rawReader);
697     delete rawReader;
698     return result;
699
700 }
701 //_____________________________________________________________________
702 TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
703                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
704                                   Char_t *type, Bool_t force)
705 {
706     //
707     // return pointer to Q histogram
708     // if force is true create a new histogram if it doesn't exist allready
709     //
710     if ( !force || arr->UncheckedAt(sector) )
711         return (TH2S*)arr->UncheckedAt(sector);
712
713     // if we are forced and histogram doesn't yes exist create it
714     Char_t name[255], title[255];
715
716     sprintf(name,"hCalib%s%.2d",type,sector);
717     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
718
719     // new histogram with Q calib information. One value for each pad!
720     TH2S* hist = new TH2S(name,title,
721                           nbinsY, ymin, ymax,
722                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
723     hist->SetDirectory(0);
724     arr->AddAt(hist,sector);
725     return hist;
726 }
727 //_____________________________________________________________________
728 TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force) /*FOLD00*/
729 {
730     //
731     // return pointer to T0 histogram
732     // if force is true create a new histogram if it doesn't exist allready
733     //
734     TObjArray *arr = &fHistoT0Array;
735     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
736 }
737 //_____________________________________________________________________
738 TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force) /*FOLD00*/
739 {
740     //
741     // return pointer to Q histogram
742     // if force is true create a new histogram if it doesn't exist allready
743     //
744     TObjArray *arr = &fHistoQArray;
745     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
746 }
747 //_____________________________________________________________________
748 TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force) /*FOLD00*/
749 {
750     //
751     // return pointer to Q histogram
752     // if force is true create a new histogram if it doesn't exist allready
753     //
754     TObjArray *arr = &fHistoRMSArray;
755     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
756 }
757 //_____________________________________________________________________
758 TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force) /*FOLD00*/
759 {
760     //
761     // return pointer to Pad Info from 'arr' for the current event and sector
762     // if force is true create it if it doesn't exist allready
763     //
764     if ( !force || arr->UncheckedAt(sector) )
765         return (TVectorF*)arr->UncheckedAt(sector);
766
767     TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
768     arr->AddAt(vect,sector);
769     return vect;
770 }
771 //_____________________________________________________________________
772 TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force) /*FOLD00*/
773 {
774     //
775     // return pointer to Pad Times Array for the current event and sector
776     // if force is true create it if it doesn't exist allready
777     //
778     TObjArray *arr = &fPadTimesArrayEvent;
779     return GetPadInfoEvent(sector,arr,force);
780 }
781 //_____________________________________________________________________
782 TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force) /*FOLD00*/
783 {
784     //
785     // return pointer to Pad Q Array for the current event and sector
786     // if force is true create it if it doesn't exist allready
787     // for debugging purposes only
788     //
789
790     TObjArray *arr = &fPadQArrayEvent;
791     return GetPadInfoEvent(sector,arr,force);
792 }
793 //_____________________________________________________________________
794 TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force) /*FOLD00*/
795 {
796     //
797     // return pointer to Pad RMS Array for the current event and sector
798     // if force is true create it if it doesn't exist allready
799     // for debugging purposes only
800     //
801     TObjArray *arr = &fPadRMSArrayEvent;
802     return GetPadInfoEvent(sector,arr,force);
803 }
804 //_____________________________________________________________________
805 TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force) /*FOLD00*/
806 {
807     //
808     // return pointer to Pad RMS Array for the current event and sector
809     // if force is true create it if it doesn't exist allready
810     // for debugging purposes only
811     //
812     TObjArray *arr = &fPadPedestalArrayEvent;
813     return GetPadInfoEvent(sector,arr,force);
814 }
815 //_____________________________________________________________________
816 AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
817 {
818     //
819     // return pointer to ROC Calibration
820     // if force is true create a new histogram if it doesn't exist allready
821     //
822     if ( !force || arr->UncheckedAt(sector) )
823         return (AliTPCCalROC*)arr->UncheckedAt(sector);
824
825     // if we are forced and histogram doesn't yes exist create it
826
827     // new AliTPCCalROC for T0 information. One value for each pad!
828     AliTPCCalROC *croc = new AliTPCCalROC(sector);
829     arr->AddAt(croc,sector);
830     return croc;
831 }
832 //_____________________________________________________________________
833 AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force) /*FOLD00*/
834 {
835     //
836     // return pointer to Carge ROC Calibration
837     // if force is true create a new histogram if it doesn't exist allready
838     //
839     TObjArray *arr = &fCalRocArrayT0;
840     return GetCalRoc(sector, arr, force);
841 }
842 //_____________________________________________________________________
843 AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force) /*FOLD00*/
844 {
845     //
846     // return pointer to T0 ROC Calibration
847     // if force is true create a new histogram if it doesn't exist allready
848     //
849     TObjArray *arr = &fCalRocArrayQ;
850     return GetCalRoc(sector, arr, force);
851 }
852 //_____________________________________________________________________
853 AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
854 {
855     //
856     // return pointer to signal width ROC Calibration
857     // if force is true create a new histogram if it doesn't exist allready
858     //
859     TObjArray *arr = &fCalRocArrayRMS;
860     return GetCalRoc(sector, arr, force);
861 }
862 //_____________________________________________________________________
863 AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
864 {
865     //
866     // return pointer to Outliers
867     // if force is true create a new histogram if it doesn't exist allready
868     //
869     TObjArray *arr = &fCalRocArrayOutliers;
870     return GetCalRoc(sector, arr, force);
871 }
872 //_____________________________________________________________________
873 void AliTPCCalibPulser::ResetEvent() /*FOLD00*/
874 {
875     //
876     //  Reset global counters  -- Should be called before each event is processed
877     //
878     fLastSector=-1;
879     fCurrentSector=-1;
880     fCurrentRow=-1;
881     fCurrentChannel=-1;
882
883     ResetPad();
884
885     fPadTimesArrayEvent.Delete();
886     fPadQArrayEvent.Delete();
887     fPadRMSArrayEvent.Delete();
888     fPadPedestalArrayEvent.Delete();
889
890     for ( Int_t i=0; i<72; ++i ){
891         fVTime0Offset[i]=0;
892         fVTime0OffsetCounter[i]=0;
893     }
894 }
895 //_____________________________________________________________________
896 void AliTPCCalibPulser::ResetPad() /*FOLD00*/
897 {
898     //
899     //  Reset pad infos -- Should be called after a pad has been processed
900     //
901     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
902         fPadSignal[i] = 0;
903     fMaxTimeBin = -1;
904     fMaxPadSignal = -1;
905     fPadPedestal  = -1;
906     fPadNoise     = -1;
907 }
908 //_____________________________________________________________________
909 void AliTPCCalibPulser::Merge(AliTPCCalibPulser *sig)
910 {
911     //
912     //  Merge reference histograms of sig to the current AliTPCCalibPulser
913     //
914
915     //merge histograms
916     for (Int_t iSec=0; iSec<72; ++iSec){
917         TH2S *hRefQmerge   = sig->GetHistoQ(iSec);
918         TH2S *hRefT0merge  = sig->GetHistoT0(iSec);
919         TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
920
921
922         if ( hRefQmerge ){
923             TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
924             TH2S *hRefQ   = GetHistoQ(iSec);
925             if ( hRefQ ) hRefQ->Add(hRefQmerge);
926             else {
927                 TH2S *hist = new TH2S(*hRefQmerge);
928                 hist->SetDirectory(0);
929                 fHistoQArray.AddAt(hist, iSec);
930             }
931             hRefQmerge->SetDirectory(dir);
932         }
933         if ( hRefT0merge ){
934             TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
935             TH2S *hRefT0  = GetHistoT0(iSec);
936             if ( hRefT0 ) hRefT0->Add(hRefT0merge);
937             else {
938                 TH2S *hist = new TH2S(*hRefT0merge);
939                 hist->SetDirectory(0);
940                 fHistoT0Array.AddAt(hist, iSec);
941             }
942             hRefT0merge->SetDirectory(dir);
943         }
944         if ( hRefRMSmerge ){
945             TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
946             TH2S *hRefRMS = GetHistoRMS(iSec);
947             if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
948             else {
949                 TH2S *hist = new TH2S(*hRefRMSmerge);
950                 hist->SetDirectory(0);
951                 fHistoRMSArray.AddAt(hist, iSec);
952             }
953             hRefRMSmerge->SetDirectory(dir);
954         }
955
956     }
957 }
958 //_____________________________________________________________________
959 void AliTPCCalibPulser::Analyse() /*FOLD00*/
960 {
961     //
962     //  Calculate calibration constants
963     //
964
965     TVectorD paramQ(3);
966     TVectorD paramT0(3);
967     TVectorD paramRMS(3);
968     TMatrixD dummy(3,3);
969
970     for (Int_t iSec=0; iSec<72; ++iSec){
971         TH2S *hT0 = GetHistoT0(iSec);
972         if (!hT0 ) continue;
973
974         AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
975         AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
976         AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
977         AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
978
979         TH2S *hQ   = GetHistoQ(iSec);
980         TH2S *hRMS = GetHistoRMS(iSec);
981
982         Short_t *array_hQ   = hQ->GetArray();
983         Short_t *array_hT0  = hT0->GetArray();
984         Short_t *array_hRMS = hRMS->GetArray();
985
986         UInt_t nChannels = fROC->GetNChannels(iSec);
987
988         //debug
989         Int_t row=0;
990         Int_t pad=0;
991         Int_t padc=0;
992         //! debug
993
994         for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
995
996
997             Float_t cogTime0 = -1000;
998             Float_t cogQ     = -1000;
999             Float_t cogRMS   = -1000;
1000             Float_t cogOut   = 0;
1001
1002
1003             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1004             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1005             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1006
1007 /*
1008             AliMathBase::FitGaus(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1009             AliMathBase::FitGaus(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1010             AliMathBase::FitGaus(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1011             cogQ     = paramQ[1];
1012             cogTime0 = paramT0[1];
1013             cogRMS   = paramRMS[1];
1014 */
1015             cogQ     = AliMathBase::GetCOG(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1016             cogTime0 = AliMathBase::GetCOG(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1017             cogRMS   = AliMathBase::GetCOG(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1018
1019
1020
1021             /*
1022             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1023                 cogOut = 1;
1024                 cogTime0 = 0;
1025                 cogQ     = 0;
1026                 cogRMS   = 0;
1027             }
1028 */
1029             rocQ->SetValue(iChannel, cogQ*cogQ);
1030             rocT0->SetValue(iChannel, cogTime0);
1031             rocRMS->SetValue(iChannel, cogRMS);
1032             rocOut->SetValue(iChannel, cogOut);
1033
1034
1035             //debug
1036             if ( fDebugLevel > 0 ){
1037                 if ( !fDebugStreamer ) {
1038                         //debug stream
1039                     TDirectory *backup = gDirectory;
1040                     fDebugStreamer = new TTreeSRedirector("deb2.root");
1041                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1042                 }
1043
1044                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1045                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1046                 padc = pad-(fROC->GetNPads(iSec,row)/2);
1047
1048                 (*fDebugStreamer) << "DataEnd" <<
1049                     "Sector="  << iSec      <<
1050                     "Pad="     << pad       <<
1051                     "PadC="    << padc      <<
1052                     "Row="     << row       <<
1053                     "PadSec="  << iChannel   <<
1054                     "Q="       << cogQ      <<
1055                     "T0="      << cogTime0  <<
1056                     "RMS="     << cogRMS    <<
1057                     "\n";
1058             }
1059             //! debug
1060
1061         }
1062
1063     }
1064     delete fDebugStreamer;
1065     fDebugStreamer = 0x0;
1066 }
1067 //_____________________________________________________________________
1068 void AliTPCCalibPulser::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1069 {
1070     //
1071     //  Write class to file
1072     //
1073
1074     TString sDir(dir);
1075     TString option;
1076
1077     if ( append )
1078         option = "update";
1079     else
1080         option = "recreate";
1081
1082     TDirectory *backup = gDirectory;
1083     TFile f(filename,option.Data());
1084     f.cd();
1085     if ( !sDir.IsNull() ){
1086         f.mkdir(sDir.Data());
1087         f.cd(sDir);
1088     }
1089     this->Write();
1090     f.Close();
1091
1092     if ( backup ) backup->cd();
1093 }
1094 //_____________________________________________________________________
1095 //_________________________  Test Functions ___________________________
1096 //_____________________________________________________________________
1097 TObjArray* AliTPCCalibPulser::TestBinning()
1098 {
1099     //
1100     //  Function to test the binning of the reference histograms
1101     //  type: T0, Q or RMS
1102     //  mode: 0 - number of filled bins per channel
1103     //        1 - number of empty bins between filled bins in one ROC
1104     //  returns TObjArray with the test histograms type*2+mode:
1105     //  position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1106
1107
1108     TObjArray *histArray = new TObjArray(6);
1109     const Char_t *type[] = {"T0","Q","RMS"};
1110     Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1111
1112     for (Int_t itype = 0; itype<3; ++itype){
1113         for (Int_t imode=0; imode<2; ++imode){
1114             Int_t icount = itype*2+imode;
1115             histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1116                                       Form("Test Binning of '%s', mode - %d",type[itype],imode),
1117                                       72,0,72),
1118                              icount);
1119         }
1120     }
1121
1122
1123     TH2S *hRef=0x0;
1124     Short_t *array=0x0;
1125     for (Int_t itype = 0; itype<3; ++itype){
1126         for (Int_t iSec=0; iSec<72; ++iSec){
1127             if ( itype == 0 ) hRef = GetHistoT0(iSec);
1128             if ( itype == 1 ) hRef = GetHistoQ(iSec);
1129             if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1130             if ( hRef == 0x0 ) continue;
1131             array = (hRef->GetArray());
1132             UInt_t nChannels = fROC->GetNChannels(iSec);
1133
1134             Int_t nempty=0;
1135             for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1136                 Int_t nfilled=0;
1137                 Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1138                 Int_t c1 = 0;
1139                 Int_t c2 = 0;
1140                 for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1141                     if ( array[offset+iBin]>0 ) {
1142                         nfilled++;
1143                         if ( c1 && c2 ) nempty++;
1144                         else c1 = 1;
1145                     }
1146                     else if ( c1 ) c2 = 1;
1147
1148                 }
1149                 ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1150             }
1151             ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);
1152         }
1153     }
1154     return histArray;
1155 }