[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.                  *
 **************************************************************************/


/* $Id$ */


//Root includes
#include <TObjArray.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TString.h>
#include <TMath.h>
#include <TF1.h>
#include <TRandom.h>
#include <TDirectory.h>
#include <TFile.h>
//AliRoot includes
#include "AliRawReader.h"
#include "AliRawReaderRoot.h"
#include "AliRawReaderDate.h"
#include "AliTPCRawStream.h"
#include "AliTPCCalROC.h"
#include "AliTPCROC.h"
#include "AliMathBase.h"
#include "TTreeStream.h"

//date
#include "event.h"

//header file
#include "AliTPCCalibPedestal.h"


///////////////////////////////////////////////////////////////////////////////////////
//          Implementation of the TPC pedestal and noise calibration
//
//   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
// 
// 
// *************************************************************************************
// *                                Class Description                                  *
// *************************************************************************************
//
// Working principle:
// ------------------
// Raw pedestal data is processed by calling one of the ProcessEvent(...) functions
// (see below). These in the end call the Update(...) function, where the data is filled
// into histograms.
//
// For each ROC one TH2F histo (ROC channel vs. ADC channel) is created when
// it is filled for the first time (GetHistoPedestal(ROC,kTRUE)). All histos are stored in the
// TObjArray fHistoPedestalArray.
//
// For a fast filling of the histogram the corresponding bin number of the channel and ADC channel
// is computed by hand and the histogram array is accessed directly via its pointer.
// ATTENTION: Doing so the the entry counter of the histogram is not increased
//            this means that e.g. the colz draw option gives an empty plot unless
//	    calling 'histo->SetEntries(1)' before drawing.
//
// After accumulating the desired statistics the Analyse() function has to be called.
// Whithin this function the pedestal and noise values are calculated for each pad, using
// the fast gaus fit function  AliMathBase::FitGaus(...), and the calibration
// storage classes (AliTPCCalROC) are filled for each ROC.
// The calibration information is stored in the TObjArrays fCalRocArrayPedestal and fCalRocArrayRMS;
//
//
//
// User interface for filling data:
// --------------------------------
//
// To Fill information one of the following functions can be used:
//
// Bool_t ProcessEvent(eventHeaderStruct *event);
//   - process Date event
//   - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
//
// Bool_t ProcessEvent(AliRawReader *rawReader);
//  - process AliRawReader event
//   - use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)
//
// Bool_t ProcessEvent(AliTPCRawStream *rawStream);
//   - process event from AliTPCRawStream
//   - call Update function for signal filling
//
// Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
//              iPad,  const Int_t iTimeBin, const Float_t signal);
//   - directly  fill signal information (sector, row, pad, time bin, pad)
//     to the reference histograms
//
// It is also possible to merge two independently taken calibrations using the function
//
// void Merge(AliTPCCalibPedestal *ped)
//   - copy histograms in 'ped' if the do not exist in this instance
//   - Add histograms in 'ped' to the histograms in this instance if the allready exist
//   - After merging call Analyse again!
//
//
//
// -- example: filling data using root raw data:
// void fillPedestal(Char_t *filename)
// {
//    rawReader = new AliRawReaderRoot(fileName);
//    if ( !rawReader ) return;
//    AliTPCCalibPedestal *calib = new AliTPCCalibPedestal;
//    while (rawReader->NextEvent()){
//      calib->ProcessEvent(rawReader);
//    }
//    calib->Analyse();
//    calib->DumpToFile("PedestalData.root");
//    delete rawReader;
//    delete calib;
// }
//
//
// What kind of information is stored and how to retrieve them:
// ------------------------------------------------------------
//
// - Accessing the 'Reference Histograms' (pedestal distribution histograms):
//
// TH2F *GetHistoPedestal(Int_t sector);
//
// - Accessing the calibration storage objects:
//
// AliTPCCalROC *GetCalRocPedestal(Int_t sector);  - for the pedestal values
// AliTPCCalROC *GetCalRocNoise(Int_t sector);     - for the Noise values
//
// example for visualisation:
// if the file "PedestalData.root" was created using the above example one could do the following:
//
// TFile filePedestal("PedestalData.root")
// AliTPCCalibPedestal *ped = (AliTPCCalibPedestal*)filePedestal->Get("AliTPCCalibPedestal");
// ped->GetCalRocPedestal(0)->Draw("colz");
// ped->GetCalRocRMS(0)->Draw("colz");
//
// or use the AliTPCCalPad functionality:
// AliTPCCalPad padPedestal(ped->GetCalPadPedestal());
// AliTPCCalPad padNoise(ped->GetCalPadRMS());
// padPedestal->MakeHisto2D()->Draw("colz");  //Draw A-Side Pedestal Information
// padNoise->MakeHisto2D()->Draw("colz");  //Draw A-Side Noise Information
//
/*
 example: fill pedestal with gausschen noise
 AliTPCCalibPedestal ped;
 ped.TestEvent();
 ped.Analyse();
 //Draw output;
 TCanvas* c1 = new TCanvas;
 c1->Divide(1,2);
 c1->cd(1);
 ped.GetHistoPedestal(0)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
 ped.GetHistoPedestal(0)->Draw("colz");
 c1->cd(2);
 ped.GetHistoPedestal(36)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
 ped.GetHistoPedestal(36)->Draw("colz");
 TCanvas* c2 = new TCanvas;
 c2->Divide(2,2);
 c2->cd(1);
 ped.GetCalRocPedestal(0)->Draw("colz");
 c2->cd(2);
 ped.GetCalRocRMS(0)->Draw("colz");
 c2->cd(3);
 ped.GetCalRocPedestal(36)->Draw("colz");
 c2->cd(4);
 ped.GetCalRocRMS(36)->Draw("colz");


*/


ClassImp(AliTPCCalibPedestal)

AliTPCCalibPedestal::AliTPCCalibPedestal() : /*FOLD00*/
  TObject(),
  fFirstTimeBin(60),
  fLastTimeBin(1000),
  fAdcMin(1),
  fAdcMax(100),
  fOldRCUformat(kTRUE),
  fROC(AliTPCROC::Instance()),
  fCalRocArrayPedestal(72),
  fCalRocArrayRMS(72),
  fHistoPedestalArray(72)
{
    //
    // default constructor
    //
}
//_____________________________________________________________________
AliTPCCalibPedestal::AliTPCCalibPedestal(const AliTPCCalibPedestal &ped) : /*FOLD00*/
  TObject(ped),
  fFirstTimeBin(ped.GetFirstTimeBin()),
  fLastTimeBin(ped.GetLastTimeBin()),
  fAdcMin(ped.GetAdcMin()),
  fAdcMax(ped.GetAdcMax()),
  fOldRCUformat(kTRUE),
  fROC(AliTPCROC::Instance()),
  fCalRocArrayPedestal(72),
  fCalRocArrayRMS(72),
  fHistoPedestalArray(72)
{
    //
    // copy constructor
    //
    for (Int_t iSec = 0; iSec < 72; ++iSec){
	const AliTPCCalROC *calPed = (AliTPCCalROC*)ped.fCalRocArrayPedestal.UncheckedAt(iSec);
	const AliTPCCalROC *calRMS = (AliTPCCalROC*)ped.fCalRocArrayRMS.UncheckedAt(iSec);
	const TH2F         *hPed   = (TH2F*)ped.fHistoPedestalArray.UncheckedAt(iSec);

	if ( calPed != 0x0 ) fCalRocArrayPedestal.AddAt(new AliTPCCalROC(*calPed), iSec);
	if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);

	if ( hPed != 0x0 ){
	    TH2F *hNew = new TH2F(*hPed);
	    hNew->SetDirectory(0);
	    fHistoPedestalArray.AddAt(hNew,iSec);
	}
    }
}
//_____________________________________________________________________
AliTPCCalibPedestal& AliTPCCalibPedestal::operator = (const  AliTPCCalibPedestal &source)
{
  //
  // assignment operator
  //
  if (&source == this) return *this;
  new (this) AliTPCCalibPedestal(source);

  return *this;
}
//_____________________________________________________________________
AliTPCCalibPedestal::~AliTPCCalibPedestal() /*FOLD00*/
{
  //
  // destructor
  //

    fCalRocArrayPedestal.Delete();
    fCalRocArrayRMS.Delete();
    fHistoPedestalArray.Delete();
    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 
    //

    //return if we are out of the specified time bin or adc range
    if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
    if ( ((Int_t)csignal>fAdcMax) || ((Int_t)csignal<fAdcMin)  ) return 0;

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

    // fast filling methode.
    // 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(AliTPCRawStream *rawStream)
{
  //
  // Event Processing loop - AliTPCRawStream
  //

  rawStream->SetOldRCUFormat(fOldRCUformat);

  Bool_t withInput = kFALSE;

  while (rawStream->Next()) {

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

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

  return withInput;
}
//_____________________________________________________________________
Bool_t AliTPCCalibPedestal::ProcessEvent(AliRawReader *rawReader)
{
  //
  //  Event processing loop - AliRawReader
  //


  AliTPCRawStream rawStream(rawReader);

  rawReader->Select("TPC");

  return ProcessEvent(&rawStream);
}
//_____________________________________________________________________
Bool_t AliTPCCalibPedestal::ProcessEvent(eventHeaderStruct *event)
{
  //
  //  process date event
  //
    AliRawReader *rawReader = new AliRawReaderDate((void*)event);
    Bool_t result=ProcessEvent(rawReader);
    delete rawReader;
    return result;
}
//_____________________________________________________________________
Bool_t AliTPCCalibPedestal::TestEvent() /*FOLD00*/
{
  //
  //  Test event loop
  // fill one oroc and one iroc with random gaus
  //

    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+3+gRandom->Gaus(0,.7));
		    if ( signal>0 )Update(iSec,iRow,iPad,iTimeBin,signal);
		}
	    }
	}
    }
    return kTRUE;
}
//_____________________________________________________________________
TH2F* 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 (TH2F*)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;ADC channel;Channel (pad)",type,sector);

    // new histogram with Q calib information. One value for each pad!
    TH2F* hist = new TH2F(name,title,
			  nbinsY, ymin, ymax,
			  fROC->GetNChannels(sector),0,fROC->GetNChannels(sector)
			 );
    hist->SetDirectory(0);
    arr->AddAt(hist,sector);
    return hist;
}
//_____________________________________________________________________
TH2F* 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 the histogram doesn't yet exist create it

    // new AliTPCCalROC for T0 information. One value for each pad!
    AliTPCCalROC *croc = new AliTPCCalROC(sector);
    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::Merge(AliTPCCalibPedestal *ped)
{
    //
    //  Merge reference histograms of sig to the current AliTPCCalibSignal
    //

    //merge histograms
    for (Int_t iSec=0; iSec<72; ++iSec){
	TH2F *hRefPedMerge   = ped->GetHistoPedestal(iSec);


	if ( hRefPedMerge ){
	    TDirectory *dir = hRefPedMerge->GetDirectory(); hRefPedMerge->SetDirectory(0);
	    TH2F *hRefPed   = GetHistoPedestal(iSec);
	    if ( hRefPed ) hRefPed->Add(hRefPedMerge);
	    else {
		TH2F *hist = new TH2F(*hRefPedMerge);
                hist->SetDirectory(0);
		fHistoPedestalArray.AddAt(hist, iSec);
	    }
	    hRefPedMerge->SetDirectory(dir);
	}
    }
}
//_____________________________________________________________________
void AliTPCCalibPedestal::Analyse() /*FOLD00*/
{
    //
    //  Calculate calibration constants
    //

    Int_t nbinsAdc = fAdcMax-fAdcMin;

    TVectorD param(3);
    TMatrixD dummy(3,3);

    Float_t *array_hP=0;


    for (Int_t iSec=0; iSec<72; ++iSec){
	TH2F *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){
            Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
	    Double_t ret = AliMathBase::FitGaus(array_hP+offset,nbinsAdc,fAdcMin,fAdcMax,&param,&dummy);
            // if the fitting failed set noise and pedestal to 0
	    if ( ret == -4 ) {
		param[1]=0;
		param[2]=0;
	    }
	    rocPedestal->SetValue(iChannel,param[1]);
            rocRMS->SetValue(iChannel,param[2]);
	}
    }
}
//_____________________________________________________________________
void AliTPCCalibPedestal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
{
    //
    //  Write class to file
    //

    TString sDir(dir);
    TString option;

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

    TDirectory *backup = gDirectory;
    TFile f(filename,option.Data());
    f.cd();
    if ( !sDir.IsNull() ){
	f.mkdir(sDir.Data());
	f.cd(sDir);
    }
    this->Write();
    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
8bc7e885	17	/* $Id$ */
8bc7e885	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"
8bc7e885	35	#include "AliTPCCalROC.h"
8bc7e885	36	#include "AliTPCROC.h"
8bc7e885	37	#include "AliMathBase.h"
8bc7e885	38	#include "TTreeStream.h"
8bc7e885	39
bc331d5b	40	//date
	41	#include "event.h"
	42
	43	//header file
	44	#include "AliTPCCalibPedestal.h"
8bc7e885	45
8bc7e885	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
8bc7e885	179
aa983e4f	180	*/
8bc7e885	181
	182
	183
aa983e4f	184	ClassImp(AliTPCCalibPedestal)
8bc7e885	185
	186	AliTPCCalibPedestal::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	{
8bc7e885	198	//
bc331d5b	199	// default constructor
8bc7e885	200	//
8bc7e885	201	}
bc331d5b	202	//_____________________________________________________________________
	203	AliTPCCalibPedestal::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	//_____________________________________________________________________
	234	AliTPCCalibPedestal& 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;
bc331d5b	243	}
8bc7e885	244	//_____________________________________________________________________
	245	AliTPCCalibPedestal::~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	//_____________________________________________________________________
	257	Int_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
8bc7e885	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	283	Bool_t AliTPCCalibPedestal::ProcessEvent(AliTPCRawStream *rawStream)
8bc7e885	284	{
8bc7e885	285	//
bc331d5b	286	// Event Processing loop - AliTPCRawStream
8bc7e885	287	//
8bc7e885	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	308	Bool_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	//_____________________________________________________________________
	322	Bool_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	333	Bool_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	356	TH2F* 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
8bc7e885	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	383	TH2F* 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	//_____________________________________________________________________
	393	AliTPCCalROC* 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	//_____________________________________________________________________
	410	AliTPCCalROC* 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	//_____________________________________________________________________
	420	AliTPCCalROC* 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	430	void 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	455	void 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	}
8bc7e885	490	}
8bc7e885	491	}
	492	//_____________________________________________________________________
	493	void 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	}