[u/mrichter/AliRoot.git] / TPC / AliTPCCalibPedestal.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


//-------------------------------------------------------
//          Implementation of the TPC pedestal and noise calibration
//
//   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
// 
// 
//-------------------------------------------------------


/* $Id$ */


//Root includes
#include <TObjArray.h>
#include <TH1F.h>
#include <TH1D.h>
#include <TH2F.h>
#include <TH2S.h>
#include <TH1S.h>
#include <TString.h>
#include <TVectorF.h>
#include <TMath.h>
#include <TF1.h>
#include <TRandom.h>
#include <TROOT.h>
#include <TDirectory.h>
#include <TSystem.h>
#include <TFile.h>


//AliRoot includes
#include "AliRawReader.h"
#include "AliRawReaderRoot.h"
#include "AliTPCRawStream.h"
#include "AliTPCCalROC.h"
#include "AliTPCCalPad.h"
#include "AliTPCROC.h"
#include "AliTPCCalibPedestal.h"
#include "AliMathBase.h"

#include "TTreeStream.h"


ClassImp(AliTPCCalibPedestal) /*FOLD00*/


AliTPCCalibPedestal::AliTPCCalibPedestal() : /*FOLD00*/
  TObject(),
  fFirstTimeBin(60),
  fLastTimeBin(1000),
  fAdcMin(1),
  fAdcMax(100),
  fROC(AliTPCROC::Instance()),
  fCalRocArrayPedestal(72),
  fCalRocArrayRMS(72),
  fHistoPedestalArray(72),
  fDebugStreamer(0),
  fDebugLevel(0)
{
    //
    // AliTPCSignal default constructor
    //
    //debug stream
    TDirectory *backup = gDirectory;
    fDebugStreamer = new TTreeSRedirector("deb2.root");
    if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer    
}

//_____________________________________________________________________
AliTPCCalibPedestal::~AliTPCCalibPedestal() /*FOLD00*/
{
  //
  // destructor
  //
    if ( fDebugStreamer ) delete fDebugStreamer;
    delete fROC;
}


