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