coding violations
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraVdriftLinearFit.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 // AliTRDCalibraVdriftLinearFit                                           //
21 //                                                                        //
22 // Does the Vdrift an ExB calibration by applying a linear fit            //
23 //                                                                        //
24 // Author:                                                                //
25 //   R. Bailhache (R.Bailhache@gsi.de)                                    //
26 //                                                                        //
27 ////////////////////////////////////////////////////////////////////////////
28
29 //Root includes
30 #include <TObjArray.h>
31 #include <TH2F.h>
32 #include <TString.h>
33 #include <TVectorD.h>
34 #include <TAxis.h>
35 #include <TLinearFitter.h>
36 #include <TMath.h>
37 #include <TTreeStream.h>
38
39 //header file
40 #include "AliTRDCalibraVdriftLinearFit.h"
41
42 ClassImp(AliTRDCalibraVdriftLinearFit) /*FOLD00*/
43
44 //_____________________________________________________________________
45 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
46   TObject(),
47   fVersion(0),
48   fLinearFitterHistoArray(540),
49   fLinearFitterPArray(540),
50   fLinearFitterEArray(540)
51 {
52     //
53     // default constructor
54     //
55 }
56 //_____________________________________________________________________
57 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /*FOLD00*/
58   TObject(ped),
59   fVersion(ped.fVersion),
60   fLinearFitterHistoArray(540),
61   fLinearFitterPArray(540),
62   fLinearFitterEArray(540)
63 {
64     //
65     // copy constructor
66     //
67   for (Int_t idet = 0; idet < 540; idet++){
68     const TVectorD     *vectorE     = (TVectorD*)ped.fLinearFitterEArray.UncheckedAt(idet);
69     const TVectorD     *vectorP     = (TVectorD*)ped.fLinearFitterPArray.UncheckedAt(idet);
70     const TH2F         *hped        = (TH2F*)ped.fLinearFitterHistoArray.UncheckedAt(idet);
71     
72     if ( vectorE != 0x0 ) fLinearFitterEArray.AddAt(new TVectorD(*vectorE), idet);
73     if ( vectorP != 0x0 ) fLinearFitterPArray.AddAt(new TVectorD(*vectorP), idet);
74     if ( hped != 0x0 ){
75       TH2F *hNew = new TH2F(*hped);
76       hNew->SetDirectory(0);
77       fLinearFitterHistoArray.AddAt(hNew,idet);
78     }
79     
80   }
81 }
82 //_____________________________________________________________________
83 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /*FOLD00*/
84   TObject(),
85   fVersion(0),
86   fLinearFitterHistoArray(540),
87   fLinearFitterPArray(540),
88   fLinearFitterEArray(540)
89 {
90   //
91   // constructor from a TObjArray
92   //
93   for (Int_t idet = 0; idet < 540; idet++){
94     const TH2F         *hped        = (TH2F*)obja.UncheckedAt(idet);
95     if ( hped != 0x0 ){
96       TH2F *hNew = new TH2F(*hped);
97       hNew->SetDirectory(0);
98       fLinearFitterHistoArray.AddAt(hNew,idet);
99     }
100   }
101 }
102 //_____________________________________________________________________
103 AliTRDCalibraVdriftLinearFit& AliTRDCalibraVdriftLinearFit::operator = (const  AliTRDCalibraVdriftLinearFit &source)
104 {
105   //
106   // assignment operator
107   //
108   if (&source == this) return *this;
109   new (this) AliTRDCalibraVdriftLinearFit(source);
110
111   return *this;
112 }
113 //_____________________________________________________________________
114 AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /*FOLD00*/
115 {
116   //
117   // destructor
118   //
119 }
120 //_____________________________________________________________________
121 void AliTRDCalibraVdriftLinearFit::Add(AliTRDCalibraVdriftLinearFit *ped)
122 {
123   //
124   // Add histo
125   //
126
127   fVersion++;
128
129   for (Int_t idet = 0; idet < 540; idet++){
130     const TH2F         *hped        = (TH2F*)ped->GetLinearFitterHisto(idet);
131     //printf("idet %d\n",idet);
132     if ( hped != 0x0 ){
133       //printf("add\n");
134       TH2F *hped1 = GetLinearFitterHisto(idet,kTRUE);
135       //printf("test2\n");
136       hped1->SetDirectory(0);
137       //printf("test4\n");
138       hped1->Add(hped);
139       //printf("test3\n");
140     }
141   }
142 }
143 //______________________________________________________________________________________
144 TH2F* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
145 {
146     //
147     // return pointer to TH2F histo 
148     // if force is true create a new histo if it doesn't exist allready
149     //
150     if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
151         return (TH2F*)fLinearFitterHistoArray.UncheckedAt(detector);
152
153     // if we are forced and TLinearFitter doesn't yes exist create it
154
155     // new TH2F
156     TString name("LFDV");
157     name += detector;
158     name += "version";
159     name +=  fVersion;
160     
161     TH2F *lfdv = new TH2F((const Char_t *)name,(const Char_t *) name
162                           ,100,-1.0,1.0,100
163                           ,-2.0,2.0);
164     lfdv->SetXTitle("tan(phi_{track})");
165     lfdv->SetYTitle("dy/dt");
166     lfdv->SetZTitle("Number of clusters");
167     lfdv->SetStats(0);
168     lfdv->SetDirectory(0);
169
170     fLinearFitterHistoArray.AddAt(lfdv,detector);
171     return lfdv;
172 }
173 //______________________________________________________________________________________
174 Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
175 {
176     //
177     // return param for this detector
178     //
179   if ( fLinearFitterPArray.UncheckedAt(detector) ){
180     const TVectorD     *vectorP     = (TVectorD*)fLinearFitterPArray.UncheckedAt(detector);
181     if(!param) param = new TVectorD(2);
182     for(Int_t k = 0; k < 2; k++){
183       (*param)[k] = (*vectorP)[k];
184     }
185     return kTRUE;
186   }
187   else return kFALSE;
188
189 }
190 //______________________________________________________________________________________
191 Bool_t AliTRDCalibraVdriftLinearFit::GetError(Int_t detector, TVectorD *error)
192 {
193     //
194     // return error for this detector 
195     //
196   if ( fLinearFitterEArray.UncheckedAt(detector) ){
197     const TVectorD     *vectorE     = (TVectorD*)fLinearFitterEArray.UncheckedAt(detector);
198     if(!error) error = new TVectorD(3);
199     for(Int_t k = 0; k < 3; k++){
200       (*error)[k] = (*vectorE)[k];
201     }
202     return kTRUE;
203   }
204   else return kFALSE;
205
206 }
207 //______________________________________________________________________________________
208 void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
209 {
210     //
211     // Fill the 2D histos for debugging
212     //
213
214   ((TH2F *) GetLinearFitterHisto(detector,kTRUE))->Fill(tnp,pars1);
215
216 }
217 //____________Functions fit Online CH2d________________________________________
218 void AliTRDCalibraVdriftLinearFit::FillPEArray()
219 {
220   //
221   // Fill fLinearFitterPArray and fLinearFitterEArray from inside
222   //
223
224   
225   Int_t *arrayI = new Int_t[540];
226   for(Int_t k = 0; k< 540; k++){
227     arrayI[k] = 0; 
228   }
229
230   // Loop over histos 
231   for(Int_t cb = 0; cb < 540; cb++){
232     const TH2F *linearfitterhisto = (TH2F*)fLinearFitterHistoArray.UncheckedAt(cb);
233     //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);    
234
235     if ( linearfitterhisto != 0 ){
236       
237       // Fill a linearfitter
238       TAxis *xaxis = linearfitterhisto->GetXaxis();
239       TAxis *yaxis = linearfitterhisto->GetYaxis();
240       TLinearFitter linearfitter = TLinearFitter(2,"pol1");
241       //printf("test\n");
242       for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
243         for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
244           if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
245             Double_t x = xaxis->GetBinCenter(ibinx+1);
246             Double_t y = yaxis->GetBinCenter(ibiny+1);
247             for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
248               linearfitter.AddPoint(&x,y);
249               arrayI[cb]++;
250             }
251           }
252         }
253       }
254       
255       //printf("Find %d entries for the detector %d\n",arrayI[cb],cb);
256
257       // Eval the linear fitter
258       if(arrayI[cb]>10){
259         TVectorD  *par  = new TVectorD(2);
260         TVectorD   pare = TVectorD(2);
261         TVectorD  *parE = new TVectorD(3);
262         linearfitter.Eval();
263         linearfitter.GetParameters(*par);
264         linearfitter.GetErrors(pare);
265         Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter.GetChisquare())/arrayI[cb]);
266         (*parE)[0] = pare[0]*ppointError;
267         (*parE)[1] = pare[1]*ppointError;
268         (*parE)[2] = (Double_t) arrayI[cb];
269         fLinearFitterPArray.AddAt(par,cb);
270         fLinearFitterEArray.AddAt(parE,cb);
271         //par->Print();
272         //parE->Print();
273       }    
274     }// if something
275   }
276    
277 }