[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / DiHadronPID / AliHistToolsDiHadronPID.cxx

/************************************************************************* 
* Copyright(c) 1998-2008, 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.                  * 
**************************************************************************/

// -----------------------------------------------------------------------
//  Tools for drawing/ manipulating histograms. (NEEDS CLEANUP!)
// -----------------------------------------------------------------------
//  Author: Misha Veldhoen (misha.veldhoen@cern.ch)

#include "AliHistToolsDiHadronPID.h"

#include <iostream>
using namespace std;

#include "TLegend.h"
#include "TMath.h"
#include "TF1.h"
#include "THn.h"
#include "TObjArray.h"
#include "TStyle.h"
#include "TArray.h"

// -----------------------------------------------------------------------
TH1F* AliHistToolsDiHadronPID::RebinVariableBinning(
	const TH1F* histIn, const Double_t* binsx, Int_t Nbinsx, Bool_t density) {

	// Rebins a histogram (hin) with a variable binning to a histogram
	// with another variable binning (binsx). If the "density" flag is set,
	// then it is expected that the bins are per unit x-axis, otherwise an
	// absolute count is assumed.
	//
	// TODO: determine over/under-flow bins.

	// Gather info from the original histogram.
	const TArrayD* binsxIn = histIn->GetXaxis()->GetXbins();
	TArrayD* binsxOut = new TArrayD(Nbinsx + 1,binsx);
	Int_t NbinsxIn = binsxIn->GetSize() - 1;
	const char* nameIn = histIn->GetName();
	const char* titleIn = histIn->GetTitle();
	const char* xAxisTitleIn = histIn->GetXaxis()->GetTitle(); 
	const char* yAxisTitleIn = histIn->GetYaxis()->GetTitle();

	// Gather info from the output histogram and create it.
	Int_t NbinsxOut = binsxOut->GetSize() - 1;
	const char* nameOut = Form("%s (rebinned)",nameIn);
	const char* titleOut = Form("%s_rebinned;%s;%s",titleIn,xAxisTitleIn,yAxisTitleIn);
	TH1F* histOut = new TH1F(nameOut,titleOut,NbinsxOut,binsxOut->GetArray());
	//histOut->Sumw2();

	// Filling the bins.
	for (Int_t iBinOut = 0; iBinOut < NbinsxOut; iBinOut++) {

		Double_t binOutLowEdge = binsxOut->At(iBinOut);
		Double_t binOutHighEdge = binsxOut->At(iBinOut+1);
		Double_t binOutWidth = binOutHighEdge - binOutLowEdge;
		//Double_t binOutCenter = binOutHighEdge - binOutWidth/2.;

		// Setting all errors of the outgoing histogram to zero (needed later):
		histOut->SetBinError(iBinOut+1,0.);

		for (Int_t iBinIn = 0; iBinIn < NbinsxIn; iBinIn++) {
			
			Double_t binInLowEdge = binsxIn->At(iBinIn);
			Double_t binInHighEdge = binsxIn->At(iBinIn+1);
			Double_t binInWidth = binInHighEdge - binInLowEdge;
			//Double_t binInCenter = binInHighEdge - binInWidth/2.;
			Double_t binInContent = histIn->GetBinContent(iBinIn+1);
			Double_t binInError2 = histIn->GetBinError(iBinIn+1)*histIn->GetBinError(iBinIn+1);

			/* -- Determining what percentage of the in-bin is included 
			   in the current out-bin -- */
			Double_t PercentageIncluded = 0.;

			// Position of the low edge of the in-bin:
			// In  1   2   3            
			// Out   |   |
			Int_t binInLowEdgePos = 0;
			if (binInLowEdge < binOutLowEdge) binInLowEdgePos = 1;
			else if (binInLowEdge > binOutHighEdge) binInLowEdgePos = 3;
			else binInLowEdgePos = 2;

			// Same for the high edge of the in-bin:
			Int_t binInHighEdgePos = 0;
			if (binInHighEdge < binOutLowEdge) binInHighEdgePos = 1;
			else if (binInHighEdge > binOutHighEdge) binInHighEdgePos = 3;
			else binInHighEdgePos = 2;

			// In-bin not included.
			// In        | |     or     | |
			// Out |   |                    |    |
			if ( binInLowEdgePos == 3 || binInHighEdgePos == 1 ) {
				//cout<<"In-bin not included."<<endl;
				PercentageIncluded = 0.;
			}

			// In-bin partially included (left side).
			// In    |    |
			// Out |  |
			else if ( binInLowEdgePos == 2 && binInHighEdgePos == 3 ) {
				//cout<<"In-bin partially included (left side)."<<endl;
				PercentageIncluded = (binOutHighEdge - binInLowEdge)/(binInHighEdge - binInLowEdge);
			}

			// In-bin partially included (middle).
			// In    |     |
			// Out     |  |
			else if ( binInLowEdgePos == 1 && binInHighEdgePos == 3 ) {
				//cout<<"In-bin partially included (middle)."<<endl;
				PercentageIncluded = (binOutHighEdge - binOutLowEdge)/(binInHighEdge - binInLowEdge);
			}

			// In-bin partially included (right side).
			// In    |     |
			// Out        |  |
			else if ( binInLowEdgePos == 1 && binInHighEdgePos == 2 ) {
				//cout<<"In-bin partially included (right side)."<<endl;
				PercentageIncluded = (binInHighEdge - binOutLowEdge)/(binInHighEdge - binInLowEdge);
			}

			// In-bin completely included.
			// In    |  |
			// Out  |    |
			else if ( binInLowEdgePos == 2 && binInHighEdgePos == 2 ) {
				//cout<<"In-bin completely included."<<endl;
				PercentageIncluded = 1.;
			}

			// Give some status:
			//cout<<"Bin Out: "<<iBinOut+1<<" ("<<binOutLowEdge<<","<<binOutHighEdge<<")  Bin In: "<<iBinIn+1<<" ("<<binInLowEdge<<","<<binInHighEdge<<"), content: "<<binInContent<<" +/- "<<TMath::Sqrt(binInError2)<<". Bin in Pos: "<<binInLowEdgePos<<" "<<binInHighEdgePos<<" pct: "<<PercentageIncluded<<endl; 

			//else cout<<"Something weird just happened."<<endl;

			Double_t binOutAddedError2 = 0.;
			binOutAddedError2 = binInError2 * PercentageIncluded * PercentageIncluded;
			Double_t binOutAddedContent = 0.;
			binOutAddedContent = binInContent * PercentageIncluded;

			if (density) {
				binOutAddedContent*=binInWidth;
			}
			
			//histOut->Fill(binOutCenter,binOutAddedContent);
			histOut->SetBinContent(iBinOut+1,histOut->GetBinContent(iBinOut+1)+binOutAddedContent);

			//cout<<histOut->GetBinError(iBinOut+1)<<endl;
			//cout<<binOutAddedError2<<endl;
			histOut->SetBinError(iBinOut+1,histOut->GetBinError(iBinOut+1)+binOutAddedError2);
			//cout<<histOut->GetBinError(iBinOut+1)<<endl;
		}

		// Once a certain binOut has been filled all the way we have to scale
		// it if the "density" flag is set.

		histOut->SetBinContent(iBinOut + 1,histOut->GetBinContent(iBinOut+1)/binOutWidth);
		histOut->SetBinError(iBinOut + 1,TMath::Sqrt(histOut->GetBinError(iBinOut+1)));
		//cout<<histOut->GetBinError(iBinOut+1)<<endl;
	}

	return histOut;

}

// -----------------------------------------------------------------------
TH1F* AliHistToolsDiHadronPID::RebinVariableBinning(const TH1F* histIn, const TH1F* histAxis, Bool_t density) {

	// Rebins histogram histIn to the x-axis of histAxis
	const TAxis* xaxis = histAxis->GetXaxis();
	Int_t nbinsx = xaxis->GetNbins();
	const Double_t* binsx = (xaxis->GetXbins())->GetArray();
	return RebinVariableBinning(histIn, const_cast<Double_t*>(binsx), nbinsx, density);

}

// -----------------------------------------------------------------------
TH1F* AliHistToolsDiHadronPID::RebinVariableBinning(const TH1F* histIn, const TAxis* xaxis, Bool_t density) {

	// Rebins histogram histIn to the x-axis of histAxis
	Int_t nbinsx = xaxis->GetNbins();
	const Double_t* binsx = (xaxis->GetXbins())->GetArray();
	return RebinVariableBinning(histIn, const_cast<Double_t*>(binsx), nbinsx, density);

}

// -----------------------------------------------------------------------
TH1F* AliHistToolsDiHadronPID::TrimHisto(const TH1F* histo, Int_t firstbin, Int_t lastbin) {

	const char* name = histo->GetName();
	const char* title = histo->GetTitle();
	if (firstbin < 0) return 0x0;
	if (lastbin > histo->GetNbinsX()) return 0x0;
	if (firstbin > lastbin) return 0x0;

	Int_t Nbins = lastbin - firstbin + 1;
	cout<<firstbin<<" "<<lastbin<<" "<<Nbins<<endl;
	//const Int_t NbinsC = Nbins;
	Double_t newaxis[41];

	for (Int_t ii = firstbin; ii < lastbin+1; ii++) {
		newaxis[ii - firstbin] = histo->GetBinLowEdge(ii); 
		//cout<<ii-firstbin<<" "<<newaxis[ii - firstbin]<<endl;
	}
	newaxis[Nbins] = histo->GetBinLowEdge(lastbin+1);

	TH1F* hout = new TH1F("hout","hout",Nbins,newaxis);
	hout->SetDirectory(0);
	hout->SetTitle(title);
	hout->SetName(name);

	for (Int_t ii = 0; ii < Nbins+1; ii++) {
		hout->SetBinContent(ii,histo->GetBinContent(firstbin+ii));
		hout->SetBinError(ii,histo->GetBinError(firstbin+ii));
	} 

	return hout;

}

// -----------------------------------------------------------------------
void AliHistToolsDiHadronPID::ConstMinusHist(TH1F* histo, Float_t cc) {

	// h -> (c-h)
	Int_t nbins = histo->GetNbinsX();
	for (Int_t iBin = 0; iBin < (nbins + 1); iBin++) {
		Float_t bincontent = histo->GetBinContent(iBin);
		histo->SetBinContent(iBin,(cc - bincontent));
	}

}

// -----------------------------------------------------------------------
TH3F* AliHistToolsDiHadronPID::MakeHist3D(const char* name, const char* title, 
	Int_t nbinsX, Double_t minX, Double_t maxX,
	Int_t nbinsY, Double_t minY, Double_t maxY,
	Int_t nbinsZ, const Double_t* zaxis) {

	const Double_t* xaxis = const_cast<Double_t*>(CreateAxis(nbinsX,minX,maxX));
	const Double_t* yaxis = const_cast<Double_t*>(CreateAxis(nbinsY,minY,maxY));

	TH3F* hout = new TH3F(name,title,nbinsX,xaxis,nbinsY,yaxis,nbinsZ,zaxis);

	return hout;

}

// -----------------------------------------------------------------------
TH2F* AliHistToolsDiHadronPID::Function2DToHist2D(const TF2* function, const TH2* grid) {

	// Creates a 2D histogram of "function" using the binning of 
	// the histogram "grid". 
	//  - We assume "grid" to have fixed binwidth.
	//  - Bins are only filled if they are within the domain of the function.
	//  - Histogram takes over the color of the function.

	// Gathering info about the axes.
	const TAxis* Xaxis = grid->GetXaxis();
	Int_t NbinsX = Xaxis->GetNbins();
	const TAxis* Yaxis = grid->GetYaxis();
	Int_t NbinsY = Yaxis->GetNbins();

	// Determining function range.
	Double_t Xmin = 0.; 
	Double_t Xmax = 0.; 
	Double_t Ymin = 0.; 
	Double_t Ymax = 0.;
	function->GetRange(Xmin, Ymin, Xmax, Ymax);

	// Creating the histogram.
	TH2F* hout = new TH2F(Form("%s_hist", function->GetName()),
		Form("%s (Hist);%s;%s", function->GetTitle(), Xaxis->GetTitle(), Yaxis->GetTitle()),
		NbinsX, Xaxis->GetBinLowEdge(1), Xaxis->GetBinUpEdge(NbinsX),
		NbinsY, Yaxis->GetBinLowEdge(1), Yaxis->GetBinUpEdge(NbinsY));
	hout->SetLineColor(function->GetLineColor());
	hout->SetMarkerColor(function->GetLineColor());
	hout->SetDirectory(0);

	// Filling the histogram.
	for (Int_t iBinX = 1; iBinX < (NbinsX + 1); iBinX++) {
		for (Int_t iBinY = 1; iBinY < (NbinsY + 1); iBinY++) {
			Double_t CoordX = Xaxis->GetBinCenter(iBinX);
			Double_t CoordY = Yaxis->GetBinCenter(iBinY);

			if ( (CoordX < Xmin) || (CoordX > Xmax) || (CoordY < Ymin) || (CoordY > Ymax) ) {continue;}

			hout->SetBinContent(iBinX, iBinY, function->Eval(CoordX, CoordY));
		}		
	}

	return hout;

}

// -----------------------------------------------------------------------
TCanvas* AliHistToolsDiHadronPID::CreateSpectraComparison(const char* name, 
	const char* title, const TH1F* h1, const TH1F* h2, Int_t markerstyle, Bool_t logy) {

	// - Creates a window comparing two histograms h1, and h2.
	// - Returns an array of pointers to the objects created
	//   in this function.

	//Int_t optstat = gStyle->GetOptStat();
	gStyle->SetOptStat(0);
	//c->UseCurrentStyle();
	//gStyle->SetOptStat(optstat);

	TH1F* spectra[2];
	TH1F* division;
	Int_t spectracolor[2] = {1,2};
	Int_t divisioncolor = 4;

	// Cloning and configuring spectra.
	spectra[0] = (TH1F*)h1->Clone();
	//spectra[0]->SetDirectory(0);
	spectra[1] = (TH1F*)h2->Clone();
	//spactra[1]->SetDirectory(0);	
	for (Int_t iSpectra = 0; iSpectra < 2; iSpectra++) {
		spectra[iSpectra]->Sumw2();
		spectra[iSpectra]->SetMarkerStyle(markerstyle);
		spectra[iSpectra]->SetLineColor(spectracolor[iSpectra]);
		spectra[iSpectra]->SetMarkerColor(spectracolor[iSpectra]);
	}

	// Creating the division.
	division = (TH1F*)spectra[0]->Clone("hDivision");
	division->Divide(spectra[1]);
	division->SetLineColor(divisioncolor);
	division->SetMarkerColor(divisioncolor);
	division->GetYaxis()->SetTitle("Ratio");

	// Creating window
	//TCanvas* c = new TCanvas(name,title,0,0,400,400);
	TCanvas* c = TCanvas::MakeDefCanvas();
	c->SetName(name);
	c->SetTitle(title);
	TPad* p1 = new TPad("p1","pad1",0,0.25,1,1);
	p1->SetRightMargin(0.05);
	p1->SetBottomMargin(0.);
	TPad* p2 = new TPad("p2","pad2",0,0,1,0.25);
	p2->SetTopMargin(0.);
	p2->SetRightMargin(0.05);
	p2->SetBottomMargin(0.3);
	p2->SetFillStyle(4000);
	p1->Draw();
	p2->Draw();

	// Determining plotting range.
	Double_t max = TMath::Max(spectra[0]->GetMaximum(),spectra[1]->GetMaximum());
	Double_t min = TMath::Min(spectra[0]->GetMinimum(),spectra[1]->GetMinimum());
	Double_t range = max-min;
	spectra[0]->SetMinimum(min-0.05*range);
	spectra[0]->SetMaximum(max+0.05*range);

	// Drawing
	p1->cd();
	if (logy) {p1->SetLogy();}
	spectra[0]->Draw("e");
	spectra[0]->GetXaxis()->SetLabelSize(0);
	spectra[0]->GetYaxis()->SetLabelSize(0.05);
	spectra[0]->GetYaxis()->SetTitleSize(0.05);
	spectra[0]->GetYaxis()->SetTitleOffset(0.8);
	//Double_t labelsize = spectra[0]->GetYaxis()->GetLabelSize();
	spectra[1]->Draw("same e");
	p2->cd();
	division->SetTitle("");
	division->GetXaxis()->SetLabelSize(0.1);
	division->GetXaxis()->SetTitleSize(0.15);
	division->GetXaxis()->SetTitleOffset(0.6);
	division->GetYaxis()->SetLabelSize(0.1);
	division->GetYaxis()->SetTitleSize(0.12);
	division->GetYaxis()->SetTitleOffset(0.25);
	division->Draw("e");

	//c->UseCurrentStyle();

	//gStyle->SetOptStat(optstat);

	TLegend* legend = new TLegend(0.75,0.8,0.95,0.95);
	legend->AddEntry(spectra[0],h1->GetTitle(),"lpe");
	legend->AddEntry(spectra[1],h2->GetTitle(),"lpe");
	p1->cd();
	legend->Draw("same");

/*
	// returning the created objects.
	TObjArray* returnobjects = new TObjArray(6);
	returnobjects->AddLast(c);
	returnobjects->AddLast(p1);
	returnobjects->AddLast(p2);
	returnobjects->AddLast(spectra[0]);
	returnobjects->AddLast(spectra[1]);
	returnobjects->AddLast(division);
*/

	return c;

}

// -----------------------------------------------------------------------
Double_t* AliHistToolsDiHadronPID::CreateAxis(Int_t nbins, Double_t min, Double_t max) {

	if (nbins <= 0) return 0x0;
	if (max < min) return 0x0;

	Double_t* axis = new Double_t[nbins + 1];
	Double_t binsize = (max - min)/((Double_t)nbins);
	for (Int_t iBin = 0; iBin < nbins + 1; iBin++) {
		axis[iBin] = min + ((Double_t)iBin) * binsize;
	}

	return axis;
}
Commit	Line	Data
97724bd1	1	/*************************************************************************
	2	* Copyright(c) 1998-2008, 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	// Tools for drawing/ manipulating histograms. (NEEDS CLEANUP!)
	18	// -----------------------------------------------------------------------
	19	// Author: Misha Veldhoen (misha.veldhoen@cern.ch)
6788af99	20
a5422983	21	#include "AliHistToolsDiHadronPID.h"
a5422983	22
6788af99	23	#include <iostream>
a5422983	24	using namespace std;
a5422983	25
6788af99	26	#include "TLegend.h"
6788af99	27	#include "TMath.h"
5c01a71f	28	#include "TF1.h"
6788af99	29	#include "THn.h"
	30	#include "TObjArray.h"
	31	#include "TStyle.h"
	32	#include "TArray.h"
	33
97724bd1	34	// -----------------------------------------------------------------------
6788af99	35	TH1F* AliHistToolsDiHadronPID::RebinVariableBinning(
fe463f34	36	const TH1F* histIn, const Double_t* binsx, Int_t Nbinsx, Bool_t density) {
6788af99	37
	38	// Rebins a histogram (hin) with a variable binning to a histogram
	39	// with another variable binning (binsx). If the "density" flag is set,
	40	// then it is expected that the bins are per unit x-axis, otherwise an
	41	// absolute count is assumed.
	42	//
	43	// TODO: determine over/under-flow bins.
	44
	45	// Gather info from the original histogram.
	46	const TArrayD* binsxIn = histIn->GetXaxis()->GetXbins();
	47	TArrayD* binsxOut = new TArrayD(Nbinsx + 1,binsx);
	48	Int_t NbinsxIn = binsxIn->GetSize() - 1;
	49	const char* nameIn = histIn->GetName();
	50	const char* titleIn = histIn->GetTitle();
	51	const char* xAxisTitleIn = histIn->GetXaxis()->GetTitle();
	52	const char* yAxisTitleIn = histIn->GetYaxis()->GetTitle();
	53
	54	// Gather info from the output histogram and create it.
	55	Int_t NbinsxOut = binsxOut->GetSize() - 1;
	56	const char* nameOut = Form("%s (rebinned)",nameIn);
	57	const char* titleOut = Form("%s_rebinned;%s;%s",titleIn,xAxisTitleIn,yAxisTitleIn);
	58	TH1F* histOut = new TH1F(nameOut,titleOut,NbinsxOut,binsxOut->GetArray());
	59	//histOut->Sumw2();
	60
	61	// Filling the bins.
	62	for (Int_t iBinOut = 0; iBinOut < NbinsxOut; iBinOut++) {
	63
	64	Double_t binOutLowEdge = binsxOut->At(iBinOut);
	65	Double_t binOutHighEdge = binsxOut->At(iBinOut+1);
	66	Double_t binOutWidth = binOutHighEdge - binOutLowEdge;
	67	//Double_t binOutCenter = binOutHighEdge - binOutWidth/2.;
	68
	69	// Setting all errors of the outgoing histogram to zero (needed later):
	70	histOut->SetBinError(iBinOut+1,0.);
	71
	72	for (Int_t iBinIn = 0; iBinIn < NbinsxIn; iBinIn++) {
	73
	74	Double_t binInLowEdge = binsxIn->At(iBinIn);
	75	Double_t binInHighEdge = binsxIn->At(iBinIn+1);
	76	Double_t binInWidth = binInHighEdge - binInLowEdge;
	77	//Double_t binInCenter = binInHighEdge - binInWidth/2.;
	78	Double_t binInContent = histIn->GetBinContent(iBinIn+1);
	79	Double_t binInError2 = histIn->GetBinError(iBinIn+1)*histIn->GetBinError(iBinIn+1);
	80
	81	/* -- Determining what percentage of the in-bin is included
	82	in the current out-bin -- */
	83	Double_t PercentageIncluded = 0.;
	84
	85	// Position of the low edge of the in-bin:
	86	// In 1 2 3
	87	// Out \| \|
	88	Int_t binInLowEdgePos = 0;
	89	if (binInLowEdge < binOutLowEdge) binInLowEdgePos = 1;
	90	else if (binInLowEdge > binOutHighEdge) binInLowEdgePos = 3;
	91	else binInLowEdgePos = 2;
	92
	93	// Same for the high edge of the in-bin:
	94	Int_t binInHighEdgePos = 0;
	95	if (binInHighEdge < binOutLowEdge) binInHighEdgePos = 1;
	96	else if (binInHighEdge > binOutHighEdge) binInHighEdgePos = 3;
	97	else binInHighEdgePos = 2;
	98
	99	// In-bin not included.
	100	// In \| \| or \| \|
101	// Out \| \| \| \|
102	if ( binInLowEdgePos == 3 \|\| binInHighEdgePos == 1 ) {
103	//cout<<"In-bin not included."<<endl;
104	PercentageIncluded = 0.;
105	}
106
107	// In-bin partially included (left side).
108	// In \| \|
109	// Out \| \|
110	else if ( binInLowEdgePos == 2 && binInHighEdgePos == 3 ) {
111	//cout<<"In-bin partially included (left side)."<<endl;
112	PercentageIncluded = (binOutHighEdge - binInLowEdge)/(binInHighEdge - binInLowEdge);
113	}
114
115	// In-bin partially included (middle).
116	// In \| \|
117	// Out \| \|
118	else if ( binInLowEdgePos == 1 && binInHighEdgePos == 3 ) {
119	//cout<<"In-bin partially included (middle)."<<endl;
120	PercentageIncluded = (binOutHighEdge - binOutLowEdge)/(binInHighEdge - binInLowEdge);
121	}
122
123	// In-bin partially included (right side).
124	// In \| \|
125	// Out \| \|
126	else if ( binInLowEdgePos == 1 && binInHighEdgePos == 2 ) {
127	//cout<<"In-bin partially included (right side)."<<endl;
128	PercentageIncluded = (binInHighEdge - binOutLowEdge)/(binInHighEdge - binInLowEdge);
129	}
130
131	// In-bin completely included.
132	// In \| \|
133	// Out \| \|
134	else if ( binInLowEdgePos == 2 && binInHighEdgePos == 2 ) {
135	//cout<<"In-bin completely included."<<endl;
136	PercentageIncluded = 1.;
137	}
138
139	// Give some status:
140	//cout<<"Bin Out: "<<iBinOut+1<<" ("<<binOutLowEdge<<","<<binOutHighEdge<<") Bin In: "<<iBinIn+1<<" ("<<binInLowEdge<<","<<binInHighEdge<<"), content: "<<binInContent<<" +/- "<<TMath::Sqrt(binInError2)<<". Bin in Pos: "<<binInLowEdgePos<<" "<<binInHighEdgePos<<" pct: "<<PercentageIncluded<<endl;
141
142	//else cout<<"Something weird just happened."<<endl;
143
144	Double_t binOutAddedError2 = 0.;
145	binOutAddedError2 = binInError2 * PercentageIncluded * PercentageIncluded;
146	Double_t binOutAddedContent = 0.;
147	binOutAddedContent = binInContent * PercentageIncluded;
148
149	if (density) {
150	binOutAddedContent*=binInWidth;
151	}
152
153	//histOut->Fill(binOutCenter,binOutAddedContent);
154	histOut->SetBinContent(iBinOut+1,histOut->GetBinContent(iBinOut+1)+binOutAddedContent);
155
156	//cout<<histOut->GetBinError(iBinOut+1)<<endl;
157	//cout<<binOutAddedError2<<endl;
158	histOut->SetBinError(iBinOut+1,histOut->GetBinError(iBinOut+1)+binOutAddedError2);
159	//cout<<histOut->GetBinError(iBinOut+1)<<endl;
160	}
161
162	// Once a certain binOut has been filled all the way we have to scale
163	// it if the "density" flag is set.
164
165	histOut->SetBinContent(iBinOut + 1,histOut->GetBinContent(iBinOut+1)/binOutWidth);
166	histOut->SetBinError(iBinOut + 1,TMath::Sqrt(histOut->GetBinError(iBinOut+1)));
167	//cout<<histOut->GetBinError(iBinOut+1)<<endl;
168	}
169
170	return histOut;
171
172	}
173
97724bd1	174	// -----------------------------------------------------------------------
fe463f34	175	TH1F* AliHistToolsDiHadronPID::RebinVariableBinning(const TH1F* histIn, const TH1F* histAxis, Bool_t density) {
5c01a71f	176
5c01a71f	177	// Rebins histogram histIn to the x-axis of histAxis
700a5a26	178	const TAxis* xaxis = histAxis->GetXaxis();
5c01a71f	179	Int_t nbinsx = xaxis->GetNbins();
	180	const Double_t* binsx = (xaxis->GetXbins())->GetArray();
	181	return RebinVariableBinning(histIn, const_cast<Double_t*>(binsx), nbinsx, density);
	182
	183	}
	184
