remove obselete clean QA macros command
[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 "TPolyLine3D.h"
50 #include "TPolyMarker3D.h"
51 #include "TObjString.h"
52
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
66 ClassImp(AliTreeDraw)
67
68
69 AliTreeDraw::AliTreeDraw():
70   fTree(0),
71   fRes(0),
72   fMean(0),
73   fPoints(0){
74   //
75   // default constructor
76   //
77 }
78
79 void  AliTreeDraw::ClearHisto(){
80   //
81   //
82   delete fRes; 
83   delete fMean;
84   fRes=0;
85   fMean=0;
86 }
87
88
89
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)
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
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)
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 ///////////////////////////////////////////////////////////////////////////////////
142 TH1F * 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 ///////////////////////////////////////////////////////////////////////////////////
170 TH1F * 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
200 Double_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
213 void 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
221 TH1F* 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
252 TH1F* 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
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
353 TH1F* 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
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
456
457
458
459 void   AliTreeDraw::GetPoints3D(const char * label, const char * chpoints, 
460                                 const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
461   //
462   // load selected points from tree
463   //
464    if (!fPoints) fPoints = new TObjArray;
465    if (tpoints->GetIndex()==0) tpoints->BuildIndex("fLabel","Label");
466    TBranch * br = tpoints->GetBranch(chpoints);
467    if (!br) return;
468    AliTrackPointArray * points = new AliTrackPointArray;
469    br->SetAddress(&points);
470
471    Int_t npoints = fTree->Draw(label,selection);
472    Float_t xyz[30000];
473    rmin*=rmin;
474    for (Int_t ii=0;ii<npoints;ii++){
475      Int_t index = (Int_t)fTree->GetV1()[ii];
476      tpoints->GetEntryWithIndex(index,index);
477      if (points->GetNPoints()<2) continue;
478      Int_t goodpoints=0;
479      for (Int_t i=0;i<points->GetNPoints(); i++){
480        xyz[goodpoints*3]   = points->GetX()[i];
481        xyz[goodpoints*3+1] = points->GetY()[i];
482        xyz[goodpoints*3+2] = points->GetZ()[i];
483        if ( points->GetX()[i]*points->GetX()[i]+points->GetY()[i]*points->GetY()[i]>rmin) goodpoints++;
484      } 
485      TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz); 
486      marker->SetMarkerColor(color);
487      marker->SetMarkerStyle(1);
488      fPoints->AddLast(marker);
489    }
490 }
491
492
493
494
495 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){
496    //
497    // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
498    // returns chi2, fitParam and covMatrix
499    // returns TString with fitted formula
500    //
501     
502    TString formulaStr(formula); 
503    TString drawStr(drawCommand);
504    TString cutStr(cuts);
505       
506    formulaStr.ReplaceAll("++", "~");
507    TObjArray* formulaTokens = formulaStr.Tokenize("~"); 
508    Int_t dim = formulaTokens->GetEntriesFast();
509    
510    fitParam.ResizeTo(dim);
511    covMatrix.ResizeTo(dim,dim);
512    
513    TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
514    fitter->StoreData(kTRUE);   
515    fitter->ClearPoints();
516    
517    Int_t entries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
518    if (entries == -1) return new TString("An ERROR has occured during fitting!");
519    Double_t **values = new Double_t*[dim+1] ; 
520    
521    for (Int_t i = 0; i < dim + 1; i++){
522       Int_t centries = 0;
523       if (i < dim) centries = fTree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
524       else  centries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
525       
526       if (entries != centries) return new TString("An ERROR has occured during fitting!");
527       values[i] = new Double_t[entries];
528       memcpy(values[i],  fTree->GetV1(), entries*sizeof(Double_t)); 
529    }
530    
531    // add points to the fitter
532    for (Int_t i = 0; i < entries; i++){
533       Double_t x[1000];
534       for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
535       fitter->AddPoint(x, values[dim][i], 1);
536    }
537
538    fitter->Eval();
539    fitter->GetParameters(fitParam);
540    fitter->GetCovarianceMatrix(covMatrix);
541    chi2 = fitter->GetChisquare();
542    chi2 = chi2;
543    
544    TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula; 
545    
546    for (Int_t iparam = 0; iparam < dim; iparam++) {
547      returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
548      if (iparam < dim-1) returnFormula.Append("+");
549    }
550    returnFormula.Append(" )");
551    delete formulaTokens;
552    delete fitter;
553    delete[] values;
554    return preturnFormula;
555 }
556