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