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