First version of some classes for comparison and performance studies (Marian)
[u/mrichter/AliRoot.git] / PWG1 / AliTreeDraw.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
17 ///////////////////////////////////////////////////////////////////////////
18 /*
19
20 Origin: marian.ivanov@cern.ch
21 Frequenlty used function for visualization 
22 marian.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"
52
53 //ALIROOT includes
54 #include "AliTrackPointArray.h"
55 #include "AliTreeDraw.h" 
56
57 #endif
58
59 //
60 //     Class for visualization and some statistacal analysis using tree
61 //     To be used in comparisons
62 //                and calib viewers based on tree    
63
64
65 ClassImp(AliTreeDraw)
66
67
68 AliTreeDraw::AliTreeDraw():
69   fTree(0),
70   fRes(0),
71   fMean(0),
72   fPoints(0){
73   //
74   // default constructor
75   //
76 }
77
78 void  AliTreeDraw::ClearHisto(){
79   //
80   //
81   delete fRes; 
82   delete fMean;
83   fRes=0;
84   fMean=0;
85 }
86
87
88
89 TH1F * AliTreeDraw::DrawXY(const char * chx, const char *chy, const char* selection, 
90                 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes)
91 {
92   //
93   Double_t* bins = CreateLogBins(nbins, minx, maxx);
94   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy);
95   char cut[1000];
96   sprintf(cut,"%s&&%s",selection,quality);
97   char expression[1000];
98   sprintf(expression,"%s:%s>>hRes2",chy,chx);
99   fTree->Draw(expression, cut, "groff");
100   TH1F* hMean=0;
101   TH1F* hRes = CreateResHisto(hRes2, &hMean);
102   AliLabelAxes(hRes, chx, chy);
103   //
104   delete hRes2;
105   delete[] bins;
106   ClearHisto();
107   fRes  = hRes;
108   fMean = hMean;
109   return hRes;
110 }
111
112
113
114 TH1F * AliTreeDraw::DrawLogXY(const char * chx, const char *chy, const char* selection, 
115                                     const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes)
116 {
117   //
118   // 
119   //
120   Double_t* bins = CreateLogBins(nbins, minx, maxx);
121   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy);
122   char cut[1000];
123   sprintf(cut,"%s&&%s",selection,quality);
124   char expression[1000];
125   sprintf(expression,"%s:%s>>hRes2",chy,chx);
126   fTree->Draw(expression, cut, "groff");
127   TH1F* hMean=0;  
128   TH1F* hRes = CreateResHisto(hRes2, &hMean);
129   AliLabelAxes(hRes, chx, chy);
130   //
131   delete hRes2;
132   delete[] bins;
133   ClearHisto();
134   fRes  = hRes;
135   fMean = hMean;
136   return hRes;
137 }
138
139 ///////////////////////////////////////////////////////////////////////////////////
140 ///////////////////////////////////////////////////////////////////////////////////
141 TH1F * AliTreeDraw::Eff(const char *variable, const char* selection, const char * quality, 
142                               Int_t nbins, Float_t min, Float_t max)
143 {
144   //
145   //
146   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max);
147   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max);
148   char inputGen[1000];  
149   sprintf(inputGen,"%s>>hGen", variable);
150   fTree->Draw(inputGen, selection, "groff");
151   char selectionRec[256];
152   sprintf(selectionRec, "%s && %s", selection, quality);
153   char inputRec[1000];  
154   sprintf(inputRec,"%s>>hRec", variable);
155   fTree->Draw(inputRec, selectionRec, "groff");
156   //
157   TH1F* hEff = CreateEffHisto(hGen, hRec);
158   AliLabelAxes(hEff, variable, "#epsilon [%]");
159   fRes = hEff;
160   delete hRec;
161   delete hGen;
162   return hEff;
163 }
164
165
166
167 ///////////////////////////////////////////////////////////////////////////////////
168 ///////////////////////////////////////////////////////////////////////////////////
169 TH1F * AliTreeDraw::EffLog(const char *variable, const char* selection, const char * quality, 
170                               Int_t nbins, Float_t min, Float_t max)
171 {
172   //
173   //
174   Double_t* bins = CreateLogBins(nbins, min, max);
175   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, bins);
176   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, bins);
177   char inputGen[1000];  
178   sprintf(inputGen,"%s>>hGen", variable);
179   fTree->Draw(inputGen, selection, "groff");
180   char selectionRec[256];
181   sprintf(selectionRec, "%s && %s", selection, quality);
182   char inputRec[1000];  
183   sprintf(inputRec,"%s>>hRec", variable);
184   fTree->Draw(inputRec, selectionRec, "groff");
185   //
186   TH1F* hEff = CreateEffHisto(hGen, hRec);
187   AliLabelAxes(hEff, variable, "#epsilon [%]");
188   fRes = hEff;
189   delete hRec;
190   delete hGen;
191   delete[] bins;
192   return hEff;
193 }
194
195
196 ///////////////////////////////////////////////////////////////////////////////////
197 ///////////////////////////////////////////////////////////////////////////////////
198
199 Double_t* AliTreeDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
200 {
201   Double_t* bins = new Double_t[nBins+1];
202   bins[0] = xMin;
203   Double_t factor = pow(xMax/xMin, 1./nBins);
204   for (Int_t i = 1; i <= nBins; i++)
205     bins[i] = factor * bins[i-1];
206   return bins;
207 }
208
209
210
211
212 void AliTreeDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
213 {
214   //
215   histo->GetXaxis()->SetTitle(xAxisTitle);
216   histo->GetYaxis()->SetTitle(yAxisTitle);
217 }
218
219
220 TH1F* AliTreeDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
221 {
222   //
223   Int_t nBins = hGen->GetNbinsX();
224   TH1F* hEff = (TH1F*) hGen->Clone("hEff");
225   hEff->SetTitle("");
226   hEff->SetStats(kFALSE);
227   hEff->SetMinimum(0.);
228   hEff->SetMaximum(110.);
229   //
230   for (Int_t iBin = 0; iBin <= nBins; iBin++) {
231     Double_t nGen = hGen->GetBinContent(iBin);
232     Double_t nRec = hRec->GetBinContent(iBin);
233     if (nGen > 0) {
234       Double_t eff = nRec/nGen;
235       hEff->SetBinContent(iBin, 100. * eff);
236       Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);      
237       //      if (error == 0) error = sqrt(0.1/nGen);
238       //
239       if (error == 0) error = 0.0001;
240       hEff->SetBinError(iBin, 100. * error);
241     } else {
242       hEff->SetBinContent(iBin, 100. * 0.5);
243       hEff->SetBinError(iBin, 100. * 0.5);
244     }
245   }
246   return hEff;
247 }
248
249
250
251 TH1F* AliTreeDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean,  Bool_t drawBinFits, 
252                      Bool_t overflowBinFits)
253 {
254   TVirtualPad* currentPad = gPad;
255   TAxis* axis = hRes2->GetXaxis();
256   Int_t nBins = axis->GetNbins();
257   TH1F* hRes, *hMean;
258   if (axis->GetXbins()->GetSize()){
259     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
260     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
261   }
262   else{
263     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
264     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
265
266   }
267   hRes->SetStats(false);
268   hRes->SetOption("E");
269   hRes->SetMinimum(0.);
270   //
271   hMean->SetStats(false);
272   hMean->SetOption("E");
273  
274   // create the fit function
275   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
276   
277   fitFunc->SetLineWidth(2);
278   fitFunc->SetFillStyle(0);
279   // create canvas for fits
280   TCanvas* canBinFits = NULL;
281   Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
282   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
283   Int_t ny = (nPads-1) / nx + 1;
284   if (drawBinFits) {
285     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
286     if (canBinFits) delete canBinFits;
287     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
288     canBinFits->Divide(nx, ny);
289   }
290
291   // loop over x bins and fit projection
292   Int_t dBin = ((overflowBinFits) ? 1 : 0);
293   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
294     if (drawBinFits) canBinFits->cd(bin + dBin);
295     TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
296     //    
297     if (hBin->GetEntries() > 5) {
298       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
299       hBin->Fit(fitFunc,"s");
300       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
301
302       if (sigma > 0.){
303         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
304         hMean->SetBinContent(bin, fitFunc->GetParameter(1));    
305       }
306       else{
307         hRes->SetBinContent(bin, 0.);
308         hMean->SetBinContent(bin,0);
309       }
310       hRes->SetBinError(bin, fitFunc->GetParError(2));
311       hMean->SetBinError(bin, fitFunc->GetParError(1));
312       
313       //
314       //
315
316     } else {
317       hRes->SetBinContent(bin, 0.);
318       hRes->SetBinError(bin, 0.);
319       hMean->SetBinContent(bin, 0.);
320       hMean->SetBinError(bin, 0.);
321     }
322     
323
324     if (drawBinFits) {
325       char name[256];
326       if (bin == 0) {
327         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
328       } else if (bin == nBins+1) {
329         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
330       } else {
331         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
332                 axis->GetTitle(), axis->GetBinUpEdge(bin));
333       }
334       canBinFits->cd(bin + dBin);
335       hBin->SetTitle(name);
336       hBin->SetStats(kTRUE);
337       hBin->DrawCopy("E");
338       canBinFits->Update();
339       canBinFits->Modified();
340       canBinFits->Update();
341     }
342     
343     delete hBin;
344   }
345
346   delete fitFunc;
347   currentPad->cd();
348   *phMean = hMean;
349   return hRes;
350 }
351
352
353
354
355 void   AliTreeDraw::GetPoints3D(const char * label, const char * chpoints, 
356                                 const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
357   //
358   // load selected points from tree
359   //
360    if (!fPoints) fPoints = new TObjArray;
361    if (tpoints->GetIndex()==0) tpoints->BuildIndex("fLabel","Label");
362    TBranch * br = tpoints->GetBranch(chpoints);
363    if (!br) return;
364    AliTrackPointArray * points = new AliTrackPointArray;
365    br->SetAddress(&points);
366
367    Int_t npoints = fTree->Draw(label,selection);
368    Float_t xyz[30000];
369    rmin*=rmin;
370    for (Int_t i=0;i<npoints;i++){
371      Int_t index = (Int_t)fTree->GetV1()[i];
372      tpoints->GetEntryWithIndex(index,index);
373      if (points->GetNPoints()<2) continue;
374      Int_t goodpoints=0;
375      for (Int_t i=0;i<points->GetNPoints(); i++){
376        xyz[goodpoints*3]   = points->GetX()[i];
377        xyz[goodpoints*3+1] = points->GetY()[i];
378        xyz[goodpoints*3+2] = points->GetZ()[i];
379        if ( points->GetX()[i]*points->GetX()[i]+points->GetY()[i]*points->GetY()[i]>rmin) goodpoints++;
380      } 
381      TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz); 
382      marker->SetMarkerColor(color);
383      marker->SetMarkerStyle(1);
384      fPoints->AddLast(marker);
385    }
386 }
387