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