Fork me on GitHub

source: svn/trunk/Utilities/Hector/src/H_RecRPObject.cc@ 19

Last change on this file since 19 was 3, checked in by Xavier Rouby, 16 years ago

first commit
File size: 10.1 KB
Line 
1/*
2---- Hector the simulator ----
3 A fast simulator of particles through generic beamlines.
4 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be
5
6 http://www.fynu.ucl.ac.be/hector.html
7
8 Centre de Physique des Particules et de Phénoménologie (CP3)
9 Université Catholique de Louvain (UCL)
10*/
11
12
13/// \file H_RecRPObject.cc
14/// \brief Class performing the reconstruction based on forward detector measurements
15///
16/// Units : angles [rad], distances [m], energies [GeV], c=[1].
17
18// local #includes
19#include "H_RecRPObject.h"
20#include "H_RomanPot.h"
21#include "H_BeamParticle.h"
22using namespace std;
23
24// reconstruction class for forward detectors.
25// Featuring the brand-new reco method from the
26// louvain group !
27
28#define MEGA 1000000.
29
30H_RecRPObject::H_RecRPObject(): emin(0), emax(-1), x1(0), x2(0), y1(0), y2(0), s1(0), s2(0),
31 txip(NOT_YET_COMPUTED), tyip(NOT_YET_COMPUTED), energy(NOT_YET_COMPUTED), q2(NOT_YET_COMPUTED), pt(NOT_YET_COMPUTED),
32 thebeam(new H_AbstractBeamLine()),
33 f_1(new TF1("f_1","[0] + [1]*x + [2]*x*x ",emin,emax)),
34 f_2(new TF1("f_2","[0] + [1]*x + [2]*x*x ",emin,emax)),
35 g_1(new TF1("g_1","[0] + [1]*x + [2]*x*x ",emin,emax)),
36 g_2(new TF1("g_2","[0] + [1]*x + [2]*x*x ",emin,emax)),
37 d_1(new TF1("d_1","[0] + [1]*x + [2]*x*x ",emin,emax)),
38 d_2(new TF1("d_2","[0] + [1]*x + [2]*x*x ",emin,emax)),
39 k_1(new TF1("k_1","[0] + [1]*x + [2]*x*x ",emin,emax)),
40 k_2(new TF1("k_2","[0] + [1]*x + [2]*x*x ",emin,emax)),
41 l_1(new TF1("l_1","[0] + [1]*x + [2]*x*x ",emin,emax)),
42 l_2(new TF1("l_2","[0] + [1]*x + [2]*x*x ",emin,emax))
43 {
44 // emax = -1; /* GeV */
45}
46
47H_RecRPObject::H_RecRPObject(float ss1, float ss2, H_AbstractBeamLine* beam) : emin(0), emax(-1), x1(0), x2(0), y1(0), y2(0), s1(ss1), s2(ss2),
48 txip(NOT_YET_COMPUTED), tyip(NOT_YET_COMPUTED), energy(NOT_YET_COMPUTED), q2(NOT_YET_COMPUTED), pt(NOT_YET_COMPUTED),
49 thebeam(beam),
50 f_1(new TF1("f_1","[0] + [1]*x + [2]*x*x ",emin,emax)),
51 f_2(new TF1("f_2","[0] + [1]*x + [2]*x*x ",emin,emax)),
52 g_1(new TF1("g_1","[0] + [1]*x + [2]*x*x ",emin,emax)),
53 g_2(new TF1("g_2","[0] + [1]*x + [2]*x*x ",emin,emax)),
54 d_1(new TF1("d_1","[0] + [1]*x + [2]*x*x ",emin,emax)),
55 d_2(new TF1("d_2","[0] + [1]*x + [2]*x*x ",emin,emax)),
56 k_1(new TF1("k_1","[0] + [1]*x + [2]*x*x ",emin,emax)),
57 k_2(new TF1("k_2","[0] + [1]*x + [2]*x*x ",emin,emax)),
58 l_1(new TF1("l_1","[0] + [1]*x + [2]*x*x ",emin,emax)),
59 l_2(new TF1("l_2","[0] + [1]*x + [2]*x*x ",emin,emax))
60 {
61 if(ss1==ss2) cout<<"<H_RecRPObject> WARNING : detectors are on same position"<<endl;
62 // emax = -1; /* GeV */
63}
64
65void H_RecRPObject::initialize() {
66 // this method sets the functions that will be used for reco later
67 // it should be used only once per beamline after the energy range was fixed.
68 // copying beamline and adding detectors
69
70 if(emax<0) {
71 cout<<"<H_RecRPObject> ERROR : energy range has to be set first !"<<endl;
72 cout<<"<H_RecRPObject> Please run setERange() or computeERange()"<<endl;
73 cout<<"<H_RecRPObject> initialization aborted"<<endl;
74 return;
75 }
76
77 if(emax<emin) {
78 cout<<"<H_RecRPObject> ERROR : maximum energy lower than minimum !"<<endl;
79 cout<<"<H_RecRPObject> Please (re-)do setERange()"<<endl;
80 cout<<"<H_RecRPObject> initialization aborted"<<endl;
81 return;
82 }
83
84 H_AbstractBeamLine * b1 = thebeam->clone();
85 H_RomanPot * rp1 = new H_RomanPot("rp1",s1,0);
86 H_RomanPot * rp2 = new H_RomanPot("rp2",s2,0);
87 b1->add(rp1);
88 b1->add(rp2);
89
90 // fitting parameters
91 const int N = 20;
92 double e_i[N];
93 double f_1i[N], f_2i[N], g_1i[N], g_2i[N], d_1i[N], d_2i[N];
94 double k_1i[N], k_2i[N], l_1i[N], l_2i[N];
95 for(int i = 0; i < N; i++) {
96 e_i[i] = emin + i * (emax - emin)/((double)N-1);
97 //
98 // the bug seems to be linked to the delete operator of the TMatrixT class in root.
99 // valgrind shows memory problems at that point.
100 // for unknwown reasons, copying the matrix gets around this bug.
101 // valgrind (related) ouptut :
102 //
103 // ==13029== Invalid read of size 4
104 // ==13029== at 0x5E75DB6: H_RecRPObject::initialize() (in /home/jdf/GGamma/Hector/lib/libHector.so)
105 // ==13029== by 0x804A2D8: intelligentreco_rpo(double, double, double, double, std::string, int) (H_IntelligentReco.cpp:314)
106 // ==13029== by 0x804A937: main (H_IntelligentReco.cpp:542)
107 // ==13029== Address 0x64AC740 is 80 bytes inside a block of size 144 free'd
108 // ==13029== at 0x4021D18: operator delete[](void*) (vg_replace_malloc.c:256)
109 // ==13029== by 0x5651F57: TMatrixT<float>::Delete_m(int, float*&) (in /home/jdf/root/5.12/lib/libMatrix.so)
110 // ==13029== by 0x565CC5D: TMatrixT<float>::~TMatrixT() (in /home/jdf/root/5.12/lib/libMatrix.so)
111 // ==13029== by 0x5E75D25: H_RecRPObject::initialize() (in /home/jdf/GGamma/Hector/lib/libHector.so)
112 //
113 //
114 TMatrix el_mattt(b1->getPartialMatrix("rp1",e_i[i],MP,QP));
115 const float *el_mat1 = el_mattt.GetMatrixArray();
116 //
117 // conclusion : first line of el_mat1 is completely messed-up if it is taken directly from
118 // the return of getpartialmatrix like this :
119 // const float *el_mat1 = (b1->getPartialMatrix("rp1",e_i[i],MP,QP)).GetMatrixArray()
120 //
121 // long-term solution (apart noticing root staff) : replacing el_mat1[i] by the equivalent el_mattt(j,k)
122 // which is anyway more transparent for the reader.
123 //
124 f_1i[i] = el_mat1[0];
125 g_1i[i] = el_mat1[1*MDIM];
126 d_1i[i] = MEGA*el_mat1[4*MDIM];
127 k_1i[i] = el_mat1[2*MDIM+2];
128 l_1i[i] = el_mat1[3*MDIM+2];
129 const float* el_mat2 = (b1->getPartialMatrix("rp2",e_i[i],MP,QP).GetMatrixArray());
130 f_2i[i] = el_mat2[0];
131 g_2i[i] = el_mat2[1*MDIM];
132 d_2i[i] = MEGA*el_mat2[4*MDIM];
133 k_2i[i] = el_mat2[2*MDIM+2];
134 l_2i[i] = el_mat2[3*MDIM+2];
135 }
136 TGraph gf_1(N,e_i,f_1i);
137 TGraph gg_1(N,e_i,g_1i);
138 TGraph gd_1(N,e_i,d_1i);
139 TGraph gf_2(N,e_i,f_2i);
140 TGraph gg_2(N,e_i,g_2i);
141 TGraph gd_2(N,e_i,d_2i);
142 TGraph gk_1(N,e_i,k_1i);
143 TGraph gl_1(N,e_i,l_1i);
144 TGraph gk_2(N,e_i,k_2i);
145 TGraph gl_2(N,e_i,l_2i);
146
147 // functions get their final shape
148 gf_1.Fit("f_1","Q");
149 gg_1.Fit("g_1","Q");
150 gd_1.Fit("d_1","Q");
151 gf_2.Fit("f_2","Q");
152 gg_2.Fit("g_2","Q");
153 gd_2.Fit("d_2","Q");
154 gk_1.Fit("k_1","Q");
155 gl_1.Fit("l_1","Q");
156 gk_2.Fit("k_2","Q");
157 gl_2.Fit("l_2","Q");
158
159 // cleaning the rest
160 delete b1;
161
162 // the end
163 return;
164}
165
166void H_RecRPObject::setDetPos(const float ss1, const float ss2) {
167 energy = NOT_YET_COMPUTED;
168 s1 = ss1;
169 s2 = ss2;
170 if(ss1==ss2) cout<<"<H_RecRPObject> WARNING : detectors are on same position"<<endl;
171 return;
172}
173
174void H_RecRPObject::setPositions(const float xx1, const float xx2, const float yy1, const float yy2) {
175 energy = NOT_YET_COMPUTED;
176 x1 = xx1;
177 x2 = xx2;
178 y1 = yy1;
179 y2 = yy2;
180 return;
181}
182
183void H_RecRPObject::setPosition_det1(const float xx1, const float yy1) {
184 energy = NOT_YET_COMPUTED;
185 x1 = xx1;
186 y1 = yy1;
187}
188
189void H_RecRPObject::setPosition_det2(const float xx2, const float yy2) {
190 energy = NOT_YET_COMPUTED;
191 x2 = xx2;
192 y2 = yy2;
193}
194
195void H_RecRPObject::setERange(const float eemin, const float eemax) {
196 energy = NOT_YET_COMPUTED;
197 emin = eemin;
198 emax = eemax;
199 f_1->SetRange(emin,emax);
200 f_2->SetRange(emin,emax);
201 g_1->SetRange(emin,emax);
202 g_2->SetRange(emin,emax);
203 d_1->SetRange(emin,emax);
204 d_2->SetRange(emin,emax);
205 k_1->SetRange(emin,emax);
206 k_2->SetRange(emin,emax);
207 l_1->SetRange(emin,emax);
208 l_2->SetRange(emin,emax);
209 return;
210}
211
212
213void H_RecRPObject::computeERange() {
214 // optional method to determine the energy range of the FIRST detector
215 // in order to refine the fits and get maximum precision.
216 energy = NOT_YET_COMPUTED;
217 H_AbstractBeamLine * b1 = thebeam->clone();
218 H_RomanPot * rp1 = new H_RomanPot("rp1",s1,0);
219 b1->add(rp1);
220 float max = 1;
221 // number of energies to check
222 const int N = 1000;
223 for(int i=0; i<N; i++) {
224 H_BeamParticle p;
225 p.setE(BE - (emin + i*(BE-emin)/((float)N)));
226 p.computePath(b1);
227 if(p.stopped(b1)) {
228 if(p.getStoppingElement()->getName()=="rp1") {
229 max = emin + i*(BE-emin)/((float)N);
230 }
231 }
232 }
233 cout<<"<H_RecRPObject> Valid energy losses run from 0 (default) to "<<max+20.<<" GeV"<<endl;
234 setERange(0,max+20.);
235 delete b1;
236 return;
237}
238
239void H_RecRPObject::computeAll() {
240 // The big game :
241 // computing E, tx, ty, Q2 and Pt and filling the variables.
242 //
243 // The root TF1 class features nice bugs, which explains the
244 // seemingly-dumb structures happening sometimes here as
245 // workarounds for these bugs. The overall thing works very
246 // well but will be cleaned later anyway.
247
248 if(energy!=NOT_YET_COMPUTED) {
249 cout<<"<H_RecRPObject> already computed variables, skipping ..."<<endl;
250 return;
251 }
252
253 TF1 par0("par0","[0]",emin,emax);
254 par0.SetParameter(0,-x1);
255 TF1 par2("par2","[0]",emin,emax);
256 par2.SetParameter(0,-y1);
257 TF1 par1("par1","[0]",emin,emax);
258 par1.SetParameter(0,-x2);
259 TF1 par3("par3","[0]",emin,emax);
260 par3.SetParameter(0,-y2);
261
262 // angle compensating method :
263 TF1 xx_E("xx_E","(g_2*(par0-d_1*x)-g_1*(par1-d_2*x))/(f_2*g_1-f_1*g_2)",emin,emax);
264 TF1 yy_E("yy_E","(par2*l_2 - par3*l_1) / (k_2*l_1 - k_1*l_2)",emin,emax);
265 TF1 xp_E("xp_E","(f_2*(par0-d_1*x)-f_1*(par1-d_2*x))/(g_2*f_1-g_1*f_2)",emin,emax);
266 TF1 yp_E("yp_E","(par2*k_2-par3*k_1)/(l_2*k_1-l_1*k_2)",emin,emax);
267 // it is possible to refine study using y info, but effect was not tested.
268 // TF1 p_xy_E("p_xy_E","(-xx_E*xx_E-yy_E*yy_E)",emin,emax);
269 TF1 p_xy_E("p_xy_E","(-xx_E*xx_E)",emin,emax);
270
271 energy = p_xy_E.GetMaximumX(emin,emax);
272 txip = xp_E.Eval(energy);
273 tyip = yp_E.Eval(energy);
274 pt = sqrt(BE*(BE-energy)*(txip*txip+tyip*tyip)/(MEGA*MEGA));
275
276 return;
277}
278
279float H_RecRPObject::getE(int a) {
280// put for backward compatibility
281 if(energy==NOT_YET_COMPUTED) { computeAll(); };
282 return energy;
283} // to be removed !!!!!
284
285float H_RecRPObject::getE() {
286 if(energy==NOT_YET_COMPUTED) { computeAll(); };
287 return energy;
288}
289
290float H_RecRPObject::getTX() {
291 if(energy==NOT_YET_COMPUTED) { computeAll(); };
292 return txip;
293}
294
295float H_RecRPObject::getTY() {
296 if(energy==NOT_YET_COMPUTED) { computeAll(); };
297 return tyip;
298}
299
300float H_RecRPObject::getQ2() {
301 if(energy==NOT_YET_COMPUTED) { computeAll(); };
302 cout<<"<H_RecRPObject::getQ2> Not implemented yet"<<endl;
303 return 0;
304}
305
306float H_RecRPObject::getPt() {
307 if(energy==NOT_YET_COMPUTED) { computeAll(); };
308 return pt;
309}
Note: See TracBrowser for help on using the repository browser.