]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/AliTreeDraw.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
modifications to TLists naming to allow merging
[u/mrichter/AliRoot.git] / PWG1 / AliTreeDraw.cxx
CommitLineData
c92725b7 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
17///////////////////////////////////////////////////////////////////////////
18/*
19
20Origin: marian.ivanov@cern.ch
21Frequenlty used function for visualization
22marian.ivanov@cern.ch
23*/
24
25#if !defined(__CINT__) || defined(__MAKECINT__)
26#include <stdio.h>
27#include <string.h>
28//ROOT includes
29#include "TROOT.h"
30#include "Rtypes.h"
31#include "TFile.h"
32#include "TTree.h"
33#include "TChain.h"
34#include "TCut.h"
35#include "TString.h"
36#include "TBenchmark.h"
37#include "TStopwatch.h"
38#include "TParticle.h"
39#include "TSystem.h"
40#include "TTimer.h"
41#include "TVector3.h"
42#include "TH1F.h"
43#include "TH2F.h"
44#include "TCanvas.h"
45#include "TPad.h"
46#include "TF1.h"
47#include "TView.h"
48#include "TView3D.h"
c92725b7 49#include "TPolyLine3D.h"
50#include "TPolyMarker3D.h"
a36eadd7 51#include "TObjString.h"
52
c92725b7 53
54//ALIROOT includes
55#include "AliTrackPointArray.h"
56#include "AliTreeDraw.h"
57
58#endif
59
60//
61// Class for visualization and some statistacal analysis using tree
62// To be used in comparisons
63// and calib viewers based on tree
64
65
66ClassImp(AliTreeDraw)
67
68
69AliTreeDraw::AliTreeDraw():
70 fTree(0),
71 fRes(0),
72 fMean(0),
73 fPoints(0){
74 //
75 // default constructor
76 //
77}
78
79void AliTreeDraw::ClearHisto(){
80 //
81 //
82 delete fRes;
83 delete fMean;
84 fRes=0;
85 fMean=0;
86}
87
88
89
90TH1F * AliTreeDraw::DrawXY(const char * chx, const char *chy, const char* selection,
91 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes)
92{
93 //
94 Double_t* bins = CreateLogBins(nbins, minx, maxx);
95 TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy);
96 char cut[1000];
97 sprintf(cut,"%s&&%s",selection,quality);
98 char expression[1000];
99 sprintf(expression,"%s:%s>>hRes2",chy,chx);
100 fTree->Draw(expression, cut, "groff");
101 TH1F* hMean=0;
102 TH1F* hRes = CreateResHisto(hRes2, &hMean);
103 AliLabelAxes(hRes, chx, chy);
104 //
105 delete hRes2;
106 delete[] bins;
107 ClearHisto();
108 fRes = hRes;
109 fMean = hMean;
110 return hRes;
111}
112
113
114
115TH1F * AliTreeDraw::DrawLogXY(const char * chx, const char *chy, const char* selection,
116 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes)
117{
118 //
119 //
120 //
121 Double_t* bins = CreateLogBins(nbins, minx, maxx);
122 TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy);
123 char cut[1000];
124 sprintf(cut,"%s&&%s",selection,quality);
125 char expression[1000];
126 sprintf(expression,"%s:%s>>hRes2",chy,chx);
127 fTree->Draw(expression, cut, "groff");
128 TH1F* hMean=0;
129 TH1F* hRes = CreateResHisto(hRes2, &hMean);
130 AliLabelAxes(hRes, chx, chy);
131 //
132 delete hRes2;
133 delete[] bins;
134 ClearHisto();
135 fRes = hRes;
136 fMean = hMean;
137 return hRes;
138}
139
140///////////////////////////////////////////////////////////////////////////////////
141///////////////////////////////////////////////////////////////////////////////////
142TH1F * AliTreeDraw::Eff(const char *variable, const char* selection, const char * quality,
143 Int_t nbins, Float_t min, Float_t max)
144{
145 //
146 //
147 TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max);
148 TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max);
149 char inputGen[1000];
150 sprintf(inputGen,"%s>>hGen", variable);
151 fTree->Draw(inputGen, selection, "groff");
152 char selectionRec[256];
153 sprintf(selectionRec, "%s && %s", selection, quality);
154 char inputRec[1000];
155 sprintf(inputRec,"%s>>hRec", variable);
156 fTree->Draw(inputRec, selectionRec, "groff");
157 //
158 TH1F* hEff = CreateEffHisto(hGen, hRec);
159 AliLabelAxes(hEff, variable, "#epsilon [%]");
160 fRes = hEff;
161 delete hRec;
162 delete hGen;
163 return hEff;
164}
165
166
167
168///////////////////////////////////////////////////////////////////////////////////
169///////////////////////////////////////////////////////////////////////////////////
170TH1F * AliTreeDraw::EffLog(const char *variable, const char* selection, const char * quality,
171 Int_t nbins, Float_t min, Float_t max)
172{
173 //
174 //
175 Double_t* bins = CreateLogBins(nbins, min, max);
176 TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, bins);
177 TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, bins);
178 char inputGen[1000];
179 sprintf(inputGen,"%s>>hGen", variable);
180 fTree->Draw(inputGen, selection, "groff");
181 char selectionRec[256];
182 sprintf(selectionRec, "%s && %s", selection, quality);
183 char inputRec[1000];
184 sprintf(inputRec,"%s>>hRec", variable);
185 fTree->Draw(inputRec, selectionRec, "groff");
186 //
187 TH1F* hEff = CreateEffHisto(hGen, hRec);
188 AliLabelAxes(hEff, variable, "#epsilon [%]");
189 fRes = hEff;
190 delete hRec;
191 delete hGen;
192 delete[] bins;
193 return hEff;
194}
195
196
197///////////////////////////////////////////////////////////////////////////////////
198///////////////////////////////////////////////////////////////////////////////////
199
200Double_t* AliTreeDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
201{
202 Double_t* bins = new Double_t[nBins+1];
203 bins[0] = xMin;
204 Double_t factor = pow(xMax/xMin, 1./nBins);
205 for (Int_t i = 1; i <= nBins; i++)
206 bins[i] = factor * bins[i-1];
207 return bins;
208}
209
210
211
212
213void AliTreeDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
214{
215 //
216 histo->GetXaxis()->SetTitle(xAxisTitle);
217 histo->GetYaxis()->SetTitle(yAxisTitle);
218}
219
220
221TH1F* AliTreeDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
222{
223 //
224 Int_t nBins = hGen->GetNbinsX();
225 TH1F* hEff = (TH1F*) hGen->Clone("hEff");
226 hEff->SetTitle("");
227 hEff->SetStats(kFALSE);
228 hEff->SetMinimum(0.);
229 hEff->SetMaximum(110.);
230 //
231 for (Int_t iBin = 0; iBin <= nBins; iBin++) {
232 Double_t nGen = hGen->GetBinContent(iBin);
233 Double_t nRec = hRec->GetBinContent(iBin);
234 if (nGen > 0) {
235 Double_t eff = nRec/nGen;
236 hEff->SetBinContent(iBin, 100. * eff);
237 Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);
238 // if (error == 0) error = sqrt(0.1/nGen);
239 //
240 if (error == 0) error = 0.0001;
241 hEff->SetBinError(iBin, 100. * error);
242 } else {
243 hEff->SetBinContent(iBin, 100. * 0.5);
244 hEff->SetBinError(iBin, 100. * 0.5);
245 }
246 }
247 return hEff;
248}
249
250
251
252TH1F* AliTreeDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean, Bool_t drawBinFits,
253 Bool_t overflowBinFits)
254{
255 TVirtualPad* currentPad = gPad;
256 TAxis* axis = hRes2->GetXaxis();
257 Int_t nBins = axis->GetNbins();
258 TH1F* hRes, *hMean;
259 if (axis->GetXbins()->GetSize()){
260 hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
261 hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
262 }
263 else{
264 hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
265 hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
266
267 }
268 hRes->SetStats(false);
269 hRes->SetOption("E");
270 hRes->SetMinimum(0.);
271 //
272 hMean->SetStats(false);
273 hMean->SetOption("E");
274
275 // create the fit function
276 TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
277
278 fitFunc->SetLineWidth(2);
279 fitFunc->SetFillStyle(0);
280 // create canvas for fits
281 TCanvas* canBinFits = NULL;
282 Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
283 Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
284 Int_t ny = (nPads-1) / nx + 1;
285 if (drawBinFits) {
286 canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
287 if (canBinFits) delete canBinFits;
288 canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
289 canBinFits->Divide(nx, ny);
290 }
291
292 // loop over x bins and fit projection
293 Int_t dBin = ((overflowBinFits) ? 1 : 0);
294 for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
295 if (drawBinFits) canBinFits->cd(bin + dBin);
296 TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
297 //
298 if (hBin->GetEntries() > 5) {
299 fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
300 hBin->Fit(fitFunc,"s");
301 Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
302
303 if (sigma > 0.){
304 hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
305 hMean->SetBinContent(bin, fitFunc->GetParameter(1));
306 }
307 else{
308 hRes->SetBinContent(bin, 0.);
309 hMean->SetBinContent(bin,0);
310 }
311 hRes->SetBinError(bin, fitFunc->GetParError(2));
312 hMean->SetBinError(bin, fitFunc->GetParError(1));
313
314 //
315 //
316
317 } else {
318 hRes->SetBinContent(bin, 0.);
319 hRes->SetBinError(bin, 0.);
320 hMean->SetBinContent(bin, 0.);
321 hMean->SetBinError(bin, 0.);
322 }
323
324
c1a02aa0 325 if (drawBinFits) {
326 char name[256];
327 if (bin == 0) {
328 sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
329 } else if (bin == nBins+1) {
330 sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
331 } else {
332 sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
333 axis->GetTitle(), axis->GetBinUpEdge(bin));
334 }
335 canBinFits->cd(bin + dBin);
336 hBin->SetTitle(name);
337 hBin->SetStats(kTRUE);
338 hBin->DrawCopy("E");
339 canBinFits->Update();
340 canBinFits->Modified();
341 canBinFits->Update();
342 }
343
344 delete hBin;
345 }
346
347 delete fitFunc;
348 currentPad->cd();
349 *phMean = hMean;
350 return hRes;
351}
352
353TH1F* AliTreeDraw::CreateResHistoI(TH2F* hRes2, TH1F **phMean, Int_t integ, Bool_t drawBinFits)
354{
355 TVirtualPad* currentPad = gPad;
356 TAxis* axis = hRes2->GetXaxis();
357 Int_t nBins = axis->GetNbins();
358 Bool_t overflowBinFits = kFALSE;
359 TH1F* hRes, *hMean;
360 if (axis->GetXbins()->GetSize()){
361 hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
362 hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
363 }
364 else{
365 hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
366 hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
367
368 }
369 hRes->SetStats(false);
370 hRes->SetOption("E");
371 hRes->SetMinimum(0.);
372 //
373 hMean->SetStats(false);
374 hMean->SetOption("E");
375
376 // create the fit function
377 TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
378
379 fitFunc->SetLineWidth(2);
380 fitFunc->SetFillStyle(0);
381 // create canvas for fits
382 TCanvas* canBinFits = NULL;
383 Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
384 Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
385 Int_t ny = (nPads-1) / nx + 1;
386 if (drawBinFits) {
387 canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
388 if (canBinFits) delete canBinFits;
389 canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
390 canBinFits->Divide(nx, ny);
391 }
392
393 // loop over x bins and fit projection
394 Int_t dBin = ((overflowBinFits) ? 1 : 0);
395 for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
396 if (drawBinFits) canBinFits->cd(bin + dBin);
397 Int_t bin0=TMath::Max(bin-integ,0);
398 Int_t bin1=TMath::Min(bin+integ,nBins);
399 TH1D* hBin = hRes2->ProjectionY("hBin", bin0, bin1);
400 //
401 if (hBin->GetEntries() > 5) {
402 fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
403 hBin->Fit(fitFunc,"s");
404 Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
405
406 if (sigma > 0.){
407 hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
408 hMean->SetBinContent(bin, fitFunc->GetParameter(1));
409 }
410 else{
411 hRes->SetBinContent(bin, 0.);
412 hMean->SetBinContent(bin,0);
413 }
414 hRes->SetBinError(bin, fitFunc->GetParError(2));
415 hMean->SetBinError(bin, fitFunc->GetParError(1));
416
417 //
418 //
419
420 } else {
421 hRes->SetBinContent(bin, 0.);
422 hRes->SetBinError(bin, 0.);
423 hMean->SetBinContent(bin, 0.);
424 hMean->SetBinError(bin, 0.);
425 }
426
427
3dd64180 428 if (drawBinFits) {
429 char name[256];
430 if (bin == 0) {
431 sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
432 } else if (bin == nBins+1) {
433 sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
434 } else {
435 sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
436 axis->GetTitle(), axis->GetBinUpEdge(bin));
437 }
438 canBinFits->cd(bin + dBin);
439 hBin->SetTitle(name);
440 hBin->SetStats(kTRUE);
441 hBin->DrawCopy("E");
442 canBinFits->Update();
443 canBinFits->Modified();
444 canBinFits->Update();
445 }
446
447 delete hBin;
448 }
449
450 delete fitFunc;
451 currentPad->cd();
452 *phMean = hMean;
453 return hRes;
454}
455
456TH1F* AliTreeDraw::CreateResHistoII(TH2F* hRes2, TH1F **phMean, Int_t integ, Bool_t drawBinFits, Int_t cut)
457{
458 TVirtualPad* currentPad = gPad;
459 TAxis* axis = hRes2->GetXaxis();
460 Int_t nBins = axis->GetNbins();
461 Bool_t overflowBinFits = kFALSE;
462 TH1F* hRes, *hMean;
463 if (axis->GetXbins()->GetSize()){
464 hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
465 hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
466 }
467 else{
468 hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
469 hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
470
471 }
472 hRes->SetStats(false);
473 hRes->SetOption("E");
474 hRes->SetMinimum(0.);
475 //
476 hMean->SetStats(false);
477 hMean->SetOption("E");
478
479 // create the fit function
480 TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
481
482 fitFunc->SetLineWidth(2);
483 fitFunc->SetFillStyle(0);
484 // create canvas for fits
485 TCanvas* canBinFits = NULL;
486 Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
487 Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
488 Int_t ny = (nPads-1) / nx + 1;
489 if (drawBinFits) {
490 canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
491 if (canBinFits) delete canBinFits;
492 canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
493 canBinFits->Divide(nx, ny);
494 }
495
496 // loop over x bins and fit projection
497 Int_t dBin = ((overflowBinFits) ? 1 : 0);
498 for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
499 if (drawBinFits) canBinFits->cd(bin + dBin);
500 Int_t bin0=TMath::Max(bin-integ,0);
501 Int_t bin1=TMath::Min(bin+integ,nBins);
502 TH1D* hBin = hRes2->ProjectionY("hBin", bin0, bin1);
503 //
504 if (hBin->GetEntries() > cut) {
505 fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
506 hBin->Fit(fitFunc,"s");
507 Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
508
509 if (sigma > 0.){
510 hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
511 hMean->SetBinContent(bin, fitFunc->GetParameter(1));
512 }
513 else{
514 hRes->SetBinContent(bin, 0.);
515 hMean->SetBinContent(bin,0);
516 }
517 hRes->SetBinError(bin, fitFunc->GetParError(2));
518 hMean->SetBinError(bin, fitFunc->GetParError(1));
519
520 //
521 //
522
523 } else {
524 hRes->SetBinContent(bin, 0.);
525 hRes->SetBinError(bin, 0.);
526 hMean->SetBinContent(bin, 0.);
527 hMean->SetBinError(bin, 0.);
528 }
529
530
c92725b7 531 if (drawBinFits) {
532 char name[256];
533 if (bin == 0) {
534 sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
535 } else if (bin == nBins+1) {
536 sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
537 } else {
538 sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
539 axis->GetTitle(), axis->GetBinUpEdge(bin));
540 }
541 canBinFits->cd(bin + dBin);
542 hBin->SetTitle(name);
543 hBin->SetStats(kTRUE);
544 hBin->DrawCopy("E");
545 canBinFits->Update();
546 canBinFits->Modified();
547 canBinFits->Update();
548 }
549
550 delete hBin;
551 }
552
553 delete fitFunc;
554 currentPad->cd();
555 *phMean = hMean;
556 return hRes;
557}
558
559
560
561
562void AliTreeDraw::GetPoints3D(const char * label, const char * chpoints,
563 const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
564 //
565 // load selected points from tree
566 //
567 if (!fPoints) fPoints = new TObjArray;
568 if (tpoints->GetIndex()==0) tpoints->BuildIndex("fLabel","Label");
569 TBranch * br = tpoints->GetBranch(chpoints);
570 if (!br) return;
571 AliTrackPointArray * points = new AliTrackPointArray;
572 br->SetAddress(&points);
573
574 Int_t npoints = fTree->Draw(label,selection);
575 Float_t xyz[30000];
576 rmin*=rmin;
3a85f1d5 577 for (Int_t ii=0;ii<npoints;ii++){
578 Int_t index = (Int_t)fTree->GetV1()[ii];
c92725b7 579 tpoints->GetEntryWithIndex(index,index);
580 if (points->GetNPoints()<2) continue;
581 Int_t goodpoints=0;
582 for (Int_t i=0;i<points->GetNPoints(); i++){
583 xyz[goodpoints*3] = points->GetX()[i];
584 xyz[goodpoints*3+1] = points->GetY()[i];
585 xyz[goodpoints*3+2] = points->GetZ()[i];
586 if ( points->GetX()[i]*points->GetX()[i]+points->GetY()[i]*points->GetY()[i]>rmin) goodpoints++;
587 }
588 TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz);
589 marker->SetMarkerColor(color);
590 marker->SetMarkerStyle(1);
591 fPoints->AddLast(marker);
592 }
593}
594
a36eadd7 595
596
597
598TString* AliTreeDraw::FitPlane(const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, TVectorD &fitParam, TMatrixD &covMatrix, Int_t start, Int_t stop){
599 //
600 // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
601 // returns chi2, fitParam and covMatrix
602 // returns TString with fitted formula
603 //
604
605 TString formulaStr(formula);
606 TString drawStr(drawCommand);
607 TString cutStr(cuts);
608
609 formulaStr.ReplaceAll("++", "~");
610 TObjArray* formulaTokens = formulaStr.Tokenize("~");
611 Int_t dim = formulaTokens->GetEntriesFast();
612
613 fitParam.ResizeTo(dim);
614 covMatrix.ResizeTo(dim,dim);
615
616 TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
617 fitter->StoreData(kTRUE);
618 fitter->ClearPoints();
619
620 Int_t entries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start);
621 if (entries == -1) return new TString("An ERROR has occured during fitting!");
622 Double_t **values = new Double_t*[dim+1] ;
623
624 for (Int_t i = 0; i < dim + 1; i++){
625 Int_t centries = 0;
626 if (i < dim) centries = fTree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
627 else centries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
628
629 if (entries != centries) return new TString("An ERROR has occured during fitting!");
630 values[i] = new Double_t[entries];
631 memcpy(values[i], fTree->GetV1(), entries*sizeof(Double_t));
632 }
633
634 // add points to the fitter
635 for (Int_t i = 0; i < entries; i++){
636 Double_t x[1000];
637 for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
638 fitter->AddPoint(x, values[dim][i], 1);
639 }
640
641 fitter->Eval();
642 fitter->GetParameters(fitParam);
643 fitter->GetCovarianceMatrix(covMatrix);
644 chi2 = fitter->GetChisquare();
645 chi2 = chi2;
646
647 TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
648
649 for (Int_t iparam = 0; iparam < dim; iparam++) {
650 returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
651 if (iparam < dim-1) returnFormula.Append("+");
652 }
653 returnFormula.Append(" )");
654 delete formulaTokens;
655 delete fitter;
656 delete[] values;
657 return preturnFormula;
658}
659