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

/** @file 
    @brief Implementation of an Resolution class */
#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>

//====================================================================
void 
AliFMDFlowResolution::Add(Double_t psiA, Double_t psiB) 
{ 
  Double_t diff    = NormalizeAngle(fOrder * (psiA - psiB));
  Double_t contrib = cos(diff);
  AliFMDFlowStat::Add(contrib);
}

//____________________________________________________________________
Double_t 
AliFMDFlowResolution::Correction(UShort_t k) const 
{ 
  Double_t e;
  return Correction(k, e); 
}

//____________________________________________________________________
Double_t 
AliFMDFlowResolution::Correction(UShort_t, Double_t& e2) const 
{ 
  e2 = fSqVar / fN;
  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 
{ 
  if (k > 4) return 0;
  Double_t delta = 0;
  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);
  return res;
}
//____________________________________________________________________
Double_t 
AliFMDFlowResolutionStar::Correction(UShort_t k) const 
{ 
  Double_t e;
  return Correction(k, e); 
}
//____________________________________________________________________
Double_t 
AliFMDFlowResolutionStar::Chi(Double_t res, UShort_t k, 
			      Double_t& delta) const 
{
  delta          = 1;
  Double_t chi   = 2;
  Double_t dr    = 0;
  for (UInt_t i = 0; i < 15; i++) { 
    if (Res(chi, k, dr) < res) chi += delta;
    else                       chi -= delta;
    delta /= 2;
  }
  return chi;
}
//____________________________________________________________________
Double_t 
AliFMDFlowResolutionStar::Res(Double_t chi, UShort_t k, Double_t& dr) const 
{ 
  // The resolution function is 
  // 
  //          sqrt(pi) x exp(-x/4) (f1(x^2/4) + f2(x^2/4))
  //   r(x) = --------------------------------------------
  //                          2 sqrt(2) 
  // 
  //        
  //        = c x (f1(y) - f2(y))
  //
  // where f1 is the modified Bessel function first kind I_{(k-1)/2}, 
  // and f2 is the modified Bessel function of the first kind
  // I_{(k+1)/2}, and 
  // 
  //          sqrt(pi) exp(-x^2/4) 
  //      c = --------------------,   y = x^2/4
  //              2 sqrt(2)
  // 
  // The derivative of the resolution function is 
  //
  //            c 
  //    r'(y) = - (4 sqrt(y) (f1'(y) - f2'(y)) - (4 y - 2)(f1(y) - f2(y)))
  //            2
  // 
  //            c                                    r(y)   
  //          = - (4 sqrt(y) (f1'(y) - f2'(y))) + --------- - sqrt(y) r(y)
  //		2             			  2 sqrt(y)       
  // 
  // Since dI_n(x)/dx = I_(n-1)(x) - n / x I_n(x), and substituting 
  // f3(y) = I_((k-3)/2)(y) 
  // 
  //            c  
  //    r'(y) = - ((4 - 2 k) f1(y) - (4 y + 2 k) f2(y) + 4 y f3(y))
  //            2   
  // 
  Double_t y  = chi * chi / 4;
  Double_t c  = sqrt(M_PI) * exp(-y) / (2 * sqrt(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);
  
  Double_t r  = c * chi * (i[2] + i[1]);
  dr        = c / 2 * ((4 - 2*k)*i[1] - (4*y + 2*k)*i[2] + 4*y*i[0]);

  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(Option_t*) 
{
  fN = 0;
  fLarge = 0;
}
//____________________________________________________________________
void 
AliFMDFlowResolutionTDR::Add(Double_t psi_a, Double_t psi_b)
{ 
  Double_t a = fabs(psi_a - psi_b);
  if (a >= .5 * M_PI) fLarge++;
  fN++;
}
//____________________________________________________________________
Double_t 
AliFMDFlowResolutionTDR::Correction(UShort_t k, Double_t& e2) const 
{ 
  // From nucl-ex/9711003v2 
  //
  //
  //   <cos n Delta phi> = 
  //
  //         sqrt(pi)
  //         -------- chi exp(-z) (I_((n-1)/2)(z) + I_((n+1)/2)(z))
  //            2
  // 
  // 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(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);
  
  Double_t r  = c * chi * (i[2] + i[1]);
  Double_t dr = c / 2 * ((4 - 2*k)*i[1] - (4*y + 2*k)*i[2] + 4*y*i[0]);
  e2        = dr * dr * echi2;
  return r;
}

//____________________________________________________________________
Double_t 
AliFMDFlowResolutionTDR::Correction(UShort_t k) const 
{ 
  Double_t e;
  return Correction(k, e); 
}
//____________________________________________________________________
Double_t 
AliFMDFlowResolutionTDR::Chi2Over2(Double_t r, Double_t& e2) const 
{
  // From nucl-ex/9711003v2 
  //
  // N equal to the total number of events, and k equal to the
  // number of events where |Psi_A - Psi_B| > pi/2, we get 
  //
  //   k   exp(-chi^2/2)                               k
  //   - = -------------     =>    chi^2 / 2 = - log 2 -
  //   N          2 		                       N
  //
  //                                               N
  //                         =>    chi^2 / 2 = log --   (*)
  //                                               k
  //            
  //                                                   2 k
  //                         =>    chi^2     = - 2 log ---
  //                                                    N
  //                                                            
  //                         =>    chi       = -/+ sgrt (- 2 log (2k/N))
  //
  // (*) this is what is returned.  
  //
  //  The error on chi is given by 
  //
  //              dchi                           1
  //    d^2chi = (------)^2 d^2(k/N) = - -------------------- d^2(k/N)
  //              d(k/N)                 4 (k/N)^2 log(2 k/N)
  // 
  // where 
  //                      k         k^2
  //               (1 - 2 -) dk^2 + --- dN^2
  //                      N	    N^2
  //   d^2(k/N) = --------------------------
  //                           N^2 
  // 
  // 
  // with dk = sqrt(k) and dN = sqrt(N), this reduces to 
  //
  //                     k      k^2        N k     k^2   k^2
  //              (1 - 2 -) k + ----- N    --- - 2 --- + ---
  //                     N      N^2         N       N     N
  //   d^2(k/N) = --------------------- = --------------------------
  //                           N^2                 N^2 
  //
  //              k N - k^2   k (N - k)   k (1 - k / N)
  //            = --------- = - ------- = - -----------
  //                 N^3      N   N^2     N      N
  //
  // Identifying r = k/N, we get 
  //
  //            
  //   d(k/N)   = sqrt(r (1 - r) / N) 
  //
  // which is the well-known result for Binomial errors.
  // Alternatively, one could compute these error using an Baysian
  // confidence level of say 68.3% (see for example 
  // http://home.fnal.gov/~paterno/images/effic.pdf). 
  // 
  // Our final expression for the error on chi then becomes 
  //
  //                      1	     
  //   d^2chi = - -------------------- r (1 - r) / N
  //		  4 r^2 log(2 r)
  //
  //                 1 - r            r - 1
  //          = - -------------- = ------------
  //              4 r N log(2 r)   4 k log(2 r)  
  if (r == 0) { 
    std::cerr << "TDR: Large/All = " << r << " <= 0!" << std::endl;
    return 0;
  }
  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;
  }
  Double_t chi2over2  = log(ratio);
  if (chi2over2 < 0) { 
    std::cerr << "TDR: log(" << ratio << ") = " << chi2over2 
	      << " < 0" << std::endl; 
    return -1;
  }
  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;
  }
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
39eefe19	1	/** @file
	2	@brief Implementation of an Resolution class */
	3	#include "flow/AliFMDFlowResolution.h"
	4	#include "flow/AliFMDFlowUtil.h"
	5	#include "flow/AliFMDFlowBessel.h"
824e71c1	6	#include <TGraphErrors.h>
	7	#include <TH1.h>
	8	#include <TH2.h>
39eefe19	9	#include <iostream>
	10	//#include <cmath>
	11
	12	//====================================================================
	13	void
	14	AliFMDFlowResolution::Add(Double_t psiA, Double_t psiB)
	15	{
	16	Double_t diff = NormalizeAngle(fOrder * (psiA - psiB));
	17	Double_t contrib = cos(diff);
	18	AliFMDFlowStat::Add(contrib);
	19	}
	20
	21	//____________________________________________________________________
	22	Double_t
	23	AliFMDFlowResolution::Correction(UShort_t k) const
	24	{
	25	Double_t e;
	26	return Correction(k, e);
	27	}
	28
	29	//____________________________________________________________________
	30	Double_t
	31	AliFMDFlowResolution::Correction(UShort_t, Double_t& e2) const
	32	{
	33	e2 = fSqVar / fN;
	34	return sqrt(2) * sqrt(fabs(fAverage));
	35	}
	36
824e71c1	37	//____________________________________________________________________
	38	void
	39	AliFMDFlowResolution::Draw(Option_t* option)
	40	{
	41	TGraph* g = new TGraph(100);
	42	for (UShort_t i = 0; i < g->GetN(); i++) {
	43	Double_t x = -1. + 2. / 100 * i;
	44	Double_t y = sqrt(2) * sqrt(fabs(x));
	45	g->SetPoint(i, x, y);
	46	}
	47	g->SetName("naive_res");
	48	g->SetTitle("Naive Resolution Function");
	49	g->GetHistogram()->SetXTitle("<cos(n(\|#Psi_{A}-#Psi_{B}\|))>");
	50	g->GetHistogram()->SetYTitle("<cos(n(#Psi-#Psi_{R}))>");
	51	g->GetHistogram()->SetDirectory(0);
	52	g->Draw(Form("lh %s", option));
	53	}
	54
	55
39eefe19	56	//====================================================================
	57	Double_t
	58	AliFMDFlowResolutionStar::Correction(UShort_t k, Double_t& e2) const
	59	{
	60	if (k > 4) return 0;
	61	Double_t delta = 0;
	62	Double_t chi = Chi(fAverage, k, delta);
	63	Double_t dr = 0;
	64	Double_t res = Res(sqrt(2) * chi, k, dr);
824e71c1	65	e2 = pow(dr * delta,2);
39eefe19	66	return res;
	67	}
	68	//____________________________________________________________________
	69	Double_t
	70	AliFMDFlowResolutionStar::Correction(UShort_t k) const
	71	{
	72	Double_t e;
	73	return Correction(k, e);
	74	}
	75	//____________________________________________________________________
	76	Double_t
	77	AliFMDFlowResolutionStar::Chi(Double_t res, UShort_t k,
	78	Double_t& delta) const
	79	{
824e71c1	80	delta = 1;
39eefe19	81	Double_t chi = 2;
	82	Double_t dr = 0;
	83	for (UInt_t i = 0; i < 15; i++) {
	84	if (Res(chi, k, dr) < res) chi += delta;
	85	else chi -= delta;
	86	delta /= 2;
	87	}
	88	return chi;
	89	}
	90	//____________________________________________________________________
	91	Double_t
824e71c1	92	AliFMDFlowResolutionStar::Res(Double_t chi, UShort_t k, Double_t& dr) const
39eefe19	93	{
	94	// The resolution function is
	95	//
	96	// sqrt(pi) x exp(-x/4) (f1(x^2/4) + f2(x^2/4))
	97	// r(x) = --------------------------------------------
	98	// 2 sqrt(2)
	99	//
	100	//
	101	// = c x (f1(y) - f2(y))
	102	//
	103	// where f1 is the modified Bessel function first kind I_{(k-1)/2},
	104	// and f2 is the modified Bessel function of the first kind
	105	// I_{(k+1)/2}, and
	106	//
	107	// sqrt(pi) exp(-x^2/4)
	108	// c = --------------------, y = x^2/4
	109	// 2 sqrt(2)
	110	//
	111	// The derivative of the resolution function is
	112	//
	113	// c
	114	// r'(y) = - (4 sqrt(y) (f1'(y) - f2'(y)) - (4 y - 2)(f1(y) - f2(y)))
	115	// 2
	116	//
	117	// c r(y)
	118	// = - (4 sqrt(y) (f1'(y) - f2'(y))) + --------- - sqrt(y) r(y)
	119	// 2 2 sqrt(y)
	120	//
	121	// Since dI_n(x)/dx = I_(n-1)(x) - n / x I_n(x), and substituting
	122	// f3(y) = I_((k-3)/2)(y)
	123	//
	124	// c
	125	// r'(y) = - ((4 - 2 k) f1(y) - (4 y + 2 k) f2(y) + 4 y f3(y))
	126	// 2
	127	//
	128	Double_t y = chi * chi / 4;
	129	Double_t c = sqrt(M_PI) * exp(-y) / (2 * sqrt(2));
	130
	131	// i[0] = I[(k-1)/2-1], i[1] = I[(k-1)/2], i[2] = I[(k-1)/2+1]
	132	Double_t i[3], di[3];
	133	AliFMDFlowBessel::Inu(Double_t(k-3)/2, Double_t(k+1)/2, y, i, di);
	134
	135	Double_t r = c * chi * (i[2] + i[1]);
	136	dr = c / 2 * ((4 - 2k)i[1] - (4y + 2k)i[2] + 4y*i[0]);
	137
	138	return r;
	139	}
	140
824e71c1	141	//____________________________________________________________________
	142	void
	143	AliFMDFlowResolutionStar::Draw(Option_t* option)
	144	{
	145	TString opt(option);
	146	opt.ToLower();
	147	Bool_t chi = opt.Contains("chi");
	148
	149	TH2* h = new TH2D(Form("star_%s_frame", (chi ? "chi" : "res")),
	150	Form("STAR %s Function for k=(1,2,4)",
	151	(chi ? "Chi" : "Resolution")),
	152	100, 0, 1, 100, 0, (chi ? 3 : 1.5));
	153	h->SetXTitle("<cos(n(\|#Psi_{A}-#Psi_{B}\|))>");
	154	h->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>"));
	155	h->SetStats(0);
	156	h->SetDirectory(0);
	157	h->Draw();
	158
	159	for (UShort_t k = 1; k <= 4; k++) {
	160	if (k == 3) continue;
	161
	162	TGraphErrors* g = new TGraphErrors(100);
	163	for (UShort_t i = 0; i < g->GetN(); i++) {
	164	Double_t e2 = 0;
	165	Double_t x = 1. / 100 * i;
	166	Double_t y = 0;
	167	Double_t c = Chi(x, k, e2);
	168	if (chi) y = c;
	169	else {
	170	Double_t dr = 0;
	171	y = Res(sqrt(2) * c, k, dr);
	172	e2 = pow(dr * e2,2);
	173	}
	174	g->SetPoint(i, x, y);
	175	g->SetPointError(i, 0, sqrt(e2));
	176	}
	177	g->SetLineColor(k);
	178	g->SetName(Form("star_%s_k%d", (chi ? "chi" : "res"), k));
	179	g->SetTitle(Form("STAR %s Function for k=%d",
	180	(chi ? "Chi" : "Resolution"), k));
	181	g->GetHistogram()->SetXTitle("<cos(n(\|#Psi_{A}-#Psi_{B}\|))>");
	182	g->GetHistogram()->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>"));
	183	g->GetHistogram()->SetDirectory(0);
	184	g->Draw(Form("l %s same", option));
	185	}
	186	}
	187
39eefe19	188	//====================================================================
39eefe19	189	void
824e71c1	190	AliFMDFlowResolutionTDR::Clear(Option_t*)
39eefe19	191	{
	192	fN = 0;
	193	fLarge = 0;
	194	}
	195	//____________________________________________________________________
	196	void
	197	AliFMDFlowResolutionTDR::Add(Double_t psi_a, Double_t psi_b)
	198	{
	199	Double_t a = fabs(psi_a - psi_b);
	200	if (a >= .5 * M_PI) fLarge++;
	201	fN++;
	202	}
	203	//____________________________________________________________________
	204	Double_t
	205	AliFMDFlowResolutionTDR::Correction(UShort_t k, Double_t& e2) const
	206	{
	207	// From nucl-ex/9711003v2
	208	//
	209	//
	210	// <cos n Delta phi> =
	211	//
	212	// sqrt(pi)
	213	// -------- chi exp(-z) (I_((n-1)/2)(z) + I_((n+1)/2)(z))
	214	// 2
	215	//
	216	// where z = chi^2 / 2
	217	//
824e71c1	218	if (fLarge == 0) {
	219	std::cerr << "TDR: K = 0" << std::endl;
	220	return -1;
	221	}
	222	if (fN == 0) {
	223	std::cerr << "TDR: N = 0" << std::endl;
	224	return -1;
	225	}
	226	Double_t r = Double_t(fLarge) / fN;
39eefe19	227	Double_t echi2 = 0;
824e71c1	228	Double_t y = Chi2Over2(r, echi2);
	229	return Res(k, y, echi2, e2);
	230	}
	231
	232
	233
	234	//____________________________________________________________________
	235	Double_t
	236	AliFMDFlowResolutionTDR::Res(UShort_t k, Double_t y, Double_t echi2,
	237	Double_t& e2) const
	238	{
	239	// y = chi^2 / 2
39eefe19	240	Double_t chi = sqrt(2 * y);
39eefe19	241	Double_t c = sqrt(M_PI) * exp(-y) / 2;
824e71c1	242
39eefe19	243	// i[0] = I[(k-1)/2-1], i[1] = I[(k-1)/2], i[2] = I[(k-1)/2+1]
	244	Double_t i[3], di[3];
	245	AliFMDFlowBessel::Inu(Double_t(k-3)/2, Double_t(k+1)/2, y, i, di);
	246
	247	Double_t r = c * chi * (i[2] + i[1]);
	248	Double_t dr = c / 2 * ((4 - 2k)i[1] - (4y + 2k)i[2] + 4y*i[0]);
	249	e2 = dr * dr * echi2;
	250	return r;
	251	}
824e71c1	252
39eefe19	253	//____________________________________________________________________
	254	Double_t
	255	AliFMDFlowResolutionTDR::Correction(UShort_t k) const
	256	{
	257	Double_t e;
	258	return Correction(k, e);
	259	}
	260	//____________________________________________________________________
	261	Double_t
824e71c1	262	AliFMDFlowResolutionTDR::Chi2Over2(Double_t r, Double_t& e2) const
39eefe19	263	{
	264	// From nucl-ex/9711003v2
	265	//
	266	// N equal to the total number of events, and k equal to the
	267	// number of events where \|Psi_A - Psi_B\| > pi/2, we get
	268	//
	269	// k exp(-chi^2/2) k
	270	// - = ------------- => chi^2 / 2 = - log 2 -
	271	// N 2 N
	272	//
	273	// N
	274	// => chi^2 / 2 = log -- (*)
	275	// k
	276	//
	277	// 2 k
	278	// => chi^2 = - 2 log ---
	279	// N
	280	//
	281	// => chi = -/+ sgrt (- 2 log (2k/N))
	282	//
	283	// (*) this is what is returned.
	284	//
	285	// The error on chi is given by
	286	//
	287	// dchi 1
	288	// d^2chi = (------)^2 d^2(k/N) = - -------------------- d^2(k/N)
	289	// d(k/N) 4 (k/N)^2 log(2 k/N)
	290	//
	291	// where
	292	// k k^2
	293	// (1 - 2 -) dk^2 + --- dN^2
	294	// N N^2
	295	// d^2(k/N) = --------------------------
	296	// N^2
	297	//
	298	//
	299	// with dk = sqrt(k) and dN = sqrt(N), this reduces to
	300	//
	301	// k k^2 N k k^2 k^2
	302	// (1 - 2 -) k + ----- N --- - 2 --- + ---
	303	// N N^2 N N N
	304	// d^2(k/N) = --------------------- = --------------------------
	305	// N^2 N^2
	306	//
	307	// k N - k^2 k (N - k) k (1 - k / N)
	308	// = --------- = - ------- = - -----------
	309	// N^3 N N^2 N N
	310	//
	311	// Identifying r = k/N, we get
	312	//
	313	//
	314	// d(k/N) = sqrt(r (1 - r) / N)
	315	//
	316	// which is the well-known result for Binomial errors.
	317	// Alternatively, one could compute these error using an Baysian
	318	// confidence level of say 68.3% (see for example
	319	// http://home.fnal.gov/~paterno/images/effic.pdf).
	320	//
	321	// Our final expression for the error on chi then becomes
	322	//
	323	// 1
	324	// d^2chi = - -------------------- r (1 - r) / N
	325	// 4 r^2 log(2 r)
	326	//
327	// 1 - r r - 1
328	// = - -------------- = ------------
329	// 4 r N log(2 r) 4 k log(2 r)
824e71c1	330	if (r == 0) {
	331	std::cerr << "TDR: Large/All = " << r << " <= 0!" << std::endl;
	332	return 0;
39eefe19	333	}
824e71c1	334	Double_t ratio = 1. / (2r); // Double_t(fN) / (2 fLarge);
39eefe19	335	if (ratio <= 0) {
	336	std::cerr << "TDR: Large/All = " << ratio << " <= 0!" << std::endl;
	337	return -1;
	338	}
	339	Double_t chi2over2 = log(ratio);
	340	if (chi2over2 < 0) {
	341	std::cerr << "TDR: log(" << ratio << ") = " << chi2over2
	342	<< " < 0" << std::endl;
	343	return -1;
	344	}
824e71c1	345	if (fLarge != 0) e2 = (r - 1) / (4 * fLarge * log(2 * r));
824e71c1	346	else e2 = (r - 1) / (4 * r * log(2 * r));
39eefe19	347	return chi2over2;
	348	}
	349
824e71c1	350	//____________________________________________________________________
	351	void
	352	AliFMDFlowResolutionTDR::Draw(Option_t* option)
	353	{
	354	TString opt(option);
	355	opt.ToLower();
	356	Bool_t chi = opt.Contains("chi");
	357
	358	TH2* h = new TH2D(Form("tdr_%s_frame", (chi ? "chi" : "res")),
	359	Form("TDR %s Function for k=(1,2,4)",
	360	(chi ? "Chi" : "Resolution")),
	361	100, 0, 1, 100, 0, (chi ? 3 : 1.5));
	362	h->SetXTitle("K/N");
	363	h->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>"));
	364	h->SetStats(0);
	365	h->Draw();
	366	h->SetDirectory(0);
	367
	368	for (UShort_t k = 1; k <= 4; k++) {
	369	if (k == 3) continue;
	370
	371	TGraphErrors* g = new TGraphErrors;
	372	Int_t i = 0;
	373	for (Double_t x = 0.02; x < 0.5; x += 0.01) {
	374	Double_t e2 = 0;
	375	Double_t y = 0;
	376	Double_t c = Chi2Over2(x, e2);
	377	if (chi) y = sqrt(2 * c);
	378	else {
	379	Double_t dr = 0;
	380	y = Res(k, c, e2, dr);
	381	e2 = dr;
	382	}
	383	g->SetPoint(i, x, y);
	384	g->SetPointError(i, 0, sqrt(fabs(e2)));
	385	i++;
	386	}
	387	g->SetLineColor(k);
	388	g->SetName(Form("tdr_%s_k%d", (chi ? "chi" : "res"), k));
	389	g->SetTitle(Form("TDR %s Function for k=%d",
	390	(chi ? "Chi" : "Resolution"), k));
	391	g->GetHistogram()->SetXTitle("K/N");
	392	g->GetHistogram()->SetYTitle((chi ? "#chi" : "<cos(n(#Psi-#Psi_{R}))>"));
	393	g->GetHistogram()->SetDirectory(0);
	394	g->Draw(Form("l %s %s", (k == 0 ? "h" : "same"), option));
	395	if (chi) break;
	396	}
	397	}
39eefe19	398
	399	//____________________________________________________________________
	400	//
	401	// EOF
	402	//