[u/mrichter/AliRoot.git] / PWG1 / AliTreeDraw.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


///////////////////////////////////////////////////////////////////////////
/*

Origin: marian.ivanov@cern.ch
Frequenlty used function for visualization 
marian.ivanov@cern.ch
*/

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <stdio.h>
#include <string.h>
//ROOT includes
#include "TROOT.h"
#include "Rtypes.h"
#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TCut.h"
#include "TString.h"
#include "TBenchmark.h"
#include "TStopwatch.h"
#include "TParticle.h"
#include "TSystem.h"
#include "TTimer.h"
#include "TVector3.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TF1.h"
#include "TView.h"
#include "TView3D.h"
#include "TGeometry.h"
#include "TPolyLine3D.h"
#include "TPolyMarker3D.h"
#include "TObjString.h"


//ALIROOT includes
#include "AliTrackPointArray.h"
#include "AliTreeDraw.h" 

#endif

//
//     Class for visualization and some statistacal analysis using tree
//     To be used in comparisons
//                and calib viewers based on tree    


ClassImp(AliTreeDraw)


AliTreeDraw::AliTreeDraw():
  fTree(0),
  fRes(0),
  fMean(0),
  fPoints(0){
  //
  // default constructor
  //
}

void  AliTreeDraw::ClearHisto(){
  //
  //
  delete fRes; 
  delete fMean;
  fRes=0;
  fMean=0;
}


TH1F * AliTreeDraw::DrawXY(const char * chx, const char *chy, const char* selection, 
		const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes)
{
  //
  Double_t* bins = CreateLogBins(nbins, minx, maxx);
  TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy);
  char cut[1000];
  sprintf(cut,"%s&&%s",selection,quality);
  char expression[1000];
  sprintf(expression,"%s:%s>>hRes2",chy,chx);
  fTree->Draw(expression, cut, "groff");
  TH1F* hMean=0;
  TH1F* hRes = CreateResHisto(hRes2, &hMean);
  AliLabelAxes(hRes, chx, chy);
  //
  delete hRes2;
  delete[] bins;
  ClearHisto();
  fRes  = hRes;
  fMean = hMean;
  return hRes;
}


TH1F * AliTreeDraw::DrawLogXY(const char * chx, const char *chy, const char* selection, 
				    const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes)
{
  //
  // 
  //
  Double_t* bins = CreateLogBins(nbins, minx, maxx);
  TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy);
  char cut[1000];
  sprintf(cut,"%s&&%s",selection,quality);
  char expression[1000];
  sprintf(expression,"%s:%s>>hRes2",chy,chx);
  fTree->Draw(expression, cut, "groff");
  TH1F* hMean=0;  
  TH1F* hRes = CreateResHisto(hRes2, &hMean);
  AliLabelAxes(hRes, chx, chy);
  //
  delete hRes2;
  delete[] bins;
  ClearHisto();
  fRes  = hRes;
  fMean = hMean;
  return hRes;
}

///////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////
TH1F * AliTreeDraw::Eff(const char *variable, const char* selection, const char * quality, 
			      Int_t nbins, Float_t min, Float_t max)
{
  //
  //
  TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max);
  TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max);
  char inputGen[1000];  
  sprintf(inputGen,"%s>>hGen", variable);
  fTree->Draw(inputGen, selection, "groff");
  char selectionRec[256];
  sprintf(selectionRec, "%s && %s", selection, quality);
  char inputRec[1000];  
  sprintf(inputRec,"%s>>hRec", variable);
  fTree->Draw(inputRec, selectionRec, "groff");
  //
  TH1F* hEff = CreateEffHisto(hGen, hRec);
  AliLabelAxes(hEff, variable, "#epsilon [%]");
  fRes = hEff;
  delete hRec;
  delete hGen;
  return hEff;
}


///////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////
TH1F * AliTreeDraw::EffLog(const char *variable, const char* selection, const char * quality, 
			      Int_t nbins, Float_t min, Float_t max)
{
  //
  //
  Double_t* bins = CreateLogBins(nbins, min, max);
  TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, bins);
  TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, bins);
  char inputGen[1000];  
  sprintf(inputGen,"%s>>hGen", variable);
  fTree->Draw(inputGen, selection, "groff");
  char selectionRec[256];
  sprintf(selectionRec, "%s && %s", selection, quality);
  char inputRec[1000];  
  sprintf(inputRec,"%s>>hRec", variable);
  fTree->Draw(inputRec, selectionRec, "groff");
  //
  TH1F* hEff = CreateEffHisto(hGen, hRec);
  AliLabelAxes(hEff, variable, "#epsilon [%]");
  fRes = hEff;
  delete hRec;
  delete hGen;
  delete[] bins;
  return hEff;
}


