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

/* Copyright (C) 2007 Christian Holm Christensen <cholm@nbi.dk>
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1 of
 * the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA
 */
/** @file 
    @brief Implementation of a 1-dimensional Flow "histogram" */
//____________________________________________________________________ 
//
// A histogram of flow bins.  The axis can by anything
// (pseudo-rapidity, transvers momentum) - there's no assumption on
// what is the basis of the histogram.  The method Event can be used
// to calculate everything in one go.   Alternatively, one can use the
// methods AddToEventPlane and AddToHarmonic.  See also the example
// TestFlow.C 
#include "flow/AliFMDFlowBinned1D.h"
#include "flow/AliFMDFlowBin.h"
#include "flow/AliFMDFlowSplitter.h"
// #include <cmath>
#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <TBrowser.h>
#include <TString.h>
#include <TH1.h>
#include <TH2.h>

//====================================================================
AliFMDFlowBinned1D::AliFMDFlowBinned1D()
  : TNamed("", ""), 
    fXAxis(),
    fN(0),
    fBins(0),
    fSplitter(0)
{
  // Default CTOR - do not use
}

//____________________________________________________________________
AliFMDFlowBinned1D::AliFMDFlowBinned1D(const char* name, 
				       const char* title, 
				       UShort_t order, 
				       UShort_t k,
				       UShort_t nxbins, 
				       Double_t* xbins, 
				       AliFMDFlowSplitter* splitter) 
  : TNamed(name,title), 
    fXAxis(nxbins, xbins),
    fBins(0), 
    fSplitter(splitter)
{
  // Constructor 
  // Parameters: 
  //   Order	Order 
  //   nxbins   Number of bins 
  //   xbins    Bin borders 
  UShort_t n = fXAxis.N();
  fN         = n;
  fBins      = new AliFMDFlowBin*[n];
  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order,k);
  if (!fSplitter) fSplitter = new AliFMDFlowShuffle;
}
//____________________________________________________________________
AliFMDFlowBinned1D::AliFMDFlowBinned1D(const char*           name, 
				       const char*           title, 
				       UShort_t              order, 
				       UShort_t              k,
				       const AliFMDFlowAxis& xaxis,
				       AliFMDFlowSplitter*   splitter)
  : TNamed(name,title), 
    fXAxis(xaxis), 
    fBins(0), 
    fSplitter(splitter)
{
  // Constructor 
  // Parameters: 
  //   Order	Order 
  //   xaxis    X axis object
  UShort_t n = fXAxis.N();
  fN         = n;
  fBins      = new AliFMDFlowBin*[n];
  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order,k);
  if (!fSplitter) fSplitter = new AliFMDFlowShuffle;
    
}
//____________________________________________________________________
AliFMDFlowBinned1D::AliFMDFlowBinned1D(const AliFMDFlowBinned1D& o)
  : TNamed(o), 
    TAttLine(o),
    TAttFill(o),
    TAttMarker(o),
    fXAxis(o.fXAxis), 
    fSplitter(0)
{
  // Copy constructor 
  // Parameters: 
  //   o   Object to copy from 
  UShort_t n = fXAxis.N();
  fSplitter  = new AliFMDFlowShuffle;
  fN         = 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)
{
  // Assignment operator
  // Parameters: 
  //   o Object to assign from 
  // 
  // Returns reference to this object 
  SetLineColor(o.GetLineColor());
  SetLineStyle(o.GetLineStyle());
  SetLineWidth(o.GetLineWidth());
  SetFillColor(o.GetFillColor());
  SetFillStyle(o.GetFillStyle());
  SetMarkerColor(o.GetMarkerColor());
  SetMarkerStyle(o.GetMarkerStyle());
  SetMarkerSize(o.GetMarkerSize());
  this->SetName(o.GetName());
  this->SetTitle(o.GetTitle());
  if (fBins) { 
    for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
    delete [] fBins;
  }
  if (fSplitter) delete fSplitter;
  fXAxis     = o.fXAxis;
  UShort_t n = fXAxis.N();
  fN         = n;
  fSplitter  = new AliFMDFlowShuffle;
  fBins      = new AliFMDFlowBin*[n];
  for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
  return *this;
}
  
//____________________________________________________________________
AliFMDFlowBinned1D::~AliFMDFlowBinned1D()
{
  // Destructor 
  // Parameters: 
  //    none 
  if (fBins) { 
    for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
    delete [] fBins;
  }
  if (fSplitter) delete fSplitter;
}

//____________________________________________________________________
AliFMDFlowBin* 
AliFMDFlowBinned1D::GetBin(UShort_t i) const
{
  // Get the ith bin 
  // Parameters: 
  //    i  Bin number 
  // 
  // Return pointer to bin, or null. 
  if (i >= fXAxis.N()) return 0;
  return fBins[i];
}
//____________________________________________________________________
AliFMDFlowBin* 
AliFMDFlowBinned1D::GetBin(Double_t x) const
{
  // Get the bin that contains x
  // Parameters: 
  //    x  X axis value
  // 
  // Return pointer to bin, or null. 
  Int_t i = fXAxis.FindBin(x);
  if (i < 0) return 0;
  UShort_t j = i;
  return GetBin(j);
}
  
