Update of the calibration code for PbPb (Raphaelle)..
[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
40 //header file
41 #include "AliTRDCalibraVdriftLinearFit.h"
42
43 ClassImp(AliTRDCalibraVdriftLinearFit) /*FOLD00*/
44
45 //_____________________________________________________________________
46 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
47   TObject(),
48   fVersion(0),
49   fLinearFitterHistoArray(540),
50   fLinearFitterPArray(540),
51   fLinearFitterEArray(540)
52 {
53   //
54   // default constructor
55   //
56 }
57 //_____________________________________________________________________
58 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /*FOLD00*/
59   TObject(ped),
60   fVersion(ped.fVersion),
61   fLinearFitterHistoArray(540),
62   fLinearFitterPArray(540),
63   fLinearFitterEArray(540)
64 {
65     //
66     // copy constructor
67     //
68   for (Int_t idet = 0; idet < 540; idet++){
69    
70     const TVectorD     *vectorE     = (TVectorD*)ped.fLinearFitterEArray.UncheckedAt(idet);
71     const TVectorD     *vectorP     = (TVectorD*)ped.fLinearFitterPArray.UncheckedAt(idet);
72     const TH2S         *hped        = (TH2S*)ped.fLinearFitterHistoArray.UncheckedAt(idet);
73     
74     if ( vectorE != 0x0 ) fLinearFitterEArray.AddAt(new TVectorD(*vectorE), idet);
75     if ( vectorP != 0x0 ) fLinearFitterPArray.AddAt(new TVectorD(*vectorP), idet);
76     if ( hped != 0x0 ){
77       TH2S *hNew = (TH2S *)hped->Clone();
78       //hNew->SetDirectory(0);
79       fLinearFitterHistoArray.AddAt(hNew,idet);
80     }
81   }
82 }
83 //_____________________________________________________________________
84 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /*FOLD00*/
85   TObject(),
86   fVersion(0),
87   fLinearFitterHistoArray(540),
88   fLinearFitterPArray(540),
89   fLinearFitterEArray(540)
90 {
91   //
92   // constructor from a TObjArray
93   //
94   for (Int_t idet = 0; idet < 540; idet++){
95     const TH2S         *hped        = (TH2S*)obja.UncheckedAt(idet);
96     if ( hped != 0x0 ){
97       TH2S *hNew = (TH2S *)hped->Clone();
98       //hNew->SetDirectory(0);
99       fLinearFitterHistoArray.AddAt(hNew,idet);
100     }
101   }
102 }
103 //_____________________________________________________________________
104 AliTRDCalibraVdriftLinearFit& AliTRDCalibraVdriftLinearFit::operator = (const  AliTRDCalibraVdriftLinearFit &source)
105 {
106   //
107   // assignment operator
108   //
109   if (&source == this) return *this;
110   new (this) AliTRDCalibraVdriftLinearFit(source);
111
112   return *this;
113 }
114 //_____________________________________________________________________
115 AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /*FOLD00*/
116 {
117   //
118   // destructor
119   //
120   fLinearFitterHistoArray.SetOwner();
121   fLinearFitterPArray.SetOwner();
122   fLinearFitterEArray.SetOwner();
123
124   fLinearFitterHistoArray.Delete();
125   fLinearFitterPArray.Delete();
126   fLinearFitterEArray.Delete();
127
128 }
129 //_____________________________________________________________________________
130 void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const
131 {
132   //
133   // Copy function
134   //
135
136   AliTRDCalibraVdriftLinearFit& target = (AliTRDCalibraVdriftLinearFit &) c;
137
138   // Copy only the histos
139   for (Int_t idet = 0; idet < 540; idet++){
140     if(fLinearFitterHistoArray.UncheckedAt(idet)){
141       TH2S *hped1 = (TH2S *)target.GetLinearFitterHisto(idet,kTRUE);
142       //hped1->SetDirectory(0);
143       hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet));
144     }
145   }
146   
147   TObject::Copy(c);
148
149 }
150 //_____________________________________________________________________________
151 Long64_t AliTRDCalibraVdriftLinearFit::Merge(const TCollection* list) 
152 {
153   // Merge list of objects (needed by PROOF)
154
155   if (!list)
156     return 0;
157   
158   if (list->IsEmpty())
159     return 1;
160   
161   TIterator* iter = list->MakeIterator();
162   TObject* obj = 0;
163   
164   // collection of generated histograms
165   Int_t count=0;
166   while((obj = iter->Next()) != 0) 
167     {
168       AliTRDCalibraVdriftLinearFit* entry = dynamic_cast<AliTRDCalibraVdriftLinearFit*>(obj);
169       if (entry == 0) continue; 
170       
171       // Copy only the histos
172       for (Int_t idet = 0; idet < 540; idet++){
173         if(entry->GetLinearFitterHisto(idet)){
174           TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
175           Double_t entriesa = hped1->GetEntries();
176           Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries();
177           if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet));
178         }
179       }
180       
181       count++;
182     }
183   
184   return count;
185 }
186 //_____________________________________________________________________
187 void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped)
188 {
189   //
190   // Add histo
191   //
192
193   fVersion++;
194
195   for (Int_t idet = 0; idet < 540; idet++){
196     const TH2S         *hped        = (TH2S*)ped->GetLinearFitterHistoNoForce(idet);
197     //printf("idet %d\n",idet);
198     if ( hped != 0x0 ){
199       //printf("add\n");
200       TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
201       Double_t entriesa = hped1->GetEntries();
202       Double_t entriesb = hped->GetEntries();
203       if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
204     }
205   }
206 }
207 //______________________________________________________________________________________
208 TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
209 {
210     //
211     // return pointer to TH2F histo 
212     // if force is true create a new histo if it doesn't exist allready
213     //
214     if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
215         return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
216
217     return GetLinearFitterHistoForce(detector);
218
219 }
220 //______________________________________________________________________________________
221 TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHistoForce(Int_t detector)
222 {
223   //
224   // return pointer to TH2F histo 
225   // if NULL create a new histo if it doesn't exist allready
226   //
227   if (fLinearFitterHistoArray.UncheckedAt(detector))
228     return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
229   
230   // if we are forced and TLinearFitter doesn't yes exist create it
231   
232   // new TH2F
233   TString name("LFDV");
234   name += detector;
235   name += "version";
236   name +=  fVersion;
237   
238   TH2S *lfdv = new TH2S((const Char_t *)name,(const Char_t *) name
239                         ,36,-0.9,0.9,48
240                         ,-1.2,1.2);
241   lfdv->SetXTitle("tan(phi_{track})");
242   lfdv->SetYTitle("dy/dt");
243   lfdv->SetZTitle("Number of clusters");
244   lfdv->SetStats(0);
245   lfdv->SetDirectory(0);
246   
247   fLinearFitterHistoArray.AddAt(lfdv,detector);
248   return lfdv;
249 }
250 //______________________________________________________________________________________
251 Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
252 {
253     //
254     // return param for this detector
255     //
256   if ( fLinearFitterPArray.UncheckedAt(detector) ){
257     const TVectorD     *vectorP     = (TVectorD*)fLinearFitterPArray.UncheckedAt(detector);
258     if(!param) param = new TVectorD(2);
259     for(Int_t k = 0; k < 2; k++){
260       (*param)[k] = (*vectorP)[k];
261     }
262     return kTRUE;
263   }
264   else return kFALSE;
265
266 }
267 //______________________________________________________________________________________
268 Bool_t AliTRDCalibraVdriftLinearFit::GetError(Int_t detector, TVectorD *error)
269 {
270     //
271     // return error for this detector 
272     //
273   if ( fLinearFitterEArray.UncheckedAt(detector) ){
274     const TVectorD     *vectorE     = (TVectorD*)fLinearFitterEArray.UncheckedAt(detector);
275     if(!error) error = new TVectorD(3);
276     for(Int_t k = 0; k < 3; k++){
277       (*error)[k] = (*vectorE)[k];
278     }
279     return kTRUE;
280   }
281   else return kFALSE;
282
283 }
284 //______________________________________________________________________________________
285 void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
286 {
287     //
288     // Fill the 2D histos for debugging
289     //
290   
291   TH2S *h = ((TH2S *) GetLinearFitterHisto(detector,kTRUE));
292   Double_t nbentries = h->GetEntries();
293   if(nbentries < 5*32767) h->Fill(tnp,pars1);
294
295 }
296 //____________Functions fit Online CH2d________________________________________
297 void AliTRDCalibraVdriftLinearFit::FillPEArray()
298 {
299   //
300   // Fill fLinearFitterPArray and fLinearFitterEArray from inside
301   //
302
303   
304   Int_t *arrayI = new Int_t[540];
305   for(Int_t k = 0; k< 540; k++){
306     arrayI[k] = 0; 
307   }
308
309   // Loop over histos 
310   for(Int_t cb = 0; cb < 540; cb++){
311     const TH2S *linearfitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
312     //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);    
313
314     if ( linearfitterhisto != 0 ){
315       
316       // Fill a linearfitter
317       TAxis *xaxis = linearfitterhisto->GetXaxis();
318       TAxis *yaxis = linearfitterhisto->GetYaxis();
319       TLinearFitter linearfitter = TLinearFitter(2,"pol1");
320       //printf("test\n");
321       for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
322         for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
323           if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
324             Double_t x = xaxis->GetBinCenter(ibinx+1);
325             Double_t y = yaxis->GetBinCenter(ibiny+1);
326             for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
327               linearfitter.AddPoint(&x,y);
328               arrayI[cb]++;
329             }
330           }
331         }
332       }
333       
334       //printf("Find %d entries for the detector %d\n",arrayI[cb],cb);
335
336       // Eval the linear fitter
337       if(arrayI[cb]>10){
338         TVectorD  *par  = new TVectorD(2);
339         TVectorD   pare = TVectorD(2);
340         TVectorD  *parE = new TVectorD(3);
341         if((linearfitter.EvalRobust(0.8)==0)) {
342         //if((linearfitter.Eval()==0)) {
343           
344           linearfitter.GetParameters(*par);
345           //linearfitter.GetErrors(pare);
346           //Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter.GetChisquare())/arrayI[cb]);
347           //(*parE)[0] = pare[0]*ppointError;
348           //(*parE)[1] = pare[1]*ppointError;
349           
350           (*parE)[0] = 0.0;
351           (*parE)[1] = 0.0;
352           (*parE)[2] = (Double_t) arrayI[cb];
353           fLinearFitterPArray.AddAt(par,cb);
354           fLinearFitterEArray.AddAt(parE,cb);
355           
356           //par->Print();
357           //parE->Print();
358         }
359       }
360       
361       //delete linearfitterhisto;
362       
363     }// if something
364
365   }
366
367   delete [] arrayI;
368    
369 }
370
371