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