Fix to the raw reader that properly decodes the parameter settings in the
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 28 Oct 2009 12:56:13 +0000 (12:56 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 28 Oct 2009 12:56:13 +0000 (12:56 +0000)
RCU trailer.

Raw writer prepared for getting the number of bits written.

Added comment in display.

Updated flow examples.

FMD/AliFMDDisplay.cxx
FMD/AliFMDRawReader.cxx
FMD/AliFMDRawWriter.cxx
FMD/AliFMDRawWriter.h
FMD/flow/TestFlow.C
FMD/flow/TestFlowSimple.C [new file with mode: 0644]

index 2122e77..b587467 100644 (file)
@@ -901,6 +901,7 @@ AliFMDDisplay::ProcessESD(UShort_t det, Char_t rng, UShort_t sec, UShort_t str,
   return kTRUE;
 }
 
+//____________________________________________________________________
 Double_t
 AliFMDDisplay::GetADCThreshold(UShort_t d, Char_t r, 
                               UShort_t s, UShort_t t) const
index 869f9a2..378b171 100644 (file)
@@ -132,22 +132,81 @@ AliFMDRawReader::NewDDL(AliAltroRawStreamV3& input, UShort_t& det)
   UInt_t ddl = input.GetDDLNumber();
   AliFMDDebug(2, ("DDL number %d", ddl));
 
-  // Get zero suppression flag
-  fZeroSuppress[ddl] = input.GetZeroSupp();
+  /* Note, previously, the ALTROCFG1 register was interpreted as 
+   * 
+   * Bits    Value    Description
+   *   0- 3      0/1   1st Baseline filter, mode 
+   *   4- 5   Over-1   2nd baseline filter, # of pre-samples
+   *   6- 9   factor   2nd baseline filter, # of post-samples 
+   *  10-          0   2nd baseline filter, enable
+   *  11-12       00   Zero suppression, glitch filter mode
+   *  13-15      001   Zero suppression, # of post samples
+   *  16-17       01   Zero suppression, # of pre  samples
+   *  18         0/1   Zero suppression, enable
+   *
+   * The interpretation used in AliAltroRawStreamerV3 - which
+   * corresponds directly to ALTRO DPCFG register - is 
+   *
+   * Bits    Value  Description
+   *   0- 3    0/1   1st Baseline filter, mode 
+   *   4         0   Polarity (if '1', then "1's inverse")
+   *   5- 6     01   Zero suppression, # of pre samples
+   *   7-10   0001   Zero suppression, # of post samples
+   *  11         0   2nd baseline filter, enable
+   *  12-13     00   Zero suppression, glitch filter mode
+   *  14-16 factor   2nd baseline filter, # of post-samples
+   *  17-18     01   2nd baseline filter, # of pre-samples 
+   *  19       0/1   Zero suppression, enable
+   *
+   *  Writing 'x' for variable values, that means we have the
+   *  following patterns for the 2 cases 
+   *
+   *    bit #  20   16   12    8    4    0
+   *     old    |0x01|0010|00xx|xxxx|xxxx|
+   *     new    |x01x|xx00|0000|1010|xxxx|
+   *
+   *  That means that we can check if bits 10-13 are '1000' or
+   *  '0000', which will tell us if the value was written with the
+   *  new or the old interpretation.    That is, we can check that 
+   *
+   *    if (((altrocfg1 >> 10) & 0x8) == 0x8) { 
+   *      // old interpretation 
+   *    }
+   *    else { 
+   *      // New interpretation 
+   *    }
+   *
+   * That means, that we should never 
+   *
+   *  - change the # of zero suppression post samples 
+   *  - Turn on 2nd baseline correction 
+   *  - Change the zero-suppression glitch filter mode
+   *
+   * This change as introduced in version 1.2 of Rcu++
+   */
+  UInt_t cfg1 = input.GetAltroCFG1();
+  if (((cfg1 >> 10) & 0x8) == 0x8) {
+    UInt_t cfg2 = input.GetAltroCFG2();
+    fZeroSuppress[ddl] = (cfg1 >> 0) & 0x1;
+    fNoiseFactor[ddl]  = (cfg1 >> 6) & 0xF;
+    fSampleRate[ddl]   = (cfg2 >> 0) & 0xF;
+  }
+  else {
+    fZeroSuppress[ddl] = input.GetZeroSupp();
+    // WARNING: We store the noise factor in the 2nd baseline
+    // filters excluded post samples, since we'll never use that
+    // mode. 
+    fNoiseFactor[ddl]  = input.GetNPostsamples();
+    // WARNING: We store the sample rate in the number of pre-trigger
+    // samples, since we'll never use that mode.
+    fSampleRate[ddl]     = input.GetNPretriggerSamples();
+  }
   AliFMDDebug(3, ("RCU @ DDL %d zero suppression: %s", 
                   ddl, (fZeroSuppress[ddl] ? "yes" : "no")));
-
-  // WARNING: We store the noise factor in the 2nd baseline
-  // filters excluded post samples, since we'll never use that
-  // mode. 
-  fNoiseFactor[ddl]  = input.GetNPostsamples();
-  AliFMDDebug(3, ("RCU @ DDL %d noise factor: %d", ddl,fNoiseFactor[ddl]));
-    
-  // WARNING: We store the sample rate in the number of pre-trigger
-  // samples, since we'll never use that mode.
-  fSampleRate[ddl]     = input.GetNPretriggerSamples();
+  AliFMDDebug(3, ("RCU @ DDL %d noise factor: %d", ddl,fNoiseFactor[ddl]));    
   AliFMDDebug(3, ("RCU @ DDL %d sample rate: %d", ddl,fSampleRate[ddl]));
 
+
   // Get Errors seen 
   Int_t nChAddrMismatch = input.GetNChAddrMismatch();
   Int_t nChLenMismatch  = input.GetNChLengthMismatch();
index acb8687..4585019 100644 (file)
@@ -140,14 +140,15 @@ AliFMDRawWriter::Exec(Option_t*)
     return;
   }
   