a5422983	185	// -----------------------------------------------------------------------
fe463f34	186	TH1F* AliHistToolsDiHadronPID::RebinVariableBinning(const TH1F* histIn, const TAxis* xaxis, Bool_t density) {
a5422983	187
	188	// Rebins histogram histIn to the x-axis of histAxis
	189	Int_t nbinsx = xaxis->GetNbins();
	190	const Double_t* binsx = (xaxis->GetXbins())->GetArray();
	191	return RebinVariableBinning(histIn, const_cast<Double_t*>(binsx), nbinsx, density);
	192
	193	}
	194
5c01a71f	195	// -----------------------------------------------------------------------
fe463f34	196	TH1F* AliHistToolsDiHadronPID::TrimHisto(const TH1F* histo, Int_t firstbin, Int_t lastbin) {
6788af99	197
	198	const char* name = histo->GetName();
	199	const char* title = histo->GetTitle();
	200	if (firstbin < 0) return 0x0;
	201	if (lastbin > histo->GetNbinsX()) return 0x0;
	202	if (firstbin > lastbin) return 0x0;
	203
	204	Int_t Nbins = lastbin - firstbin + 1;
	205	cout<<firstbin<<" "<<lastbin<<" "<<Nbins<<endl;
	206	//const Int_t NbinsC = Nbins;
	207	Double_t newaxis[41];
	208
	209	for (Int_t ii = firstbin; ii < lastbin+1; ii++) {
	210	newaxis[ii - firstbin] = histo->GetBinLowEdge(ii);
	211	//cout<<ii-firstbin<<" "<<newaxis[ii - firstbin]<<endl;
	212	}
	213	newaxis[Nbins] = histo->GetBinLowEdge(lastbin+1);
	214
	215	TH1F* hout = new TH1F("hout","hout",Nbins,newaxis);
	216	hout->SetDirectory(0);
	217	hout->SetTitle(title);
	218	hout->SetName(name);
	219
	220	for (Int_t ii = 0; ii < Nbins+1; ii++) {
	221	hout->SetBinContent(ii,histo->GetBinContent(firstbin+ii));
	222	hout->SetBinError(ii,histo->GetBinError(firstbin+ii));
	223	}
	224
	225	return hout;
	226
	227	}
