FMD/flow/TestFlow.C

   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  */
  18 //___________________________________________________________________
  19 /** @file
  20     @brief Example
  21
  22     Compile and run like
  23     @verbatim
  24     Root> gSystem->Load("libFMDflow.so");
  25     Root> gSystem->AddIncludePath("-I$ALICE_ROOT/FMD")
  26     Root> .x FMD/flow/TestFlow.C
  27     @endverbatim
  28
  29     or
  30     @verbatim
  31     $ g++ `root-config --cflags --libs` -I$ALICE_ROOT/include \
  32        -I$ALICE_ROOT/FMD $ALICE_ROOT/FMD/flow/TestFlow.C -o testflow
  33     $ ./testflow
  34     @endverbatim
  35 */
  36 #include <FMD/flow/AliFMDFlowBinned1D.h>
  37 #include <FMD/flow/AliFMDFlowTrue.h>
  38 #include <FMD/flow/AliFMDFlowUtil.h>
  39 #include <TRandom.h>
  40 #include <TCanvas.h>
  41 #include <TVirtualPad.h>
  42 #include <TArrayD.h>
  43 #include <TBrowser.h>
  44 #include <iostream>
  45 #include <TMath.h>
  46 #include <TH2.h>
  47 #include <TStyle.h>
  48 #include <TFile.h>
  49
  50 //____________________________________________________________________
  51 /** Generate events */
  52 struct Generator
  53 {
  54   /** Constructor
  55       @param psi Possibly fixed event plane (if <0 -> random)
  56       @param v1  value of v1
  57       @param v2  value of v2
  58       @param min Minimum number of observations
  59       @param max Maximum number of observations */
  60   Generator(Double_t psi=-1,
  61             Double_t v1=.05,  Double_t v2=.05,
  62             UInt_t   min=100, UInt_t   max=1000)
  63     : fPsi(psi), fV1(v1), fV2(v2), fMin(min), fMax(max)
  64   {}
  65   /** Prepare for the next event.
  66       @return Number of observations */
  67   UInt_t Prepare()
  68   {
  69     // Generate a uniform random direction
  70     fN = 0;
  71     if (fPsi >= 0) fPsiR = fPsi;
  72     else           fPsiR = gRandom->Uniform(0, 2 * TMath::Pi());
  73     return unsigned(gRandom->Uniform(fMin, fMax));
  74   }
  75   /** Create an observation */
  76   Double_t operator()()
  77   {
  78     // Generate a uniform random direction
  79     Double_t phi  =  gRandom->Uniform(0, 2 * TMath::Pi());
  80     Double_t rel  =  phi - fPsiR;
  81     Double_t dphi =  -2 * TMath::Sin(rel) * fV1;
  82     dphi          -= TMath::Sin(2 * rel) * fV2;
  83     phi           += dphi;
  84     return phi;
  85   }
  86   /** Get the event plane */
  87   Double_t Psi() const { return fPsiR; }
  88   Double_t fPsi;
  89   Double_t fPsiR;
  90   Double_t fV1;
  91   Double_t fV2;
  92   Double_t fMin;
  93   Double_t fMax;
  94   UInt_t   fN;
  95 };
  96
  97 /** Run test program */
  98 void
  99 TestFlow(UInt_t n_events=100, Int_t seg=-1, UInt_t n_max=20000,
 100          Float_t v2=0.05)
 101 {
 102   Generator           generator(-1, 0.00, v2, n_max, n_max);
 103   AliFMDFlowAxis      a(10, -5, 5);
 104   AliFMDFlowBinned1D  flow("test","Test",2, 1, a);
 105   AliFMDFlowTrue1D    real("true","true", 2, a);
 106   TH2D                hist("hist","hist",a.N(),a.Bins(),
 107                            (seg>0?seg:90),0, 2* TMath::Pi());
 108   TH1D                rela("rela","real",
 109                            (seg>0?seg:90),0, 2* TMath::Pi());
 110   Double_t            dphi = (seg <= 0 ? 0 : 2 * TMath::Pi() / seg);
 111   TArrayD             phis(n_max);
 112   TArrayD             tphis(n_max);
 113   TArrayD             xs(n_max);
 114   flow.SetLineColor(kRed+2);
 115   real.SetLineColor(kBlue+2);
 116   rela.SetXTitle("#Psi_{R}");
 117   rela.SetYTitle("#Psi_{2}");
 118
 119   std::cout << std::endl;
 120   for (UInt_t i = 0; i < n_events; i++) {
 121     std::cout << "\rEvent # " << i << std::flush;
 122     // Generate an event, get the number of objects, get the true
 123     // event plane, shuffle the phis around randomly.
 124     UInt_t   n_obs  = generator.Prepare();
 125     Double_t rpsi_r = NormalizeAngle(generator.Psi());
 126
 127     real.SetPsi(rpsi_r);
 128     std::cout << " (" << n_obs << " observations)" << std::flush;
 129     for (UInt_t j = 0; j < n_obs; j++) {
 130       if (j % 2000 == 0) std::cout << "." << std::flush;
 131       Double_t x   = gRandom->Gaus(0, 3);
 132       Double_t phi = generator();
 133       tphis[j]     = phi;
 134       if (seg >= 0) {
 135         Int_t iseg = Int_t(phi / dphi);
 136         phi        = dphi * (iseg + .5);
 137       }
 138       phis[j]      = phi;
 139       xs[j]        = x;
 140       hist.Fill(x, NormalizeAngle(phi-rpsi_r));
 141       rela.Fill(NormalizeAngle(phi-rpsi_r));
 142     }
 143     flow.Event(n_obs, phis.fArray,  xs.fArray);
 144     real.Event(n_obs, tphis.fArray, xs.fArray);
 145   }
 146   std::cout << std::endl;
 147   flow.Print("s");
 148   real.Print("s");
 149
 150   gStyle->SetPalette(1);
 151
 152   TCanvas* c = new TCanvas("C");
 153   c->SetFillColor(0);
 154   c->Divide(2,2);
 155   TVirtualPad* p = c->cd(1);
 156   // p->Divide(1,2);
 157   // p = p->cd(1);
 158   p->SetFillColor(0);
 159   flow.Draw("th:");
 160   real.Draw("th:same");
 161   // p = c->cd(1);
 162   p = c->cd(2);
 163   rela.Scale(1. / n_events);
 164   rela.Draw("");
 165   p = c->cd(3);
 166   p->SetFillColor(0);
 167   flow.Draw("tr:");
 168   real.Draw("tr:same");
 169   c->cd(4);
 170   hist.Scale(1. / n_events);
 171   hist.Draw("lego2");
 172   TBrowser* b = new TBrowser("b");
 173   b->Add(&flow);
 174   b->Add(&real);
 175
 176   TFile* file = TFile::Open("flow_test.root", "RECREATE");
 177   flow.Write();
 178   real.Write();
 179   file->Close();
 180   delete file;
 181 }
 182
 183 #ifndef __CINT__
 184 #include <TApplication.h>
 185 /** Entry point for test program */
 186 int
 187 main()
 188 {
 189   TApplication app("app", 0, 0);
 190   TestFlow();
 191   app.Run();
 192 }
 193
 194 #endif
 195
 196