[u/mrichter/AliRoot.git] / FMD / flow / AliFMDFlowBinned1D.cxx

/** @file 
    @brief Implementation of a 1-dimensional Flow "histogram" */
#include "flow/AliFMDFlowBinned1D.h"
#include "flow/AliFMDFlowBin.h"
#include <cmath>
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <TBrowser.h>
#include <TString.h>
#include <TH1.h>
#include <TH2.h>

//====================================================================
AliFMDFlowBinned1D::AliFMDFlowBinned1D(UShort_t order, 
				       UShort_t nxbins, 
				       Double_t* xbins) 
  : fXAxis(nxbins, xbins),
    fBins(0)
{
  UShort_t n = fXAxis.N();
  fBins   = new AliFMDFlowBin*[n];
  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order);
}
//____________________________________________________________________
AliFMDFlowBinned1D::AliFMDFlowBinned1D(UShort_t order, 
				       const AliFMDFlowAxis& xaxis)
  : fXAxis(xaxis)
{
  UShort_t n = fXAxis.N();
  fBins   = new AliFMDFlowBin*[n];
  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order);
}
//____________________________________________________________________
AliFMDFlowBinned1D::AliFMDFlowBinned1D(const AliFMDFlowBinned1D& o)
  : TObject(o), 
    fXAxis(o.fXAxis)
{
  UShort_t n = fXAxis.N();
  fBins   = new AliFMDFlowBin*[n];
  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
}
//____________________________________________________________________
AliFMDFlowBinned1D&
AliFMDFlowBinned1D::operator=(const AliFMDFlowBinned1D& o)
{
  if (fBins) { 
    for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
    delete [] fBins;
  }
  fXAxis     = o.fXAxis;
  UShort_t n = fXAxis.N();
  fBins   = new AliFMDFlowBin*[n];
  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
  return *this;
}
  
//____________________________________________________________________
AliFMDFlowBinned1D::~AliFMDFlowBinned1D()
{
  if (fBins) { 
    for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
    delete [] fBins;
  }
}

//____________________________________________________________________
AliFMDFlowBin* 
AliFMDFlowBinned1D::GetBin(UShort_t i) const
{
  if (i >= fXAxis.N()) return 0;
  return fBins[i];
}
//____________________________________________________________________
AliFMDFlowBin* 
AliFMDFlowBinned1D::GetBin(Double_t x) const
{
  Int_t i = fXAxis.FindBin(x);
  if (i < 0) return 0;
  UShort_t j = i;
  return GetBin(j);
}
  
//____________________________________________________________________
void 
AliFMDFlowBinned1D::Begin()
{
  for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Begin();
}
//____________________________________________________________________
void 
AliFMDFlowBinned1D::End()
{
  for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->End();
}
//____________________________________________________________________
void 
AliFMDFlowBinned1D::Finish()
{
  for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Finish();
}
//____________________________________________________________________
Bool_t 
AliFMDFlowBinned1D::AddToEventPlane(Double_t x, Double_t phi, Double_t w, Bool_t a)
{
  AliFMDFlowBin* bin = GetBin(x);
  if (!bin) return kFALSE;
  bin->AddToEventPlane(phi, w, a);
  return kTRUE;
}

//____________________________________________________________________
Bool_t 
AliFMDFlowBinned1D::AddToHarmonic(Double_t x, Double_t phi)
{
  AliFMDFlowBin* bin = GetBin(x);
  if (!bin) return kFALSE;
  bin->AddToHarmonic(phi);
  return kTRUE;
}

//____________________________________________________________________
void 
AliFMDFlowBinned1D::Event(Double_t* phis, Double_t* xs, Double_t* ws, ULong_t n)
{
  Begin();
  for (UInt_t i = 0; i < n; i++) 
    AddToEventPlane(xs[i], phis[i], (ws ? ws[i] : 1), 
		    Float_t(rand()) / RAND_MAX > 0.5);
  for (UInt_t i = 0; i < n; i++) 
    AddToHarmonic(xs[i], phis[i]);
  End();
}

