Add new class for Pedestal Calibration (M.Ivanov, J.Wiechula)
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibPedestal.cxx
CommitLineData
8bc7e885 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
30//Root includes
31#include <TObjArray.h>
32#include <TH1F.h>
33#include <TH1D.h>
34#include <TH2F.h>
35#include <TH2S.h>
36#include <TH1S.h>
37#include <TString.h>
38#include <TVectorF.h>
39#include <TMath.h>
40#include <TF1.h>
41#include <TRandom.h>
42#include <TROOT.h>
43#include <TDirectory.h>
44#include <TSystem.h>
45#include <TFile.h>
46
47
48//AliRoot includes
49#include "AliRawReader.h"
50#include "AliRawReaderRoot.h"
51#include "AliTPCRawStream.h"
52#include "AliTPCCalROC.h"
53#include "AliTPCCalPad.h"
54#include "AliTPCROC.h"
55#include "AliTPCCalibPedestal.h"
56#include "AliMathBase.h"
57
58#include "TTreeStream.h"
59
60
61
62ClassImp(AliTPCCalibPedestal) /*FOLD00*/
63
64
65
66
67
68
69
70AliTPCCalibPedestal::AliTPCCalibPedestal() : /*FOLD00*/
71 TObject(),
72 fFirstTimeBin(60),
73 fLastTimeBin(1000),
74 fAdcMin(1),
75 fAdcMax(100),
76 fROC(AliTPCROC::Instance()),
77 fCalRocArrayPedestal(72),
78 fCalRocArrayRMS(72),
79 fHistoPedestalArray(72),
80 fDebugStreamer(0),
81 fDebugLevel(0)
82{
83 //
84 // AliTPCSignal default constructor
85 //
86 //debug stream
87 TDirectory *backup = gDirectory;
88 fDebugStreamer = new TTreeSRedirector("deb2.root");
89 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
90}
91
92//_____________________________________________________________________
93AliTPCCalibPedestal::~AliTPCCalibPedestal() /*FOLD00*/
94{
95 //
96 // destructor
97 //
98 if ( fDebugStreamer ) delete fDebugStreamer;
99 delete fROC;
100}
101
102
103
104
105//_____________________________________________________________________
106Int_t AliTPCCalibPedestal::Update(const Int_t icsector, /*FOLD00*/
107 const Int_t icRow,
108 const Int_t icPad,
109 const Int_t icTimeBin,
110 const Float_t csignal)
111{
112 //
113 // Signal filling methode
114 //
115 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
116
117 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
118
119 // dirty and fast filling methode. No boundary checking!!!
120 // Attention: the entry counter of the histogram is not increased
121 // this means that e.g. the colz draw option gives an empty plot
122 Int_t bin = (iChannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
123
124 GetHistoPedestal(icsector,kTRUE)->GetArray()[bin]++;
125
126 return 0;
127}
128//_____________________________________________________________________
129Bool_t AliTPCCalibPedestal::ProcessEvent(AliRawReader *rawReader) /*FOLD00*/
130{
131 //
132 // simple event processing loop
133 //
134
135
136 AliTPCRawStream input(rawReader);
137
138 rawReader->Select("TPC");
139
140 input.SetOldRCUFormat(1);
141 printf("start event processing\n");
142
143
144 Bool_t withInput = kFALSE;
145
146 while (input.Next()) {
147
148 Int_t isector = input.GetSector(); // current sector
149 Int_t iRow = input.GetRow(); // current row
150 Int_t iPad = input.GetPad(); // current pad
151 Int_t iTimeBin = input.GetTime(); // current time bin
152 Float_t signal = input.GetSignal(); // current ADC signal
153
154 Update(isector,iRow,iPad,iTimeBin,signal);
155 withInput = kTRUE;
156 }
157
158 printf("end event processing\n");
159 if ( fDebugLevel>0 )
160 fDebugStreamer->GetFile()->Write();
161 return withInput;
162}
163//_____________________________________________________________________
164Bool_t AliTPCCalibPedestal::TestEvent() /*FOLD00*/
165{
166 //
167 // Test event loop
168 //
169
170 gRandom->SetSeed(0);
171
172 for (UInt_t iSec=0; iSec<72; iSec++){
173 if (iSec%36>0) continue;
174 for (UInt_t iRow=0; iRow < fROC->GetNRows(iSec); iRow++){
175 for (UInt_t iPad=0; iPad < fROC->GetNPads(iSec,iRow); iPad++){
176 for (UInt_t iTimeBin=0; iTimeBin<1024; iTimeBin++){
177 Float_t signal=(Int_t)(iRow+5+gRandom->Gaus(0,.7));
178 if ( signal>0 )Update(iSec,iRow,iPad,iTimeBin,signal);
179 }
180 }
181 }
182 }
183 return kTRUE;
184}
185//_____________________________________________________________________
186TH2S* AliTPCCalibPedestal::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
187 Int_t nbinsY, Float_t ymin, Float_t ymax,
188 Char_t *type, Bool_t force)
189{
190 //
191 // return pointer to Q histogram
192 // if force is true create a new histogram if it doesn't exist allready
193 //
194 if ( !force || arr->UncheckedAt(sector) )
195 return (TH2S*)arr->UncheckedAt(sector);
196
197 // if we are forced and histogram doesn't yes exist create it
198 Char_t name[255], title[255];
199
200 sprintf(name,"hCalib%s%.2d",type,sector);
201 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
202
203 // new histogram with Q calib information. One value for each pad!
204 TH2S* hist = new TH2S(name,title,
205 nbinsY, ymin, ymax,
206 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector)
207 );
208 hist->SetDirectory(0);
209 arr->AddAt(hist,sector);
210 return hist;
211}
212//_____________________________________________________________________
213TH2S* AliTPCCalibPedestal::GetHistoPedestal(Int_t sector, Bool_t force) /*FOLD00*/
214{
215 //
216 // return pointer to T0 histogram
217 // if force is true create a new histogram if it doesn't exist allready
218 //
219 TObjArray *arr = &fHistoPedestalArray;
220 return GetHisto(sector, arr, fAdcMax-fAdcMin, fAdcMin, fAdcMax, "Pedestal", force);
221}
222//_____________________________________________________________________
223AliTPCCalROC* AliTPCCalibPedestal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
224{
225 //
226 // return pointer to ROC Calibration
227 // if force is true create a new histogram if it doesn't exist allready
228 //
229 if ( !force || arr->UncheckedAt(sector) )
230 return (AliTPCCalROC*)arr->UncheckedAt(sector);
231
232 // if we are forced and histogram doesn't yes exist create it
233
234 // new AliTPCCalROC for T0 information. One value for each pad!
235 AliTPCCalROC *croc = new AliTPCCalROC(sector);
236 //init values
237 for ( UInt_t iChannel = 0; iChannel<croc->GetNchannels(); iChannel++){
238 croc->SetValue(iChannel, 0);
239 }
240 arr->AddAt(croc,sector);
241 return croc;
242}
243//_____________________________________________________________________
244AliTPCCalROC* AliTPCCalibPedestal::GetCalRocPedestal(Int_t sector, Bool_t force) /*FOLD00*/
245{
246 //
247 // return pointer to Carge ROC Calibration
248 // if force is true create a new histogram if it doesn't exist allready
249 //
250 TObjArray *arr = &fCalRocArrayPedestal;
251 return GetCalRoc(sector, arr, force);
252}
253//_____________________________________________________________________
254AliTPCCalROC* AliTPCCalibPedestal::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
255{
256 //
257 // return pointer to signal width ROC Calibration
258 // if force is true create a new histogram if it doesn't exist allready
259 //
260 TObjArray *arr = &fCalRocArrayRMS;
261 return GetCalRoc(sector, arr, force);
262}
263//_____________________________________________________________________
264void AliTPCCalibPedestal::Analyse() /*FOLD00*/
265{
266 //
267 // Calculate calibration constants
268 //
269
270 Int_t nbinsAdc = fAdcMax-fAdcMin;
271
272 TH1F *py = new TH1F("htemp_py","htemp_py", nbinsAdc, fAdcMin, fAdcMax);
273 TF1 *gaus = new TF1("fit","gaus");
274 TVectorD param(3);
275
276 Float_t *array_py=0;
277 Short_t *array_hP=0;
278
279 array_py = py->GetArray();
280
281 for (Int_t iSec=0; iSec<72; iSec++){
282 TH2S *hP = GetHistoPedestal(iSec);
283 if ( !hP ) continue;
284
285 AliTPCCalROC *rocPedestal = GetCalRocPedestal(iSec,kTRUE);
286 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
287
288 array_hP = hP->GetArray();
289 UInt_t nChannels = fROC->GetNChannels(iSec);
290
291 for (UInt_t iChannel=0; iChannel<nChannels; iChannel++){
292
293 // set bin content of py in a dirty but fast way
294 for (Int_t iAdc=0;iAdc<nbinsAdc;iAdc++)
295 array_py[iAdc+1] = (Float_t)array_hP[(iChannel+1)*(nbinsAdc+2)+(iAdc+1)];
296
297 gaus->SetParameters(0,0);
298 py->Fit(gaus,"nq");
299
300 rocPedestal->SetValue(iChannel,gaus->GetParameter(1));
301 rocRMS->SetValue(iChannel,gaus->GetParameter(2));
302
303 //AliMathBase::FitGaus(nbinsAdc,array_hP+(iChannel+1)*(nbinsAdc+2),&param)
304 //rocPedestal->SetValue(iChannel,param[1]);
305 //rocRMS->SetValue(iChannel,param[2]);
306 }
307 }
308 delete py;
309 delete gaus;
310 delete fDebugStreamer;
311 fDebugStreamer = 0x0;
312}
313//_____________________________________________________________________
314void AliTPCCalibPedestal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
315{
316 //
317 // Write class to file
318 //
319
320 TDirectory *backup = gDirectory;
321 TString sDir(dir);
322 TString option;
323
324 if ( append )
325 option = "update";
326 else
327 option = "recreate";
328
329 TFile f(filename,option.Data());
330 if ( !sDir.IsNull() ){
331 f.mkdir(sDir.Data());
332 f.cd(sDir);
333 }
334 gDirectory->WriteTObject(this);
335 f.Close();
336
337 if ( backup ) backup->cd();
338
339}