//____________________________________________________________________
void 
AliFMDFlowBinned1D::Begin()
{
  // Called at the beginning of an event
  // Parameters: 
  //   none
  for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Begin();
  fSplitter->Begin();
}
//____________________________________________________________________
void 
AliFMDFlowBinned1D::End()
{
  // Called at the end of an event
  // Parameters: 
  //   none
  for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->End();
  fSplitter->End();
}
//____________________________________________________________________
void 
AliFMDFlowBinned1D::Finish()
{
  // Called at the end of an job
  // Parameters: 
  //   none
  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)
{
  // Called to add a contribution to the event plane 
  // Parameters:
  //    x   Bin to fill into 
  //    w   Weight
  //    phi The angle phi in radians 
  //    a   If true, add to sub-event A, otherwise sub-event B
  // 
  // Return false if x falls outside the defined range, true otherwise
  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, 
				  Double_t wp, Double_t wh)
{
  // Called to add a contribution to the harmonic
  // Parameters: 
  //    x    Bin to fill into 
  //    phi  The angle phi in radians
  //    wp   weight of event plane
  //    wh   weight of harmonic

  // Return false if x falls outside the defined range, true otherwise
  AliFMDFlowBin* bin = GetBin(x);
  if (!bin) return kFALSE;
  bin->AddToHarmonic(phi, wp, wh);
  return kTRUE;
}

//____________________________________________________________________
UShort_t
AliFMDFlowBinned1D::Order() const 
{ 
  return GetBin(UShort_t(0))->Order(); 
}
//____________________________________________________________________
UShort_t
AliFMDFlowBinned1D::PsiOrder() const 
{ 
  return GetBin(UShort_t(0))->Order(); 
}

//____________________________________________________________________
void 
AliFMDFlowBinned1D::Event(ULong_t   n,  Double_t* phis, Double_t* xs, 
			  Double_t* wp, Double_t* wh)
{
  // Process a full event. 
  // Parameters: 
  //   phis   List of n phi=[0,2pi] angles 
  //   xs     List of n x values. 
  //   ws     Weights
  //   n      Size of phis and xs
  Begin();
  fSplitter->Event(phis, xs, n);
  for (ULong_t i = 0; i < n; i++) 
    AddToEventPlane(xs[i], phis[i], (wp ? wp[i] : 1), fSplitter->Select(i));
  for (ULong_t i = 0; i < n; i++) 
    AddToHarmonic(xs[i], phis[i], (wp ? wp[i] : 1), (wh ? wh[i] : 1));
  End();
}

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