5c01a71f	228
5c01a71f	229	// -----------------------------------------------------------------------
fe463f34	230	void AliHistToolsDiHadronPID::ConstMinusHist(TH1F* histo, Float_t cc) {
5c01a71f	231
	232	// h -> (c-h)
	233	Int_t nbins = histo->GetNbinsX();
	234	for (Int_t iBin = 0; iBin < (nbins + 1); iBin++) {
	235	Float_t bincontent = histo->GetBinContent(iBin);
	236	histo->SetBinContent(iBin,(cc - bincontent));
	237	}
	238
	239	}
	240
	241	// -----------------------------------------------------------------------
	242	TH3F* AliHistToolsDiHadronPID::MakeHist3D(const char* name, const char* title,
fe463f34	243	Int_t nbinsX, Double_t minX, Double_t maxX,
	244	Int_t nbinsY, Double_t minY, Double_t maxY,
	245	Int_t nbinsZ, const Double_t* zaxis) {
5c01a71f	246
	247	const Double_t* xaxis = const_cast<Double_t*>(CreateAxis(nbinsX,minX,maxX));
	248	const Double_t* yaxis = const_cast<Double_t*>(CreateAxis(nbinsY,minY,maxY));
	249
	250	TH3F* hout = new TH3F(name,title,nbinsX,xaxis,nbinsY,yaxis,nbinsZ,zaxis);
	251
	252	return hout;
	253
	254	}
	255
	256	// -----------------------------------------------------------------------
	257	TH2F* AliHistToolsDiHadronPID::Function2DToHist2D(const TF2* function, const TH2* grid) {
	258
	259	// Creates a 2D histogram of "function" using the binning of
	260	// the histogram "grid".
	261	// - We assume "grid" to have fixed binwidth.
	262	// - Bins are only filled if they are within the domain of the function.
	263	// - Histogram takes over the color of the function.
	264
	265	// Gathering info about the axes.
700a5a26	266	const TAxis* Xaxis = grid->GetXaxis();
5c01a71f	267	Int_t NbinsX = Xaxis->GetNbins();
700a5a26	268	const TAxis* Yaxis = grid->GetYaxis();
5c01a71f	269	Int_t NbinsY = Yaxis->GetNbins();
	270
	271	// Determining function range.
	272	Double_t Xmin = 0.;
	273	Double_t Xmax = 0.;
	274	Double_t Ymin = 0.;
	275	Double_t Ymax = 0.;
	276	function->GetRange(Xmin, Ymin, Xmax, Ymax);
	277
	278	// Creating the histogram.
	279	TH2F* hout = new TH2F(Form("%s_hist", function->GetName()),
	280	Form("%s (Hist);%s;%s", function->GetTitle(), Xaxis->GetTitle(), Yaxis->GetTitle()),
	281	NbinsX, Xaxis->GetBinLowEdge(1), Xaxis->GetBinUpEdge(NbinsX),
	282	NbinsY, Yaxis->GetBinLowEdge(1), Yaxis->GetBinUpEdge(NbinsY));
	283	hout->SetLineColor(function->GetLineColor());
	284	hout->SetMarkerColor(function->GetLineColor());
	285	hout->SetDirectory(0);
	286
	287	// Filling the histogram.
	288	for (Int_t iBinX = 1; iBinX < (NbinsX + 1); iBinX++) {
	289	for (Int_t iBinY = 1; iBinY < (NbinsY + 1); iBinY++) {
	290	Double_t CoordX = Xaxis->GetBinCenter(iBinX);
	291	Double_t CoordY = Yaxis->GetBinCenter(iBinY);
	292
	293	if ( (CoordX < Xmin) \|\| (CoordX > Xmax) \|\| (CoordY < Ymin) \|\| (CoordY > Ymax) ) {continue;}
	294
	295	hout->SetBinContent(iBinX, iBinY, function->Eval(CoordX, CoordY));
	296	}
	297	}
	298
	299	return hout;
	300
	301	}
	302
	303	// -----------------------------------------------------------------------
