TPC/AliTPCCalPad.cxx

   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 //                                                                           //
  20 //  TPC calibration class for parameters which are saved per pad             //
  21 //  Each AliTPCCalPad consists of 72 AliTPCCalROC-objects                    //
  22 //                                                                           //
  23 ///////////////////////////////////////////////////////////////////////////////
  24
  25 #include "AliTPCCalPad.h"
  26 #include "AliTPCCalROC.h"
  27 #include <TObjArray.h>
  28 #include <TAxis.h>
  29 #include <TGraph.h>
  30 #include <TGraph2D.h>
  31 #include <TH2F.h>
  32 #include "TTreeStream.h"
  33 #include "TFile.h"
  34 #include "TKey.h"
  35 #include <iostream>
  36
  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
  72   for (Int_t isec = 0; isec < kNsec; isec++) {
  73          fROC[isec] = 0;
  74      if (c.fROC[isec])
  75        fROC[isec] = new AliTPCCalROC(*(c.fROC[isec]));
  76   }
  77 }
  78
  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
  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 }
 135
 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();
 144    if (!fROC[sector]) fROC[sector] = new AliTPCCalROC(sector);
 145    for (UInt_t ichannel = 0; ichannel < roc->GetNchannels(); ichannel++)
 146       fROC[sector]->SetValue(ichannel, roc->GetValue(ichannel));
 147 }
 148
 149
 150
 151 //_____________________________________________________________________________
 152 void AliTPCCalPad::Add(Float_t c1)
 153 {
 154     //
 155     // add constant c1 to all channels of all ROCs
 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 {
 168   //
 169   // multiply each channel of all ROCs with c1
 170   //
 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 {
 181   //
 182   // multiply AliTPCCalPad 'pad' by c1 and add each channel to the coresponing channel in all ROCs
 183   //  - pad by pad -
 184   //
 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 {
 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++) {
 200         if (fROC[isec]){
 201             fROC[isec]->Multiply(pad->GetCalROC(isec));
 202         }
 203     }
 204 }
 205
 206 //_____________________________________________________________________________
 207 void AliTPCCalPad::Divide(const AliTPCCalPad * pad)
 208 {
 209   //
 210   // divide each channel of all ROCs by the coresponding channel of 'pad'
 211   //     - pad by pad -
 212   //
 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 //_____________________________________________________________________________
 221 TGraph  *  AliTPCCalPad::MakeGraph(Int_t type, Float_t ratio){
 222   //
 223   //   type=1 - mean
 224   //        2 - median
 225   //        3 - LTM
 226   //
 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 }
 255
 256 //_____________________________________________________________________________
 257 Double_t AliTPCCalPad::GetMeanRMS(Double_t &rms)
 258 {
 259     //
 260     // Calculates mean and RMS of all ROCs
 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(sum2*n1-mean*mean));
 280     return mean;
 281 }
 282
 283
 284 //_____________________________________________________________________________
 285 Double_t AliTPCCalPad::GetMean(AliTPCCalPad* outlierPad)
 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++) {
 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        }
 300     }
 301     return TMath::Mean(n,arr);
 302 }
 303
 304 //_____________________________________________________________________________
 305 Double_t AliTPCCalPad::GetRMS(AliTPCCalPad* outlierPad)
 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++) {
 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        }
 320     }
 321     return TMath::Mean(n,arr);
 322 }
 323
 324 //_____________________________________________________________________________
 325 Double_t AliTPCCalPad::GetMedian(AliTPCCalPad* outlierPad)
 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++) {
 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        }
 340     }
 341     return TMath::Mean(n,arr);
 342 }
 343
 344 //_____________________________________________________________________________
 345 Double_t AliTPCCalPad::GetLTM(Double_t *sigma, Double_t fraction, AliTPCCalPad* outlierPad)
 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;
 359        AliTPCCalROC* outlierROC = 0;
 360        if (outlierPad) outlierROC = outlierPad->GetCalROC(isec);
 361             arrm[n] = calRoc->GetLTM(sTemp,fraction, outlierROC);
 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
 375   //
 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 //_____________________________________________________________________________
 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
 425   //
 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
 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 {
 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++){
 470       if (printCurrentSector) std::cout << "LocalFit in sector " << isec << "\r" << std::flush;
 471       if (PadOutliers)
 472          pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, PadOutliers->GetCalROC(isec), robust, chi2Threshold, robustFraction));
 473       else
 474          pad->SetCalROC(GetCalROC(isec)->LocalFit(rowRadius, padRadius, 0, robust, chi2Threshold, robustFraction));
 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 }
 642
 643