+
+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,"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, 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 chi2Threshold, Double_t robustFraction) {
+ //
+ // AliTPCCalROC::GetNeighbourhoodValue - smoothing - PRIVATE
+ // in this function the fit for LocalFit is done
+ //
+
+ fitterQ->ClearPoints();
+ TVectorD fitParam(6);
+ Int_t npoints=0;
+ Double_t xx[6];
+ Float_t dlx, dly;
+ Float_t lPad[3] = {0};
+ Float_t localXY[3] = {0};
+
+ AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
+ tpcROCinstance->GetPositionLocal(fSector, row, pad, lPad); // calculate position lPad by pad and row number
+
+ 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; // 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[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[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
+ }
+ fitterQ->Eval();
+ fitterQ->GetParameters(fitParam);
+ Float_t chi2Q = 0;
+ if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
+ //if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
+ if (robust && chi2Q > chi2Threshold) {
+ //std::cout << "robust fitter called... " << std::endl;
+ 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 rmin = row - rRadius;
+ UInt_t rmax = row + rRadius;
+
+ // if window goes out of ROC
+ if (rmin < 0) {
+ rmax = rmax - rmin;
+ rmin = 0;
+ }
+ if (rmax >= GetNrows()) {
+ rmin = rmin - (rmax - GetNrows()+1);
+ rmax = GetNrows() - 1;
+ if (rmin < 0 ) rmin = 0; // if the window is bigger than the ROC
+ }
+
+ Int_t pmin, pmax;
+ Int_t i = 0;
+
+ for (UInt_t r = rmin; r <= rmax; r++) {
+ pmin = pad - pRadius;
+ pmax = pad + pRadius;
+ if (pmin < 0) {
+ pmax = pmax - pmin;
+ pmin = 0;
+ }
+ 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
+ }
+ for (Int_t p = pmin; p <= pmax; 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;
+ (*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, 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
+ //
+ TLinearFitter* fitterG = 0;
+ Double_t xx[6];
+
+ if (fitType == 1) {
+ fitterG = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
+ 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;
+
+ Float_t dlx, dly;
+ Float_t centerPad[3] = {0};
+ Float_t localXY[3] = {0};
+
+ AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
+ tpcROCinstance->GetPositionLocal(fSector, GetNrows()/2, GetNPads(GetNrows()/2)/2, centerPad); // calculate center of ROC
+
+ // loop over all channels and read data into fitterG
+ 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 > chi2Threshold) {
+ // std::cout << "robust fitter called... " << std::endl;
+ fitterG->EvalRobust(robustFraction);
+ fitterG->GetParameters(fitParam);
+ }
+ delete fitterG;
+}
+
+
+AliTPCCalROC* AliTPCCalROC::CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sector){
+ //
+ // Create ROC with global fit parameters
+ // 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};
+ Float_t localXY[3] = {0};
+ AliTPCCalROC * ROCfitted = new AliTPCCalROC(sector);
+ AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
+ tpcROCinstance->GetPositionLocal(sector, ROCfitted->GetNrows()/2, ROCfitted->GetNPads(ROCfitted->GetNrows()/2)/2, centerPad); // calculate center of ROC
+ Int_t fitType = 1;
+ if (fitParam.GetNoElements() == 6) fitType = 1;
+ else fitType = 0;
+ Double_t value = 0;
+ if (fitType == 1) { // parabolic fit
+ 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 = 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);
+ }
+ }
+ }
+ else { // linear fit
+ 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 = 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;
+}
+
+