//____________________________________________________________________
TH1*
AliFMDFlowBinned1D::MakeHistogram(UInt_t which, UInt_t what)
{
  // Make a histogram of some of the stuff in this object 

  // Some strings 
  const char* names[]   = { "Bare",  "Naive", "STAR", "TDR" };
  const char* hNames[]  = { "flow",  "resolution", "counts" }; 
  const char* hTitles[] = { "Flow",  "Resolution", "Counts" }; 
  const char* hUnits[]  = { "v_{%d}","<cos(%d(#Psi_{%d}-#Psi_{R}))>", "N" };
  AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::kNone, 
				     AliFMDFlowBin::kNaive,
				     AliFMDFlowBin::kStar,
				     AliFMDFlowBin::kTdr };
  // Figure out how many things to draw
  UShort_t nm = 0;
  Short_t  sm = -1;
  for (UShort_t i = 0; i < 4; i++) { 
    if ((which & (1 << i)) != 0) {
      nm++;
      sm = i;
    }
  }
  if (what == AliFMDFlowBin::kCounts) { nm = 1; sm = 0; }

  TH1* h = 0;
  if (nm > 1) { 
    // Make 2D histogram 
    h = new TH2D(Form("%s_%s", GetName(), hNames[what]), 
		 Form("%s %s", GetTitle(), hTitles[what]), 
		 fXAxis.N(), fXAxis.Bins(), nm, 0, nm);

    // Set titles and such
    h->SetXTitle("x");
    h->SetYTitle("method");
    switch (what) { 
    case AliFMDFlowBin::kHarmonic:   
      h->SetZTitle(Form(hUnits[what], Order())); break;
    case AliFMDFlowBin::kResolution: 
      h->SetZTitle(Form(hUnits[what], Order(), PsiOrder())); break;
    default:
      h->SetZTitle(hUnits[what]); break;
    }
    h->GetYaxis()->SetNdivisions(nm+1, kFALSE);

    // Set Bin labels for the methods 
    UInt_t j = 0;
    for (UShort_t i = 0; i < 4; i++) 
      if (which & (1 << i)) h->GetYaxis()->SetBinLabel(++j, names[i]);
  }
  else {
    TString name(what == AliFMDFlowBin::kCounts ? 
		 Form("%s_%s", GetName(), hNames[what]) : 
		 Form("%s_%s_%s", GetName(), hNames[what], names[sm]));
    TString title(what == AliFMDFlowBin::kCounts ?
		  Form("%s %s", GetTitle(), hTitles[what]) : 
		  Form("%s %s %s", GetTitle(), hTitles[what], names[sm]));
    h = new TH1D(name.Data(), title.Data(), fXAxis.N(), fXAxis.Bins());
    h->SetXTitle("x");
    switch (what) { 
    case AliFMDFlowBin::kHarmonic:   
      h->SetYTitle(Form(hUnits[what], Order())); break;
    case AliFMDFlowBin::kResolution: 
      h->SetYTitle(Form(hUnits[what], Order(), PsiOrder())); break;
    default:
      h->SetYTitle(hUnits[what]); break;
    }
  }

  for (UShort_t i = 0; i < fXAxis.N(); i++) { 
    Double_t       x   = fXAxis.BinCenter(i);
    AliFMDFlowBin* bin = GetBin(x);
    Double_t       y=0, e2=0, dummy;
    if (bin->Counts() <= 0) continue;
    if (nm == 1) { 
      switch (what) { 
      case AliFMDFlowBin::kHarmonic:   
	if (bin->Correction(dummy, types[sm]) < .01) continue;
	y = bin->Value(e2, types[sm]); break;
      case AliFMDFlowBin::kResolution: 
	y = bin->Correction(e2, types[sm]); break;
      case AliFMDFlowBin::kCounts:     
	y = bin->Counts(); break;
      }
      h->SetBinContent(i+1, y);
      h->SetBinError(i+1, sqrt(e2));
      continue;
    }
    UInt_t j = 0;
    for (UShort_t k = 0; k < 4; k++)  { 
      if (!(which & (1 << k))) continue;
      switch (what) { 
      case AliFMDFlowBin::kHarmonic:   
	if (bin->Correction(dummy,types[k]) < .01) continue;
	y = bin->Value(e2, types[k]); break;
      case AliFMDFlowBin::kResolution: 
	y = bin->Correction(e2, types[k]); break;
      case AliFMDFlowBin::kCounts:     
	y = bin->Counts(); break;
      }
      h->SetBinContent(i+1, j+1, y);
      h->SetBinError(i+1, j+1, sqrt(e2));
      j++;
    }
  }
  h->SetLineColor(GetLineColor());
  h->SetLineWidth(GetLineWidth());
  h->SetLineStyle(GetLineStyle());
  h->SetFillColor(GetFillColor());
  h->SetFillStyle(GetFillStyle());
  h->SetMarkerColor(GetMarkerColor());
  h->SetMarkerStyle(GetMarkerStyle());
  h->SetMarkerSize(GetMarkerSize());
  h->SetDirectory(0);
  return h;
}
  

//____________________________________________________________________
void 
AliFMDFlowBinned1D::Draw(Option_t* option)
{
  // Draw the distribution of the harmonics or the event plane
  // resolution. 
  // Parameters: 
  //    option     String of options 
  // 
  // Options:
  //    b          Draw bare averages of cos(n(phi-Psi))
  //    n          Draw harmonic scaled by naive correction
  //    s          Draw harmonic scaled by STAR correction 
  //    t          Draw harmonic scaled by TDR correction (*)
  //    r          Draw resolution rather than harmonic 
  //    c          Draw counts rather than harmonic 
  //    h          Draw harmonics (*)
  //   
  // (*) This is the default value 
  // 
  // If more than one of b, n, s, or t is given, a 2D histogram is
  // drawn. 
  // 
  // Only one of r, c, h can be specified.  The first specified wins.
  //
  TString opt(option);
  opt.ToLower();
  TString dopt = "";
  Int_t idx = opt.Index(":");
  if (idx != kNPOS) {
    dopt = opt(idx+1, opt.Length()-idx-1);
    opt.Remove(idx, opt.Length()-idx);
  }
  UInt_t which = 0;
  if (opt.Contains("b")) which |= AliFMDFlowBin::kNone;
  if (opt.Contains("n")) which |= AliFMDFlowBin::kNaive;
  if (opt.Contains("s")) which |= AliFMDFlowBin::kStar;  
  if (opt.Contains("t")) which |= AliFMDFlowBin::kTdr;
  
  UInt_t what = AliFMDFlowBin::kHarmonic;
  if (opt.Contains("c")) what = AliFMDFlowBin::kCounts;
  if (opt.Contains("r")) what = AliFMDFlowBin::kResolution;

  TH1* h = MakeHistogram(which, what);
  h->Draw(dopt.Data());
}

  
  
