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