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