Move from THnSparseS to TH2S and changes due to shared pads and baseline subtraction
[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.Delete();
121   fLinearFitterPArray.Delete();
122   fLinearFitterEArray.Delete();
123
124 }
125 //_____________________________________________________________________________
126 void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const
127 {
128   //
129   // Copy function
130   //
131
132   AliTRDCalibraVdriftLinearFit& target = (AliTRDCalibraVdriftLinearFit &) c;
133
134   // Copy only the histos
135   for (Int_t idet = 0; idet < 540; idet++){
136     if(fLinearFitterHistoArray.UncheckedAt(idet)){
137       TH2S *hped1 = (TH2S *)target.GetLinearFitterHisto(idet,kTRUE);
138       //hped1->SetDirectory(0);
139       hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet));
140     }
141   }
142   
143   TObject::Copy(c);
144
145 }
146 //_____________________________________________________________________________
147 Long64_t AliTRDCalibraVdriftLinearFit::Merge(const 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           TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
171           Double_t entriesa = hped1->GetEntries();
172           Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries();
173           if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet));
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 TH2S         *hped        = (TH2S*)ped->GetLinearFitterHisto(idet);
193     //printf("idet %d\n",idet);
194     if ( hped != 0x0 ){
195       //printf("add\n");
196       TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
197       Double_t entriesa = hped1->GetEntries();
198       Double_t entriesb = hped->GetEntries();
199       if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
200     }
201   }
202 }
203 //______________________________________________________________________________________
204 TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
205 {
206     //
207     // return pointer to TH2F histo 
208     // if force is true create a new histo if it doesn't exist allready
209     //
210     if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
211         return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
212
213     // if we are forced and TLinearFitter doesn't yes exist create it
214
215     // new TH2F
216     TString name("LFDV");
217     name += detector;
218     name += "version";
219     name +=  fVersion;
220
221     TH2S *lfdv = new TH2S((const Char_t *)name,(const Char_t *) name
222                   ,36,-0.9,0.9,48
223                   ,-1.2,1.2);
224     lfdv->SetXTitle("tan(phi_{track})");
225     lfdv->SetYTitle("dy/dt");
226     lfdv->SetZTitle("Number of clusters");
227     lfdv->SetStats(0);
228     lfdv->SetDirectory(0);
229
230     fLinearFitterHistoArray.AddAt(lfdv,detector);
231     return lfdv;
232 }
233 //______________________________________________________________________________________
234 Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
235 {
236     //
237     // return param for this detector
238     //
239   if ( fLinearFitterPArray.UncheckedAt(detector) ){
240     const TVectorD     *vectorP     = (TVectorD*)fLinearFitterPArray.UncheckedAt(detector);
241     if(!param) param = new TVectorD(2);
242     for(Int_t k = 0; k < 2; k++){
243       (*param)[k] = (*vectorP)[k];
244     }
245     return kTRUE;
246   }
247   else return kFALSE;
248
249 }
250 //______________________________________________________________________________________
251 Bool_t AliTRDCalibraVdriftLinearFit::GetError(Int_t detector, TVectorD *error)
252 {
253     //
254     // return error for this detector 
255     //
256   if ( fLinearFitterEArray.UncheckedAt(detector) ){
257     const TVectorD     *vectorE     = (TVectorD*)fLinearFitterEArray.UncheckedAt(detector);
258     if(!error) error = new TVectorD(3);
259     for(Int_t k = 0; k < 3; k++){
260       (*error)[k] = (*vectorE)[k];
261     }
262     return kTRUE;
263   }
264   else return kFALSE;
265
266 }
267 //______________________________________________________________________________________
268 void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
269 {
270     //
271     // Fill the 2D histos for debugging
272     //
273   
274   TH2S *h = ((TH2S *) GetLinearFitterHisto(detector,kTRUE));
275   Double_t nbentries = h->GetEntries();
276   if(nbentries < 5*32767) h->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 TH2S *linearfitterhisto = (TH2S*)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
337       delete linearfitterhisto;
338       
339     }// if something
340   }
341    
342 }
343