///////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////

Double_t* AliTreeDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
{
  Double_t* bins = new Double_t[nBins+1];
  bins[0] = xMin;
  Double_t factor = pow(xMax/xMin, 1./nBins);
  for (Int_t i = 1; i <= nBins; i++)
    bins[i] = factor * bins[i-1];
  return bins;
}


void AliTreeDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
{
  //
  histo->GetXaxis()->SetTitle(xAxisTitle);
  histo->GetYaxis()->SetTitle(yAxisTitle);
}


TH1F* AliTreeDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
{
  //
  Int_t nBins = hGen->GetNbinsX();
  TH1F* hEff = (TH1F*) hGen->Clone("hEff");
  hEff->SetTitle("");
  hEff->SetStats(kFALSE);
  hEff->SetMinimum(0.);
  hEff->SetMaximum(110.);
  //
  for (Int_t iBin = 0; iBin <= nBins; iBin++) {
    Double_t nGen = hGen->GetBinContent(iBin);
    Double_t nRec = hRec->GetBinContent(iBin);
    if (nGen > 0) {
      Double_t eff = nRec/nGen;
      hEff->SetBinContent(iBin, 100. * eff);
      Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);      
      //      if (error == 0) error = sqrt(0.1/nGen);
      //
      if (error == 0) error = 0.0001;
      hEff->SetBinError(iBin, 100. * error);
    } else {
      hEff->SetBinContent(iBin, 100. * 0.5);
      hEff->SetBinError(iBin, 100. * 0.5);
    }
  }
  return hEff;
}


TH1F* AliTreeDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean,  Bool_t drawBinFits, 
		     Bool_t overflowBinFits)
{
  TVirtualPad* currentPad = gPad;
  TAxis* axis = hRes2->GetXaxis();
  Int_t nBins = axis->GetNbins();
  TH1F* hRes, *hMean;
  if (axis->GetXbins()->GetSize()){
    hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
    hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
  }
  else{
    hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
    hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());

  }
  hRes->SetStats(false);
  hRes->SetOption("E");
  hRes->SetMinimum(0.);
  //
  hMean->SetStats(false);
  hMean->SetOption("E");
 
  // create the fit function
  TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
  
  fitFunc->SetLineWidth(2);
  fitFunc->SetFillStyle(0);
  // create canvas for fits
  TCanvas* canBinFits = NULL;
  Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
  Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
  Int_t ny = (nPads-1) / nx + 1;
  if (drawBinFits) {
    canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
    if (canBinFits) delete canBinFits;
    canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
    canBinFits->Divide(nx, ny);
  }

  // loop over x bins and fit projection
  Int_t dBin = ((overflowBinFits) ? 1 : 0);
  for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
    if (drawBinFits) canBinFits->cd(bin + dBin);
    TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
    //    
    if (hBin->GetEntries() > 5) {
      fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
      hBin->Fit(fitFunc,"s");
      Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));

      if (sigma > 0.){
	hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
	hMean->SetBinContent(bin, fitFunc->GetParameter(1));	
      }
      else{
	hRes->SetBinContent(bin, 0.);
	hMean->SetBinContent(bin,0);
      }
      hRes->SetBinError(bin, fitFunc->GetParError(2));
      hMean->SetBinError(bin, fitFunc->GetParError(1));
      
      //
      //

    } else {
      hRes->SetBinContent(bin, 0.);
      hRes->SetBinError(bin, 0.);
      hMean->SetBinContent(bin, 0.);
      hMean->SetBinError(bin, 0.);
    }
    

    if (drawBinFits) {
      char name[256];
      if (bin == 0) {
	sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
      } else if (bin == nBins+1) {
	sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
      } else {
	sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
		axis->GetTitle(), axis->GetBinUpEdge(bin));
      }
      canBinFits->cd(bin + dBin);
      hBin->SetTitle(name);
      hBin->SetStats(kTRUE);
      hBin->DrawCopy("E");
      canBinFits->Update();
      canBinFits->Modified();
      canBinFits->Update();
    }
    
    delete hBin;
  }

  delete fitFunc;
  currentPad->cd();
  *phMean = hMean;
  return hRes;
}


