Added drawing capabilities
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 19 Sep 2007 13:12:50 +0000 (13:12 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 19 Sep 2007 13:12:50 +0000 (13:12 +0000)
FMD/flow/AliFMDFlowAxis.cxx
FMD/flow/AliFMDFlowAxis.h
FMD/flow/AliFMDFlowBin.cxx
FMD/flow/AliFMDFlowBin.h
FMD/flow/AliFMDFlowBinned1D.cxx
FMD/flow/AliFMDFlowBinned1D.h
FMD/flow/AliFMDFlowResolution.cxx
FMD/flow/AliFMDFlowResolution.h
FMD/flow/AliFMDFlowStat.h
FMD/flow/AliFMDFlowTrue.h

index 51e1a4f..dc35464 100644 (file)
@@ -2,6 +2,8 @@
     @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) 
@@ -91,6 +93,19 @@ AliFMDFlowAxis::BinUpper(UShort_t i) const
 }
 
 //____________________________________________________________________
+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
 //
index 81a5f6f..045bc02 100644 (file)
@@ -57,6 +57,8 @@ public:
   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;
index eeb7842..d7edbd9 100644 (file)
@@ -4,6 +4,7 @@
 #include <cmath>
 #include <iostream>
 #include <iomanip>
+#include <TBrowser.h>
 
 //====================================================================
 void 
