[u/mrichter/AliRoot.git] / TPC / AliTPCCalPad.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//  TPC calibration class for parameters which are saved per pad             //
//  Each AliTPCCalPad consists of 72 AliTPCCalROC-objects                    //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#include "AliTPCCalPad.h"
#include "AliTPCCalROC.h"
#include <TObjArray.h>
#include <TAxis.h>
#include <TGraph.h>
#include <TGraph2D.h>
#include <TH2F.h>
#include "TTreeStream.h"
#include "TFile.h"
#include "TKey.h"
#include <iostream>

ClassImp(AliTPCCalPad)

//_____________________________________________________________________________
AliTPCCalPad::AliTPCCalPad():TNamed()
{
  //
  // AliTPCCalPad default constructor
  //

  for (Int_t isec = 0; isec < kNsec; isec++) {
    fROC[isec] = 0;
  }

}

//_____________________________________________________________________________
AliTPCCalPad::AliTPCCalPad(const Text_t *name, const Text_t *title)
                :TNamed(name,title)
{
  //
  // AliTPCCalPad constructor
  //
  for (Int_t isec = 0; isec < kNsec; isec++) {
    fROC[isec] = new AliTPCCalROC(isec);
  }
}


//_____________________________________________________________________________
AliTPCCalPad::AliTPCCalPad(const AliTPCCalPad &c):TNamed(c)
{
  //
  // AliTPCCalPad copy constructor
  //

  for (Int_t isec = 0; isec < kNsec; isec++) {
         fROC[isec] = 0;
     if (c.fROC[isec])
       fROC[isec] = new AliTPCCalROC(*(c.fROC[isec]));
  }
}

//_____________________________________________________________________________
AliTPCCalPad::AliTPCCalPad(TObjArray * array):TNamed()
{
  //
  // AliTPCCalPad default constructor
  //

  for (Int_t isec = 0; isec < kNsec; isec++) {
    fROC[isec] = (AliTPCCalROC *)array->At(isec);
  }

}


///_____________________________________________________________________________
AliTPCCalPad::~AliTPCCalPad()
{
  //
  // AliTPCCalPad destructor
  //

  for (Int_t isec = 0; isec < kNsec; isec++) {
    if (fROC[isec]) {
      delete fROC[isec];
      fROC[isec] = 0;
    }
  }

}

//_____________________________________________________________________________
AliTPCCalPad &AliTPCCalPad::operator=(const AliTPCCalPad &c)
{
  //
  // Assignment operator
  //

  if (this != &c) ((AliTPCCalPad &) c).Copy(*this);
  return *this;

}

//_____________________________________________________________________________
void AliTPCCalPad::Copy(TObject &c) const
{
  //
  // Copy function
  //

  for (Int_t isec = 0; isec < kNsec; isec++) {
    if (fROC[isec]) {
      fROC[isec]->Copy(*((AliTPCCalPad &) c).fROC[isec]);
    }
  }
  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 c1 to all channels of all ROCs
    //

    for (Int_t isec = 0; isec < kNsec; isec++) {
	if (fROC[isec]){
	    fROC[isec]->Add(c1);
	}
    }
}

//_____________________________________________________________________________
void AliTPCCalPad::Multiply(Float_t c1)
{
  //
  // multiply each channel of all ROCs with c1
  //    
    for (Int_t isec = 0; isec < kNsec; isec++) {
	if (fROC[isec]){
	    fROC[isec]->Multiply(c1);
	}
    }
}

//_____________________________________________________________________________
void AliTPCCalPad::Add(const AliTPCCalPad * pad, Double_t c1)
{
  //
  // 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]){
	    fROC[isec]->Add(pad->GetCalROC(isec),c1);
	}
    }
}

//_____________________________________________________________________________
void AliTPCCalPad::Multiply(const AliTPCCalPad * pad)
{
  //
  // 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));
	}
    }
}

//_____________________________________________________________________________
void AliTPCCalPad::Divide(const AliTPCCalPad * pad)
{
  //
  // 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));
	}
    }
}

