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