-  TClonesArray* digits = new TClonesArray("AliFMDDigit", 1000);
   fFMD->SetTreeAddress();
-  TBranch* digitBranch = digitTree->GetBranch(fFMD->GetName());
-  if (!digitBranch) {
-    AliError(Form("no branch for %s", fFMD->GetName()));
-    return;
-  }
-  digitBranch->SetAddress(&digits);
+  TClonesArray* digits = fFMD->Digits(); 
+  // new TClonesArray("AliFMDDigit", 1000);
+  // TBranch* digitBranch = digitTree->GetBranch(fFMD->GetName());
+  // if (!digitBranch) {
+  //   AliError(Form("no branch for %s", fFMD->GetName()));
+  //   return;
+  // }
+  // digitBranch->SetAddress(&digits);
   
   Int_t nEvents = Int_t(digitTree->GetEntries());
   AliFMDDebug(5, ("Got a total of %5d events from tree", nEvents));
@@ -164,12 +165,12 @@ AliFMDRawWriter::Exec(Option_t*)
 
 #if 1
 //____________________________________________________________________
-void
+Long_t
 AliFMDRawWriter::WriteDigits(TClonesArray* digits)
 {
   // WRite an array of digits to disk file 
   Int_t nDigits = digits->GetEntries();
-  if (nDigits < 1) return;
+  if (nDigits < 1) return 0;
   AliFMDDebug(5, ("Got a total of %5d digits from tree", nDigits));
 
   AliFMDParameters* pars = AliFMDParameters::Instance();
@@ -192,9 +193,10 @@ AliFMDRawWriter::WriteDigits(TClonesArray* digits)
   // The Altro buffer 
   AliAltroBufferV3* altro = 0;
   
-  Int_t totalWords = 0;
-  Int_t nCounts    = 0;
-  
+  Int_t  totalWords = 0;
+  Int_t  nCounts    = 0;
+  Long_t nBits      = 0;
+
   // Loop over the digits in the event.  Note, that we assume the
   // the digits are in order in the branch.   If they were not, we'd
   // have to cache all channels before we could write the data to
@@ -240,7 +242,8 @@ AliFMDRawWriter::WriteDigits(TClonesArray* digits)
                         time, prevaddr, nWords));
        totalWords += nWords;
        ZeroSuppress(data.fArray, nWords, peds.fArray, noise.fArray, threshold);