//____________________________________________________________________
void 
AliFMDFlowBinned1D::Browse(TBrowser* b)
{
  b->Add(&fXAxis, "xaxis");
  for (UInt_t i = 0; i < fXAxis.N(); i++) 
    b->Add(fBins[i], Form("bin_%03d", i));
}

//____________________________________________________________________
void 
AliFMDFlowBinned1D::Draw(Option_t* option)
{
  TString opt(option);
  opt.ToLower();
  const char* names[] = { "Bare", "Naive", "STAR", "TDR" };
  const AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::none, 
					   AliFMDFlowBin::naive, 
					   AliFMDFlowBin::star, 
					   AliFMDFlowBin::tdr };
  Bool_t meths[] = { opt.Contains("b"), 
		     opt.Contains("n"),
		     opt.Contains("s"), 
		     opt.Contains("t") };
  Bool_t res     = opt.Contains("r");
  UShort_t nm = 0;
  Short_t  sm = -1;
  for (UShort_t i = 0; i < 4; i++) { if (meths[i]) { nm++; sm = i; } }
  TH1* h = 0;
  if (nm > 1) { 
    h = new TH2D((res ? "res" : "flow"), (res ? "Resolution" : "Flow"), 
		 fXAxis.N(), fXAxis.Bins(), nm, 0, nm);
    h->SetXTitle("x");
    h->SetYTitle("method");
    h->GetYaxis()->SetNdivisions(nm+1, kFALSE);
    h->SetZTitle((res ? "<cos(n(#Psi_{m}-#Psi_{R}))>" : "v"));
    UInt_t j = 0;
    for (UShort_t i = 0; i < 4; i++) 
      if (meths[i]) h->GetYaxis()->SetBinLabel(++j, names[i]);
  }
  else {
    h = new TH1D(Form("%s_%s", (res ? "res" : "flow"), names[sm]),
		 Form("%s_%s", (res ? "Resolution" : "Flow"), names[sm]),
		 fXAxis.N(), fXAxis.Bins());
    h->SetXTitle("x");
    h->SetYTitle((res ? "<cos(n(#Psi_{m}-#Psi_{R}))>" : "v"));
  }

  for (UShort_t i = 0; i < fXAxis.N(); i++) { 
    Double_t       x   = fXAxis.BinCenter(i);
    AliFMDFlowBin* bin = GetBin(x);
    Double_t       v, e2;
    if (nm == 1) { 
      if (res) v = bin->Correction(e2, types[sm]);
      else     v = bin->Value(e2, types[sm]);
      h->SetBinContent(i+1, v);
      h->SetBinError(i+1, sqrt(e2));
      continue;
    }
    UInt_t j = 0;
    for (UShort_t k = 0; k < 4; k++)  { 
      if (!meths[k]) continue;
      if (res) v = bin->Correction(e2, types[k]);
      else     v = bin->Value(e2, types[k]);
      h->SetBinContent(i+1, j+1, v);
      h->SetBinError(i+1, j+1, sqrt(e2));
      j++;
    }
  }
  h->Draw(option);
}

  
  
//____________________________________________________________________
void 
AliFMDFlowBinned1D::Print(Option_t* option) const
{
  TString opt(option);
  opt.ToLower();
  Bool_t det = opt.Contains("d");
  Bool_t sum = opt.Contains("s");
  if (det) { 
    for (UShort_t i = 0; i < fXAxis.N(); i++) { 
      Double_t x = fXAxis.BinCenter(i);
      std::streamsize         old_p = std::cout.precision(3);
      std::ios_base::fmtflags old_f = std::cout.setf(std::ios_base::fixed, 
						     std::ios_base::floatfield);
      std::cout << "x=" << std::setw(5) << x << std::endl;
      fBins[i]->Print();
      std::cout.precision(old_p);
      std::cout.setf(old_f, std::ios_base::floatfield);
    }
  }
  
  if (sum) { 
    UInt_t       nType = 4;
    const char*  names[] = { "Bare",    "Naive",    "STAR",    "TDR" };
    AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::none, 
				       AliFMDFlowBin::naive, 
				       AliFMDFlowBin::star, 
				       AliFMDFlowBin::tdr };
    std::cout << "    x";
    for (UInt_t i = 0; i < nType; i++) 
      std::cout << " | " << std::setw(6+6+5) << names[i];
    std::cout << "\n-----" << std::setfill('-');
    for (UInt_t i = 0; i < nType; i++) 
      std::cout << "-+-" <<  std::setw(6+6+5) << "-";
    std::cout << std::setfill(' ') << std::endl;
    
    std::streamsize         old_p = std::cout.precision(2);
    std::ios_base::fmtflags old_f = std::cout.setf(std::ios_base::fixed, 
						   std::ios_base::floatfield);
    for (UShort_t i = 0; i < fXAxis.N(); i++) { 
      Double_t x = fXAxis.BinCenter(i);
      std::cout << std::setprecision(2) << std::setw(5) << x << std::flush;
      for (UShort_t j = 0; j < nType; j++) { 
	Double_t e2v;
	Double_t v  = fBins[i]->Value(e2v, types[j]);
	std::cout << std::setprecision(3)    << " | "
		  << std::setw(6) << 100 * v << " +/- " 
		  << std::setw(6) << 100 * sqrt(e2v);
      }
      std::cout << std::endl;
    }
    std::cout.precision(old_p);
    std::cout.setf(old_f, std::ios_base::floatfield);
  }
}

    
//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
39eefe19	1	/** @file
	2	@brief Implementation of a 1-dimensional Flow "histogram" */
	3	#include "flow/AliFMDFlowBinned1D.h"
	4	#include "flow/AliFMDFlowBin.h"
	5	#include <cmath>
	6	#include <cstdlib>
	7	#include <iostream>
	8	#include <iomanip>
	9	#include <TBrowser.h>
	10	#include <TString.h>
824e71c1	11	#include <TH1.h>
824e71c1	12	#include <TH2.h>
39eefe19	13
	14	//====================================================================
	15	AliFMDFlowBinned1D::AliFMDFlowBinned1D(UShort_t order,
	16	UShort_t nxbins,
	17	Double_t* xbins)
	18	: fXAxis(nxbins, xbins),
	19	fBins(0)
	20	{
	21	UShort_t n = fXAxis.N();
	22	fBins = new AliFMDFlowBin*[n];
	23	for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order);
	24	}
	25	//____________________________________________________________________
	26	AliFMDFlowBinned1D::AliFMDFlowBinned1D(UShort_t order,
	27	const AliFMDFlowAxis& xaxis)
	28	: fXAxis(xaxis)
	29	{
	30	UShort_t n = fXAxis.N();
	31	fBins = new AliFMDFlowBin*[n];
	32	for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order);
	33	}
	34	//____________________________________________________________________
	35	AliFMDFlowBinned1D::AliFMDFlowBinned1D(const AliFMDFlowBinned1D& o)
2219d590	36	: TObject(o),
2219d590	37	fXAxis(o.fXAxis)
39eefe19	38	{
	39	UShort_t n = fXAxis.N();
	40	fBins = new AliFMDFlowBin*[n];
	41	for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
	42	}
	43	//____________________________________________________________________
	44	AliFMDFlowBinned1D&
	45	AliFMDFlowBinned1D::operator=(const AliFMDFlowBinned1D& o)
	46	{
	47	if (fBins) {
	48	for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
	49	delete [] fBins;
	50	}
	51	fXAxis = o.fXAxis;
	52	UShort_t n = fXAxis.N();
	53	fBins = new AliFMDFlowBin*[n];
	54	for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
	55	return *this;
	56	}
	57
	58	//____________________________________________________________________
	59	AliFMDFlowBinned1D::~AliFMDFlowBinned1D()
	60	{
	61	if (fBins) {
	62	for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
	63	delete [] fBins;
	64	}
	65	}
	66
	67	//____________________________________________________________________
	68	AliFMDFlowBin*
	69	AliFMDFlowBinned1D::GetBin(UShort_t i) const
	70	{
	71	if (i >= fXAxis.N()) return 0;
	72	return fBins[i];
	73	}
	74	//____________________________________________________________________
	75	AliFMDFlowBin*
	76	AliFMDFlowBinned1D::GetBin(Double_t x) const
	77	{
	78	Int_t i = fXAxis.FindBin(x);
	79	if (i < 0) return 0;
	80	UShort_t j = i;
	81	return GetBin(j);
	82	}
	83
	84	//____________________________________________________________________
	85	void
	86	AliFMDFlowBinned1D::Begin()
	87	{
	88	for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Begin();
	89	}
	90	//____________________________________________________________________
	91	void
	92	AliFMDFlowBinned1D::End()
	93	{
	94	for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->End();
	95	}
	96	//____________________________________________________________________
	97	void
	98	AliFMDFlowBinned1D::Finish()
	99	{
	100	for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Finish();
	101	}
102	//____________________________________________________________________
103	Bool_t
104	AliFMDFlowBinned1D::AddToEventPlane(Double_t x, Double_t phi, Double_t w, Bool_t a)
105	{
106	AliFMDFlowBin* bin = GetBin(x);
107	if (!bin) return kFALSE;
108	bin->AddToEventPlane(phi, w, a);
109	return kTRUE;
110	}
111
112	//____________________________________________________________________
113	Bool_t
114	AliFMDFlowBinned1D::AddToHarmonic(Double_t x, Double_t phi)
115	{
116	AliFMDFlowBin* bin = GetBin(x);
117	if (!bin) return kFALSE;
118	bin->AddToHarmonic(phi);
119	return kTRUE;
120	}
121
122	//____________________________________________________________________
123	void
124	AliFMDFlowBinned1D::Event(Double_t* phis, Double_t* xs, Double_t* ws, ULong_t n)
125	{
126	Begin();
127	for (UInt_t i = 0; i < n; i++)
128	AddToEventPlane(xs[i], phis[i], (ws ? ws[i] : 1),
129	Float_t(rand()) / RAND_MAX > 0.5);
130	for (UInt_t i = 0; i < n; i++)
131	AddToHarmonic(xs[i], phis[i]);
132	End();
133	}
134
135	//____________________________________________________________________
136	void
137	AliFMDFlowBinned1D::Browse(TBrowser* b)
138	{
139	b->Add(&fXAxis, "xaxis");
140	for (UInt_t i = 0; i < fXAxis.N(); i++)
141	b->Add(fBins[i], Form("bin_%03d", i));
142	}
143
824e71c1	144	//____________________________________________________________________
	145	void
	146	AliFMDFlowBinned1D::Draw(Option_t* option)
	147	{
	148	TString opt(option);
	149	opt.ToLower();
	150	const char* names[] = { "Bare", "Naive", "STAR", "TDR" };
	151	const AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::none,
	152	AliFMDFlowBin::naive,
	153	AliFMDFlowBin::star,
	154	AliFMDFlowBin::tdr };
	155	Bool_t meths[] = { opt.Contains("b"),
	156	opt.Contains("n"),
	157	opt.Contains("s"),
	158	opt.Contains("t") };
	159	Bool_t res = opt.Contains("r");
	160	UShort_t nm = 0;
	161	Short_t sm = -1;
	162	for (UShort_t i = 0; i < 4; i++) { if (meths[i]) { nm++; sm = i; } }
	163	TH1* h = 0;
	164	if (nm > 1) {
	165	h = new TH2D((res ? "res" : "flow"), (res ? "Resolution" : "Flow"),
	166	fXAxis.N(), fXAxis.Bins(), nm, 0, nm);
	167	h->SetXTitle("x");
	168	h->SetYTitle("method");
	169	h->GetYaxis()->SetNdivisions(nm+1, kFALSE);
	170	h->SetZTitle((res ? "<cos(n(#Psi_{m}-#Psi_{R}))>" : "v"));
	171	UInt_t j = 0;
	172	for (UShort_t i = 0; i < 4; i++)
	173	if (meths[i]) h->GetYaxis()->SetBinLabel(++j, names[i]);
	174	}
	175	else {
	176	h = new TH1D(Form("%s_%s", (res ? "res" : "flow"), names[sm]),
	177	Form("%s_%s", (res ? "Resolution" : "Flow"), names[sm]),
	178	fXAxis.N(), fXAxis.Bins());
	179	h->SetXTitle("x");
	180	h->SetYTitle((res ? "<cos(n(#Psi_{m}-#Psi_{R}))>" : "v"));
	181	}
	182
	183	for (UShort_t i = 0; i < fXAxis.N(); i++) {
	184	Double_t x = fXAxis.BinCenter(i);
	185	AliFMDFlowBin* bin = GetBin(x);
	186	Double_t v, e2;
	187	if (nm == 1) {
	188	if (res) v = bin->Correction(e2, types[sm]);
	189	else v = bin->Value(e2, types[sm]);
	190	h->SetBinContent(i+1, v);
	191	h->SetBinError(i+1, sqrt(e2));
	192	continue;
	193	}
	194	UInt_t j = 0;
	195	for (UShort_t k = 0; k < 4; k++) {
	196	if (!meths[k]) continue;
	197	if (res) v = bin->Correction(e2, types[k]);
	198	else v = bin->Value(e2, types[k]);
	199	h->SetBinContent(i+1, j+1, v);
	200	h->SetBinError(i+1, j+1, sqrt(e2));
	201	j++;
	202	}
	203	}
	204	h->Draw(option);
	205	}
	206
	207
