@brief Implementation of an Axis in a Flow "histogram" */
#include "flow/AliFMDFlowAxis.h"
#include <cmath>
+#include <iostream>
+#include <iomanip>
//====================================================================
AliFMDFlowAxis::AliFMDFlowAxis(UShort_t n, Double_t* bins)
return fBins[i+1];
}
+//____________________________________________________________________
+void
+AliFMDFlowAxis::Print(Option_t*) const
+{
+ std::ios_base::fmtflags old_flags = std::cout.setf(std::ios_base::fixed,
+ std::ios_base::floatfield);
+ for (UShort_t i = 0; i < fN; i++)
+ std::cout << std::setw(5) << BinLower(i) << " - "
+ << std::setw(5) << BinUpper(i) << std::endl;
+ std::cout.setf(old_flags, std::ios_base::floatfield);
+}
+
+
//____________________________________________________________________
//
// EOF
Double_t* Bins() const { return fBins; }
/** Get the number of bins */
UShort_t N() const { return fN; }
+ /** Print the axis */
+ void Print(Option_t* option="") const; //*MENU*
protected:
/** Number of bins */
UShort_t fN;
#include <cmath>
#include <iostream>
#include <iomanip>
+#include <TBrowser.h>
//====================================================================
void
AliFMDFlowBin::Finish()
{}
+//____________________________________________________________________
+void
+AliFMDFlowBin::Browse(TBrowser* b)
+{
+ b->Add(&fPsi, "Full event plane");
+ b->Add(&fPsiA, "Sub-event A event plane");
+ b->Add(&fPsiB, "Sub-event A event plane");
+ b->Add(&fRes, "Naive resolution");
+ b->Add(&fResStar, "STAR resolution");
+ b->Add(&fResTdr, "TDR resolution");
+ b->Add(&fHarmonic, "Harmonic");
+}
+
//____________________________________________________________________
void
AliFMDFlowBin::Print(Option_t*) const
#include <flow/AliFMDFlowResolution.h>
#include <TObject.h>
+//Forward declaration
+class TBrowser;
+
/** @defgroup c_binned Binned flow
@brief This group contains code for binned flow analysis. Two
kinds of "histograms" are defined - a 1 dimensional and a 2
@return the value of the harmonic */
virtual Double_t Correction(Double_t& e2, CorType t=naive) const;
/** Print summary to standard output */
- virtual void Print(Option_t* option="") const;
-
+ virtual void Print(Option_t* option="") const; //*MENU*
+ /** Return true */
+ virtual Bool_t IsFolder() const { return kTRUE; }
+ /** Browse this item */
+ virtual void Browse(TBrowser* b);
/** Get the event plane angle */
virtual Double_t Psi() const { return fPsi.Psi(); }
/** Get the sub-event A plane angle */
#include <iomanip>
#include <TBrowser.h>
#include <TString.h>
+#include <TH1.h>
+#include <TH2.h>
//====================================================================
AliFMDFlowBinned1D::AliFMDFlowBinned1D(UShort_t order,
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
@return The flow object in bin @a i or 0 if out of range */
virtual AliFMDFlowBin* GetBin(UShort_t i) const;
/** Print to standard out */
- virtual void Print(Option_t* option="s") const;
+ virtual void Print(Option_t* option="s") const; //*MENU*
+ /** Draw as a histogram
+ @param option Option string.
+ - s Draw STAR method.
+ - t Draw TDR method
+ - n Draw Naive method
+ - b Draw bare method
+ - r Draw resolution rather than harmonic.
+ */
+ virtual void Draw(Option_t* option="stnb"); //*MENU*
/** Whether this is to be considered a folder */
Bool_t IsFolder() const { return kTRUE; }
/** Browse this object */
#include "flow/AliFMDFlowResolution.h"
#include "flow/AliFMDFlowUtil.h"
#include "flow/AliFMDFlowBessel.h"
+#include <TGraphErrors.h>
+#include <TH1.h>
+#include <TH2.h>
#include <iostream>
//#include <cmath>
return sqrt(2) * sqrt(fabs(fAverage));
}
+//____________________________________________________________________
+void
+AliFMDFlowResolution::Draw(Option_t* option)
+{
+ TGraph* g = new TGraph(100);
+ for (UShort_t i = 0; i < g->GetN(); i++) {
+ Double_t x = -1. + 2. / 100 * i;
+ Double_t y = sqrt(2) * sqrt(fabs(x));
+ g->SetPoint(i, x, y);
+ }
+ g->SetName("naive_res");
+ g->SetTitle("Naive Resolution Function");
+ g->GetHistogram()->SetXTitle("<cos(n(|#Psi_{A}-#Psi_{B}|))>");
+ g->GetHistogram()->SetYTitle("<cos(n(#Psi-#Psi_{R}))>");
+ g->GetHistogram()->SetDirectory(0);
+ g->Draw(Form("lh %s", option));
+}
+
+
//====================================================================
Double_t
AliFMDFlowResolutionStar::Correction(UShort_t k, Double_t& e2) const
Double_t chi = Chi(fAverage, k, delta);
Double_t dr = 0;
Double_t res = Res(sqrt(2) * chi, k, dr);
- e2 = pow(dr * delta,2);
+ e2 = pow(dr * delta,2);
return res;
}
//____________________________________________________________________
AliFMDFlowResolutionStar::Chi(Double_t res, UShort_t k,
Double_t& delta) const
{
- delta = 1;
+ delta = 1;
Double_t chi = 2;
Double_t dr = 0;
for (UInt_t i = 0; i < 15; i++) {
}
//____________________________________________________________________
Double_t
-AliFMDFlowResolutionStar::Res(Double_t chi, UShort_t k,
- Double_t& dr) const
+AliFMDFlowResolutionStar::Res(Double_t chi, UShort_t k, Double_t& dr) const
{
// The resolution function is
//
return r;
}
+//____________________________________________________________________
+void
+AliFMDFlowResolutionStar::Draw(Option_t* option)
+{
+ TString opt(option);
+ opt.ToLower();
+ Bool_t chi = opt.Contains("chi");
+
+ TH2* h = new TH2D(Form("star_%s_frame", (chi ? "chi" : "res")),
+ Form("STAR %s Function for k=(1,2,4)",
+ (chi ? "Chi" : "Resolution")),
+ 100, 0, 1, 100, 0, (chi ? 3 : 1.5));
+ h->SetXTitle("<cos(n(|#Psi_{A}-#Psi_{B}|))>");
+ h->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>"));
+ h->SetStats(0);
+ h->SetDirectory(0);
+ h->Draw();
+
+ for (UShort_t k = 1; k <= 4; k++) {
+ if (k == 3) continue;
+
+ TGraphErrors* g = new TGraphErrors(100);
+ for (UShort_t i = 0; i < g->GetN(); i++) {
+ Double_t e2 = 0;
+ Double_t x = 1. / 100 * i;
+ Double_t y = 0;
+ Double_t c = Chi(x, k, e2);
+ if (chi) y = c;
+ else {
+ Double_t dr = 0;
+ y = Res(sqrt(2) * c, k, dr);
+ e2 = pow(dr * e2,2);
+ }
+ g->SetPoint(i, x, y);
+ g->SetPointError(i, 0, sqrt(e2));
+ }
+ g->SetLineColor(k);
+ g->SetName(Form("star_%s_k%d", (chi ? "chi" : "res"), k));
+ g->SetTitle(Form("STAR %s Function for k=%d",
+ (chi ? "Chi" : "Resolution"), k));
+ g->GetHistogram()->SetXTitle("<cos(n(|#Psi_{A}-#Psi_{B}|))>");
+ g->GetHistogram()->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>"));
+ g->GetHistogram()->SetDirectory(0);
+ g->Draw(Form("l %s same", option));
+ }
+}
+
//====================================================================
void
-AliFMDFlowResolutionTDR::Clear()
+AliFMDFlowResolutionTDR::Clear(Option_t*)
{
fN = 0;
fLarge = 0;
//
// where z = chi^2 / 2
//
+ if (fLarge == 0) {
+ std::cerr << "TDR: K = 0" << std::endl;
+ return -1;
+ }
+ if (fN == 0) {
+ std::cerr << "TDR: N = 0" << std::endl;
+ return -1;
+ }
+ Double_t r = Double_t(fLarge) / fN;
Double_t echi2 = 0;
- Double_t y = Chi2Over2(echi2);
+ Double_t y = Chi2Over2(r, echi2);
+ return Res(k, y, echi2, e2);
+}
+
+
+
+//____________________________________________________________________
+Double_t
+AliFMDFlowResolutionTDR::Res(UShort_t k, Double_t y, Double_t echi2,
+ Double_t& e2) const
+{
+ // y = chi^2 / 2
Double_t chi = sqrt(2 * y);
Double_t c = sqrt(M_PI) * exp(-y) / 2;
-
+
// i[0] = I[(k-1)/2-1], i[1] = I[(k-1)/2], i[2] = I[(k-1)/2+1]
Double_t i[3], di[3];
AliFMDFlowBessel::Inu(Double_t(k-3)/2, Double_t(k+1)/2, y, i, di);
e2 = dr * dr * echi2;
return r;
}
+
//____________________________________________________________________
Double_t
AliFMDFlowResolutionTDR::Correction(UShort_t k) const
}
//____________________________________________________________________
Double_t
-AliFMDFlowResolutionTDR::Chi2Over2(Double_t& e2) const
+AliFMDFlowResolutionTDR::Chi2Over2(Double_t r, Double_t& e2) const
{
// From nucl-ex/9711003v2
//
// 1 - r r - 1
// = - -------------- = ------------
// 4 r N log(2 r) 4 k log(2 r)
-
- if (fLarge == 0) {
- std::cerr << "TDR: Large = 0" << std::endl;
- return -1;
+ if (r == 0) {
+ std::cerr << "TDR: Large/All = " << r << " <= 0!" << std::endl;
+ return 0;
}
- Double_t ratio = Double_t(fN) / (2 * fLarge);
+ Double_t ratio = 1. / (2*r); // Double_t(fN) / (2 * fLarge);
if (ratio <= 0) {
std::cerr << "TDR: Large/All = " << ratio << " <= 0!" << std::endl;
return -1;
<< " < 0" << std::endl;
return -1;
}
- Double_t r = Double_t(fLarge) / fN;
- e2 = (r - 1) / (4 * fLarge * log(2 * r));
+ if (fLarge != 0) e2 = (r - 1) / (4 * fLarge * log(2 * r));
+ else e2 = (r - 1) / (4 * r * log(2 * r));
return chi2over2;
}
+//____________________________________________________________________
+void
+AliFMDFlowResolutionTDR::Draw(Option_t* option)
+{
+ TString opt(option);
+ opt.ToLower();
+ Bool_t chi = opt.Contains("chi");
+
+ TH2* h = new TH2D(Form("tdr_%s_frame", (chi ? "chi" : "res")),
+ Form("TDR %s Function for k=(1,2,4)",
+ (chi ? "Chi" : "Resolution")),
+ 100, 0, 1, 100, 0, (chi ? 3 : 1.5));
+ h->SetXTitle("K/N");
+ h->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>"));
+ h->SetStats(0);
+ h->Draw();
+ h->SetDirectory(0);
+
+ for (UShort_t k = 1; k <= 4; k++) {
+ if (k == 3) continue;
+
+ TGraphErrors* g = new TGraphErrors;
+ Int_t i = 0;
+ for (Double_t x = 0.02; x < 0.5; x += 0.01) {
+ Double_t e2 = 0;
+ Double_t y = 0;
+ Double_t c = Chi2Over2(x, e2);
+ if (chi) y = sqrt(2 * c);
+ else {
+ Double_t dr = 0;
+ y = Res(k, c, e2, dr);
+ e2 = dr;
+ }
+ g->SetPoint(i, x, y);
+ g->SetPointError(i, 0, sqrt(fabs(e2)));
+ i++;
+ }
+ g->SetLineColor(k);
+ g->SetName(Form("tdr_%s_k%d", (chi ? "chi" : "res"), k));
+ g->SetTitle(Form("TDR %s Function for k=%d",
+ (chi ? "Chi" : "Resolution"), k));
+ g->GetHistogram()->SetXTitle("K/N");
+ g->GetHistogram()->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>"));
+ g->GetHistogram()->SetDirectory(0);
+ g->Draw(Form("l %s %s", (k == 0 ? "h" : "same"), option));
+ if (chi) break;
+ }
+}
//____________________________________________________________________
//
virtual Double_t Correction(UShort_t k) const;
/** Get the harmnic order */
UShort_t Order() const { return fOrder; }
+ /** Draw this corrrection function
+ @param option Options passed to drawing */
+ virtual void Draw(Option_t* option=""); //*MENU*
protected:
/** Order */
UShort_t fOrder;
which is taken to be @f$ \delta\chi@f$
@return @f$\chi@f$ */
virtual Double_t Chi(Double_t res, UShort_t k, Double_t& delta) const;
+ /** Draw this corrrection function
+ @param option Options passed to drawing */
+ virtual void Draw(Option_t* option=""); //*MENU*
protected:
/** Calculate resolution
@param chi @f$ \chi@f$
: AliFMDFlowResolution(n), fLarge(0) {}
/** DTOR */
~AliFMDFlowResolutionTDR() {}
- virtual void Clear();
+ virtual void Clear(Option_t* option="");
/** add a data point */
virtual void Add(Double_t psi_a, Double_t psi_b);
/** Get the correction for harmonic strength of order @a k
/** Get @f$ \chi^2/2@f$
@param e2 The square error on the correction
@return @f$ \chi^2/2@f$ */
- virtual Double_t Chi2Over2(Double_t& e2) const;
+ virtual Double_t Chi2Over2(Double_t r, Double_t& e2) const;
+ /** Draw this corrrection function
+ @param option Options passed to drawing */
+ virtual void Draw(Option_t* option=""); //*MENU*
protected:
+ /** Calculate resolution
+ @param k Order factor
+ @param y @f$ \chi^2/2@f$
+ @param echi @f$\delta\chi@f$
+ @param dr On return, the derivative of @f$ R(\chi)@f$
+ @return
+ @f[
+ \frac{\sqrt{\pi/2}}{2}\chi e^{-\chi^2/2}
+ (I_{\frac{(k-1)}{2}}(\chi^2/2)+ I_{\frac{(k+1)}{2}}(\chi^2/2))
+ @f] */
+ Double_t Res(UShort_t k, Double_t y, Double_t echi2, Double_t& e2) const;
/** Number of events with large diviation */
ULong_t fLarge;
/** Define for ROOT I/O */
\end{array}\right.
@f]
*/
-class AliFMDFlowStat
+class AliFMDFlowStat : public TObject
{
public:
/** Constructor */
/** Destructor */
virtual ~AliFMDFlowStat() {}
/** Reset */
- void Clear() { fAverage = fSqVar = 0; fN = 0; }
+ void Clear(Option_t* option="");
/** Add another data point */
void Add(Double_t x);
/** Get the average */
ClassDef(AliFMDFlowStat,1);
};
+//__________________________________________________________________
+inline void
+AliFMDFlowStat::Clear(Option_t*)
+{
+ fAverage = fSqVar = 0;
+ fN = 0;
+}
+
//__________________________________________________________________
inline void
AliFMDFlowStat::Add(Double_t x)
@return the value of the harmonic */
Double_t Correction(Double_t& e2, CorType t=naive) const;
/** Print to standard out. */
- void Print(Option_t* option="") const;
+ void Print(Option_t* option="s") const;
/** The well-known event plane */
Double_t fPsiR;
/** True resolution */
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff