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