git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix coverity defects
[u/mrichter/AliRoot.git] / TRD / Cal / AliTRDCalDet.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// TRD calibration class for parameters which saved per detector //
21// //
22///////////////////////////////////////////////////////////////////////////////
23
4c87f26e 24#include <TMath.h>
25#include <TH1F.h>
26#include <TH2F.h>
27#include <TStyle.h>
28
7754cd1f 29#include "AliTRDCalDet.h"
4c87f26e 30#include "AliTRDgeometry.h"
31#include "AliMathBase.h"
32#include "AliTRDpadPlane.h"
7754cd1f 33
34ClassImp(AliTRDCalDet)
35
36//_____________________________________________________________________________
37AliTRDCalDet::AliTRDCalDet():TNamed()
38{
39 //
40 // AliTRDCalDet default constructor
41 //
42
43 for (Int_t idet = 0; idet < kNdet; idet++) {
44 fData[idet] = 0;
45 }
46
47}
48
49//_____________________________________________________________________________
50AliTRDCalDet::AliTRDCalDet(const Text_t *name, const Text_t *title)
51 :TNamed(name,title)
52{
53 //
54 // AliTRDCalDet constructor
55 //
56
57 for (Int_t idet = 0; idet < kNdet; idet++) {
58 fData[idet] = 0;
59 }
60
61}
62
7754cd1f 63//_____________________________________________________________________________
64AliTRDCalDet::AliTRDCalDet(const AliTRDCalDet &c):TNamed(c)
65{
66 //
67 // AliTRDCalDet copy constructor
68 //
69
70 ((AliTRDCalDet &) c).Copy(*this);
71
72}
73
74///_____________________________________________________________________________
75AliTRDCalDet::~AliTRDCalDet()
76{
77 //
78 // AliTRDCalDet destructor
79 //
80
81}
82
83//_____________________________________________________________________________
84AliTRDCalDet &AliTRDCalDet::operator=(const AliTRDCalDet &c)
85{
86 //
87 // Assignment operator
88 //
89
90 if (this != &c) ((AliTRDCalDet &) c).Copy(*this);
91 return *this;
92
93}
94
95//_____________________________________________________________________________
96void AliTRDCalDet::Copy(TObject &c) const
97{
98 //
99 // Copy function
100 //
101
102 for (Int_t idet = 0; idet < kNdet; idet++) {
103 ((AliTRDCalDet &) c).fData[idet] = fData[idet];
104 }
105
106 TObject::Copy(c);
107
108}
109
4c87f26e 110//___________________________________________________________________________________
2e32a5ae 111Double_t AliTRDCalDet::GetMean(AliTRDCalDet * const outlierDet) const
4c87f26e 112{
113 //
114 // Calculate the mean
115 //
116
117 if (!outlierDet) return TMath::Mean(kNdet,fData);
118 Double_t *ddata = new Double_t[kNdet];
2e32a5ae 119 Int_t nPoints = 0;
4c87f26e 120 for (Int_t i=0;i<kNdet;i++) {
121 if (!(outlierDet->GetValue(i))) {
2e32a5ae 122 ddata[nPoints]= fData[nPoints];
123 nPoints++;
4c87f26e 124 }
125 }
2e32a5ae 126 Double_t mean = TMath::Mean(nPoints,ddata);
4c87f26e 127 delete [] ddata;
128 return mean;
129}
130
131//_______________________________________________________________________________________
2e32a5ae 132Double_t AliTRDCalDet::GetMedian(AliTRDCalDet * const outlierDet) const
4c87f26e 133{
134 //
135 // Calculate the median
136 //
137
138 if (!outlierDet) return (Double_t) TMath::Median(kNdet,fData);
139 Double_t *ddata = new Double_t[kNdet];
2e32a5ae 140 Int_t nPoints = 0;
4c87f26e 141 for (Int_t i=0;i<kNdet;i++) {
142 if (!(outlierDet->GetValue(i))) {
2e32a5ae 143 ddata[nPoints]= fData[nPoints];
144 nPoints++;
4c87f26e 145 }
146 }
2e32a5ae 147 Double_t mean = TMath::Median(nPoints,ddata);
4c87f26e 148 delete [] ddata;
149 return mean;
150
151}
152
153//____________________________________________________________________________________________
2e32a5ae 154Double_t AliTRDCalDet::GetRMS(AliTRDCalDet * const outlierDet) const
4c87f26e 155{
156 //
157 // Calculate the RMS
158 //
159
160 if (!outlierDet) return TMath::RMS(kNdet,fData);
161 Double_t *ddata = new Double_t[kNdet];
2e32a5ae 162 Int_t nPoints = 0;
4c87f26e 163 for (Int_t i=0;i<kNdet;i++) {
164 if (!(outlierDet->GetValue(i))) {
2e32a5ae 165 ddata[nPoints]= fData[nPoints];
166 nPoints++;
4c87f26e 167 }
168 }
2e32a5ae 169 Double_t mean = TMath::RMS(nPoints,ddata);
4c87f26e 170 delete [] ddata;
171 return mean;
172}
173
a5dcf618 174//____________________________________________________________________________________________
175Double_t AliTRDCalDet::GetRMSRobust(Double_t robust) const
176{
177 //
178 // Calculate the robust RMS
179 //
180
181 // sorted
5c85bc96 182 Int_t *index = new Int_t[kNdet+1];
a5dcf618 183 TMath::Sort((Int_t)kNdet,fData,index);
184
185 // reject
186 Double_t reject = (Int_t) (kNdet*(1-robust)/2.0);
187 if(reject <= 0.0) reject = 0.0;
188 if(reject >= kNdet) reject = 0.0;
189 //printf("Rejecter %f\n",reject);
190
191 Double_t *ddata = new Double_t[kNdet];
192 Int_t nPoints = 0;
193 for (Int_t i=0;i<kNdet;i++) {
194 Bool_t rej = kFALSE;
195 for(Int_t k = 0; k < reject; k++){
196 if(i==index[k]) rej = kTRUE;
197 if(i==index[kNdet-(k+1)]) rej = kTRUE;
198 }
199 if(!rej){
200 ddata[nPoints]= fData[i];
201 nPoints++;
202 }
203 }
204 //printf("Number of points %d\n",nPoints);
205 Double_t mean = TMath::RMS(nPoints,ddata);
206 delete [] ddata;
207 delete [] index;
208 return mean;
209}
210
4c87f26e 211//______________________________________________________________________________________________
2e32a5ae 212Double_t AliTRDCalDet::GetLTM(Double_t *sigma
213 , Double_t fraction
214 , AliTRDCalDet * const outlierDet)
4c87f26e 215{
216 //
217 // Calculate LTM mean and sigma
218 //
219
220 Double_t *ddata = new Double_t[kNdet];
221 Double_t mean=0, lsigma=0;
2e32a5ae 222 UInt_t nPoints = 0;
4c87f26e 223 for (Int_t i=0;i<kNdet;i++) {
224 if (!outlierDet || !(outlierDet->GetValue(i))) {
2e32a5ae 225 ddata[nPoints]= fData[nPoints];
226 nPoints++;
4c87f26e 227 }
228 }
2e32a5ae 229 Int_t hh = TMath::Min(TMath::Nint(fraction *nPoints), Int_t(nPoints));
230 AliMathBase::EvaluateUni(nPoints,ddata, mean, lsigma, hh);
4c87f26e 231 if (sigma) *sigma=lsigma;
232 delete [] ddata;
233 return mean;
234}
235
236//_________________________________________________________________________
237TH1F * AliTRDCalDet::MakeHisto1Distribution(Float_t min, Float_t max,Int_t type)
238{
239 //
240 // make 1D histo
241 // type -1 = user defined range
242 // 0 = nsigma cut nsigma=min
243 // 1 = delta cut around median delta=min
244 //
245
246 if (type>=0){
247 if (type==0){
248 // nsigma range
249 Float_t mean = GetMean();
250 Float_t sigma = GetRMS();
251 Float_t nsigma = TMath::Abs(min);
252 min = mean-nsigma*sigma;
253 max = mean+nsigma*sigma;
254 }
255 if (type==1){
256 // fixed range
257 Float_t mean = GetMedian();
258 Float_t delta = min;
259 min = mean-delta;
260 max = mean+delta;
261 }
262 if (type==2){
263 //
264 // LTM mean +- nsigma
265 //
266 Double_t sigma;
267 Float_t mean = GetLTM(&sigma,max);
268 sigma*=min;
269 min = mean-sigma;
270 max = mean+sigma;
271 }
272 }
273 char name[1000];
02f3bfcc 274 snprintf(name,1000,"%s CalDet 1Distribution",GetTitle());
4c87f26e 275 TH1F * his = new TH1F(name,name,100, min,max);
276 for (Int_t idet=0; idet<kNdet; idet++){
277 his->Fill(GetValue(idet));
278 }
279 return his;
280}
281
282//________________________________________________________________________________
283TH1F * AliTRDCalDet::MakeHisto1DAsFunctionOfDet(Float_t min, Float_t max,Int_t type)
284{
285 //
286 // make 1D histo
287 // type -1 = user defined range
288 // 0 = nsigma cut nsigma=min
289 // 1 = delta cut around median delta=min
290 //
291
292 if (type>=0){
293 if (type==0){
294 // nsigma range
295 Float_t mean = GetMean();
296 Float_t sigma = GetRMS();
297 Float_t nsigma = TMath::Abs(min);
298 min = mean-nsigma*sigma;
299 max = mean+nsigma*sigma;
300 }
301 if (type==1){
302 // fixed range
303 Float_t mean = GetMedian();
304 Float_t delta = min;
305 min = mean-delta;
306 max = mean+delta;
307 }
308 if (type==2){
309 Double_t sigma;
310 Float_t mean = GetLTM(&sigma,max);
311 sigma*=min;
312 min = mean-sigma;
313 max = mean+sigma;
314
315 }
316 }
317
318 char name[1000];
a987273c 319 snprintf(name,1000,"%s CalDet as function of det",GetTitle());
4c87f26e 320 TH1F * his = new TH1F(name,name,kNdet, 0, kNdet);
321 for(Int_t det = 0; det< kNdet; det++){
322 his->Fill(det+0.5,GetValue(det));
323 }
324
325 his->SetMaximum(max);
326 his->SetMinimum(min);
327 return his;
328}
329
330//_____________________________________________________________________________
331TH2F *AliTRDCalDet::MakeHisto2DCh(Int_t ch, Float_t min, Float_t max, Int_t type)
332{
333 //
334 // Make 2D graph
335 // ch - chamber
336 // type -1 = user defined range
337 // 0 = nsigma cut nsigma=min
338 // 1 = delta cut around median delta=min
339 //
340
341 gStyle->SetPalette(1);
342 if (type>=0){
343 if (type==0){
344 // nsigma range
345 Float_t mean = GetMean();
346 Float_t sigma = GetRMS();
347 Float_t nsigma = TMath::Abs(min);
348 min = mean-nsigma*sigma;
349 max = mean+nsigma*sigma;
350 }
351 if (type==1){
352 // fixed range
353 Float_t mean = GetMedian();
354 Float_t delta = min;
355 min = mean-delta;
356 max = mean+delta;
357 }
358 if (type==2){
359 Double_t sigma;
360 Float_t mean = GetLTM(&sigma,max);
361 sigma*=min;
362 min = mean-sigma;
363 max = mean+sigma;
364
365 }
366 }
367
015cd5ba 368 AliTRDgeometry *trdGeo = new AliTRDgeometry();
4c87f26e 369
370 Double_t poslocal[3] = {0.0,0.0,0.0};
371 Double_t posglobal[3] = {0.0,0.0,0.0};
372
373 char name[1000];
a987273c 374 snprintf(name,1000,"%s CalDet 2D ch %d",GetTitle(),ch);
4c87f26e 375 TH2F * his = new TH2F(name, name, 400,-400.0,400.0,400,-400.0,400.0);
376
377 // Where we begin
378 Int_t offsetch = 6*ch;
379
380
381 for (Int_t isec = 0; isec < kNsect; isec++){
382 for(Int_t ipl = 0; ipl < kNplan; ipl++){
383 Int_t det = offsetch+isec*30+ipl;
015cd5ba 384 AliTRDpadPlane *padPlane = trdGeo->GetPadPlane(ipl,ch);
4c87f26e 385 for (Int_t icol=0; icol<padPlane->GetNcols(); icol++){
015cd5ba 386 poslocal[0] = trdGeo->GetTime0(ipl);
4c87f26e 387 poslocal[2] = padPlane->GetRowPos(0);
388 poslocal[1] = padPlane->GetColPos(icol);
015cd5ba 389 trdGeo->RotateBack(det,poslocal,posglobal);
4c87f26e 390 Int_t binx = 1+TMath::Nint((posglobal[0]+400.0)*0.5);
391 Int_t biny = 1+TMath::Nint((posglobal[1]+400.0)*0.5);
392 his->SetBinContent(binx,biny,fData[det]);
393 }
394 }
395 }
396 his->SetXTitle("x (cm)");
397 his->SetYTitle("y (cm)");
398 his->SetStats(0);
399 his->SetMaximum(max);
400 his->SetMinimum(min);
015cd5ba 401 delete trdGeo;
4c87f26e 402 return his;
403}
404
405//_____________________________________________________________________________
406TH2F *AliTRDCalDet::MakeHisto2DSmPl(Int_t sm, Int_t pl, Float_t min, Float_t max, Int_t type)
407{
408 //
409 // Make 2D graph
410 // sm - supermodule number
411 // pl - plane number
412 // type -1 = user defined range
413 // 0 = nsigma cut nsigma=min
414 // 1 = delta cut around median delta=min
415 //
416
417 gStyle->SetPalette(1);
418 if (type>=0){
419 if (type==0){
420 // nsigma range
421 Float_t mean = GetMean();
422 Float_t sigma = GetRMS();
423 Float_t nsigma = TMath::Abs(min);
424 min = mean-nsigma*sigma;
425 max = mean+nsigma*sigma;
426 }
427 if (type==1){
428 // fixed range
429 Float_t mean = GetMedian();
430 Float_t delta = min;
431 min = mean-delta;
432 max = mean+delta;
433 }
434 if (type==2){
435 Double_t sigma;
436 Float_t mean = GetLTM(&sigma,max);
437 sigma*=min;
438 min = mean-sigma;
439 max = mean+sigma;
440
441 }
442 }
443
015cd5ba 444 AliTRDgeometry *trdGeo = new AliTRDgeometry();
445 AliTRDpadPlane *padPlane0 = trdGeo->GetPadPlane(pl,0);
4c87f26e 446 Double_t row0 = padPlane0->GetRow0();
447 Double_t col0 = padPlane0->GetCol0();
448
4c87f26e 449 char name[1000];
a987273c 450 snprintf(name,1000,"%s CalDet 2D sm %d and pl %d",GetTitle(),sm,pl);
4c87f26e 451 TH2F * his = new TH2F( name, name, 5, -TMath::Abs(row0), TMath::Abs(row0)
452 , 4,-2*TMath::Abs(col0),2*TMath::Abs(col0));
453
454 // Where we begin
455 Int_t offsetsmpl = 30*sm+pl;
456
457
458 for (Int_t k = 0; k < kNcham; k++){
459 Int_t det = offsetsmpl+k*6;
460 Int_t kb = kNcham-1-k;
461 his->SetBinContent(kb+1,2,fData[det]);
462 his->SetBinContent(kb+1,3,fData[det]);
463 }
464 his->SetXTitle("z (cm)");
465 his->SetYTitle("y (cm)");
466 his->SetStats(0);
467 his->SetMaximum(max);
468 his->SetMinimum(min);
015cd5ba 469 delete trdGeo;
4c87f26e 470 return his;
471}
472
473//_____________________________________________________________________________
474void AliTRDCalDet::Add(Float_t c1)
475{
476 //
477 // Add constant for all detectors
478 //
479
480 for (Int_t idet = 0; idet < kNdet; idet++) {
481 fData[idet] += c1;
482 }
483}
484
485//_____________________________________________________________________________
486void AliTRDCalDet::Multiply(Float_t c1)
487{
488 //
489 // multiply constant for all detectors
490 //
491 for (Int_t idet = 0; idet < kNdet; idet++) {
492 fData[idet] *= c1;
493 }
494}
495
496//_____________________________________________________________________________
497void AliTRDCalDet::Add(const AliTRDCalDet * calDet, Double_t c1)
498{
499 //
500 // add caldet channel by channel multiplied by c1
501 //
502 for (Int_t idet = 0; idet < kNdet; idet++) {
503 fData[idet] += calDet->GetValue(idet)*c1;
504 }
505}
506
507//_____________________________________________________________________________
508void AliTRDCalDet::Multiply(const AliTRDCalDet * calDet)
509{
510 //
511 // multiply caldet channel by channel
512 //
513 for (Int_t idet = 0; idet < kNdet; idet++) {
514 fData[idet] *= calDet->GetValue(idet);
515 }
516}
517
518//_____________________________________________________________________________
519void AliTRDCalDet::Divide(const AliTRDCalDet * calDet)
520{
521 //
522 // divide caldet channel by channel
523 //
524 Float_t kEpsilon=0.00000000000000001;
525 for (Int_t idet = 0; idet < kNdet; idet++) {
526 if (TMath::Abs(calDet->GetValue(idet))>kEpsilon){
527 fData[idet] /= calDet->GetValue(idet);
528 }
529 }
530}
adaea377 531//_____________________________________________________________________________
532Double_t AliTRDCalDet::CalcMean(Bool_t wghtPads)
533{
534 // Calculate the mean value after rejection of the chambers not calibrated
535 // wghPads = kTRUE weighted with the number of pads in case of a AliTRDCalPad created (t0)
536 //
537
538 Int_t iSM;
539 Double_t sum = 0.0;
540 Int_t ndet = 0;
541 Double_t meanALL = 0.0;
542 Double_t meanWP = 0.0;
543 Double_t pads=0.0;
544 Double_t padsALL=(144*16*24+144*12*6)*18;
5c85bc96 545 Double_t *meanSM = new Double_t[18];
546 Double_t *meanSMWP = new Double_t[18];
adaea377 547
548 Int_t det = 0;
549 while(det < 540) {
550 Float_t val= fData[det];
551 iSM = (Int_t)(det / (6*5));
552 pads=(((Int_t) (det % (6 * 5)) / 6) == 2) ? 144*12 : 144*16;
553 meanALL+=val/540.;
554 meanSM[iSM]+=val/30.;
555 meanWP+=val*(pads/padsALL);
556 meanSMWP[iSM]+=val*(pads/(padsALL/18.));
557
558 /*
559 printf(" det %d val %.3f meanALL %.5f meanWP %.5f meanSM[%d] %.5f meanSMWP[%d] %.5f \n",
560 det,
561 val,
562 meanALL,
563 meanWP,
564 iSM,
565 meanSM[iSM],
566 iSM,
567 meanSMWP[iSM]);
568 */
569
570 det++;
571 }
a5dcf618 572
adaea377 573 // debug
574 /*
575 printf(" ALL %.5f \n",meanALL);
576 printf(" WP %.5f \n",meanWP);
577 for(Int_t i=0; i<18; i++) printf(" SM %02d %.5f \n",i,meanSM[i]);
578 for(Int_t i=0; i<18; i++) printf(" SM %02d %.5f \n",i,meanSMWP[i]);
579 */
580
581 det=0;
582 while(det < 540) {
583 Float_t val= fData[det];
584 if (( (!wghtPads) &&
585 (TMath::Abs(val - meanALL) > 0.0001) &&
586 (TMath::Abs(val - meanSM[(Int_t)(det / (6*5))]) > 0.0001) ) ||
587 ( (wghtPads) &&
588 (TMath::Abs(val - meanWP) > 0.0001) &&
589 (TMath::Abs(val - meanSMWP[(Int_t)(det / (6*5) )]) > 0.0001) )
590 ) {
591 sum+=val;
592 ndet++;
593 }
594 det++;
595 }
596
5c85bc96 597 delete []meanSM;
598 delete []meanSMWP;
599
adaea377 600 return (sum/ndet);
601}