//____________________________________________________________________
void 
AliFMDFlowBinned1D::Print(Option_t* option) const
{
  // Print information to standard output. 
  // Parameters: 
  //     option    String of options 
  // 
  // Options: 
  //     d         Print details 
  //     s         Print summary 
  // 
  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         oldP = std::cout.precision(3);
      std::ios_base::fmtflags oldF = 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(oldP);
      std::cout.setf(oldF, std::ios_base::floatfield);
    }
  }
  
  if (sum) { 
    UInt_t       nType = 4;
    const char*  names[] = { "Bare",    "Naive",    "STAR",    "TDR" };
    AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::kNone, 
				       AliFMDFlowBin::kNaive, 
				       AliFMDFlowBin::kStar, 
				       AliFMDFlowBin::kTdr };
    std::cout << GetName() << " - " << GetTitle() << "\n"
	      << "    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         oldP = std::cout.precision(2);
    std::ios_base::fmtflags oldF = 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(oldP);
    std::cout.setf(oldF, std::ios_base::floatfield);
  }
}

    
//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
97e94238	1	/* Copyright (C) 2007 Christian Holm Christensen <cholm@nbi.dk>
	2	*
	3	* This library is free software; you can redistribute it and/or
	4	* modify it under the terms of the GNU Lesser General Public License
	5	* as published by the Free Software Foundation; either version 2.1 of
	6	* the License, or (at your option) any later version.
	7	*
	8	* This library is distributed in the hope that it will be useful, but
	9	* WITHOUT ANY WARRANTY; without even the implied warranty of
	10	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	11	* Lesser General Public License for more details.
	12	*
	13	* You should have received a copy of the GNU Lesser General Public
	14	* License along with this library; if not, write to the Free Software
	15	* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
	16	* USA
	17	*/
39eefe19	18	/** @file
39eefe19	19	@brief Implementation of a 1-dimensional Flow "histogram" */
97e94238	20	//____________________________________________________________________
	21	//
	22	// A histogram of flow bins. The axis can by anything
	23	// (pseudo-rapidity, transvers momentum) - there's no assumption on
	24	// what is the basis of the histogram. The method Event can be used
	25	// to calculate everything in one go. Alternatively, one can use the
	26	// methods AddToEventPlane and AddToHarmonic. See also the example
	27	// TestFlow.C
39eefe19	28	#include "flow/AliFMDFlowBinned1D.h"
39eefe19	29	#include "flow/AliFMDFlowBin.h"
9b98d361	30	#include "flow/AliFMDFlowSplitter.h"
9b98d361	31	// #include <cmath>
39eefe19	32	#include <cstdlib>
	33	#include <iostream>
	34	#include <iomanip>
	35	#include <TBrowser.h>
	36	#include <TString.h>
824e71c1	37	#include <TH1.h>
824e71c1	38	#include <TH2.h>
39eefe19	39
39eefe19	40	//====================================================================
9b98d361	41	AliFMDFlowBinned1D::AliFMDFlowBinned1D()
	42	: TNamed("", ""),
	43	fXAxis(),
	44	fN(0),
	45	fBins(0),
	46	fSplitter(0)
	47	{
	48	// Default CTOR - do not use
	49	}
	50
	51	//____________________________________________________________________
	52	AliFMDFlowBinned1D::AliFMDFlowBinned1D(const char* name,
	53	const char* title,
	54	UShort_t order,
	55	UShort_t k,
39eefe19	56	UShort_t nxbins,
9b98d361	57	Double_t* xbins,
	58	AliFMDFlowSplitter* splitter)
	59	: TNamed(name,title),
	60	fXAxis(nxbins, xbins),
	61	fBins(0),
	62	fSplitter(splitter)
39eefe19	63	{
97e94238	64	// Constructor
	65	// Parameters:
	66	// Order Order
	67	// nxbins Number of bins
	68	// xbins Bin borders
39eefe19	69	UShort_t n = fXAxis.N();
9b98d361	70	fN = n;
	71	fBins = new AliFMDFlowBin*[n];
	72	for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order,k);
	73	if (!fSplitter) fSplitter = new AliFMDFlowShuffle;
39eefe19	74	}
39eefe19	75	//____________________________________________________________________
9b98d361	76	AliFMDFlowBinned1D::AliFMDFlowBinned1D(const char* name,
	77	const char* title,
	78	UShort_t order,
	79	UShort_t k,
	80	const AliFMDFlowAxis& xaxis,
	81	AliFMDFlowSplitter* splitter)
	82	: TNamed(name,title),
	83	fXAxis(xaxis),
	84	fBins(0),
	85	fSplitter(splitter)
