X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCCalPad.cxx;h=fe68f65c30f59224c2b720144148b0b95feb77f4;hb=03190ad77fc63ad66fe07de295dbafd8b9e06f04;hp=20537a08fae2b134f7b01f1fd533f99566cd56ae;hpb=e6c51e6e8ecb4344e040fc8426201593f22a4e4f;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCCalPad.cxx b/TPC/AliTPCCalPad.cxx index 20537a08fae..fe68f65c30f 100644 --- a/TPC/AliTPCCalPad.cxx +++ b/TPC/AliTPCCalPad.cxx @@ -17,7 +17,8 @@ /////////////////////////////////////////////////////////////////////////////// // // -// TPC calibration class for parameters which saved per pad // +// TPC calibration class for parameters which are saved per pad // +// Each AliTPCCalPad consists of 72 AliTPCCalROC-objects // // // /////////////////////////////////////////////////////////////////////////////// @@ -29,6 +30,13 @@ #include #include #include "TTreeStream.h" +#include "TFile.h" +#include "TKey.h" +#include +#include +#include +#include +#include ClassImp(AliTPCCalPad) @@ -66,12 +74,10 @@ AliTPCCalPad::AliTPCCalPad(const AliTPCCalPad &c):TNamed(c) // for (Int_t isec = 0; isec < kNsec; isec++) { - fROC[isec] = 0; - if (c.fROC[isec]) - fROC[isec] = new AliTPCCalROC(*(c.fROC[isec])); + fROC[isec] = 0; + if (c.fROC[isec]) + fROC[isec] = new AliTPCCalROC(*(c.fROC[isec])); } - - } //_____________________________________________________________________________ @@ -131,11 +137,26 @@ void AliTPCCalPad::Copy(TObject &c) const TObject::Copy(c); } + +void AliTPCCalPad::SetCalROC(AliTPCCalROC* roc, Int_t sector){ + // + // Set AliTPCCalROC copies values from 'roc' + // if sector == -1 the sector specified in 'roc' is used + // else sector specified in 'roc' is ignored and specified sector is filled + // + if (sector == -1) sector = roc->GetSector(); + if (!fROC[sector]) fROC[sector] = new AliTPCCalROC(sector); + for (UInt_t ichannel = 0; ichannel < roc->GetNchannels(); ichannel++) + fROC[sector]->SetValue(ichannel, roc->GetValue(ichannel)); +} + + + //_____________________________________________________________________________ void AliTPCCalPad::Add(Float_t c1) { // - // add constant for all channels of all ROCs + // add constant c1 to all channels of all ROCs // for (Int_t isec = 0; isec < kNsec; isec++) { @@ -148,9 +169,9 @@ void AliTPCCalPad::Add(Float_t c1) //_____________________________________________________________________________ void AliTPCCalPad::Multiply(Float_t c1) { - // - // multiply constant for all channels of all ROCs - // + // + // multiply each channel of all ROCs with c1 + // for (Int_t isec = 0; isec < kNsec; isec++) { if (fROC[isec]){ fROC[isec]->Multiply(c1); @@ -161,11 +182,12 @@ void AliTPCCalPad::Multiply(Float_t c1) //_____________________________________________________________________________ void AliTPCCalPad::Add(const AliTPCCalPad * pad, Double_t c1) { - // - // add calpad channel by channel multiplied by c1 - all ROCs - // + // + // multiply AliTPCCalPad 'pad' by c1 and add each channel to the coresponing channel in all ROCs + // - pad by pad - + // for (Int_t isec = 0; isec < kNsec; isec++) { - if (fROC[isec]){ + if (fROC[isec] && pad->GetCalROC(isec)){ fROC[isec]->Add(pad->GetCalROC(isec),c1); } } @@ -174,10 +196,11 @@ void AliTPCCalPad::Add(const AliTPCCalPad * pad, Double_t c1) //_____________________________________________________________________________ void AliTPCCalPad::Multiply(const AliTPCCalPad * pad) { - // - // multiply calpad channel by channel - all ROCs - // - for (Int_t isec = 0; isec < kNsec; isec++) { + // + // multiply each channel of all ROCs with the coresponding channel of 'pad' + // - pad by pad - + // + for (Int_t isec = 0; isec < kNsec; isec++) { if (fROC[isec]){ fROC[isec]->Multiply(pad->GetCalROC(isec)); } @@ -187,9 +210,10 @@ void AliTPCCalPad::Multiply(const AliTPCCalPad * pad) //_____________________________________________________________________________ void AliTPCCalPad::Divide(const AliTPCCalPad * pad) { - // - // divide calpad channel by channel - all ROCs - // + // + // divide each channel of all ROCs by the coresponding channel of 'pad' + // - pad by pad - + // for (Int_t isec = 0; isec < kNsec; isec++) { if (fROC[isec]){ fROC[isec]->Divide(pad->GetCalROC(isec)); @@ -203,6 +227,7 @@ TGraph * AliTPCCalPad::MakeGraph(Int_t type, Float_t ratio){ // type=1 - mean // 2 - median // 3 - LTM + // Int_t npoints = 0; for (Int_t i=0;i<72;i++) if (fROC[i]) npoints++; TGraph * graph = new TGraph(npoints); @@ -236,7 +261,7 @@ TGraph * AliTPCCalPad::MakeGraph(Int_t type, Float_t ratio){ Double_t AliTPCCalPad::GetMeanRMS(Double_t &rms) { // - // Calculate mean an RMS of all rocs + // Calculates mean and RMS of all ROCs // Double_t sum = 0, sum2 = 0, n=0, val=0; for (Int_t isec = 0; isec < kNsec; isec++) { @@ -261,7 +286,7 @@ Double_t AliTPCCalPad::GetMeanRMS(Double_t &rms) //_____________________________________________________________________________ -Double_t AliTPCCalPad::GetMean() +Double_t AliTPCCalPad::GetMean(AliTPCCalPad* outlierPad) { // // return mean of the mean of all ROCs @@ -269,17 +294,19 @@ Double_t AliTPCCalPad::GetMean() Double_t arr[kNsec]; Int_t n=0; for (Int_t isec = 0; isec < kNsec; isec++) { - AliTPCCalROC *calRoc = fROC[isec]; - if ( calRoc ){ - arr[n] = calRoc->GetMean(); - n++; - } + AliTPCCalROC *calRoc = fROC[isec]; + if ( calRoc ){ + AliTPCCalROC* outlierROC = 0; + if (outlierPad) outlierROC = outlierPad->GetCalROC(isec); + arr[n] = calRoc->GetMean(outlierROC); + n++; + } } return TMath::Mean(n,arr); } //_____________________________________________________________________________ -Double_t AliTPCCalPad::GetRMS() +Double_t AliTPCCalPad::GetRMS(AliTPCCalPad* outlierPad) { // // return mean of the RMS of all ROCs @@ -287,17 +314,19 @@ Double_t AliTPCCalPad::GetRMS() Double_t arr[kNsec]; Int_t n=0; for (Int_t isec = 0; isec < kNsec; isec++) { - AliTPCCalROC *calRoc = fROC[isec]; - if ( calRoc ){ - arr[n] = calRoc->GetRMS(); - n++; - } + AliTPCCalROC *calRoc = fROC[isec]; + if ( calRoc ){ + AliTPCCalROC* outlierROC = 0; + if (outlierPad) outlierROC = outlierPad->GetCalROC(isec); + arr[n] = calRoc->GetRMS(outlierROC); + n++; + } } return TMath::Mean(n,arr); } //_____________________________________________________________________________ -Double_t AliTPCCalPad::GetMedian() +Double_t AliTPCCalPad::GetMedian(AliTPCCalPad* outlierPad) { // // return mean of the median of all ROCs @@ -305,17 +334,19 @@ Double_t AliTPCCalPad::GetMedian() Double_t arr[kNsec]; Int_t n=0; for (Int_t isec = 0; isec < kNsec; isec++) { - AliTPCCalROC *calRoc = fROC[isec]; - if ( calRoc ){ - arr[n] = calRoc->GetMedian(); - n++; - } + AliTPCCalROC *calRoc = fROC[isec]; + if ( calRoc ){ + AliTPCCalROC* outlierROC = 0; + if (outlierPad) outlierROC = outlierPad->GetCalROC(isec); + arr[n] = calRoc->GetMedian(outlierROC); + n++; + } } return TMath::Mean(n,arr); } //_____________________________________________________________________________ -Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction) +Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction, AliTPCCalPad* outlierPad) { // // return mean of the LTM and sigma of all ROCs @@ -329,7 +360,9 @@ Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction) AliTPCCalROC *calRoc = fROC[isec]; if ( calRoc ){ if ( sigma ) sTemp=arrs+n; - arrm[n] = calRoc->GetLTM(sTemp,fraction); + AliTPCCalROC* outlierROC = 0; + if (outlierPad) outlierROC = outlierPad->GetCalROC(isec); + arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC); n++; } } @@ -343,6 +376,7 @@ TH1F * AliTPCCalPad::MakeHisto1D(Float_t min, Float_t max,Int_t type){ // make 1D histo // type -1 = user defined range // 0 = nsigma cut nsigma=min + // if (type>=0){ if (type==0){ // nsigma range @@ -392,6 +426,7 @@ TH2F *AliTPCCalPad::MakeHisto2D(Int_t side){ // Make 2D graph // side - specify the side A = 0 C = 1 // type - used types of determination of boundaries in z + // Float_t kEpsilon = 0.000000000001; TH2F * his = new TH2F(GetName(), GetName(), 250,-250,250,250,-250,250); AliTPCROC * roc = AliTPCROC::Instance(); @@ -418,102 +453,368 @@ TH2F *AliTPCCalPad::MakeHisto2D(Int_t side){ } +AliTPCCalPad* AliTPCCalPad::LocalFit(const char* padName, Int_t rowRadius, Int_t padRadius, AliTPCCalPad* PadOutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction, Bool_t printCurrentSector) const { + // + // Loops over all AliTPCCalROCs and performs a localFit in each ROC + // AliTPCCalPad with fit-data is returned + // rowRadius and padRadius specifies a window around a given pad in one sector. + // 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 + // + // + AliTPCCalPad* pad = new AliTPCCalPad(padName, padName); + for (Int_t isec = 0; isec < 72; isec++){ + if (printCurrentSector) std::cout << "LocalFit in sector " << isec << "\r" << std::flush; + if (PadOutliers) + pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, PadOutliers->GetCalROC(isec), robust, chi2Threshold, robustFraction)); + else + pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, 0, robust, chi2Threshold, robustFraction)); + } + return pad; +} -void AliTPCCalPad::MakeTree(const char * fileName, TObjArray * array) { +AliTPCCalPad* AliTPCCalPad::GlobalFit(const char* padName, AliTPCCalPad* PadOutliers, Bool_t robust, Int_t fitType, Double_t chi2Threshold, Double_t robustFraction, Double_t err, TObjArray *fitParArr, TObjArray *fitCovArr){ + // + // Loops over all AliTPCCalROCs and performs a globalFit in each ROC + // AliTPCCalPad with fit-data is returned + // chi2Threshold: Threshold for chi2 when EvalRobust is called + // robustFraction: Fraction of data that will be used in EvalRobust + // chi2Threshold: Threshold for chi2 when EvalRobust is called + // robustFraction: Fraction of data that will be used in EvalRobust + // err: error of the data points + // if fitParArr and/or fitCovArr is given, write fitParameters and/or covariance Matrices into the array + // + AliTPCCalPad* pad = new AliTPCCalPad(padName, padName); + TVectorD fitParam(0); + TMatrixD covMatrix(0,0); + Float_t chi2 = 0; + for (Int_t isec = 0; isec < 72; isec++){ + if (PadOutliers) + GetCalROC(isec)->GlobalFit(PadOutliers->GetCalROC(isec), robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err); + else + GetCalROC(isec)->GlobalFit(0, robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction, err); + + AliTPCCalROC *roc=AliTPCCalROC::CreateGlobalFitCalROC(fitParam, isec); + pad->SetCalROC(roc); + delete roc; + if ( fitParArr ) fitParArr->AddAtAndExpand(new TVectorD(fitParam), isec); + if ( fitCovArr ) fitCovArr->AddAtAndExpand(new TMatrixD(covMatrix), isec); + } + return pad; +} +//_____________________________________________________________________________ +TObjArray* AliTPCCalPad::CreateFormulaArray(const char *fitFormula) +{ // - // Write tree with all available information + // create an array of TFormulas for the each parameter of the fit function // - TTreeSRedirector cstream(fileName); - AliTPCROC* tpcROCinstance = AliTPCROC::Instance(); - - TString* names = new TString[array->GetEntries()]; - for (Int_t ivalue = 0; ivalue < array->GetEntries(); ivalue++) - names[ivalue].Append(((AliTPCCalPad*)array->At(ivalue))->GetName()); - for (UInt_t isector = 0; isector < tpcROCinstance->GetNSectors(); isector++) { - // - // get statistic for given sector - // - TVectorF median(array->GetEntries()); - TVectorF mean(array->GetEntries()); - TVectorF rms(array->GetEntries()); - TVectorF ltm(array->GetEntries()); - - TVectorF *vectorArray = new TVectorF[array->GetEntries()]; - for (Int_t ivalue = 0; ivalue < array->GetEntries(); ivalue++) - vectorArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(isector)); - - for (Int_t ivalue = 0; ivalue < array->GetEntries(); ivalue++) { - AliTPCCalPad* calPad = (AliTPCCalPad*) array->At(ivalue); - AliTPCCalROC* calROC = calPad->GetCalROC(isector); - median[ivalue] = calROC->GetMedian(); - mean[ivalue] = calROC->GetMean(); - rms[ivalue] = calROC->GetRMS(); - ltm[ivalue] = calROC->GetLTM(); + // split fit string in single parameters + // find dimension of the fit: + TString fitString(fitFormula); + fitString.ReplaceAll("++","#"); + fitString.ReplaceAll(" ",""); + TObjArray *arrFitParams = fitString.Tokenize("#"); + Int_t ndim = arrFitParams->GetEntries(); + //create array of TFormulas to evaluate the parameters + TObjArray *arrFitFormulas = new TObjArray(ndim); + arrFitFormulas->SetOwner(kTRUE); + for (Int_t idim=0;idimAt(idim))->GetString(); + s.ReplaceAll("gx","[0]"); + s.ReplaceAll("gy","[1]"); + s.ReplaceAll("lx","[2]"); + s.ReplaceAll("ly","[3]"); + s.ReplaceAll("sector","[4]"); + arrFitFormulas->AddAt(new TFormula(Form("param%02d",idim),s.Data()),idim); + } + delete arrFitParams; + + return arrFitFormulas; +} +//_____________________________________________________________________________ +void AliTPCCalPad::EvalFormulaArray(const TObjArray &arrFitFormulas, TVectorD &results, + const Int_t sec, const Int_t row, const Int_t pad) +{ + // + // evaluate the fit formulas + // + Int_t ndim=arrFitFormulas.GetEntries(); + results.ResizeTo(ndim); + + AliTPCROC* tpcROCinstance = AliTPCROC::Instance(); // to calculate the pad's position + Float_t localXYZ[3]; + Float_t globalXYZ[3]; + tpcROCinstance->GetPositionLocal(sec, row, pad, localXYZ); + tpcROCinstance->GetPositionGlobal(sec, row, pad, globalXYZ); + //calculate parameter values + for (Int_t idim=0;idimSetParameters(globalXYZ[0],globalXYZ[1],localXYZ[0],localXYZ[1],sec); + results[idim]=f->Eval(0); + } +} +//_____________________________________________________________________________ +void AliTPCCalPad::GlobalSidesFit(const AliTPCCalPad* PadOutliers, const char* fitFormula, TVectorD &fitParamSideA, TVectorD &fitParamSideC,TMatrixD &covMatrixSideA, TMatrixD &covMatrixSideC, Float_t & chi2SideA, Float_t & chi2SideC, AliTPCCalPad *pointError, Bool_t robust, Double_t robustFraction){ + // + // Performs a fit on both sides. + // Valid information for the fitFormula are the variables + // - gx, gy, lx ,ly: meaning global x, global y, local x, local y value of the padName + // - sector: the sector number. + // eg. a formula might look 'gy' or '(sector<36) ++ gy' or 'gx ++ gy' or 'gx ++ gy ++ lx ++ lx^2' and so on + // + // PadOutliers - 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 + // + + TObjArray* arrFitFormulas=CreateFormulaArray(fitFormula); + Int_t ndim = arrFitFormulas->GetEntries(); + //resize output data arrays + fitParamSideA.ResizeTo(ndim+1); + fitParamSideC.ResizeTo(ndim+1); + covMatrixSideA.ResizeTo(ndim+1,ndim+1); + covMatrixSideC.ResizeTo(ndim+1,ndim+1); + // create linear fitter for A- and C- Side + TLinearFitter* fitterGA = new TLinearFitter(ndim+1,Form("hyp%d",ndim)); + TLinearFitter* fitterGC = new TLinearFitter(ndim+1,Form("hyp%d",ndim)); + fitterGA->StoreData(kTRUE); + fitterGC->StoreData(kTRUE); + //parameter values + TVectorD parValues(ndim); + + AliTPCCalROC *rocErr=0x0; + + for (UInt_t isec = 0; isecGetCalROC(isec); + AliTPCCalROC *rocData=GetCalROC(isec); + if (pointError) rocErr=pointError->GetCalROC(isec); + if (!rocData) continue; + for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) { + for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) { + //check for outliers + if (rocOut && rocOut->GetValue(irow,ipad)) continue; + //calculate parameter values + EvalFormulaArray(*arrFitFormulas,parValues,isec,irow,ipad); + //get value + Float_t value=rocData->GetValue(irow,ipad); + //point error + Int_t err=1; + if (rocErr) { + err=TMath::Nint(rocErr->GetValue(irow,ipad)); + if (err==0) err=1; + } + //add points to the fitters + if (isec/18%2==0){ + fitterGA->AddPoint(parValues.GetMatrixArray(),value,err); + }else{ + fitterGC->AddPoint(parValues.GetMatrixArray(),value,err); + } } - - // - // fill vectors of variable per pad - // - TVectorF *posArray = new TVectorF[6]; - for (Int_t ivalue = 0; ivalue < 6; ivalue++) - posArray[ivalue].ResizeTo(tpcROCinstance->GetNChannels(isector)); - - Float_t posG[3] = {0}; - Float_t posL[3] = {0}; - Int_t ichannel = 0; - for (UInt_t irow = 0; irow < tpcROCinstance->GetNRows(isector); irow++) { - for (UInt_t ipad = 0; ipad < tpcROCinstance->GetNPads(isector, irow); ipad++) { - tpcROCinstance->GetPositionLocal(isector, irow, ipad, posL); - tpcROCinstance->GetPositionGlobal(isector, irow, ipad, posG); - posArray[0][ichannel] = irow; - posArray[1][ichannel] = ipad; - posArray[2][ichannel] = posL[0]; - posArray[3][ichannel] = posL[1]; - posArray[4][ichannel] = posG[0]; - posArray[5][ichannel] = posG[1]; - - // loop over array containing AliTPCCalPads - for (Int_t ivalue = 0; ivalue < array->GetEntries(); ivalue++) { - AliTPCCalPad* calPad = (AliTPCCalPad*) array->At(ivalue); - (vectorArray[ivalue])[ichannel] = calPad->GetCalROC(isector)->GetValue(irow, ipad); + } + } + if (robust){ + fitterGA->EvalRobust(robustFraction); + fitterGC->EvalRobust(robustFraction); + } else { + fitterGA->Eval(); + fitterGC->Eval(); + } + chi2SideA=fitterGA->GetChisquare()/(fitterGA->GetNpoints()-(ndim+1)); + chi2SideC=fitterGC->GetChisquare()/(fitterGC->GetNpoints()-(ndim+1)); + fitterGA->GetParameters(fitParamSideA); + fitterGC->GetParameters(fitParamSideC); + fitterGA->GetCovarianceMatrix(covMatrixSideA); + fitterGC->GetCovarianceMatrix(covMatrixSideC); + + delete arrFitFormulas; + delete fitterGA; + delete fitterGC; + +} +// +AliTPCCalPad *AliTPCCalPad::CreateCalPadFit(const char* fitFormula, const TVectorD &fitParamSideA, const TVectorD &fitParamSideC) +{ + // + // + // + TObjArray *arrFitFormulas=CreateFormulaArray(fitFormula); + Int_t ndim = arrFitFormulas->GetEntries(); + //check if dimension of fit formula and fit parameters agree + if (ndim!=fitParamSideA.GetNrows()||ndim!=fitParamSideC.GetNrows()){ + printf("AliTPCCalPad::CreateCalPadFit: Dimensions of fit formula and fit Parameters does not match!"); + return 0; + } + //create cal pad + AliTPCCalPad *pad=new AliTPCCalPad("fitResultPad",Form("Fit result: %s",fitFormula)); + //fill cal pad with fit results if requested + for (UInt_t isec = 0; isecGetCalROC(isec); + for (UInt_t irow = 0; irow < roc->GetNrows(); irow++) { + for (UInt_t ipad = 0; ipad < roc->GetNPads(irow); ipad++) { + const TVectorD *fitPar=0; + TVectorD fitResArray; + if (isec/18%2==0){ + fitPar=&fitParamSideA; + }else{ + fitPar=&fitParamSideC; + } + EvalFormulaArray(*arrFitFormulas,fitResArray, isec, irow, ipad); + for (Int_t idim=0;idimSetValue(irow,ipad,fitResArray.Sum()); + } + } + } + delete arrFitFormulas; + return pad; +} +/* +void AliTPCCalPad::GlobalSidesFit(const AliTPCCalPad* PadOutliers, TVectorD &fitParamSideA, TVectorD &fitParamSideC,TMatrixD &covMatrixSideA, TMatrixD &covMatrixSideC, Float_t & chi2SideA, Float_t & chi2SideC, Int_t fitType, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction){ + // + // Makes a GlobalFit over each side and return fit-parameters, covariance and chi2 for each side + // fitType == 0: fit plane function + // fitType == 1: fit parabolic function + // PadOutliers - 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 + // + TLinearFitter* fitterGA = 0; + TLinearFitter* fitterGC = 0; + + if (fitType == 1) { + fitterGA = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5"); + fitterGC = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5"); + } + else { + fitterGA = new TLinearFitter(3,"x0++x1++x2"); + fitterGC = new TLinearFitter(3,"x0++x1++x2"); + } + fitterGA->StoreData(kTRUE); + fitterGC->StoreData(kTRUE); + fitterGA->ClearPoints(); + fitterGC->ClearPoints(); + Double_t xx[6]; + Int_t npointsA=0; + Int_t npointsC=0; + + Float_t localXY[3] = {0}; // pad's position, needed to get the pad's position + Float_t lx, ly; // pads position + + AliTPCROC* tpcROCinstance = AliTPCROC::Instance(); // to calculate the pad's position + + // loop over all sectors and pads and read data into fitterGA and fitterGC + if (fitType == 1) { + // parabolic fit + fitParamSideA.ResizeTo(6); + fitParamSideC.ResizeTo(6); + covMatrixSideA.ResizeTo(6,6); + covMatrixSideC.ResizeTo(6,6); + for (UInt_t isec = 0; isec<72; isec++){ + for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) { + for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) { + // fill fitterG + tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY); // calculate position localXY by sector, pad and row number + lx = localXY[0]; + ly = localXY[1]; + xx[0] = 1; + xx[1] = lx; + xx[2] = ly; + xx[3] = lx*lx; + xx[4] = ly*ly; + xx[5] = lx*ly; + if (!PadOutliers || PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) { + // if given pad is no outlier, add it to TLinearFitter, decide to which of both +// sector 0 - 17: IROC, A +// sector 18 - 35: IROC, C +// sector 36 - 53: OROC, A +// sector 54 - 71: CROC, C + if (isec <= 17 || (isec >= 36 && isec <= 53)) { // Side A + npointsA++; + fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1); + } + else { // side C + npointsC++; + fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1); + } } - ichannel++; } } - - cstream << "calPads" << - "sector=" << isector; - - for (Int_t ivalue = 0; ivalue < array->GetEntries(); ivalue++) { - cstream << "calPads" << - (Char_t*)((names[ivalue] + "Median=").Data()) << median[ivalue] << - (Char_t*)((names[ivalue] + "Mean=").Data()) << mean[ivalue] << - (Char_t*)((names[ivalue] + "RMS=").Data()) << rms[ivalue] << - (Char_t*)((names[ivalue] + "LTM=").Data()) << ltm[ivalue]; - } - - for (Int_t ivalue = 0; ivalue < array->GetEntries(); ivalue++) { - cstream << "calPads" << - (Char_t*)((names[ivalue] + ".=").Data()) << &vectorArray[ivalue]; + } + } + else { + // linear fit + fitParamSideA.ResizeTo(3); + fitParamSideC.ResizeTo(3); + covMatrixSideA.ResizeTo(3,3); + covMatrixSideC.ResizeTo(3,3); + + for (UInt_t isec = 0; isec<72; isec++){ + for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) { + for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) { + // fill fitterG + tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY); // calculate position localXY by sector, pad and row number + lx = localXY[0]; + ly = localXY[1]; + xx[0] = 1; + xx[1] = lx; + xx[2] = ly; + if (!PadOutliers || PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) { + // if given pad is no outlier, add it to TLinearFitter, decide to which of both +// sector 0 - 17: IROC, A +// sector 18 - 35: IROC, C +// sector 36 - 53: OROC, A +// sector 54 - 71: CROC, C + if (isec <= 17 || (isec >= 36 && isec <= 53)) { + // Side A + npointsA++; + fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1); + } + else { + // side C + npointsC++; + fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1); + } + } + } } - - cstream << "calPads" << - "row.=" << &posArray[0] << - "pad.=" << &posArray[1] << - "lx.=" << &posArray[2] << - "ly.=" << &posArray[3] << - "gx.=" << &posArray[4] << - "gy.=" << &posArray[5]; - - cstream << "calPads" << - "\n"; - - delete[] posArray; - delete[] vectorArray; - } - delete[] names; + } + } + + fitterGA->Eval(); + fitterGC->Eval(); + fitterGA->GetParameters(fitParamSideA); + fitterGC->GetParameters(fitParamSideC); + fitterGA->GetCovarianceMatrix(covMatrixSideA); + fitterGC->GetCovarianceMatrix(covMatrixSideC); + if (fitType == 1){ + chi2SideA = fitterGA->GetChisquare()/(npointsA-6.); + chi2SideC = fitterGC->GetChisquare()/(npointsC-6.); + } + else { + chi2SideA = fitterGA->GetChisquare()/(npointsA-3.); + chi2SideC = fitterGC->GetChisquare()/(npointsC-3.); + } + if (robust && chi2SideA > chi2Threshold) { + // std::cout << "robust fitter called... " << std::endl; + fitterGA->EvalRobust(robustFraction); + fitterGA->GetParameters(fitParamSideA); + } + if (robust && chi2SideC > chi2Threshold) { + // std::cout << "robust fitter called... " << std::endl; + fitterGC->EvalRobust(robustFraction); + fitterGC->GetParameters(fitParamSideC); + } + delete fitterGA; + delete fitterGC; } - +*/