Algebraic function added (Marian)
[u/mrichter/AliRoot.git] / TPC / AliTPCCalROC.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
17 ///////////////////////////////////////////////////////////////////////////////
18 //                                                                           //
19 //  Calibration base class for a single ROC                                  //
20 //  Contains one float value per pad                                         //
21 //     mapping of the pads taken form AliTPCROC                              //
22 //                                                                           //
23 ///////////////////////////////////////////////////////////////////////////////
24
25 #include "AliTPCCalROC.h"
26 #include "TMath.h"
27 #include "TClass.h"
28 #include "TFile.h"
29 #include "TH1F.h"
30 #include "TH2F.h"
31 #include "AliMathBase.h"
32 ClassImp(AliTPCCalROC)
33
34
35 //_____________________________________________________________________________
36 AliTPCCalROC::AliTPCCalROC()
37              :TObject(),
38               fSector(0),
39               fNChannels(0),
40               fNRows(0),
41               fIndexes(0),
42               fData(0)
43 {
44   //
45   // Default constructor
46   //
47
48 }
49
50 //_____________________________________________________________________________
51 AliTPCCalROC::AliTPCCalROC(UInt_t  sector)
52              :TObject(),
53               fSector(0),
54               fNChannels(0),
55               fNRows(0),
56               fIndexes(0),
57               fData(0)
58 {
59   //
60   // Constructor that initializes a given sector
61   //
62   fSector = sector;
63   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
64   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
65   fIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
66   fData = new Float_t[fNChannels];
67   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = 0.;
68 }
69
70 //_____________________________________________________________________________
71 AliTPCCalROC::AliTPCCalROC(const AliTPCCalROC &c)
72              :TObject(c),
73               fSector(0),
74               fNChannels(0),
75               fNRows(0),
76               fIndexes(0),
77               fData(0)
78 {
79   //
80   // AliTPCCalROC copy constructor
81   //
82   fSector = c.fSector;
83   fNChannels    =  AliTPCROC::Instance()->GetNChannels(fSector);
84   fNRows        =  AliTPCROC::Instance()->GetNRows(fSector);
85   fIndexes      =  AliTPCROC::Instance()->GetRowIndexes(fSector);
86   //
87   fData   = new Float_t[fNChannels];
88   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata] = c.fData[idata];
89 }
90 //____________________________________________________________________________
91 AliTPCCalROC & AliTPCCalROC::operator =(const AliTPCCalROC & param)
92 {
93   //
94   // assignment operator - dummy
95   //
96   fData=param.fData;
97   return (*this);
98 }
99
100
101 //_____________________________________________________________________________
102 AliTPCCalROC::~AliTPCCalROC()
103 {
104   //
105   // AliTPCCalROC destructor
106   //
107   if (fData) {
108     delete [] fData;
109     fData = 0;
110   }
111 }
112
113
114
115 void AliTPCCalROC::Streamer(TBuffer &R__b)
116 {
117    // Stream an object of class AliTPCCalROC.
118    if (R__b.IsReading()) {
119       AliTPCCalROC::Class()->ReadBuffer(R__b, this);
120       fIndexes =  AliTPCROC::Instance()->GetRowIndexes(fSector);
121    } else {
122       AliTPCCalROC::Class()->WriteBuffer(R__b,this);
123    }
124 }
125
126 //  //
127 //   // algebra
128 //   void Add(Float_t c1);
129 //   void Multiply(Float_t c1);
130 //   void Add(const AliTPCCalROC * roc, Double_t c1 = 1);
131 //   void Divide(const AliTPCCalROC * roc);   
132
133 void AliTPCCalROC::Add(Float_t c1){
134   //
135   // add constant
136   //
137   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata]+=c1;
138 }
139 void AliTPCCalROC::Multiply(Float_t c1){
140   //
141   // add constant
142   //
143   for (UInt_t  idata = 0; idata< fNChannels; idata++) fData[idata]*=c1;
144 }
145
146 void AliTPCCalROC::Add(const AliTPCCalROC * roc, Double_t c1){
147   //
148   // add values 
149   //
150   for (UInt_t  idata = 0; idata< fNChannels; idata++){
151     fData[idata]+=roc->fData[idata]*c1;
152   }
153 }
154
155
156 void AliTPCCalROC::Multiply(const AliTPCCalROC*  roc) {
157   //
158   // multiply values - per by pad
159   //
160   for (UInt_t  idata = 0; idata< fNChannels; idata++){
161     fData[idata]*=roc->fData[idata];
162   }
163 }
164
165
166 void AliTPCCalROC::Divide(const AliTPCCalROC*  roc) {
167   //
168   // divide values 
169   //
170   Float_t kEpsilon=0.00000000000000001;
171   for (UInt_t  idata = 0; idata< fNChannels; idata++){
172     if (TMath::Abs(roc->fData[idata])>kEpsilon)
173       fData[idata]/=roc->fData[idata];
174   }
175 }
176
177
178
179
180 Double_t AliTPCCalROC::GetLTM(Double_t *sigma, Double_t fraction){
181   //
182   //  Calculate LTM mean and sigma
183   //
184   Double_t *ddata = new Double_t[fNChannels];
185   Double_t mean=0, lsigma=0;
186   Int_t hh = TMath::Min(TMath::Nint(fraction *fNChannels), Int_t(fNChannels));
187   for (UInt_t i=0;i<fNChannels;i++) ddata[i]= fData[i];
188   AliMathBase::EvaluateUni(UInt_t(fNChannels),ddata, mean, lsigma, hh);
189   if (sigma) *sigma=lsigma;
190   delete [] ddata;
191   return mean;
192 }
193
194 TH1F * AliTPCCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type){
195   //
196   // make 1D histo
197   // type -1 = user defined range
198   //       0 = nsigma cut nsigma=min
199   //       1 = delta cut around median delta=min
200   if (type>=0){
201     if (type==0){
202       // nsigma range
203       Float_t mean  = GetMean();
204       Float_t sigma = GetRMS();
205       Float_t nsigma = TMath::Abs(min);
206       min = mean-nsigma*sigma;
207       max = mean+nsigma*sigma;
208     }
209     if (type==1){
210       // fixed range
211       Float_t mean   = GetMedian();
212       Float_t  delta = min;
213       min = mean-delta;
214       max = mean+delta;
215     }
216     if (type==2){
217       //
218       // LTM mean +- nsigma
219       //
220       Double_t sigma;
221       Float_t mean  = GetLTM(&sigma,max);
222       sigma*=min;
223       min = mean-sigma;
224       max = mean+sigma;
225     }
226   }
227   char  name[1000];
228   sprintf(name,"%s ROC 1D%d",GetTitle(),fSector);
229   TH1F * his = new TH1F(name,name,100, min,max);
230   for (UInt_t irow=0; irow<fNRows; irow++){
231     UInt_t npads = (Int_t)GetNPads(irow);
232     for (UInt_t ipad=0; ipad<=npads; ipad++){
233       his->Fill(GetValue(irow,ipad));
234     }
235   }
236   return his;
237 }
238
239
240
241 TH2F * AliTPCCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type){
242   //
243   // make 2D histo
244   // type -1 = user defined range
245   //       0 = nsigma cut nsigma=min
246   //       1 = delta cut around median delta=min
247   if (type>=0){
248     if (type==0){
249       // nsigma range
250       Float_t mean  = GetMean();
251       Float_t sigma = GetRMS();
252       Float_t nsigma = TMath::Abs(min);
253       min = mean-nsigma*sigma;
254       max = mean+nsigma*sigma;
255     }
256     if (type==1){
257       // fixed range
258       Float_t mean   = GetMedian();
259       Float_t  delta = min;
260       min = mean-delta;
261       max = mean+delta;
262     }
263     if (type==2){
264       Double_t sigma;
265       Float_t mean  = GetLTM(&sigma,max);
266       sigma*=min;
267       min = mean-sigma;
268       max = mean+sigma;
269
270     }
271   }
272   UInt_t maxPad = 0;
273   for (UInt_t irow=0; irow<fNRows; irow++){
274     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
275   }
276   char  name[1000];
277   sprintf(name,"%s ROC%d",GetTitle(),fSector);
278   TH2F * his = new TH2F(name,name,fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
279   for (UInt_t irow=0; irow<fNRows; irow++){
280     UInt_t npads = (Int_t)GetNPads(irow);
281     for (UInt_t ipad=0; ipad<=npads; ipad++){
282       his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,GetValue(irow,ipad));
283     }
284   }
285   his->SetMaximum(max);
286   his->SetMinimum(min);
287   return his;
288 }
289
290 TH2F * AliTPCCalROC::MakeHistoOutliers(Float_t delta, Float_t fraction, Int_t type){
291   //
292   // Make Histogram with outliers
293   // mode = 0 - sigma cut used
294   // mode = 1 - absolute cut used
295   // fraction - fraction of values used to define sigma
296   // delta - in mode 0 - nsigma cut
297   //            mode 1 - delta cut
298   Double_t sigma;
299   Float_t mean  = GetLTM(&sigma,fraction);  
300   if (type==0) delta*=sigma; 
301   UInt_t maxPad = 0;
302   for (UInt_t irow=0; irow<fNRows; irow++){
303     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
304   }
305
306   char  name[1000];
307   sprintf(name,"%s ROC Outliers%d",GetTitle(),fSector);
308   TH2F * his = new TH2F(name,name,fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
309   for (UInt_t irow=0; irow<fNRows; irow++){
310     UInt_t npads = (Int_t)GetNPads(irow);
311     for (UInt_t ipad=0; ipad<=npads; ipad++){
312       if (TMath::Abs(GetValue(irow,ipad)-mean)>delta)
313         his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,1);
314     }
315   }
316   return his;
317 }
318
319
320
321 void AliTPCCalROC::Draw(Option_t* opt){
322   //
323   // create histogram with values and draw it
324   //
325   TH1 * his=0; 
326   TString option=opt;
327   option.ToUpper();
328   if (option.Contains("1D")){
329     his = MakeHisto1D();
330   }
331   else{
332     his = MakeHisto2D();
333   }
334   his->Draw(option);
335 }
336
337
338
339
340
341 void AliTPCCalROC::Test(){
342   //
343   // example function to show functionality and tes AliTPCCalROC
344   //
345   AliTPCCalROC  roc0(0);  
346   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
347     for (UInt_t ipad = 0; ipad <roc0.GetNPads(irow); ipad++){
348       Float_t value  = irow+ipad/1000.;
349       roc0.SetValue(irow,ipad,value);
350     }
351   }
352   //
353   AliTPCCalROC roc1(roc0);
354   for (UInt_t irow = 0; irow <roc1.GetNrows(); irow++){
355     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
356       Float_t value  = irow+ipad/1000.;
357       if (roc1.GetValue(irow,ipad)!=value){
358         printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
359       }
360     }
361   }  
362   TFile f("calcTest.root","recreate");
363   roc0.Write("Roc0");
364   AliTPCCalROC * roc2 = (AliTPCCalROC*)f.Get("Roc0");
365   f.Close();
366   //
367   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
368     if (roc0.GetNPads(irow)!=roc2->GetNPads(irow))
369       printf("NPads - Read/Write error\trow=%d\n",irow);
370     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
371       Float_t value  = irow+ipad/1000.;
372       if (roc2->GetValue(irow,ipad)!=value){
373         printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
374       }
375     }
376   }   
377 }
378