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