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 -L$ALICE_ROOT/lib/tgt_${ALICE_TARGET} \
  33        -lFMDflow $ALICE_ROOT/FMD/flow/TestFlow.C -o testflow
  34     $ ./testflow
  35     @endverbatim
  36 */
  37 #include <flow/AliFMDFlowBinned1D.h>
  38 #include <flow/AliFMDFlowTrue.h>
  39 #include <flow/AliFMDFlowUtil.h>
  40 #include <TRandom.h>
  41 #include <TCanvas.h>
  42 #include <TVirtualPad.h>
  43 #include <TLegend.h>
  44 #include <TArrayD.h>
  45 #include <TBrowser.h>
  46 #include <iostream>
  47 #include <TMath.h>
  48 #include <TH2.h>
  49 #include <TStyle.h>
  50 #include <TFile.h>
  51
  52 //____________________________________________________________________
  53 /** Generate events */
  54 struct Generator
  55 {
  56   /** Constructor
  57       @param psi Possibly fixed event plane (if <0 -> random)
  58       @param v1  value of v1
  59       @param v2  value of v2
  60       @param min Minimum number of observations
  61       @param max Maximum number of observations */
  62   Generator(Double_t psi=-1,
  63             Double_t v1=.05,  Double_t v2=.05,
  64             UInt_t   min=100, UInt_t   max=1000)
  65     : fPsi(psi), fV1(v1), fV2(v2), fMin(min), fMax(max)
  66   {}
  67   /** Prepare for the next event.
  68       @return Number of observations */
  69   UInt_t Prepare()
  70   {
  71     // Generate a uniform random direction
  72     fN = 0;
  73     if (fPsi >= 0) fPsiR = fPsi;
  74     else           fPsiR = gRandom->Uniform(0, 2 * TMath::Pi());
  75     return unsigned(gRandom->Uniform(fMin, fMax));
  76   }
  77   /** Create an observation */
  78   Double_t operator()()
  79   {
  80     // Generate a uniform random direction
  81     Double_t phi  =  gRandom->Uniform(0, 2 * TMath::Pi());
  82     Double_t rel  =  phi - fPsiR;
  83     Double_t dphi =  -2 * TMath::Sin(rel) * fV1;
  84     dphi          -= TMath::Sin(2 * rel) * fV2;
  85     phi           += dphi;
  86     return phi;
  87   }
  88   /** Get the event plane */
  89   Double_t Psi() const { return fPsiR; }
  90   Double_t fPsi;
  91   Double_t fPsiR;
  92   Double_t fV1;
  93   Double_t fV2;
  94   Double_t fMin;
  95   Double_t fMax;
  96   UInt_t   fN;
  97 };
  98
  99 struct Tester
 100 {
 101   Tester(Float_t v2=0.05, Int_t seg=-1, UInt_t n_max=20000)
 102     : fGenerator(-1, 0.00, v2, n_max, n_max),
 103       fAxis(10, -5, 5),
 104       fFlow("flow", "Analysed", 2, 1, fAxis),
 105       fReal("true", "Truth", 2, fAxis),
 106       fHist("hist", "Histogram", fAxis.N(), fAxis.Bins(),
 107             (seg > 0 ? seg : 90), 0, 2*TMath::Pi()),
 108       fRela("rela", "Truth", (seg > 0 ? seg : 90), 0, 2*TMath::Pi()),
 109       fDPhi(seg <= 0 ? 0 : 2 * TMath::Pi() / seg),
 110       fPhis(n_max),
 111       fTruePhis(n_max),
 112       fXs(n_max)
 113   {
 114     fFlow.SetLineColor(kRed+2);
 115     fReal.SetLineColor(kBlue+2);
 116     fRela.SetXTitle("#Psi_{R}");
 117     fRela.SetYTitle("#Psi_{2}");
 118   }
 119   void Run(UInt_t n_events)
 120   {
 121     for (UInt_t i = 0; i < n_events; i++) Event(i);
 122     std::cout << std::endl;
 123     PrintResult();
 124     DrawResult(n_events);
 125     StoreResult();
 126   }
 127   void Event(Int_t i)
 128   {
 129     std::cout << "\rEvent # " << i << std::flush;
 130     // Generate an event, get the number of objects, get the true
 131     // event plane, shuffle the phis around randomly.
 132     UInt_t   n_obs  = fGenerator.Prepare();
 133     Double_t rpsi_r = NormalizeAngle(fGenerator.Psi());
 134
 135     fReal.SetPsi(rpsi_r);
 136     std::cout << " (" << n_obs << " observations)" << std::flush;
 137     for (UInt_t j = 0; j < n_obs; j++) {
 138       Observation(j, rpsi_r);
 139     }
 140     fFlow.Event(n_obs, fPhis.fArray,     fXs.fArray);
 141     fReal.Event(n_obs, fTruePhis.fArray, fXs.fArray);
 142   }
 143   void Observation(Int_t j, Double_t rpsi_r)
 144   {
 145     if (j % 2000 == 0) std::cout << "." << std::flush;
 146     Double_t x   = gRandom->Gaus(0, 3);
 147     Double_t phi = fGenerator();
 148     fTruePhis[j] = phi;
 149     if (fDPhi != 0) {
 150       Int_t iseg = Int_t(phi / fDPhi);
 151       phi        = fDPhi * (iseg + .5);
 152     }
 153     fPhis[j]     = phi;
 154     fXs[j]       = x;
 155     fHist.Fill(x, NormalizeAngle(phi-rpsi_r));
 156     fRela.Fill(NormalizeAngle(phi-rpsi_r));
 157   }
 158   void PrintResult()
 159   {
 160     fFlow.Print("s");
 161     fReal.Print("s");
 162   }
 163   void DrawResult(Int_t n_events, Bool_t showHists=kFALSE)
 164   {
 165     gStyle->SetPalette(1);
 166     gStyle->SetOptStat(0);
 167     gStyle->SetOptTitle(0);
 168
 169     TCanvas* c = new TCanvas("C", "TestFlow", 600, (showHists ? 800 : 600));
 170     c->SetFillColor(0);
 171     c->SetBorderSize(0);
 172     c->SetBorderMode(0);
 173     Double_t x2 = (showHists ? 0.5 : 1);
 174     TPad* pflow = new TPad("pflow", "PFlow", 0.0, 0.5, x2, 1.0, 0, 0, 0);
 175     TPad* pres  = new TPad("pres",  "PRes",  0.0, 0.0, x2, 0.5, 0, 0, 0);
 176     TPad* phist = 0;
 177     TPad* p2d   = 0;
 178     if (showHists) {
 179       phist = new TPad("phist", "PHist", x2, 0.5, 1.0, 1.0, 0, 0, 0);
 180       p2d   = new TPad("p2d",   "P2d",   x2, 0.0, 1.0, 0.5, 0, 0, 0);
 181     }
 182     pflow->SetTopMargin(0.05);
 183     pflow->SetBottomMargin(0.00);
 184     pflow->SetRightMargin(0.05);
 185     pres->SetTopMargin(0.00);
 186     pres->SetRightMargin(0.05);
 187
 188     c->cd();
 189     pflow->Draw();
 190     pflow->cd();
 191     TH1* fth =fFlow.MakeHistogram(AliFMDFlowBin::kTdr,AliFMDFlowBin::kHarmonic);
 192     TH1* tth =fReal.MakeHistogram(AliFMDFlowBin::kTdr,AliFMDFlowBin::kHarmonic);
 193     Double_t min = TMath::Min(fth->GetMinimum(), tth->GetMinimum());
 194     Double_t max = TMath::Max(fth->GetMaximum(), tth->GetMaximum());
 195     fth->SetMinimum(0.9 * min);
 196     tth->SetMinimum(0.9 * min);
 197     fth->SetMaximum(1.2 * max);
 198     tth->SetMaximum(1.2 * max);
 199     fth->Draw();
 200     tth->Draw("same");
 201     TLegend* l = pflow->BuildLegend(0.5, 0.7, 0.94, 0.94);
 202     l->SetFillColor(0);
 203     l->SetBorderSize(0);
 204
 205     c->cd();
 206     pres->Draw();
 207     pres->cd();
 208     TH1* ftr = fFlow.MakeHistogram(AliFMDFlowBin::kTdr,
 209                                    AliFMDFlowBin::kResolution);
 210     TH1* ttr = fReal.MakeHistogram(AliFMDFlowBin::kTdr,
 211                                    AliFMDFlowBin::kResolution);
 212     min      = TMath::Min(ftr->GetMinimum(), ttr->GetMinimum());
 213     max      = TMath::Max(ftr->GetMaximum(), ttr->GetMaximum());
 214     ftr->SetMinimum(0.8 * min);
 215     ttr->SetMinimum(0.8 * min);
 216     ftr->SetMaximum(1.2 * max);
 217     ttr->SetMaximum(1.2 * max);
 218     ftr->Draw();
 219     ttr->Draw("same");
 220     l = pres->BuildLegend(0.5, 0.7, 0.94, 0.94);
 221     l->SetFillColor(0);
 222     l->SetBorderSize(0);
 223
 224     if (p2d) {
 225       p2d->Draw();
 226       p2d->cd();
 227       p2d->SetTopMargin(0.05);
 228       p2d->SetRightMargin(0.05);
 229       p2d->SetFillColor(0);
 230       fRela.Scale(1. / n_events);
 231       fRela.DrawCopy("");
 232     }
 233
 234     if (phist) {
 235       phist->Draw();
 236       phist->cd();
 237       phist->SetTopMargin(0.05);
 238       phist->SetRightMargin(0.05);
 239       fHist.Scale(1. / n_events);
 240       fHist.DrawCopy("lego2");
 241     }
 242     c->Modified();
 243     c->Update();
 244     c->cd();
 245   }
 246   void BrowseResult()
 247   {
 248     TBrowser* b = new TBrowser("b");
 249     b->Add(&fFlow);
 250     b->Add(&fReal);
 251   }
 252   void StoreResult()
 253   {
 254     TFile* file = TFile::Open("flow_test.root", "RECREATE");
 255     fFlow.Write();
 256     fReal.Write();
 257     file->Close();
 258     delete file;
 259   }
 260   Generator          fGenerator;   // Event plane generator
 261   AliFMDFlowAxis     fAxis;        // Bins
 262   AliFMDFlowBinned1D fFlow;        // Our histogram
 263   AliFMDFlowTrue1D   fReal;        // Histogram of true values
 264   TH2D               fHist;        // Histogram of data
 265   TH1D               fRela;        // Histogram of resolution
 266   Double_t           fDPhi;        // Phi segmentation
 267   TArrayD            fPhis;        // Cache of phis
 268   TArrayD            fTruePhis;    // Cache of true phis
 269   TArrayD            fXs;          // Cache of X's
 270 };
 271 //____________________________________________________________________
 272 void
 273 TestFlow(UInt_t n_events=100, Int_t seg=-1, UInt_t n_max=20000,
 274          Float_t v2=0.05)
 275 {
 276   std::cout << "Flow:                   " << v2 << "\n"
 277             << "# of events:            " << n_events << "\n"
 278             << "# of phi segments:      ";
 279   if (seg < 0) std::cout << "none\n";
 280   else         std::cout << seg << "\n";
 281   std::cout << "Maximum # observations: " << n_max << std::endl;
 282
 283   Tester t(v2, seg, n_max);
 284   t.Run(n_events);
 285 }
 286
 287 #ifndef __CINT__
 288 #include <sstream>
 289 #include <TApplication.h>
 290 /** Entry point for test program */
 291 void
 292 usage(std::ostream& o, const char* argv0)
 293 {
 294   o << "Usage: " << argv0 << " [OPTIONS]\n\n"
 295     << "Options:\n"
 296     << "   -n N_EVENT          Set the number of events\n"
 297     << "   -s N_SEGMENTS       Set number of phi segments\n"
 298     << "   -o MAX_OBSERVATIONS Set maximum number of observations per event\n"
 299     << "   -v FLOW             Set the flow\n"
 300     << "   -h                  Show this help\n"
 301     << std::endl;
 302 }
 303
 304 template <typename T>
 305 void str2val(const char* str, T& val)
 306 {
 307   std::stringstream s(str);
 308   s >> val;
 309 }
 310
 311 int
 312 main(int argc, char** argv)
 313 {
 314   UInt_t  n_events = 100;
 315   Int_t   seg      = -1;
 316   UInt_t  n_max    = 20000;
 317   Float_t v2      = 0.05;
 318   for (int i = 1; i < argc; i++) {
 319     if (argv[i][0] == '-') {
 320       switch (argv[i][1]) {
 321       case 'h': usage(std::cout, argv[0]); return 0;
 322       case 'n': str2val(argv[++i], n_events); break;
 323       case 'o': str2val(argv[++i], n_max); break;
 324       case 's': str2val(argv[++i], seg); break;
 325       case 'v': str2val(argv[++i], v2); break;
 326       default:
 327         std::cerr << "Unknown option: " << argv[i] << std::endl;
 328         return 1;
 329       }
 330     }
 331   }
 332   TApplication app("app", 0, 0);
 333   TestFlow(n_events, seg, n_max, v2);
 334   app.Run();
 335 }
 336
 337 #endif
 338
 339