]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCCalROC.cxx
Updated flags for low flux case (A. Dainese)
[u/mrichter/AliRoot.git] / TPC / AliTPCCalROC.cxx
index 609b306898fd299ff90e971af1c392cc52cb09b1..eb760541144b8096faa8424295e508dd54a247e4 100644 (file)
@@ -31,6 +31,7 @@
 #include "TH1F.h"
 #include "TH2F.h"
 #include "TArrayI.h"
+
 //
 //
 #include "AliTPCCalROC.h"
@@ -191,17 +192,20 @@ 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;
 }
@@ -210,59 +214,70 @@ 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){
   //
   //  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){
@@ -681,7 +696,7 @@ void AliTPCCalROC::GetNeighbourhood(TArrayI* &rowArray, TArrayI* &padArray, Int_
       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
@@ -702,7 +717,7 @@ void AliTPCCalROC::GetNeighbourhood(TArrayI* &rowArray, TArrayI* &padArray, Int_
 
 
 
-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){
+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!
@@ -711,14 +726,20 @@ void AliTPCCalROC::GlobalFit(const AliTPCCalROC* ROCoutliers, Bool_t robust, TVe
   // 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) 
+  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;
@@ -731,58 +752,45 @@ void AliTPCCalROC::GlobalFit(const AliTPCCalROC* ROCoutliers, Bool_t robust, TVe
   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 > chi2Threshold) {
-    //    std::cout << "robust fitter called... " << std::endl;
-    fitterG->EvalRobust(robustFraction);
+  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;
 }
 
@@ -807,8 +815,8 @@ AliTPCCalROC* AliTPCCalROC::CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sect
     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);
       }
@@ -818,8 +826,8 @@ AliTPCCalROC* AliTPCCalROC::CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sect
     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);
       }
@@ -828,4 +836,3 @@ AliTPCCalROC* AliTPCCalROC::CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sect
   return ROCfitted;
 }
 
-