Corrected UInt_t <-> Int_t conversion
[u/mrichter/AliRoot.git] / TPC / AliTPCCalROC.cxx
index 609b306..38df7d3 100644 (file)
@@ -31,6 +31,7 @@
 #include "TH1F.h"
 #include "TH2F.h"
 #include "TArrayI.h"
+
 //
 //
 #include "AliTPCCalROC.h"
@@ -47,7 +48,7 @@ AliTPCCalROC::AliTPCCalROC()
              fSector(0),
              fNChannels(0),
              fNRows(0),
-             fIndexes(0),
+             fkIndexes(0),
              fData(0)
 {
   //
@@ -62,7 +63,7 @@ AliTPCCalROC::AliTPCCalROC(UInt_t  sector)
              fSector(0),
              fNChannels(0),
              fNRows(0),
-             fIndexes(0),
+             fkIndexes(0),
              fData(0)
 {
   //
@@ -71,7 +72,7 @@ AliTPCCalROC::AliTPCCalROC(UInt_t  sector)
   fSector = sector;
   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
-  fIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
+  fkIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
   fData = new Float_t[fNChannels];
   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = 0.;
 }
@@ -82,7 +83,7 @@ AliTPCCalROC::AliTPCCalROC(const AliTPCCalROC &c)
              fSector(0),
              fNChannels(0),
              fNRows(0),
-             fIndexes(0),
+             fkIndexes(0),
              fData(0)
 {
   //
@@ -91,7 +92,7 @@ AliTPCCalROC::AliTPCCalROC(const AliTPCCalROC &c)
   fSector = c.fSector;
   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
-  fIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
+  fkIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
   //
   fData   = new Float_t[fNChannels];
   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = c.fData[idata];
@@ -128,7 +129,7 @@ void AliTPCCalROC::Streamer(TBuffer &R__b)
    //
    if (R__b.IsReading()) {
       AliTPCCalROC::Class()->ReadBuffer(R__b, this);
-      fIndexes =  AliTPCROC::Instance()->GetRowIndexes(fSector);
+      fkIndexes =  AliTPCROC::Instance()->GetRowIndexes(fSector);
    } else {
       AliTPCCalROC::Class()->WriteBuffer(R__b,this);
    }
@@ -187,82 +188,96 @@ void AliTPCCalROC::Divide(const AliTPCCalROC*  roc) {
   }
 }
 
-Double_t AliTPCCalROC::GetMean(AliTPCCalROC* outlierROC) {
+Double_t AliTPCCalROC::GetMean(AliTPCCalROC *const outlierROC) const {
    //
    //  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;
+   Int_t nPoints = 0;
    for (UInt_t i=0;i<fNChannels;i++) {
       if (!(outlierROC->GetValue(i))) {
-         ddata[NPoints]= fData[NPoints];
-         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) {
+Double_t AliTPCCalROC::GetMedian(AliTPCCalROC *const outlierROC) const {
    //
    //  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;
+   Int_t nPoints = 0;
    for (UInt_t i=0;i<fNChannels;i++) {
       if (!(outlierROC->GetValue(i))) {
-         ddata[NPoints]= fData[NPoints];
-         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) {
+Double_t AliTPCCalROC::GetRMS(AliTPCCalROC *const outlierROC) const {
    //
    //  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;
+   Int_t nPoints = 0;
    for (UInt_t i=0;i<fNChannels;i++) {
       if (!(outlierROC->GetValue(i))) {
-         ddata[NPoints]= fData[NPoints];
-         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){
+Double_t AliTPCCalROC::GetLTM(Double_t *const sigma, Double_t fraction, AliTPCCalROC *const 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;
+  UInt_t nPoints = 0;
   for (UInt_t i=0;i<fNChannels;i++) {
      if (!outlierROC || !(outlierROC->GetValue(i))) {
-        ddata[NPoints]= fData[NPoints];
-        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){
@@ -299,9 +314,8 @@ TH1F * AliTPCCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type){
       max = mean+sigma;
     }
   }
-  char  name[1000];
-  sprintf(name,"%s ROC 1D%d",GetTitle(),fSector);
-  TH1F * his = new TH1F(name,name,100, min,max);
+  TString name=Form("%s ROC 1D%d",GetTitle(),fSector);
+  TH1F * his = new TH1F(name.Data(),name.Data(),100, min,max);
   for (UInt_t irow=0; irow<fNRows; irow++){
     UInt_t npads = (Int_t)GetNPads(irow);
     for (UInt_t ipad=0; ipad<npads; ipad++){
@@ -349,9 +363,9 @@ TH2F * AliTPCCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type){
   for (UInt_t irow=0; irow<fNRows; irow++){
     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
   }
-  char  name[1000];
-  sprintf(name,"%s ROC%d",GetTitle(),fSector);
-  TH2F * his = new TH2F(name,name,fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
+
+  TString name=Form("%s ROC%d",GetTitle(),fSector);
+  TH2F * his = new TH2F(name.Data(),name.Data(),fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
   for (UInt_t irow=0; irow<fNRows; irow++){
     UInt_t npads = (Int_t)GetNPads(irow);
     for (UInt_t ipad=0; ipad<npads; ipad++){
@@ -380,9 +394,8 @@ TH2F * AliTPCCalROC::MakeHistoOutliers(Float_t delta, Float_t fraction, Int_t ty
     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
   }
 
-  char  name[1000];
-  sprintf(name,"%s ROC Outliers%d",GetTitle(),fSector);
-  TH2F * his = new TH2F(name,name,fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
+  TString name=Form("%s ROC Outliers%d",GetTitle(),fSector);
+  TH2F * his = new TH2F(name.Data(),name.Data(),fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
   for (UInt_t irow=0; irow<fNRows; irow++){
     UInt_t npads = (Int_t)GetNPads(irow);
     for (UInt_t ipad=0; ipad<npads; ipad++){
@@ -566,20 +579,20 @@ AliTPCCalROC * AliTPCCalROC::LocalFit(Int_t rowRadius, Int_t padRadius, AliTPCCa
   // chi2Threshold: Threshold for chi2 when EvalRobust is called
   // robustFraction: Fraction of data that will be used in EvalRobust
   //
-  AliTPCCalROC * ROCfitted = new AliTPCCalROC(fSector);
+  AliTPCCalROC * xROCfitted = 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));
+      xROCfitted->SetValue(row, pad, GetNeighbourhoodValue(&fitterQ, row, pad, rowRadius, padRadius, ROCoutliers, robust, chi2Threshold, robustFraction));
   }
-  return ROCfitted;
+  return xROCfitted;
 }
 
 
-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) {
+Double_t AliTPCCalROC::GetNeighbourhoodValue(TLinearFitter* fitterQ, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius, AliTPCCalROC *const ROCoutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction) {
   //
   // AliTPCCalROC::GetNeighbourhoodValue - smoothing - PRIVATE
   // in this function the fit for LocalFit is done
@@ -681,7 +694,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 +715,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 +724,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 +750,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;
 }
 
@@ -796,36 +802,35 @@ AliTPCCalROC* AliTPCCalROC::CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sect
   Float_t dlx, dly;
   Float_t centerPad[3] = {0};
   Float_t localXY[3] = {0};
-  AliTPCCalROC * ROCfitted = new AliTPCCalROC(sector);
+  AliTPCCalROC * xROCfitted = new AliTPCCalROC(sector);
   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
-  tpcROCinstance->GetPositionLocal(sector, ROCfitted->GetNrows()/2, ROCfitted->GetNPads(ROCfitted->GetNrows()/2)/2, centerPad);  // calculate center of ROC 
+  tpcROCinstance->GetPositionLocal(sector, xROCfitted->GetNrows()/2, xROCfitted->GetNPads(xROCfitted->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++) {
+    for (UInt_t irow = 0; irow < xROCfitted->GetNrows(); irow++) {
+      for (UInt_t ipad = 0; ipad < xROCfitted->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);
+       xROCfitted->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++) {
+    for (UInt_t irow = 0; irow < xROCfitted->GetNrows(); irow++) {
+      for (UInt_t ipad = 0; ipad < xROCfitted->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);
+       xROCfitted->SetValue(irow, ipad, value);
       }
     }   
   }
-  return ROCfitted;
+  return xROCfitted;
 }
 
-