//_____________________________________________________________________________
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);
  npoints=0;   
  for (Int_t isec=0;isec<72;isec++){
    if (!fROC[isec]) continue;
    if (type==0)  graph->SetPoint(npoints,isec,fROC[isec]->GetMean());      
    if (type==1)  graph->SetPoint(npoints,isec,fROC[isec]->GetMedian());
    if (type==2)  graph->SetPoint(npoints,isec,fROC[isec]->GetLTM(0,ratio));    
    npoints++;
  }

  graph->GetXaxis()->SetTitle("Sector"); 
  if (type==0) {
    graph->GetYaxis()->SetTitle("Mean");   
    graph->SetMarkerStyle(22);    
  }
  if (type==1) {
    graph->GetYaxis()->SetTitle("Median");   
    graph->SetMarkerStyle(22);    
  }
  if (type==2) {
      graph->GetYaxis()->SetTitle(Form("Mean%f",ratio));      
      graph->SetMarkerStyle(24);
  }

  return graph;
}

//_____________________________________________________________________________
Double_t AliTPCCalPad::GetMeanRMS(Double_t &rms)
{
    //
    // 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++) {
        AliTPCCalROC *calRoc = fROC[isec];
	if ( calRoc ){
	    for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++){
		for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++){
		    val = calRoc->GetValue(irow,ipad);
		    sum+=val;
		    sum2+=val*val;
                    n++;
		}
	    }

	}
    }
    Double_t n1 = 1./n;
    Double_t mean = sum*n1;
    rms  = TMath::Sqrt(TMath::Abs(sum2*n1-mean*mean));
    return mean;
}


//_____________________________________________________________________________
Double_t AliTPCCalPad::GetMean(AliTPCCalPad* outlierPad)
{
    //
    // return mean of the mean of all ROCs
    //
    Double_t arr[kNsec];
    Int_t n=0;
    for (Int_t isec = 0; isec < kNsec; isec++) {
       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(AliTPCCalPad* outlierPad)
{
    //
    // return mean of the RMS of all ROCs
    //
    Double_t arr[kNsec];
    Int_t n=0;
    for (Int_t isec = 0; isec < kNsec; isec++) {
       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(AliTPCCalPad* outlierPad)
{
    //
    // return mean of the median of all ROCs
    //
    Double_t arr[kNsec];
    Int_t n=0;
    for (Int_t isec = 0; isec < kNsec; isec++) {
       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, AliTPCCalPad* outlierPad)
{
    //
    // return mean of the LTM and sigma of all ROCs
    //
    Double_t arrm[kNsec];
    Double_t arrs[kNsec];
    Double_t *sTemp=0x0;
    Int_t n=0;

    for (Int_t isec = 0; isec < kNsec; isec++) {
        AliTPCCalROC *calRoc = fROC[isec];
	if ( calRoc ){
	    if ( sigma ) sTemp=arrs+n;
       AliTPCCalROC* outlierROC = 0;
       if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
	    arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC);
            n++;
	}
    }
    if ( sigma ) *sigma = TMath::Mean(n,arrs);
    return TMath::Mean(n,arrm);
}

//_____________________________________________________________________________
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
      Float_t mean  = GetMean();
      Float_t sigma = GetRMS();
      Float_t nsigma = TMath::Abs(min);
      min = mean-nsigma*sigma;
      max = mean+nsigma*sigma;
    }
    if (type==1){
      // fixed range
      Float_t mean   = GetMedian();
      Float_t  delta = min;
      min = mean-delta;
      max = mean+delta;
    }
    if (type==2){
      //
      // LTM mean +- nsigma
      //
      Double_t sigma;
      Float_t mean  = GetLTM(&sigma,max);
      sigma*=min;
      min = mean-sigma;
      max = mean+sigma;
    }
  }
  char  name[1000];
  sprintf(name,"%s Pad 1D",GetTitle());
  TH1F * his = new TH1F(name,name,100, min,max);
    for (Int_t isec = 0; isec < kNsec; isec++) {
	if (fROC[isec]){
	    for (UInt_t irow=0; irow<fROC[isec]->GetNrows(); irow++){
		UInt_t npads = (Int_t)fROC[isec]->GetNPads(irow);
		for (UInt_t ipad=0; ipad<npads; ipad++){
		    his->Fill(fROC[isec]->GetValue(irow,ipad));
		}
	    }
	}
    }
  return his;
}

//_____________________________________________________________________________
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(); 
  for (Int_t isec=0; isec<72; isec++){
    if (side==0 && isec%36>=18) continue;
    if (side>0 && isec%36<18) continue;
    if (fROC[isec]){
      AliTPCCalROC * calRoc = fROC[isec];
      for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++)
	for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++)
	  if (TMath::Abs(calRoc->GetValue(irow,ipad))>kEpsilon){
	    Float_t xyz[3];
	    roc->GetPositionGlobal(isec,irow,ipad,xyz);
	    Int_t binx = 1+TMath::Nint((xyz[0]+250.)*0.5);
	    Int_t biny = 1+TMath::Nint((xyz[1]+250.)*0.5);
	    Float_t value = calRoc->GetValue(irow,ipad);	    
	    his->SetBinContent(binx,biny,value);
	  }
    }
  }
  his->SetXTitle("x (cm)");
  his->SetYTitle("y (cm)");
  return his;
}


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;
}


