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