a39a7358cf4bbb5517d2424f245df230f2229f77
[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 //-------------------------------------------------------
18 //          Implementation of the TPC pedestal and noise calibration
19 //
20 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
21 // 
22 // 
23 //-------------------------------------------------------
24
25
26 /* $Id$ */
27
28 /*
29  example: fill pedestal with gausschen noise
30  AliTPCCalibPedestal ped;
31  ped.TestEvent();
32  ped.Analyse();
33  //Draw output;
34  TCanvas* c1 = new TCanvas;
35  c1->Divide(1,2);
36  c1->cd(1);
37  ped.GetHistoPedestal(0)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
38  ped.GetHistoPedestal(0)->Draw("colz");
39  c1->cd(2);
40  ped.GetHistoPedestal(36)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
41  ped.GetHistoPedestal(36)->Draw("colz");
42  TCanvas* c2 = new TCanvas;
43  c2->Divide(2,2);
44  c2->cd(1);
45  ped.GetCalRocPedestal(0)->Draw("colz");
46  c2->cd(2);
47  ped.GetCalRocRMS(0)->Draw("colz");
48  c2->cd(3);
49  ped.GetCalRocPedestal(36)->Draw("colz");
50  c2->cd(4);
51  ped.GetCalRocRMS(36)->Draw("colz");
52
53
54 */
55
56 //Root includes
57 #include <TObjArray.h>
58 #include <TH1F.h>
59 #include <TH2F.h>
60 #include <TString.h>
61 #include <TMath.h>
62 #include <TF1.h>
63 #include <TRandom.h>
64 #include <TDirectory.h>
65 #include <TFile.h>
66 //AliRoot includes
67 #include "AliRawReader.h"
68 #include "AliRawReaderRoot.h"
69 #include "AliRawReaderDate.h"
70 #include "AliTPCRawStream.h"
71 #include "AliTPCCalROC.h"
72 #include "AliTPCROC.h"
73 #include "AliMathBase.h"
74 #include "TTreeStream.h"
75
76 //date
77 #include "event.h"
78
79 //header file
80 #include "AliTPCCalibPedestal.h"
81
82
83 ClassImp(AliTPCCalibPedestal) /*FOLD00*/
84
85
86
87
88
89
90
91 AliTPCCalibPedestal::AliTPCCalibPedestal() : /*FOLD00*/
92   TObject(),
93   fFirstTimeBin(60),
94   fLastTimeBin(1000),
95   fAdcMin(1),
96   fAdcMax(100),
97   fROC(AliTPCROC::Instance()),
98   fCalRocArrayPedestal(72),
99   fCalRocArrayRMS(72),
100   fHistoPedestalArray(72)
101 {
102     //
103     // default constructor
104     //
105 }
106 //_____________________________________________________________________
107 AliTPCCalibPedestal::AliTPCCalibPedestal(const AliTPCCalibPedestal &ped) : /*FOLD00*/
108   TObject(ped),
109   fFirstTimeBin(ped.GetFirstTimeBin()),
110   fLastTimeBin(ped.GetLastTimeBin()),
111   fAdcMin(ped.GetAdcMin()),
112   fAdcMax(ped.GetAdcMax()),
113   fROC(AliTPCROC::Instance()),
114   fCalRocArrayPedestal(72),
115   fCalRocArrayRMS(72),
116   fHistoPedestalArray(72)
117 {
118     //
119     // copy constructor
120     //
121     for (Int_t iSec = 0; iSec < 72; iSec++){
122         const AliTPCCalROC *calPed = (AliTPCCalROC*)ped.fCalRocArrayPedestal.UncheckedAt(iSec);
123         const AliTPCCalROC *calRMS = (AliTPCCalROC*)ped.fCalRocArrayRMS.UncheckedAt(iSec);
124         const TH2F         *hPed   = (TH2F*)ped.fHistoPedestalArray.UncheckedAt(iSec);
125
126         if ( calPed != 0x0 ) fCalRocArrayPedestal.AddAt(new AliTPCCalROC(*calPed), iSec);
127         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
128
129         if ( hPed != 0x0 ){
130             TH2F *hNew = new TH2F(*hPed);
131             hNew->SetDirectory(0);
132             fHistoPedestalArray.AddAt(hNew,iSec);
133         }
134     }
135 }
136 //_____________________________________________________________________
137 AliTPCCalibPedestal& AliTPCCalibPedestal::operator = (const  AliTPCCalibPedestal &source)
138 {
139   //
140   // assignment operator
141   //
142   if (&source == this) return *this;
143   new (this) AliTPCCalibPedestal(source);
144
145   return *this;
146 }
147 //_____________________________________________________________________
148 AliTPCCalibPedestal::~AliTPCCalibPedestal() /*FOLD00*/
149 {
150   //
151   // destructor
152   //
153   delete fROC;
154 }
155
156
157
158
159 //_____________________________________________________________________
160 Int_t AliTPCCalibPedestal::Update(const Int_t icsector, /*FOLD00*/
161                                 const Int_t icRow,
162                                 const Int_t icPad,
163                                 const Int_t icTimeBin,
164                                 const Float_t csignal)
165 {
166     //
167     // Signal filling methode 
168     //
169     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
170
171     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
172
173     // fast filling methode.
174     // Attention: the entry counter of the histogram is not increased
175     //            this means that e.g. the colz draw option gives an empty plot
176     Int_t bin = 0;
177     if ( !(((Int_t)csignal>fAdcMax ) || ((Int_t)csignal<fAdcMin)) )
178         bin = (iChannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
179
180     GetHistoPedestal(icsector,kTRUE)->GetArray()[bin]++;
181
182     return 0;
183 }
184 //_____________________________________________________________________
185 Bool_t AliTPCCalibPedestal::ProcessEvent(AliTPCRawStream *rawStream)
186 {
187   //
188   // Event Processing loop - AliTPCRawStream
189   //
190
191   rawStream->SetOldRCUFormat(1);
192
193   Bool_t withInput = kFALSE;
194
195   while (rawStream->Next()) {
196
197       Int_t isector  = rawStream->GetSector();                       //  current sector
198       Int_t iRow     = rawStream->GetRow();                          //  current row
199       Int_t iPad     = rawStream->GetPad();                          //  current pad
200       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
201       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
202
203       Update(isector,iRow,iPad,iTimeBin,signal);
204       withInput = kTRUE;
205   }
206
207   return withInput;
208 }
209 //_____________________________________________________________________
210 Bool_t AliTPCCalibPedestal::ProcessEvent(AliRawReader *rawReader)
211 {
212   //
213   //  Event processing loop - AliRawReader
214   //
215
216
217   AliTPCRawStream rawStream(rawReader);
218
219   rawReader->Select("TPC");
220
221   return ProcessEvent(&rawStream);
222 }
223 //_____________________________________________________________________
224 Bool_t AliTPCCalibPedestal::ProcessEvent(eventHeaderStruct *event)
225 {
226   //
227   //  process date event
228   //
229     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
230     Bool_t result=ProcessEvent(rawReader);
231     delete rawReader;
232     return result;
233 }
234 //_____________________________________________________________________
235 Bool_t AliTPCCalibPedestal::TestEvent() /*FOLD00*/
236 {
237   //
238   //  Test event loop
239   // fill one oroc and one iroc with random gaus
240   //
241
242     gRandom->SetSeed(0);
243
244     for (UInt_t iSec=0; iSec<72; iSec++){
245         if (iSec%36>0) continue;
246         for (UInt_t iRow=0; iRow < fROC->GetNRows(iSec); iRow++){
247             for (UInt_t iPad=0; iPad < fROC->GetNPads(iSec,iRow); iPad++){
248                 for (UInt_t iTimeBin=0; iTimeBin<1024; iTimeBin++){
249                     Float_t signal=(Int_t)(iRow+3+gRandom->Gaus(0,.7));
250                     if ( signal>0 )Update(iSec,iRow,iPad,iTimeBin,signal);
251                 }
252             }
253         }
254     }
255     return kTRUE;
256 }
257 //_____________________________________________________________________
258 TH2F* AliTPCCalibPedestal::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
259                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
260                                   Char_t *type, Bool_t force)
261 {
262     //
263     // return pointer to Q histogram
264     // if force is true create a new histogram if it doesn't exist allready
265     //
266     if ( !force || arr->UncheckedAt(sector) )
267         return (TH2F*)arr->UncheckedAt(sector);
268
269     // if we are forced and histogram doesn't yes exist create it
270     Char_t name[255], title[255];
271
272     sprintf(name,"hCalib%s%.2d",type,sector);
273     sprintf(title,"%s calibration histogram sector %.2d;ADC channel;Channel (pad)",type,sector);
274
275     // new histogram with Q calib information. One value for each pad!
276     TH2F* hist = new TH2F(name,title,
277                           nbinsY, ymin, ymax,
278                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector)
279                          );
280     hist->SetDirectory(0);
281     arr->AddAt(hist,sector);
282     return hist;
283 }
284 //_____________________________________________________________________
285 TH2F* AliTPCCalibPedestal::GetHistoPedestal(Int_t sector, Bool_t force) /*FOLD00*/
286 {
287     //
288     // return pointer to T0 histogram
289     // if force is true create a new histogram if it doesn't exist allready
290     //
291     TObjArray *arr = &fHistoPedestalArray;
292     return GetHisto(sector, arr, fAdcMax-fAdcMin, fAdcMin, fAdcMax, "Pedestal", force);
293 }
294 //_____________________________________________________________________
295 AliTPCCalROC* AliTPCCalibPedestal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
296 {
297     //
298     // return pointer to ROC Calibration
299     // if force is true create a new histogram if it doesn't exist allready
300     //
301     if ( !force || arr->UncheckedAt(sector) )
302         return (AliTPCCalROC*)arr->UncheckedAt(sector);
303
304     // if we are forced and histogram doesn't yes exist create it
305
306     // new AliTPCCalROC for T0 information. One value for each pad!
307     AliTPCCalROC *croc = new AliTPCCalROC(sector);
308     //init values
309     for ( UInt_t iChannel = 0; iChannel<croc->GetNchannels(); iChannel++){
310         croc->SetValue(iChannel, 0);
311     }
312     arr->AddAt(croc,sector);
313     return croc;
314 }
315 //_____________________________________________________________________
316 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocPedestal(Int_t sector, Bool_t force) /*FOLD00*/
317 {
318     //
319     // return pointer to Carge ROC Calibration
320     // if force is true create a new histogram if it doesn't exist allready
321     //
322     TObjArray *arr = &fCalRocArrayPedestal;
323     return GetCalRoc(sector, arr, force);
324 }
325 //_____________________________________________________________________
326 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
327 {
328     //
329     // return pointer to signal width ROC Calibration
330     // if force is true create a new histogram if it doesn't exist allready
331     //
332     TObjArray *arr = &fCalRocArrayRMS;
333     return GetCalRoc(sector, arr, force);
334 }
335 //_____________________________________________________________________
336 void AliTPCCalibPedestal::Analyse() /*FOLD00*/
337 {
338     //
339     //  Calculate calibration constants
340     //
341
342     Int_t nbinsAdc = fAdcMax-fAdcMin;
343
344     TVectorD param(3);
345     TMatrixD dummy(3,3);
346
347     Float_t *array_hP=0;
348
349
350     for (Int_t iSec=0; iSec<72; iSec++){
351         TH2F *hP = GetHistoPedestal(iSec);
352         if ( !hP ) continue;
353
354         AliTPCCalROC *rocPedestal = GetCalRocPedestal(iSec,kTRUE);
355         AliTPCCalROC *rocRMS      = GetCalRocRMS(iSec,kTRUE);
356
357         array_hP = hP->GetArray();
358         UInt_t nChannels = fROC->GetNChannels(iSec);
359
360         for (UInt_t iChannel=0; iChannel<nChannels; iChannel++){
361             Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
362             Double_t ret = AliMathBase::FitGaus(array_hP+offset,nbinsAdc,fAdcMin,fAdcMax,&param,&dummy);
363             // if the fitting failed set noise and pedestal to 0
364             if ( ret == -4 ) {
365                 param[1]=0;
366                 param[2]=0;
367             }
368             rocPedestal->SetValue(iChannel,param[1]);
369             rocRMS->SetValue(iChannel,param[2]);
370         }
371     }
372 }
373 //_____________________________________________________________________
374 void AliTPCCalibPedestal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
375 {
376     //
377     //  Write class to file
378     //
379
380     TString sDir(dir);
381     TString option;
382
383     if ( append )
384         option = "update";
385     else
386         option = "recreate";
387
388     TDirectory *backup = gDirectory;
389     TFile f(filename,option.Data());
390     f.cd();
391     if ( !sDir.IsNull() ){
392         f.mkdir(sDir.Data());
393         f.cd(sDir);
394     }
395     this->Write();
396     f.Close();
397
398     if ( backup ) backup->cd();
399 }