X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=TPC%2FAliTPCCalROC.cxx;h=609b306898fd299ff90e971af1c392cc52cb09b1;hb=f11b30713680a73c0edf325133c3943e12c65b99;hp=831c3e11489d66db078392c86d5eaecd9aa03628;hpb=92e56aebddc7d84b1e6ff17234de408af0f9863e;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCCalROC.cxx b/TPC/AliTPCCalROC.cxx index 831c3e11489..609b306898f 100644 --- a/TPC/AliTPCCalROC.cxx +++ b/TPC/AliTPCCalROC.cxx @@ -16,8 +16,8 @@ /////////////////////////////////////////////////////////////////////////////// // // -// 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 // // // /////////////////////////////////////////////////////////////////////////////// @@ -36,12 +36,14 @@ #include "AliTPCCalROC.h" #include "AliMathBase.h" +#include "TRandom3.h" // only needed by the AliTPCCalROCTest() method + ClassImp(AliTPCCalROC) //_____________________________________________________________________________ AliTPCCalROC::AliTPCCalROC() - :TObject(), + :TNamed(), fSector(0), fNChannels(0), fNRows(0), @@ -56,7 +58,7 @@ AliTPCCalROC::AliTPCCalROC() //_____________________________________________________________________________ AliTPCCalROC::AliTPCCalROC(UInt_t sector) - :TObject(), + :TNamed(), fSector(0), fNChannels(0), fNRows(0), @@ -76,7 +78,7 @@ AliTPCCalROC::AliTPCCalROC(UInt_t sector) //_____________________________________________________________________________ AliTPCCalROC::AliTPCCalROC(const AliTPCCalROC &c) - :TObject(c), + :TNamed(c), fSector(0), fNChannels(0), fNRows(0), @@ -121,7 +123,9 @@ AliTPCCalROC::~AliTPCCalROC() 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); @@ -130,29 +134,29 @@ void AliTPCCalROC::Streamer(TBuffer &R__b) } } -// // -// // 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; @@ -162,7 +166,8 @@ void AliTPCCalROC::Add(const AliTPCCalROC * roc, Double_t 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]; @@ -172,7 +177,8 @@ void AliTPCCalROC::Multiply(const AliTPCCalROC* roc) { 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++){ @@ -181,18 +187,79 @@ void AliTPCCalROC::Divide(const AliTPCCalROC* roc) { } } +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 + // + if (!outlierROC) return TMath::Mean(fNChannels, fData); + Double_t *ddata = new Double_t[fNChannels]; + Int_t NPoints = 0; + for (UInt_t i=0;iGetValue(i))) { + ddata[NPoints]= fData[NPoints]; + NPoints++; + } + } + Double_t 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 + // + if (!outlierROC) return TMath::Median(fNChannels, fData); + Double_t *ddata = new Double_t[fNChannels]; + Int_t NPoints = 0; + for (UInt_t i=0;iGetValue(i))) { + ddata[NPoints]= fData[NPoints]; + NPoints++; + } + } + Double_t mean = TMath::Median(NPoints, ddata); + delete [] ddata; + return mean; +} +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 + // + if (!outlierROC) return TMath::RMS(fNChannels, fData); + Double_t *ddata = new Double_t[fNChannels]; + Int_t NPoints = 0; + for (UInt_t i=0;iGetValue(i))) { + ddata[NPoints]= fData[NPoints]; + NPoints++; + } + } + Double_t mean = TMath::RMS(NPoints, ddata); + delete [] ddata; + return mean; +} -Double_t AliTPCCalROC::GetLTM(Double_t *sigma, Double_t fraction){ +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 // Double_t *ddata = new Double_t[fNChannels]; Double_t mean=0, lsigma=0; - Int_t hh = TMath::Min(TMath::Nint(fraction *fNChannels), Int_t(fNChannels)); - for (UInt_t i=0;iGetValue(i))) { + ddata[NPoints]= fData[NPoints]; + NPoints++; + } + } + Int_t hh = TMath::Min(TMath::Nint(fraction *NPoints), Int_t(NPoints)); + AliMathBase::EvaluateUni(NPoints,ddata, mean, lsigma, hh); if (sigma) *sigma=lsigma; delete [] ddata; return mean; @@ -204,6 +271,7 @@ 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 @@ -251,6 +319,7 @@ TH2F * AliTPCCalROC::MakeHisto2D(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 @@ -302,6 +371,7 @@ TH2F * AliTPCCalROC::MakeHistoOutliers(Float_t delta, Float_t fraction, Int_t ty // 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; @@ -345,10 +415,14 @@ void AliTPCCalROC::Draw(Option_t* opt){ -void AliTPCCalROC::Test(){ +void AliTPCCalROC::Test() { // - // example function to show functionality and tes AliTPCCalROC + // example function to show functionality and test AliTPCCalROC // + + Float_t kEpsilon=0.00001; + + // create CalROC with defined values AliTPCCalROC roc0(0); for (UInt_t irow = 0; irow GetNPads(irow)) printf("NPads - Read/Write error\trow=%d\n",irow); for (UInt_t ipad = 0; ipad GetValue(irow,ipad)!=value){ - printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad); + printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad); + } + } + } + + // + // Algebra Tests + // + + // Add constant + AliTPCCalROC roc3(roc0); + roc3.Add(1.5); + for (UInt_t irow = 0; irow kEpsilon){ + printf("Add constant - error\trow=%d\tpad=%d\n",irow,ipad); + } + } + } + + // Add another CalROC + AliTPCCalROC roc4(roc0); + roc4.Add(&roc0, -1.5); + for (UInt_t irow = 0; irow kEpsilon){ + printf("Add CalROC - error\trow=%d\tpad=%d\n",irow,ipad); + } + } + } + + // Multiply with constant + AliTPCCalROC roc5(roc0); + roc5.Multiply(-1.4); + for (UInt_t irow = 0; irow kEpsilon){ + printf("Multiply with constant - error\trow=%d\tpad=%d\n",irow,ipad); + } + } + } + + // Multiply another CalROC + AliTPCCalROC roc6(roc0); + roc6.Multiply(&roc0); + for (UInt_t irow = 0; irow kEpsilon*100){ + printf("Multiply with CalROC - error\trow=%d\tpad=%d\n",irow,ipad); } } - } + } + + + // Divide by CalROC + AliTPCCalROC roc7(roc0); + roc7.Divide(&roc0); + for (UInt_t irow = 0; irow kEpsilon){ + printf("Multiply with CalROC - error\trow=%d\tpad=%d\n",irow,ipad); + } + } + } + + // + // Statistics Test + // + + // create CalROC with defined values + TRandom3 rnd(0); + AliTPCCalROC sroc0(0); + for (UInt_t ichannel = 0; ichannel < sroc0.GetNchannels(); ichannel++){ + Float_t value = rnd.Gaus(10., 2.); + sroc0.SetValue(ichannel,value); + } + + printf("Mean (should be close to 10): %f\n", sroc0.GetMean()); + printf("RMS (should be close to 2): %f\n", sroc0.GetRMS()); + printf("Median (should be close to 10): %f\n", sroc0.GetMedian()); + printf("LTM (should be close to 10): %f\n", sroc0.GetLTM()); + + //AliTPCCalROC* sroc1 = sroc0.LocalFit(4, 8); + + //delete sroc1; + +// std::cout << TMath::Abs(roc5.GetValue(irow,ipad)-value) << std::endl; } -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); @@ -427,27 +598,34 @@ Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row, TArrayI *neighbourhoodRows = 0; TArrayI *neighbourhoodPads = 0; + + //std::cerr << "Trying to get neighbourhood for row " << row << ", pad " << pad << std::endl; GetNeighbourhood(neighbourhoodRows, neighbourhoodPads, row, pad, rRadius, pRadius); + //std::cerr << "Got neighbourhood for row " << row << ", pad " << pad << std::endl; Int_t r, p; for (Int_t i=0; i < (2*rRadius+1)*(2*pRadius+1); i++) { r = neighbourhoodRows->At(i); p = neighbourhoodPads->At(i); - if (r == -1 || p == -1) continue; + if (r == -1 || p == -1) continue; // window is bigger than ROC 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); + if (!ROCoutliers || ROCoutliers->GetValue(r,p) != 1) { + fitterQ->AddPoint(&xx[1], GetValue(r, p), 1); npoints++; } } + + delete neighbourhoodRows; + delete neighbourhoodPads; + if (npoints < 0.5 * ((2*rRadius+1)*(2*pRadius+1)) ) { // std::cerr << "Too few data points for fitting @ row " << row << ", pad " << pad << " in sector " << fSector << std::endl; return 0.; // for diagnostic @@ -455,20 +633,16 @@ Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row, 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]; - delete neighbourhoodRows; - delete neighbourhoodPads; - //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; } @@ -477,12 +651,11 @@ Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row, 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; @@ -514,29 +687,31 @@ void AliTPCCalROC::GetNeighbourhood(TArrayI* &rowArray, TArrayI* &padArray, Int_ if (pmin < 0 ) pmin = 0; // if the window is bigger than the ROC } for (Int_t p = pmin; p <= pmax; p++) { - rowArrayTemp[i] = r; - padArrayTemp[i] = p; + (*rowArray)[i] = r; + (*padArray)[i] = p; i++; } } for (Int_t j = i; j < rowArray->GetSize(); j++){ // unused padArray-entries, in the case that the window is bigger than the ROC //std::cout << "trying to write -1" << std::endl; - rowArrayTemp[j] = -1; - padArrayTemp[j] = -1; + (*rowArray)[j] = -1; + (*padArray)[j] = -1; //std::cout << "writing -1" << std::endl; - } + } } -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){ + // + // 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 // - // 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]; @@ -571,7 +746,7 @@ void AliTPCCalROC::GlobalFit(const AliTPCCalROC* ROCoutliers, Bool_t robust, TVe xx[3] = dlx*dlx; xx[4] = dly*dly; xx[5] = dlx*dly; - if (ROCoutliers && ROCoutliers->GetValue(irow, ipad) != 1) { + if (!ROCoutliers || ROCoutliers->GetValue(irow, ipad) != 1) { npoints++; fitterG->AddPoint(xx, GetValue(irow, ipad), 1); } @@ -590,7 +765,7 @@ void AliTPCCalROC::GlobalFit(const AliTPCCalROC* ROCoutliers, Bool_t robust, TVe xx[0] = 1; xx[1] = dlx; xx[2] = dly; - if (ROCoutliers && ROCoutliers->GetValue(irow, ipad) != 1) { + if (!ROCoutliers || ROCoutliers->GetValue(irow, ipad) != 1) { npoints++; fitterG->AddPoint(xx, GetValue(irow, ipad), 1); } @@ -603,23 +778,20 @@ void AliTPCCalROC::GlobalFit(const AliTPCCalROC* ROCoutliers, Bool_t robust, TVe if (fitType == 1) chi2 = fitterG->GetChisquare()/(npoints-6.); else chi2 = fitterG->GetChisquare()/(npoints-3.); - if (robust && chi2 > 5) { + if (robust && chi2 > chi2Threshold) { // std::cout << "robust fitter called... " << std::endl; - fitterG->EvalRobust(0.7); + fitterG->EvalRobust(robustFraction); fitterG->GetParameters(fitParam); } 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};