39eefe19	86	{
97e94238	87	// Constructor
	88	// Parameters:
	89	// Order Order
	90	// xaxis X axis object
39eefe19	91	UShort_t n = fXAxis.N();
9b98d361	92	fN = n;
	93	fBins = new AliFMDFlowBin*[n];
	94	for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(order,k);
	95	if (!fSplitter) fSplitter = new AliFMDFlowShuffle;
	96
39eefe19	97	}
	98	//____________________________________________________________________
	99	AliFMDFlowBinned1D::AliFMDFlowBinned1D(const AliFMDFlowBinned1D& o)
9b98d361	100	: TNamed(o),
	101	TAttLine(o),
	102	TAttFill(o),
	103	TAttMarker(o),
	104	fXAxis(o.fXAxis),
	105	fSplitter(0)
39eefe19	106	{
97e94238	107	// Copy constructor
	108	// Parameters:
	109	// o Object to copy from
39eefe19	110	UShort_t n = fXAxis.N();
9b98d361	111	fSplitter = new AliFMDFlowShuffle;
	112	fN = n;
	113	fBins = new AliFMDFlowBin*[n];
39eefe19	114	for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
	115	}
	116	//____________________________________________________________________
	117	AliFMDFlowBinned1D&
	118	AliFMDFlowBinned1D::operator=(const AliFMDFlowBinned1D& o)
	119	{
97e94238	120	// Assignment operator
	121	// Parameters:
	122	// o Object to assign from
	123	//
	124	// Returns reference to this object
9b98d361	125	SetLineColor(o.GetLineColor());
	126	SetLineStyle(o.GetLineStyle());
	127	SetLineWidth(o.GetLineWidth());
	128	SetFillColor(o.GetFillColor());
	129	SetFillStyle(o.GetFillStyle());
	130	SetMarkerColor(o.GetMarkerColor());
	131	SetMarkerStyle(o.GetMarkerStyle());
	132	SetMarkerSize(o.GetMarkerSize());
	133	this->SetName(o.GetName());
	134	this->SetTitle(o.GetTitle());
39eefe19	135	if (fBins) {
	136	for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
	137	delete [] fBins;
	138	}
9b98d361	139	if (fSplitter) delete fSplitter;
39eefe19	140	fXAxis = o.fXAxis;
39eefe19	141	UShort_t n = fXAxis.N();
9b98d361	142	fN = n;
	143	fSplitter = new AliFMDFlowShuffle;
	144	fBins = new AliFMDFlowBin*[n];
39eefe19	145	for (UInt_t i = 0; i < n; i++) fBins[i]= new AliFMDFlowBin(*(o.fBins[i]));
	146	return *this;
	147	}
	148
	149	//____________________________________________________________________
	150	AliFMDFlowBinned1D::~AliFMDFlowBinned1D()
	151	{
97e94238	152	// Destructor
	153	// Parameters:
	154	// none
39eefe19	155	if (fBins) {
	156	for (UInt_t i = 0; i < fXAxis.N(); i++) delete fBins[i];
	157	delete [] fBins;
	158	}
9b98d361	159	if (fSplitter) delete fSplitter;
39eefe19	160	}
	161
	162	//____________________________________________________________________
	163	AliFMDFlowBin*
	164	AliFMDFlowBinned1D::GetBin(UShort_t i) const
	165	{
97e94238	166	// Get the ith bin
	167	// Parameters:
	168	// i Bin number
	169	//
	170	// Return pointer to bin, or null.
39eefe19	171	if (i >= fXAxis.N()) return 0;
	172	return fBins[i];
	173	}
	174	//____________________________________________________________________
	175	AliFMDFlowBin*
	176	AliFMDFlowBinned1D::GetBin(Double_t x) const
	177	{
97e94238	178	// Get the bin that contains x
	179	// Parameters:
	180	// x X axis value
	181	//
	182	// Return pointer to bin, or null.
39eefe19	183	Int_t i = fXAxis.FindBin(x);
	184	if (i < 0) return 0;
	185	UShort_t j = i;
	186	return GetBin(j);
	187	}
	188
	189	//____________________________________________________________________
	190	void
	191	AliFMDFlowBinned1D::Begin()
	192	{
97e94238	193	// Called at the beginning of an event
	194	// Parameters:
	195	// none
39eefe19	196	for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Begin();
9b98d361	197	fSplitter->Begin();
39eefe19	198	}
	199	//____________________________________________________________________
	200	void
	201	AliFMDFlowBinned1D::End()
	202	{
97e94238	203	// Called at the end of an event
	204	// Parameters:
	205	// none
39eefe19	206	for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->End();
9b98d361	207	fSplitter->End();
39eefe19	208	}
	209	//____________________________________________________________________
	210	void
	211	AliFMDFlowBinned1D::Finish()
	212	{
97e94238	213	// Called at the end of an job
	214	// Parameters:
	215	// none
39eefe19	216	for (UInt_t i = 0; i < fXAxis.N(); i++) fBins[i]->Finish();
	217	}
	218	//____________________________________________________________________
	219	Bool_t