AliTPCCalPad* AliTPCCalPad::GlobalFit(const char* padName, AliTPCCalPad* PadOutliers, Bool_t robust, Int_t fitType, Double_t chi2Threshold, Double_t robustFraction){
   //
   // 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
   //
   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);
      else 
         GetCalROC(isec)->GlobalFit(0, robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction);
      pad->SetCalROC(AliTPCCalROC::CreateGlobalFitCalROC(fitParam, isec));
   }
   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);  
               }
            }
         }
      }
    }
  }
  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);  
               }
            }
         }
      }
    }
  }    
  
  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;
}
Commit	Line	Data
07627591	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	///////////////////////////////////////////////////////////////////////////////
	19	// //
a6d2bd0c	20	// TPC calibration class for parameters which are saved per pad //
a6d2bd0c	21	// Each AliTPCCalPad consists of 72 AliTPCCalROC-objects //
07627591	22	// //
	23	///////////////////////////////////////////////////////////////////////////////
	24
	25	#include "AliTPCCalPad.h"
	26	#include "AliTPCCalROC.h"
184bcc16	27	#include <TObjArray.h>
	28	#include <TAxis.h>
	29	#include <TGraph.h>
200be8a6	30	#include <TGraph2D.h>
200be8a6	31	#include <TH2F.h>
586007f3	32	#include "TTreeStream.h"
ca5dca67	33	#include "TFile.h"
ca5dca67	34	#include "TKey.h"
72d0ab7e	35	#include <iostream>
586007f3	36
07627591	37	ClassImp(AliTPCCalPad)
	38
	39	//_____________________________________________________________________________
	40	AliTPCCalPad::AliTPCCalPad():TNamed()
	41	{
	42	//
	43	// AliTPCCalPad default constructor
	44	//
	45
	46	for (Int_t isec = 0; isec < kNsec; isec++) {
	47	fROC[isec] = 0;
	48	}
	49
	50	}
	51
	52	//_____________________________________________________________________________
	53	AliTPCCalPad::AliTPCCalPad(const Text_t name, const Text_t title)
	54	:TNamed(name,title)
	55	{
	56	//
	57	// AliTPCCalPad constructor
	58	//
	59	for (Int_t isec = 0; isec < kNsec; isec++) {
	60	fROC[isec] = new AliTPCCalROC(isec);
	61	}
	62	}
	63
	64
	65	//_____________________________________________________________________________
	66	AliTPCCalPad::AliTPCCalPad(const AliTPCCalPad &c):TNamed(c)
	67	{
	68	//
	69	// AliTPCCalPad copy constructor
	70	//
	71
e6c51e6e	72	for (Int_t isec = 0; isec < kNsec; isec++) {
7fb8d357	73	fROC[isec] = 0;
	74	if (c.fROC[isec])
	75	fROC[isec] = new AliTPCCalROC(*(c.fROC[isec]));
e6c51e6e	76	}
07627591	77	}
07627591	78
184bcc16	79	//_____________________________________________________________________________
	80	AliTPCCalPad::AliTPCCalPad(TObjArray * array):TNamed()
	81	{
	82	//
	83	// AliTPCCalPad default constructor
	84	//
	85
	86	for (Int_t isec = 0; isec < kNsec; isec++) {
	87	fROC[isec] = (AliTPCCalROC *)array->At(isec);
	88	}
	89
	90	}
	91
	92
