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