97e94238	220	AliFMDFlowBinned1D::AddToEventPlane(Double_t x, Double_t phi,
97e94238	221	Double_t w, Bool_t a)
39eefe19	222	{
97e94238	223	// Called to add a contribution to the event plane
	224	// Parameters:
	225	// x Bin to fill into
	226	// w Weight
	227	// phi The angle phi in radians
	228	// a If true, add to sub-event A, otherwise sub-event B
	229	//
	230	// Return false if x falls outside the defined range, true otherwise
39eefe19	231	AliFMDFlowBin* bin = GetBin(x);
	232	if (!bin) return kFALSE;
	233	bin->AddToEventPlane(phi, w, a);
	234	return kTRUE;
	235	}
	236
	237	//____________________________________________________________________
	238	Bool_t
9b98d361	239	AliFMDFlowBinned1D::AddToHarmonic(Double_t x, Double_t phi,
9b98d361	240	Double_t wp, Double_t wh)
39eefe19	241	{
97e94238	242	// Called to add a contribution to the harmonic
97e94238	243	// Parameters:
9b98d361	244	// x Bin to fill into
	245	// phi The angle phi in radians
	246	// wp weight of event plane
	247	// wh weight of harmonic
	248
97e94238	249	// Return false if x falls outside the defined range, true otherwise
39eefe19	250	AliFMDFlowBin* bin = GetBin(x);
39eefe19	251	if (!bin) return kFALSE;
9b98d361	252	bin->AddToHarmonic(phi, wp, wh);
39eefe19	253	return kTRUE;
	254	}
	255
9b98d361	256	//____________________________________________________________________
	257	UShort_t
	258	AliFMDFlowBinned1D::Order() const
	259	{
	260	return GetBin(UShort_t(0))->Order();
	261	}
	262	//____________________________________________________________________
	263	UShort_t
	264	AliFMDFlowBinned1D::PsiOrder() const
	265	{
	266	return GetBin(UShort_t(0))->Order();
	267	}
	268
39eefe19	269	//____________________________________________________________________
39eefe19	270	void
9b98d361	271	AliFMDFlowBinned1D::Event(ULong_t n, Double_t* phis, Double_t* xs,
9b98d361	272	Double_t* wp, Double_t* wh)
39eefe19	273	{
97e94238	274	// Process a full event.
	275	// Parameters:
	276	// phis List of n phi=[0,2pi] angles
	277	// xs List of n x values.
	278	// ws Weights
	279	// n Size of phis and xs
39eefe19	280	Begin();
9b98d361	281	fSplitter->Event(phis, xs, n);
	282	for (ULong_t i = 0; i < n; i++)
	283	AddToEventPlane(xs[i], phis[i], (wp ? wp[i] : 1), fSplitter->Select(i));
	284	for (ULong_t i = 0; i < n; i++)
	285	AddToHarmonic(xs[i], phis[i], (wp ? wp[i] : 1), (wh ? wh[i] : 1));
39eefe19	286	End();
	287	}
	288
	289	//____________________________________________________________________
	290	void
	291	AliFMDFlowBinned1D::Browse(TBrowser* b)
	292	{
97e94238	293	// Browse this object.
39eefe19	294	b->Add(&fXAxis, "xaxis");
9b98d361	295	if (fSplitter) b->Add(fSplitter, "Splitter");
39eefe19	296	for (UInt_t i = 0; i < fXAxis.N(); i++)
	297	b->Add(fBins[i], Form("bin_%03d", i));
	298	}
	299
824e71c1	300	//____________________________________________________________________
9b98d361	301	TH1*
9b98d361	302	AliFMDFlowBinned1D::MakeHistogram(UInt_t which, UInt_t what)
824e71c1	303	{
9b98d361	304	// Make a histogram of some of the stuff in this object
	305
	306	// Some strings
	307	const char* names[] = { "Bare", "Naive", "STAR", "TDR" };
	308	const char* hNames[] = { "flow", "resolution", "counts" };
	309	const char* hTitles[] = { "Flow", "Resolution", "Counts" };
	310	const char* hUnits[] = { "v_{%d}","<cos(%d(#Psi_{%d}-#Psi_{R}))>", "N" };
	311	AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::kNone,
	312	AliFMDFlowBin::kNaive,
	313	AliFMDFlowBin::kStar,
	314	AliFMDFlowBin::kTdr };
	315	// Figure out how many things to draw
824e71c1	316	UShort_t nm = 0;
824e71c1	317	Short_t sm = -1;
9b98d361	318	for (UShort_t i = 0; i < 4; i++) {
	319	if ((which & (1 << i)) != 0) {
	320	nm++;
	321	sm = i;
	322	}
	323	}
	324	if (what == AliFMDFlowBin::kCounts) { nm = 1; sm = 0; }
	325