07627591	93	///_____________________________________________________________________________
	94	AliTPCCalPad::~AliTPCCalPad()
	95	{
	96	//
	97	// AliTPCCalPad destructor
	98	//
	99
	100	for (Int_t isec = 0; isec < kNsec; isec++) {
	101	if (fROC[isec]) {
	102	delete fROC[isec];
	103	fROC[isec] = 0;
	104	}
	105	}
	106
	107	}
	108
	109	//_____________________________________________________________________________
	110	AliTPCCalPad &AliTPCCalPad::operator=(const AliTPCCalPad &c)
	111	{
	112	//
	113	// Assignment operator
	114	//
	115
	116	if (this != &c) ((AliTPCCalPad &) c).Copy(*this);
	117	return *this;
	118
	119	}
	120
	121	//_____________________________________________________________________________
	122	void AliTPCCalPad::Copy(TObject &c) const
	123	{
	124	//
	125	// Copy function
	126	//
	127
	128	for (Int_t isec = 0; isec < kNsec; isec++) {
	129	if (fROC[isec]) {
	130	fROC[isec]->Copy(*((AliTPCCalPad &) c).fROC[isec]);
	131	}
	132	}
	133	TObject::Copy(c);
	134	}
184bcc16	135
a6d2bd0c	136
	137	void AliTPCCalPad::SetCalROC(AliTPCCalROC* roc, Int_t sector){
	138	//
	139	// Set AliTPCCalROC copies values from 'roc'
	140	// if sector == -1 the sector specified in 'roc' is used
	141	// else sector specified in 'roc' is ignored and specified sector is filled
	142	//
	143	if (sector == -1) sector = roc->GetSector();
72d0ab7e	144	if (!fROC[sector]) fROC[sector] = new AliTPCCalROC(sector);
a6d2bd0c	145	for (UInt_t ichannel = 0; ichannel < roc->GetNchannels(); ichannel++)
	146	fROC[sector]->SetValue(ichannel, roc->GetValue(ichannel));
	147	}
	148
	149
	150
90127643	151	//_____________________________________________________________________________
	152	void AliTPCCalPad::Add(Float_t c1)
	153	{
	154	//
a6d2bd0c	155	// add constant c1 to all channels of all ROCs
90127643	156	//
	157
	158	for (Int_t isec = 0; isec < kNsec; isec++) {
	159	if (fROC[isec]){
	160	fROC[isec]->Add(c1);
	161	}
	162	}
	163	}
	164
	165	//_____________________________________________________________________________
	166	void AliTPCCalPad::Multiply(Float_t c1)
	167	{
a6d2bd0c	168	//
	169	// multiply each channel of all ROCs with c1
	170	//
90127643	171	for (Int_t isec = 0; isec < kNsec; isec++) {
	172	if (fROC[isec]){
	173	fROC[isec]->Multiply(c1);
	174	}
	175	}
	176	}
	177
	178	//_____________________________________________________________________________
	179	void AliTPCCalPad::Add(const AliTPCCalPad * pad, Double_t c1)
	180	{
a6d2bd0c	181	//
	182	// multiply AliTPCCalPad 'pad' by c1 and add each channel to the coresponing channel in all ROCs
	183	// - pad by pad -
	184	//
90127643	185	for (Int_t isec = 0; isec < kNsec; isec++) {
	186	if (fROC[isec]){
	187	fROC[isec]->Add(pad->GetCalROC(isec),c1);
	188	}
	189	}
	190	}
	191
	192	//_____________________________________________________________________________
	193	void AliTPCCalPad::Multiply(const AliTPCCalPad * pad)
	194	{
a6d2bd0c	195	//
	196	// multiply each channel of all ROCs with the coresponding channel of 'pad'
	197	// - pad by pad -
	198	//
	199	for (Int_t isec = 0; isec < kNsec; isec++) {
90127643	200	if (fROC[isec]){
	201	fROC[isec]->Multiply(pad->GetCalROC(isec));
	202	}
	203	}
	204	}
	205
	206	//_____________________________________________________________________________
	207	void AliTPCCalPad::Divide(const AliTPCCalPad * pad)
	208	{
a6d2bd0c	209	//
	210	// divide each channel of all ROCs by the coresponding channel of 'pad'
	211	// - pad by pad -
	212	//
90127643	213	for (Int_t isec = 0; isec < kNsec; isec++) {
	214	if (fROC[isec]){
	215	fROC[isec]->Divide(pad->GetCalROC(isec));
	216	}
	217	}
	218	}
	219
	220	//_____________________________________________________________________________
