1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////
20 Origin: marian.ivanov@cern.ch
21 Frequenlty used function for visualization
25 #if !defined(__CINT__) || defined(__MAKECINT__)
36 #include "TBenchmark.h"
37 #include "TStopwatch.h"
38 #include "TParticle.h"
49 #include "TPolyLine3D.h"
50 #include "TPolyMarker3D.h"
51 #include "TObjString.h"
55 #include "AliTrackPointArray.h"
56 #include "AliTreeDraw.h"
61 // Class for visualization and some statistacal analysis using tree
62 // To be used in comparisons
63 // and calib viewers based on tree
69 AliTreeDraw::AliTreeDraw():
75 // default constructor
79 void AliTreeDraw::ClearHisto(){
90 TH1F * 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)
94 Double_t* bins = CreateLogBins(nbins, minx, maxx);
95 TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy);
97 snprintf(cut,1000,"%s&&%s",selection,quality);
98 char expression[1000];
99 snprintf(expression,1000,"%s:%s>>hRes2",chy,chx);
100 fTree->Draw(expression, cut, "groff");
102 TH1F* hRes = CreateResHisto(hRes2, &hMean);
103 AliLabelAxes(hRes, chx, chy);
115 TH1F * 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)
121 Double_t* bins = CreateLogBins(nbins, minx, maxx);
122 TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy);
124 snprintf(cut,1000,"%s&&%s",selection,quality);
125 char expression[1000];
126 snprintf(expression,1000,"%s:%s>>hRes2",chy,chx);
127 fTree->Draw(expression, cut, "groff");
129 TH1F* hRes = CreateResHisto(hRes2, &hMean);
130 AliLabelAxes(hRes, chx, chy);
140 ///////////////////////////////////////////////////////////////////////////////////
141 ///////////////////////////////////////////////////////////////////////////////////
142 TH1F * AliTreeDraw::Eff(const char *variable, const char* selection, const char * quality,
143 Int_t nbins, Float_t min, Float_t max)
147 TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max);
148 TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max);
150 snprintf(inputGen,1000,"%s>>hGen", variable);
151 fTree->Draw(inputGen, selection, "groff");
152 char selectionRec[256];
153 snprintf(selectionRec,256, "%s && %s", selection, quality);
155 snprintf(inputRec,1000,"%s>>hRec", variable);
156 fTree->Draw(inputRec, selectionRec, "groff");
158 TH1F* hEff = CreateEffHisto(hGen, hRec);
159 AliLabelAxes(hEff, variable, "#epsilon [%]");
168 ///////////////////////////////////////////////////////////////////////////////////
169 ///////////////////////////////////////////////////////////////////////////////////
170 TH1F * AliTreeDraw::EffLog(const char *variable, const char* selection, const char * quality,
171 Int_t nbins, Float_t min, Float_t max)
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);
179 snprintf(inputGen,1000,"%s>>hGen", variable);
180 fTree->Draw(inputGen, selection, "groff");
181 char selectionRec[256];
182 snprintf(selectionRec,256, "%s && %s", selection, quality);
184 snprintf(inputRec,1000,"%s>>hRec", variable);
185 fTree->Draw(inputRec, selectionRec, "groff");
187 TH1F* hEff = CreateEffHisto(hGen, hRec);
188 AliLabelAxes(hEff, variable, "#epsilon [%]");
197 ///////////////////////////////////////////////////////////////////////////////////
198 ///////////////////////////////////////////////////////////////////////////////////
200 Double_t* AliTreeDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
202 Double_t* bins = new Double_t[nBins+1];
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];
213 void AliTreeDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
216 histo->GetXaxis()->SetTitle(xAxisTitle);
217 histo->GetYaxis()->SetTitle(yAxisTitle);
221 TH1F* AliTreeDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
224 Int_t nBins = hGen->GetNbinsX();
225 TH1F* hEff = (TH1F*) hGen->Clone("hEff");
227 hEff->SetStats(kFALSE);
228 hEff->SetMinimum(0.);
229 hEff->SetMaximum(110.);
231 for (Int_t iBin = 0; iBin <= nBins; iBin++) {
232 Double_t nGen = hGen->GetBinContent(iBin);
233 Double_t nRec = hRec->GetBinContent(iBin);
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);
240 if (error == 0) error = 0.0001;
241 hEff->SetBinError(iBin, 100. * error);
243 hEff->SetBinContent(iBin, 100. * 0.5);
244 hEff->SetBinError(iBin, 100. * 0.5);
252 TH1F* AliTreeDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean, Bool_t drawBinFits,
253 Bool_t overflowBinFits)
255 TVirtualPad* currentPad = gPad;
256 TAxis* axis = hRes2->GetXaxis();
257 Int_t nBins = axis->GetNbins();
259 if (axis->GetXbins()->GetSize()){
260 hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
261 hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
264 hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
265 hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
268 hRes->SetStats(false);
269 hRes->SetOption("E");
270 hRes->SetMinimum(0.);
272 hMean->SetStats(false);
273 hMean->SetOption("E");
275 // create the fit function
276 TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
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;
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);
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);
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));
304 hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
305 hMean->SetBinContent(bin, fitFunc->GetParameter(1));
308 hRes->SetBinContent(bin, 0.);
309 hMean->SetBinContent(bin,0);
311 hRes->SetBinError(bin, fitFunc->GetParError(2));
312 hMean->SetBinError(bin, fitFunc->GetParError(1));
318 hRes->SetBinContent(bin, 0.);
319 hRes->SetBinError(bin, 0.);
320 hMean->SetBinContent(bin, 0.);
321 hMean->SetBinError(bin, 0.);
328 snprintf(name,256, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
329 } else if (bin == nBins+1) {
330 snprintf(name,256, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
332 snprintf(name,256, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
333 axis->GetTitle(), axis->GetBinUpEdge(bin));
335 canBinFits->cd(bin + dBin);
336 hBin->SetTitle(name);
337 hBin->SetStats(kTRUE);
339 canBinFits->Update();
340 canBinFits->Modified();
341 canBinFits->Update();
353 TH1F* AliTreeDraw::CreateResHistoI(TH2F* hRes2, TH1F **phMean, Int_t integ, Bool_t drawBinFits)
355 TVirtualPad* currentPad = gPad;
356 TAxis* axis = hRes2->GetXaxis();
357 Int_t nBins = axis->GetNbins();
358 //Bool_t overflowBinFits = kFALSE;
361 if (axis->GetXbins()->GetSize()){
362 hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
363 hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
366 hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
367 hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
370 hRes->SetStats(false);
371 hRes->SetOption("E");
372 hRes->SetMinimum(0.);
374 hMean->SetStats(false);
375 hMean->SetOption("E");
377 // create the fit function
378 TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
380 fitFunc->SetLineWidth(2);
381 fitFunc->SetFillStyle(0);
382 // create canvas for fits
383 TCanvas* canBinFits = NULL;
384 //Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
386 Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
387 Int_t ny = (nPads-1) / nx + 1;
389 canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
390 if (canBinFits) delete canBinFits;
391 canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
392 canBinFits->Divide(nx, ny);
395 // loop over x bins and fit projection
396 //Int_t dBin = ((overflowBinFits) ? 1 : 0);
398 for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
399 if (drawBinFits) canBinFits->cd(bin + dBin);
400 Int_t bin0=TMath::Max(bin-integ,0);
401 Int_t bin1=TMath::Min(bin+integ,nBins);
402 TH1D* hBin = hRes2->ProjectionY("hBin", bin0, bin1);
404 if (hBin->GetEntries() > 5) {
405 fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
406 hBin->Fit(fitFunc,"s");
407 Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
410 hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
411 hMean->SetBinContent(bin, fitFunc->GetParameter(1));
414 hRes->SetBinContent(bin, 0.);
415 hMean->SetBinContent(bin,0);
417 hRes->SetBinError(bin, fitFunc->GetParError(2));
418 hMean->SetBinError(bin, fitFunc->GetParError(1));
424 hRes->SetBinContent(bin, 0.);
425 hRes->SetBinError(bin, 0.);
426 hMean->SetBinContent(bin, 0.);
427 hMean->SetBinError(bin, 0.);
434 snprintf(name,256, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
435 } else if (bin == nBins+1) {
436 snprintf(name,256, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
438 snprintf(name,256, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
439 axis->GetTitle(), axis->GetBinUpEdge(bin));
441 canBinFits->cd(bin + dBin);
442 hBin->SetTitle(name);
443 hBin->SetStats(kTRUE);
445 canBinFits->Update();
446 canBinFits->Modified();
447 canBinFits->Update();
459 TH1F* AliTreeDraw::CreateResHistoII(TH2F* hRes2, TH1F **phMean, Int_t integ, Bool_t drawBinFits, Int_t cut)
461 TVirtualPad* currentPad = gPad;
462 TAxis* axis = hRes2->GetXaxis();
463 Int_t nBins = axis->GetNbins();
464 //Bool_t overflowBinFits = kFALSE;
466 if (axis->GetXbins()->GetSize()){
467 hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
468 hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
471 hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
472 hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
475 hRes->SetStats(false);
476 hRes->SetOption("E");
477 hRes->SetMinimum(0.);
479 hMean->SetStats(false);
480 hMean->SetOption("E");
482 // create the fit function
483 TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
484 fitFunc->SetParLimits(2,0.0001,100.);
486 fitFunc->SetLineWidth(2);
487 fitFunc->SetFillStyle(0);
488 // create canvas for fits
489 TCanvas* canBinFits = NULL;
490 //Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
492 Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
493 Int_t ny = (nPads-1) / nx + 1;
495 canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
496 if (canBinFits) delete canBinFits;
497 canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
498 canBinFits->Divide(nx, ny);
501 // loop over x bins and fit projection
502 //Int_t dBin = ((overflowBinFits) ? 1 : 0);
504 for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
505 if (drawBinFits) canBinFits->cd(bin + dBin);
506 Int_t bin0=TMath::Max(bin-integ,0);
507 Int_t bin1=TMath::Min(bin+integ,nBins);
508 TH1D* hBin = hRes2->ProjectionY("hBin", bin0, bin1);
510 if (hBin->GetEntries() > cut) {
511 fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
512 hBin->Fit(fitFunc,"s");
513 Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
516 hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
517 hMean->SetBinContent(bin, fitFunc->GetParameter(1));
520 hRes->SetBinContent(bin, 0.);
521 hMean->SetBinContent(bin,0);
523 hRes->SetBinError(bin, fitFunc->GetParError(2));
524 hMean->SetBinError(bin, fitFunc->GetParError(1));
530 hRes->SetBinContent(bin, 0.);
531 hRes->SetBinError(bin, 0.);
532 hMean->SetBinContent(bin, 0.);
533 hMean->SetBinError(bin, 0.);
540 snprintf(name,256, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
541 } else if (bin == nBins+1) {
542 snprintf(name,256, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
544 snprintf(name,256, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
545 axis->GetTitle(), axis->GetBinUpEdge(bin));
547 canBinFits->cd(bin + dBin);
548 hBin->SetTitle(name);
549 hBin->SetStats(kTRUE);
551 canBinFits->Update();
552 canBinFits->Modified();
553 canBinFits->Update();
568 void AliTreeDraw::GetPoints3D(const char * label, const char * chpoints,
569 const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
571 // load selected points from tree
573 if (!fPoints) fPoints = new TObjArray;
574 if (tpoints->GetIndex()==0) tpoints->BuildIndex("fLabel","Label");
575 TBranch * br = tpoints->GetBranch(chpoints);
577 AliTrackPointArray * points = new AliTrackPointArray;
578 br->SetAddress(&points);
580 Int_t npoints = fTree->Draw(label,selection);
583 for (Int_t ii=0;ii<npoints;ii++){
584 Int_t index = (Int_t)fTree->GetV1()[ii];
585 tpoints->GetEntryWithIndex(index,index);
586 if (points->GetNPoints()<2) continue;
588 for (Int_t i=0;i<points->GetNPoints(); i++){
589 xyz[goodpoints*3] = points->GetX()[i];
590 xyz[goodpoints*3+1] = points->GetY()[i];
591 xyz[goodpoints*3+2] = points->GetZ()[i];
592 if ( points->GetX()[i]*points->GetX()[i]+points->GetY()[i]*points->GetY()[i]>rmin) goodpoints++;
594 TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz);
595 marker->SetMarkerColor(color);
596 marker->SetMarkerStyle(1);
597 fPoints->AddLast(marker);
604 TString* AliTreeDraw::FitPlane(const char* drawCommand, const char* formula, const char* cuts, Double_t & chi2, TVectorD &fitParam, TMatrixD &covMatrix, Int_t start, Int_t stop){
606 // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
607 // returns chi2, fitParam and covMatrix
608 // returns TString with fitted formula
611 TString formulaStr(formula);
612 TString drawStr(drawCommand);
613 TString cutStr(cuts);
615 formulaStr.ReplaceAll("++", "~");
616 TObjArray* formulaTokens = formulaStr.Tokenize("~");
617 Int_t dim = formulaTokens->GetEntriesFast();
619 fitParam.ResizeTo(dim);
620 covMatrix.ResizeTo(dim,dim);
622 TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
623 fitter->StoreData(kTRUE);
624 fitter->ClearPoints();
626 Int_t entries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start);
627 if (entries == -1) return new TString("An ERROR has occured during fitting!");
628 Double_t **values = new Double_t*[dim+1] ;
630 for (Int_t i = 0; i < dim + 1; i++) {
634 for (Int_t i = 0; i < dim + 1; i++) {
636 if (i < dim) centries = fTree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
637 else centries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
639 if (entries != centries) {
640 for (Int_t j = 0; j < dim + 1; j++) {
641 if(values[j]) delete values[j];
644 return new TString("An ERROR has occured during fitting!");
647 values[i] = new Double_t[entries];
648 memcpy(values[i], fTree->GetV1(), entries*sizeof(Double_t));
652 // add points to the fitter
653 for (Int_t i = 0; i < entries; i++) {
655 for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
656 fitter->AddPoint(x, values[dim][i], 1);
660 fitter->GetParameters(fitParam);
661 fitter->GetCovarianceMatrix(covMatrix);
662 chi2 = fitter->GetChisquare();
664 TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula;
666 for (Int_t iparam = 0; iparam < dim; iparam++) {
667 returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
668 if (iparam < dim-1) returnFormula.Append("+");
670 returnFormula.Append(" )");
671 delete formulaTokens;
674 for (Int_t i = 0; i < dim + 1; i++) {
678 return preturnFormula;