Line data Source code
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 : /* $Id$ */
17 :
18 : ///////////////////////////////////////////////////////////////////////////////
19 : // //
20 : // Calibration base class for a single ROC //
21 : // Contains one UShort_t value per pad //
22 : // However, values are set and get as float, there are stored internally as //
23 : // (UShort_t) value * 10000 //
24 : // //
25 : ///////////////////////////////////////////////////////////////////////////////
26 :
27 : #include <TStyle.h>
28 : #include <TMath.h>
29 : #include <TH2F.h>
30 : #include <TH1F.h>
31 :
32 : #include "AliMathBase.h"
33 :
34 : #include "AliTRDCalROC.h"
35 :
36 48 : ClassImp(AliTRDCalROC)
37 :
38 : //_____________________________________________________________________________
39 : AliTRDCalROC::AliTRDCalROC()
40 7020 : :TObject()
41 7020 : ,fPla(0)
42 7020 : ,fCha(0)
43 7020 : ,fNrows(0)
44 7020 : ,fNcols(0)
45 7020 : ,fNchannels(0)
46 7020 : ,fData(0)
47 35100 : {
48 : //
49 : // Default constructor
50 : //
51 :
52 14040 : }
53 :
54 : //_____________________________________________________________________________
55 : AliTRDCalROC::AliTRDCalROC(Int_t p, Int_t c)
56 0 : :TObject()
57 0 : ,fPla(p)
58 0 : ,fCha(c)
59 0 : ,fNrows(0)
60 0 : ,fNcols(144)
61 0 : ,fNchannels(0)
62 0 : ,fData(0)
63 0 : {
64 : //
65 : // Constructor that initializes a given pad plane type
66 : //
67 :
68 : //
69 : // The pad plane parameter
70 : //
71 0 : switch (p) {
72 : case 0:
73 0 : if (c == 2) {
74 : // L0C0 type
75 0 : fNrows = 12;
76 0 : }
77 : else {
78 : // L0C1 type
79 0 : fNrows = 16;
80 : }
81 : break;
82 : case 1:
83 0 : if (c == 2) {
84 : // L1C0 type
85 0 : fNrows = 12;
86 0 : }
87 : else {
88 : // L1C1 type
89 0 : fNrows = 16;
90 : }
91 : break;
92 : case 2:
93 0 : if (c == 2) {
94 : // L2C0 type
95 0 : fNrows = 12;
96 0 : }
97 : else {
98 : // L2C1 type
99 0 : fNrows = 16;
100 : }
101 : break;
102 : case 3:
103 0 : if (c == 2) {
104 : // L3C0 type
105 0 : fNrows = 12;
106 0 : }
107 : else {
108 : // L3C1 type
109 0 : fNrows = 16;
110 : }
111 : break;
112 : case 4:
113 0 : if (c == 2) {
114 : // L4C0 type
115 0 : fNrows = 12;
116 0 : }
117 : else {
118 : // L4C1 type
119 0 : fNrows = 16;
120 : }
121 : break;
122 : case 5:
123 0 : if (c == 2) {
124 : // L5C0 type
125 0 : fNrows = 12;
126 0 : }
127 : else {
128 : // L5C1 type
129 0 : fNrows = 16;
130 : }
131 : break;
132 : };
133 :
134 0 : fNchannels = fNrows * fNcols;
135 0 : if (fNchannels != 0) {
136 0 : fData = new UShort_t[fNchannels];
137 0 : }
138 :
139 0 : for (Int_t i=0; i<fNchannels; ++i) {
140 0 : fData[i] = 0;
141 : }
142 :
143 0 : }
144 :
145 : //_____________________________________________________________________________
146 : AliTRDCalROC::AliTRDCalROC(const AliTRDCalROC &c)
147 0 : :TObject(c)
148 0 : ,fPla(c.fPla)
149 0 : ,fCha(c.fCha)
150 0 : ,fNrows(c.fNrows)
151 0 : ,fNcols(c.fNcols)
152 0 : ,fNchannels(c.fNchannels)
153 0 : ,fData(0)
154 0 : {
155 : //
156 : // AliTRDCalROC copy constructor
157 : //
158 :
159 : Int_t iBin = 0;
160 :
161 0 : fData = new UShort_t[fNchannels];
162 0 : for (iBin = 0; iBin < fNchannels; iBin++) {
163 0 : fData[iBin] = ((AliTRDCalROC &) c).fData[iBin];
164 : }
165 :
166 0 : }
167 :
168 : //_____________________________________________________________________________
169 : AliTRDCalROC::~AliTRDCalROC()
170 0 : {
171 : //
172 : // AliTRDCalROC destructor
173 : //
174 :
175 0 : if (fData) {
176 0 : delete [] fData;
177 0 : fData = 0;
178 0 : }
179 :
180 0 : }
181 :
182 : //_____________________________________________________________________________
183 : AliTRDCalROC &AliTRDCalROC::operator=(const AliTRDCalROC &c)
184 : {
185 : //
186 : // Assignment operator
187 : //
188 :
189 0 : if (this == &c) {
190 0 : return *this;
191 : }
192 :
193 0 : fPla = c.fPla;
194 0 : fCha = c.fCha;
195 0 : fNrows = c.fNrows;
196 0 : fNcols = c.fNcols;
197 0 : fNchannels = c.fNchannels;
198 :
199 0 : if (fData) {
200 0 : delete [] fData;
201 : }
202 0 : fData = new UShort_t[fNchannels];
203 0 : for (Int_t iBin = 0; iBin < fNchannels; iBin++) {
204 0 : fData[iBin] = ((AliTRDCalROC &) c).fData[iBin];
205 : }
206 :
207 0 : return *this;
208 :
209 0 : }
210 :
211 : //_____________________________________________________________________________
212 : void AliTRDCalROC::Copy(TObject &c) const
213 : {
214 : //
215 : // Copy function
216 : //
217 :
218 0 : ((AliTRDCalROC &) c).fPla = fPla;
219 0 : ((AliTRDCalROC &) c).fCha = fCha;
220 :
221 0 : ((AliTRDCalROC &) c).fNrows = fNrows;
222 0 : ((AliTRDCalROC &) c).fNcols = fNcols;
223 :
224 : Int_t iBin = 0;
225 :
226 0 : ((AliTRDCalROC &) c).fNchannels = fNchannels;
227 :
228 0 : if (((AliTRDCalROC &) c).fData) delete [] ((AliTRDCalROC &) c).fData;
229 0 : ((AliTRDCalROC &) c).fData = new UShort_t[fNchannels];
230 0 : for (iBin = 0; iBin < fNchannels; iBin++) {
231 0 : ((AliTRDCalROC &) c).fData[iBin] = fData[iBin];
232 : }
233 :
234 0 : TObject::Copy(c);
235 :
236 0 : }
237 :
238 : //___________________________________________________________________________________
239 : Double_t AliTRDCalROC::GetMean(AliTRDCalROC* const outlierROC) const
240 : {
241 : //
242 : // Calculate the mean
243 : //
244 :
245 0 : Double_t *ddata = new Double_t[fNchannels];
246 : Int_t nPoints = 0;
247 0 : for (Int_t i=0;i<fNchannels;i++) {
248 0 : if ((!outlierROC) || (!(outlierROC->GetValue(i)))) {
249 : //if(fData[i] > 0.000000000000001){
250 0 : ddata[nPoints]= (Double_t) fData[i]/10000;
251 0 : nPoints++;
252 : //}
253 0 : }
254 : }
255 0 : Double_t mean = TMath::Mean(nPoints,ddata);
256 0 : delete [] ddata;
257 0 : return mean;
258 : }
259 : //___________________________________________________________________________________
260 : Double_t AliTRDCalROC::GetMeanNotNull() const
261 : {
262 : //
263 : // Calculate the mean rejecting null value
264 : //
265 :
266 0 : Double_t *ddata = new Double_t[fNchannels];
267 : Int_t nPoints = 0;
268 0 : for (Int_t i=0;i<fNchannels;i++) {
269 :
270 0 : if(fData[i] > 0.000000000000001){
271 0 : ddata[nPoints]= (Double_t) fData[i]/10000;
272 0 : nPoints++;
273 0 : }
274 :
275 : }
276 0 : if(nPoints<1) return -1.;
277 0 : Double_t mean = TMath::Mean(nPoints,ddata);
278 0 : delete [] ddata;
279 : return mean;
280 0 : }
281 :
282 : //_______________________________________________________________________________________
283 : Double_t AliTRDCalROC::GetMedian(AliTRDCalROC* const outlierROC) const
284 : {
285 : //
286 : // Calculate the median
287 : //
288 :
289 0 : Double_t *ddata = new Double_t[fNchannels];
290 : Int_t nPoints = 0;
291 0 : for (Int_t i=0;i<fNchannels;i++) {
292 0 : if ((!outlierROC) || (!(outlierROC->GetValue(i)))) {
293 0 : if(fData[i] > 0.000000000000001){
294 0 : ddata[nPoints]= (Double_t) fData[i]/10000;
295 0 : nPoints++;
296 0 : }
297 : }
298 : }
299 0 : Double_t mean = TMath::Median(nPoints,ddata);
300 0 : delete [] ddata;
301 0 : return mean;
302 : }
303 :
304 : //____________________________________________________________________________________________
305 : Double_t AliTRDCalROC::GetRMS(AliTRDCalROC* const outlierROC) const
306 : {
307 : //
308 : // Calculate the RMS
309 : //
310 :
311 0 : Double_t *ddata = new Double_t[fNchannels];
312 : Int_t nPoints = 0;
313 0 : for (Int_t i=0;i<fNchannels;i++) {
314 0 : if ((!outlierROC) || (!(outlierROC->GetValue(i)))) {
315 : //if(fData[i] > 0.000000000000001){
316 0 : ddata[nPoints]= (Double_t)fData[i]/10000;
317 0 : nPoints++;
318 : //}
319 0 : }
320 : }
321 0 : Double_t mean = TMath::RMS(nPoints,ddata);
322 0 : delete [] ddata;
323 0 : return mean;
324 : }
325 : //____________________________________________________________________________________________
326 : Double_t AliTRDCalROC::GetRMSNotNull() const
327 : {
328 : //
329 : // Calculate the RMS
330 : //
331 :
332 0 : Double_t *ddata = new Double_t[fNchannels];
333 : Int_t nPoints = 0;
334 0 : for (Int_t i=0;i<fNchannels;i++) {
335 0 : if(fData[i] > 0.000000000000001){
336 0 : ddata[nPoints]= (Double_t)fData[i]/10000;
337 0 : nPoints++;
338 0 : }
339 : }
340 0 : if(nPoints<1) return -1.;
341 0 : Double_t mean = TMath::RMS(nPoints,ddata);
342 0 : delete [] ddata;
343 : return mean;
344 0 : }
345 : //______________________________________________________________________________________________
346 : Double_t AliTRDCalROC::GetLTM(Double_t *sigma, Double_t fraction, AliTRDCalROC* const outlierROC)
347 : {
348 : //
349 : // Calculate LTM mean and sigma
350 : //
351 :
352 0 : Double_t *ddata = new Double_t[fNchannels];
353 0 : Double_t mean=0, lsigma=0;
354 : UInt_t nPoints = 0;
355 0 : for (Int_t i=0;i<fNchannels;i++) {
356 0 : if (!outlierROC || !(outlierROC->GetValue(i))) {
357 0 : if(fData[i] > 0.000000000000001){
358 0 : ddata[nPoints]= (Double_t) fData[i]/10000;
359 0 : nPoints++;
360 0 : }
361 : }
362 : }
363 0 : Int_t hh = TMath::Min(TMath::Nint(fraction *nPoints), Int_t(nPoints));
364 0 : AliMathBase::EvaluateUni(nPoints,ddata, mean, lsigma, hh);
365 0 : if (sigma) *sigma=lsigma;
366 0 : delete [] ddata;
367 0 : return mean;
368 0 : }
369 :
370 : //___________________________________________________________________________________
371 : Bool_t AliTRDCalROC::Add(Float_t c1)
372 : {
373 : //
374 : // add constant
375 : //
376 :
377 : Bool_t result = kTRUE;
378 0 : for (Int_t idata = 0; idata< fNchannels; idata++) {
379 0 : if(((GetValue(idata)+c1) <= 6.5535) && ((GetValue(idata)+c1) >= 0.0)) SetValue(idata,GetValue(idata)+c1);
380 : else {
381 0 : SetValue(idata,0.0);
382 : result = kFALSE;
383 : }
384 : }
385 0 : return result;
386 : }
387 :
388 : //_______________________________________________________________________________________
389 : Bool_t AliTRDCalROC::Multiply(Float_t c1)
390 : {
391 : //
392 : // multiply constant
393 : //
394 :
395 : Bool_t result = kTRUE;
396 0 : if(c1 < 0) return kFALSE;
397 0 : for (Int_t idata = 0; idata< fNchannels; idata++) {
398 0 : if((GetValue(idata)*c1) <= 6.5535) SetValue(idata,GetValue(idata)*c1);
399 : else {
400 0 : SetValue(idata,0.0);
401 : result = kFALSE;
402 : }
403 : }
404 0 : return result;
405 0 : }
406 :
407 : //____________________________________________________________________________________________
408 : Bool_t AliTRDCalROC::Add(const AliTRDCalROC * roc, Double_t c1)
409 : {
410 : //
411 : // add values
412 : //
413 :
414 : Bool_t result = kTRUE;
415 0 : for (Int_t idata = 0; idata< fNchannels; idata++){
416 0 : if(((GetValue(idata)+roc->GetValue(idata)*c1) <= 6.5535) &&
417 0 : ((GetValue(idata)+roc->GetValue(idata)*c1) >= 0.0))
418 0 : SetValue(idata,GetValue(idata)+roc->GetValue(idata)*c1);
419 : else {
420 0 : SetValue(idata,0.0);
421 : result = kFALSE;
422 : }
423 : }
424 0 : return result;
425 : }
426 :
427 : //____________________________________________________________________________________________
428 : Bool_t AliTRDCalROC::Multiply(const AliTRDCalROC* roc)
429 : {
430 : //
431 : // multiply values - per by pad
432 : //
433 :
434 : Bool_t result = kTRUE;
435 0 : for (Int_t idata = 0; idata< fNchannels; idata++){
436 0 : if((GetValue(idata)*roc->GetValue(idata)) <= 6.5535)
437 0 : SetValue(idata,GetValue(idata)*roc->GetValue(idata));
438 : else {
439 0 : SetValue(idata,0.0);
440 : result = kFALSE;
441 : }
442 : }
443 0 : return result;
444 : }
445 :
446 : //______________________________________________________________________________________________
447 : Bool_t AliTRDCalROC::Divide(const AliTRDCalROC* roc)
448 : {
449 : //
450 : // divide values
451 : //
452 :
453 : Bool_t result = kTRUE;
454 : Float_t kEpsilon=0.00000000000000001;
455 0 : for (Int_t idata = 0; idata< fNchannels; idata++){
456 0 : if (TMath::Abs(roc->GetValue(idata))>kEpsilon){
457 0 : if((GetValue(idata)/roc->GetValue(idata)) <= 6.5535)
458 0 : SetValue(idata,GetValue(idata)/roc->GetValue(idata));
459 : else {
460 0 : SetValue(idata,0.0);
461 : result = kFALSE;
462 : }
463 : }
464 : else result = kFALSE;
465 : }
466 0 : return result;
467 : }
468 : //______________________________________________________________________________________________
469 : Bool_t AliTRDCalROC::Unfold()
470 : {
471 : //
472 : // Compute the mean value per pad col
473 : // Divide with this value each pad col
474 : // This is for the noise study
475 : // Return kFALSE if one or more of the pad col was not normalised
476 : //
477 :
478 : Bool_t result = kTRUE;
479 : Float_t kEpsilon=0.00000000000000001;
480 0 : Double_t mmeannotnull = GetMeanNotNull();
481 0 : Double_t rmsnotnull = GetRMSNotNull();
482 : //printf("mmeannotnull %f and rmsnotnull %f\n",mmeannotnull,rmsnotnull);
483 0 : if((mmeannotnull<0.) || (rmsnotnull<0.)) return kFALSE;
484 :
485 :
486 :
487 : // calcul the mean value per col
488 0 : for(Int_t icol = 0; icol < fNcols; icol++){
489 :
490 : Float_t mean = 0.0;
491 : Float_t nb = 0.0;
492 :
493 0 : for(Int_t irow = 0; irow < fNrows; irow++){
494 0 : if(TMath::Abs(GetValue(icol,irow)-mmeannotnull) < 5*rmsnotnull){
495 0 : mean += GetValue(icol,irow);
496 0 : nb += 1.0;
497 0 : }
498 : }
499 :
500 0 : if(nb > kEpsilon) {
501 :
502 0 : mean = mean/nb;
503 :
504 0 : if(mean > kEpsilon){
505 0 : for(Int_t irow = 0; irow < fNrows; irow++){
506 0 : Float_t value = GetValue(icol,irow);
507 0 : SetValue(icol,irow,(Float_t)(value/mean));
508 : }
509 0 : }
510 : else result = kFALSE;
511 :
512 : }
513 :
514 : else result = kFALSE;
515 :
516 :
517 : }
518 :
519 :
520 0 : return result;
521 0 : }
522 : //__________________________________________________________________________________
523 : TH2F * AliTRDCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type, Float_t mu)
524 : {
525 : //
526 : // make 2D histo
527 : // type -1 = user defined range
528 : // 0 = nsigma cut nsigma=min
529 : // 1 = delta cut around median delta=min
530 :
531 : Float_t kEpsilonr = 0.005;
532 0 : gStyle->SetPalette(1);
533 :
534 0 : if (type>=0){
535 0 : if (type==0){
536 : // nsigma range
537 0 : Float_t mean = GetMean();
538 0 : Float_t sigma = GetRMS();
539 0 : Float_t nsigma = TMath::Abs(min);
540 0 : sigma *= nsigma;
541 0 : if(sigma < kEpsilonr) sigma = kEpsilonr;
542 0 : min = mean-sigma;
543 0 : max = mean+sigma;
544 0 : }
545 0 : if (type==1){
546 : // fixed range
547 0 : Float_t mean = GetMedian();
548 : Float_t delta = min;
549 0 : min = mean-delta;
550 0 : max = mean+delta;
551 0 : }
552 0 : if (type==2){
553 0 : Double_t sigma;
554 0 : Float_t mean = GetLTM(&sigma,max);
555 0 : sigma*=min;
556 0 : if(sigma < kEpsilonr) sigma = kEpsilonr;
557 0 : min = mean-sigma;
558 0 : max = mean+sigma;
559 :
560 0 : }
561 : }
562 :
563 0 : char name[1000];
564 0 : snprintf(name,1000,"%s 2D Plane %d Chamber %d",GetTitle(),fPla, fCha);
565 0 : TH2F * his = new TH2F(name,name,fNrows,0, fNrows, fNcols, 0, fNcols);
566 0 : for (Int_t irow=0; irow<fNrows; irow++){
567 0 : for (Int_t icol=0; icol<fNcols; icol++){
568 0 : his->Fill(irow+0.5,icol+0.5,GetValue(icol,irow)*mu);
569 : }
570 : }
571 0 : his->SetStats(0);
572 0 : his->SetMaximum(max);
573 0 : his->SetMinimum(min);
574 0 : return his;
575 0 : }
576 :
577 : //_______________________________________________________________________________________
578 : TH1F * AliTRDCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type, Float_t mu)
579 : {
580 : //
581 : // make 1D histo
582 : // type -1 = user defined range
583 : // 0 = nsigma cut nsigma=min
584 : // 1 = delta cut around median delta=min
585 :
586 : Float_t kEpsilonr = 0.005;
587 :
588 0 : if (type>=0){
589 0 : if (type==0){
590 : // nsigma range
591 0 : Float_t mean = GetMean();
592 0 : Float_t sigma = GetRMS();
593 0 : Float_t nsigma = TMath::Abs(min);
594 0 : sigma *= nsigma;
595 0 : if(sigma < kEpsilonr) sigma = kEpsilonr;
596 0 : min = mean-sigma;
597 0 : max = mean+sigma;
598 0 : }
599 0 : if (type==1){
600 : // fixed range
601 0 : Float_t mean = GetMedian();
602 : Float_t delta = min;
603 0 : min = mean-delta;
604 0 : max = mean+delta;
605 0 : }
606 0 : if (type==2){
607 : //
608 : // LTM mean +- nsigma
609 : //
610 0 : Double_t sigma;
611 0 : Float_t mean = GetLTM(&sigma,max);
612 0 : sigma*=min;
613 0 : if(sigma < kEpsilonr) sigma = kEpsilonr;
614 0 : min = mean-sigma;
615 0 : max = mean+sigma;
616 0 : }
617 : }
618 0 : char name[1000];
619 0 : snprintf(name,1000,"%s 1D Plane %d Chamber %d",GetTitle(),fPla, fCha);
620 0 : TH1F * his = new TH1F(name,name,100, min,max);
621 0 : for (Int_t irow=0; irow<fNrows; irow++){
622 0 : for (Int_t icol=0; icol<fNcols; icol++){
623 0 : his->Fill(GetValue(icol,irow)*mu);
624 : }
625 : }
626 0 : return his;
627 0 : }
|
