Fork me on GitHub

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

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

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