void   AliTreeDraw::GetPoints3D(const char * label, const char * chpoints, 
				const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
  //
  // load selected points from tree
  //
   if (!fPoints) fPoints = new TObjArray;
   if (tpoints->GetIndex()==0) tpoints->BuildIndex("fLabel","Label");
   TBranch * br = tpoints->GetBranch(chpoints);
   if (!br) return;
   AliTrackPointArray * points = new AliTrackPointArray;
   br->SetAddress(&points);

   Int_t npoints = fTree->Draw(label,selection);
   Float_t xyz[30000];
   rmin*=rmin;
   for (Int_t i=0;i<npoints;i++){
     Int_t index = (Int_t)fTree->GetV1()[i];
     tpoints->GetEntryWithIndex(index,index);
     if (points->GetNPoints()<2) continue;
     Int_t goodpoints=0;
     for (Int_t i=0;i<points->GetNPoints(); i++){
       xyz[goodpoints*3]   = points->GetX()[i];
       xyz[goodpoints*3+1] = points->GetY()[i];
       xyz[goodpoints*3+2] = points->GetZ()[i];
       if ( points->GetX()[i]*points->GetX()[i]+points->GetY()[i]*points->GetY()[i]>rmin) goodpoints++;
     } 
     TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz); 
     marker->SetMarkerColor(color);
     marker->SetMarkerStyle(1);
     fPoints->AddLast(marker);
   }
}


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){
   //
   // fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
   // returns chi2, fitParam and covMatrix
   // returns TString with fitted formula
   //
    
   TString formulaStr(formula); 
   TString drawStr(drawCommand);
   TString cutStr(cuts);
      
   formulaStr.ReplaceAll("++", "~");
   TObjArray* formulaTokens = formulaStr.Tokenize("~"); 
   Int_t dim = formulaTokens->GetEntriesFast();
   
   fitParam.ResizeTo(dim);
   covMatrix.ResizeTo(dim,dim);
   
   TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
   fitter->StoreData(kTRUE);   
   fitter->ClearPoints();
   
   Int_t entries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff",  stop-start, start);
   if (entries == -1) return new TString("An ERROR has occured during fitting!");
   Double_t **values = new Double_t*[dim+1] ; 
   
   for (Int_t i = 0; i < dim + 1; i++){
      Int_t centries = 0;
      if (i < dim) centries = fTree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
      else  centries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
      
      if (entries != centries) return new TString("An ERROR has occured during fitting!");
      values[i] = new Double_t[entries];
      memcpy(values[i],  fTree->GetV1(), entries*sizeof(Double_t)); 
   }
   
   // add points to the fitter
   for (Int_t i = 0; i < entries; i++){
      Double_t x[1000];
      for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
      fitter->AddPoint(x, values[dim][i], 1);
   }

   fitter->Eval();
   fitter->GetParameters(fitParam);
   fitter->GetCovarianceMatrix(covMatrix);
   chi2 = fitter->GetChisquare();
   chi2 = chi2;
   
   TString *preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = *preturnFormula; 
   
   for (Int_t iparam = 0; iparam < dim; iparam++) {
     returnFormula.Append(Form("%s*(%f)",((TObjString*)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
     if (iparam < dim-1) returnFormula.Append("+");
   }
   returnFormula.Append(" )");
   delete formulaTokens;
   delete fitter;
   delete[] values;
   return preturnFormula;
}
Commit	Line	Data
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
	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"
a36eadd7	52	#include "TObjString.h"
a36eadd7	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
	67	ClassImp(AliTreeDraw)
	68
	69
	70	AliTreeDraw::AliTreeDraw():
	71	fTree(0),
	72	fRes(0),
	73	fMean(0),
	74	fPoints(0){
	75	//
	76	// default constructor
	77	//
	78	}
	79
	80	void AliTreeDraw::ClearHisto(){
	81	//
	82	//
	83	delete fRes;
	84	delete fMean;
	85	fRes=0;
	86	fMean=0;
	87	}
	88
	89
	90
	91	TH1F * 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
	116	TH1F * 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	///////////////////////////////////////////////////////////////////////////////////
143	TH1F * 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	///////////////////////////////////////////////////////////////////////////////////
171	TH1F * 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
201	Double_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
214	void 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
222	TH1F* 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
253	TH1F* 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
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
354
355
356
357	void AliTreeDraw::GetPoints3D(const char * label, const char * chpoints,
358	const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
359	//
360	// load selected points from tree
361	//
362	if (!fPoints) fPoints = new TObjArray;
363	if (tpoints->GetIndex()==0) tpoints->BuildIndex("fLabel","Label");
364	TBranch * br = tpoints->GetBranch(chpoints);
365	if (!br) return;
366	AliTrackPointArray * points = new AliTrackPointArray;
367	br->SetAddress(&points);
368
369	Int_t npoints = fTree->Draw(label,selection);
370	Float_t xyz[30000];
371	rmin*=rmin;
372	for (Int_t i=0;i<npoints;i++){
373	Int_t index = (Int_t)fTree->GetV1()[i];
374	tpoints->GetEntryWithIndex(index,index);
375	if (points->GetNPoints()<2) continue;
376	Int_t goodpoints=0;
377	for (Int_t i=0;i<points->GetNPoints(); i++){
378	xyz[goodpoints*3] = points->GetX()[i];
379	xyz[goodpoints*3+1] = points->GetY()[i];
380	xyz[goodpoints*3+2] = points->GetZ()[i];
381	if ( points->GetX()[i]points->GetX()[i]+points->GetY()[i]points->GetY()[i]>rmin) goodpoints++;
382	}
383	TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz);
384	marker->SetMarkerColor(color);
385	marker->SetMarkerStyle(1);
386	fPoints->AddLast(marker);
387	}
388	}
389
a36eadd7	390
	391
	392
	393	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){
	394	//
	395	// fit an arbitrary function, specified by formula into the data, specified by drawCommand and cuts
	396	// returns chi2, fitParam and covMatrix
	397	// returns TString with fitted formula
	398	//
	399
	400	TString formulaStr(formula);
	401	TString drawStr(drawCommand);
	402	TString cutStr(cuts);
	403
	404	formulaStr.ReplaceAll("++", "~");
	405	TObjArray* formulaTokens = formulaStr.Tokenize("~");
	406	Int_t dim = formulaTokens->GetEntriesFast();
	407
	408	fitParam.ResizeTo(dim);
	409	covMatrix.ResizeTo(dim,dim);
	410
	411	TLinearFitter* fitter = new TLinearFitter(dim+1, Form("hyp%d",dim));
	412	fitter->StoreData(kTRUE);
	413	fitter->ClearPoints();
	414
	415	Int_t entries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start, start);
	416	if (entries == -1) return new TString("An ERROR has occured during fitting!");
	417	Double_t *values = new Double_t[dim+1] ;
	418
	419	for (Int_t i = 0; i < dim + 1; i++){
	420	Int_t centries = 0;
	421	if (i < dim) centries = fTree->Draw(((TObjString*)formulaTokens->At(i))->GetName(), cutStr.Data(), "goff", stop-start,start);
	422	else centries = fTree->Draw(drawStr.Data(), cutStr.Data(), "goff", stop-start,start);
	423
	424	if (entries != centries) return new TString("An ERROR has occured during fitting!");
	425	values[i] = new Double_t[entries];
	426	memcpy(values[i], fTree->GetV1(), entries*sizeof(Double_t));
	427	}
	428
	429	// add points to the fitter
	430	for (Int_t i = 0; i < entries; i++){
	431	Double_t x[1000];
	432	for (Int_t j=0; j<dim;j++) x[j]=values[j][i];
	433	fitter->AddPoint(x, values[dim][i], 1);
	434	}
	435
	436	fitter->Eval();
	437	fitter->GetParameters(fitParam);
	438	fitter->GetCovarianceMatrix(covMatrix);
	439	chi2 = fitter->GetChisquare();
	440	chi2 = chi2;
	441
	442	TString preturnFormula = new TString(Form("( %f+",fitParam[0])), &returnFormula = preturnFormula;
	443
	444	for (Int_t iparam = 0; iparam < dim; iparam++) {
	445	returnFormula.Append(Form("%s(%f)",((TObjString)formulaTokens->At(iparam))->GetName(),fitParam[iparam+1]));
	446	if (iparam < dim-1) returnFormula.Append("+");
	447	}
	448	returnFormula.Append(" )");
	449	delete formulaTokens;
	450	delete fitter;
	451	delete[] values;
	452	return preturnFormula;
	453	}
454