AliRoot Core  edcc906 (edcc906)
AliTPCCalROC.cxx
Go to the documentation of this file.
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 
23 
24 // ROOT includes
25 #include "TMath.h"
26 #include "TClass.h"
27 #include "TFile.h"
28 #include "TH1F.h"
29 #include "TH2F.h"
30 #include "TArrayI.h"
31 
32 //
33 //
34 #include "AliTPCCalROC.h"
35 #include "AliMathBase.h"
36 
37 #include "TRandom3.h" // only needed by the AliTPCCalROCTest() method
38 
40 ClassImp(AliTPCCalROC)
42 
43 
44 //_____________________________________________________________________________
46  :TNamed(),
47  fSector(0),
48  fNChannels(0),
49  fNRows(0),
50  fkIndexes(0),
51  fData(0)
52 {
53  //
54  // Default constructor
55  //
56 
57 }
58 
59 //_____________________________________________________________________________
61  :TNamed(),
62  fSector(0),
63  fNChannels(0),
64  fNRows(0),
65  fkIndexes(0),
66  fData(0)
67 {
69 
70  fSector = sector;
74  fData = new Float_t[fNChannels];
75  Reset();
76 
77 }
78 
79 //_____________________________________________________________________________
81  :TNamed(c),
82  fSector(0),
83  fNChannels(0),
84  fNRows(0),
85  fkIndexes(0),
86  fData(0)
87 {
89 
90  fSector = c.fSector;
94  //
95  fData = new Float_t[fNChannels];
96  for (UInt_t idata = 0; idata< fNChannels; idata++) fData[idata] = c.fData[idata];
97 }
98 //____________________________________________________________________________
100 {
102 
103  if (this == &param) return (*this);
104  fSector = param.fSector;
108  //
109  if (fData) delete [] fData;
110  fData = new Float_t[fNChannels];
111  for (UInt_t idata = 0; idata< fNChannels; idata++) fData[idata] = param.fData[idata];
112  return (*this);
113 }
114 
115 
116 //_____________________________________________________________________________
118 {
120 
121  if (fData) {
122  delete [] fData;
123  fData = 0;
124  }
125 }
126 
127 
128 
129 void AliTPCCalROC::Streamer(TBuffer &R__b)
130 {
132 
133  if (R__b.IsReading()) {
134  AliTPCCalROC::Class()->ReadBuffer(R__b, this);
136  } else {
137  AliTPCCalROC::Class()->WriteBuffer(R__b,this);
138  }
139 }
140 
141 void AliTPCCalROC::Print(Option_t* /*option*/) const{
142  //
143  // Print summary content
144  //
145  printf("|ROC%d|\t%f|\t%f|\t%f|\n",fSector, GetMean(),GetMedian(),GetRMS());
146 }
147 
148 
149 
150 //
151 
152 Bool_t AliTPCCalROC::MedianFilter(Int_t deltaRow, Int_t deltaPad, AliTPCCalROC* outlierROC, Bool_t doEdge){
154 
155  Float_t *newBuffer=new Float_t[fNChannels] ;
156  Double_t *cacheBuffer=new Double_t[fNChannels];
157  //
158  for (Int_t iRow=0; iRow< Int_t(fNRows); iRow++){
159  Int_t nPads=GetNPads(iRow); // number of rows in current row
160  for (Int_t iPad=0; iPad<nPads; iPad++){
161  Double_t value=GetValue(iRow,iPad);
162  Int_t counter=0;
163  //
164  for (Int_t dRow=-deltaRow; dRow<=deltaRow; dRow++){
165  Int_t jRow=iRow+dRow; //take the row - mirror if ouside of range
166  Float_t sign0=1.;
167  if (jRow<0) sign0=-1.;
168  if (UInt_t(jRow)>=fNRows) sign0=-1.;
169  Int_t jPads= GetNPads(iRow+sign0*dRow);
170  Int_t offset=(nPads-jPads)/2.;
171  //
172  for (Int_t dPad=-deltaPad; dPad<=deltaPad; dPad++){
173  Float_t sign=sign0;
174  Int_t jPad=iPad+dPad+offset; //take the pad - mirror if ouside of range
175  Int_t kRow=jRow;
176  if (jPad<0) sign=-1;
177  if (jPad>=jPads) sign=-1;
178  if (sign<0){
179  kRow=iRow-dRow;
180  jPad=iPad-dPad+offset;
181  if (!doEdge) continue;
182  }
183  if (IsInRange(UInt_t(kRow),jPad)){
184  Bool_t isOutlier=(outlierROC==NULL)?kFALSE:outlierROC->GetValue(kRow,jPad)>0;
185  if (!isOutlier){
186  cacheBuffer[counter]=sign*(GetValue(kRow,jPad)-value);
187  counter++;
188  }
189  }
190  }
191  }
192  newBuffer[fkIndexes[iRow]+iPad]=0.;
193  if (counter>1) newBuffer[fkIndexes[iRow]+iPad] = TMath::Median(counter, cacheBuffer)+value;
194  }
195  }
196  memcpy(fData, newBuffer,GetNchannels()*sizeof(Float_t));
197  delete []newBuffer;
198  delete []cacheBuffer;
199  return kTRUE;
200 }
201 
202 
203 
204 Bool_t AliTPCCalROC::LTMFilter(Int_t deltaRow, Int_t deltaPad, Float_t fraction, Int_t type, AliTPCCalROC* outlierROC, Bool_t doEdge){
207 
208  if (fraction<0 || fraction>1) return kFALSE;
209  Float_t *newBuffer=new Float_t[fNChannels] ;
210  Double_t *cacheBuffer=new Double_t[fNChannels];
211  //
212  for (Int_t iRow=0; iRow< Int_t(fNRows); iRow++){
213  Int_t nPads=GetNPads(iRow); // number of rows in current row
214  for (Int_t iPad=0; iPad<nPads; iPad++){
215  Double_t value=GetValue(iRow,iPad);
216  Int_t counter=0;
217  //
218  for (Int_t dRow=-deltaRow; dRow<=deltaRow; dRow++){
219  Int_t jRow=iRow+dRow; //take the row - mirror if ouside of range
220  Float_t sign0=1.;
221  if (jRow<0) sign0=-1.;
222  if (UInt_t(jRow)>=fNRows) sign0=-1.;
223  Int_t jPads= GetNPads(iRow+sign0*dRow);
224  Int_t offset=(nPads-jPads)/2.;
225  //
226  for (Int_t dPad=-deltaPad; dPad<=deltaPad; dPad++){
227  Float_t sign=sign0;
228  Int_t jPad=iPad+dPad+offset; //take the pad - mirror if ouside of range
229  Int_t kRow=jRow;
230  if (jPad<0) sign=-1;
231  if (jPad>=jPads) sign=-1;
232  if (sign<0){
233  if (!doEdge) continue;
234  kRow=iRow-dRow;
235  jPad=iPad-dPad+offset;
236  }
237  if (IsInRange(UInt_t(kRow),jPad)){
238  Bool_t isOutlier=(outlierROC==NULL)?kFALSE:outlierROC->GetValue(kRow,jPad)>0;
239  if (!isOutlier){
240  cacheBuffer[counter]=sign*(GetValue(kRow,jPad)-value);
241  counter++;
242  }
243  }
244  }
245  }
246  Double_t mean=0,rms=0;
247  if (TMath::Nint(fraction*Double_t(counter))>1 ) AliMathBase::EvaluateUni(counter,cacheBuffer,mean,rms,TMath::Nint(fraction*Double_t(counter)));
248  mean+=value;
249  newBuffer[fkIndexes[iRow]+iPad] = (type==0)? mean:rms;
250  }
251  }
252  memcpy(fData, newBuffer,GetNchannels()*sizeof(Float_t));
253  delete []newBuffer;
254  delete []cacheBuffer;
255  return kTRUE;
256 }
257 
258 Bool_t AliTPCCalROC::Convolute(Double_t sigmaPad, Double_t sigmaRow, AliTPCCalROC*outlierROC, TF1 */*fpad*/, TF1 */*frow*/){
261 
262  Float_t *newBuffer=new Float_t[fNChannels] ;
263  //
264  for (Int_t iRow=0; iRow< Int_t(fNRows); iRow++){
265  Int_t nPads=GetNPads(iRow); // number of rows in current row
266  for (Int_t iPad=0; iPad<nPads; iPad++){
267  Int_t jRow0=TMath::Max(TMath::Nint(iRow-sigmaRow*4.),0);
268  Int_t jRow1=TMath::Min(TMath::Nint(iRow+sigmaRow*4.),Int_t(fNRows));
269  Int_t jPad0=TMath::Max(TMath::Nint(iPad-sigmaPad*4.),0);
270  Int_t jPad1=TMath::Min(TMath::Nint(iRow+sigmaPad*4.),Int_t(nPads));
271  //
272  Double_t sumW=0;
273  Double_t sumCal=0;
274  for (Int_t jRow=jRow0; jRow<=jRow1; jRow++){
275  for (Int_t jPad=jPad0; jPad<=jPad1; jPad++){
276  if (!IsInRange(jRow,jPad)) continue;
277  Bool_t isOutlier=(outlierROC==NULL)?kFALSE:outlierROC->GetValue(jRow,jPad)>0;
278  if (!isOutlier){
279  Double_t weight= TMath::Gaus(jPad,iPad,sigmaPad)*TMath::Gaus(jRow,iRow,sigmaRow);
280  sumCal+=weight*GetValue(jRow,jPad);
281  sumW+=weight;
282  }
283  }
284  }
285  if (sumW>0){
286  sumCal/=sumW;
287  newBuffer[fkIndexes[iRow]+iPad]=sumCal;
288  }
289  }
290  }
291  memcpy(fData, newBuffer,GetNchannels()*sizeof(Float_t));
292  delete []newBuffer;
293  return kTRUE;
294 }
295 
296 
297 //
298 
299 
300 
301 
302 
303 // algebra fuctions:
304 
305 void AliTPCCalROC::Add(Float_t c1){
307 
308  for (UInt_t idata = 0; idata< fNChannels; idata++) fData[idata]+=c1;
309 }
310 
311 
312 void AliTPCCalROC::Multiply(Float_t c1){
314 
315  for (UInt_t idata = 0; idata< fNChannels; idata++) fData[idata]*=c1;
316 }
317 
318 
319 void AliTPCCalROC::Add(const AliTPCCalROC * roc, Double_t c1){
322 
323  if (!roc) return;
324  for (UInt_t idata = 0; idata< fNChannels; idata++){
325  fData[idata]+=roc->fData[idata]*c1;
326  }
327 }
328 
329 
333 
334  if (!roc) return;
335  for (UInt_t idata = 0; idata< fNChannels; idata++){
336  fData[idata]*=roc->fData[idata];
337  }
338 }
339 
340 
344 
345  if (!roc) return;
346  Float_t kEpsilon=0.00000000000000001;
347  for (UInt_t idata = 0; idata< fNChannels; idata++){
348  if (TMath::Abs(roc->fData[idata])>kEpsilon)
349  fData[idata]/=roc->fData[idata];
350  }
351 }
352 
354 {
356 
357  memset(fData,0,sizeof(Float_t)*fNChannels); // set all values to 0
358 }
359 
360 //____________________________________________________________
361 Double_t AliTPCCalROC::GetStats(EStatType statType, AliTPCCalROC *const outlierROC, EPadType padType) const
362 {
367 
368  if(fSector<36) padType=kAll;
369  if (!outlierROC) {
370  if(padType==kAll) {
371  if (statType==kMean) return TMath::Mean (fNChannels, fData);
372  else if(statType==kMedian) return TMath::Median (fNChannels, fData);
373  else if(statType==kRMS) return TMath::RMS (fNChannels, fData);
374  else if(statType==kMinElement) return TMath::MinElement(fNChannels, fData);
375  else if(statType==kMaxElement) return TMath::MaxElement(fNChannels, fData);
376  }
377  else if(padType==kOROCmedium) {
378  if (statType==kMean) return TMath::Mean (fkIndexes[64], fData);
379  else if(statType==kMedian) return TMath::Median (fkIndexes[64], fData);
380  else if(statType==kRMS) return TMath::RMS (fkIndexes[64], fData);
381  else if(statType==kMinElement) return TMath::MinElement(fkIndexes[64], fData);
382  else if(statType==kMaxElement) return TMath::MaxElement(fkIndexes[64], fData);
383  }
384  else if(padType==kOROClong) {
385  const Int_t offset=fkIndexes[64];
386  if (statType==kMean) return TMath::Mean (fNChannels-offset, fData+offset);
387  else if(statType==kMedian) return TMath::Median (fNChannels-offset, fData+offset);
388  else if(statType==kRMS) return TMath::RMS (fNChannels-offset, fData+offset);
389  else if(statType==kMinElement) return TMath::MinElement(fNChannels-offset, fData+offset);
390  else if(statType==kMaxElement) return TMath::MaxElement(fNChannels-offset, fData+offset);
391  }
392  }//end of no outlierROC
393 
394  Float_t ddata[fNChannels];
395 
396  Int_t nPoints = 0;
397  UInt_t indexMin = 0;
398  UInt_t indexMax = fNChannels;
399  if(padType==kOROCmedium) indexMax=fkIndexes[64];
400  else if(padType==kOROClong) indexMin=fkIndexes[64];
401 
402  for (UInt_t i=indexMin;i<indexMax;i++) {
403  if (outlierROC->GetValue(i)>1e-20) continue;
404  ddata[nPoints]= fData[i];
405  nPoints++;
406  }
407  Double_t value = 0;
408  if(nPoints>0){
409  if (statType==kMean) value = TMath::Mean (nPoints, ddata);
410  else if(statType==kMedian) value = TMath::Median (nPoints, ddata);
411  else if(statType==kRMS) value = TMath::RMS (nPoints, ddata);
412  else if(statType==kMinElement) value = TMath::MinElement(nPoints, ddata);
413  else if(statType==kMaxElement) value = TMath::MaxElement(nPoints, ddata);
414  }
415  return value;
416 }
417 
418 //_____________________________________________________________________________
419 Double_t AliTPCCalROC::GetMean(AliTPCCalROC *const outlierROC,EPadType padType) const
420 {
424 
425  return GetStats(kMean,outlierROC,padType);
426 }
427 
428 //_____________________________________________________________________________
429 Double_t AliTPCCalROC::GetMedian(AliTPCCalROC *const outlierROC,EPadType padType) const
430 {
434 
435  return GetStats(kMedian,outlierROC,padType);
436 }
437 
438 //_____________________________________________________________________________
439 Double_t AliTPCCalROC::GetRMS(AliTPCCalROC *const outlierROC,EPadType padType) const
440 {
444 
445  return GetStats(kRMS,outlierROC,padType);
446 }
447 
448 //_____________________________________________________________________________
449 Double_t AliTPCCalROC::GetMinElement(AliTPCCalROC *const outlierROC,EPadType padType) const
450 {
454 
455  return GetStats(kMinElement,outlierROC,padType);
456 }
457 
458 //_____________________________________________________________________________
459 Double_t AliTPCCalROC::GetMaxElement(AliTPCCalROC *const outlierROC,EPadType padType) const
460 {
464 
465  return GetStats(kMaxElement,outlierROC,padType);
466 }
467 
468 //_____________________________________________________________________________
469 Double_t AliTPCCalROC::GetLTM(Double_t *const sigma, Double_t fraction, AliTPCCalROC *const outlierROC, EPadType padType){
474 
475  if(fSector<36) padType=kAll;
476  Double_t ddata[fNChannels];
477  UInt_t nPoints = 0;
478  UInt_t indexMin = 0;
479  UInt_t indexMax = fNChannels;
480  if(padType==kOROCmedium) indexMax=fkIndexes[64];
481  else if(padType==kOROClong) indexMin=fkIndexes[64];
482  for (UInt_t i=indexMin;i<indexMax;i++) {
483  if (outlierROC && (outlierROC->GetValue(i) >1e-20)) continue;
484  ddata[nPoints]= fData[i];
485  nPoints++;
486  }
487 
488  Double_t ltm =0, lsigma=0;
489  if(nPoints>0) {
490  Int_t hh = TMath::Min(TMath::Nint(fraction *nPoints), Int_t(nPoints));
491  AliMathBase::EvaluateUni(nPoints,ddata, ltm, lsigma, hh);
492  if (sigma) *sigma=lsigma;
493  }
494 
495  return ltm;
496 }
497 
498 //_____________________________________________________________________________
499 TH1F * AliTPCCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type){
504 
505  if (type>=0){
506  if (type==0){
507  // nsigma range
508  Float_t mean = GetMean();
509  Float_t sigma = GetRMS();
510  Float_t nsigma = TMath::Abs(min);
511  min = mean-nsigma*sigma;
512  max = mean+nsigma*sigma;
513  }
514  if (type==1){
515  // fixed range
516  Float_t mean = GetMedian();
517  Float_t delta = min;
518  min = mean-delta;
519  max = mean+delta;
520  }
521  if (type==2){
522  //
523  // LTM mean +- nsigma
524  //
525  Double_t sigma;
526  Float_t mean = GetLTM(&sigma,max);
527  sigma*=min;
528  min = mean-sigma;
529  max = mean+sigma;
530  }
531  }
532  TString name=Form("%s ROC 1D%d",GetTitle(),fSector);
533  TH1F * his = new TH1F(name.Data(),name.Data(),100, min,max);
534  for (UInt_t irow=0; irow<fNRows; irow++){
535  UInt_t npads = (Int_t)GetNPads(irow);
536  for (UInt_t ipad=0; ipad<npads; ipad++){
537  his->Fill(GetValue(irow,ipad));
538  }
539  }
540  return his;
541 }
542 
543 
544 
545 TH2F * AliTPCCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type){
550 
551  if (type>=0){
552  if (type==0){
553  // nsigma range
554  Float_t mean = GetMean();
555  Float_t sigma = GetRMS();
556  Float_t nsigma = TMath::Abs(min);
557  min = mean-nsigma*sigma;
558  max = mean+nsigma*sigma;
559  }
560  if (type==1){
561  // fixed range
562  Float_t mean = GetMedian();
563  Float_t delta = min;
564  min = mean-delta;
565  max = mean+delta;
566  }
567  if (type==2){
568  Double_t sigma;
569  Float_t mean = GetLTM(&sigma,max);
570  sigma*=min;
571  min = mean-sigma;
572  max = mean+sigma;
573 
574  }
575  }
576  UInt_t maxPad = 0;
577  for (UInt_t irow=0; irow<fNRows; irow++){
578  if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
579  }
580 
581  TString name=Form("%s ROC%d",GetTitle(),fSector);
582  TH2F * his = new TH2F(name.Data(),name.Data(),fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
583  for (UInt_t irow=0; irow<fNRows; irow++){
584  UInt_t npads = (Int_t)GetNPads(irow);
585  for (UInt_t ipad=0; ipad<npads; ipad++){
586  his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,GetValue(irow,ipad));
587  }
588  }
589  his->SetMaximum(max);
590  his->SetMinimum(min);
591  return his;
592 }
593 
594 TH2F * AliTPCCalROC::MakeHistoOutliers(Float_t delta, Float_t fraction, Int_t type){
601 
602  Double_t sigma;
603  Float_t mean = GetLTM(&sigma,fraction);
604  if (type==0) delta*=sigma;
605  UInt_t maxPad = 0;
606  for (UInt_t irow=0; irow<fNRows; irow++){
607  if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
608  }
609 
610  TString name=Form("%s ROC Outliers%d",GetTitle(),fSector);
611  TH2F * his = new TH2F(name.Data(),name.Data(),fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
612  for (UInt_t irow=0; irow<fNRows; irow++){
613  UInt_t npads = (Int_t)GetNPads(irow);
614  for (UInt_t ipad=0; ipad<npads; ipad++){
615  if (TMath::Abs(GetValue(irow,ipad)-mean)>delta)
616  his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,1);
617  }
618  }
619  return his;
620 }
621 
622 
623 
624 void AliTPCCalROC::Draw(Option_t* opt){
626 
627  TH1 * his=0;
628  TString option=opt;
629  option.ToUpper();
630  if (option.Contains("1D")){
631  his = MakeHisto1D();
632  }
633  else{
634  his = MakeHisto2D();
635  }
636  his->Draw(option);
637 }
638 
639 
640 
641 
642 
645 
646  Float_t kEpsilon=0.00001;
647 
648  // create CalROC with defined values
649  AliTPCCalROC roc0(0);
650  for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
651  for (UInt_t ipad = 0; ipad <roc0.GetNPads(irow); ipad++){
652  Float_t value = irow+ipad/1000.;
653  roc0.SetValue(irow,ipad,value);
654  }
655  }
656 
657  // copy CalROC, readout values and compare to original
658  AliTPCCalROC roc1(roc0);
659  for (UInt_t irow = 0; irow <roc1.GetNrows(); irow++){
660  for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
661  Float_t value = irow+ipad/1000.;
662  if (roc1.GetValue(irow,ipad)!=value){
663  printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
664  }
665  }
666  }
667 
668  // write original CalROC to file
669  TFile f("calcTest.root","recreate");
670  roc0.Write("Roc0");
671  AliTPCCalROC * roc2 = (AliTPCCalROC*)f.Get("Roc0");
672  f.Close();
673 
674  // read CalROC from file and compare to original CalROC
675  for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
676  if (roc0.GetNPads(irow)!=roc2->GetNPads(irow))
677  printf("NPads - Read/Write error\trow=%d\n",irow);
678  for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
679  Float_t value = irow+ipad/1000.;
680  if (roc2->GetValue(irow,ipad)!=value){
681  printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
682  }
683  }
684  }
685 
686  //
687  // Algebra Tests
688  //
689 
690  // Add constant
691  AliTPCCalROC roc3(roc0);
692  roc3.Add(1.5);
693  for (UInt_t irow = 0; irow <roc3.GetNrows(); irow++){
694  for (UInt_t ipad = 0; ipad <roc3.GetNPads(irow); ipad++){
695  Float_t value = irow+ipad/1000. + 1.5;
696  if (TMath::Abs(roc3.GetValue(irow,ipad)-value) > kEpsilon){
697  printf("Add constant - error\trow=%d\tpad=%d\n",irow,ipad);
698  }
699  }
700  }
701 
702  // Add another CalROC
703  AliTPCCalROC roc4(roc0);
704  roc4.Add(&roc0, -1.5);
705  for (UInt_t irow = 0; irow <roc4.GetNrows(); irow++){
706  for (UInt_t ipad = 0; ipad <roc4.GetNPads(irow); ipad++){
707  Float_t value = irow+ipad/1000. - 1.5 * (irow+ipad/1000.);
708  if (TMath::Abs(roc4.GetValue(irow,ipad)-value) > kEpsilon){
709  printf("Add CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
710  }
711  }
712  }
713 
714  // Multiply with constant
715  AliTPCCalROC roc5(roc0);
716  roc5.Multiply(-1.4);
717  for (UInt_t irow = 0; irow <roc5.GetNrows(); irow++){
718  for (UInt_t ipad = 0; ipad <roc5.GetNPads(irow); ipad++){
719  Float_t value = (irow+ipad/1000.) * (-1.4);
720  if (TMath::Abs(roc5.GetValue(irow,ipad)-value) > kEpsilon){
721  printf("Multiply with constant - error\trow=%d\tpad=%d\n",irow,ipad);
722  }
723  }
724  }
725 
726  // Multiply another CalROC
727  AliTPCCalROC roc6(roc0);
728  roc6.Multiply(&roc0);
729  for (UInt_t irow = 0; irow <roc6.GetNrows(); irow++){
730  for (UInt_t ipad = 0; ipad <roc6.GetNPads(irow); ipad++){
731  Float_t value = (irow+ipad/1000.) * (irow+ipad/1000.);
732  if (TMath::Abs(roc6.GetValue(irow,ipad)-value) > kEpsilon*100){
733  printf("Multiply with CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
734  }
735  }
736  }
737 
738 
739  // Divide by CalROC
740  AliTPCCalROC roc7(roc0);
741  roc7.Divide(&roc0);
742  for (UInt_t irow = 0; irow <roc7.GetNrows(); irow++){
743  for (UInt_t ipad = 0; ipad <roc7.GetNPads(irow); ipad++){
744  Float_t value = 1.;
745  if (irow+ipad == 0) value = roc0.GetValue(irow,ipad);
746  if (TMath::Abs(roc7.GetValue(irow,ipad)-value) > kEpsilon){
747  printf("Multiply with CalROC - error\trow=%d\tpad=%d\n",irow,ipad);
748  }
749  }
750  }
751 
752  //
753  // Statistics Test
754  //
755 
756  // create CalROC with defined values
757  TRandom3 rnd(0);
758  AliTPCCalROC sroc0(0);
759  for (UInt_t ichannel = 0; ichannel < sroc0.GetNchannels(); ichannel++){
760  Float_t value = rnd.Gaus(10., 2.);
761  sroc0.SetValue(ichannel,value);
762  }
763 
764  printf("Mean (should be close to 10): %f\n", sroc0.GetMean());
765  printf("RMS (should be close to 2): %f\n", sroc0.GetRMS());
766  printf("Median (should be close to 10): %f\n", sroc0.GetMedian());
767  printf("LTM (should be close to 10): %f\n", sroc0.GetLTM());
768 
769  //AliTPCCalROC* sroc1 = sroc0.LocalFit(4, 8);
770 
771  //delete sroc1;
772 
773 // std::cout << TMath::Abs(roc5.GetValue(irow,ipad)-value) << std::endl;
774 }
775 
776 
777 AliTPCCalROC * AliTPCCalROC::LocalFit(Int_t rowRadius, Int_t padRadius, AliTPCCalROC* ROCoutliers, Bool_t robust, Double_t chi2Threshold, Double_t robustFraction) {
790 
791  AliTPCCalROC * xROCfitted = new AliTPCCalROC(fSector);
792  TLinearFitter fitterQ(6,"hyp5");
793  // TLinearFitter fitterQ(6,"x0++x1++x2++x3++x4++x5");
794  fitterQ.StoreData(kTRUE);
795  for (UInt_t row=0; row < GetNrows(); row++) {
796  //std::cout << "Entering row " << row << " of " << GetNrows() << " @ sector "<< fSector << " for local fitting... "<< std::endl;
797  for (UInt_t pad=0; pad < GetNPads(row); pad++)
798  xROCfitted->SetValue(row, pad, GetNeighbourhoodValue(&fitterQ, row, pad, rowRadius, padRadius, ROCoutliers, robust, chi2Threshold, robustFraction));
799  }
800  return xROCfitted;
801 }
802 
803 
804 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) {
807 
808  fitterQ->ClearPoints();
809  TVectorD fitParam(6);
810  Int_t npoints=0;
811  Double_t xx[6];
812  Float_t dlx, dly;
813  Float_t lPad[3] = {0};
814  Float_t localXY[3] = {0};
815 
816  AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
817  tpcROCinstance->GetPositionLocal(fSector, row, pad, lPad); // calculate position lPad by pad and row number
818 
819  TArrayI *neighbourhoodRows = 0;
820  TArrayI *neighbourhoodPads = 0;
821 
822  //std::cerr << "Trying to get neighbourhood for row " << row << ", pad " << pad << std::endl;
823  GetNeighbourhood(neighbourhoodRows, neighbourhoodPads, row, pad, rRadius, pRadius);
824  //std::cerr << "Got neighbourhood for row " << row << ", pad " << pad << std::endl;
825 
826  Int_t r, p;
827  for (Int_t i=0; i < (2*rRadius+1)*(2*pRadius+1); i++) {
828  r = neighbourhoodRows->At(i);
829  p = neighbourhoodPads->At(i);
830  if (r == -1 || p == -1) continue; // window is bigger than ROC
831  tpcROCinstance->GetPositionLocal(fSector, r, p, localXY); // calculate position localXY by pad and row number
832  dlx = lPad[0] - localXY[0];
833  dly = lPad[1] - localXY[1];
834  //xx[0] = 1;
835  xx[1] = dlx;
836  xx[2] = dly;
837  xx[3] = dlx*dlx;
838  xx[4] = dly*dly;
839  xx[5] = dlx*dly;
840  if (!ROCoutliers || ROCoutliers->GetValue(r,p) != 1) {
841  fitterQ->AddPoint(&xx[1], GetValue(r, p), 1);
842  npoints++;
843  }
844  }
845 
846  delete neighbourhoodRows;
847  delete neighbourhoodPads;
848 
849  if (npoints < 0.5 * ((2*rRadius+1)*(2*pRadius+1)) ) {
850  // std::cerr << "Too few data points for fitting @ row " << row << ", pad " << pad << " in sector " << fSector << std::endl;
851  return 0.; // for diagnostic
852  }
853  fitterQ->Eval();
854  fitterQ->GetParameters(fitParam);
855  Float_t chi2Q = 0;
856  if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
857  //if (robust) chi2Q = fitterQ->GetChisquare()/(npoints-6.);
858  if (robust && chi2Q > chi2Threshold) {
859  //std::cout << "robust fitter called... " << std::endl;
860  fitterQ->EvalRobust(robustFraction);
861  fitterQ->GetParameters(fitParam);
862  }
863  Double_t value = fitParam[0];
864 
865  //if (value < 0) std::cerr << "negative fit-value " << value << " in sector "<< this->fSector << " @ row: " << row << " and pad: " << pad << ", with fitter Chi2 = " << chi2Q << std::endl;
866  return value;
867 }
868 
869 
870 
871 
872 void AliTPCCalROC::GetNeighbourhood(TArrayI* &rowArray, TArrayI* &padArray, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius) {
875 
876  rowArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
877  padArray = new TArrayI((2*rRadius+1)*(2*pRadius+1));
878 
879  Int_t rmin = row - rRadius;
880  UInt_t rmax = row + rRadius;
881 
882  // if window goes out of ROC
883  if (rmin < 0) {
884  rmax = rmax - rmin;
885  rmin = 0;
886  }
887  if (rmax >= GetNrows()) {
888  rmin = rmin - (rmax - GetNrows()+1);
889  rmax = GetNrows() - 1;
890  if (rmin < 0 ) rmin = 0; // if the window is bigger than the ROC
891  }
892 
893  Int_t pmin, pmax;
894  Int_t i = 0;
895 
896  for (UInt_t r = rmin; r <= rmax; r++) {
897  pmin = pad - pRadius;
898  pmax = pad + pRadius;
899  if (pmin < 0) {
900  pmax = pmax - pmin;
901  pmin = 0;
902  }
903  if (pmax >= (Int_t)GetNPads(r)) {
904  pmin = pmin - (pmax - GetNPads(r)+1);
905  pmax = GetNPads(r) - 1;
906  if (pmin < 0 ) pmin = 0; // if the window is bigger than the ROC
907  }
908  for (Int_t p = pmin; p <= pmax; p++) {
909  (*rowArray)[i] = r;
910  (*padArray)[i] = p;
911  i++;
912  }
913  }
914  for (Int_t j = i; j < rowArray->GetSize(); j++){ // unused padArray-entries, in the case that the window is bigger than the ROC
915  //std::cout << "trying to write -1" << std::endl;
916  (*rowArray)[j] = -1;
917  (*padArray)[j] = -1;
918  //std::cout << "writing -1" << std::endl;
919  }
920 }
921 
922 
923 
924 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, EPadType padType){
934 
935  if(fSector<36) padType=kAll;
936 
937  TLinearFitter* fitterG = 0;
938  Double_t xx[6];
939 
940  if (fitType == 1) {
941  fitterG = new TLinearFitter (6,"x0++x1++x2++x3++x4++x5");
942  fitParam.ResizeTo(6);
943  covMatrix.ResizeTo(6,6);
944  } else {
945  fitterG = new TLinearFitter(3,"x0++x1++x2");
946  fitParam.ResizeTo(3);
947  covMatrix.ResizeTo(3,3);
948  }
949  fitterG->StoreData(kTRUE);
950  fitterG->ClearPoints();
951  Int_t npoints=0;
952 
953  Float_t dlx, dly;
954  Float_t centerPad[3] = {0};
955  Float_t localXY[3] = {0};
956 
957  AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
958  if(padType==kAll)
959  tpcROCinstance->GetPositionLocal(fSector, GetNrows()/2, GetNPads(GetNrows()/2)/2, centerPad); // calculate center of ROC
960  else if(padType==kOROCmedium)
961  tpcROCinstance->GetPositionLocal(fSector, 64/2, GetNPads(64/2)/2, centerPad); // calculate center of ROC medium pads
962  else if(padType==kOROClong)
963  tpcROCinstance->GetPositionLocal(fSector, (GetNrows()-64)/2+64, GetNPads((GetNrows()-64)/2+64)/2, centerPad); // calculate center of ROC long pads
964 
965  UInt_t irowMin = 0;
966  UInt_t irowMax = GetNrows();
967  if(padType==kOROCmedium) irowMax = 64;
968  else if(padType==kOROClong) irowMin = 64;
969 
970  // loop over all channels and read data into fitterG
971  for (UInt_t irow = irowMin; irow < irowMax; irow++) {
972  for (UInt_t ipad = 0; ipad < GetNPads(irow); ipad++) {
973  // fill fitterG
974  if (ROCoutliers && ROCoutliers->GetValue(irow, ipad) != 0) continue;
975  tpcROCinstance->GetPositionLocal(fSector, irow, ipad, localXY); // calculate position localXY by pad and row number
976  dlx = localXY[0] - centerPad[0];
977  dly = localXY[1] - centerPad[1];
978  xx[0] = 1;
979  xx[1] = dlx;
980  xx[2] = dly;
981  xx[3] = dlx*dlx;
982  xx[4] = dly*dly;
983  xx[5] = dlx*dly;
984  npoints++;
985  fitterG->AddPoint(xx, GetValue(irow, ipad), err);
986  }
987  }
988  if(npoints>10) { // make sure there is something to fit
989  fitterG->Eval();
990  fitterG->GetParameters(fitParam);
991  fitterG->GetCovarianceMatrix(covMatrix);
992  if (fitType == 1)
993  chi2 = fitterG->GetChisquare()/(npoints-6.);
994  else chi2 = fitterG->GetChisquare()/(npoints-3.);
995  if (robust && chi2 > chi2Threshold) {
996  // std::cout << "robust fitter called... " << std::endl;
997  fitterG->EvalRobust(robustFraction);
998  fitterG->GetParameters(fitParam);
999  }
1000  } else {
1001  // set parameteres to 0
1002  Int_t nParameters = 3;
1003  if (fitType == 1)
1004  nParameters = 6;
1005 
1006  for(Int_t i = 0; i < nParameters; i++)
1007  fitParam[i] = 0;
1008  }
1009 
1010  delete fitterG;
1011 }
1012 
1013 AliTPCCalROC* AliTPCCalROC::CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sector, EPadType padType, AliTPCCalROC *oldTPCCalROC ){
1020 
1021  if(sector<36) padType=kAll;
1022 
1023  Float_t dlx, dly;
1024  Float_t centerPad[3] = {0};
1025  Float_t localXY[3] = {0};
1026  AliTPCCalROC * xROCfitted = NULL;
1027  if(oldTPCCalROC!=0 && padType!=kAll) xROCfitted = oldTPCCalROC;
1028  else xROCfitted = new AliTPCCalROC(sector);
1029  AliTPCROC* tpcROCinstance = AliTPCROC::Instance();
1030  if(padType==kAll)
1031  tpcROCinstance->GetPositionLocal(sector, xROCfitted->GetNrows()/2, xROCfitted->GetNPads(xROCfitted->GetNrows()/2)/2, centerPad); // calculate center of ROC
1032  else if(padType==kOROCmedium)
1033  tpcROCinstance->GetPositionLocal(sector, 64/2, xROCfitted->GetNPads(64/2)/2, centerPad); // calculate center of ROC medium pads
1034  else if(padType==kOROClong)
1035  tpcROCinstance->GetPositionLocal(sector, (xROCfitted->GetNrows()-64)/2+64, xROCfitted->GetNPads((xROCfitted->GetNrows()-64)/2+64)/2, centerPad); // calculate center of ROC long pads
1036 
1037  UInt_t irowMin = 0;
1038  UInt_t irowMax = xROCfitted->GetNrows();
1039  if(padType==kOROCmedium) irowMax = 64;
1040  else if(padType==kOROClong) irowMin = 64;
1041 
1042  Int_t fitType = 1;
1043  if (fitParam.GetNoElements() == 6) fitType = 1;
1044  else fitType = 0;
1045  Double_t value = 0;
1046  if (fitType == 1) { // parabolic fit
1047  for (UInt_t irow = irowMin; irow < irowMax; irow++) {
1048  for (UInt_t ipad = 0; ipad < xROCfitted->GetNPads(irow); ipad++) {
1049  tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY); // calculate position localXY by pad and row number
1050  dlx = localXY[0] - centerPad[0];
1051  dly = localXY[1] - centerPad[1];
1052  value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly + fitParam[3]*dlx*dlx + fitParam[4]*dly*dly + fitParam[5]*dlx*dly;
1053  xROCfitted->SetValue(irow, ipad, value);
1054  }
1055  }
1056  }
1057  else { // linear fit
1058  for (UInt_t irow = irowMin; irow < irowMax; irow++) {
1059  for (UInt_t ipad = 0; ipad < xROCfitted->GetNPads(irow); ipad++) {
1060  tpcROCinstance->GetPositionLocal(sector, irow, ipad, localXY); // calculate position localXY by pad and row number
1061  dlx = localXY[0] - centerPad[0];
1062  dly = localXY[1] - centerPad[1];
1063  value = fitParam[0] + fitParam[1]*dlx + fitParam[2]*dly;
1064  xROCfitted->SetValue(irow, ipad, value);
1065  }
1066  }
1067  }
1068  return xROCfitted;
1069 }
1070 
UInt_t fSector
sector number
Definition: AliTPCCalROC.h:80
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
const Double_t kEpsilon
TH1F * MakeHisto1D(Float_t min=4, Float_t max=-4, Int_t type=0)
void Add(Float_t c1)
UInt_t fNChannels
number of channels
Definition: AliTPCCalROC.h:81
Double_t GetMinElement(AliTPCCalROC *const outlierROC=0, EPadType padType=kAll) const
Double_t GetMedian(AliTPCCalROC *const outlierROC=0, EPadType padType=kAll) const
Double_t GetRMS(AliTPCCalROC *const outlierROC=0, EPadType padType=kAll) const
virtual ~AliTPCCalROC()
UInt_t GetNRows(UInt_t sector) const
Definition: AliTPCROC.h:28
TH2F * MakeHistoOutliers(Float_t delta=4, Float_t fraction=0.7, Int_t mode=0)
Float_t GetValue(UInt_t row, UInt_t pad) const
Definition: AliTPCCalROC.h:38
void Multiply(Float_t c1)
AliTPCCalROC * LocalFit(Int_t rowRadius, Int_t padRadius, AliTPCCalROC *ROCoutliers=0, Bool_t robust=kFALSE, Double_t chi2Threshold=5, Double_t robustFraction=0.7)
Float_t p[]
Definition: kNNTest.C:133
UInt_t GetNChannels(UInt_t sector) const
Definition: AliTPCROC.h:29
void GetPositionLocal(UInt_t sector, UInt_t row, UInt_t pad, Float_t *pos)
Definition: AliTPCROC.cxx:414
void GlobalFit(const AliTPCCalROC *ROCoutliers, Bool_t robust, TVectorD &fitParam, TMatrixD &covMatrix, Float_t &chi2, Int_t fitType=1, Double_t chi2Threshold=5, Double_t robustFraction=0.7, Double_t err=1, EPadType padType=kAll)
static void Test()
npoints
Definition: driftITSTPC.C:85
void Divide(const AliTPCCalROC *roc)
TH2F * MakeHisto2D(Float_t min=4, Float_t max=-4, Int_t type=0)
UInt_t GetNPads(UInt_t row) const
Definition: AliTPCCalROC.h:37
UInt_t fNRows
number of rows
Definition: AliTPCCalROC.h:82
const UInt_t * fkIndexes
Definition: AliTPCCalROC.h:83
Double_t chi2
Definition: AnalyzeLaser.C:7
Double_t 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)
Bool_t IsInRange(UInt_t row, UInt_t pad)
Definition: AliTPCCalROC.h:36
Double_t GetMean(AliTPCCalROC *const outlierROC=0, EPadType padType=kAll) const
const UInt_t * GetRowIndexes(UInt_t sector) const
Definition: AliTPCROC.h:33
Geometry class for a single ROC.
Definition: AliTPCROC.h:14
void SetValue(UInt_t row, UInt_t pad, Float_t vd)
Definition: AliTPCCalROC.h:40
virtual void Print(Option_t *option="") const
TPC calibration base class for one ROC.
Definition: AliTPCCalROC.h:20
Float_t * fData
Data.
Definition: AliTPCCalROC.h:85
AliTPCCalROC & operator=(const AliTPCCalROC &param)
TF1 * f
Definition: interpolTest.C:21
UInt_t GetNchannels() const
Definition: AliTPCCalROC.h:34
void GetNeighbourhood(TArrayI *&rowArray, TArrayI *&padArray, Int_t row, Int_t pad, Int_t rRadius, Int_t pRadius)
Double_t GetStats(EStatType statType, AliTPCCalROC *const outlierROC=0, EPadType padType=kAll) const
Bool_t MedianFilter(Int_t deltaRow, Int_t deltaPad, AliTPCCalROC *outlierROC=0, Bool_t doEdge=kTRUE)
virtual void Draw(Option_t *option="")
Bool_t Convolute(Double_t sigmaPad, Double_t sigmaRow, AliTPCCalROC *outlierPad=0, TF1 *fpad=0, TF1 *frow=0)
static AliTPCROC * Instance()
Definition: AliTPCROC.cxx:34
class TVectorT< Double_t > TVectorD
static AliTPCCalROC * CreateGlobalFitCalROC(TVectorD &fitParam, Int_t sector, EPadType padType=kAll, AliTPCCalROC *oldTPCCalROC=0)
Double_t GetLTM(Double_t *const sigma=0, Double_t fraction=0.9, AliTPCCalROC *const outlierROC=0, EPadType padType=kAll)
Bool_t LTMFilter(Int_t deltaRow, Int_t deltaPad, Float_t fraction, Int_t type, AliTPCCalROC *outlierROC=0, Bool_t doEdge=kTRUE)
Double_t GetMaxElement(AliTPCCalROC *const outlierROC=0, EPadType padType=kAll) const
static void EvaluateUni(const Int_t nvectors, const Double_t *data, Double_t &mean, Double_t &sigma, const Int_t hSub)
Definition: AliMathBase.cxx:62
class TMatrixT< Double_t > TMatrixD
UInt_t GetNrows() const
Definition: AliTPCCalROC.h:33