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