184bcc16	221	TGraph * AliTPCCalPad::MakeGraph(Int_t type, Float_t ratio){
	222	//
	223	// type=1 - mean
	224	// 2 - median
	225	// 3 - LTM
a6d2bd0c	226	//
184bcc16	227	Int_t npoints = 0;
	228	for (Int_t i=0;i<72;i++) if (fROC[i]) npoints++;
	229	TGraph * graph = new TGraph(npoints);
	230	npoints=0;
	231	for (Int_t isec=0;isec<72;isec++){
	232	if (!fROC[isec]) continue;
	233	if (type==0) graph->SetPoint(npoints,isec,fROC[isec]->GetMean());
	234	if (type==1) graph->SetPoint(npoints,isec,fROC[isec]->GetMedian());
	235	if (type==2) graph->SetPoint(npoints,isec,fROC[isec]->GetLTM(0,ratio));
	236	npoints++;
	237	}
	238
	239	graph->GetXaxis()->SetTitle("Sector");
	240	if (type==0) {
	241	graph->GetYaxis()->SetTitle("Mean");
	242	graph->SetMarkerStyle(22);
	243	}
	244	if (type==1) {
	245	graph->GetYaxis()->SetTitle("Median");
	246	graph->SetMarkerStyle(22);
	247	}
	248	if (type==2) {
	249	graph->GetYaxis()->SetTitle(Form("Mean%f",ratio));
	250	graph->SetMarkerStyle(24);
	251	}
	252
	253	return graph;
	254	}
200be8a6	255
90127643	256	//_____________________________________________________________________________
	257	Double_t AliTPCCalPad::GetMeanRMS(Double_t &rms)
	258	{
	259	//
a6d2bd0c	260	// Calculates mean and RMS of all ROCs
90127643	261	//
	262	Double_t sum = 0, sum2 = 0, n=0, val=0;
	263	for (Int_t isec = 0; isec < kNsec; isec++) {
	264	AliTPCCalROC *calRoc = fROC[isec];
	265	if ( calRoc ){
	266	for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++){
	267	for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++){
	268	val = calRoc->GetValue(irow,ipad);
	269	sum+=val;
	270	sum2+=val*val;
	271	n++;
	272	}
	273	}
	274
	275	}
	276	}
	277	Double_t n1 = 1./n;
	278	Double_t mean = sum*n1;
	279	rms = TMath::Sqrt(TMath::Abs(sum2n1-meanmean));
	280	return mean;
	281	}
	282
	283
	284	//_____________________________________________________________________________
ca5dca67	285	Double_t AliTPCCalPad::GetMean(AliTPCCalPad* outlierPad)
90127643	286	{
	287	//
	288	// return mean of the mean of all ROCs
	289	//
	290	Double_t arr[kNsec];
	291	Int_t n=0;
	292	for (Int_t isec = 0; isec < kNsec; isec++) {
ca5dca67	293	AliTPCCalROC *calRoc = fROC[isec];
	294	if ( calRoc ){
	295	AliTPCCalROC* outlierROC = 0;
	296	if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
	297	arr[n] = calRoc->GetMean(outlierROC);
	298	n++;
	299	}
90127643	300	}
	301	return TMath::Mean(n,arr);
	302	}
	303
	304	//_____________________________________________________________________________
ca5dca67	305	Double_t AliTPCCalPad::GetRMS(AliTPCCalPad* outlierPad)
90127643	306	{
	307	//
	308	// return mean of the RMS of all ROCs
	309	//
	310	Double_t arr[kNsec];
	311	Int_t n=0;
	312	for (Int_t isec = 0; isec < kNsec; isec++) {
ca5dca67	313	AliTPCCalROC *calRoc = fROC[isec];
	314	if ( calRoc ){
	315	AliTPCCalROC* outlierROC = 0;
	316	if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
	317	arr[n] = calRoc->GetRMS(outlierROC);
	318	n++;
	319	}
90127643	320	}
	321	return TMath::Mean(n,arr);
	322	}
	323
	324	//_____________________________________________________________________________
ca5dca67	325	Double_t AliTPCCalPad::GetMedian(AliTPCCalPad* outlierPad)
90127643	326	{
	327	//
	328	// return mean of the median of all ROCs
	329	//
	330	Double_t arr[kNsec];
	331	Int_t n=0;
	332	for (Int_t isec = 0; isec < kNsec; isec++) {
ca5dca67	333	AliTPCCalROC *calRoc = fROC[isec];
	334	if ( calRoc ){
	335	AliTPCCalROC* outlierROC = 0;
	336	if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
	337	arr[n] = calRoc->GetMedian(outlierROC);
	338	n++;
	339	}
90127643	340	}
	341	return TMath::Mean(n,arr);
	342	}
	343
	344	//_____________________________________________________________________________
