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