a/ AliTRDCalibTask.cxx .h: one histo more to quantify the event selection if any...
[u/mrichter/AliRoot.git] / TRD / Cal / AliTRDCalDet.cxx
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
24 #include <TMath.h>
25 #include <TH1F.h>
26 #include <TH2F.h>
27 #include <TStyle.h>
28
29 #include "AliTRDCalDet.h"
30 #include "AliTRDgeometry.h"
31 #include "AliMathBase.h"
32 #include "AliTRDpadPlane.h"
33
34 ClassImp(AliTRDCalDet)
35
36 //_____________________________________________________________________________
37 AliTRDCalDet::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 //_____________________________________________________________________________
50 AliTRDCalDet::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
63 //_____________________________________________________________________________
64 AliTRDCalDet::AliTRDCalDet(const AliTRDCalDet &c):TNamed(c)
65 {
66   //
67   // AliTRDCalDet copy constructor
68   //
69
70   ((AliTRDCalDet &) c).Copy(*this);
71
72 }
73
74 ///_____________________________________________________________________________
75 AliTRDCalDet::~AliTRDCalDet()
76 {
77   //
78   // AliTRDCalDet destructor
79   //
80
81 }
82
83 //_____________________________________________________________________________
84 AliTRDCalDet &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 //_____________________________________________________________________________
96 void 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
110 //___________________________________________________________________________________
111 Double_t AliTRDCalDet::GetMean(AliTRDCalDet * const outlierDet) const
112 {
113   //
114   // Calculate the mean
115   //
116
117   if (!outlierDet) return TMath::Mean(kNdet,fData);
118   Double_t *ddata = new Double_t[kNdet];
119   Int_t nPoints = 0;
120   for (Int_t i=0;i<kNdet;i++) {
121     if (!(outlierDet->GetValue(i))) {
122       ddata[nPoints]= fData[nPoints];
123       nPoints++;
124     }
125   }
126   Double_t mean = TMath::Mean(nPoints,ddata);
127   delete [] ddata;
128   return mean;
129 }
130
131 //_______________________________________________________________________________________
132 Double_t AliTRDCalDet::GetMedian(AliTRDCalDet * const outlierDet) const
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];
140   Int_t nPoints = 0;
141   for (Int_t i=0;i<kNdet;i++) {
142     if (!(outlierDet->GetValue(i))) {
143       ddata[nPoints]= fData[nPoints];
144       nPoints++;
145     }
146   }
147   Double_t mean = TMath::Median(nPoints,ddata);
148   delete [] ddata;
149   return mean;
150
151 }
152
153 //____________________________________________________________________________________________
154 Double_t AliTRDCalDet::GetRMS(AliTRDCalDet * const outlierDet) const
155 {
156   //
157   // Calculate the RMS
158   //
159
160   if (!outlierDet) return TMath::RMS(kNdet,fData);
161   Double_t *ddata = new Double_t[kNdet];
162   Int_t nPoints = 0;
163   for (Int_t i=0;i<kNdet;i++) {
164     if (!(outlierDet->GetValue(i))) {
165       ddata[nPoints]= fData[nPoints];
166       nPoints++;
167     }
168   }
169   Double_t mean = TMath::RMS(nPoints,ddata);
170   delete [] ddata;
171   return mean;
172 }
173
174 //____________________________________________________________________________________________
175 Double_t AliTRDCalDet::GetRMSRobust(Double_t robust) const
176 {
177   //
178   // Calculate the robust RMS
179   //
180
181   // sorted
182   Int_t *index = new Int_t[kNdet];
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
211 //______________________________________________________________________________________________
212 Double_t AliTRDCalDet::GetLTM(Double_t *sigma
213                             , Double_t fraction
214                             , AliTRDCalDet * const outlierDet)
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;
222   UInt_t nPoints = 0;
223   for (Int_t i=0;i<kNdet;i++) {
224      if (!outlierDet || !(outlierDet->GetValue(i))) {
225         ddata[nPoints]= fData[nPoints];
226         nPoints++;
227      }
228   }
229   Int_t hh = TMath::Min(TMath::Nint(fraction *nPoints), Int_t(nPoints));
230   AliMathBase::EvaluateUni(nPoints,ddata, mean, lsigma, hh);
231   if (sigma) *sigma=lsigma;
232   delete [] ddata;
233   return mean;
234 }
235
236 //_________________________________________________________________________
237 TH1F * 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];
274   snprintf(name,1000,"%s CalDet 1Distribution",GetTitle());
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 //________________________________________________________________________________
283 TH1F * 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];
319   snprintf(name,1000,"%s CalDet as function of det",GetTitle());
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 //_____________________________________________________________________________
331 TH2F *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     
368   AliTRDgeometry *trdGeo = new AliTRDgeometry();
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];
374   snprintf(name,1000,"%s CalDet 2D ch %d",GetTitle(),ch);
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;
384       AliTRDpadPlane *padPlane = trdGeo->GetPadPlane(ipl,ch);
385       for (Int_t icol=0; icol<padPlane->GetNcols(); icol++){
386         poslocal[0] = trdGeo->GetTime0(ipl);
387         poslocal[2] = padPlane->GetRowPos(0);
388         poslocal[1] = padPlane->GetColPos(icol);
389         trdGeo->RotateBack(det,poslocal,posglobal);
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);
401   delete trdGeo;
402   return his;
403 }
404
405 //_____________________________________________________________________________
406 TH2F *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      
444   AliTRDgeometry *trdGeo = new AliTRDgeometry();
445   AliTRDpadPlane *padPlane0 = trdGeo->GetPadPlane(pl,0);
446   Double_t row0    = padPlane0->GetRow0();
447   Double_t col0    = padPlane0->GetCol0();
448
449   char  name[1000];
450   snprintf(name,1000,"%s CalDet 2D sm %d and pl %d",GetTitle(),sm,pl);
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);
469   delete trdGeo;
470   return his;
471 }
472
473 //_____________________________________________________________________________
474 void 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 //_____________________________________________________________________________
486 void 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 //_____________________________________________________________________________
497 void 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 //_____________________________________________________________________________
508 void 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 //_____________________________________________________________________________
519 void 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 }
531
532