Fork me on GitHub

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

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

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