-       if (altro) altro->WriteChannel(prevaddr,nWords,data.fArray,threshold);
+       if (altro) 
+         /*nBits+=*/altro->WriteChannel(prevaddr,nWords,data.fArray,threshold);
        data.Reset(-1);
        peds.Reset(0);
        noise.Reset(0);
@@ -255,8 +258,8 @@ AliFMDRawWriter::WriteDigits(TClonesArray* digits)
          // When the first argument is false, we write the real
          // header. 
          AliFMDDebug(15, ("Closing output"));
-         altro->Flush();
-         altro->WriteDataHeader(kFALSE, kFALSE);
+         /* nBits += */ altro->Flush();
+         /* nBits += */ altro->WriteDataHeader(kFALSE, kFALSE);
          delete altro;
          altro = 0;
        }
@@ -300,13 +303,15 @@ AliFMDRawWriter::WriteDigits(TClonesArray* digits)
   // already 
   if (altro) {
     ZeroSuppress(data.fArray, nWords, peds.fArray, noise.fArray, threshold);
-    if (nWords > 0) altro->WriteChannel(prevaddr,nWords,data.fArray,threshold);
-    altro->Flush();
-    altro->WriteDataHeader(kFALSE, kFALSE);
+    if (nWords > 0) 
+      /* nBits += */ altro->WriteChannel(prevaddr,nWords,data.fArray,threshold);
+    /* nBits += */ altro->Flush();
+    /* nBits += */ altro->WriteDataHeader(kFALSE, kFALSE);
     delete altro;
   }
-  AliFMDDebug(5, ("Wrote a total of %d words for %d counts", 
-                 nWords, nCounts));
+  AliFMDDebug(5, ("Wrote a total of %d words in %d bytes for %d counts", 
+                 nWords, nBits / 8, nCounts));
+  return nBits;
 }
 //____________________________________________________________________
 void
index d1c309f..e2c2947 100644 (file)
@@ -56,7 +56,7 @@ public:
   /** Write an array of AliFMDDigit objects as raw ALTRO data. 
       @param digits Array of AliFMDDigit objects to convert to raw
       ALTRO data. */
-  virtual void WriteDigits(TClonesArray* digits);
+  virtual Long_t WriteDigits(TClonesArray* digits);
   /** Do zero-suppression of channel data. 
       @param data      Contain @a nWords of valid data.  On input, it 
                        contains the full channel data.  On output it
index 9f52947..38729a1 100644 (file)
     or 
     @verbatim 
     $ g++ `root-config --cflags --libs` -I$ALICE_ROOT/include \
-       -I$ALICE_ROOT/FMD $ALICE_ROOT/FMD/flow/TestFlow.C -o testflow
+       -I$ALICE_ROOT/FMD -L$ALICE_ROOT/lib/tgt_${ALICE_TARGET} \
+       -lFMDflow $ALICE_ROOT/FMD/flow/TestFlow.C -o testflow
     $ ./testflow 
     @endverbatim 
 */
-#include <FMD/flow/AliFMDFlowBinned1D.h>
-#include <FMD/flow/AliFMDFlowTrue.h>
-#include <FMD/flow/AliFMDFlowUtil.h>
+#include <flow/AliFMDFlowBinned1D.h>
+#include <flow/AliFMDFlowTrue.h>
+#include <flow/AliFMDFlowUtil.h>
 #include <TRandom.h>
 #include <TCanvas.h>
 #include <TVirtualPad.h>
+#include <TLegend.h>
 #include <TArrayD.h>
 #include <TBrowser.h>
 #include <iostream>
@@ -94,100 +96,241 @@ struct Generator
   UInt_t   fN;
 };
 
-/** Run test program */
-void
-TestFlow(UInt_t n_events=100, Int_t seg=-1, UInt_t n_max=20000,
-        Float_t v2=0.05)
+struct Tester 
 {
-  Generator           generator(-1, 0.00, v2, n_max, n_max);
-  AliFMDFlowAxis      a(10, -5, 5);
-  AliFMDFlowBinned1D  flow("test","Test",2, 1, a);
-  AliFMDFlowTrue1D    real("true","true", 2, a);
-  TH2D                hist("hist","hist",a.N(),a.Bins(),
-                          (seg>0?seg:90),0, 2* TMath::Pi());
-  TH1D                rela("rela","real",
-                          (seg>0?seg:90),0, 2* TMath::Pi());
-  Double_t            dphi = (seg <= 0 ? 0 : 2 * TMath::Pi() / seg);
-  TArrayD             phis(n_max);
-  TArrayD             tphis(n_max);
-  TArrayD             xs(n_max);
-  flow.SetLineColor(kRed+2);
-  real.SetLineColor(kBlue+2);
-  rela.SetXTitle("#Psi_{R}");
-  rela.SetYTitle("#Psi_{2}");
-
-  std::cout << std::endl;
-  for (UInt_t i = 0; i < n_events; i++) {
+  Tester(Float_t v2=0.05, Int_t seg=-1, UInt_t n_max=20000)
+    : fGenerator(-1, 0.00, v2, n_max, n_max),
+      fAxis(10, -5, 5), 
+      fFlow("flow", "Analysed", 2, 1, fAxis), 
+      fReal("true", "Truth", 2, fAxis), 
+      fHist("hist", "Histogram", fAxis.N(), fAxis.Bins(),
+           (seg > 0 ? seg : 90), 0, 2*TMath::Pi()),
+      fRela("rela", "Truth", (seg > 0 ? seg : 90), 0, 2*TMath::Pi()),
+      fDPhi(seg <= 0 ? 0 : 2 * TMath::Pi() / seg), 
+      fPhis(n_max), 
+      fTruePhis(n_max), 
+      fXs(n_max)
+  {
+    fFlow.SetLineColor(kRed+2);
+    fReal.SetLineColor(kBlue+2);
+    fRela.SetXTitle("#Psi_{R}");
+    fRela.SetYTitle("#Psi_{2}");
+  }
+  void Run(UInt_t n_events) 
+  {
+    for (UInt_t i = 0; i < n_events; i++) Event(i);
+    std::cout << std::endl;
+    PrintResult();
+    DrawResult(n_events);
+    StoreResult();
+  }
+  void Event(Int_t i)
+  {
     std::cout << "\rEvent # " << i << std::flush;
     // Generate an event, get the number of objects, get the true
     // event plane, shuffle the phis around randomly. 
-    UInt_t   n_obs  = generator.Prepare();
-    Double_t rpsi_r = NormalizeAngle(generator.Psi());
+    UInt_t   n_obs  = fGenerator.Prepare();
+    Double_t rpsi_r = NormalizeAngle(fGenerator.Psi());
     
-    real.SetPsi(rpsi_r);
+    fReal.SetPsi(rpsi_r);
     std::cout << " (" << n_obs << " observations)" << std::flush;
     for (UInt_t j = 0; j < n_obs; j++) {
-      if (j % 2000 == 0) std::cout << "." << std::flush;
-      Double_t x   = gRandom->Gaus(0, 3);
-      Double_t phi = generator();
-      tphis[j]     = phi;
-      if (seg >= 0) { 
-       Int_t iseg = Int_t(phi / dphi);
-       phi        = dphi * (iseg + .5);
-      }
-      phis[j]      = phi;
-      xs[j]        = x;
-      hist.Fill(x, NormalizeAngle(phi-rpsi_r));
-      rela.Fill(NormalizeAngle(phi-rpsi_r));
+      Observation(j, rpsi_r);
+    }
+    fFlow.Event(n_obs, fPhis.fArray,     fXs.fArray);
+    fReal.Event(n_obs, fTruePhis.fArray, fXs.fArray);
+  }
+  void Observation(Int_t j, Double_t rpsi_r)
+  {
+    if (j % 2000 == 0) std::cout << "." << std::flush;
+    Double_t x   = gRandom->Gaus(0, 3);
+    Double_t phi = fGenerator();
+    fTruePhis[j] = phi;
+    if (fDPhi != 0) { 
+      Int_t iseg = Int_t(phi / fDPhi);
+      phi        = fDPhi * (iseg + .5);
     }
-    flow.Event(n_obs, phis.fArray,  xs.fArray);
-    real.Event(n_obs, tphis.fArray, xs.fArray);
+    fPhis[j]     = phi;
+    fXs[j]       = x;
+    fHist.Fill(x, NormalizeAngle(phi-rpsi_r));
+    fRela.Fill(NormalizeAngle(phi-rpsi_r));
+  }
+  void PrintResult()
+  {
+    fFlow.Print("s");
+    fReal.Print("s");
   }
-  std::cout << std::endl;
-  flow.Print("s");
-  real.Print("s");
+  void DrawResult(Int_t n_events, Bool_t showHists=kFALSE)
+  {
+    gStyle->SetPalette(1);
+    gStyle->SetOptStat(0);
+    gStyle->SetOptTitle(0);
+    
+    TCanvas* c = new TCanvas("C", "TestFlow", 600, (showHists ? 800 : 600));
+    c->SetFillColor(0);
+    c->SetBorderSize(0);
+    c->SetBorderMode(0);
+    Double_t x2 = (showHists ? 0.5 : 1);
+    TPad* pflow = new TPad("pflow", "PFlow", 0.0, 0.5, x2, 1.0, 0, 0, 0);
+    TPad* pres  = new TPad("pres",  "PRes",  0.0, 0.0, x2, 0.5, 0, 0, 0);
+    TPad* phist = 0;
+    TPad* p2d   = 0;
+    if (showHists) {
+      phist = new TPad("phist", "PHist", x2, 0.5, 1.0, 1.0, 0, 0, 0);
+      p2d   = new TPad("p2d",   "P2d",   x2, 0.0, 1.0, 0.5, 0, 0, 0);
+    }
+    pflow->SetTopMargin(0.05);
+    pflow->SetBottomMargin(0.00);
+    pflow->SetRightMargin(0.05);
+    pres->SetTopMargin(0.00);
+    pres->SetRightMargin(0.05);
 
-  gStyle->SetPalette(1);
-  
-  TCanvas* c = new TCanvas("C");
-  c->SetFillColor(0);
-  c->Divide(2,2);
-  TVirtualPad* p = c->cd(1);
-  // p->Divide(1,2);
-  // p = p->cd(1);
-  p->SetFillColor(0);
-  flow.Draw("th:");
-  real.Draw("th:same");
-  // p = c->cd(1);
-  p = c->cd(2);
-  rela.Scale(1. / n_events);
-  rela.Draw("");
-  p = c->cd(3);
-  p->SetFillColor(0);
-  flow.Draw("tr:");
-  real.Draw("tr:same");
-  c->cd(4);
-  hist.Scale(1. / n_events);
-  hist.Draw("lego2");
-  TBrowser* b = new TBrowser("b");
-  b->Add(&flow);
-  b->Add(&real);
+    c->cd();
+    pflow->Draw();
+    pflow->cd();
+    TH1* fth =fFlow.MakeHistogram(AliFMDFlowBin::kTdr,AliFMDFlowBin::kHarmonic);
+    TH1* tth =fReal.MakeHistogram(AliFMDFlowBin::kTdr,AliFMDFlowBin::kHarmonic);
+    Double_t min = TMath::Min(fth->GetMinimum(), tth->GetMinimum());
+    Double_t max = TMath::Max(fth->GetMaximum(), tth->GetMaximum());
+    fth->SetMinimum(0.9 * min); 
+    tth->SetMinimum(0.9 * min); 
+    fth->SetMaximum(1.2 * max); 
+    tth->SetMaximum(1.2 * max); 
+    fth->Draw();
+    tth->Draw("same");
+    TLegend* l = pflow->BuildLegend(0.5, 0.7, 0.94, 0.94);
+    l->SetFillColor(0);
+    l->SetBorderSize(0);
+
+    c->cd();
+    pres->Draw();
+    pres->cd();
+    TH1* ftr = fFlow.MakeHistogram(AliFMDFlowBin::kTdr,
+                                  AliFMDFlowBin::kResolution);
+    TH1* ttr = fReal.MakeHistogram(AliFMDFlowBin::kTdr,
+                                  AliFMDFlowBin::kResolution);
+    min      = TMath::Min(ftr->GetMinimum(), ttr->GetMinimum());
+    max      = TMath::Max(ftr->GetMaximum(), ttr->GetMaximum());
+    ftr->SetMinimum(0.8 * min); 
+    ttr->SetMinimum(0.8 * min); 
+    ftr->SetMaximum(1.2 * max); 
+    ttr->SetMaximum(1.2 * max); 
+    ftr->Draw();
+    ttr->Draw("same");
+    l = pres->BuildLegend(0.5, 0.7, 0.94, 0.94);
+    l->SetFillColor(0);
+    l->SetBorderSize(0);
+
+    if (p2d) {
+      p2d->Draw();
+      p2d->cd();
+      p2d->SetTopMargin(0.05);
+      p2d->SetRightMargin(0.05);
+      p2d->SetFillColor(0);
+      fRela.Scale(1. / n_events);
+      fRela.DrawCopy("");
+    }
 
-  TFile* file = TFile::Open("flow_test.root", "RECREATE");
-  flow.Write();
-  real.Write();
-  file->Close();
-  delete file;
+    if (phist) { 
+      phist->Draw();
+      phist->cd();
+      phist->SetTopMargin(0.05);
+      phist->SetRightMargin(0.05);
+      fHist.Scale(1. / n_events);
+      fHist.DrawCopy("lego2");
+    }
+    c->Modified();
+    c->Update();
+    c->cd();
+  }
+  void BrowseResult()
+  {
+    TBrowser* b = new TBrowser("b");
+    b->Add(&fFlow);
+    b->Add(&fReal);
+  }
+  void StoreResult() 
+  {
+    TFile* file = TFile::Open("flow_test.root", "RECREATE");
+    fFlow.Write();
+    fReal.Write();
+    file->Close();
+    delete file;
+  }
+  Generator          fGenerator;   // Event plane generator
+  AliFMDFlowAxis     fAxis;        // Bins
+  AliFMDFlowBinned1D fFlow;        // Our histogram
+  AliFMDFlowTrue1D   fReal;        // Histogram of true values
+  TH2D               fHist;        // Histogram of data
+  TH1D               fRela;        // Histogram of resolution
+  Double_t           fDPhi;        // Phi segmentation
+  TArrayD            fPhis;        // Cache of phis
+  TArrayD            fTruePhis;    // Cache of true phis
+  TArrayD            fXs;          // Cache of X's
+};
+//____________________________________________________________________
+void
+TestFlow(UInt_t n_events=100, Int_t seg=-1, UInt_t n_max=20000,
+        Float_t v2=0.05)
+{
+  std::cout << "Flow:                   " << v2 << "\n"
+           << "# of events:            " << n_events << "\n" 
+           << "# of phi segments:      ";
+  if (seg < 0) std::cout << "none\n"; 
+  else         std::cout << seg << "\n";
+  std::cout << "Maximum # observations: " << n_max << std::endl;
+                
+  Tester t(v2, seg, n_max);
+  t.Run(n_events);
 }
 
 #ifndef __CINT__
+#include <sstream>
 #include <TApplication.h>
 /** Entry point for test program */
+void 
+usage(std::ostream& o, const char* argv0)
+{
+  o << "Usage: " << argv0 << " [OPTIONS]\n\n"
+    << "Options:\n"
+    << "   -n N_EVENT          Set the number of events\n"
+    << "   -s N_SEGMENTS       Set number of phi segments\n" 
+    << "   -o MAX_OBSERVATIONS Set maximum number of observations per event\n"
+    << "   -v FLOW             Set the flow\n"
+    << "   -h                  Show this help\n"
+    << std::endl;
+}
+
+template <typename T>
+void str2val(const char* str, T& val)
+{
+  std::stringstream s(str);
+  s >> val;
+}
+
 int
-main()
+main(int argc, char** argv)
 {
+  UInt_t  n_events = 100;
+  Int_t   seg      = -1;
+  UInt_t  n_max    = 20000;
+  Float_t v2      = 0.05;
+  for (int i = 1; i < argc; i++) { 
+    if (argv[i][0] == '-') { 
+      switch (argv[i][1]) { 
+      case 'h': usage(std::cout, argv[0]); return 0;
+      case 'n': str2val(argv[++i], n_events); break;
+      case 'o': str2val(argv[++i], n_max); break;
+      case 's': str2val(argv[++i], seg); break;
+      case 'v': str2val(argv[++i], v2); break;
+      default:
+       std::cerr << "Unknown option: " << argv[i] << std::endl;
+       return 1;
+      }
+    }
+  }
   TApplication app("app", 0, 0);
-  TestFlow();
+  TestFlow(n_events, seg, n_max, v2);
   app.Run();
 }
 
diff --git a/FMD/flow/TestFlowSimple.C b/FMD/flow/TestFlowSimple.C
new file mode 100644 (file)
index 0000000..ec88564
--- /dev/null
@@ -0,0 +1,31 @@
+Double_t GeneratePhi(Double_t psiR, Double_t v1, Double_t v2)
+{
+  Double_t phi  =  gRandom->Uniform(0, 2 * TMath::Pi());
+  Double_t rel  =  phi - psiR;
+  Double_t dphi =  -2 * TMath::Sin(rel) * v1;
+  dphi          -= TMath::Sin(2 * rel) * v2;
+  phi           += dphi;
+  return phi;
+}
+
+void TestFlowSimple()
+{
+  gSystem->Load("libFMDflow.so");
+  AliFMDFlowAxis      axis(10, -5, 5);
+  AliFMDFlowBinned1D  flow("flow", "analysed", 2, 1, axis);
+  TArrayD             phis(20000);
+  TArrayD             xs(20000);
+
+  for (Int_t i = 0; i < 100; i++) { 
+    Double_t psiR = gRandom->Uniform(0, 2*TMath::Pi());
+    Int_t    nObs = gRandom->Integer(20000);
+    
+    for (Int_t j = 0; j < nObs; j++) { 
+      xs[j]   = gRandom->Uniform(-5, 5);
+      phis[j] = GeneratePhi(psiR, 0, 0.05);
+    }
+    flow.Event(nObs, phis.fArray, xs.fArray);
+  }
+  flow.Draw("ht");
+}
+