1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
27 #include <TDirectory.h>
31 #include "AliRawReader.h"
32 #include "AliRawReaderRoot.h"
33 #include "AliRawReaderDate.h"
34 #include "AliTPCCalROC.h"
35 #include "AliTPCROC.h"
36 #include "AliMathBase.h"
37 #include "TTreeStream.h"
43 #include "AliTPCCalibPedestal.h"
46 ///////////////////////////////////////////////////////////////////////////////////////
47 // Implementation of the TPC pedestal and noise calibration
49 // Origin: Jens Wiechula, Marian Ivanov J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
52 // *************************************************************************************
53 // * Class Description *
54 // *************************************************************************************
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
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.
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.
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;
80 // User interface for filling data:
81 // --------------------------------
83 // To Fill information one of the following functions can be used:
85 // Bool_t ProcessEvent(eventHeaderStruct *event);
86 // - process Date event
87 // - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
89 // Bool_t ProcessEvent(AliRawReader *rawReader);
90 // - process AliRawReader event
91 // - use AliTPCRawStreamV3 to loop over data and call ProcessEvent(AliTPCRawStreamV3 *rawStream)
93 // Bool_t ProcessEvent(AliTPCRawStreamV3 *rawStream);
94 // - process event from AliTPCRawStreamV3
95 // - call Update function for signal filling
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
102 // It is also possible to merge two independently taken calibrations using the function
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!
111 // -- example: filling data using root raw data:
112 // void fillPedestal(Char_t *filename)
114 // rawReader = new AliRawReaderRoot(fileName);
115 // if ( !rawReader ) return;
116 // AliTPCCalibPedestal *calib = new AliTPCCalibPedestal;
117 // while (rawReader->NextEvent()){
118 // calib->ProcessEvent(rawReader);
121 // calib->DumpToFile("PedestalData.root");
127 // What kind of information is stored and how to retrieve them:
128 // ------------------------------------------------------------
130 // - Accessing the 'Reference Histograms' (pedestal distribution histograms):
132 // TH2F *GetHistoPedestal(Int_t sector);
134 // - Accessing the calibration storage objects:
136 // AliTPCCalROC *GetCalRocPedestal(Int_t sector); - for the pedestal values, mean from gaus fit
137 // AliTPCCalROC *GetCalRocSigma(Int_t sector); - for the Noise values, sigma from guas fit
138 // AliTPCCalROC *GetCalRocMean(Int_t sector); - for the pedestal values, truncated mean
139 // AliTPCCalROC *GetCalRocRMS(Int_t sector); - for the Noise values, rms from truncated mean
141 // example for visualisation:
142 // if the file "PedestalData.root" was created using the above example one could do the following:
144 // TFile filePedestal("PedestalData.root")
145 // AliTPCCalibPedestal *ped = (AliTPCCalibPedestal*)filePedestal->Get("AliTPCCalibPedestal");
146 // ped->GetCalRocPedestal(0)->Draw("colz");
147 // ped->GetCalRocRMS(0)->Draw("colz");
149 // or use the AliTPCCalPad functionality:
150 // AliTPCCalPad padPedestal(ped->GetCalPadPedestal());
151 // AliTPCCalPad padNoise(ped->GetCalPadRMS());
152 // padPedestal->MakeHisto2D()->Draw("colz"); //Draw A-Side Pedestal Information
153 // padNoise->MakeHisto2D()->Draw("colz"); //Draw A-Side Noise Information
156 example: fill pedestal with gausschen noise
157 AliTPCCalibPedestal ped;
161 TCanvas* c1 = new TCanvas;
164 ped.GetHistoPedestal(0)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
165 ped.GetHistoPedestal(0)->Draw("colz");
167 ped.GetHistoPedestal(36)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
168 ped.GetHistoPedestal(36)->Draw("colz");
169 TCanvas* c2 = new TCanvas;
172 ped.GetCalRocPedestal(0)->Draw("colz");
174 ped.GetCalRocRMS(0)->Draw("colz");
176 ped.GetCalRocPedestal(36)->Draw("colz");
178 ped.GetCalRocRMS(36)->Draw("colz");
181 // Time dependent pedestals:
183 // If wished there is the possibility to calculate for each channel and time bin
184 // the mean pedestal [pedestals(t)]. This is done by
186 // 1) setting SetTimeAnalysis(kTRUE),
187 // 2) processing the data by looping over the events using ProcessEvent(..)
188 // 3) calling the Analyse() and AnalyseTime(nevents) functions (providing nevents)
189 // 4) getting the pedestals(t) using TArrayF **timePed = calibPedestal.GetTimePedestals();
190 // 5) looking at values using timePed[row][pad].At(timebin)
192 // This functionality is intended to be used on an LDC bu the detector algorithm
193 // (TPCPEDESTALda) to generate a data set used for configuration of the pattern
194 // memory for baseline subtraction in the ALTROs. Later the information should also
195 // be stored as reference data.
199 ClassImp(AliTPCCalibPedestal)
201 AliTPCCalibPedestal::AliTPCCalibPedestal() :
202 AliTPCCalibRawBase(),
207 fTimeAnalysis(kFALSE),
208 fCalRocArrayPedestal(72),
209 fCalRocArraySigma(72),
210 fHistoPedestalArray(72),
212 fCalRocArrayMean(72),
216 // default constructor
218 SetNameTitle("AliTPCCalibPedestal","AliTPCCalibPedestal");
224 //_____________________________________________________________________
225 AliTPCCalibPedestal::AliTPCCalibPedestal(const AliTPCCalibPedestal &ped) :
226 AliTPCCalibRawBase(ped),
227 fAdcMin(ped.GetAdcMin()),
228 fAdcMax(ped.GetAdcMax()),
229 fAnaMeanDown(ped.fAnaMeanDown),
230 fAnaMeanUp(ped.fAnaMeanUp),
231 fTimeAnalysis(ped.fTimeAnalysis),
232 fCalRocArrayPedestal(72),
233 fCalRocArraySigma(72),
234 fHistoPedestalArray(72),
235 fTimeSignal(ped.fTimeSignal),
236 fCalRocArrayMean(72),
242 for (Int_t iSec = 0; iSec < 72; ++iSec){
243 const AliTPCCalROC *calPed = (AliTPCCalROC*)ped.fCalRocArrayPedestal.UncheckedAt(iSec);
244 const AliTPCCalROC *calRMS = (AliTPCCalROC*)ped.fCalRocArrayRMS.UncheckedAt(iSec);
245 const TH2F *hPed = (TH2F*)ped.fHistoPedestalArray.UncheckedAt(iSec);
247 if ( calPed != 0x0 ) fCalRocArrayPedestal.AddAt(new AliTPCCalROC(*calPed), iSec);
248 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
251 TH2F *hNew = new TH2F(*hPed);
252 hNew->SetDirectory(0);
253 fHistoPedestalArray.AddAt(hNew,iSec);
257 AliTPCCalibPedestal::AliTPCCalibPedestal(const TMap *config):
258 AliTPCCalibRawBase(),
263 fTimeAnalysis(kFALSE),
264 fCalRocArrayPedestal(72),
265 fCalRocArraySigma(72),
266 fHistoPedestalArray(72),
268 fCalRocArrayMean(72),
272 // This constructor uses a TMap for setting some parametes
274 SetNameTitle("AliTPCCalibPedestal","AliTPCCalibPedestal");
277 if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
278 if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
279 if (config->GetValue("AdcMin")) fAdcMin = ((TObjString*)config->GetValue("AdcMin"))->GetString().Atoi();
280 if (config->GetValue("AdcMax")) fAdcMax = ((TObjString*)config->GetValue("AdcMax"))->GetString().Atoi();
281 if (config->GetValue("TimeAnalysis")) SetTimeAnalysis(((TObjString*)config->GetValue("TimeAnalysis"))->GetString().Atoi());
285 //_____________________________________________________________________
286 AliTPCCalibPedestal& AliTPCCalibPedestal::operator = (const AliTPCCalibPedestal &source)
289 // assignment operator
291 if (&source == this) return *this;
292 new (this) AliTPCCalibPedestal(source);
298 //_____________________________________________________________________
299 AliTPCCalibPedestal::~AliTPCCalibPedestal()
305 fCalRocArrayPedestal.Delete();
306 fCalRocArrayRMS.Delete();
307 fCalRocArraySigma.Delete();
308 fHistoPedestalArray.Delete();
311 for (Int_t i = 0; i < 159; i++) {
312 delete [] fTimeSignal[i];
315 delete [] fTimeSignal;
319 // do not delete fMapping, because we do not own it.
324 //_____________________________________________________________________
325 void AliTPCCalibPedestal::SetTimeAnalysis(Bool_t time)
328 // Use time dependent analysis: Pedestals are analysed as a function
329 // of the drift time. There is one mean value generated for each time
330 // bin and each channel. It can be used as reference data and for
331 // configuration of the ALTRO pattern memory for baseline subtraction.
333 // ATTENTION: Use only on LDC in TPCPEDESTALda! On a LDC we get data
334 // only from one sector. For the full TPC we would need a lot of
335 // memory (36*159*140*1024*4bytes = 3.3GB)!
338 fTimeAnalysis = time;
340 if ( !fTimeAnalysis ) return;
342 // prepare array for one sector (159*140*1024*4bytes = 92MB):
343 fTimeSignal = new TArrayF*[159];
344 for (Int_t i = 0; i < 159; i++) { // padrows
345 fTimeSignal[i] = new TArrayF[140];
346 for (Int_t j = 0; j < 140; j++) { // pads per row
347 fTimeSignal[i][j].Set(1024);
348 for (Int_t k = 0; k < 1024; k++) { // time bins per pad
349 fTimeSignal[i][j].AddAt(0., k);
356 //_____________________________________________________________________
357 Int_t AliTPCCalibPedestal::Update(const Int_t icsector,
360 const Int_t icTimeBin,
361 const Float_t csignal)
364 // Signal filling method
366 if (icRow<0) return 0;
367 if (icPad<0) return 0;
368 if (icTimeBin<0) return 0;
370 // Time dependent pedestals
371 if ( fTimeAnalysis ) {
372 if ( icsector < 36 ) // IROC
373 fTimeSignal[icRow][icPad].AddAt(fTimeSignal[icRow][icPad].At(icTimeBin)+csignal, icTimeBin);
375 fTimeSignal[icRow+63][icPad].AddAt(fTimeSignal[icRow+63][icPad].At(icTimeBin)+csignal, icTimeBin);
377 //return if we are out of the specified time bin or adc range
378 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
379 if ( ((Int_t)csignal>fAdcMax) || ((Int_t)csignal<fAdcMin) ) return 0;
381 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
383 // fast filling method
384 // Attention: the entry counter of the histogram is not increased
385 // this means that e.g. the colz draw option gives an empty plot
386 Int_t bin = (iChannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
388 GetHistoPedestal(icsector,kTRUE)->GetArray()[bin]++;
394 //_____________________________________________________________________
395 Bool_t AliTPCCalibPedestal::TestEvent()
399 // fill one oroc and one iroc with random gaus
404 for (UInt_t iSec=0; iSec<72; ++iSec){
405 if (iSec%36>0) continue;
406 for (UInt_t iRow=0; iRow < fROC->GetNRows(iSec); ++iRow){
407 for (UInt_t iPad=0; iPad < fROC->GetNPads(iSec,iRow); ++iPad){
408 for (UInt_t iTimeBin=0; iTimeBin<1024; ++iTimeBin){
409 Float_t signal=(Int_t)(iRow+3+gRandom->Gaus(0,.7));
410 if ( signal>0 )Update(iSec,iRow,iPad,iTimeBin,signal);
419 //_____________________________________________________________________
420 TH2F* AliTPCCalibPedestal::GetHisto(Int_t sector, TObjArray *arr,
421 Int_t nbinsY, Float_t ymin, Float_t ymax,
422 const Char_t *type, Bool_t force)
425 // return pointer to Q histogram
426 // if force is true create a new histogram if it doesn't exist allready
428 if ( !force || arr->UncheckedAt(sector) )
429 return (TH2F*)arr->UncheckedAt(sector);
431 // if we are forced and histogram doesn't yes exist create it
432 // new histogram with Q calib information. One value for each pad!
433 TH2F* hist = new TH2F(Form("hCalib%s%.2d",type,sector),
434 Form("%s calibration histogram sector %.2d;ADC channel;Channel (pad)",type,sector),
436 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector)
438 hist->SetDirectory(0);
439 arr->AddAt(hist,sector);
444 //_____________________________________________________________________
445 TH2F* AliTPCCalibPedestal::GetHistoPedestal(Int_t sector, Bool_t force)
448 // return pointer to T0 histogram
449 // if force is true create a new histogram if it doesn't exist allready
451 TObjArray *arr = &fHistoPedestalArray;
452 return GetHisto(sector, arr, fAdcMax-fAdcMin, fAdcMin, fAdcMax, "Pedestal", force);
456 //_____________________________________________________________________
457 AliTPCCalROC* AliTPCCalibPedestal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force)
460 // return pointer to ROC Calibration
461 // if force is true create a new histogram if it doesn't exist allready
463 if ( !force || arr->UncheckedAt(sector) )
464 return (AliTPCCalROC*)arr->UncheckedAt(sector);
466 // if we are forced and the histogram doesn't yet exist create it
468 // new AliTPCCalROC for T0 information. One value for each pad!
469 AliTPCCalROC *croc = new AliTPCCalROC(sector);
470 arr->AddAt(croc,sector);
475 //_____________________________________________________________________
476 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocPedestal(Int_t sector, Bool_t force)
479 // return pointer to ROC with Pedestal data
480 // if force is true create a new histogram if it doesn't exist allready
482 TObjArray *arr = &fCalRocArrayPedestal;
483 return GetCalRoc(sector, arr, force);
487 //_____________________________________________________________________
488 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocSigma(Int_t sector, Bool_t force)
491 // return pointer to ROC with signal witdth in sigma
492 // if force is true create a new histogram if it doesn't exist allready
494 TObjArray *arr = &fCalRocArraySigma;
495 return GetCalRoc(sector, arr, force);
497 //_____________________________________________________________________
498 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocMean(Int_t sector, Bool_t force)
501 // return pointer to ROC with signal mean information
502 // if force is true create a new histogram if it doesn't exist allready
504 TObjArray *arr = &fCalRocArrayMean;
505 return GetCalRoc(sector, arr, force);
508 //_____________________________________________________________________
509 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocRMS(Int_t sector, Bool_t force)
512 // return pointer to signal width ROC Calibration
513 // if force is true create a new histogram if it doesn't exist allready
515 TObjArray *arr = &fCalRocArrayRMS;
516 return GetCalRoc(sector, arr, force);
520 //_____________________________________________________________________
521 void AliTPCCalibPedestal::Merge(AliTPCCalibPedestal * const ped)
524 // Merge reference histograms of sig to the current AliTPCCalibPedestal
528 for (Int_t iSec=0; iSec<72; ++iSec){
529 TH2F *hRefPedMerge = ped->GetHistoPedestal(iSec);
532 TDirectory *dir = hRefPedMerge->GetDirectory(); hRefPedMerge->SetDirectory(0);
533 TH2F *hRefPed = GetHistoPedestal(iSec);
534 if ( hRefPed ) hRefPed->Add(hRefPedMerge);
536 TH2F *hist = new TH2F(*hRefPedMerge);
537 hist->SetDirectory(0);
538 fHistoPedestalArray.AddAt(hist, iSec);
540 hRefPedMerge->SetDirectory(dir);
549 //_____________________________________________________________________
550 Long64_t AliTPCCalibPedestal::Merge(TCollection * const list)
553 // Merge all objects of this type in list
559 AliTPCCalibPedestal *ce=0;
562 while ( (o=next()) ){
563 ce=dynamic_cast<AliTPCCalibPedestal*>(o);
573 //_____________________________________________________________________
574 void AliTPCCalibPedestal::Analyse()
577 // Calculate calibration constants
580 Int_t nbinsAdc = fAdcMax-fAdcMin;
585 TH1F *hChannel=new TH1F("hChannel","hChannel",nbinsAdc,fAdcMin,fAdcMax);
589 for (Int_t iSec=0; iSec<72; ++iSec){
590 TH2F *hP = GetHistoPedestal(iSec);
593 AliTPCCalROC *rocPedestal = GetCalRocPedestal(iSec,kTRUE);
594 AliTPCCalROC *rocSigma = GetCalRocSigma(iSec,kTRUE);
595 AliTPCCalROC *rocMean = GetCalRocMean(iSec,kTRUE);
596 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
598 arrayhP = hP->GetArray();
599 UInt_t nChannels = fROC->GetNChannels(iSec);
601 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
602 Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
603 //calculate mean and sigma using a gaus fit
605 AliMathBase::FitGaus(arrayhP+offset,nbinsAdc,fAdcMin,fAdcMax,¶m,&dummy);
606 // if the fitting failed set noise and pedestal to 0
607 // is now done in AliMathBase::FitGaus !
608 // if ( ret == -4 ) {
612 if ( param[1]<fAdcMin || param[1]>fAdcMax ){
616 rocPedestal->SetValue(iChannel,param[1]);
617 rocSigma->SetValue(iChannel,param[2]);
618 //calculate mean and RMS using a truncated means
619 hChannel->Set(nbinsAdc+2,arrayhP+offset-1);
620 hChannel->SetEntries(param[3]);
623 if ( param[3]>0 ) AliMathBase::TruncatedMean(hChannel,¶m,fAnaMeanDown,fAnaMeanUp);
624 rocMean->SetValue(iChannel,param[1]);
625 rocRMS->SetValue(iChannel,param[2]);
632 //_____________________________________________________________________
633 void AliTPCCalibPedestal::AnalyseTime(Int_t nevents)
636 // Calculate for each channel and time bin the mean pedestal. This
637 // is used on LDC by TPCPEDESTALda to generate data used for configuration
638 // of the pattern memory for baseline subtraction in the ALTROs.
641 if ( nevents <= 0 ) return;
642 if ( fTimeAnalysis ) {
643 for (Int_t i = 0; i < 159; i++) { // padrows
644 for (Int_t j = 0; j < 140; j++) { // pads per row
645 for (Int_t k = 0; k < 1024; k++) { // time bins per pad
646 fTimeSignal[i][j].AddAt(fTimeSignal[i][j].At(k)/(Float_t)nevents, k);