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