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
  49 //____________________________________________________________________
  50 /** Generate events */
  51 struct Generator
  52 {
  53   /** Constructor
  54       @param psi Possibly fixed event plane (if <0 -> random)
  55       @param v1  value of v1
  56       @param v2  value of v2
  57       @param min Minimum number of observations
  58       @param max Maximum number of observations */
  59   Generator(Double_t psi=-1,
  60             Double_t v1=.05,  Double_t v2=.05,
  61             UInt_t   min=100, UInt_t   max=1000)
  62     : fPsi(psi), fV1(v1), fV2(v2), fMin(min), fMax(max)
  63   {}
  64   /** Prepare for the next event.
  65       @return Number of observations */
  66   UInt_t Prepare()
  67   {
  68     // Generate a uniform random direction
  69     fN = 0;
  70     if (fPsi >= 0) fPsiR = fPsi;
  71     else           fPsiR = gRandom->Uniform(0, 2 * TMath::Pi());
  72     return unsigned(gRandom->Uniform(fMin, fMax));
  73   }
  74   /** Create an observation */
  75   Double_t operator()()
  76   {
  77     // Generate a uniform random direction
  78     Double_t phi  =  gRandom->Uniform(0, 2 * TMath::Pi());
  79     Double_t rel  =  phi - fPsiR;
  80     Double_t dphi =  -2 * TMath::Sin(rel) * fV1;
  81     dphi          -= TMath::Sin(2 * rel) * fV2;
  82     phi           += dphi;
  83     return phi;
  84   }
  85   /** Get the event plane */
  86   Double_t Psi() const { return fPsiR; }
  87   Double_t fPsi;
  88   Double_t fPsiR;
  89   Double_t fV1;
  90   Double_t fV2;
  91   Double_t fMin;
  92   Double_t fMax;
  93   UInt_t   fN;
  94 };
  95
  96 /** Run test program */
  97 void
  98 TestFlow(UInt_t n_events=100, Int_t seg=-1, UInt_t n_max=20000)
  99 {
 100   Generator           generator(-1, 0.00, 0.05, n_max, n_max);
 101   AliFMDFlowAxis      a(10, -5, 5);
 102   AliFMDFlowBinned1D* flow = new AliFMDFlowBinned1D(2, a);
 103   AliFMDFlowTrue1D*   real = new AliFMDFlowTrue1D(2, a);
 104   TH2D*               hist = new TH2D("hist","hist",a.N(),a.Bins(),
 105                                       (seg>0?seg:90),0, 2* TMath::Pi());
 106   Double_t            dphi = (seg <= 0 ? 0 : 2 * TMath::Pi() / seg);
 107   TArrayD             phis;
 108   TArrayD             tphis;
 109   TArrayD             xs;
 110
 111   std::cout << std::endl;
 112   for (UInt_t i = 0; i < n_events; i++) {
 113     std::cout << "\rEvent # " << i << std::flush;
 114     // Generate an event, get the number of objects, get the true
 115     // event plane, shuffle the phis around randomly.
 116     UInt_t   n_obs  = generator.Prepare();
 117     Double_t rpsi_r = NormalizeAngle(generator.Psi());
 118     phis.Set(n_obs);
 119     phis.Reset(0);
 120     tphis.Set(n_obs);
 121     tphis.Reset(0);
 122     xs.Set(n_obs);
 123     xs.Reset(0);
 124
 125     real->SetPsi(rpsi_r);
 126     std::cout << " (" << n_obs << " observations)" << std::flush;
 127     for (UInt_t j = 0; j < n_obs; j++) {
 128       if (j % 2000 == 0) std::cout << "." << std::flush;
 129       Double_t x   = gRandom->Gaus(0, 3);
 130       Double_t phi = generator();
 131       tphis[j]     = phi;
 132       if (seg >= 0) {
 133         Int_t iseg = Int_t(phi / dphi);
 134         phi        = dphi * (iseg + .5);
 135       }
 136       hist->Fill(x, phi);
 137       phis[j]      = phi;
 138       xs[j]        = x;
 139     }
 140     flow->Event(phis.fArray,  xs.fArray, 0, phis.fN);
 141     real->Event(tphis.fArray, xs.fArray, 0, tphis.fN);
 142   }
 143   std::cout << std::endl;
 144   flow->Print("s");
 145   real->Print("s");
 146
 147   gStyle->SetPalette(1);
 148
 149   TCanvas* c = new TCanvas("C");
 150   c->SetFillColor(0);
 151   c->Divide(2,1);
 152   TVirtualPad* p = c->cd(1);
 153   p->Divide(1,2);
 154   p = p->cd(1);
 155   p->SetFillColor(0);
 156   flow->Draw("bnst colz");
 157   p = c->cd(1);
 158   p = p->cd(2);
 159   p->SetFillColor(0);
 160   flow->Draw("bnstr colz");
 161   c->cd(2);
 162   hist->Draw("colz");
 163   new TBrowser("b", flow);
 164 }
 165
 166 #ifndef __CINT__
 167 #include <TApplication.h>
 168 /** Entry point for test program */
 169 int
 170 main()
 171 {
 172   TApplication app("app", 0, 0);
 173   TestFlow();
 174   app.Run();
 175 }
 176
 177 #endif
 178
 179