///////////////////////////////////////////////////////////////////////////////
// //
-// Calibration base class for a single ROC //
-// Contains one float value per pad //
+// Calibration base class for a single ROC //
+// Contains one float value per pad //
// mapping of the pads taken form AliTPCROC //
// //
///////////////////////////////////////////////////////////////////////////////
#include "TH1F.h"
#include "TH2F.h"
#include "TArrayI.h"
+
//
//
#include "AliTPCCalROC.h"
//_____________________________________________________________________________
AliTPCCalROC::AliTPCCalROC()
- :TObject(),
+ :TNamed(),
fSector(0),
fNChannels(0),
fNRows(0),
//_____________________________________________________________________________
AliTPCCalROC::AliTPCCalROC(UInt_t sector)
- :TObject(),
+ :TNamed(),
fSector(0),
fNChannels(0),
fNRows(0),
//_____________________________________________________________________________
AliTPCCalROC::AliTPCCalROC(const AliTPCCalROC &c)
- :TObject(c),
+ :TNamed(c),
fSector(0),
fNChannels(0),
fNRows(0),
void AliTPCCalROC::Streamer(TBuffer &R__b)
{
+ //
// Stream an object of class AliTPCCalROC.
+ //
if (R__b.IsReading()) {
AliTPCCalROC::Class()->ReadBuffer(R__b, this);
fIndexes = AliTPCROC::Instance()->GetRowIndexes(fSector);
}
}
-// //
-// // algebra
-// void Add(Float_t c1);
-// void Multiply(Float_t c1);
-// void Add(const AliTPCCalROC * roc, Double_t c1 = 1);
-// void Divide(const AliTPCCalROC * roc);
+
+// algebra fuctions:
void AliTPCCalROC::Add(Float_t c1){
//
- // add constant
+ // add c1 to each channel of the ROC
//
for (UInt_t idata = 0; idata< fNChannels; idata++) fData[idata]+=c1;
}
+
+
void AliTPCCalROC::Multiply(Float_t c1){
//
- // add constant
+ // multiply each channel of the ROC with c1
//
for (UInt_t idata = 0; idata< fNChannels; idata++) fData[idata]*=c1;
}
+
void AliTPCCalROC::Add(const AliTPCCalROC * roc, Double_t c1){
//
- // add values
+ // multiply AliTPCCalROC roc by c1 and add each channel to the coresponing channel in the ROC
+ // - pad by pad
//
for (UInt_t idata = 0; idata< fNChannels; idata++){
fData[idata]+=roc->fData[idata]*c1;
void AliTPCCalROC::Multiply(const AliTPCCalROC* roc) {
//
- // multiply values - per by pad
+ // multiply each channel of the ROC with the coresponding channel of 'roc'
+ // - pad by pad -
//
for (UInt_t idata = 0; idata< fNChannels; idata++){
fData[idata]*=roc->fData[idata];
void AliTPCCalROC::Divide(const AliTPCCalROC* roc) {
//
- // divide values
+ // divide each channel of the ROC by the coresponding channel of 'roc'
+ // - pad by pad -
//
Float_t kEpsilon=0.00000000000000001;
for (UInt_t idata = 0; idata< fNChannels; idata++){
}
Double_t AliTPCCalROC::GetMean(AliTPCCalROC* outlierROC) {
+ //
+ // returns the mean value of the ROC
+ // pads with value != 0 in outlierROC are not used for the calculation
+ // return 0 if no data is accepted by the outlier cuts
+ //
if (!outlierROC) return TMath::Mean(fNChannels, fData);
Double_t *ddata = new Double_t[fNChannels];
Int_t NPoints = 0;
for (UInt_t i=0;i<fNChannels;i++) {
if (!(outlierROC->GetValue(i))) {
- ddata[NPoints]= fData[NPoints];
+ ddata[NPoints]= fData[i];
NPoints++;
}
}
- Double_t mean = TMath::Mean(NPoints, ddata);
+ Double_t mean = 0;
+ if(NPoints>0)
+ mean = TMath::Mean(NPoints, ddata);
delete [] ddata;
return mean;
}
Double_t AliTPCCalROC::GetMedian(AliTPCCalROC* outlierROC) {
+ //
+ // returns the median value of the ROC
+ // pads with value != 0 in outlierROC are not used for the calculation
+ // return 0 if no data is accepted by the outlier cuts
+ //
if (!outlierROC) return TMath::Median(fNChannels, fData);
Double_t *ddata = new Double_t[fNChannels];
Int_t NPoints = 0;
for (UInt_t i=0;i<fNChannels;i++) {
if (!(outlierROC->GetValue(i))) {
- ddata[NPoints]= fData[NPoints];
+ ddata[NPoints]= fData[i];
NPoints++;
}
}
- Double_t mean = TMath::Median(NPoints, ddata);
+ Double_t median = 0;
+ if(NPoints>0)
+ median = TMath::Median(NPoints, ddata);
delete [] ddata;
- return mean;
+ return median;
}
Double_t AliTPCCalROC::GetRMS(AliTPCCalROC* outlierROC) {
+ //
+ // returns the RMS value of the ROC
+ // pads with value != 0 in outlierROC are not used for the calculation
+ // return 0 if no data is accepted by the outlier cuts
+ //
if (!outlierROC) return TMath::RMS(fNChannels, fData);
Double_t *ddata = new Double_t[fNChannels];
Int_t NPoints = 0;
for (UInt_t i=0;i<fNChannels;i++) {
if (!(outlierROC->GetValue(i))) {
- ddata[NPoints]= fData[NPoints];
+ ddata[NPoints]= fData[i];
NPoints++;
}
}
- Double_t mean = TMath::RMS(NPoints, ddata);
+ Double_t rms = 0;
+ if(NPoints>0)
+ rms = TMath::RMS(NPoints, ddata);
delete [] ddata;
- return mean;
+ return rms;
}
Double_t AliTPCCalROC::GetLTM(Double_t *sigma, Double_t fraction, AliTPCCalROC* outlierROC){
//
- // Calculate LTM mean and sigma
+ // returns the LTM and sigma
+ // pads with value != 0 in outlierROC are not used for the calculation
+ // return 0 if no data is accepted by the outlier cuts
//
Double_t *ddata = new Double_t[fNChannels];
- Double_t mean=0, lsigma=0;
UInt_t NPoints = 0;
for (UInt_t i=0;i<fNChannels;i++) {
if (!outlierROC || !(outlierROC->GetValue(i))) {
- ddata[NPoints]= fData[NPoints];
+ ddata[NPoints]= fData[i];
NPoints++;
}
}
- Int_t hh = TMath::Min(TMath::Nint(fraction *NPoints), Int_t(NPoints));
- AliMathBase::EvaluateUni(NPoints,ddata, mean, lsigma, hh);
- if (sigma) *sigma=lsigma;
+
+ Double_t ltm =0, lsigma=0;
+ if(NPoints>0) {
+ Int_t hh = TMath::Min(TMath::Nint(fraction *NPoints), Int_t(NPoints));
+ AliMathBase::EvaluateUni(NPoints,ddata, ltm, lsigma, hh);
+ if (sigma) *sigma=lsigma;
+ }
+
delete [] ddata;
- return mean;
+ return ltm;
}
TH1F * AliTPCCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type){
// type -1 = user defined range
// 0 = nsigma cut nsigma=min
// 1 = delta cut around median delta=min
+ //
if (type>=0){
if (type==0){
// nsigma range
// type -1 = user defined range
// 0 = nsigma cut nsigma=min
// 1 = delta cut around median delta=min
+ //
if (type>=0){
if (type==0){
// nsigma range
// fraction - fraction of values used to define sigma
// delta - in mode 0 - nsigma cut
// mode 1 - delta cut
+ //
Double_t sigma;
Float_t mean = GetLTM(&sigma,fraction);
if (type==0) delta*=sigma;
}
-AliTPCCalROC * AliTPCCalROC::LocalFit(Int_t rowRadius, Int_t padRadius, AliTPCCalROC* ROCoutliers, Bool_t robust) {
+AliTPCCalROC * AliTPCCalROC::LocalFit(Int_t rowRadius, Int_t padRadius, AliTPCCalROC* ROCoutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction) {
//
// MakeLocalFit - smoothing
+ // returns a AliTPCCalROC with smoothed data
+ // rowRadius and padRadius specifies a window around a given pad.
+ // The data of this window are fitted with a parabolic function.
+ // This function is evaluated at the pad's position.
+ // At the edges the window is shifted, so that the specified pad is not anymore in the center of the window.
// rowRadius - radius - rows to be used for smoothing
// padradius - radius - pads to be used for smoothing
// ROCoutlier - map of outliers - pads not to be used for local smoothing
// robust - robust method of fitting - (much slower)
-
+ // chi2Threshold: Threshold for chi2 when EvalRobust is called
+ // robustFraction: Fraction of data that will be used in EvalRobust
+ //
AliTPCCalROC * ROCfitted = new AliTPCCalROC(fSector);
- TLinearFitter fitterQ(6,"x0++x1++x2++x3++x4++x5");
+ TLinearFitter fitterQ(6,"hyp5");
+ // TLinearFitter fitterQ(6,"x0++x1++x2++x3++x4++x5");
fitterQ.StoreData(kTRUE);
for (UInt_t row=0; row < GetNrows(); row++) {
//std::cout << "Entering row " << row << " of " << GetNrows() << " @ sector "<< fSector << " for local fitting... "<< std::endl;
for (UInt_t pad=0; pad < GetNPads(row); pad++)
- ROCfitted->SetValue(row, pad, GetNeighbourhoodValue(&fitterQ, row, pad, rowRadius, padRadius, ROCoutliers, robust));
+ ROCfitted->SetValue(row, pad, GetNeighbourhoodValue(&fitterQ, row, pad, rowRadius, padRadius, ROCoutliers, robust, chi2Threshold, robustFraction));
}
return ROCfitted;
}
-Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius, AliTPCCalROC* ROCoutliers, Bool_t robust) {
+Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius, AliTPCCalROC* ROCoutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction) {
+ //
+ // AliTPCCalROC::GetNeighbourhoodValue - smoothing - PRIVATE
+ // in this function the fit for LocalFit is done
//
- // AliTPCCalROC::GetNeighbourhoodValue - smoothing (PRIVATE)
- // rowRadius - radius - rows to be used for smoothing
- // padradius - radius - pads to be used for smoothing
- // ROCoutlier - map of outliers - pads not to be used for local smoothing
- // robust - robust method of fitting - (much slower)
-
-
fitterQ->ClearPoints();
TVectorD fitParam(6);
tpcROCinstance->GetPositionLocal(fSector, r, p, localXY); // calculate position localXY by pad and row number
dlx = lPad[0] - localXY[0];
dly = lPad[1] - localXY[1];
- xx[0] = 1;
+ //xx[0] = 1;
xx[1] = dlx;
xx[2] = dly;
xx[3] = dlx*dlx;
xx[4] = dly*dly;
xx[5] = dlx*dly;
if (!ROCoutliers || ROCoutliers->GetValue(r,p) != 1) {
- fitterQ->AddPoint(xx, GetValue(r, p), 1);
+ fitterQ->AddPoint(&xx[1], GetValue(r, p), 1);
npoints++;
}
}
fitterQ->Eval();
fitterQ->GetParameters(fitParam);
Float_t chi2Q = 0;
- chi2Q = fitterQ->GetChisquare()/(npoints-6.);
+ if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
//if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
- if (robust && chi2Q > 5) {
+ if (robust && chi2Q > chi2Threshold) {
//std::cout << "robust fitter called... " << std::endl;
- fitterQ->EvalRobust(0.7);
+ fitterQ->EvalRobust(robustFraction);
fitterQ->GetParameters(fitParam);
}
Double_t value = fitParam[0];
//if (value < 0) std::cerr << "negative fit-value " << value << " in sector "<< this->fSector << " @ row: " << row << " and pad: " << pad << ", with fitter Chi2 = " << chi2Q << std::endl;
-
return value;
}
void AliTPCCalROC::GetNeighbourhood(TArrayI* &rowArray, TArrayI* &padArray, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius) {
//
- //
+ // AliTPCCalROC::GetNeighbourhood - PRIVATE
+ // in this function the window for LocalFit is determined
//
rowArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
padArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
- //Int_t* rowArrayTemp = rowArray->GetArray();
- //Int_t* padArrayTemp = padArray->GetArray();
Int_t rmin = row - rRadius;
UInt_t rmax = row + rRadius;
pmax = pmax - pmin;
pmin = 0;
}
- if (pmax >= GetNPads(r)) {
+ if (pmax >= (Int_t)GetNPads(r)) {
pmin = pmin - (pmax - GetNPads(r)+1);
pmax = GetNPads(r) - 1;
if (pmin < 0 ) pmin = 0; // if the window is bigger than the ROC
-void AliTPCCalROC::GlobalFit(const AliTPCCalROC* ROCoutliers, Bool_t robust, TVectorD &fitParam, TMatrixD &covMatrix, Float_t & chi2, Int_t fitType){
+void AliTPCCalROC::GlobalFit(const AliTPCCalROC* ROCoutliers, Bool_t robust, TVectorD &fitParam, TMatrixD &covMatrix, Float_t & chi2, Int_t fitType, Double_t chi2Threshold, Double_t robustFraction, Double_t err){
+ //
+ // Makes a GlobalFit for the given secotr and return fit-parameters, covariance and chi2
+ // The origin of the fit function is the center of the ROC!
+ // fitType == 0: fit plane function
+ // fitType == 1: fit parabolic function
+ // ROCoutliers - pads with value !=0 are not used in fitting procedure
+ // chi2Threshold: Threshold for chi2 when EvalRobust is called
+ // robustFraction: Fraction of data that will be used in EvalRobust
+ // err: error of the data points
//
- // Makes global fit
- // do GlobalFit for given Secotr and return fit-parameters, covariance, and whatever
- // fitType == 0: fit plane
- // fitType == 1: fit parabolic
- // ROCoutliers - pads signed=1 - not used in fitting procedure
-
TLinearFitter* fitterG = 0;
Double_t xx[6];
- if (fitType == 1)
+ if (fitType == 1) {
fitterG = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
- else
+ fitParam.ResizeTo(6);
+ covMatrix.ResizeTo(6,6);
+ } else {
fitterG = new TLinearFitter(3,"x0++x1++x2");
+ fitParam.ResizeTo(3);
+ covMatrix.ResizeTo(3,3);
+ }
fitterG->StoreData(kTRUE);
fitterG->ClearPoints();
Int_t npoints=0;
tpcROCinstance->GetPositionLocal(fSector, GetNrows()/2, GetNPads(GetNrows()/2)/2, centerPad); // calculate center of ROC
// loop over all channels and read data into fitterG
- if (fitType == 1) { // parabolic fit
- fitParam.ResizeTo(6);
- covMatrix.ResizeTo(6,6);
- for (UInt_t irow = 0; irow < GetNrows(); irow++) {
- for (UInt_t ipad = 0; ipad < GetNPads(irow); ipad++) {
- // fill fitterG
- tpcROCinstance->GetPositionLocal(fSector, irow, ipad, localXY); // calculate position localXY by pad and row number
- dlx = centerPad[0] - localXY[0];
- dly = centerPad[1] - localXY[1];
- xx[0] = 1;
- xx[1] = dlx;
- xx[2] = dly;
- xx[3] = dlx*dlx;
- xx[4] = dly*dly;
- xx[5] = dlx*dly;
- if (!ROCoutliers || ROCoutliers->GetValue(irow, ipad) != 1) {
- npoints++;
- fitterG->AddPoint(xx, GetValue(irow, ipad), 1);
- }
- }
- }
- }
- else { // linear fit
- fitParam.ResizeTo(3);
- covMatrix.ResizeTo(3,3);
- for (UInt_t irow = 0; irow < GetNrows(); irow++) {
- for (UInt_t ipad = 0; ipad < GetNPads(irow); ipad++) {
- // fill fitterG
- tpcROCinstance->GetPositionLocal(fSector, irow, ipad, localXY); // calculate position localXY by pad and row number
- dlx = centerPad[0] - localXY[0];
- dly = centerPad[1] - localXY[1];
- xx[0] = 1;
- xx[1] = dlx;
- xx[2] = dly;
- if (!ROCoutliers || ROCoutliers->GetValue(irow, ipad) != 1) {
- npoints++;
- fitterG->AddPoint(xx, GetValue(irow, ipad), 1);
- }
- }
+ for (UInt_t irow = 0; irow < GetNrows(); irow++) {
+ for (UInt_t ipad = 0; ipad < GetNPads(irow); ipad++) {
+ // fill fitterG
+ if (ROCoutliers && ROCoutliers->GetValue(irow, ipad) != 0) continue;
+ tpcROCinstance->GetPositionLocal(fSector, irow, ipad, localXY); // calculate position localXY by pad and row number
+ dlx = localXY[0] - centerPad[0];
+ dly = localXY[1] - centerPad[1];
+ xx[0] = 1;
+ xx[1] = dlx;
+ xx[2] = dly;
+ xx[3] = dlx*dlx;
+ xx[4] = dly*dly;
+ xx[5] = dlx*dly;
+ npoints++;
+ fitterG->AddPoint(xx, GetValue(irow, ipad), err);
}
}
- fitterG->Eval();
- fitterG->GetParameters(fitParam);
- fitterG->GetCovarianceMatrix(covMatrix);
- if (fitType == 1)
- chi2 = fitterG->GetChisquare()/(npoints-6.);
- else chi2 = fitterG->GetChisquare()/(npoints-3.);
- if (robust && chi2 > 5) {
- // std::cout << "robust fitter called... " << std::endl;
- fitterG->EvalRobust(0.7);
+ if(npoints>10) { // make sure there is something to fit
+ fitterG->Eval();
fitterG->GetParameters(fitParam);
+ fitterG->GetCovarianceMatrix(covMatrix);
+ if (fitType == 1)
+ chi2 = fitterG->GetChisquare()/(npoints-6.);
+ else chi2 = fitterG->GetChisquare()/(npoints-3.);
+ if (robust && chi2 > chi2Threshold) {
+ // std::cout << "robust fitter called... " << std::endl;
+ fitterG->EvalRobust(robustFraction);
+ fitterG->GetParameters(fitParam);
+ }
+ } else {
+ // set parameteres to 0
+ Int_t nParameters = 3;
+ if (fitType == 1)
+ nParameters = 6;
+
+ for(Int_t i = 0; i < nParameters; i++)
+ fitParam[i] = 0;
}
+
delete fitterG;
}
-//
AliTPCCalROC* AliTPCCalROC::CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sector){
- //
//
// Create ROC with global fit parameters
- // fitType == 0: fit plane
- // fitType == 1: fit parabolic
- // loop over all channels and write fit values into ROCfitted
+ // The origin of the fit function is the center of the ROC!
+ // loop over all channels, write fit values into new ROC and return it
//
Float_t dlx, dly;
Float_t centerPad[3] = {0};
for (UInt_t irow = 0; irow < ROCfitted->GetNrows(); irow++) {
for (UInt_t ipad = 0; ipad < ROCfitted->GetNPads(irow); ipad++) {
tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY); // calculate position localXY by pad and row number
- dlx = centerPad[0] - localXY[0];
- dly = centerPad[1] - localXY[1];
+ dlx = localXY[0] - centerPad[0];
+ dly = localXY[1] - centerPad[1];
value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
ROCfitted->SetValue(irow, ipad, value);
}
for (UInt_t irow = 0; irow < ROCfitted->GetNrows(); irow++) {
for (UInt_t ipad = 0; ipad < ROCfitted->GetNPads(irow); ipad++) {
tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY); // calculate position localXY by pad and row number
- dlx = centerPad[0] - localXY[0];
- dly = centerPad[1] - localXY[1];
+ dlx = localXY[0] - centerPad[0];
+ dly = localXY[1] - centerPad[1];
value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly;
ROCfitted->SetValue(irow, ipad, value);
}
return ROCfitted;
}
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff