Fix Coverity
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraVdriftLinearFit.cxx
CommitLineData
3a0f6479 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
0bc7827a 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
3a0f6479 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>
661c7be6 37#include <TStyle.h>
38#include <TCanvas.h>
3a0f6479 39#include <TTreeStream.h>
661c7be6 40#include <TGraphErrors.h>
41#include <TDirectory.h>
42#include <TTreeStream.h>
43#include <TF1.h>
3a0f6479 44
45//header file
46#include "AliTRDCalibraVdriftLinearFit.h"
47
48ClassImp(AliTRDCalibraVdriftLinearFit) /*FOLD00*/
49
50//_____________________________________________________________________
51AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
52 TObject(),
53 fVersion(0),
b2277aa2 54 fNameCalibUsed(""),
3a0f6479 55 fLinearFitterHistoArray(540),
56 fLinearFitterPArray(540),
2a1a7b36 57 fLinearFitterEArray(540),
661c7be6 58 fRobustFit(kTRUE),
59 fDebugStreamer(0x0),
60 fDebugLevel(0),
61 fSeeDetector(0)
3a0f6479 62{
ba94744a 63 //
64 // default constructor
65 //
3a0f6479 66}
3a0f6479 67//_____________________________________________________________________
68AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /*FOLD00*/
69 TObject(ped),
70 fVersion(ped.fVersion),
b2277aa2 71 fNameCalibUsed(ped.fNameCalibUsed),
3a0f6479 72 fLinearFitterHistoArray(540),
73 fLinearFitterPArray(540),
2a1a7b36 74 fLinearFitterEArray(540),
661c7be6 75 fRobustFit(kTRUE),
76 fDebugStreamer(0x0),
77 fDebugLevel(0),
78 fSeeDetector(0)
3a0f6479 79{
80 //
81 // copy constructor
82 //
83 for (Int_t idet = 0; idet < 540; idet++){
ba94744a 84
3a0f6479 85 const TVectorD *vectorE = (TVectorD*)ped.fLinearFitterEArray.UncheckedAt(idet);
86 const TVectorD *vectorP = (TVectorD*)ped.fLinearFitterPArray.UncheckedAt(idet);
0ca1720e 87 const TH2S *hped = (TH2S*)ped.fLinearFitterHistoArray.UncheckedAt(idet);
3a0f6479 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 ){
0ca1720e 92 TH2S *hNew = (TH2S *)hped->Clone();
ba94744a 93 //hNew->SetDirectory(0);
3a0f6479 94 fLinearFitterHistoArray.AddAt(hNew,idet);
95 }
3a0f6479 96 }
97}
98//_____________________________________________________________________
d0569428 99AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /*FOLD00*/
100 TObject(),
101 fVersion(0),
b2277aa2 102 fNameCalibUsed(""),
d0569428 103 fLinearFitterHistoArray(540),
104 fLinearFitterPArray(540),
2a1a7b36 105 fLinearFitterEArray(540),
661c7be6 106 fRobustFit(kTRUE),
107 fDebugStreamer(0x0),
108 fDebugLevel(0),
109 fSeeDetector(0)
d0569428 110{
111 //
112 // constructor from a TObjArray
113 //
114 for (Int_t idet = 0; idet < 540; idet++){
0ca1720e 115 const TH2S *hped = (TH2S*)obja.UncheckedAt(idet);
d0569428 116 if ( hped != 0x0 ){
0ca1720e 117 TH2S *hNew = (TH2S *)hped->Clone();
ba94744a 118 //hNew->SetDirectory(0);
d0569428 119 fLinearFitterHistoArray.AddAt(hNew,idet);
120 }
121 }
122}
123//_____________________________________________________________________
3a0f6479 124AliTRDCalibraVdriftLinearFit& 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//_____________________________________________________________________
135AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /*FOLD00*/
136{
137 //
138 // destructor
139 //
f5cda287 140 fLinearFitterHistoArray.SetOwner();
141 fLinearFitterPArray.SetOwner();
142 fLinearFitterEArray.SetOwner();
143
1ca79a00 144 fLinearFitterHistoArray.Delete();
145 fLinearFitterPArray.Delete();
146 fLinearFitterEArray.Delete();
147
661c7be6 148 if ( fDebugStreamer ) delete fDebugStreamer;
149
3a0f6479 150}
64942b85 151//_____________________________________________________________________________
152void 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)){
0ca1720e 163 TH2S *hped1 = (TH2S *)target.GetLinearFitterHisto(idet,kTRUE);
ba94744a 164 //hped1->SetDirectory(0);
0ca1720e 165 hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet));
64942b85 166 }
64942b85 167 }
ba94744a 168
64942b85 169 TObject::Copy(c);
170
171}
172//_____________________________________________________________________________
6bbdc11a 173Long64_t AliTRDCalibraVdriftLinearFit::Merge(const TCollection* list)
64942b85 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)
ba94744a 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)){
0ca1720e 196 TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
ba94744a 197 Double_t entriesa = hped1->GetEntries();
0ca1720e 198 Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries();
ba94744a 199 if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet));
200 }
64942b85 201 }
202
ba94744a 203 count++;
64942b85 204 }
64942b85 205
206 return count;
207}
3a0f6479 208//_____________________________________________________________________
fe4ee353 209void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped)
3a0f6479 210{
211 //
212 // Add histo
213 //
214
215 fVersion++;
216
217 for (Int_t idet = 0; idet < 540; idet++){
fe4ee353 218 const TH2S *hped = (TH2S*)ped->GetLinearFitterHistoNoForce(idet);
3a0f6479 219 //printf("idet %d\n",idet);
220 if ( hped != 0x0 ){
221 //printf("add\n");
0ca1720e 222 TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
ba94744a 223 Double_t entriesa = hped1->GetEntries();
224 Double_t entriesb = hped->GetEntries();
225 if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
3a0f6479 226 }
227 }
228}
229//______________________________________________________________________________________
0ca1720e 230TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
3a0f6479 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) )
0ca1720e 237 return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
3a0f6479 238
fe4ee353 239 return GetLinearFitterHistoForce(detector);
240
241}
242//______________________________________________________________________________________
243TH2S* 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;
3a0f6479 271}
272//______________________________________________________________________________________
273Bool_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//______________________________________________________________________________________
290Bool_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//______________________________________________________________________________________
307void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
308{
309 //
310 // Fill the 2D histos for debugging
311 //
ba94744a 312
0ca1720e 313 TH2S *h = ((TH2S *) GetLinearFitterHisto(detector,kTRUE));
ba94744a 314 Double_t nbentries = h->GetEntries();
0ca1720e 315 if(nbentries < 5*32767) h->Fill(tnp,pars1);
3a0f6479 316
317}
318//____________Functions fit Online CH2d________________________________________
319void 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++){
0ca1720e 333 const TH2S *linearfitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
3a0f6479 334 //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);
335
0ca1720e 336 if ( linearfitterhisto != 0 ){
3a0f6479 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");
a5dcf618 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;
3a0f6479 347 for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
348 for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
8a606eb5 349 if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
3a0f6479 350 Double_t x = xaxis->GetBinCenter(ibinx+1);
351 Double_t y = yaxis->GetBinCenter(ibiny+1);
a5dcf618 352
3a0f6479 353 for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
a5dcf618 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 }
3a0f6479 364 }
a5dcf618 365
3a0f6479 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);
a5dcf618 377 //printf("Fit\n");
2a1a7b36 378 if(((fRobustFit) && (linearfitter.EvalRobust(0.8)==0)) || ((!fRobustFit) && (linearfitter.Eval()==0))) {
a5dcf618 379 //if((linearfitter.Eval()==0)) {
380 //printf("Take the param\n");
c1105918 381 linearfitter.GetParameters(*par);
a5dcf618 382 //printf("Done\n");
c1105918 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 }
a5dcf618 397 //printf("Finish\n");
ba94744a 398 }
c1105918 399
07f9b7f2 400 //delete linearfitterhisto;
ba94744a 401
3a0f6479 402 }// if something
fdc15553 403
3a0f6479 404 }
a987273c 405
406 delete [] arrayI;
3a0f6479 407
d1e1a064 408}
c1105918 409
661c7be6 410//____________Functions fit Online CH2d________________________________________
411void 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);
ce42378b 429 if (!gg) continue;
661c7be6 430 // Number of points of the TGraphErrors
431 if(gg->GetN() < 20) {
432 if(gg) delete gg;
433 continue;
434 }
435 //printf("det %d, number of points %d, nEntries %d\n",cb,gg->GetN(),nEntries);
436 gg->Fit(f1,"Q0");
437
438 TVectorD *par = new TVectorD(2);
439 TVectorD *parE = new TVectorD(3);
440 (*parE)[0] = 0.0;
441 (*parE)[1] = 0.0;
442 (*parE)[2] = (Double_t) nEntries;
443 (*par)[0] = f1->GetParameter(0);
444 (*par)[1] = f1->GetParameter(1);
445 fLinearFitterPArray.AddAt(par,cb);
446 fLinearFitterEArray.AddAt(parE,cb);
447
448 if(fDebugLevel==0) {
449 if(gg) delete gg;
450 }
451 else {
452 if(cb==fSeeDetector) {
453 gStyle->SetPalette(1);
454 gStyle->SetOptStat(1111);
455 gStyle->SetPadBorderMode(0);
456 gStyle->SetCanvasColor(10);
457 gStyle->SetPadLeftMargin(0.13);
458 gStyle->SetPadRightMargin(0.01);
459 TCanvas *stat = new TCanvas("stat","",50,50,600,800);
460 stat->cd(1);
461 fitterhisto->Draw("colz");
462 gg->Draw("P");
463 f1->Draw("same");
464 break;
465 }
466 }
467 }
468 }
469 if(fDebugLevel==0) delete f1;
470}
471
472//_________Helper function__________________________________________________
473TGraphErrors* AliTRDCalibraVdriftLinearFit::DrawMS(const TH2 *const h2, Int_t &nEntries)
474{
475 TF1 fg("fg", "gaus", -10., 30.);
476 TGraphErrors *gp = new TGraphErrors();
c1105918 477
661c7be6 478 TAxis *ax(h2->GetXaxis());
479 TAxis *ay(h2->GetYaxis());
480 TH1D *h1(NULL);
481 for(Int_t ipt(0), jpt(1), ig(0); ipt<ax->GetNbins(); ipt++, jpt++){
482 h1 = h2->ProjectionY("py", jpt, jpt);
483 fg.SetParameter(1, h1->GetMean());
484 //Float_t x(ax->GetBinCenter(jpt));
485 Int_t n=Int_t(h1->Integral(1, h1->GetNbinsX()));
486 nEntries+=n;
487 if(n < 15){
488 //Warning("drawMS()", Form("reject x[%d]=%f on n=%d", jpt, x, n));
489 continue;
490 }
491 h1->Fit(&fg, "WWQ0");
492 if(fg.GetNDF()<2){
493 //Warning("drawMS()", Form("reject x[%d]=%f on NDF=%d", jpt, x, fg.GetNDF()));
494 continue;
495 }
64c7fdf8 496 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;
661c7be6 497 gp->SetPoint(ig, ax->GetBinCenter(jpt), fg.GetParameter(1));
498 gp->SetPointError(ig, 0, TMath::Sqrt(pow(fg.GetParError(1),2) + (1/pow(fg.GetParameter(0),2))));
499 ig++;
500 }
501 delete h1;
502 return gp;
503}