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