ca5dca67	345	Double_t AliTPCCalPad::GetLTM(Double_t sigma, Double_t fraction, AliTPCCalPad outlierPad)
90127643	346	{
	347	//
	348	// return mean of the LTM and sigma of all ROCs
	349	//
	350	Double_t arrm[kNsec];
	351	Double_t arrs[kNsec];
	352	Double_t *sTemp=0x0;
	353	Int_t n=0;
	354
	355	for (Int_t isec = 0; isec < kNsec; isec++) {
	356	AliTPCCalROC *calRoc = fROC[isec];
	357	if ( calRoc ){
	358	if ( sigma ) sTemp=arrs+n;
ca5dca67	359	AliTPCCalROC* outlierROC = 0;
	360	if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
	361	arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC);
90127643	362	n++;
	363	}
	364	}
	365	if ( sigma ) *sigma = TMath::Mean(n,arrs);
	366	return TMath::Mean(n,arrm);
	367	}
	368
	369	//_____________________________________________________________________________
	370	TH1F * AliTPCCalPad::MakeHisto1D(Float_t min, Float_t max,Int_t type){
	371	//
	372	// make 1D histo
	373	// type -1 = user defined range
	374	// 0 = nsigma cut nsigma=min
a6d2bd0c	375	//
90127643	376	if (type>=0){
	377	if (type==0){
	378	// nsigma range
	379	Float_t mean = GetMean();
	380	Float_t sigma = GetRMS();
	381	Float_t nsigma = TMath::Abs(min);
	382	min = mean-nsigma*sigma;
	383	max = mean+nsigma*sigma;
	384	}
	385	if (type==1){
	386	// fixed range
	387	Float_t mean = GetMedian();
	388	Float_t delta = min;
	389	min = mean-delta;
	390	max = mean+delta;
	391	}
	392	if (type==2){
	393	//
	394	// LTM mean +- nsigma
	395	//
	396	Double_t sigma;
	397	Float_t mean = GetLTM(&sigma,max);
	398	sigma*=min;
	399	min = mean-sigma;
	400	max = mean+sigma;
	401	}
	402	}
	403	char name[1000];
	404	sprintf(name,"%s Pad 1D",GetTitle());
	405	TH1F * his = new TH1F(name,name,100, min,max);
	406	for (Int_t isec = 0; isec < kNsec; isec++) {
	407	if (fROC[isec]){
	408	for (UInt_t irow=0; irow<fROC[isec]->GetNrows(); irow++){
	409	UInt_t npads = (Int_t)fROC[isec]->GetNPads(irow);
	410	for (UInt_t ipad=0; ipad<npads; ipad++){
	411	his->Fill(fROC[isec]->GetValue(irow,ipad));
	412	}
	413	}
	414	}
	415	}
	416	return his;
	417	}
	418
	419	//_____________________________________________________________________________
200be8a6	420	TH2F *AliTPCCalPad::MakeHisto2D(Int_t side){
	421	//
	422	// Make 2D graph
	423	// side - specify the side A = 0 C = 1
	424	// type - used types of determination of boundaries in z
a6d2bd0c	425	//
200be8a6	426	Float_t kEpsilon = 0.000000000001;
	427	TH2F * his = new TH2F(GetName(), GetName(), 250,-250,250,250,-250,250);
	428	AliTPCROC * roc = AliTPCROC::Instance();
	429	for (Int_t isec=0; isec<72; isec++){
	430	if (side==0 && isec%36>=18) continue;
	431	if (side>0 && isec%36<18) continue;
	432	if (fROC[isec]){
	433	AliTPCCalROC * calRoc = fROC[isec];
	434	for (UInt_t irow=0; irow<calRoc->GetNrows(); irow++)
	435	for (UInt_t ipad=0; ipad<calRoc->GetNPads(irow); ipad++)
	436	if (TMath::Abs(calRoc->GetValue(irow,ipad))>kEpsilon){
	437	Float_t xyz[3];
	438	roc->GetPositionGlobal(isec,irow,ipad,xyz);
	439	Int_t binx = 1+TMath::Nint((xyz[0]+250.)*0.5);
	440	Int_t biny = 1+TMath::Nint((xyz[1]+250.)*0.5);
	441	Float_t value = calRoc->GetValue(irow,ipad);
	442	his->SetBinContent(binx,biny,value);
	443	}
	444	}
	445	}
	446	his->SetXTitle("x (cm)");
	447	his->SetYTitle("y (cm)");
	448	return his;
	449	}
	450
	451
