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