]> git.uio.no Git - u/mrichter/AliRoot.git/blob - TPC/AliTPCCalibPulser.cxx
- The default value of useFastDecoder in ctor changed to kFALSE,
[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 "AliTPCRawStreamFast.h"
185 #include "AliTPCCalROC.h"
186 #include "AliTPCCalPad.h"
187 #include "AliTPCROC.h"
188 #include "AliTPCParam.h"
189 #include "AliTPCCalibPulser.h"
190 #include "AliTPCcalibDB.h"
191 #include "AliMathBase.h"
192 #include "TTreeStream.h"
193
194 //date
195 #include "event.h"
196
197
198
199
200 ClassImp(AliTPCCalibPulser)
201
202 AliTPCCalibPulser::AliTPCCalibPulser() :
203     TObject(),
204     fFirstTimeBin(60),
205     fLastTimeBin(120),
206     fNbinsT0(200),
207     fXminT0(-2),
208     fXmaxT0(2),
209     fNbinsQ(200),
210     fXminQ(1),
211     fXmaxQ(40),
212     fNbinsRMS(100),
213     fXminRMS(0.1),
214     fXmaxRMS(5.1),
215     fLastSector(-1),
216     fOldRCUformat(kTRUE),
217     fROC(AliTPCROC::Instance()),
218     fParam(new AliTPCParam),
219     fPedestalTPC(0x0),
220     fPadNoiseTPC(0x0),
221     fPedestalROC(0x0),
222     fPadNoiseROC(0x0),
223     fCalRocArrayT0(72),
224     fCalRocArrayQ(72),
225     fCalRocArrayRMS(72),
226     fCalRocArrayOutliers(72),
227     fHistoQArray(72),
228     fHistoT0Array(72),
229     fHistoRMSArray(72),
230     fPadTimesArrayEvent(72),
231     fPadQArrayEvent(72),
232     fPadRMSArrayEvent(72),
233     fPadPedestalArrayEvent(72),
234     fCurrentChannel(-1),
235     fCurrentSector(-1),
236     fCurrentRow(-1),
237     fMaxPadSignal(-1),
238     fMaxTimeBin(-1),
239     fPadSignal(1024),
240     fPadPedestal(0),
241     fPadNoise(0),
242     fVTime0Offset(72),
243     fVTime0OffsetCounter(72),
244 //    fEvent(-1),
245     fDebugStreamer(0x0),
246     fDebugLevel(0)
247 {
248     //
249     // AliTPCSignal default constructor
250     //
251
252 }
253 //_____________________________________________________________________
254 AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
255     TObject(sig),
256     fFirstTimeBin(sig.fFirstTimeBin),
257     fLastTimeBin(sig.fLastTimeBin),
258     fNbinsT0(sig.fNbinsT0),
259     fXminT0(sig.fXminT0),
260     fXmaxT0(sig.fXmaxT0),
261     fNbinsQ(sig.fNbinsQ),
262     fXminQ(sig.fXminQ),
263     fXmaxQ(sig.fXmaxQ),
264     fNbinsRMS(sig.fNbinsRMS),
265     fXminRMS(sig.fXminRMS),
266     fXmaxRMS(sig.fXmaxRMS),
267     fLastSector(-1),
268     fOldRCUformat(kTRUE),
269     fROC(AliTPCROC::Instance()),
270     fParam(new AliTPCParam),
271     fPedestalTPC(0x0),
272     fPadNoiseTPC(0x0),
273     fPedestalROC(0x0),
274     fPadNoiseROC(0x0),
275     fCalRocArrayT0(72),
276     fCalRocArrayQ(72),
277     fCalRocArrayRMS(72),
278     fCalRocArrayOutliers(72),
279     fHistoQArray(72),
280     fHistoT0Array(72),
281     fHistoRMSArray(72),
282     fPadTimesArrayEvent(72),
283     fPadQArrayEvent(72),
284     fPadRMSArrayEvent(72),
285     fPadPedestalArrayEvent(72),
286     fCurrentChannel(-1),
287     fCurrentSector(-1),
288     fCurrentRow(-1),
289     fMaxPadSignal(-1),
290     fMaxTimeBin(-1),
291     fPadSignal(1024),
292     fPadPedestal(0),
293     fPadNoise(0),
294     fVTime0Offset(72),
295     fVTime0OffsetCounter(72),
296 //    fEvent(-1),
297     fDebugStreamer(0x0),
298     fDebugLevel(sig.fDebugLevel)
299 {
300     //
301     // AliTPCSignal default constructor
302     //
303
304     for (Int_t iSec = 0; iSec < 72; ++iSec){
305         const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
306         const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
307         const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
308         const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
309
310         const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
311         const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
312         const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
313
314         if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
315         if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
316         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
317         if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
318
319         if ( hQ != 0x0 ){
320             TH2S *hNew = new TH2S(*hQ);
321             hNew->SetDirectory(0);
322             fHistoQArray.AddAt(hNew,iSec);
323         }
324         if ( hT0 != 0x0 ){
325             TH2S *hNew = new TH2S(*hT0);
326             hNew->SetDirectory(0);
327             fHistoQArray.AddAt(hNew,iSec);
328         }
329         if ( hRMS != 0x0 ){
330             TH2S *hNew = new TH2S(*hRMS);
331             hNew->SetDirectory(0);
332             fHistoQArray.AddAt(hNew,iSec);
333         }
334     }
335
336 }
337 //_____________________________________________________________________
338 AliTPCCalibPulser& AliTPCCalibPulser::operator = (const  AliTPCCalibPulser &source)
339 {
340   //
341   // assignment operator
342   //
343   if (&source == this) return *this;
344   new (this) AliTPCCalibPulser(source);
345
346   return *this;
347 }
348 //_____________________________________________________________________
349 AliTPCCalibPulser::~AliTPCCalibPulser()
350 {
351     //
352     // destructor
353     //
354
355     fCalRocArrayT0.Delete();
356     fCalRocArrayQ.Delete();
357     fCalRocArrayRMS.Delete();
358     fCalRocArrayOutliers.Delete();
359
360     fHistoQArray.Delete();
361     fHistoT0Array.Delete();
362     fHistoRMSArray.Delete();
363
364     fPadTimesArrayEvent.Delete();
365     fPadQArrayEvent.Delete();
366     fPadRMSArrayEvent.Delete();
367     fPadPedestalArrayEvent.Delete();
368
369     if ( fDebugStreamer) delete fDebugStreamer;
370     delete fROC;
371     delete fParam;
372 }
373 //_____________________________________________________________________
374 Int_t AliTPCCalibPulser::Update(const Int_t icsector,
375                                 const Int_t icRow,
376                                 const Int_t icPad,
377                                 const Int_t icTimeBin,
378                                 const Float_t csignal)
379 {
380     //
381     // Signal filling methode on the fly pedestal and time offset correction if necessary.
382     // no extra analysis necessary. Assumes knowledge of the signal shape!
383     // assumes that it is looped over consecutive time bins of one pad
384     //
385     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
386
387     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
388
389     //init first pad and sector in this event
390     if ( fCurrentChannel == -1 ) {
391         fCurrentChannel = iChannel;
392         fCurrentSector  = icsector;
393         fCurrentRow     = icRow;
394     }
395
396     //process last pad if we change to a new one
397     if ( iChannel != fCurrentChannel ){
398         ProcessPad();
399         fCurrentChannel = iChannel;
400         fCurrentSector  = icsector;
401         fCurrentRow     = icRow;
402     }
403
404     //fill signals for current pad
405     fPadSignal[icTimeBin]=csignal;
406     if ( csignal > fMaxPadSignal ){
407         fMaxPadSignal = csignal;
408         fMaxTimeBin   = icTimeBin;
409     }
410     return 0;
411 }
412 //_____________________________________________________________________
413 void AliTPCCalibPulser::FindPedestal(Float_t part)
414 {
415     //
416     // find pedestal and noise for the current pad. Use either database or
417     // truncated mean with part*100%
418     //
419     Bool_t noPedestal = kTRUE;;
420     if (fPedestalTPC&&fPadNoiseTPC){
421         //use pedestal database
422         //only load new pedestals if the sector has changed
423         if ( fCurrentSector!=fLastSector ){
424             fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
425             fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
426             fLastSector=fCurrentSector;
427         }
428
429         if ( fPedestalROC&&fPadNoiseROC ){
430             fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
431             fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
432             noPedestal   = kFALSE;
433         }
434
435     }
436
437     //if we are not running with pedestal database, or for the current sector there is no information
438     //available, calculate the pedestal and noise on the fly
439     if ( noPedestal ) {
440         const Int_t kPedMax = 100;  //maximum pedestal value
441         Float_t  max    =  0;
442         Float_t  maxPos =  0;
443         Int_t    median =  -1;
444         Int_t    count0 =  0;
445         Int_t    count1 =  0;
446         //
447         Float_t padSignal=0;
448         //
449         UShort_t histo[kPedMax];
450         memset(histo,0,kPedMax*sizeof(UShort_t));
451
452         for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
453             padSignal = fPadSignal.GetMatrixArray()[i];
454             if (padSignal<=0) continue;
455             if (padSignal>max && i>10) {
456                 max = padSignal;
457                 maxPos = i;
458             }
459             if (padSignal>kPedMax-1) continue;
460             histo[Int_t(padSignal+0.5)]++;
461             count0++;
462         }
463             //
464         for (Int_t i=1; i<kPedMax; ++i){
465             if (count1<count0*0.5) median=i;
466             count1+=histo[i];
467         }
468         // truncated mean
469         //
470         Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
471         //
472         for (Int_t idelta=1; idelta<10; ++idelta){
473             if (median-idelta<=0) continue;
474             if (median+idelta>kPedMax) continue;
475             if (count<part*count1){
476                 count+=histo[median-idelta];
477                 mean +=histo[median-idelta]*(median-idelta);
478                 rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
479                 count+=histo[median+idelta];
480                 mean +=histo[median+idelta]*(median+idelta);
481                 rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
482             }
483         }
484         fPadPedestal = 0;
485         fPadNoise    = 0;
486         if ( count > 0 ) {
487             mean/=count;
488             rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
489             fPadPedestal = mean;
490             fPadNoise    = rms;
491         } 
492     }
493 }
494 //_____________________________________________________________________
495 void AliTPCCalibPulser::FindPulserSignal(TVectorD &param, Float_t &qSum)
496 {
497     //
498     //  Find position, signal width and height of the CE signal (last signal)
499     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
500     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
501     //
502
503     Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
504     Int_t   cemaxpos       = fMaxTimeBin;
505     Float_t ceSumThreshold = 8.*fPadNoise;  // threshold for the signal sum
506     const Int_t    kCemin  = 2;             // range for the analysis of the ce signal +- channels from the peak
507     const Int_t    kCemax  = 7;
508
509     if (cemaxpos!=0){
510         ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
511         for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
512             Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
513             if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
514                 ceTime+=signal*(i+0.5);
515                 ceRMS +=signal*(i+0.5)*(i+0.5);
516                 ceQsum+=signal;
517             }
518         }
519     }
520     if (ceQmax&&ceQsum>ceSumThreshold) {
521         ceTime/=ceQsum;
522         ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
523         fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
524         fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
525
526         //Normalise Q to the pad area
527         Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
528
529         ceQsum/=norm;
530     } else {
531         ceQmax=0;
532         ceTime=0;
533         ceRMS =0;
534         ceQsum=0;
535     }
536     param[0] = ceQmax;
537     param[1] = ceTime;
538     param[2] = ceRMS;
539     qSum     = ceQsum;
540 }
541 //_____________________________________________________________________
542 void AliTPCCalibPulser::ProcessPad()
543 {
544     //
545     //  Process data of current pad
546     //
547
548     FindPedestal();
549     TVectorD param(3);
550     Float_t  qSum;
551     FindPulserSignal(param, qSum);
552
553     Double_t meanT  = param[1];
554     Double_t sigmaT = param[2];
555
556     //Fill Event T0 counter
557     (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
558
559     //Fill Q histogram
560     GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
561
562     //Fill RMS histogram
563     GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
564
565
566     //Fill debugging info
567     if ( fDebugLevel>0 ){
568         (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
569         (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
570         (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
571     }
572     ResetPad();
573 }
574 //_____________________________________________________________________
575 void AliTPCCalibPulser::EndEvent()
576 {
577     //
578     //  Process data of current event
579     //
580
581     //check if last pad has allready been processed, if not do so
582     if ( fMaxTimeBin>-1 ) ProcessPad();
583
584     //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
585     for ( Int_t iSec = 0; iSec<72; ++iSec ){
586         TVectorF *vTimes = GetPadTimesEvent(iSec);
587         if ( !vTimes ) continue;
588
589         for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
590             Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
591             Float_t time  = (*vTimes).GetMatrixArray()[iChannel];
592
593             GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
594
595
596             //Debug start
597             if ( fDebugLevel>0 ){
598                 if ( !fDebugStreamer ) {
599                         //debug stream
600                     TDirectory *backup = gDirectory;
601                     fDebugStreamer = new TTreeSRedirector("deb2.root");
602                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
603                 }
604
605                 Int_t row=0;
606                 Int_t pad=0;
607                 Int_t padc=0;
608
609                 Float_t q   = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
610                 Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
611
612                 UInt_t channel=iChannel;
613                 Int_t sector=iSec;
614
615                 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
616                 pad = channel-fROC->GetRowIndexes(sector)[row];
617                 padc = pad-(fROC->GetNPads(sector,row)/2);
618
619                 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
620                                     Form("hSignalD%d.%d.%d",sector,row,pad),
621                                     fLastTimeBin-fFirstTimeBin,
622                                     fFirstTimeBin,fLastTimeBin);
623                 h1->SetDirectory(0);
624
625                 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
626                     h1->Fill(i,fPadSignal(i));
627
628                 (*fDebugStreamer) << "DataPad" <<
629 //                  "Event=" << fEvent <<
630                     "Sector="<< sector <<
631                     "Row="   << row<<
632                     "Pad="   << pad <<
633                     "PadC="  << padc <<
634                     "PadSec="<< channel <<
635                     "Time0="  << time0 <<
636                     "Time="  << time <<
637                     "RMS="   << rms <<
638                     "Sum="   << q <<
639                     "hist.=" << h1 <<
640                     "\n";
641
642                 delete h1;
643             }
644             //Debug end
645         }
646     }
647 }
648 //_____________________________________________________________________
649 Bool_t AliTPCCalibPulser::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
650 {
651   //
652   // Event Processing loop - AliTPCRawStream
653   //
654   ResetEvent();
655
656   Bool_t withInput = kFALSE;
657
658   while ( rawStreamFast->NextDDL() ){
659       while ( rawStreamFast->NextChannel() ){
660           Int_t isector  = rawStreamFast->GetSector();                       //  current sector
661           Int_t iRow     = rawStreamFast->GetRow();                          //  current row
662           Int_t iPad     = rawStreamFast->GetPad();                          //  current pad
663           Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
664           Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
665
666           while ( rawStreamFast->NextBunch() ){
667               for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
668                   Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
669                   Update(isector,iRow,iPad,iTimeBin+1,signal);
670                   withInput = kTRUE;
671               }
672           }
673       }
674   }
675   if (withInput){
676       EndEvent();
677   }
678   return withInput;
679 }
680 //_____________________________________________________________________
681 Bool_t AliTPCCalibPulser::ProcessEventFast(AliRawReader *rawReader)
682 {
683   //
684   //  Event processing loop - AliRawReader
685   //
686   AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader);
687   Bool_t res=ProcessEventFast(rawStreamFast);
688   delete rawStreamFast;
689   return res;
690 }
691 //_____________________________________________________________________
692 Bool_t AliTPCCalibPulser::ProcessEvent(AliTPCRawStream *rawStream)
693 {
694   //
695   // Event Processing loop - AliTPCRawStream
696   //
697
698   rawStream->SetOldRCUFormat(fOldRCUformat);
699
700   ResetEvent();
701
702   Bool_t withInput = kFALSE;
703
704   while (rawStream->Next()) {
705       Int_t isector  = rawStream->GetSector();                       //  current sector
706       Int_t iRow     = rawStream->GetRow();                          //  current row
707       Int_t iPad     = rawStream->GetPad();                          //  current pad
708       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
709       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
710
711       Update(isector,iRow,iPad,iTimeBin,signal);
712       withInput = kTRUE;
713   }
714   if (withInput){
715       EndEvent();
716   }
717   return withInput;
718 }
719 //_____________________________________________________________________
720 Bool_t AliTPCCalibPulser::ProcessEvent(AliRawReader *rawReader)
721 {
722   //
723   //  Event processing loop - AliRawReader
724   //
725
726
727   AliTPCRawStream rawStream(rawReader);
728
729   rawReader->Select("TPC");
730
731   return ProcessEvent(&rawStream);
732 }
733 //_____________________________________________________________________
734 Bool_t AliTPCCalibPulser::ProcessEvent(eventHeaderStruct *event)
735 {
736   //
737   //  Event processing loop - date event
738   //
739     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
740     Bool_t result=ProcessEvent(rawReader);
741     delete rawReader;
742     return result;
743
744 }
745 //_____________________________________________________________________
746 TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
747                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
748                                   Char_t *type, Bool_t force)
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     if ( !force || arr->UncheckedAt(sector) )
755         return (TH2S*)arr->UncheckedAt(sector);
756
757     // if we are forced and histogram doesn't yes exist create it
758     Char_t name[255], title[255];
759
760     sprintf(name,"hCalib%s%.2d",type,sector);
761     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
762
763     // new histogram with Q calib information. One value for each pad!
764     TH2S* hist = new TH2S(name,title,
765                           nbinsY, ymin, ymax,
766                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
767     hist->SetDirectory(0);
768     arr->AddAt(hist,sector);
769     return hist;
770 }
771 //_____________________________________________________________________
772 TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
773 {
774     //
775     // return pointer to T0 histogram
776     // if force is true create a new histogram if it doesn't exist allready
777     //
778     TObjArray *arr = &fHistoT0Array;
779     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
780 }
781 //_____________________________________________________________________
782 TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
783 {
784     //
785     // return pointer to Q histogram
786     // if force is true create a new histogram if it doesn't exist allready
787     //
788     TObjArray *arr = &fHistoQArray;
789     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
790 }
791 //_____________________________________________________________________
792 TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
793 {
794     //
795     // return pointer to Q histogram
796     // if force is true create a new histogram if it doesn't exist allready
797     //
798     TObjArray *arr = &fHistoRMSArray;
799     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
800 }
801 //_____________________________________________________________________
802 TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
803 {
804     //
805     // return pointer to Pad Info from 'arr' for the current event and sector
806     // if force is true create it if it doesn't exist allready
807     //
808     if ( !force || arr->UncheckedAt(sector) )
809         return (TVectorF*)arr->UncheckedAt(sector);
810
811     TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
812     arr->AddAt(vect,sector);
813     return vect;
814 }
815 //_____________________________________________________________________
816 TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
817 {
818     //
819     // return pointer to Pad Times Array for the current event and sector
820     // if force is true create it if it doesn't exist allready
821     //
822     TObjArray *arr = &fPadTimesArrayEvent;
823     return GetPadInfoEvent(sector,arr,force);
824 }
825 //_____________________________________________________________________
826 TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
827 {
828     //
829     // return pointer to Pad Q Array for the current event and sector
830     // if force is true create it if it doesn't exist allready
831     // for debugging purposes only
832     //
833
834     TObjArray *arr = &fPadQArrayEvent;
835     return GetPadInfoEvent(sector,arr,force);
836 }
837 //_____________________________________________________________________
838 TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
839 {
840     //
841     // return pointer to Pad RMS Array for the current event and sector
842     // if force is true create it if it doesn't exist allready
843     // for debugging purposes only
844     //
845     TObjArray *arr = &fPadRMSArrayEvent;
846     return GetPadInfoEvent(sector,arr,force);
847 }
848 //_____________________________________________________________________
849 TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
850 {
851     //
852     // return pointer to Pad RMS Array for the current event and sector
853     // if force is true create it if it doesn't exist allready
854     // for debugging purposes only
855     //
856     TObjArray *arr = &fPadPedestalArrayEvent;
857     return GetPadInfoEvent(sector,arr,force);
858 }
859 //_____________________________________________________________________
860 AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
861 {
862     //
863     // return pointer to ROC Calibration
864     // if force is true create a new histogram if it doesn't exist allready
865     //
866     if ( !force || arr->UncheckedAt(sector) )
867         return (AliTPCCalROC*)arr->UncheckedAt(sector);
868
869     // if we are forced and histogram doesn't yes exist create it
870
871     // new AliTPCCalROC for T0 information. One value for each pad!
872     AliTPCCalROC *croc = new AliTPCCalROC(sector);
873     arr->AddAt(croc,sector);
874     return croc;
875 }
876 //_____________________________________________________________________
877 AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
878 {
879     //
880     // return pointer to Carge ROC Calibration
881     // if force is true create a new histogram if it doesn't exist allready
882     //
883     TObjArray *arr = &fCalRocArrayT0;
884     return GetCalRoc(sector, arr, force);
885 }
886 //_____________________________________________________________________
887 AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
888 {
889     //
890     // return pointer to T0 ROC Calibration
891     // if force is true create a new histogram if it doesn't exist allready
892     //
893     TObjArray *arr = &fCalRocArrayQ;
894     return GetCalRoc(sector, arr, force);
895 }
896 //_____________________________________________________________________
897 AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
898 {
899     //
900     // return pointer to signal width ROC Calibration
901     // if force is true create a new histogram if it doesn't exist allready
902     //
903     TObjArray *arr = &fCalRocArrayRMS;
904     return GetCalRoc(sector, arr, force);
905 }
906 //_____________________________________________________________________
907 AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
908 {
909     //
910     // return pointer to Outliers
911     // if force is true create a new histogram if it doesn't exist allready
912     //
913     TObjArray *arr = &fCalRocArrayOutliers;
914     return GetCalRoc(sector, arr, force);
915 }
916 //_____________________________________________________________________
917 void AliTPCCalibPulser::ResetEvent()
918 {
919     //
920     //  Reset global counters  -- Should be called before each event is processed
921     //
922     fLastSector=-1;
923     fCurrentSector=-1;
924     fCurrentRow=-1;
925     fCurrentChannel=-1;
926
927     ResetPad();
928
929     fPadTimesArrayEvent.Delete();
930     fPadQArrayEvent.Delete();
931     fPadRMSArrayEvent.Delete();
932     fPadPedestalArrayEvent.Delete();
933
934     for ( Int_t i=0; i<72; ++i ){
935         fVTime0Offset[i]=0;
936         fVTime0OffsetCounter[i]=0;
937     }
938 }
939 //_____________________________________________________________________
940 void AliTPCCalibPulser::ResetPad()
941 {
942     //
943     //  Reset pad infos -- Should be called after a pad has been processed
944     //
945     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
946         fPadSignal[i] = 0;
947     fMaxTimeBin = -1;
948     fMaxPadSignal = -1;
949     fPadPedestal  = -1;
950     fPadNoise     = -1;
951 }
952 //_____________________________________________________________________
953 void AliTPCCalibPulser::Merge(AliTPCCalibPulser *sig)
954 {
955     //
956     //  Merge reference histograms of sig to the current AliTPCCalibPulser
957     //
958
959     //merge histograms
960     for (Int_t iSec=0; iSec<72; ++iSec){
961         TH2S *hRefQmerge   = sig->GetHistoQ(iSec);
962         TH2S *hRefT0merge  = sig->GetHistoT0(iSec);
963         TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
964
965
966         if ( hRefQmerge ){
967             TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
968             TH2S *hRefQ   = GetHistoQ(iSec);
969             if ( hRefQ ) hRefQ->Add(hRefQmerge);
970             else {
971                 TH2S *hist = new TH2S(*hRefQmerge);
972                 hist->SetDirectory(0);
973                 fHistoQArray.AddAt(hist, iSec);
974             }
975             hRefQmerge->SetDirectory(dir);
976         }
977         if ( hRefT0merge ){
978             TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
979             TH2S *hRefT0  = GetHistoT0(iSec);
980             if ( hRefT0 ) hRefT0->Add(hRefT0merge);
981             else {
982                 TH2S *hist = new TH2S(*hRefT0merge);
983                 hist->SetDirectory(0);
984                 fHistoT0Array.AddAt(hist, iSec);
985             }
986             hRefT0merge->SetDirectory(dir);
987         }
988         if ( hRefRMSmerge ){
989             TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
990             TH2S *hRefRMS = GetHistoRMS(iSec);
991             if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
992             else {
993                 TH2S *hist = new TH2S(*hRefRMSmerge);
994                 hist->SetDirectory(0);
995                 fHistoRMSArray.AddAt(hist, iSec);
996             }
997             hRefRMSmerge->SetDirectory(dir);
998         }
999
1000     }
1001 }
1002 //_____________________________________________________________________
1003 void AliTPCCalibPulser::Analyse()
1004 {
1005     //
1006     //  Calculate calibration constants
1007     //
1008
1009     TVectorD paramQ(3);
1010     TVectorD paramT0(3);
1011     TVectorD paramRMS(3);
1012     TMatrixD dummy(3,3);
1013
1014     for (Int_t iSec=0; iSec<72; ++iSec){
1015         TH2S *hT0 = GetHistoT0(iSec);
1016         if (!hT0 ) continue;
1017
1018         AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
1019         AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
1020         AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1021         AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1022
1023         TH2S *hQ   = GetHistoQ(iSec);
1024         TH2S *hRMS = GetHistoRMS(iSec);
1025
1026         Short_t *arrayhQ   = hQ->GetArray();
1027         Short_t *arrayhT0  = hT0->GetArray();
1028         Short_t *arrayhRMS = hRMS->GetArray();
1029
1030         UInt_t nChannels = fROC->GetNChannels(iSec);
1031
1032         //debug
1033         Int_t row=0;
1034         Int_t pad=0;
1035         Int_t padc=0;
1036         //! debug
1037
1038         for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1039
1040             Float_t cogTime0 = -1000;
1041             Float_t cogQ     = -1000;
1042             Float_t cogRMS   = -1000;
1043             Float_t cogOut   = 0;
1044
1045             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1046             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1047             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1048 /*
1049             AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1050             AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1051             AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1052             cogQ     = paramQ[1];
1053             cogTime0 = paramT0[1];
1054             cogRMS   = paramRMS[1];
1055 */
1056             cogQ     = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1057             cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1058             cogRMS   = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1059
1060             /*
1061             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1062                 cogOut = 1;
1063                 cogTime0 = 0;
1064                 cogQ     = 0;
1065                 cogRMS   = 0;
1066             }
1067 */
1068             rocQ->SetValue(iChannel, cogQ*cogQ);
1069             rocT0->SetValue(iChannel, cogTime0);
1070             rocRMS->SetValue(iChannel, cogRMS);
1071             rocOut->SetValue(iChannel, cogOut);
1072
1073
1074             //debug
1075             if ( fDebugLevel > 0 ){
1076                 if ( !fDebugStreamer ) {
1077                         //debug stream
1078                     TDirectory *backup = gDirectory;
1079                     fDebugStreamer = new TTreeSRedirector("deb2.root");
1080                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1081                 }
1082
1083                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1084                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1085                 padc = pad-(fROC->GetNPads(iSec,row)/2);
1086
1087                 (*fDebugStreamer) << "DataEnd" <<
1088                     "Sector="  << iSec      <<
1089                     "Pad="     << pad       <<
1090                     "PadC="    << padc      <<
1091                     "Row="     << row       <<
1092                     "PadSec="  << iChannel   <<
1093                     "Q="       << cogQ      <<
1094                     "T0="      << cogTime0  <<
1095                     "RMS="     << cogRMS    <<
1096                     "\n";
1097             }
1098             //! debug
1099         }
1100     }
1101     delete fDebugStreamer;
1102     fDebugStreamer = 0x0;
1103 }
1104 //_____________________________________________________________________
1105 void AliTPCCalibPulser::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1106 {
1107     //
1108     //  Write class to file
1109     //
1110
1111     TString sDir(dir);
1112     TString option;
1113
1114     if ( append )
1115         option = "update";
1116     else
1117         option = "recreate";
1118
1119     TDirectory *backup = gDirectory;
1120     TFile f(filename,option.Data());
1121     f.cd();
1122     if ( !sDir.IsNull() ){
1123         f.mkdir(sDir.Data());
1124         f.cd(sDir);
1125     }
1126     this->Write();
1127     f.Close();
1128
1129     if ( backup ) backup->cd();
1130 }
1131 //_____________________________________________________________________
1132 //_________________________  Test Functions ___________________________
1133 //_____________________________________________________________________
1134 TObjArray* AliTPCCalibPulser::TestBinning()
1135 {
1136     //
1137     //  Function to test the binning of the reference histograms
1138     //  type: T0, Q or RMS
1139     //  mode: 0 - number of filled bins per channel
1140     //        1 - number of empty bins between filled bins in one ROC
1141     //  returns TObjArray with the test histograms type*2+mode:
1142     //  position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1143
1144
1145     TObjArray *histArray = new TObjArray(6);
1146     const Char_t *type[] = {"T0","Q","RMS"};
1147     Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1148
1149     for (Int_t itype = 0; itype<3; ++itype){
1150         for (Int_t imode=0; imode<2; ++imode){
1151             Int_t icount = itype*2+imode;
1152             histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1153                                       Form("Test Binning of '%s', mode - %d",type[itype],imode),
1154                                       72,0,72),
1155                              icount);
1156         }
1157     }
1158
1159
1160     TH2S *hRef=0x0;
1161     Short_t *array=0x0;
1162     for (Int_t itype = 0; itype<3; ++itype){
1163         for (Int_t iSec=0; iSec<72; ++iSec){
1164             if ( itype == 0 ) hRef = GetHistoT0(iSec);
1165             if ( itype == 1 ) hRef = GetHistoQ(iSec);
1166             if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1167             if ( hRef == 0x0 ) continue;
1168             array = (hRef->GetArray());
1169             UInt_t nChannels = fROC->GetNChannels(iSec);
1170
1171             Int_t nempty=0;
1172             for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1173                 Int_t nfilled=0;
1174                 Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1175                 Int_t c1 = 0;
1176                 Int_t c2 = 0;
1177                 for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1178                     if ( array[offset+iBin]>0 ) {
1179                         nfilled++;
1180                         if ( c1 && c2 ) nempty++;
1181                         else c1 = 1;
1182                     }
1183                     else if ( c1 ) c2 = 1;
1184
1185                 }
1186                 ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1187             }
1188             ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);
1189         }
1190     }
1191     return histArray;
1192 }