824e71c1	326	TH1* h = 0;
824e71c1	327	if (nm > 1) {
9b98d361	328	// Make 2D histogram
	329	h = new TH2D(Form("%s_%s", GetName(), hNames[what]),
	330	Form("%s %s", GetTitle(), hTitles[what]),
824e71c1	331	fXAxis.N(), fXAxis.Bins(), nm, 0, nm);
9b98d361	332
9b98d361	333	// Set titles and such
824e71c1	334	h->SetXTitle("x");
824e71c1	335	h->SetYTitle("method");
9b98d361	336	switch (what) {
	337	case AliFMDFlowBin::kHarmonic:
	338	h->SetZTitle(Form(hUnits[what], Order())); break;
	339	case AliFMDFlowBin::kResolution:
	340	h->SetZTitle(Form(hUnits[what], Order(), PsiOrder())); break;
	341	default:
	342	h->SetZTitle(hUnits[what]); break;
	343	}
824e71c1	344	h->GetYaxis()->SetNdivisions(nm+1, kFALSE);
9b98d361	345
9b98d361	346	// Set Bin labels for the methods
824e71c1	347	UInt_t j = 0;
824e71c1	348	for (UShort_t i = 0; i < 4; i++)
9b98d361	349	if (which & (1 << i)) h->GetYaxis()->SetBinLabel(++j, names[i]);
824e71c1	350	}
824e71c1	351	else {
9b98d361	352	TString name(what == AliFMDFlowBin::kCounts ?
	353	Form("%s_%s", GetName(), hNames[what]) :
	354	Form("%s_%s_%s", GetName(), hNames[what], names[sm]));
	355	TString title(what == AliFMDFlowBin::kCounts ?
	356	Form("%s %s", GetTitle(), hTitles[what]) :
	357	Form("%s %s %s", GetTitle(), hTitles[what], names[sm]));
	358	h = new TH1D(name.Data(), title.Data(), fXAxis.N(), fXAxis.Bins());
824e71c1	359	h->SetXTitle("x");
9b98d361	360	switch (what) {
	361	case AliFMDFlowBin::kHarmonic:
	362	h->SetYTitle(Form(hUnits[what], Order())); break;
	363	case AliFMDFlowBin::kResolution:
	364	h->SetYTitle(Form(hUnits[what], Order(), PsiOrder())); break;
	365	default:
	366	h->SetYTitle(hUnits[what]); break;
	367	}
824e71c1	368	}
	369
	370	for (UShort_t i = 0; i < fXAxis.N(); i++) {
	371	Double_t x = fXAxis.BinCenter(i);
	372	AliFMDFlowBin* bin = GetBin(x);
9b98d361	373	Double_t y=0, e2=0, dummy;
9b98d361	374	if (bin->Counts() <= 0) continue;
824e71c1	375	if (nm == 1) {
9b98d361	376	switch (what) {
	377	case AliFMDFlowBin::kHarmonic:
	378	if (bin->Correction(dummy, types[sm]) < .01) continue;
	379	y = bin->Value(e2, types[sm]); break;
	380	case AliFMDFlowBin::kResolution:
	381	y = bin->Correction(e2, types[sm]); break;
	382	case AliFMDFlowBin::kCounts:
	383	y = bin->Counts(); break;
	384	}
	385	h->SetBinContent(i+1, y);
824e71c1	386	h->SetBinError(i+1, sqrt(e2));
	387	continue;
	388	}
	389	UInt_t j = 0;
	390	for (UShort_t k = 0; k < 4; k++) {
9b98d361	391	if (!(which & (1 << k))) continue;
	392	switch (what) {
	393	case AliFMDFlowBin::kHarmonic:
	394	if (bin->Correction(dummy,types[k]) < .01) continue;
	395	y = bin->Value(e2, types[k]); break;
	396	case AliFMDFlowBin::kResolution:
	397	y = bin->Correction(e2, types[k]); break;
	398	case AliFMDFlowBin::kCounts:
	399	y = bin->Counts(); break;
	400	}
	401	h->SetBinContent(i+1, j+1, y);
824e71c1	402	h->SetBinError(i+1, j+1, sqrt(e2));
	403	j++;
	404	}
	405	}