//_____________________________________________________________________
Int_t AliTPCCalibPedestal::Update(const Int_t icsector, /*FOLD00*/
				const Int_t icRow,
				const Int_t icPad,
				const Int_t icTimeBin,
				const Float_t csignal)
{
    //
    // Signal filling methode 
    //
    if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;

    Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector

    // dirty and fast filling methode. No boundary checking!!!
    // Attention: the entry counter of the histogram is not increased
    //            this means that e.g. the colz draw option gives an empty plot
    Int_t bin = (iChannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);

    GetHistoPedestal(icsector,kTRUE)->GetArray()[bin]++;

    return 0;
}
//_____________________________________________________________________
Bool_t AliTPCCalibPedestal::ProcessEvent(AliRawReader *rawReader) /*FOLD00*/
{
  //
  //  simple event processing loop
  //


    AliTPCRawStream input(rawReader);

  rawReader->Select("TPC");

  input.SetOldRCUFormat(1);
  printf("start event processing\n");


  Bool_t withInput = kFALSE;

  while (input.Next()) {

      Int_t isector  = input.GetSector();                       //  current sector
      Int_t iRow     = input.GetRow();                          //  current row
      Int_t iPad     = input.GetPad();                          //  current pad
      Int_t iTimeBin = input.GetTime();                         //  current time bin
      Float_t signal = input.GetSignal();                       //  current ADC signal

      Update(isector,iRow,iPad,iTimeBin,signal);
      withInput = kTRUE;
  }

  printf("end event processing\n");
  if ( fDebugLevel>0 )
      fDebugStreamer->GetFile()->Write();
  return withInput;
}
//_____________________________________________________________________
Bool_t AliTPCCalibPedestal::TestEvent() /*FOLD00*/
{
  //
  //  Test event loop
  //

    gRandom->SetSeed(0);

    for (UInt_t iSec=0; iSec<72; iSec++){
        if (iSec%36>0) continue;
	for (UInt_t iRow=0; iRow < fROC->GetNRows(iSec); iRow++){
	    for (UInt_t iPad=0; iPad < fROC->GetNPads(iSec,iRow); iPad++){
		for (UInt_t iTimeBin=0; iTimeBin<1024; iTimeBin++){
		    Float_t signal=(Int_t)(iRow+5+gRandom->Gaus(0,.7));
		    if ( signal>0 )Update(iSec,iRow,iPad,iTimeBin,signal);
		}
	    }
	}
    }
    return kTRUE;
}
//_____________________________________________________________________
TH2S* AliTPCCalibPedestal::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
				  Int_t nbinsY, Float_t ymin, Float_t ymax,
				  Char_t *type, Bool_t force)
{
    //
    // return pointer to Q histogram
    // if force is true create a new histogram if it doesn't exist allready
    //
    if ( !force || arr->UncheckedAt(sector) )
	return (TH2S*)arr->UncheckedAt(sector);

    // if we are forced and histogram doesn't yes exist create it
    Char_t name[255], title[255];

    sprintf(name,"hCalib%s%.2d",type,sector);
    sprintf(title,"%s calibration histogram sector %.2d",type,sector);

    // new histogram with Q calib information. One value for each pad!
    TH2S* hist = new TH2S(name,title,
			  nbinsY, ymin, ymax,
			  fROC->GetNChannels(sector),0,fROC->GetNChannels(sector)
			 );
    hist->SetDirectory(0);
    arr->AddAt(hist,sector);
    return hist;
}
//_____________________________________________________________________
TH2S* AliTPCCalibPedestal::GetHistoPedestal(Int_t sector, Bool_t force) /*FOLD00*/
{
    //
    // return pointer to T0 histogram
    // if force is true create a new histogram if it doesn't exist allready
    //
    TObjArray *arr = &fHistoPedestalArray;
    return GetHisto(sector, arr, fAdcMax-fAdcMin, fAdcMin, fAdcMax, "Pedestal", force);
}
//_____________________________________________________________________
AliTPCCalROC* AliTPCCalibPedestal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
{
    //
    // return pointer to ROC Calibration
    // if force is true create a new histogram if it doesn't exist allready
    //
    if ( !force || arr->UncheckedAt(sector) )
	return (AliTPCCalROC*)arr->UncheckedAt(sector);

    // if we are forced and histogram doesn't yes exist create it

    // new AliTPCCalROC for T0 information. One value for each pad!
    AliTPCCalROC *croc = new AliTPCCalROC(sector);
    //init values
    for ( UInt_t iChannel = 0; iChannel<croc->GetNchannels(); iChannel++){
	croc->SetValue(iChannel, 0);
    }
    arr->AddAt(croc,sector);
    return croc;
}
//_____________________________________________________________________
AliTPCCalROC* AliTPCCalibPedestal::GetCalRocPedestal(Int_t sector, Bool_t force) /*FOLD00*/
{
    //
    // return pointer to Carge ROC Calibration
    // if force is true create a new histogram if it doesn't exist allready
    //
    TObjArray *arr = &fCalRocArrayPedestal;
    return GetCalRoc(sector, arr, force);
}
//_____________________________________________________________________
AliTPCCalROC* AliTPCCalibPedestal::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
{
    //
    // return pointer to signal width ROC Calibration
    // if force is true create a new histogram if it doesn't exist allready
    //
    TObjArray *arr = &fCalRocArrayRMS;
    return GetCalRoc(sector, arr, force);
}
//_____________________________________________________________________
void AliTPCCalibPedestal::Analyse() /*FOLD00*/
{
    //
    //  Calculate calibration constants
    //

    Int_t nbinsAdc = fAdcMax-fAdcMin;

    TH1F *py = new TH1F("htemp_py","htemp_py", nbinsAdc, fAdcMin, fAdcMax);
    TF1 *gaus = new TF1("fit","gaus");
    TVectorD param(3);

    Float_t *array_py=0;
    Short_t *array_hP=0;

    array_py = py->GetArray();

    for (Int_t iSec=0; iSec<72; iSec++){
	TH2S *hP = GetHistoPedestal(iSec);
        if ( !hP ) continue;

	AliTPCCalROC *rocPedestal = GetCalRocPedestal(iSec,kTRUE);
	AliTPCCalROC *rocRMS      = GetCalRocRMS(iSec,kTRUE);

	array_hP = hP->GetArray();
        UInt_t nChannels = fROC->GetNChannels(iSec);

	for (UInt_t iChannel=0; iChannel<nChannels; iChannel++){

            // set bin content of py in a dirty but fast way
	    for (Int_t iAdc=0;iAdc<nbinsAdc;iAdc++)
		array_py[iAdc+1] = (Float_t)array_hP[(iChannel+1)*(nbinsAdc+2)+(iAdc+1)];

	    gaus->SetParameters(0,0);
	    py->Fit(gaus,"nq");

	    rocPedestal->SetValue(iChannel,gaus->GetParameter(1));
            rocRMS->SetValue(iChannel,gaus->GetParameter(2));

	    //AliMathBase::FitGaus(nbinsAdc,array_hP+(iChannel+1)*(nbinsAdc+2),&param)
	    //rocPedestal->SetValue(iChannel,param[1]);
            //rocRMS->SetValue(iChannel,param[2]);
	}
    }
    delete py;
    delete gaus;
    delete fDebugStreamer;
    fDebugStreamer = 0x0;
}
//_____________________________________________________________________
void AliTPCCalibPedestal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
{
    //
    //  Write class to file
    //

    TDirectory *backup = gDirectory;
    TString sDir(dir);
    TString option;

    if ( append )
	option = "update";
    else
        option = "recreate";

    TFile f(filename,option.Data());
    if ( !sDir.IsNull() ){
	f.mkdir(sDir.Data());
	f.cd(sDir);
    }
    gDirectory->WriteTObject(this);
    f.Close();

    if ( backup ) backup->cd();

}
Commit	Line	Data
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
	62	ClassImp(AliTPCCalibPedestal) /FOLD00/
	63
	64
65
66
67
68
69
70	AliTPCCalibPedestal::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	//_____________________________________________________________________
93	AliTPCCalibPedestal::~AliTPCCalibPedestal() /FOLD00/
94	{
95	//
96	// destructor
97	//
98	if ( fDebugStreamer ) delete fDebugStreamer;
99	delete fROC;
100	}
101
102
103
104
105	//_____________________________________________________________________
106	Int_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	//_____________________________________________________________________
129	Bool_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	//_____________________________________________________________________
164	Bool_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	//_____________________________________________________________________
186	TH2S* 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	//_____________________________________________________________________
213	TH2S* 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	//_____________________________________________________________________
223	AliTPCCalROC* 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	//_____________________________________________________________________
244	AliTPCCalROC* 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	//_____________________________________________________________________
254	AliTPCCalROC* 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	//_____________________________________________________________________
264	void 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	//_____________________________________________________________________
314	void 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	}