208
39eefe19	209	//____________________________________________________________________
	210	void
	211	AliFMDFlowBinned1D::Print(Option_t* option) const
	212	{
	213	TString opt(option);
	214	opt.ToLower();
	215	Bool_t det = opt.Contains("d");
	216	Bool_t sum = opt.Contains("s");
	217	if (det) {
	218	for (UShort_t i = 0; i < fXAxis.N(); i++) {
	219	Double_t x = fXAxis.BinCenter(i);
	220	std::streamsize old_p = std::cout.precision(3);
	221	std::ios_base::fmtflags old_f = std::cout.setf(std::ios_base::fixed,
	222	std::ios_base::floatfield);
	223	std::cout << "x=" << std::setw(5) << x << std::endl;
	224	fBins[i]->Print();
	225	std::cout.precision(old_p);
	226	std::cout.setf(old_f, std::ios_base::floatfield);
	227	}
	228	}
	229
	230	if (sum) {
	231	UInt_t nType = 4;
	232	const char* names[] = { "Bare", "Naive", "STAR", "TDR" };
	233	AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::none,
	234	AliFMDFlowBin::naive,
	235	AliFMDFlowBin::star,
	236	AliFMDFlowBin::tdr };
	237	std::cout << " x";
	238	for (UInt_t i = 0; i < nType; i++)
	239	std::cout << " \| " << std::setw(6+6+5) << names[i];
	240	std::cout << "\n-----" << std::setfill('-');
	241	for (UInt_t i = 0; i < nType; i++)
	242	std::cout << "-+-" << std::setw(6+6+5) << "-";
	243	std::cout << std::setfill(' ') << std::endl;
	244
	245	std::streamsize old_p = std::cout.precision(2);
	246	std::ios_base::fmtflags old_f = std::cout.setf(std::ios_base::fixed,
	247	std::ios_base::floatfield);
	248	for (UShort_t i = 0; i < fXAxis.N(); i++) {
	249	Double_t x = fXAxis.BinCenter(i);
	250	std::cout << std::setprecision(2) << std::setw(5) << x << std::flush;
	251	for (UShort_t j = 0; j < nType; j++) {
	252	Double_t e2v;
	253	Double_t v = fBins[i]->Value(e2v, types[j]);
	254	std::cout << std::setprecision(3) << " \| "
	255	<< std::setw(6) << 100 * v << " +/- "
	256	<< std::setw(6) << 100 * sqrt(e2v);
	257	}
	258	std::cout << std::endl;
	259	}
	260	std::cout.precision(old_p);
	261	std::cout.setf(old_f, std::ios_base::floatfield);
	262	}
	263	}
	264
	265
	266	//____________________________________________________________________
	267	//
	268	// EOF
	269	//