Corrected UInt_t <-> Int_t conversion
[u/mrichter/AliRoot.git] / TPC / AliTPCCalROC.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
17 ///////////////////////////////////////////////////////////////////////////////
18 //                                                                           //
19 //     Calibration base class for a single ROC                               //
20 //     Contains one float value per pad                                      //
21 //     mapping of the pads taken form AliTPCROC                              //
22 //                                                                           //
23 ///////////////////////////////////////////////////////////////////////////////
24
25 //
26 // ROOT includes 
27 //
28 #include "TMath.h"
29 #include "TClass.h"
30 #include "TFile.h"
31 #include "TH1F.h"
32 #include "TH2F.h"
33 #include "TArrayI.h"
34
35 //
36 //
37 #include "AliTPCCalROC.h"
38 #include "AliMathBase.h"
39
40 #include "TRandom3.h"      // only needed by the AliTPCCalROCTest() method
41
42 ClassImp(AliTPCCalROC)
43
44
45 //_____________________________________________________________________________
46 AliTPCCalROC::AliTPCCalROC()
47              :TNamed(),
48               fSector(0),
49               fNChannels(0),
50               fNRows(0),
51               fkIndexes(0),
52               fData(0)
53 {
54   //
55   // Default constructor
56   //
57
58 }
59
60 //_____________________________________________________________________________
61 AliTPCCalROC::AliTPCCalROC(UInt_t  sector)
62              :TNamed(),
63               fSector(0),
64               fNChannels(0),
65               fNRows(0),
66               fkIndexes(0),
67               fData(0)
68 {
69   //
70   // Constructor that initializes a given sector
71   //
72   fSector = sector;
73   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
74   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
75   fkIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
76   fData = new Float_t[fNChannels];
77   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = 0.;
78 }
79
80 //_____________________________________________________________________________
81 AliTPCCalROC::AliTPCCalROC(const AliTPCCalROC &c)
82              :TNamed(c),
83               fSector(0),
84               fNChannels(0),
85               fNRows(0),
86               fkIndexes(0),
87               fData(0)
88 {
89   //
90   // AliTPCCalROC copy constructor
91   //
92   fSector = c.fSector;
93   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
94   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
95   fkIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
96   //
97   fData   = new Float_t[fNChannels];
98   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = c.fData[idata];
99 }
100 //____________________________________________________________________________
101 AliTPCCalROC & AliTPCCalROC::operator =(const AliTPCCalROC & param)
102 {
103   //
104   // assignment operator - dummy
105   //
106   fData=param.fData;
107   return (*this);
108 }
109
110
111 //_____________________________________________________________________________
112 AliTPCCalROC::~AliTPCCalROC()
113 {
114   //
115   // AliTPCCalROC destructor
116   //
117   if (fData) {
118     delete [] fData;
119     fData = 0;
120   }
121 }
122
123
124
125 void AliTPCCalROC::Streamer(TBuffer &R__b)
126 {
127    //
128    // Stream an object of class AliTPCCalROC.
129    //
130    if (R__b.IsReading()) {
131       AliTPCCalROC::Class()->ReadBuffer(R__b, this);
132       fkIndexes =  AliTPCROC::Instance()->GetRowIndexes(fSector);
133    } else {
134       AliTPCCalROC::Class()->WriteBuffer(R__b,this);
135    }
136 }
137
138
139 // algebra fuctions:
140
141 void AliTPCCalROC::Add(Float_t c1){
142   //
143   // add c1 to each channel of the ROC  
144   //
145   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata]+=c1;
146 }
147
148
149 void AliTPCCalROC::Multiply(Float_t c1){
150   //
151   // multiply each channel of the ROC with c1
152   //
153   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata]*=c1;
154 }
155
156
157 void AliTPCCalROC::Add(const AliTPCCalROC * roc, Double_t c1){
158   //
159   // multiply AliTPCCalROC roc by c1 and add each channel to the coresponing channel in the ROC
160   //  - pad by pad 
161   //
162   for (UInt_t  idata = 0; idata< fNChannels; idata++){
163     fData[idata]+=roc->fData[idata]*c1;
164   }
165 }
166
167
168 void AliTPCCalROC::Multiply(const AliTPCCalROC*  roc) {
169   //
170   // multiply each channel of the ROC with the coresponding channel of 'roc'
171   //     - pad by pad -
172   //
173   for (UInt_t  idata = 0; idata< fNChannels; idata++){
174     fData[idata]*=roc->fData[idata];
175   }
176 }
177
178
179 void AliTPCCalROC::Divide(const AliTPCCalROC*  roc) {
180   //
181   // divide each channel of the ROC by the coresponding channel of 'roc'
182   //     - pad by pad -
183   //
184   Float_t kEpsilon=0.00000000000000001;
185   for (UInt_t  idata = 0; idata< fNChannels; idata++){
186     if (TMath::Abs(roc->fData[idata])>kEpsilon)
187       fData[idata]/=roc->fData[idata];
188   }
189 }
190
191 Double_t AliTPCCalROC::GetMean(AliTPCCalROC *const outlierROC) const {
192    //
193    //  returns the mean value of the ROC
194    //  pads with value != 0 in outlierROC are not used for the calculation
195    //  return 0 if no data is accepted by the outlier cuts 
196    //
197    if (!outlierROC) return TMath::Mean(fNChannels, fData);
198    Double_t *ddata = new Double_t[fNChannels];
199    Int_t nPoints = 0;
200    for (UInt_t i=0;i<fNChannels;i++) {
201       if (!(outlierROC->GetValue(i))) {
202          ddata[nPoints]= fData[i];
203          nPoints++;
204       }
205    }
206    Double_t mean = 0;
207    if(nPoints>0)
208      mean = TMath::Mean(nPoints, ddata);
209    delete [] ddata;
210    return mean;
211 }
212
213 Double_t AliTPCCalROC::GetMedian(AliTPCCalROC *const outlierROC) const {
214    //
215    //  returns the median value of the ROC
216    //  pads with value != 0 in outlierROC are not used for the calculation
217    //  return 0 if no data is accepted by the outlier cuts 
218    //
219    if (!outlierROC) return TMath::Median(fNChannels, fData);
220    Double_t *ddata = new Double_t[fNChannels];
221    Int_t nPoints = 0;
222    for (UInt_t i=0;i<fNChannels;i++) {
223       if (!(outlierROC->GetValue(i))) {
224          ddata[nPoints]= fData[i];
225          nPoints++;
226       }
227    }
228    Double_t median = 0;
229    if(nPoints>0)
230      median = TMath::Median(nPoints, ddata);
231    delete [] ddata;
232    return median;
233 }
234
235 Double_t AliTPCCalROC::GetRMS(AliTPCCalROC *const outlierROC) const {
236    //
237    //  returns the RMS value of the ROC
238    //  pads with value != 0 in outlierROC are not used for the calculation
239    //  return 0 if no data is accepted by the outlier cuts 
240    //
241    if (!outlierROC) return TMath::RMS(fNChannels, fData);
242    Double_t *ddata = new Double_t[fNChannels];
243    Int_t nPoints = 0;
244    for (UInt_t i=0;i<fNChannels;i++) {
245       if (!(outlierROC->GetValue(i))) {
246          ddata[nPoints]= fData[i];
247          nPoints++;
248       }
249    }
250    Double_t rms = 0;
251    if(nPoints>0)
252      rms = TMath::RMS(nPoints, ddata);
253    delete [] ddata;
254    return rms;
255 }
256
257 Double_t AliTPCCalROC::GetLTM(Double_t *const sigma, Double_t fraction, AliTPCCalROC *const outlierROC){
258   //
259   //  returns the LTM and sigma
260   //  pads with value != 0 in outlierROC are not used for the calculation
261    //  return 0 if no data is accepted by the outlier cuts 
262   //
263   Double_t *ddata = new Double_t[fNChannels];
264   UInt_t nPoints = 0;
265   for (UInt_t i=0;i<fNChannels;i++) {
266      if (!outlierROC || !(outlierROC->GetValue(i))) {
267         ddata[nPoints]= fData[i];
268         nPoints++;
269      }
270   }
271
272   Double_t ltm =0, lsigma=0;
273   if(nPoints>0) {
274     Int_t hh = TMath::Min(TMath::Nint(fraction *nPoints), Int_t(nPoints));
275     AliMathBase::EvaluateUni(nPoints,ddata, ltm, lsigma, hh);
276     if (sigma) *sigma=lsigma;
277   }
278   
279   delete [] ddata;
280   return ltm;
281 }
282
283 TH1F * AliTPCCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type){
284   //
285   // make 1D histo
286   // type -1 = user defined range
287   //       0 = nsigma cut nsigma=min
288   //       1 = delta cut around median delta=min
289   //
290   if (type>=0){
291     if (type==0){
292       // nsigma range
293       Float_t mean  = GetMean();
294       Float_t sigma = GetRMS();
295       Float_t nsigma = TMath::Abs(min);
296       min = mean-nsigma*sigma;
297       max = mean+nsigma*sigma;
298     }
299     if (type==1){
300       // fixed range
301       Float_t mean   = GetMedian();
302       Float_t  delta = min;
303       min = mean-delta;
304       max = mean+delta;
305     }
306     if (type==2){
307       //
308       // LTM mean +- nsigma
309       //
310       Double_t sigma;
311       Float_t mean  = GetLTM(&sigma,max);
312       sigma*=min;
313       min = mean-sigma;
314       max = mean+sigma;
315     }
316   }
317   TString name=Form("%s ROC 1D%d",GetTitle(),fSector);
318   TH1F * his = new TH1F(name.Data(),name.Data(),100, min,max);
319   for (UInt_t irow=0; irow<fNRows; irow++){
320     UInt_t npads = (Int_t)GetNPads(irow);
321     for (UInt_t ipad=0; ipad<npads; ipad++){
322       his->Fill(GetValue(irow,ipad));
323     }
324   }
325   return his;
326 }
327
328
329
330 TH2F * AliTPCCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type){
331   //
332   // make 2D histo
333   // type -1 = user defined range
334   //       0 = nsigma cut nsigma=min
335   //       1 = delta cut around median delta=min
336   //
337   if (type>=0){
338     if (type==0){
339       // nsigma range
340       Float_t mean  = GetMean();
341       Float_t sigma = GetRMS();
342       Float_t nsigma = TMath::Abs(min);
343       min = mean-nsigma*sigma;
344       max = mean+nsigma*sigma;
345     }
346     if (type==1){
347       // fixed range
348       Float_t mean   = GetMedian();
349       Float_t  delta = min;
350       min = mean-delta;
351       max = mean+delta;
352     }
353     if (type==2){
354       Double_t sigma;
355       Float_t mean  = GetLTM(&sigma,max);
356       sigma*=min;
357       min = mean-sigma;
358       max = mean+sigma;
359
360     }
361   }
362   UInt_t maxPad = 0;
363   for (UInt_t irow=0; irow<fNRows; irow++){
364     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
365   }
366
367   TString name=Form("%s ROC%d",GetTitle(),fSector);
368   TH2F * his = new TH2F(name.Data(),name.Data(),fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
369   for (UInt_t irow=0; irow<fNRows; irow++){
370     UInt_t npads = (Int_t)GetNPads(irow);
371     for (UInt_t ipad=0; ipad<npads; ipad++){
372       his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,GetValue(irow,ipad));
373     }
374   }
375   his->SetMaximum(max);
376   his->SetMinimum(min);
377   return his;
378 }
379
380 TH2F * AliTPCCalROC::MakeHistoOutliers(Float_t delta, Float_t fraction, Int_t type){
381   //
382   // Make Histogram with outliers
383   // mode = 0 - sigma cut used
384   // mode = 1 - absolute cut used
385   // fraction - fraction of values used to define sigma
386   // delta - in mode 0 - nsigma cut
387   //            mode 1 - delta cut
388   //
389   Double_t sigma;
390   Float_t mean  = GetLTM(&sigma,fraction);  
391   if (type==0) delta*=sigma; 
392   UInt_t maxPad = 0;
393   for (UInt_t irow=0; irow<fNRows; irow++){
394     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
395   }
396
397   TString name=Form("%s ROC Outliers%d",GetTitle(),fSector);
398   TH2F * his = new TH2F(name.Data(),name.Data(),fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
399   for (UInt_t irow=0; irow<fNRows; irow++){
400     UInt_t npads = (Int_t)GetNPads(irow);
401     for (UInt_t ipad=0; ipad<npads; ipad++){
402       if (TMath::Abs(GetValue(irow,ipad)-mean)>delta)
403         his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,1);
404     }
405   }
406   return his;
407 }
408
409
410
411 void AliTPCCalROC::Draw(Option_t* opt){
412   //
413   // create histogram with values and draw it
414   //
415   TH1 * his=0; 
416   TString option=opt;
417   option.ToUpper();
418   if (option.Contains("1D")){
419     his = MakeHisto1D();
420   }
421   else{
422     his = MakeHisto2D();
423   }
424   his->Draw(option);
425 }
426
427
428
429
430
431 void AliTPCCalROC::Test() {
432   //
433   // example function to show functionality and test AliTPCCalROC
434   //
435
436   Float_t kEpsilon=0.00001;
437   
438   // create CalROC with defined values
439   AliTPCCalROC  roc0(0);  
440   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
441     for (UInt_t ipad = 0; ipad <roc0.GetNPads(irow); ipad++){
442       Float_t value  = irow+ipad/1000.;
443       roc0.SetValue(irow,ipad,value);
444     }
445   }
446   
447   // copy CalROC, readout values and compare to original
448   AliTPCCalROC roc1(roc0);
449   for (UInt_t irow = 0; irow <roc1.GetNrows(); irow++){
450     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
451       Float_t value  = irow+ipad/1000.;
452       if (roc1.GetValue(irow,ipad)!=value){
453         printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
454       }
455     }
456   }
457
458   // write original CalROC to file
459   TFile f("calcTest.root","recreate");
460   roc0.Write("Roc0");
461   AliTPCCalROC * roc2 = (AliTPCCalROC*)f.Get("Roc0");
462   f.Close();
463   
464   // read CalROC from file and compare to original CalROC
465   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
466     if (roc0.GetNPads(irow)!=roc2->GetNPads(irow))
467       printf("NPads - Read/Write error\trow=%d\n",irow);
468     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
469       Float_t value  = irow+ipad/1000.;
470       if (roc2->GetValue(irow,ipad)!=value){
471         printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
472       }
473     }
474   }
475
476   // 
477   // Algebra Tests
478   //
479   
480   // Add constant
481   AliTPCCalROC roc3(roc0);
482   roc3.Add(1.5);
483   for (UInt_t irow = 0; irow <roc3.GetNrows(); irow++){
484     for (UInt_t ipad = 0; ipad <roc3.GetNPads(irow); ipad++){
485       Float_t value  = irow+ipad/1000. + 1.5;
486       if (TMath::Abs(roc3.GetValue(irow,ipad)-value) > kEpsilon){
487         printf("Add constant - error\trow=%d\tpad=%d\n",irow,ipad);
488       }
489     }
490   }
491
492   // Add another CalROC
493   AliTPCCalROC roc4(roc0);
494   roc4.Add(&roc0, -1.5);
495   for (UInt_t irow = 0; irow <roc4.GetNrows(); irow++){
496     for (UInt_t ipad = 0; ipad <roc4.GetNPads(irow); ipad++){
497       Float_t value  = irow+ipad/1000. - 1.5 * (irow+ipad/1000.);
498       if (TMath::Abs(roc4.GetValue(irow,ipad)-value) > kEpsilon){
499         printf("Add CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
500       }
501     }
502   }
503   
504   // Multiply with constant
505   AliTPCCalROC roc5(roc0);
506   roc5.Multiply(-1.4);
507   for (UInt_t irow = 0; irow <roc5.GetNrows(); irow++){
508     for (UInt_t ipad = 0; ipad <roc5.GetNPads(irow); ipad++){
509       Float_t value  = (irow+ipad/1000.) * (-1.4);
510       if (TMath::Abs(roc5.GetValue(irow,ipad)-value) > kEpsilon){
511         printf("Multiply with constant - error\trow=%d\tpad=%d\n",irow,ipad);
512       }
513     }
514   }
515
516   // Multiply another CalROC
517   AliTPCCalROC roc6(roc0);
518   roc6.Multiply(&roc0);
519   for (UInt_t irow = 0; irow <roc6.GetNrows(); irow++){
520     for (UInt_t ipad = 0; ipad <roc6.GetNPads(irow); ipad++){
521       Float_t value  = (irow+ipad/1000.) * (irow+ipad/1000.);
522       if (TMath::Abs(roc6.GetValue(irow,ipad)-value) > kEpsilon*100){
523         printf("Multiply with CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
524       }
525     }
526   }
527
528
529   // Divide by CalROC
530   AliTPCCalROC roc7(roc0);
531   roc7.Divide(&roc0);
532   for (UInt_t irow = 0; irow <roc7.GetNrows(); irow++){
533     for (UInt_t ipad = 0; ipad <roc7.GetNPads(irow); ipad++){
534       Float_t value  = 1.;
535       if (irow+ipad == 0) value = roc0.GetValue(irow,ipad);
536       if (TMath::Abs(roc7.GetValue(irow,ipad)-value) > kEpsilon){
537         printf("Multiply with CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
538       }
539     }
540   }
541
542   //
543   // Statistics Test
544   //
545   
546   // create CalROC with defined values
547   TRandom3 rnd(0);
548   AliTPCCalROC sroc0(0);
549   for (UInt_t ichannel = 0; ichannel < sroc0.GetNchannels(); ichannel++){
550     Float_t value  = rnd.Gaus(10., 2.);
551     sroc0.SetValue(ichannel,value);
552   }
553
554   printf("Mean   (should be close to 10): %f\n", sroc0.GetMean());
555   printf("RMS    (should be close to  2): %f\n", sroc0.GetRMS());
556   printf("Median (should be close to 10): %f\n", sroc0.GetMedian());
557   printf("LTM    (should be close to 10): %f\n", sroc0.GetLTM());
558
559   //AliTPCCalROC* sroc1 = sroc0.LocalFit(4, 8);
560   
561   //delete sroc1;
562
563 //        std::cout << TMath::Abs(roc5.GetValue(irow,ipad)-value) << std::endl;
564 }
565
566
567 AliTPCCalROC * AliTPCCalROC::LocalFit(Int_t rowRadius, Int_t padRadius, AliTPCCalROC* ROCoutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction) {
568   //
569   // MakeLocalFit - smoothing
570   // returns a AliTPCCalROC with smoothed data
571   // rowRadius and padRadius specifies a window around a given pad. 
572   // The data of this window are fitted with a parabolic function. 
573   // This function is evaluated at the pad's position.
574   // At the edges the window is shifted, so that the specified pad is not anymore in the center of the window. 
575   // rowRadius  -  radius - rows to be used for smoothing
576   // padradius  -  radius - pads to be used for smoothing
577   // ROCoutlier -  map of outliers - pads not to be used for local smoothing
578   // robust     -  robust method of fitting  - (much slower)
579   // chi2Threshold: Threshold for chi2 when EvalRobust is called
580   // robustFraction: Fraction of data that will be used in EvalRobust
581   //
582   AliTPCCalROC * xROCfitted = new AliTPCCalROC(fSector);
583   TLinearFitter fitterQ(6,"hyp5");
584    // TLinearFitter fitterQ(6,"x0++x1++x2++x3++x4++x5");
585   fitterQ.StoreData(kTRUE);
586   for (UInt_t row=0; row < GetNrows(); row++) {
587     //std::cout << "Entering row " << row << " of " << GetNrows() << " @ sector "<< fSector << " for local fitting... "<< std::endl;
588     for (UInt_t pad=0; pad < GetNPads(row); pad++)
589       xROCfitted->SetValue(row, pad, GetNeighbourhoodValue(&fitterQ, row, pad, rowRadius, padRadius, ROCoutliers, robust, chi2Threshold, robustFraction));
590   }
591   return xROCfitted;
592 }
593
594
595 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) {
596   //
597   // AliTPCCalROC::GetNeighbourhoodValue - smoothing - PRIVATE
598   // in this function the fit for LocalFit is done
599   //
600
601   fitterQ->ClearPoints();
602   TVectorD fitParam(6);
603   Int_t    npoints=0;
604   Double_t xx[6];
605   Float_t dlx, dly;
606   Float_t lPad[3] = {0};
607   Float_t localXY[3] = {0};
608   
609   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
610   tpcROCinstance->GetPositionLocal(fSector, row, pad, lPad);  // calculate position lPad by pad and row number
611   
612   TArrayI *neighbourhoodRows = 0;
613   TArrayI *neighbourhoodPads = 0;
614   
615   //std::cerr << "Trying to get neighbourhood for row " << row << ", pad " << pad << std::endl;
616   GetNeighbourhood(neighbourhoodRows, neighbourhoodPads, row, pad, rRadius, pRadius);
617   //std::cerr << "Got neighbourhood for row " << row << ", pad " << pad << std::endl;
618   
619   Int_t r, p;
620   for (Int_t i=0; i < (2*rRadius+1)*(2*pRadius+1); i++) {
621     r = neighbourhoodRows->At(i);
622     p = neighbourhoodPads->At(i);
623     if (r == -1 || p == -1) continue;   // window is bigger than ROC
624     tpcROCinstance->GetPositionLocal(fSector, r, p, localXY);   // calculate position localXY by pad and row number
625     dlx = lPad[0] - localXY[0];
626     dly = lPad[1] - localXY[1];
627     //xx[0] = 1;
628     xx[1] = dlx;
629     xx[2] = dly;
630     xx[3] = dlx*dlx;
631     xx[4] = dly*dly;
632     xx[5] = dlx*dly;
633     if (!ROCoutliers || ROCoutliers->GetValue(r,p) != 1) {
634       fitterQ->AddPoint(&xx[1], GetValue(r, p), 1);
635       npoints++;
636     }
637   }
638   
639   delete neighbourhoodRows;
640   delete neighbourhoodPads;
641   
642   if (npoints < 0.5 * ((2*rRadius+1)*(2*pRadius+1)) ) {
643     // std::cerr << "Too few data points for fitting @ row " << row << ", pad " << pad << " in sector " << fSector << std::endl;
644     return 0.;  // for diagnostic
645   }
646   fitterQ->Eval();
647   fitterQ->GetParameters(fitParam);
648   Float_t chi2Q = 0;
649   if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
650   //if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
651   if (robust && chi2Q > chi2Threshold) {
652     //std::cout << "robust fitter called... " << std::endl;
653     fitterQ->EvalRobust(robustFraction);
654     fitterQ->GetParameters(fitParam);
655   }
656   Double_t value = fitParam[0];
657   
658   //if (value < 0) std::cerr << "negative fit-value " << value << " in sector "<< this->fSector << " @ row: " << row << " and pad: " << pad << ", with fitter Chi2 = " << chi2Q <<  std::endl;
659   return value;
660 }
661
662
663
664
665 void AliTPCCalROC::GetNeighbourhood(TArrayI* &rowArray, TArrayI* &padArray, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius) {
666   //
667   // AliTPCCalROC::GetNeighbourhood - PRIVATE
668   // in this function the window for LocalFit is determined
669   //
670   rowArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
671   padArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
672   
673   Int_t rmin = row - rRadius;
674   UInt_t rmax = row + rRadius;
675   
676   // if window goes out of ROC
677   if (rmin < 0) {
678     rmax = rmax - rmin;
679     rmin = 0;
680   }
681   if (rmax >= GetNrows()) {
682     rmin = rmin - (rmax - GetNrows()+1);
683     rmax = GetNrows() - 1;
684       if (rmin  < 0 ) rmin = 0; // if the window is bigger than the ROC
685   }
686   
687   Int_t pmin, pmax;
688   Int_t i = 0;
689   
690   for (UInt_t r = rmin; r <= rmax; r++) {
691     pmin = pad - pRadius;
692     pmax = pad + pRadius;
693     if (pmin < 0) {
694       pmax = pmax - pmin;
695       pmin = 0;
696     }
697     if (pmax >= (Int_t)GetNPads(r)) {
698       pmin = pmin - (pmax - GetNPads(r)+1);
699       pmax = GetNPads(r) - 1;
700       if (pmin  < 0 ) pmin = 0; // if the window is bigger than the ROC
701     }
702     for (Int_t p = pmin; p <= pmax; p++) {
703       (*rowArray)[i] = r;
704       (*padArray)[i] = p;
705       i++;
706     }
707   }
708   for (Int_t j = i; j < rowArray->GetSize(); j++){  // unused padArray-entries, in the case that the window is bigger than the ROC
709     //std::cout << "trying to write -1" << std::endl;
710     (*rowArray)[j] = -1;
711     (*padArray)[j] = -1;
712     //std::cout << "writing -1" << std::endl;
713   }
714 }
715
716
717
718 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){
719   //
720   // Makes a  GlobalFit for the given secotr and return fit-parameters, covariance and chi2
721   // The origin of the fit function is the center of the ROC!
722   // fitType == 0: fit plane function
723   // fitType == 1: fit parabolic function
724   // ROCoutliers - pads with value !=0 are not used in fitting procedure
725   // chi2Threshold: Threshold for chi2 when EvalRobust is called
726   // robustFraction: Fraction of data that will be used in EvalRobust
727   // err: error of the data points
728   //
729   TLinearFitter* fitterG = 0;
730   Double_t xx[6];
731   
732   if (fitType  == 1) {
733     fitterG = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
734     fitParam.ResizeTo(6);
735     covMatrix.ResizeTo(6,6);
736   } else {
737     fitterG = new TLinearFitter(3,"x0++x1++x2");
738     fitParam.ResizeTo(3);
739     covMatrix.ResizeTo(3,3);
740   }
741   fitterG->StoreData(kTRUE);   
742   fitterG->ClearPoints();
743   Int_t    npoints=0;
744   
745   Float_t dlx, dly;
746   Float_t centerPad[3] = {0};
747   Float_t localXY[3] = {0};
748   
749   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
750   tpcROCinstance->GetPositionLocal(fSector, GetNrows()/2, GetNPads(GetNrows()/2)/2, centerPad);  // calculate center of ROC 
751   
752   // loop over all channels and read data into fitterG
753   for (UInt_t irow = 0; irow < GetNrows(); irow++) {
754     for (UInt_t ipad = 0; ipad < GetNPads(irow); ipad++) {
755       // fill fitterG
756       if (ROCoutliers && ROCoutliers->GetValue(irow, ipad) != 0) continue;
757       tpcROCinstance->GetPositionLocal(fSector, irow, ipad, localXY);   // calculate position localXY by pad and row number
758       dlx = localXY[0] - centerPad[0];
759       dly = localXY[1] - centerPad[1];
760       xx[0] = 1;
761       xx[1] = dlx;
762       xx[2] = dly;
763       xx[3] = dlx*dlx;
764       xx[4] = dly*dly;
765       xx[5] = dlx*dly;
766       npoints++;
767       fitterG->AddPoint(xx, GetValue(irow, ipad), err);
768     }
769   }
770   if(npoints>10) { // make sure there is something to fit
771     fitterG->Eval();
772     fitterG->GetParameters(fitParam);
773     fitterG->GetCovarianceMatrix(covMatrix);
774     if (fitType == 1)
775       chi2 = fitterG->GetChisquare()/(npoints-6.);
776     else chi2 = fitterG->GetChisquare()/(npoints-3.);
777     if (robust && chi2 > chi2Threshold) {
778       //    std::cout << "robust fitter called... " << std::endl;
779       fitterG->EvalRobust(robustFraction);
780       fitterG->GetParameters(fitParam);
781     }
782   } else {
783     // set parameteres to 0
784     Int_t nParameters = 3;
785     if (fitType  == 1)
786       nParameters = 6;
787
788     for(Int_t i = 0; i < nParameters; i++)
789       fitParam[i] = 0;
790   }
791   
792   delete fitterG;
793 }
794
795
796 AliTPCCalROC* AliTPCCalROC::CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sector){
797   //
798   // Create ROC with global fit parameters
799   // The origin of the fit function is the center of the ROC!
800   // loop over all channels, write fit values into new ROC and return it
801   //
802   Float_t dlx, dly;
803   Float_t centerPad[3] = {0};
804   Float_t localXY[3] = {0};
805   AliTPCCalROC * xROCfitted = new AliTPCCalROC(sector);
806   AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
807   tpcROCinstance->GetPositionLocal(sector, xROCfitted->GetNrows()/2, xROCfitted->GetNPads(xROCfitted->GetNrows()/2)/2, centerPad);  // calculate center of ROC 
808   Int_t fitType = 1;
809   if (fitParam.GetNoElements() == 6) fitType = 1;
810   else fitType = 0;
811   Double_t value = 0;
812   if (fitType == 1) { // parabolic fit
813     for (UInt_t irow = 0; irow < xROCfitted->GetNrows(); irow++) {
814       for (UInt_t ipad = 0; ipad < xROCfitted->GetNPads(irow); ipad++) {
815         tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY);   // calculate position localXY by pad and row number
816         dlx = localXY[0] - centerPad[0];
817         dly = localXY[1] - centerPad[1];
818         value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
819         xROCfitted->SetValue(irow, ipad, value);
820       }
821     }   
822   }
823   else {  // linear fit
824     for (UInt_t irow = 0; irow < xROCfitted->GetNrows(); irow++) {
825       for (UInt_t ipad = 0; ipad < xROCfitted->GetNPads(irow); ipad++) {
826         tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY);   // calculate position localXY by pad and row number
827         dlx = localXY[0] - centerPad[0];
828         dly = localXY[1] - centerPad[1];
829         value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly;
830         xROCfitted->SetValue(irow, ipad, value);
831       }
832     }   
833   }
834   return xROCfitted;
835 }
836