1 /* Copyright (C) 2007 Christian Holm Christensen <cholm@nbi.dk>
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.
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.
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
18 //___________________________________________________________________
24 Root> gSystem->Load("libFMDflow.so");
25 Root> gSystem->AddIncludePath("-I$ALICE_ROOT/FMD")
26 Root> .x FMD/flow/TestFlow.C
31 $ g++ `root-config --cflags --libs` -I$ALICE_ROOT/include \
32 -I$ALICE_ROOT/FMD $ALICE_ROOT/FMD/flow/TestFlow.C -o testflow
36 #include <FMD/flow/AliFMDFlowBinned1D.h>
37 #include <FMD/flow/AliFMDFlowTrue.h>
38 #include <FMD/flow/AliFMDFlowUtil.h>
41 #include <TVirtualPad.h>
50 //____________________________________________________________________
51 /** Generate events */
55 @param psi Possibly fixed event plane (if <0 -> random)
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)
65 /** Prepare for the next event.
66 @return Number of observations */
69 // Generate a uniform random direction
71 if (fPsi >= 0) fPsiR = fPsi;
72 else fPsiR = gRandom->Uniform(0, 2 * TMath::Pi());
73 return unsigned(gRandom->Uniform(fMin, fMax));
75 /** Create an observation */
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;
86 /** Get the event plane */
87 Double_t Psi() const { return fPsiR; }
97 /** Run test program */
99 TestFlow(UInt_t n_events=100, Int_t seg=-1, UInt_t n_max=20000,
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);
112 TArrayD tphis(n_max);
114 flow.SetLineColor(kRed+2);
115 real.SetLineColor(kBlue+2);
116 rela.SetXTitle("#Psi_{R}");
117 rela.SetYTitle("#Psi_{2}");
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());
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();
135 Int_t iseg = Int_t(phi / dphi);
136 phi = dphi * (iseg + .5);
140 hist.Fill(x, NormalizeAngle(phi-rpsi_r));
141 rela.Fill(NormalizeAngle(phi-rpsi_r));
143 flow.Event(n_obs, phis.fArray, xs.fArray);
144 real.Event(n_obs, tphis.fArray, xs.fArray);
146 std::cout << std::endl;
150 gStyle->SetPalette(1);
152 TCanvas* c = new TCanvas("C");
155 TVirtualPad* p = c->cd(1);
160 real.Draw("th:same");
163 rela.Scale(1. / n_events);
168 real.Draw("tr:same");
170 hist.Scale(1. / n_events);
172 TBrowser* b = new TBrowser("b");
176 TFile* file = TFile::Open("flow_test.root", "RECREATE");
184 #include <TApplication.h>
185 /** Entry point for test program */
189 TApplication app("app", 0, 0);