142a831e9a136362e9afc0ca79851c4b861e5640
[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 <TObjArray.h>
22 #include <TH1F.h>
23 #include <TH2F.h>
24 #include <TString.h>
25 #include <TMath.h>
26 #include <TF1.h>
27 #include <TRandom.h>
28 #include <TDirectory.h>
29 #include <TFile.h>
30 //AliRoot includes
31 #include "AliRawReader.h"
32 #include "AliRawReaderRoot.h"
33 #include "AliRawReaderDate.h"
34 #include "AliTPCRawStream.h"
35 #include "AliTPCCalROC.h"
36 #include "AliTPCROC.h"
37 #include "AliMathBase.h"
38 #include "TTreeStream.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
138 // AliTPCCalROC *GetCalRocNoise(Int_t sector);     - for the Noise values
139 //
140 // example for visualisation:
141 // if the file "PedestalData.root" was created using the above example one could do the following:
142 //
143 // TFile filePedestal("PedestalData.root")
144 // AliTPCCalibPedestal *ped = (AliTPCCalibPedestal*)filePedestal->Get("AliTPCCalibPedestal");
145 // ped->GetCalRocPedestal(0)->Draw("colz");
146 // ped->GetCalRocRMS(0)->Draw("colz");
147 //
148 // or use the AliTPCCalPad functionality:
149 // AliTPCCalPad padPedestal(ped->GetCalPadPedestal());
150 // AliTPCCalPad padNoise(ped->GetCalPadRMS());
151 // padPedestal->MakeHisto2D()->Draw("colz");  //Draw A-Side Pedestal Information
152 // padNoise->MakeHisto2D()->Draw("colz");  //Draw A-Side Noise Information
153 //
154 /*
155  example: fill pedestal with gausschen noise
156  AliTPCCalibPedestal ped;
157  ped.TestEvent();
158  ped.Analyse();
159  //Draw output;
160  TCanvas* c1 = new TCanvas;
161  c1->Divide(1,2);
162  c1->cd(1);
163  ped.GetHistoPedestal(0)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
164  ped.GetHistoPedestal(0)->Draw("colz");
165  c1->cd(2);
166  ped.GetHistoPedestal(36)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
167  ped.GetHistoPedestal(36)->Draw("colz");
168  TCanvas* c2 = new TCanvas;
169  c2->Divide(2,2);
170  c2->cd(1);
171  ped.GetCalRocPedestal(0)->Draw("colz");
172  c2->cd(2);
173  ped.GetCalRocRMS(0)->Draw("colz");
174  c2->cd(3);
175  ped.GetCalRocPedestal(36)->Draw("colz");
176  c2->cd(4);
177  ped.GetCalRocRMS(36)->Draw("colz");
178
179
180 */
181
182
183
184 ClassImp(AliTPCCalibPedestal)
185
186 AliTPCCalibPedestal::AliTPCCalibPedestal() : /*FOLD00*/
187   TObject(),
188   fFirstTimeBin(60),
189   fLastTimeBin(1000),
190   fAdcMin(1),
191   fAdcMax(100),
192   fOldRCUformat(kTRUE),
193   fROC(AliTPCROC::Instance()),
194   fCalRocArrayPedestal(72),
195   fCalRocArrayRMS(72),
196   fHistoPedestalArray(72)
197 {
198     //
199     // default constructor
200     //
201 }
202 //_____________________________________________________________________
203 AliTPCCalibPedestal::AliTPCCalibPedestal(const AliTPCCalibPedestal &ped) : /*FOLD00*/
204   TObject(ped),
205   fFirstTimeBin(ped.GetFirstTimeBin()),
206   fLastTimeBin(ped.GetLastTimeBin()),
207   fAdcMin(ped.GetAdcMin()),
208   fAdcMax(ped.GetAdcMax()),
209   fOldRCUformat(kTRUE),
210   fROC(AliTPCROC::Instance()),
211   fCalRocArrayPedestal(72),
212   fCalRocArrayRMS(72),
213   fHistoPedestalArray(72)
214 {
215     //
216     // copy constructor
217     //
218     for (Int_t iSec = 0; iSec < 72; ++iSec){
219         const AliTPCCalROC *calPed = (AliTPCCalROC*)ped.fCalRocArrayPedestal.UncheckedAt(iSec);
220         const AliTPCCalROC *calRMS = (AliTPCCalROC*)ped.fCalRocArrayRMS.UncheckedAt(iSec);
221         const TH2F         *hPed   = (TH2F*)ped.fHistoPedestalArray.UncheckedAt(iSec);
222
223         if ( calPed != 0x0 ) fCalRocArrayPedestal.AddAt(new AliTPCCalROC(*calPed), iSec);
224         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
225
226         if ( hPed != 0x0 ){
227             TH2F *hNew = new TH2F(*hPed);
228             hNew->SetDirectory(0);
229             fHistoPedestalArray.AddAt(hNew,iSec);
230         }
231     }
232 }
233 //_____________________________________________________________________
234 AliTPCCalibPedestal& AliTPCCalibPedestal::operator = (const  AliTPCCalibPedestal &source)
235 {
236   //
237   // assignment operator
238   //
239   if (&source == this) return *this;
240   new (this) AliTPCCalibPedestal(source);
241
242   return *this;
243 }
244 //_____________________________________________________________________
245 AliTPCCalibPedestal::~AliTPCCalibPedestal() /*FOLD00*/
246 {
247   //
248   // destructor
249   //
250
251     fCalRocArrayPedestal.Delete();
252     fCalRocArrayRMS.Delete();
253     fHistoPedestalArray.Delete();
254     delete fROC;
255 }
256 //_____________________________________________________________________
257 Int_t AliTPCCalibPedestal::Update(const Int_t icsector, /*FOLD00*/
258                                 const Int_t icRow,
259                                 const Int_t icPad,
260                                 const Int_t icTimeBin,
261                                 const Float_t csignal)
262 {
263     //
264     // Signal filling methode 
265     //
266
267     //return if we are out of the specified time bin or adc range
268     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
269     if ( ((Int_t)csignal>fAdcMax) || ((Int_t)csignal<fAdcMin)  ) return 0;
270
271     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
272
273     // fast filling methode.
274     // Attention: the entry counter of the histogram is not increased
275     //            this means that e.g. the colz draw option gives an empty plot
276     Int_t bin = (iChannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
277
278     GetHistoPedestal(icsector,kTRUE)->GetArray()[bin]++;
279
280     return 0;
281 }
282 //_____________________________________________________________________
283 Bool_t AliTPCCalibPedestal::ProcessEvent(AliTPCRawStream *rawStream)
284 {
285   //
286   // Event Processing loop - AliTPCRawStream
287   //
288
289   rawStream->SetOldRCUFormat(fOldRCUformat);
290
291   Bool_t withInput = kFALSE;
292
293   while (rawStream->Next()) {
294
295       Int_t isector  = rawStream->GetSector();                       //  current sector
296       Int_t iRow     = rawStream->GetRow();                          //  current row
297       Int_t iPad     = rawStream->GetPad();                          //  current pad
298       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
299       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
300
301       Update(isector,iRow,iPad,iTimeBin,signal);
302       withInput = kTRUE;
303   }
304
305   return withInput;
306 }
307 //_____________________________________________________________________
308 Bool_t AliTPCCalibPedestal::ProcessEvent(AliRawReader *rawReader)
309 {
310   //
311   //  Event processing loop - AliRawReader
312   //
313
314
315   AliTPCRawStream rawStream(rawReader);
316
317   rawReader->Select("TPC");
318
319   return ProcessEvent(&rawStream);
320 }
321 //_____________________________________________________________________
322 Bool_t AliTPCCalibPedestal::ProcessEvent(eventHeaderStruct *event)
323 {
324   //
325   //  process date event
326   //
327     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
328     Bool_t result=ProcessEvent(rawReader);
329     delete rawReader;
330     return result;
331 }
332 //_____________________________________________________________________
333 Bool_t AliTPCCalibPedestal::TestEvent() /*FOLD00*/
334 {
335   //
336   //  Test event loop
337   // fill one oroc and one iroc with random gaus
338   //
339
340     gRandom->SetSeed(0);
341
342     for (UInt_t iSec=0; iSec<72; ++iSec){
343         if (iSec%36>0) continue;
344         for (UInt_t iRow=0; iRow < fROC->GetNRows(iSec); ++iRow){
345             for (UInt_t iPad=0; iPad < fROC->GetNPads(iSec,iRow); ++iPad){
346                 for (UInt_t iTimeBin=0; iTimeBin<1024; ++iTimeBin){
347                     Float_t signal=(Int_t)(iRow+3+gRandom->Gaus(0,.7));
348                     if ( signal>0 )Update(iSec,iRow,iPad,iTimeBin,signal);
349                 }
350             }
351         }
352     }
353     return kTRUE;
354 }
355 //_____________________________________________________________________
356 TH2F* AliTPCCalibPedestal::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
357                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
358                                   Char_t *type, Bool_t force)
359 {
360     //
361     // return pointer to Q histogram
362     // if force is true create a new histogram if it doesn't exist allready
363     //
364     if ( !force || arr->UncheckedAt(sector) )
365         return (TH2F*)arr->UncheckedAt(sector);
366
367     // if we are forced and histogram doesn't yes exist create it
368     Char_t name[255], title[255];
369
370     sprintf(name,"hCalib%s%.2d",type,sector);
371     sprintf(title,"%s calibration histogram sector %.2d;ADC channel;Channel (pad)",type,sector);
372
373     // new histogram with Q calib information. One value for each pad!
374     TH2F* hist = new TH2F(name,title,
375                           nbinsY, ymin, ymax,
376                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector)
377                          );
378     hist->SetDirectory(0);
379     arr->AddAt(hist,sector);
380     return hist;
381 }
382 //_____________________________________________________________________
383 TH2F* AliTPCCalibPedestal::GetHistoPedestal(Int_t sector, Bool_t force) /*FOLD00*/
384 {
385     //
386     // return pointer to T0 histogram
387     // if force is true create a new histogram if it doesn't exist allready
388     //
389     TObjArray *arr = &fHistoPedestalArray;
390     return GetHisto(sector, arr, fAdcMax-fAdcMin, fAdcMin, fAdcMax, "Pedestal", force);
391 }
392 //_____________________________________________________________________
393 AliTPCCalROC* AliTPCCalibPedestal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
394 {
395     //
396     // return pointer to ROC Calibration
397     // if force is true create a new histogram if it doesn't exist allready
398     //
399     if ( !force || arr->UncheckedAt(sector) )
400         return (AliTPCCalROC*)arr->UncheckedAt(sector);
401
402     // if we are forced and the histogram doesn't yet exist create it
403
404     // new AliTPCCalROC for T0 information. One value for each pad!
405     AliTPCCalROC *croc = new AliTPCCalROC(sector);
406     arr->AddAt(croc,sector);
407     return croc;
408 }
409 //_____________________________________________________________________
410 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocPedestal(Int_t sector, Bool_t force) /*FOLD00*/
411 {
412     //
413     // return pointer to Carge ROC Calibration
414     // if force is true create a new histogram if it doesn't exist allready
415     //
416     TObjArray *arr = &fCalRocArrayPedestal;
417     return GetCalRoc(sector, arr, force);
418 }
419 //_____________________________________________________________________
420 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
421 {
422     //
423     // return pointer to signal width ROC Calibration
424     // if force is true create a new histogram if it doesn't exist allready
425     //
426     TObjArray *arr = &fCalRocArrayRMS;
427     return GetCalRoc(sector, arr, force);
428 }
429 //_____________________________________________________________________
430 void AliTPCCalibPedestal::Merge(AliTPCCalibPedestal *ped)
431 {
432     //
433     //  Merge reference histograms of sig to the current AliTPCCalibSignal
434     //
435
436     //merge histograms
437     for (Int_t iSec=0; iSec<72; ++iSec){
438         TH2F *hRefPedMerge   = ped->GetHistoPedestal(iSec);
439
440
441         if ( hRefPedMerge ){
442             TDirectory *dir = hRefPedMerge->GetDirectory(); hRefPedMerge->SetDirectory(0);
443             TH2F *hRefPed   = GetHistoPedestal(iSec);
444             if ( hRefPed ) hRefPed->Add(hRefPedMerge);
445             else {
446                 TH2F *hist = new TH2F(*hRefPedMerge);
447                 hist->SetDirectory(0);
448                 fHistoPedestalArray.AddAt(hist, iSec);
449             }
450             hRefPedMerge->SetDirectory(dir);
451         }
452     }
453 }
454 //_____________________________________________________________________
455 void AliTPCCalibPedestal::Analyse() /*FOLD00*/
456 {
457     //
458     //  Calculate calibration constants
459     //
460
461     Int_t nbinsAdc = fAdcMax-fAdcMin;
462
463     TVectorD param(3);
464     TMatrixD dummy(3,3);
465
466     Float_t *array_hP=0;
467
468
469     for (Int_t iSec=0; iSec<72; ++iSec){
470         TH2F *hP = GetHistoPedestal(iSec);
471         if ( !hP ) continue;
472
473         AliTPCCalROC *rocPedestal = GetCalRocPedestal(iSec,kTRUE);
474         AliTPCCalROC *rocRMS      = GetCalRocRMS(iSec,kTRUE);
475
476         array_hP = hP->GetArray();
477         UInt_t nChannels = fROC->GetNChannels(iSec);
478
479         for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
480             Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
481             Double_t ret = AliMathBase::FitGaus(array_hP+offset,nbinsAdc,fAdcMin,fAdcMax,&param,&dummy);
482             // if the fitting failed set noise and pedestal to 0
483             if ( ret == -4 ) {
484                 param[1]=0;
485                 param[2]=0;
486             }
487             rocPedestal->SetValue(iChannel,param[1]);
488             rocRMS->SetValue(iChannel,param[2]);
489         }
490     }
491 }
492 //_____________________________________________________________________
493 void AliTPCCalibPedestal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
494 {
495     //
496     //  Write class to file
497     //
498
499     TString sDir(dir);
500     TString option;
501
502     if ( append )
503         option = "update";
504     else
505         option = "recreate";
506
507     TDirectory *backup = gDirectory;
508     TFile f(filename,option.Data());
509     f.cd();
510     if ( !sDir.IsNull() ){
511         f.mkdir(sDir.Data());
512         f.cd(sDir);
513     }
514     this->Write();
515     f.Close();
516
517     if ( backup ) backup->cd();
518 }