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