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