72d0ab7e	452	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 {
a6d2bd0c	453	//
	454	// Loops over all AliTPCCalROCs and performs a localFit in each ROC
	455	// AliTPCCalPad with fit-data is returned
	456	// rowRadius and padRadius specifies a window around a given pad in one sector.
	457	// The data of this window are fitted with a parabolic function.
	458	// This function is evaluated at the pad's position.
	459	// At the edges the window is shifted, so that the specified pad is not anymore in the center of the window.
	460	// rowRadius - radius - rows to be used for smoothing
	461	// padradius - radius - pads to be used for smoothing
	462	// ROCoutlier - map of outliers - pads not to be used for local smoothing
	463	// robust - robust method of fitting - (much slower)
	464	// chi2Threshold: Threshold for chi2 when EvalRobust is called
	465	// robustFraction: Fraction of data that will be used in EvalRobust
	466	//
	467	//
	468	AliTPCCalPad* pad = new AliTPCCalPad(padName, padName);
	469	for (Int_t isec = 0; isec < 72; isec++){
72d0ab7e	470	if (printCurrentSector) std::cout << "LocalFit in sector " << isec << "\r" << std::flush;
a6d2bd0c	471	if (PadOutliers)
72d0ab7e	472	pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, PadOutliers->GetCalROC(isec), robust, chi2Threshold, robustFraction));
a6d2bd0c	473	else
72d0ab7e	474	pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, 0, robust, chi2Threshold, robustFraction));
a6d2bd0c	475	}
	476	return pad;
	477	}
	478
	479
	480	AliTPCCalPad* AliTPCCalPad::GlobalFit(const char* padName, AliTPCCalPad* PadOutliers, Bool_t robust, Int_t fitType, Double_t chi2Threshold, Double_t robustFraction){
	481	//
	482	// Loops over all AliTPCCalROCs and performs a globalFit in each ROC
	483	// AliTPCCalPad with fit-data is returned
	484	// chi2Threshold: Threshold for chi2 when EvalRobust is called
	485	// robustFraction: Fraction of data that will be used in EvalRobust
	486	// chi2Threshold: Threshold for chi2 when EvalRobust is called
	487	// robustFraction: Fraction of data that will be used in EvalRobust
	488	//
	489	AliTPCCalPad* pad = new AliTPCCalPad(padName, padName);
	490	TVectorD fitParam(0);
	491	TMatrixD covMatrix(0,0);
	492	Float_t chi2 = 0;
	493	for (Int_t isec = 0; isec < 72; isec++){
	494	if (PadOutliers)
	495	GetCalROC(isec)->GlobalFit(PadOutliers->GetCalROC(isec), robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction);
	496	else
	497	GetCalROC(isec)->GlobalFit(0, robust, fitParam, covMatrix, chi2, fitType, chi2Threshold, robustFraction);
	498	pad->SetCalROC(AliTPCCalROC::CreateGlobalFitCalROC(fitParam, isec));
	499	}
	500	return pad;
	501	}
	502
	503
	504	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){
	505	//
	506	// Makes a GlobalFit over each side and return fit-parameters, covariance and chi2 for each side
	507	// fitType == 0: fit plane function
	508	// fitType == 1: fit parabolic function
	509	// PadOutliers - pads with value !=0 are not used in fitting procedure
	510	// chi2Threshold: Threshold for chi2 when EvalRobust is called
	511	// robustFraction: Fraction of data that will be used in EvalRobust
	512	//
	513	TLinearFitter* fitterGA = 0;
	514	TLinearFitter* fitterGC = 0;
	515
	516	if (fitType == 1) {
	517	fitterGA = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
	518	fitterGC = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
	519	}
	520	else {
	521	fitterGA = new TLinearFitter(3,"x0++x1++x2");
	522	fitterGC = new TLinearFitter(3,"x0++x1++x2");
	523	}
	524	fitterGA->StoreData(kTRUE);
	525	fitterGC->StoreData(kTRUE);
	526	fitterGA->ClearPoints();
	527	fitterGC->ClearPoints();
	528	Double_t xx[6];
	529	Int_t npointsA=0;
	530	Int_t npointsC=0;
	531
	532	Float_t localXY[3] = {0}; // pad's position, needed to get the pad's position
	533	Float_t lx, ly; // pads position
	534
	535	AliTPCROC* tpcROCinstance = AliTPCROC::Instance(); // to calculate the pad's position
	536
	537	// loop over all sectors and pads and read data into fitterGA and fitterGC
	538	if (fitType == 1) {
539	// parabolic fit
540	fitParamSideA.ResizeTo(6);
541	fitParamSideC.ResizeTo(6);
542	covMatrixSideA.ResizeTo(6,6);
543	covMatrixSideC.ResizeTo(6,6);
544	for (UInt_t isec = 0; isec<72; isec++){
545	for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
546	for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
547	// fill fitterG
548	tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY); // calculate position localXY by sector, pad and row number
549	lx = localXY[0];
550	ly = localXY[1];
551	xx[0] = 1;
552	xx[1] = lx;
553	xx[2] = ly;
554	xx[3] = lx*lx;
555	xx[4] = ly*ly;
556	xx[5] = lx*ly;
557	if (!PadOutliers \|\| PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) {
558	// if given pad is no outlier, add it to TLinearFitter, decide to which of both
559	// sector 0 - 17: IROC, A
560	// sector 18 - 35: IROC, C
561	// sector 36 - 53: OROC, A
562	// sector 54 - 71: CROC, C
563	if (isec <= 17 \|\| (isec >= 36 && isec <= 53)) { // Side A
564	npointsA++;
565	fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
566	}
567	else { // side C
568	npointsC++;
569	fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
570	}
571	}
572	}
573	}
574	}
575	}
576	else {
577	// linear fit
578	fitParamSideA.ResizeTo(3);
579	fitParamSideC.ResizeTo(3);
580	covMatrixSideA.ResizeTo(3,3);
581	covMatrixSideC.ResizeTo(3,3);
582
583	for (UInt_t isec = 0; isec<72; isec++){
584	for (UInt_t irow = 0; irow < GetCalROC(isec)->GetNrows(); irow++) {
585	for (UInt_t ipad = 0; ipad < GetCalROC(isec)->GetNPads(irow); ipad++) {
586	// fill fitterG
587	tpcROCinstance->GetPositionLocal(isec, irow, ipad, localXY); // calculate position localXY by sector, pad and row number
588	lx = localXY[0];
589	ly = localXY[1];
590	xx[0] = 1;
591	xx[1] = lx;
592	xx[2] = ly;
593	if (!PadOutliers \|\| PadOutliers->GetCalROC(isec)->GetValue(irow, ipad) != 1) {
594	// if given pad is no outlier, add it to TLinearFitter, decide to which of both
595	// sector 0 - 17: IROC, A
596	// sector 18 - 35: IROC, C
597	// sector 36 - 53: OROC, A
598	// sector 54 - 71: CROC, C
599	if (isec <= 17 \|\| (isec >= 36 && isec <= 53)) {
600	// Side A
601	npointsA++;
602	fitterGA->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
603	}
604	else {
605	// side C
606	npointsC++;
607	fitterGC->AddPoint(xx, GetCalROC(isec)->GetValue(irow, ipad), 1);
608	}
609	}
610	}
611	}
612	}
613	}
614
615	fitterGA->Eval();
616	fitterGC->Eval();
617	fitterGA->GetParameters(fitParamSideA);
618	fitterGC->GetParameters(fitParamSideC);
619	fitterGA->GetCovarianceMatrix(covMatrixSideA);
620	fitterGC->GetCovarianceMatrix(covMatrixSideC);
621	if (fitType == 1){
622	chi2SideA = fitterGA->GetChisquare()/(npointsA-6.);
623	chi2SideC = fitterGC->GetChisquare()/(npointsC-6.);
624	}
625	else {
626	chi2SideA = fitterGA->GetChisquare()/(npointsA-3.);
627	chi2SideC = fitterGC->GetChisquare()/(npointsC-3.);
628	}
629	if (robust && chi2SideA > chi2Threshold) {
630	// std::cout << "robust fitter called... " << std::endl;
631	fitterGA->EvalRobust(robustFraction);
632	fitterGA->GetParameters(fitParamSideA);
633	}
634	if (robust && chi2SideC > chi2Threshold) {
635	// std::cout << "robust fitter called... " << std::endl;
636	fitterGC->EvalRobust(robustFraction);
637	fitterGC->GetParameters(fitParamSideC);
638	}
639	delete fitterGA;
640	delete fitterGC;
641	}
200be8a6	642
586007f3	643