changes from fzhou
[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 <TStyle.h>
38 #include <TCanvas.h>
39 #include <TTreeStream.h>
40 #include <TGraphErrors.h>
41 #include <TDirectory.h>
42 #include <TTreeStream.h>
43 #include <TF1.h>
44
45 //header file
46 #include "AliTRDCalibraVdriftLinearFit.h"
47
48 ClassImp(AliTRDCalibraVdriftLinearFit) /*FOLD00*/
49
50 //_____________________________________________________________________
51 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
52   TObject(),
53   fVersion(0),
54   fNameCalibUsed(""),
55   fLinearFitterHistoArray(540),
56   fLinearFitterPArray(540),
57   fLinearFitterEArray(540),
58   fRobustFit(kTRUE),
59   fMinNpointsFit(11),
60   fDebugStreamer(0x0),
61   fDebugLevel(0),
62   fSeeDetector(0)
63 {
64   //
65   // default constructor
66   //
67 }
68 //_____________________________________________________________________
69 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /*FOLD00*/
70   TObject(ped),
71   fVersion(ped.fVersion),
72   fNameCalibUsed(ped.fNameCalibUsed),
73   fLinearFitterHistoArray(540),
74   fLinearFitterPArray(540),
75   fLinearFitterEArray(540),
76   fRobustFit(kTRUE),
77   fMinNpointsFit(10),
78   fDebugStreamer(0x0),
79   fDebugLevel(0),
80   fSeeDetector(0)
81 {
82     //
83     // copy constructor
84     //
85   for (Int_t idet = 0; idet < 540; idet++){
86    
87     const TVectorD     *vectorE     = (TVectorD*)ped.fLinearFitterEArray.UncheckedAt(idet);
88     const TVectorD     *vectorP     = (TVectorD*)ped.fLinearFitterPArray.UncheckedAt(idet);
89     const TH2S         *hped        = (TH2S*)ped.fLinearFitterHistoArray.UncheckedAt(idet);
90     
91     if ( vectorE != 0x0 ) fLinearFitterEArray.AddAt(new TVectorD(*vectorE), idet);
92     if ( vectorP != 0x0 ) fLinearFitterPArray.AddAt(new TVectorD(*vectorP), idet);
93     if ( hped != 0x0 ){
94       TH2S *hNew = (TH2S *)hped->Clone();
95       //hNew->SetDirectory(0);
96       fLinearFitterHistoArray.AddAt(hNew,idet);
97     }
98   }
99 }
100 //_____________________________________________________________________
101 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /*FOLD00*/
102   TObject(),
103   fVersion(0),
104   fNameCalibUsed(""),
105   fLinearFitterHistoArray(540),
106   fLinearFitterPArray(540),
107   fLinearFitterEArray(540),
108   fRobustFit(kTRUE),
109   fMinNpointsFit(10),
110   fDebugStreamer(0x0),
111   fDebugLevel(0),
112   fSeeDetector(0)
113 {
114   //
115   // constructor from a TObjArray
116   //
117   for (Int_t idet = 0; idet < 540; idet++){
118     const TH2S         *hped        = (TH2S*)obja.UncheckedAt(idet);
119     if ( hped != 0x0 ){
120       TH2S *hNew = (TH2S *)hped->Clone();
121       //hNew->SetDirectory(0);
122       fLinearFitterHistoArray.AddAt(hNew,idet);
123     }
124   }
125 }
126 //_____________________________________________________________________
127 AliTRDCalibraVdriftLinearFit& AliTRDCalibraVdriftLinearFit::operator = (const  AliTRDCalibraVdriftLinearFit &source)
128 {
129   //
130   // assignment operator
131   //
132   if (&source == this) return *this;
133   new (this) AliTRDCalibraVdriftLinearFit(source);
134
135   return *this;
136 }
137 //_____________________________________________________________________
138 AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /*FOLD00*/
139 {
140   //
141   // destructor
142   //
143   fLinearFitterHistoArray.SetOwner();
144   fLinearFitterPArray.SetOwner();
145   fLinearFitterEArray.SetOwner();
146
147   fLinearFitterHistoArray.Delete();
148   fLinearFitterPArray.Delete();
149   fLinearFitterEArray.Delete();
150
151   if ( fDebugStreamer ) delete fDebugStreamer;
152
153 }
154 //_____________________________________________________________________________
155 void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const
156 {
157   //
158   // Copy function
159   //
160
161   AliTRDCalibraVdriftLinearFit& target = (AliTRDCalibraVdriftLinearFit &) c;
162
163   // Copy only the histos
164   for (Int_t idet = 0; idet < 540; idet++){
165     if(fLinearFitterHistoArray.UncheckedAt(idet)){
166       TH2S *hped1 = (TH2S *)target.GetLinearFitterHisto(idet,kTRUE);
167       //hped1->SetDirectory(0);
168       hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet));
169     }
170   }
171   
172   TObject::Copy(c);
173
174 }
175 //_____________________________________________________________________________
176 Long64_t AliTRDCalibraVdriftLinearFit::Merge(const TCollection* list) 
177 {
178   // Merge list of objects (needed by PROOF)
179
180   if (!list)
181     return 0;
182   
183   if (list->IsEmpty())
184     return 1;
185   
186   TIterator* iter = list->MakeIterator();
187   TObject* obj = 0;
188   
189   // collection of generated histograms
190   Int_t count=0;
191   while((obj = iter->Next()) != 0) 
192     {
193       AliTRDCalibraVdriftLinearFit* entry = dynamic_cast<AliTRDCalibraVdriftLinearFit*>(obj);
194       if (entry == 0) continue; 
195       
196       // Copy only the histos
197       for (Int_t idet = 0; idet < 540; idet++){
198         if(entry->GetLinearFitterHisto(idet)){
199           TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
200           Double_t entriesa = hped1->GetEntries();
201           Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries();
202           if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet));
203         }
204       }
205       
206       count++;
207     }
208   
209   return count;
210 }
211 //_____________________________________________________________________
212 void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped)
213 {
214   //
215   // Add histo
216   //
217
218   fVersion++;
219
220   for (Int_t idet = 0; idet < 540; idet++){
221     const TH2S         *hped        = (TH2S*)ped->GetLinearFitterHistoNoForce(idet);
222     //printf("idet %d\n",idet);
223     if ( hped != 0x0 ){
224       //printf("add\n");
225       TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
226       Double_t entriesa = hped1->GetEntries();
227       Double_t entriesb = hped->GetEntries();
228       if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
229     }
230   }
231 }
232 //______________________________________________________________________________________
233 TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
234 {
235     //
236     // return pointer to TH2F histo 
237     // if force is true create a new histo if it doesn't exist allready
238     //
239     if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
240         return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
241
242     return GetLinearFitterHistoForce(detector);
243
244 }
245 //______________________________________________________________________________________
246 TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHistoForce(Int_t detector)
247 {
248   //
249   // return pointer to TH2F histo 
250   // if NULL create a new histo if it doesn't exist allready
251   //
252   if (fLinearFitterHistoArray.UncheckedAt(detector))
253     return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
254   
255   // if we are forced and TLinearFitter doesn't yes exist create it
256   
257   // new TH2F
258   TString name("LFDV");
259   name += detector;
260   name += "version";
261   name +=  fVersion;
262   
263   TH2S *lfdv = new TH2S((const Char_t *)name,(const Char_t *) name
264                         ,36,-0.9,0.9,48
265                         ,-1.2,1.2);
266   lfdv->SetXTitle("tan(phi_{track})");
267   lfdv->SetYTitle("dy/dt");
268   lfdv->SetZTitle("Number of clusters");
269   lfdv->SetStats(0);
270   lfdv->SetDirectory(0);
271   
272   fLinearFitterHistoArray.AddAt(lfdv,detector);
273   return lfdv;
274 }
275 //______________________________________________________________________________________
276 Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
277 {
278     //
279     // return param for this detector
280     //
281   if ( fLinearFitterPArray.UncheckedAt(detector) ){
282     const TVectorD     *vectorP     = (TVectorD*)fLinearFitterPArray.UncheckedAt(detector);
283     if(!param) param = new TVectorD(2);
284     for(Int_t k = 0; k < 2; k++){
285       (*param)[k] = (*vectorP)[k];
286     }
287     return kTRUE;
288   }
289   else return kFALSE;
290
291 }
292 //______________________________________________________________________________________
293 Bool_t AliTRDCalibraVdriftLinearFit::GetError(Int_t detector, TVectorD *error)
294 {
295     //
296     // return error for this detector 
297     //
298   if ( fLinearFitterEArray.UncheckedAt(detector) ){
299     const TVectorD     *vectorE     = (TVectorD*)fLinearFitterEArray.UncheckedAt(detector);
300     if(!error) error = new TVectorD(3);
301     for(Int_t k = 0; k < 3; k++){
302       (*error)[k] = (*vectorE)[k];
303     }
304     return kTRUE;
305   }
306   else return kFALSE;
307
308 }
309 //______________________________________________________________________________________
310 void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
311 {
312     //
313     // Fill the 2D histos for debugging
314     //
315   
316   TH2S *h = ((TH2S *) GetLinearFitterHisto(detector,kTRUE));
317   Double_t nbentries = h->GetEntries();
318   if(nbentries < 5*32767) h->Fill(tnp,pars1);
319
320 }
321 //____________Functions fit Online CH2d________________________________________
322 void AliTRDCalibraVdriftLinearFit::FillPEArray()
323 {
324   //
325   // Fill fLinearFitterPArray and fLinearFitterEArray from inside
326   //
327
328   
329   Int_t *arrayI = new Int_t[540];
330   for(Int_t k = 0; k< 540; k++){
331     arrayI[k] = 0; 
332   }
333
334   // Loop over histos 
335   for(Int_t cb = 0; cb < 540; cb++){
336     const TH2S *linearfitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
337     //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);    
338
339     if ( linearfitterhisto != 0 ){
340       
341       // Fill a linearfitter
342       TAxis *xaxis = linearfitterhisto->GetXaxis();
343       TAxis *yaxis = linearfitterhisto->GetYaxis();
344       TLinearFitter linearfitter = TLinearFitter(2,"pol1");
345       //printf("test\n");
346       Double_t integral = linearfitterhisto->Integral();
347       //printf("Integral is %f\n",integral);
348       Bool_t securitybreaking = kFALSE;
349       if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
350       for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
351         for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
352           if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
353             Double_t x = xaxis->GetBinCenter(ibinx+1);
354             Double_t y = yaxis->GetBinCenter(ibiny+1);
355             
356             for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
357               if(!securitybreaking){
358                 linearfitter.AddPoint(&x,y);
359                 arrayI[cb]++;
360               }
361               else {
362                 if(arrayI[cb]< 1198){
363                   linearfitter.AddPoint(&x,y);
364                   arrayI[cb]++; 
365                 }
366               }
367             }
368             
369           }
370         }
371       }
372       
373       //printf("Find %d entries for the detector %d\n",arrayI[cb],cb);
374
375       // Eval the linear fitter
376       if(arrayI[cb]>10){
377         TVectorD  *par  = new TVectorD(2);
378         TVectorD   pare = TVectorD(2);
379         TVectorD  *parE = new TVectorD(3);
380         //printf("Fit\n");
381         if(((fRobustFit) && (linearfitter.EvalRobust(0.8)==0)) || ((!fRobustFit) && (linearfitter.Eval()==0))) {
382           //if((linearfitter.Eval()==0)) {
383           //printf("Take the param\n");
384           linearfitter.GetParameters(*par);
385           //printf("Done\n");
386           //linearfitter.GetErrors(pare);
387           //Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter.GetChisquare())/arrayI[cb]);
388           //(*parE)[0] = pare[0]*ppointError;
389           //(*parE)[1] = pare[1]*ppointError;
390           
391           (*parE)[0] = 0.0;
392           (*parE)[1] = 0.0;
393           (*parE)[2] = (Double_t) arrayI[cb];
394           fLinearFitterPArray.AddAt(par,cb);
395           fLinearFitterEArray.AddAt(parE,cb);
396           
397           //par->Print();
398           //parE->Print();
399         }
400         //printf("Finish\n");
401       }
402       
403       //delete linearfitterhisto;
404       
405     }// if something
406
407   }
408
409   delete [] arrayI;
410    
411 }
412
413 //____________Functions fit Online CH2d________________________________________
414 void AliTRDCalibraVdriftLinearFit::FillPEArray2()
415 {
416   //
417   // Fill fFitterPArray and fFitterEArray from inside
418   //
419
420   // Loop over histos 
421   TF1 *f1 = new TF1("f1","[0]+[1]*x",-1,1);
422       
423   for(Int_t cb = 0; cb < 540; cb++){
424     //const TH2S *fitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
425     TH2S *fitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
426     //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) fitterhisto);    
427   
428     if ( fitterhisto != 0 ){
429       
430       Int_t nEntries=0;
431       TGraphErrors *gg=DrawMS(fitterhisto,nEntries);
432       if (!gg) continue;
433       // Number of points of the TGraphErrors
434       //printf("Number of points %d for detector %d\n",gg->GetN(),cb);
435       if(gg->GetN() <  fMinNpointsFit) {
436         if(gg) delete gg;
437         continue;       
438       }
439       //printf("det %d, number of points %d, nEntries %d\n",cb,gg->GetN(),nEntries);
440       gg->Fit(f1,"Q0");
441       
442       TVectorD  *par  = new TVectorD(2);
443       TVectorD  *parE = new TVectorD(3);
444       (*parE)[0] = 0.0;
445       (*parE)[1] = 0.0;
446       (*parE)[2] = (Double_t) nEntries;
447       (*par)[0] = f1->GetParameter(0);
448       (*par)[1] = f1->GetParameter(1);
449       fLinearFitterPArray.AddAt(par,cb);
450       fLinearFitterEArray.AddAt(parE,cb);
451       //printf("Value %f for detector %d\n",(*par)[1],cb);
452             
453       if(fDebugLevel==0) {
454         if(gg) delete gg;
455       }
456       else {
457         if(cb==fSeeDetector) {
458           gStyle->SetPalette(1);
459           gStyle->SetOptStat(1111);
460           gStyle->SetPadBorderMode(0);
461           gStyle->SetCanvasColor(10);
462           gStyle->SetPadLeftMargin(0.13);
463           gStyle->SetPadRightMargin(0.01);
464           TCanvas *stat = new TCanvas("stat","",50,50,600,800);
465           stat->cd(1);
466           fitterhisto->Draw("colz");
467           gg->Draw("P");
468           f1->Draw("same");
469           break;
470         }
471       } 
472     }
473   }
474   if(fDebugLevel==0) delete f1;
475 }
476
477 //_________Helper function__________________________________________________
478 TGraphErrors* AliTRDCalibraVdriftLinearFit::DrawMS(const TH2 *const h2, Int_t &nEntries)
479 {
480   //
481   // Debug function
482   //
483
484   TF1 fg("fg", "gaus", -10., 30.);
485   TGraphErrors *gp = new TGraphErrors();
486
487   TAxis *ax(h2->GetXaxis());
488   TAxis *ay(h2->GetYaxis());
489   TH1D *h1(NULL);
490   for(Int_t ipt(0), jpt(1), ig(0); ipt<ax->GetNbins(); ipt++, jpt++){
491     h1 = h2->ProjectionY("py", jpt, jpt);
492     fg.SetParameter(1, h1->GetMean());
493     //Float_t x(ax->GetBinCenter(jpt)); 
494     Int_t n=Int_t(h1->Integral(1, h1->GetNbinsX()));
495     nEntries+=n;
496     if(n < 15){
497       //Warning("drawMS()", Form("reject x[%d]=%f on n=%d", jpt, x, n));
498       continue;
499     }
500     h1->Fit(&fg, "WWQ0");
501     if(fg.GetNDF()<2){
502       //Warning("drawMS()", Form("reject x[%d]=%f on NDF=%d", jpt, x, fg.GetNDF()));
503       continue;
504     }
505     if(((fg.GetParameter(1)+fg.GetParameter(2)/2)>ay->GetXmax()) || ((fg.GetParameter(1)-fg.GetParameter(2)/2)<ay->GetXmin()) || (TMath::Abs(fg.GetParameter(0))< 0.00001)) continue;
506     gp->SetPoint(ig, ax->GetBinCenter(jpt), fg.GetParameter(1));
507     gp->SetPointError(ig, 0, TMath::Sqrt(pow(fg.GetParError(1),2) + (1/pow(fg.GetParameter(0),2))));
508     ig++;
509   }
510   delete h1;
511   return gp;
512 }