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