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