Ading 2 new types of analysis (Jens)
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibPedestal.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 /* $Id$ */
18
19
20 //Root includes
21 #include <TH1F.h>
22 #include <TH2F.h>
23 #include <TString.h>
24 #include <TMath.h>
25 #include <TF1.h>
26 #include <TRandom.h>
27 #include <TDirectory.h>
28 #include <TFile.h>
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 "AliTPCROC.h"
36 #include "AliMathBase.h"
37 #include "TTreeStream.h"
38 #include "AliTPCRawStreamFast.h"
39
40 //date
41 #include "event.h"
42
43 //header file
44 #include "AliTPCCalibPedestal.h"
45
46
47 ///////////////////////////////////////////////////////////////////////////////////////
48 //          Implementation of the TPC pedestal and noise 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 // Working principle:
58 // ------------------
59 // Raw pedestal data is processed by calling one of the ProcessEvent(...) functions
60 // (see below). These in the end call the Update(...) function, where the data is filled
61 // into histograms.
62 //
63 // For each ROC one TH2F histo (ROC channel vs. ADC channel) is created when
64 // it is filled for the first time (GetHistoPedestal(ROC,kTRUE)). All histos are stored in the
65 // TObjArray fHistoPedestalArray.
66 //
67 // For a fast filling of the histogram the corresponding bin number of the channel and ADC channel
68 // is computed by hand and the histogram array is accessed directly via its pointer.
69 // ATTENTION: Doing so the the entry counter of the histogram is not increased
70 //            this means that e.g. the colz draw option gives an empty plot unless
71 //          calling 'histo->SetEntries(1)' before drawing.
72 //
73 // After accumulating the desired statistics the Analyse() function has to be called.
74 // Whithin this function the pedestal and noise values are calculated for each pad, using
75 // the fast gaus fit function  AliMathBase::FitGaus(...), and the calibration
76 // storage classes (AliTPCCalROC) are filled for each ROC.
77 // The calibration information is stored in the TObjArrays fCalRocArrayPedestal and fCalRocArrayRMS;
78 //
79 //
80 //
81 // User interface for filling data:
82 // --------------------------------
83 //
84 // To Fill information one of the following functions can be used:
85 //
86 // Bool_t ProcessEvent(eventHeaderStruct *event);
87 //   - process Date event
88 //   - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
89 //
90 // Bool_t ProcessEvent(AliRawReader *rawReader);
91 //  - process AliRawReader event
92 //   - use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)
93 //
94 // Bool_t ProcessEvent(AliTPCRawStream *rawStream);
95 //   - process event from AliTPCRawStream
96 //   - call Update function for signal filling
97 //
98 // Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
99 //              iPad,  const Int_t iTimeBin, const Float_t signal);
100 //   - directly  fill signal information (sector, row, pad, time bin, pad)
101 //     to the reference histograms
102 //
103 // It is also possible to merge two independently taken calibrations using the function
104 //
105 // void Merge(AliTPCCalibPedestal *ped)
106 //   - copy histograms in 'ped' if the do not exist in this instance
107 //   - Add histograms in 'ped' to the histograms in this instance if the allready exist
108 //   - After merging call Analyse again!
109 //
110 //
111 //
112 // -- example: filling data using root raw data:
113 // void fillPedestal(Char_t *filename)
114 // {
115 //    rawReader = new AliRawReaderRoot(fileName);
116 //    if ( !rawReader ) return;
117 //    AliTPCCalibPedestal *calib = new AliTPCCalibPedestal;
118 //    while (rawReader->NextEvent()){
119 //      calib->ProcessEvent(rawReader);
120 //    }
121 //    calib->Analyse();
122 //    calib->DumpToFile("PedestalData.root");
123 //    delete rawReader;
124 //    delete calib;
125 // }
126 //
127 //
128 // What kind of information is stored and how to retrieve them:
129 // ------------------------------------------------------------
130 //
131 // - Accessing the 'Reference Histograms' (pedestal distribution histograms):
132 //
133 // TH2F *GetHistoPedestal(Int_t sector);
134 //
135 // - Accessing the calibration storage objects:
136 //
137 // AliTPCCalROC *GetCalRocPedestal(Int_t sector);  - for the pedestal values, mean from gaus fit
138 // AliTPCCalROC *GetCalRocSigma(Int_t sector);     - for the Noise values, sigma from guas fit
139 // AliTPCCalROC *GetCalRocMean(Int_t sector);  - for the pedestal values, truncated mean
140 // AliTPCCalROC *GetCalRocRMS(Int_t sector);     - for the Noise values, rms from truncated mean
141 //
142 // example for visualisation:
143 // if the file "PedestalData.root" was created using the above example one could do the following:
144 //
145 // TFile filePedestal("PedestalData.root")
146 // AliTPCCalibPedestal *ped = (AliTPCCalibPedestal*)filePedestal->Get("AliTPCCalibPedestal");
147 // ped->GetCalRocPedestal(0)->Draw("colz");
148 // ped->GetCalRocRMS(0)->Draw("colz");
149 //
150 // or use the AliTPCCalPad functionality:
151 // AliTPCCalPad padPedestal(ped->GetCalPadPedestal());
152 // AliTPCCalPad padNoise(ped->GetCalPadRMS());
153 // padPedestal->MakeHisto2D()->Draw("colz");  //Draw A-Side Pedestal Information
154 // padNoise->MakeHisto2D()->Draw("colz");  //Draw A-Side Noise Information
155 //
156 /*
157  example: fill pedestal with gausschen noise
158  AliTPCCalibPedestal ped;
159  ped.TestEvent();
160  ped.Analyse();
161  //Draw output;
162  TCanvas* c1 = new TCanvas;
163  c1->Divide(1,2);
164  c1->cd(1);
165  ped.GetHistoPedestal(0)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
166  ped.GetHistoPedestal(0)->Draw("colz");
167  c1->cd(2);
168  ped.GetHistoPedestal(36)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
169  ped.GetHistoPedestal(36)->Draw("colz");
170  TCanvas* c2 = new TCanvas;
171  c2->Divide(2,2);
172  c2->cd(1);
173  ped.GetCalRocPedestal(0)->Draw("colz");
174  c2->cd(2);
175  ped.GetCalRocRMS(0)->Draw("colz");
176  c2->cd(3);
177  ped.GetCalRocPedestal(36)->Draw("colz");
178  c2->cd(4);
179  ped.GetCalRocRMS(36)->Draw("colz");
180 */
181 //
182 // Time dependent pedestals:
183 //
184 // If wished there is the possibility to calculate for each channel and time bin
185 // the mean pedestal [pedestals(t)]. This is done by
186 //
187 // 1) setting SetTimeAnalysis(kTRUE),
188 // 2) processing the data by looping over the events using ProcessEvent(..)
189 // 3) calling the Analyse() and AnalyseTime(nevents) functions (providing nevents)
190 // 4) getting the pedestals(t) using   TArrayF **timePed = calibPedestal.GetTimePedestals();
191 // 5) looking at values using   timePed[row][pad].At(timebin)
192 //
193 // This functionality is intended to be used on an LDC bu the detector algorithm
194 // (TPCPEDESTALda) to generate a data set used for configuration of the pattern
195 // memory for baseline subtraction in the ALTROs. Later the information should also
196 // be stored as reference data.
197 //
198
199
200 ClassImp(AliTPCCalibPedestal)
201
202 AliTPCCalibPedestal::AliTPCCalibPedestal() : 
203   TObject(),
204   fFirstTimeBin(60),
205   fLastTimeBin(1000),
206   fAdcMin(1),
207   fAdcMax(100),
208   fAnaMeanDown(0.),
209   fAnaMeanUp(1.),
210   fOldRCUformat(kTRUE),
211   fTimeAnalysis(kFALSE),
212   fROC(AliTPCROC::Instance()),
213   fMapping(NULL),
214   fCalRocArrayPedestal(72),
215   fCalRocArraySigma(72),
216   fHistoPedestalArray(72),
217   fTimeSignal(NULL),
218   fCalRocArrayMean(72),
219   fCalRocArrayRMS(72)
220 {
221   //
222   // default constructor
223   //
224 }
225
226
227 //_____________________________________________________________________
228 AliTPCCalibPedestal::AliTPCCalibPedestal(const AliTPCCalibPedestal &ped) : 
229   TObject(ped),
230   fFirstTimeBin(ped.GetFirstTimeBin()),
231   fLastTimeBin(ped.GetLastTimeBin()),
232   fAdcMin(ped.GetAdcMin()),
233   fAdcMax(ped.GetAdcMax()),
234   fAnaMeanDown(ped.fAnaMeanDown),
235   fAnaMeanUp(ped.fAnaMeanUp),
236   fOldRCUformat(ped.fOldRCUformat),
237   fTimeAnalysis(ped.fTimeAnalysis),
238   fROC(AliTPCROC::Instance()),
239   fMapping(NULL),
240   fCalRocArrayPedestal(72),
241   fCalRocArraySigma(72),
242   fHistoPedestalArray(72),
243   fTimeSignal(ped.fTimeSignal),
244   fCalRocArrayMean(72),
245   fCalRocArrayRMS(72)
246 {
247   //
248   // copy constructor
249   //
250   for (Int_t iSec = 0; iSec < 72; ++iSec){
251     const AliTPCCalROC *calPed = (AliTPCCalROC*)ped.fCalRocArrayPedestal.UncheckedAt(iSec);
252     const AliTPCCalROC *calRMS = (AliTPCCalROC*)ped.fCalRocArrayRMS.UncheckedAt(iSec);
253     const TH2F         *hPed   = (TH2F*)ped.fHistoPedestalArray.UncheckedAt(iSec);
254     
255     if ( calPed != 0x0 ) fCalRocArrayPedestal.AddAt(new AliTPCCalROC(*calPed), iSec);
256     if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
257     
258     if ( hPed != 0x0 ){
259       TH2F *hNew = new TH2F(*hPed);
260       hNew->SetDirectory(0);
261       fHistoPedestalArray.AddAt(hNew,iSec);
262     }
263   }
264 }
265
266
267 //_____________________________________________________________________
268 AliTPCCalibPedestal& AliTPCCalibPedestal::operator = (const  AliTPCCalibPedestal &source)
269 {
270   //
271   // assignment operator
272   //
273   if (&source == this) return *this;
274   new (this) AliTPCCalibPedestal(source);
275
276   return *this;
277 }
278
279
280 //_____________________________________________________________________
281 AliTPCCalibPedestal::~AliTPCCalibPedestal() 
282 {
283   //
284   // destructor
285   //
286
287   fCalRocArrayPedestal.Delete();
288   fCalRocArrayRMS.Delete();
289   fHistoPedestalArray.Delete();
290
291   if ( fTimeSignal ) {
292     for (Int_t i = 0; i < 159; i++) {
293       delete [] fTimeSignal[i];
294       fTimeSignal[i] = 0;
295     }
296     delete [] fTimeSignal;
297     fTimeSignal = 0;
298   }
299
300   // do not delete fMapping, because we do not own it.
301
302 }
303
304
305 //_____________________________________________________________________
306 void AliTPCCalibPedestal::SetTimeAnalysis(Bool_t time)
307 {
308   //
309   // Use time dependent analysis: Pedestals are analysed as a function
310   // of the drift time. There is one mean value generated for each time
311   // bin and each channel. It can be used as reference data and for
312   // configuration of the ALTRO pattern memory for baseline subtraction.
313   //
314   // ATTENTION: Use only on LDC in TPCPEDESTALda! On a LDC we get data
315   // only from one sector. For the full TPC we would need a lot of
316   // memory (36*159*140*1024*4bytes = 3.3GB)!
317   //
318
319   fTimeAnalysis = time;
320
321   if ( !fTimeAnalysis ) return;
322
323   // prepare array for one sector (159*140*1024*4bytes = 92MB):
324   fTimeSignal = new TArrayF*[159];
325   for (Int_t i = 0; i < 159; i++) {  // padrows
326     fTimeSignal[i] = new TArrayF[140];
327     for (Int_t j = 0; j < 140; j++) {  // pads per row
328       fTimeSignal[i][j].Set(1024);
329       for (Int_t k = 0; k < 1024; k++) {  // time bins per pad
330         fTimeSignal[i][j].AddAt(0., k);
331       }
332     }      
333   }
334 }
335
336
337 //_____________________________________________________________________
338 Int_t AliTPCCalibPedestal::Update(const Int_t icsector, 
339                                   const Int_t icRow,
340                                   const Int_t icPad,
341                                   const Int_t icTimeBin,
342                                   const Float_t csignal)
343 {
344   //
345   // Signal filling method
346   //
347   if (icRow<0) return 0;
348   if (icPad<0) return 0;
349   if (icTimeBin<0) return 0;
350  
351   // Time dependent pedestals
352   if ( fTimeAnalysis ) {
353     if ( icsector < 36 ) // IROC
354       fTimeSignal[icRow][icPad].AddAt(fTimeSignal[icRow][icPad].At(icTimeBin)+csignal, icTimeBin);
355     else 
356       fTimeSignal[icRow+63][icPad].AddAt(fTimeSignal[icRow+63][icPad].At(icTimeBin)+csignal, icTimeBin);
357   }
358   //return if we are out of the specified time bin or adc range
359   if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
360   if ( ((Int_t)csignal>fAdcMax) || ((Int_t)csignal<fAdcMin)  ) return 0;
361
362   Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
363
364   // fast filling method
365   // Attention: the entry counter of the histogram is not increased
366   //            this means that e.g. the colz draw option gives an empty plot
367   Int_t bin = (iChannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
368
369   GetHistoPedestal(icsector,kTRUE)->GetArray()[bin]++;
370
371   return 0;
372 }
373
374
375 //_____________________________________________________________________
376 Bool_t AliTPCCalibPedestal::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
377 {
378   //
379   // Event Processing loop - AliTPCRawStream
380   //
381   Bool_t withInput = kFALSE;
382
383   while ( rawStreamFast->NextDDL() ){
384       while ( rawStreamFast->NextChannel() ){
385           Int_t isector  = rawStreamFast->GetSector();                       //  current sector
386           Int_t iRow     = rawStreamFast->GetRow();                          //  current row
387           Int_t iPad     = rawStreamFast->GetPad();                          //  current pad
388
389           while ( rawStreamFast->NextBunch() ){
390             Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
391             Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
392             for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
393                   Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
394                   Update(isector,iRow,iPad,iTimeBin+1,signal);
395                   withInput = kTRUE;
396               }
397           }
398       }
399   }
400
401   return withInput;
402 }
403 //_____________________________________________________________________
404 Bool_t AliTPCCalibPedestal::ProcessEventFast(AliRawReader *rawReader)
405 {
406   //
407   //  Event processing loop - AliRawReader
408   //
409   AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
410   Bool_t res=ProcessEventFast(rawStreamFast);
411   delete rawStreamFast;
412   return res;
413 }
414
415 //_____________________________________________________________________
416 Bool_t AliTPCCalibPedestal::ProcessEvent(AliTPCRawStream *rawStream)
417 {
418   //
419   // Event Processing loop - AliTPCRawStream
420   //
421
422   rawStream->SetOldRCUFormat(fOldRCUformat);
423
424   Bool_t withInput = kFALSE;
425
426   while (rawStream->Next()) {
427
428     Int_t iSector  = rawStream->GetSector();      //  current ROC
429     Int_t iRow     = rawStream->GetRow();         //  current row
430     Int_t iPad     = rawStream->GetPad();         //  current pad
431     Int_t iTimeBin = rawStream->GetTime();        //  current time bin
432     Float_t signal = rawStream->GetSignal();      //  current ADC signal
433     
434     Update(iSector,iRow,iPad,iTimeBin,signal);
435     withInput = kTRUE;
436   }
437
438   return withInput;
439 }
440
441
442 //_____________________________________________________________________
443 Bool_t AliTPCCalibPedestal::ProcessEvent(AliRawReader *rawReader)
444 {
445   //
446   //  Event processing loop - AliRawReader
447   //
448
449   // if fMapping is NULL the rawstream will crate its own mapping
450   AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
451   rawReader->Select("TPC");
452   return ProcessEvent(&rawStream);
453 }
454
455
456 //_____________________________________________________________________
457 Bool_t AliTPCCalibPedestal::ProcessEvent(eventHeaderStruct *event)
458 {
459   //
460   //  process date event
461   //
462
463   AliRawReader *rawReader = new AliRawReaderDate((void*)event);
464   Bool_t result=ProcessEvent(rawReader);
465   delete rawReader;
466   return result;
467 }
468
469
470 //_____________________________________________________________________
471 Bool_t AliTPCCalibPedestal::TestEvent() 
472 {
473   //
474   //  Test event loop
475   // fill one oroc and one iroc with random gaus
476   //
477
478     gRandom->SetSeed(0);
479
480     for (UInt_t iSec=0; iSec<72; ++iSec){
481         if (iSec%36>0) continue;
482         for (UInt_t iRow=0; iRow < fROC->GetNRows(iSec); ++iRow){
483             for (UInt_t iPad=0; iPad < fROC->GetNPads(iSec,iRow); ++iPad){
484                 for (UInt_t iTimeBin=0; iTimeBin<1024; ++iTimeBin){
485                     Float_t signal=(Int_t)(iRow+3+gRandom->Gaus(0,.7));
486                     if ( signal>0 )Update(iSec,iRow,iPad,iTimeBin,signal);
487                 }
488             }
489         }
490     }
491     return kTRUE;
492 }
493
494
495 //_____________________________________________________________________
496 TH2F* AliTPCCalibPedestal::GetHisto(Int_t sector, TObjArray *arr, 
497                                     Int_t nbinsY, Float_t ymin, Float_t ymax,
498                                     Char_t *type, Bool_t force)
499 {
500     //
501     // return pointer to Q histogram
502     // if force is true create a new histogram if it doesn't exist allready
503     //
504     if ( !force || arr->UncheckedAt(sector) )
505         return (TH2F*)arr->UncheckedAt(sector);
506
507     // if we are forced and histogram doesn't yes exist create it
508     Char_t name[255], title[255];
509
510     sprintf(name,"hCalib%s%.2d",type,sector);
511     sprintf(title,"%s calibration histogram sector %.2d;ADC channel;Channel (pad)",type,sector);
512
513     // new histogram with Q calib information. One value for each pad!
514     TH2F* hist = new TH2F(name,title,
515                           nbinsY, ymin, ymax,
516                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector)
517                          );
518     hist->SetDirectory(0);
519     arr->AddAt(hist,sector);
520     return hist;
521 }
522
523
524 //_____________________________________________________________________
525 TH2F* AliTPCCalibPedestal::GetHistoPedestal(Int_t sector, Bool_t force) 
526 {
527     //
528     // return pointer to T0 histogram
529     // if force is true create a new histogram if it doesn't exist allready
530     //
531     TObjArray *arr = &fHistoPedestalArray;
532     return GetHisto(sector, arr, fAdcMax-fAdcMin, fAdcMin, fAdcMax, "Pedestal", force);
533 }
534
535
536 //_____________________________________________________________________
537 AliTPCCalROC* AliTPCCalibPedestal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) 
538 {
539     //
540     // return pointer to ROC Calibration
541     // if force is true create a new histogram if it doesn't exist allready
542     //
543     if ( !force || arr->UncheckedAt(sector) )
544         return (AliTPCCalROC*)arr->UncheckedAt(sector);
545
546     // if we are forced and the histogram doesn't yet exist create it
547
548     // new AliTPCCalROC for T0 information. One value for each pad!
549     AliTPCCalROC *croc = new AliTPCCalROC(sector);
550     arr->AddAt(croc,sector);
551     return croc;
552 }
553
554
555 //_____________________________________________________________________
556 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocPedestal(Int_t sector, Bool_t force) 
557 {
558     //
559     // return pointer to ROC with Pedestal data
560     // if force is true create a new histogram if it doesn't exist allready
561     //
562     TObjArray *arr = &fCalRocArrayPedestal;
563     return GetCalRoc(sector, arr, force);
564 }
565
566
567 //_____________________________________________________________________
568 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocSigma(Int_t sector, Bool_t force) 
569 {
570     //
571     // return pointer to  ROC with signal witdth in sigma
572     // if force is true create a new histogram if it doesn't exist allready
573     //
574     TObjArray *arr = &fCalRocArraySigma;
575     return GetCalRoc(sector, arr, force);
576 }
577 //_____________________________________________________________________
578 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocMean(Int_t sector, Bool_t force)
579 {
580   //
581     // return pointer to ROC with signal mean information
582     // if force is true create a new histogram if it doesn't exist allready
583   //
584   TObjArray *arr = &fCalRocArrayMean;
585   return GetCalRoc(sector, arr, force);
586 }
587
588 //_____________________________________________________________________
589 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocRMS(Int_t sector, Bool_t force) 
590 {
591   //
592     // return pointer to signal width ROC Calibration
593     // if force is true create a new histogram if it doesn't exist allready
594   //
595   TObjArray *arr = &fCalRocArrayRMS;
596   return GetCalRoc(sector, arr, force);
597 }
598
599
600 //_____________________________________________________________________
601 void AliTPCCalibPedestal::Merge(AliTPCCalibPedestal *ped)
602 {
603   //
604   //  Merge reference histograms of sig to the current AliTPCCalibSignal
605   //
606
607   // merge histograms
608   for (Int_t iSec=0; iSec<72; ++iSec){
609     TH2F *hRefPedMerge   = ped->GetHistoPedestal(iSec);
610
611     if ( hRefPedMerge ){
612       TDirectory *dir = hRefPedMerge->GetDirectory(); hRefPedMerge->SetDirectory(0);
613       TH2F *hRefPed   = GetHistoPedestal(iSec);
614       if ( hRefPed ) hRefPed->Add(hRefPedMerge);
615       else {
616         TH2F *hist = new TH2F(*hRefPedMerge);
617         hist->SetDirectory(0);
618         fHistoPedestalArray.AddAt(hist, iSec);
619       }
620       hRefPedMerge->SetDirectory(dir);
621     }
622   }
623
624   // merge array
625   // ...
626
627 }
628
629
630 //_____________________________________________________________________
631 void AliTPCCalibPedestal::Analyse() 
632 {
633   //
634   //  Calculate calibration constants
635   //
636
637   Int_t nbinsAdc = fAdcMax-fAdcMin;
638
639   TVectorD param(4);
640   TMatrixD dummy(3,3);
641
642   TH1F *hChannel=new TH1F("hChannel","hChannel",nbinsAdc,fAdcMin,fAdcMax);
643   
644   Float_t *array_hP=0;  
645
646   for (Int_t iSec=0; iSec<72; ++iSec){
647     TH2F *hP = GetHistoPedestal(iSec);
648     if ( !hP ) continue;
649
650     AliTPCCalROC *rocPedestal = GetCalRocPedestal(iSec,kTRUE);
651     AliTPCCalROC *rocSigma    = GetCalRocSigma(iSec,kTRUE);
652     AliTPCCalROC *rocMean     = GetCalRocMean(iSec,kTRUE);
653     AliTPCCalROC *rocRMS      = GetCalRocRMS(iSec,kTRUE);
654
655     array_hP = hP->GetArray();
656     UInt_t nChannels = fROC->GetNChannels(iSec);
657
658     for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
659       Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
660       //calculate mean and sigma using a gaus fit
661       Double_t ret = AliMathBase::FitGaus(array_hP+offset,nbinsAdc,fAdcMin,fAdcMax,&param,&dummy);
662       // if the fitting failed set noise and pedestal to 0
663       // is now done in AliMathBase::FitGaus !
664 //       if ( ret == -4 ) {
665 //      param[1]=0;
666 //      param[2]=0;
667 //       }
668       rocPedestal->SetValue(iChannel,param[1]);
669       rocSigma->SetValue(iChannel,param[2]);
670       //calculate mean and RMS using a truncated means
671       hChannel->Set(nbinsAdc+2,array_hP+offset-1);
672       hChannel->SetEntries(param[3]);
673       param[1]=0;
674       param[2]=0;
675       if ( param[3]>0 ) AliMathBase::TruncatedMean(hChannel,&param,fAnaMeanDown,fAnaMeanUp);
676       rocMean->SetValue(iChannel,param[1]);
677       rocRMS->SetValue(iChannel,param[2]);
678     }
679   }
680   delete hChannel;
681 }
682
683
684 //_____________________________________________________________________
685 void AliTPCCalibPedestal::AnalyseTime(Int_t nevents)
686 {
687   //
688   // Calculate for each channel and time bin the mean pedestal. This
689   // is used on LDC by TPCPEDESTALda to generate data used for configuration
690   // of the pattern memory for baseline subtraction in the ALTROs.
691   //
692
693   if ( nevents <= 0 ) return;
694   if ( fTimeAnalysis ) {
695     for (Int_t i = 0; i < 159; i++) {  // padrows
696       for (Int_t j = 0; j < 140; j++) {  // pads per row
697         for (Int_t k = 0; k < 1024; k++) {  // time bins per pad
698           fTimeSignal[i][j].AddAt(fTimeSignal[i][j].At(k)/(Float_t)nevents, k);
699         }
700       }
701     }
702   }
703 }
704
705
706 //_____________________________________________________________________
707 void AliTPCCalibPedestal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) 
708 {
709   //
710   //  Write class to file
711   //
712
713   TString sDir(dir);
714   TString option;
715
716   if ( append )
717     option = "update";
718   else
719     option = "recreate";
720
721   TDirectory *backup = gDirectory;
722   TFile f(filename,option.Data());
723   f.cd();
724   if ( !sDir.IsNull() ){
725     f.mkdir(sDir.Data());
726     f.cd(sDir);
727   }
728   this->Write();
729   f.Close();
730
731   if ( backup ) backup->cd();
732 }