Small update (number of bins)
[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   fLinearFitterHistoArray.Delete();
120   fLinearFitterPArray.Delete();
121   fLinearFitterEArray.Delete();
122
123 }
124 //_____________________________________________________________________________
125 void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const
126 {
127   //
128   // Copy function
129   //
130
131   AliTRDCalibraVdriftLinearFit& target = (AliTRDCalibraVdriftLinearFit &) c;
132
133   // Copy only the histos
134   for (Int_t idet = 0; idet < 540; idet++){
135     if(fLinearFitterHistoArray.UncheckedAt(idet)){
136       TH2F *hped1 = target.GetLinearFitterHisto(idet,kTRUE);
137       hped1->SetDirectory(0);
138       hped1->Add((const TH1F *)fLinearFitterHistoArray.UncheckedAt(idet));
139     }
140   
141   }
142
143   TObject::Copy(c);
144
145 }
146 //_____________________________________________________________________________
147 Long64_t AliTRDCalibraVdriftLinearFit::Merge(TCollection* list) 
148 {
149   // Merge list of objects (needed by PROOF)
150
151   if (!list)
152     return 0;
153   
154   if (list->IsEmpty())
155     return 1;
156   
157   TIterator* iter = list->MakeIterator();
158   TObject* obj = 0;
159   
160   // collection of generated histograms
161   Int_t count=0;
162   while((obj = iter->Next()) != 0) 
163   {
164     AliTRDCalibraVdriftLinearFit* entry = dynamic_cast<AliTRDCalibraVdriftLinearFit*>(obj);
165     if (entry == 0) continue; 
166
167     // Copy only the histos
168     for (Int_t idet = 0; idet < 540; idet++){
169       if(entry->GetLinearFitterHisto(idet)){
170         TH2F *hped1 = GetLinearFitterHisto(idet,kTRUE);
171         hped1->SetDirectory(0);
172         hped1->Add(entry->GetLinearFitterHisto(idet));
173       }
174       
175     }
176     
177     count++;
178   }
179   
180   return count;
181 }
182 //_____________________________________________________________________
183 void AliTRDCalibraVdriftLinearFit::Add(AliTRDCalibraVdriftLinearFit *ped)
184 {
185   //
186   // Add histo
187   //
188
189   fVersion++;
190
191   for (Int_t idet = 0; idet < 540; idet++){
192     const TH2F         *hped        = (TH2F*)ped->GetLinearFitterHisto(idet);
193     //printf("idet %d\n",idet);
194     if ( hped != 0x0 ){
195       //printf("add\n");
196       TH2F *hped1 = GetLinearFitterHisto(idet,kTRUE);
197       //printf("test2\n");
198       hped1->SetDirectory(0);
199       //printf("test4\n");
200       hped1->Add(hped);
201       //printf("test3\n");
202     }
203   }
204 }
205 //______________________________________________________________________________________
206 TH2F* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
207 {
208     //
209     // return pointer to TH2F histo 
210     // if force is true create a new histo if it doesn't exist allready
211     //
212     if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
213         return (TH2F*)fLinearFitterHistoArray.UncheckedAt(detector);
214
215     // if we are forced and TLinearFitter doesn't yes exist create it
216
217     // new TH2F
218     TString name("LFDV");
219     name += detector;
220     name += "version";
221     name +=  fVersion;
222     
223     TH2F *lfdv = new TH2F((const Char_t *)name,(const Char_t *) name
224                           ,40,-1.0,1.0,50
225                           ,-2.0,2.0);
226     lfdv->SetXTitle("tan(phi_{track})");
227     lfdv->SetYTitle("dy/dt");
228     lfdv->SetZTitle("Number of clusters");
229     lfdv->SetStats(0);
230     lfdv->SetDirectory(0);
231
232     fLinearFitterHistoArray.AddAt(lfdv,detector);
233     return lfdv;
234 }
235 //______________________________________________________________________________________
236 Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
237 {
238     //
239     // return param for this detector
240     //
241   if ( fLinearFitterPArray.UncheckedAt(detector) ){
242     const TVectorD     *vectorP     = (TVectorD*)fLinearFitterPArray.UncheckedAt(detector);
243     if(!param) param = new TVectorD(2);
244     for(Int_t k = 0; k < 2; k++){
245       (*param)[k] = (*vectorP)[k];
246     }
247     return kTRUE;
248   }
249   else return kFALSE;
250
251 }
252 //______________________________________________________________________________________
253 Bool_t AliTRDCalibraVdriftLinearFit::GetError(Int_t detector, TVectorD *error)
254 {
255     //
256     // return error for this detector 
257     //
258   if ( fLinearFitterEArray.UncheckedAt(detector) ){
259     const TVectorD     *vectorE     = (TVectorD*)fLinearFitterEArray.UncheckedAt(detector);
260     if(!error) error = new TVectorD(3);
261     for(Int_t k = 0; k < 3; k++){
262       (*error)[k] = (*vectorE)[k];
263     }
264     return kTRUE;
265   }
266   else return kFALSE;
267
268 }
269 //______________________________________________________________________________________
270 void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
271 {
272     //
273     // Fill the 2D histos for debugging
274     //
275
276   ((TH2F *) GetLinearFitterHisto(detector,kTRUE))->Fill(tnp,pars1);
277
278 }
279 //____________Functions fit Online CH2d________________________________________
280 void AliTRDCalibraVdriftLinearFit::FillPEArray()
281 {
282   //
283   // Fill fLinearFitterPArray and fLinearFitterEArray from inside
284   //
285
286   
287   Int_t *arrayI = new Int_t[540];
288   for(Int_t k = 0; k< 540; k++){
289     arrayI[k] = 0; 
290   }
291
292   // Loop over histos 
293   for(Int_t cb = 0; cb < 540; cb++){
294     const TH2F *linearfitterhisto = (TH2F*)fLinearFitterHistoArray.UncheckedAt(cb);
295     //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);    
296
297     if ( linearfitterhisto != 0 ){
298       
299       // Fill a linearfitter
300       TAxis *xaxis = linearfitterhisto->GetXaxis();
301       TAxis *yaxis = linearfitterhisto->GetYaxis();
302       TLinearFitter linearfitter = TLinearFitter(2,"pol1");
303       //printf("test\n");
304       for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
305         for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
306           if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
307             Double_t x = xaxis->GetBinCenter(ibinx+1);
308             Double_t y = yaxis->GetBinCenter(ibiny+1);
309             for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
310               linearfitter.AddPoint(&x,y);
311               arrayI[cb]++;
312             }
313           }
314         }
315       }
316       
317       //printf("Find %d entries for the detector %d\n",arrayI[cb],cb);
318
319       // Eval the linear fitter
320       if(arrayI[cb]>10){
321         TVectorD  *par  = new TVectorD(2);
322         TVectorD   pare = TVectorD(2);
323         TVectorD  *parE = new TVectorD(3);
324         linearfitter.Eval();
325         linearfitter.GetParameters(*par);
326         linearfitter.GetErrors(pare);
327         Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter.GetChisquare())/arrayI[cb]);
328         (*parE)[0] = pare[0]*ppointError;
329         (*parE)[1] = pare[1]*ppointError;
330         (*parE)[2] = (Double_t) arrayI[cb];
331         fLinearFitterPArray.AddAt(par,cb);
332         fLinearFitterEArray.AddAt(parE,cb);
333         //par->Print();
334         //parE->Print();
335       }    
336     }// if something
337   }
338    
339 }
340