Adding statistical function (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 Double_t AliTPCCalROC::GetLTM(Double_t *sigma, Double_t fraction){
128   //
129   //  Calculate LTM mean and sigma
130   //
131   Double_t *ddata = new Double_t[fNChannels];
132   Double_t mean=0, lsigma=0;
133   Int_t hh = TMath::Min(TMath::Nint(fraction *fNChannels), Int_t(fNChannels));
134   for (UInt_t i=0;i<fNChannels;i++) ddata[i]= fData[i];
135   AliMathBase::EvaluateUni(UInt_t(fNChannels),ddata, mean, lsigma, hh);
136   if (sigma) *sigma=lsigma;
137   delete [] ddata;
138   return mean;
139 }
140
141 TH1F * AliTPCCalROC::MakeHisto1D(Float_t min, Float_t max,Int_t type){
142   //
143   // make 1D histo
144   // type -1 = user defined range
145   //       0 = nsigma cut nsigma=min
146   //       1 = delta cut around median delta=min
147   if (type>=0){
148     if (type==0){
149       // nsigma range
150       Float_t mean  = GetMean();
151       Float_t sigma = GetRMS();
152       Float_t nsigma = TMath::Abs(min);
153       min = mean-nsigma*sigma;
154       max = mean+nsigma*sigma;
155     }
156     if (type==1){
157       // fixed range
158       Float_t mean   = GetMedian();
159       Float_t  delta = min;
160       min = mean-delta;
161       max = mean+delta;
162     }
163     if (type==2){
164       //
165       // LTM mean +- nsigma
166       //
167       Double_t sigma;
168       Float_t mean  = GetLTM(&sigma,max);
169       sigma*=min;
170       min = mean-sigma;
171       max = mean+sigma;
172     }
173   }
174   char  name[1000];
175   sprintf(name,"%s ROC 1D%d",GetTitle(),fSector);
176   TH1F * his = new TH1F(name,name,100, min,max);
177   for (UInt_t irow=0; irow<fNRows; irow++){
178     UInt_t npads = (Int_t)GetNPads(irow);
179     for (UInt_t ipad=0; ipad<=npads; ipad++){
180       his->Fill(GetValue(irow,ipad));
181     }
182   }
183   return his;
184 }
185
186
187
188 TH2F * AliTPCCalROC::MakeHisto2D(Float_t min, Float_t max,Int_t type){
189   //
190   // make 2D histo
191   // type -1 = user defined range
192   //       0 = nsigma cut nsigma=min
193   //       1 = delta cut around median delta=min
194   if (type>=0){
195     if (type==0){
196       // nsigma range
197       Float_t mean  = GetMean();
198       Float_t sigma = GetRMS();
199       Float_t nsigma = TMath::Abs(min);
200       min = mean-nsigma*sigma;
201       max = mean+nsigma*sigma;
202     }
203     if (type==1){
204       // fixed range
205       Float_t mean   = GetMedian();
206       Float_t  delta = min;
207       min = mean-delta;
208       max = mean+delta;
209     }
210     if (type==2){
211       Double_t sigma;
212       Float_t mean  = GetLTM(&sigma,max);
213       sigma*=min;
214       min = mean-sigma;
215       max = mean+sigma;
216
217     }
218   }
219   UInt_t maxPad = 0;
220   for (UInt_t irow=0; irow<fNRows; irow++){
221     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
222   }
223   char  name[1000];
224   sprintf(name,"%s ROC%d",GetTitle(),fSector);
225   TH2F * his = new TH2F(name,name,fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
226   for (UInt_t irow=0; irow<fNRows; irow++){
227     UInt_t npads = (Int_t)GetNPads(irow);
228     for (UInt_t ipad=0; ipad<=npads; ipad++){
229       his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,GetValue(irow,ipad));
230     }
231   }
232   his->SetMaximum(max);
233   his->SetMinimum(min);
234   return his;
235 }
236
237 TH2F * AliTPCCalROC::MakeHistoOutliers(Float_t delta, Float_t fraction, Int_t type){
238   //
239   // Make Histogram with outliers
240   // mode = 0 - sigma cut used
241   // mode = 1 - absolute cut used
242   // fraction - fraction of values used to define sigma
243   // delta - in mode 0 - nsigma cut
244   //            mode 1 - delta cut
245   Double_t sigma;
246   Float_t mean  = GetLTM(&sigma,fraction);  
247   if (type==0) delta*=sigma; 
248   UInt_t maxPad = 0;
249   for (UInt_t irow=0; irow<fNRows; irow++){
250     if (GetNPads(irow)>maxPad) maxPad = GetNPads(irow);
251   }
252
253   char  name[1000];
254   sprintf(name,"%s ROC Outliers%d",GetTitle(),fSector);
255   TH2F * his = new TH2F(name,name,fNRows+10,-5, fNRows+5, maxPad+10, -(Int_t(maxPad/2))-5, maxPad/2+5);
256   for (UInt_t irow=0; irow<fNRows; irow++){
257     UInt_t npads = (Int_t)GetNPads(irow);
258     for (UInt_t ipad=0; ipad<=npads; ipad++){
259       if (TMath::Abs(GetValue(irow,ipad)-mean)>delta)
260         his->Fill(irow+0.5,Int_t(ipad)-Int_t(npads/2)+0.5,1);
261     }
262   }
263   return his;
264 }
265
266
267
268 void AliTPCCalROC::Draw(Option_t* opt){
269   //
270   // create histogram with values and draw it
271   //
272   TH1 * his=0; 
273   TString option=opt;
274   option.ToUpper();
275   if (option.Contains("1D")){
276     his = MakeHisto1D();
277   }
278   else{
279     his = MakeHisto2D();
280   }
281   his->Draw(option);
282 }
283
284
285
286
287
288 void AliTPCCalROC::Test(){
289   //
290   // example function to show functionality and tes AliTPCCalROC
291   //
292   AliTPCCalROC  roc0(0);  
293   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
294     for (UInt_t ipad = 0; ipad <roc0.GetNPads(irow); ipad++){
295       Float_t value  = irow+ipad/1000.;
296       roc0.SetValue(irow,ipad,value);
297     }
298   }
299   //
300   AliTPCCalROC roc1(roc0);
301   for (UInt_t irow = 0; irow <roc1.GetNrows(); irow++){
302     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
303       Float_t value  = irow+ipad/1000.;
304       if (roc1.GetValue(irow,ipad)!=value){
305         printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
306       }
307     }
308   }  
309   TFile f("calcTest.root","recreate");
310   roc0.Write("Roc0");
311   AliTPCCalROC * roc2 = (AliTPCCalROC*)f.Get("Roc0");
312   f.Close();
313   //
314   for (UInt_t irow = 0; irow <roc0.GetNrows(); irow++){
315     if (roc0.GetNPads(irow)!=roc2->GetNPads(irow))
316       printf("NPads - Read/Write error\trow=%d\n",irow);
317     for (UInt_t ipad = 0; ipad <roc1.GetNPads(irow); ipad++){
318       Float_t value  = irow+ipad/1000.;
319       if (roc2->GetValue(irow,ipad)!=value){
320         printf("Read/Write error\trow=%d\tpad=%d\n",irow,ipad);
321       }
322     }
323   }   
324 }
325