@@ -114,6 +115,19 @@ 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
 {
index 67205ea..f40b747 100644 (file)
@@ -8,6 +8,9 @@
 #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
@@ -100,8 +103,11 @@ public:
       @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 */
index b28fc1f..9e5bc6d 100644 (file)
@@ -8,6 +8,8 @@
 #include <iomanip>
 #include <TBrowser.h>
 #include <TString.h>
+#include <TH1.h>
+#include <TH2.h>
 
 //====================================================================
 AliFMDFlowBinned1D::AliFMDFlowBinned1D(UShort_t order, 
@@ -141,6 +143,71 @@ AliFMDFlowBinned1D::Browse(TBrowser* b)
 
 //____________________________________________________________________
 void 
+AliFMDFlowBinned1D::Draw(Option_t* option)
+{
+  TString opt(option);
+  opt.ToLower();
+  const char* names[] = { "Bare", "Naive", "STAR", "TDR" };
+  const AliFMDFlowBin::CorType types[] = { AliFMDFlowBin::none, 
+                                          AliFMDFlowBin::naive, 
+                                          AliFMDFlowBin::star, 
+                                          AliFMDFlowBin::tdr };
+  Bool_t meths[] = { opt.Contains("b"), 
+                    opt.Contains("n"),
+                    opt.Contains("s"), 
+                    opt.Contains("t") };
+  Bool_t res     = opt.Contains("r");
+  UShort_t nm = 0;
+  Short_t  sm = -1;
+  for (UShort_t i = 0; i < 4; i++) { if (meths[i]) { nm++; sm = i; } }
+  TH1* h = 0;
+  if (nm > 1) { 
+    h = new TH2D((res ? "res" : "flow"), (res ? "Resolution" : "Flow"), 
+                fXAxis.N(), fXAxis.Bins(), nm, 0, nm);
+    h->SetXTitle("x");
+    h->SetYTitle("method");
+    h->GetYaxis()->SetNdivisions(nm+1, kFALSE);
+    h->SetZTitle((res ? "<cos(n(#Psi_{m}-#Psi_{R}))>" : "v"));
+    UInt_t j = 0;
+    for (UShort_t i = 0; i < 4; i++) 
+      if (meths[i]) h->GetYaxis()->SetBinLabel(++j, names[i]);
+  }
+  else {
+    h = new TH1D(Form("%s_%s", (res ? "res" : "flow"), names[sm]),
+                Form("%s_%s", (res ? "Resolution" : "Flow"), names[sm]),
+                fXAxis.N(), fXAxis.Bins());
+    h->SetXTitle("x");
+    h->SetYTitle((res ? "<cos(n(#Psi_{m}-#Psi_{R}))>" : "v"));
+  }
+
+  for (UShort_t i = 0; i < fXAxis.N(); i++) { 
+    Double_t       x   = fXAxis.BinCenter(i);
+    AliFMDFlowBin* bin = GetBin(x);
+    Double_t       v, e2;
+    if (nm == 1) { 
+      if (res) v = bin->Correction(e2, types[sm]);
+      else     v = bin->Value(e2, types[sm]);
+      h->SetBinContent(i+1, v);
+      h->SetBinError(i+1, sqrt(e2));
+      continue;
+    }
+    UInt_t j = 0;
+    for (UShort_t k = 0; k < 4; k++)  { 
+      if (!meths[k]) continue;
+      if (res) v = bin->Correction(e2, types[k]);
+      else     v = bin->Value(e2, types[k]);
+      h->SetBinContent(i+1, j+1, v);
+      h->SetBinError(i+1, j+1, sqrt(e2));
+      j++;
+    }
+  }
+  h->Draw(option);
+}
+
+  
+  
+//____________________________________________________________________
+void 
 AliFMDFlowBinned1D::Print(Option_t* option) const
 {
   TString opt(option);
index a145f61..e65d812 100644 (file)
@@ -75,7 +75,16 @@ public:
       @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 */ 
index 57fdc0e..6a327bd 100644 (file)
@@ -3,6 +3,9 @@
 #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>
 
@@ -31,6 +34,25 @@ AliFMDFlowResolution::Correction(UShort_t, Double_t& e2) const
   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 
@@ -40,7 +62,7 @@ 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;
 }
 //____________________________________________________________________
@@ -55,7 +77,7 @@ Double_t
 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++) { 
@@ -67,8 +89,7 @@ AliFMDFlowResolutionStar::Chi(Double_t res, UShort_t k,
 }
 //____________________________________________________________________
 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 
   // 
@@ -117,9 +138,56 @@ AliFMDFlowResolutionStar::Res(Double_t chi, UShort_t k,
   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;
@@ -147,11 +215,31 @@ AliFMDFlowResolutionTDR::Correction(UShort_t k, Double_t& e2) const
   // 
   // 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);
@@ -161,6 +249,7 @@ AliFMDFlowResolutionTDR::Correction(UShort_t k, Double_t& e2) const
   e2        = dr * dr * echi2;
   return r;
 }
+
 //____________________________________________________________________
 Double_t 
 AliFMDFlowResolutionTDR::Correction(UShort_t k) const 
@@ -170,7 +259,7 @@ 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 
   //
@@ -238,12 +327,11 @@ AliFMDFlowResolutionTDR::Chi2Over2(Double_t& e2) const
   //                 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;
@@ -254,11 +342,59 @@ AliFMDFlowResolutionTDR::Chi2Over2(Double_t& e2) const
              << " < 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;
+  }
+}
 
 //____________________________________________________________________
 //
index e3954ff..021231e 100644 (file)
@@ -69,6 +69,9 @@ public:
   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;
@@ -146,6 +149,9 @@ public:
       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$
@@ -208,7 +214,7 @@ public:
     : 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 
@@ -223,8 +229,22 @@ public:
   /** 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 */
index 223de72..5b1cb0f 100644 (file)
@@ -30,7 +30,7 @@
     \end{array}\right.
     @f] 
 */
-class AliFMDFlowStat
+class AliFMDFlowStat : public TObject
 {
 public:
   /** Constructor */
@@ -38,7 +38,7 @@ public:
   /** 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 */
@@ -60,6 +60,14 @@ protected:
 
 //__________________________________________________________________
 inline void 
+AliFMDFlowStat::Clear(Option_t*) 
+{ 
+  fAverage = fSqVar = 0; 
+  fN = 0; 
+}
+
+//__________________________________________________________________
+inline void 
 AliFMDFlowStat::Add(Double_t x) 
 { 
   if (TMath::IsNaN(x) || !TMath::Finite(x)) return;
index eef57b7..1da887e 100644 (file)
@@ -49,7 +49,7 @@ struct AliFMDFlowTrueBin : public AliFMDFlowBin
       @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 */