07d62e30	304	TCanvas* AliHistToolsDiHadronPID::CreateSpectraComparison(const char* name,
fe463f34	305	const char* title, const TH1F* h1, const TH1F* h2, Int_t markerstyle, Bool_t logy) {
5c01a71f	306
	307	// - Creates a window comparing two histograms h1, and h2.
	308	// - Returns an array of pointers to the objects created
	309	// in this function.
	310
	311	//Int_t optstat = gStyle->GetOptStat();
	312	gStyle->SetOptStat(0);
	313	//c->UseCurrentStyle();
	314	//gStyle->SetOptStat(optstat);
	315
	316	TH1F* spectra[2];
	317	TH1F* division;
	318	Int_t spectracolor[2] = {1,2};
	319	Int_t divisioncolor = 4;
	320
	321	// Cloning and configuring spectra.
	322	spectra[0] = (TH1F*)h1->Clone();
	323	//spectra[0]->SetDirectory(0);
	324	spectra[1] = (TH1F*)h2->Clone();
	325	//spactra[1]->SetDirectory(0);
	326	for (Int_t iSpectra = 0; iSpectra < 2; iSpectra++) {
	327	spectra[iSpectra]->Sumw2();
	328	spectra[iSpectra]->SetMarkerStyle(markerstyle);
	329	spectra[iSpectra]->SetLineColor(spectracolor[iSpectra]);
	330	spectra[iSpectra]->SetMarkerColor(spectracolor[iSpectra]);
	331	}
	332
	333	// Creating the division.
	334	division = (TH1F*)spectra[0]->Clone("hDivision");
	335	division->Divide(spectra[1]);
	336	division->SetLineColor(divisioncolor);
	337	division->SetMarkerColor(divisioncolor);
	338	division->GetYaxis()->SetTitle("Ratio");
	339
	340	// Creating window
	341	//TCanvas* c = new TCanvas(name,title,0,0,400,400);
	342	TCanvas* c = TCanvas::MakeDefCanvas();
	343	c->SetName(name);
	344	c->SetTitle(title);
	345	TPad* p1 = new TPad("p1","pad1",0,0.25,1,1);
	346	p1->SetRightMargin(0.05);
	347	p1->SetBottomMargin(0.);
	348	TPad* p2 = new TPad("p2","pad2",0,0,1,0.25);
	349	p2->SetTopMargin(0.);
	350	p2->SetRightMargin(0.05);
	351	p2->SetBottomMargin(0.3);
	352	p2->SetFillStyle(4000);
	353	p1->Draw();
	354	p2->Draw();
	355
	356	// Determining plotting range.
	357	Double_t max = TMath::Max(spectra[0]->GetMaximum(),spectra[1]->GetMaximum());
	358	Double_t min = TMath::Min(spectra[0]->GetMinimum(),spectra[1]->GetMinimum());
	359	Double_t range = max-min;
	360	spectra[0]->SetMinimum(min-0.05*range);
	361	spectra[0]->SetMaximum(max+0.05*range);
	362
	363	// Drawing
	364	p1->cd();
	365	if (logy) {p1->SetLogy();}
	366	spectra[0]->Draw("e");
	367	spectra[0]->GetXaxis()->SetLabelSize(0);
	368	spectra[0]->GetYaxis()->SetLabelSize(0.05);
	369	spectra[0]->GetYaxis()->SetTitleSize(0.05);
a5422983	370	spectra[0]->GetYaxis()->SetTitleOffset(0.8);
5c01a71f	371	//Double_t labelsize = spectra[0]->GetYaxis()->GetLabelSize();
	372	spectra[1]->Draw("same e");
	373	p2->cd();
	374	division->SetTitle("");
	375	division->GetXaxis()->SetLabelSize(0.1);
	376	division->GetXaxis()->SetTitleSize(0.15);
	377	division->GetXaxis()->SetTitleOffset(0.6);
	378	division->GetYaxis()->SetLabelSize(0.1);
	379	division->GetYaxis()->SetTitleSize(0.12);
a5422983	380	division->GetYaxis()->SetTitleOffset(0.25);
5c01a71f	381	division->Draw("e");
	382
	383	//c->UseCurrentStyle();
	384
	385	//gStyle->SetOptStat(optstat);
	386
	387	TLegend* legend = new TLegend(0.75,0.8,0.95,0.95);
	388	legend->AddEntry(spectra[0],h1->GetTitle(),"lpe");
	389	legend->AddEntry(spectra[1],h2->GetTitle(),"lpe");
	390	p1->cd();
	391	legend->Draw("same");
	392
07d62e30	393	/*
5c01a71f	394	// returning the created objects.
	395	TObjArray* returnobjects = new TObjArray(6);
	396	returnobjects->AddLast(c);
	397	returnobjects->AddLast(p1);
	398	returnobjects->AddLast(p2);
	399	returnobjects->AddLast(spectra[0]);
	400	returnobjects->AddLast(spectra[1]);
	401	returnobjects->AddLast(division);
07d62e30	402	*/
5c01a71f	403
07d62e30	404	return c;
5c01a71f	405
	406	}
	407
	408	// -----------------------------------------------------------------------
fe463f34	409	Double_t* AliHistToolsDiHadronPID::CreateAxis(Int_t nbins, Double_t min, Double_t max) {
5c01a71f	410
	411	if (nbins <= 0) return 0x0;
	412	if (max < min) return 0x0;
	413
	414	Double_t* axis = new Double_t[nbins + 1];
	415	Double_t binsize = (max - min)/((Double_t)nbins);
	416	for (Int_t iBin = 0; iBin < nbins + 1; iBin++) {
	417	axis[iBin] = min + ((Double_t)iBin) * binsize;
	418	}
	419
	420	return axis;
	421	}
	422
	423
	424
	425
	426
	427