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