Additional include files (Matthias)
[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
8bc7e885 17/* $Id$ */
18
8bc7e885 19
20//Root includes
8bc7e885 21#include <TH1F.h>
8bc7e885 22#include <TH2F.h>
8bc7e885 23#include <TString.h>
8bc7e885 24#include <TMath.h>
25#include <TF1.h>
26#include <TRandom.h>
8bc7e885 27#include <TDirectory.h>
8bc7e885 28#include <TFile.h>
8bc7e885 29//AliRoot includes
30#include "AliRawReader.h"
31#include "AliRawReaderRoot.h"
bc331d5b 32#include "AliRawReaderDate.h"
8bc7e885 33#include "AliTPCRawStream.h"
34#include "AliTPCCalROC.h"
8bc7e885 35#include "AliTPCROC.h"
8bc7e885 36#include "AliMathBase.h"
8bc7e885 37#include "TTreeStream.h"
38
bc331d5b 39//date
40#include "event.h"
41
42//header file
43#include "AliTPCCalibPedestal.h"
8bc7e885 44
45
aa983e4f 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");
8bc7e885 177
178
aa983e4f 179*/
8bc7e885 180
181
182
aa983e4f 183ClassImp(AliTPCCalibPedestal)
8bc7e885 184
185AliTPCCalibPedestal::AliTPCCalibPedestal() : /*FOLD00*/
186 TObject(),
187 fFirstTimeBin(60),
188 fLastTimeBin(1000),
189 fAdcMin(1),
190 fAdcMax(100),
aa983e4f 191 fOldRCUformat(kTRUE),
8bc7e885 192 fROC(AliTPCROC::Instance()),
193 fCalRocArrayPedestal(72),
194 fCalRocArrayRMS(72),
bc331d5b 195 fHistoPedestalArray(72)
8bc7e885 196{
197 //
bc331d5b 198 // default constructor
8bc7e885 199 //
8bc7e885 200}
bc331d5b 201//_____________________________________________________________________
202AliTPCCalibPedestal::AliTPCCalibPedestal(const AliTPCCalibPedestal &ped) : /*FOLD00*/
203 TObject(ped),
204 fFirstTimeBin(ped.GetFirstTimeBin()),
205 fLastTimeBin(ped.GetLastTimeBin()),
206 fAdcMin(ped.GetAdcMin()),
207 fAdcMax(ped.GetAdcMax()),
aa983e4f 208 fOldRCUformat(kTRUE),
bc331d5b 209 fROC(AliTPCROC::Instance()),
210 fCalRocArrayPedestal(72),
211 fCalRocArrayRMS(72),
212 fHistoPedestalArray(72)
213{
214 //
215 // copy constructor
216 //
aa983e4f 217 for (Int_t iSec = 0; iSec < 72; ++iSec){
bc331d5b 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//_____________________________________________________________________
233AliTPCCalibPedestal& AliTPCCalibPedestal::operator = (const AliTPCCalibPedestal &source)
234{
235 //
236 // assignment operator
237 //
238 if (&source == this) return *this;
239 new (this) AliTPCCalibPedestal(source);
8bc7e885 240
bc331d5b 241 return *this;
242}
8bc7e885 243//_____________________________________________________________________
244AliTPCCalibPedestal::~AliTPCCalibPedestal() /*FOLD00*/
245{
246 //
247 // destructor
248 //
8bc7e885 249
aa983e4f 250 fCalRocArrayPedestal.Delete();
251 fCalRocArrayRMS.Delete();
252 fHistoPedestalArray.Delete();
253 delete fROC;
254}
8bc7e885 255//_____________________________________________________________________
256Int_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 //
aa983e4f 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;
8bc7e885 269
270 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
271
bc331d5b 272 // fast filling methode.
8bc7e885 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
aa983e4f 275 Int_t bin = (iChannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
8bc7e885 276
277 GetHistoPedestal(icsector,kTRUE)->GetArray()[bin]++;
278
279 return 0;
280}
281//_____________________________________________________________________
bc331d5b 282Bool_t AliTPCCalibPedestal::ProcessEvent(AliTPCRawStream *rawStream)
8bc7e885 283{
284 //
bc331d5b 285 // Event Processing loop - AliTPCRawStream
8bc7e885 286 //
287
aa983e4f 288 rawStream->SetOldRCUFormat(fOldRCUformat);
8bc7e885 289
290 Bool_t withInput = kFALSE;
291
bc331d5b 292 while (rawStream->Next()) {
8bc7e885 293
bc331d5b 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
8bc7e885 299
300 Update(isector,iRow,iPad,iTimeBin,signal);
301 withInput = kTRUE;
302 }
303
8bc7e885 304 return withInput;
305}
306//_____________________________________________________________________
bc331d5b 307Bool_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//_____________________________________________________________________
321Bool_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//_____________________________________________________________________
8bc7e885 332Bool_t AliTPCCalibPedestal::TestEvent() /*FOLD00*/
333{
334 //
335 // Test event loop
bc331d5b 336 // fill one oroc and one iroc with random gaus
8bc7e885 337 //
338
339 gRandom->SetSeed(0);
340
aa983e4f 341 for (UInt_t iSec=0; iSec<72; ++iSec){
8bc7e885 342 if (iSec%36>0) continue;
aa983e4f 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){
bc331d5b 346 Float_t signal=(Int_t)(iRow+3+gRandom->Gaus(0,.7));
8bc7e885 347 if ( signal>0 )Update(iSec,iRow,iPad,iTimeBin,signal);
348 }
349 }
350 }
351 }
352 return kTRUE;
353}
354//_____________________________________________________________________
bc331d5b 355TH2F* AliTPCCalibPedestal::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
8bc7e885 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) )
bc331d5b 364 return (TH2F*)arr->UncheckedAt(sector);
8bc7e885 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);
bc331d5b 370 sprintf(title,"%s calibration histogram sector %.2d;ADC channel;Channel (pad)",type,sector);
8bc7e885 371
372 // new histogram with Q calib information. One value for each pad!
bc331d5b 373 TH2F* hist = new TH2F(name,title,
8bc7e885 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//_____________________________________________________________________
bc331d5b 382TH2F* AliTPCCalibPedestal::GetHistoPedestal(Int_t sector, Bool_t force) /*FOLD00*/
8bc7e885 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//_____________________________________________________________________
392AliTPCCalROC* 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
aa983e4f 401 // if we are forced and the histogram doesn't yet exist create it
8bc7e885 402
403 // new AliTPCCalROC for T0 information. One value for each pad!
404 AliTPCCalROC *croc = new AliTPCCalROC(sector);
8bc7e885 405 arr->AddAt(croc,sector);
406 return croc;
407}
408//_____________________________________________________________________
409AliTPCCalROC* 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//_____________________________________________________________________
419AliTPCCalROC* 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//_____________________________________________________________________
aa983e4f 429void 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//_____________________________________________________________________
8bc7e885 454void AliTPCCalibPedestal::Analyse() /*FOLD00*/
455{
456 //
457 // Calculate calibration constants
458 //
459
460 Int_t nbinsAdc = fAdcMax-fAdcMin;
461
8bc7e885 462 TVectorD param(3);
bc331d5b 463 TMatrixD dummy(3,3);
8bc7e885 464
bc331d5b 465 Float_t *array_hP=0;
8bc7e885 466
8bc7e885 467
aa983e4f 468 for (Int_t iSec=0; iSec<72; ++iSec){
bc331d5b 469 TH2F *hP = GetHistoPedestal(iSec);
8bc7e885 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
aa983e4f 478 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
bc331d5b 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]);
8bc7e885 488 }
489 }
8bc7e885 490}
491//_____________________________________________________________________
492void AliTPCCalibPedestal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
493{
494 //
495 // Write class to file
496 //
497
8bc7e885 498 TString sDir(dir);
499 TString option;
500
501 if ( append )
502 option = "update";
503 else
504 option = "recreate";
505
bc331d5b 506 TDirectory *backup = gDirectory;
8bc7e885 507 TFile f(filename,option.Data());
bc331d5b 508 f.cd();
8bc7e885 509 if ( !sDir.IsNull() ){
510 f.mkdir(sDir.Data());
511 f.cd(sDir);
512 }
bc331d5b 513 this->Write();
8bc7e885 514 f.Close();
515
516 if ( backup ) backup->cd();
8bc7e885 517}