9b98d361	406	h->SetLineColor(GetLineColor());
	407	h->SetLineWidth(GetLineWidth());
	408	h->SetLineStyle(GetLineStyle());
	409	h->SetFillColor(GetFillColor());
	410	h->SetFillStyle(GetFillStyle());
	411	h->SetMarkerColor(GetMarkerColor());
	412	h->SetMarkerStyle(GetMarkerStyle());
	413	h->SetMarkerSize(GetMarkerSize());
	414	h->SetDirectory(0);
	415	return h;
	416	}
	417
	418
	419	//____________________________________________________________________
	420	void
	421	AliFMDFlowBinned1D::Draw(Option_t* option)
	422	{
	423	// Draw the distribution of the harmonics or the event plane
	424	// resolution.
	425	// Parameters:
	426	// option String of options
	427	//
	428	// Options:
	429	// b Draw bare averages of cos(n(phi-Psi))
	430	// n Draw harmonic scaled by naive correction
	431	// s Draw harmonic scaled by STAR correction
	432	// t Draw harmonic scaled by TDR correction (*)
	433	// r Draw resolution rather than harmonic
	434	// c Draw counts rather than harmonic
	435	// h Draw harmonics (*)
	436	//
	437	// (*) This is the default value
	438	//
	439	// If more than one of b, n, s, or t is given, a 2D histogram is
	440	// drawn.
	441	//
	442	// Only one of r, c, h can be specified. The first specified wins.
	443	//
	444	TString opt(option);
	445	opt.ToLower();
	446	TString dopt = "";
	447	Int_t idx = opt.Index(":");
	448	if (idx != kNPOS) {
	449	dopt = opt(idx+1, opt.Length()-idx-1);
	450	opt.Remove(idx, opt.Length()-idx);
	451	}
	452	UInt_t which = 0;
	453	if (opt.Contains("b")) which \|= AliFMDFlowBin::kNone;
	454	if (opt.Contains("n")) which \|= AliFMDFlowBin::kNaive;
	455	if (opt.Contains("s")) which \|= AliFMDFlowBin::kStar;
	456	if (opt.Contains("t")) which \|= AliFMDFlowBin::kTdr;
	457
	458	UInt_t what = AliFMDFlowBin::kHarmonic;
	459	if (opt.Contains("c")) what = AliFMDFlowBin::kCounts;
	460	if (opt.Contains("r")) what = AliFMDFlowBin::kResolution;
	461
	462	TH1* h = MakeHistogram(which, what);
	463	h->Draw(dopt.Data());
824e71c1	464	}
	465
	466
	467
39eefe19	468	//____________________________________________________________________
	469	void
	470	AliFMDFlowBinned1D::Print(Option_t* option) const
	471	{
97e94238	472	// Print information to standard output.
	473	// Parameters:
	474	// option String of options
	475	//
	476	// Options:
	477	// d Print details
	478	// s Print summary
	479	//
39eefe19	480	TString opt(option);
	481	opt.ToLower();
	482	Bool_t det = opt.Contains("d");
	483	Bool_t sum = opt.Contains("s");
	484	if (det) {
	485	for (UShort_t i = 0; i < fXAxis.N(); i++) {
	486	Double_t x = fXAxis.BinCenter(i);
97e94238	487	std::streamsize oldP = std::cout.precision(3);
97e94238	488	std::ios_base::fmtflags oldF = std::cout.setf(std::ios_base::fixed,
39eefe19	489	std::ios_base::floatfield);
	490	std::cout << "x=" << std::setw(5) << x << std::endl;
	491	fBins[i]->Print();
97e94238	492	std::cout.precision(oldP);
97e94238	493	std::cout.setf(oldF, std::ios_base::floatfield);
39eefe19	494	}
	495	}
	496
	497	if (sum) {
	498	UInt_t nType = 4;
	499	const char* names[] = { "Bare", "Naive", "STAR", "TDR" };
97e94238	500	AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::kNone,
	501	AliFMDFlowBin::kNaive,
	502	AliFMDFlowBin::kStar,
	503	AliFMDFlowBin::kTdr };
9b98d361	504	std::cout << GetName() << " - " << GetTitle() << "\n"
9b98d361	505	<< " x";
39eefe19	506	for (UInt_t i = 0; i < nType; i++)
	507	std::cout << " \| " << std::setw(6+6+5) << names[i];
	508	std::cout << "\n-----" << std::setfill('-');
	509	for (UInt_t i = 0; i < nType; i++)
	510	std::cout << "-+-" << std::setw(6+6+5) << "-";
	511	std::cout << std::setfill(' ') << std::endl;
	512
97e94238	513	std::streamsize oldP = std::cout.precision(2);
	514	std::ios_base::fmtflags oldF = std::cout.setf(std::ios_base::fixed,
	515	std::ios_base::floatfield);
39eefe19	516	for (UShort_t i = 0; i < fXAxis.N(); i++) {
	517	Double_t x = fXAxis.BinCenter(i);
	518	std::cout << std::setprecision(2) << std::setw(5) << x << std::flush;
	519	for (UShort_t j = 0; j < nType; j++) {
	520	Double_t e2v;
	521	Double_t v = fBins[i]->Value(e2v, types[j]);
	522	std::cout << std::setprecision(3) << " \| "
	523	<< std::setw(6) << 100 * v << " +/- "
	524	<< std::setw(6) << 100 * sqrt(e2v);
	525	}
	526	std::cout << std::endl;
	527	}
97e94238	528	std::cout.precision(oldP);
97e94238	529	std::cout.setf(oldF, std::ios_base::floatfield);
39eefe19	530	}
	531	}
	532
	533
	534	//____________________________________________________________________
	535	//
	536	// EOF
	537	//