Additional include files (Matthias)
[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
39 //date
40 #include "event.h"
41
42 //header file
43 #include "AliTPCCalibPedestal.h"
44
45
46 ///////////////////////////////////////////////////////////////////////////////////////
47 //          Implementation of the TPC pedestal and noise calibration
48 //
49 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
50 // 
51 // 
52 // *************************************************************************************
53 // *                                Class Description                                  *
54 // *************************************************************************************
55 //
56 // Working principle:
57 // ------------------
58 // Raw pedestal data is processed by calling one of the ProcessEvent(...) functions
59 // (see below). These in the end call the Update(...) function, where the data is filled
60 // into histograms.
61 //
62 // For each ROC one TH2F histo (ROC channel vs. ADC channel) is created when
63 // it is filled for the first time (GetHistoPedestal(ROC,kTRUE)). All histos are stored in the
64 // TObjArray fHistoPedestalArray.
65 //
66 // For a fast filling of the histogram the corresponding bin number of the channel and ADC channel
67 // is computed by hand and the histogram array is accessed directly via its pointer.
68 // ATTENTION: Doing so the the entry counter of the histogram is not increased
69 //            this means that e.g. the colz draw option gives an empty plot unless
70 //          calling 'histo->SetEntries(1)' before drawing.
71 //
72 // After accumulating the desired statistics the Analyse() function has to be called.
73 // Whithin this function the pedestal and noise values are calculated for each pad, using
74 // the fast gaus fit function  AliMathBase::FitGaus(...), and the calibration
75 // storage classes (AliTPCCalROC) are filled for each ROC.
76 // The calibration information is stored in the TObjArrays fCalRocArrayPedestal and fCalRocArrayRMS;
77 //
78 //
79 //
80 // User interface for filling data:
81 // --------------------------------
82 //
83 // To Fill information one of the following functions can be used:
84 //
85 // Bool_t ProcessEvent(eventHeaderStruct *event);
86 //   - process Date event
87 //   - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
88 //
89 // Bool_t ProcessEvent(AliRawReader *rawReader);
90 //  - process AliRawReader event
91 //   - use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)
92 //
93 // Bool_t ProcessEvent(AliTPCRawStream *rawStream);
94 //   - process event from AliTPCRawStream
95 //   - call Update function for signal filling
96 //
97 // Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
98 //              iPad,  const Int_t iTimeBin, const Float_t signal);
99 //   - directly  fill signal information (sector, row, pad, time bin, pad)
100 //     to the reference histograms
101 //
102 // It is also possible to merge two independently taken calibrations using the function
103 //
104 // void Merge(AliTPCCalibPedestal *ped)
105 //   - copy histograms in 'ped' if the do not exist in this instance
106 //   - Add histograms in 'ped' to the histograms in this instance if the allready exist
107 //   - After merging call Analyse again!
108 //
109 //
110 //
111 // -- example: filling data using root raw data:
112 // void fillPedestal(Char_t *filename)
113 // {
114 //    rawReader = new AliRawReaderRoot(fileName);
115 //    if ( !rawReader ) return;
116 //    AliTPCCalibPedestal *calib = new AliTPCCalibPedestal;
117 //    while (rawReader->NextEvent()){
118 //      calib->ProcessEvent(rawReader);
119 //    }
120 //    calib->Analyse();
121 //    calib->DumpToFile("PedestalData.root");
122 //    delete rawReader;
123 //    delete calib;
124 // }
125 //
126 //
127 // What kind of information is stored and how to retrieve them:
128 // ------------------------------------------------------------
129 //
130 // - Accessing the 'Reference Histograms' (pedestal distribution histograms):
131 //
132 // TH2F *GetHistoPedestal(Int_t sector);
133 //
134 // - Accessing the calibration storage objects:
135 //
136 // AliTPCCalROC *GetCalRocPedestal(Int_t sector);  - for the pedestal values
137 // AliTPCCalROC *GetCalRocNoise(Int_t sector);     - for the Noise values
138 //
139 // example for visualisation:
140 // if the file "PedestalData.root" was created using the above example one could do the following:
141 //
142 // TFile filePedestal("PedestalData.root")
143 // AliTPCCalibPedestal *ped = (AliTPCCalibPedestal*)filePedestal->Get("AliTPCCalibPedestal");
144 // ped->GetCalRocPedestal(0)->Draw("colz");
145 // ped->GetCalRocRMS(0)->Draw("colz");
146 //
147 // or use the AliTPCCalPad functionality:
148 // AliTPCCalPad padPedestal(ped->GetCalPadPedestal());
149 // AliTPCCalPad padNoise(ped->GetCalPadRMS());
150 // padPedestal->MakeHisto2D()->Draw("colz");  //Draw A-Side Pedestal Information
151 // padNoise->MakeHisto2D()->Draw("colz");  //Draw A-Side Noise Information
152 //
153 /*
154  example: fill pedestal with gausschen noise
155  AliTPCCalibPedestal ped;
156  ped.TestEvent();
157  ped.Analyse();
158  //Draw output;
159  TCanvas* c1 = new TCanvas;
160  c1->Divide(1,2);
161  c1->cd(1);
162  ped.GetHistoPedestal(0)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
163  ped.GetHistoPedestal(0)->Draw("colz");
164  c1->cd(2);
165  ped.GetHistoPedestal(36)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
166  ped.GetHistoPedestal(36)->Draw("colz");
167  TCanvas* c2 = new TCanvas;
168  c2->Divide(2,2);
169  c2->cd(1);
170  ped.GetCalRocPedestal(0)->Draw("colz");
171  c2->cd(2);
172  ped.GetCalRocRMS(0)->Draw("colz");
173  c2->cd(3);
174  ped.GetCalRocPedestal(36)->Draw("colz");
175  c2->cd(4);
176  ped.GetCalRocRMS(36)->Draw("colz");
177
178
179 */
180
181
182
183 ClassImp(AliTPCCalibPedestal)
184
185 AliTPCCalibPedestal::AliTPCCalibPedestal() : /*FOLD00*/
186   TObject(),
187   fFirstTimeBin(60),
188   fLastTimeBin(1000),
189   fAdcMin(1),
190   fAdcMax(100),
191   fOldRCUformat(kTRUE),
192   fROC(AliTPCROC::Instance()),
193   fCalRocArrayPedestal(72),
194   fCalRocArrayRMS(72),
195   fHistoPedestalArray(72)
196 {
197     //
198     // default constructor
199     //
200 }
201 //_____________________________________________________________________
202 AliTPCCalibPedestal::AliTPCCalibPedestal(const AliTPCCalibPedestal &ped) : /*FOLD00*/
203   TObject(ped),
204   fFirstTimeBin(ped.GetFirstTimeBin()),
205   fLastTimeBin(ped.GetLastTimeBin()),
206   fAdcMin(ped.GetAdcMin()),
207   fAdcMax(ped.GetAdcMax()),
208   fOldRCUformat(kTRUE),
209   fROC(AliTPCROC::Instance()),
210   fCalRocArrayPedestal(72),
211   fCalRocArrayRMS(72),
212   fHistoPedestalArray(72)
213 {
214     //
215     // copy constructor
216     //
217     for (Int_t iSec = 0; iSec < 72; ++iSec){
218         const AliTPCCalROC *calPed = (AliTPCCalROC*)ped.fCalRocArrayPedestal.UncheckedAt(iSec);
219         const AliTPCCalROC *calRMS = (AliTPCCalROC*)ped.fCalRocArrayRMS.UncheckedAt(iSec);
220         const TH2F         *hPed   = (TH2F*)ped.fHistoPedestalArray.UncheckedAt(iSec);
221
222         if ( calPed != 0x0 ) fCalRocArrayPedestal.AddAt(new AliTPCCalROC(*calPed), iSec);
223         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
224
225         if ( hPed != 0x0 ){
226             TH2F *hNew = new TH2F(*hPed);
227             hNew->SetDirectory(0);
228             fHistoPedestalArray.AddAt(hNew,iSec);
229         }
230     }
231 }
232 //_____________________________________________________________________
233 AliTPCCalibPedestal& AliTPCCalibPedestal::operator = (const  AliTPCCalibPedestal &source)
234 {
235   //
236   // assignment operator
237   //
238   if (&source == this) return *this;
239   new (this) AliTPCCalibPedestal(source);
240
241   return *this;
242 }
243 //_____________________________________________________________________
244 AliTPCCalibPedestal::~AliTPCCalibPedestal() /*FOLD00*/
245 {
246   //
247   // destructor
248   //
249
250     fCalRocArrayPedestal.Delete();
251     fCalRocArrayRMS.Delete();
252     fHistoPedestalArray.Delete();
253     delete fROC;
254 }
255 //_____________________________________________________________________
256 Int_t AliTPCCalibPedestal::Update(const Int_t icsector, /*FOLD00*/
257                                 const Int_t icRow,
258                                 const Int_t icPad,
259                                 const Int_t icTimeBin,
260                                 const Float_t csignal)
261 {
262     //
263     // Signal filling methode 
264     //
265
266     //return if we are out of the specified time bin or adc range
267     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
268     if ( ((Int_t)csignal>fAdcMax) || ((Int_t)csignal<fAdcMin)  ) return 0;
269
270     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
271
272     // fast filling methode.
273     // Attention: the entry counter of the histogram is not increased
274     //            this means that e.g. the colz draw option gives an empty plot
275     Int_t bin = (iChannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
276
277     GetHistoPedestal(icsector,kTRUE)->GetArray()[bin]++;
278
279     return 0;
280 }
281 //_____________________________________________________________________
282 Bool_t AliTPCCalibPedestal::ProcessEvent(AliTPCRawStream *rawStream)
283 {
284   //
285   // Event Processing loop - AliTPCRawStream
286   //
287
288   rawStream->SetOldRCUFormat(fOldRCUformat);
289
290   Bool_t withInput = kFALSE;
291
292   while (rawStream->Next()) {
293
294       Int_t isector  = rawStream->GetSector();                       //  current sector
295       Int_t iRow     = rawStream->GetRow();                          //  current row
296       Int_t iPad     = rawStream->GetPad();                          //  current pad
297       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
298       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
299
300       Update(isector,iRow,iPad,iTimeBin,signal);
301       withInput = kTRUE;
302   }
303
304   return withInput;
305 }
306 //_____________________________________________________________________
307 Bool_t AliTPCCalibPedestal::ProcessEvent(AliRawReader *rawReader)
308 {
309   //
310   //  Event processing loop - AliRawReader
311   //
312
313
314   AliTPCRawStream rawStream(rawReader);
315
316   rawReader->Select("TPC");
317
318   return ProcessEvent(&rawStream);
319 }
320 //_____________________________________________________________________
321 Bool_t AliTPCCalibPedestal::ProcessEvent(eventHeaderStruct *event)
322 {
323   //
324   //  process date event
325   //
326     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
327     Bool_t result=ProcessEvent(rawReader);
328     delete rawReader;
329     return result;
330 }
331 //_____________________________________________________________________
332 Bool_t AliTPCCalibPedestal::TestEvent() /*FOLD00*/
333 {
334   //
335   //  Test event loop
336   // fill one oroc and one iroc with random gaus
337   //
338
339     gRandom->SetSeed(0);
340
341     for (UInt_t iSec=0; iSec<72; ++iSec){
342         if (iSec%36>0) continue;
343         for (UInt_t iRow=0; iRow < fROC->GetNRows(iSec); ++iRow){
344             for (UInt_t iPad=0; iPad < fROC->GetNPads(iSec,iRow); ++iPad){
345                 for (UInt_t iTimeBin=0; iTimeBin<1024; ++iTimeBin){
346                     Float_t signal=(Int_t)(iRow+3+gRandom->Gaus(0,.7));
347                     if ( signal>0 )Update(iSec,iRow,iPad,iTimeBin,signal);
348                 }
349             }
350         }
351     }
352     return kTRUE;
353 }
354 //_____________________________________________________________________
355 TH2F* AliTPCCalibPedestal::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
356                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
357                                   Char_t *type, Bool_t force)
358 {
359     //
360     // return pointer to Q histogram
361     // if force is true create a new histogram if it doesn't exist allready
362     //
363     if ( !force || arr->UncheckedAt(sector) )
364         return (TH2F*)arr->UncheckedAt(sector);
365
366     // if we are forced and histogram doesn't yes exist create it
367     Char_t name[255], title[255];
368
369     sprintf(name,"hCalib%s%.2d",type,sector);
370     sprintf(title,"%s calibration histogram sector %.2d;ADC channel;Channel (pad)",type,sector);
371
372     // new histogram with Q calib information. One value for each pad!
373     TH2F* hist = new TH2F(name,title,
374                           nbinsY, ymin, ymax,
375                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector)
376                          );
377     hist->SetDirectory(0);
378     arr->AddAt(hist,sector);
379     return hist;
380 }
381 //_____________________________________________________________________
382 TH2F* AliTPCCalibPedestal::GetHistoPedestal(Int_t sector, Bool_t force) /*FOLD00*/
383 {
384     //
385     // return pointer to T0 histogram
386     // if force is true create a new histogram if it doesn't exist allready
387     //
388     TObjArray *arr = &fHistoPedestalArray;
389     return GetHisto(sector, arr, fAdcMax-fAdcMin, fAdcMin, fAdcMax, "Pedestal", force);
390 }
391 //_____________________________________________________________________
392 AliTPCCalROC* AliTPCCalibPedestal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
393 {
394     //
395     // return pointer to ROC Calibration
396     // if force is true create a new histogram if it doesn't exist allready
397     //
398     if ( !force || arr->UncheckedAt(sector) )
399         return (AliTPCCalROC*)arr->UncheckedAt(sector);
400
401     // if we are forced and the histogram doesn't yet exist create it
402
403     // new AliTPCCalROC for T0 information. One value for each pad!
404     AliTPCCalROC *croc = new AliTPCCalROC(sector);
405     arr->AddAt(croc,sector);
406     return croc;
407 }
408 //_____________________________________________________________________
409 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocPedestal(Int_t sector, Bool_t force) /*FOLD00*/
410 {
411     //
412     // return pointer to Carge ROC Calibration
413     // if force is true create a new histogram if it doesn't exist allready
414     //
415     TObjArray *arr = &fCalRocArrayPedestal;
416     return GetCalRoc(sector, arr, force);
417 }
418 //_____________________________________________________________________
419 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
420 {
421     //
422     // return pointer to signal width ROC Calibration
423     // if force is true create a new histogram if it doesn't exist allready
424     //
425     TObjArray *arr = &fCalRocArrayRMS;
426     return GetCalRoc(sector, arr, force);
427 }
428 //_____________________________________________________________________
429 void AliTPCCalibPedestal::Merge(AliTPCCalibPedestal *ped)
430 {
431     //
432     //  Merge reference histograms of sig to the current AliTPCCalibSignal
433     //
434
435     //merge histograms
436     for (Int_t iSec=0; iSec<72; ++iSec){
437         TH2F *hRefPedMerge   = ped->GetHistoPedestal(iSec);
438
439
440         if ( hRefPedMerge ){
441             TDirectory *dir = hRefPedMerge->GetDirectory(); hRefPedMerge->SetDirectory(0);
442             TH2F *hRefPed   = GetHistoPedestal(iSec);
443             if ( hRefPed ) hRefPed->Add(hRefPedMerge);
444             else {
445                 TH2F *hist = new TH2F(*hRefPedMerge);
446                 hist->SetDirectory(0);
447                 fHistoPedestalArray.AddAt(hist, iSec);
448             }
449             hRefPedMerge->SetDirectory(dir);
450         }
451     }
452 }
453 //_____________________________________________________________________
454 void AliTPCCalibPedestal::Analyse() /*FOLD00*/
455 {
456     //
457     //  Calculate calibration constants
458     //
459
460     Int_t nbinsAdc = fAdcMax-fAdcMin;
461
462     TVectorD param(3);
463     TMatrixD dummy(3,3);
464
465     Float_t *array_hP=0;
466
467
468     for (Int_t iSec=0; iSec<72; ++iSec){
469         TH2F *hP = GetHistoPedestal(iSec);
470         if ( !hP ) continue;
471
472         AliTPCCalROC *rocPedestal = GetCalRocPedestal(iSec,kTRUE);
473         AliTPCCalROC *rocRMS      = GetCalRocRMS(iSec,kTRUE);
474
475         array_hP = hP->GetArray();
476         UInt_t nChannels = fROC->GetNChannels(iSec);
477
478         for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
479             Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
480             Double_t ret = AliMathBase::FitGaus(array_hP+offset,nbinsAdc,fAdcMin,fAdcMax,&param,&dummy);
481             // if the fitting failed set noise and pedestal to 0
482             if ( ret == -4 ) {
483                 param[1]=0;
484                 param[2]=0;
485             }
486             rocPedestal->SetValue(iChannel,param[1]);
487             rocRMS->SetValue(iChannel,param[2]);
488         }
489     }
490 }
491 //_____________________________________________________________________
492 void AliTPCCalibPedestal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
493 {
494     //
495     //  Write class to file
496     //
497
498     TString sDir(dir);
499     TString option;
500
501     if ( append )
502         option = "update";
503     else
504         option = "recreate";
505
506     TDirectory *backup = gDirectory;
507     TFile f(filename,option.Data());
508     f.cd();
509     if ( !sDir.IsNull() ){
510         f.mkdir(sDir.Data());
511         f.cd(sDir);
512     }
513     this->Write();
514     f.Close();
515
516     if ( backup ) backup->cd();
517 }