Fork me on GitHub

source: git/external/Hector/H_RecRPObject.cc@ b5bcacb

ImprovedOutputFile
Last change on this file since b5bcacb was 3c40083, checked in by pavel <pavel@…>, 11 years ago

switch to a more stable Hector version
  • Property mode set to 100644
File size: 10.8 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/// \file H_RecRPObject.cc
13/// \brief Classes aiming at reconstruction particle properties
14
15// C++ #includes
16#include <iostream>
17#include <iomanip>
18
19// ROOT #includes
20//#include "TGraph.h"
21//#include "TF1.h"
22//#include "TCanvas.h"
23
24// local #includes
25#include "H_RecRPObject.h"
26#include "H_RomanPot.h"
27#include "H_BeamParticle.h"
28using namespace std;
29
30H_RecRPObject::H_RecRPObject() {
31 x1 = 0.;
32 x2 = 0.;
33 y1 = 0.;
34 y2 = 0.;
35 s1 = 0.;
36 s2 = 0.;
37 corr1_TM = 0;
38 corr2_TM = 0;
39 corr1_AM = 0;
40 corr2_AM = 0;
41 thx = NOT_YET_COMPUTED;
42 thy = NOT_YET_COMPUTED;
43 x0 = NOT_YET_COMPUTED;
44 y0 = NOT_YET_COMPUTED;
45 energy = NOT_YET_COMPUTED;
46 virtuality = NOT_YET_COMPUTED;
47 matrp1 = new TMatrix(MDIM,MDIM);
48 matrp2 = new TMatrix(MDIM,MDIM);
49 thebeam = new H_AbstractBeamLine();
50}
51
52H_RecRPObject::H_RecRPObject(const float S1, const float S2, const H_AbstractBeamLine& beamline) {
53 x1 = 0;
54 x2 = 0;
55 y1 = 0;
56 y2 = 0;
57 s1 = S1;
58 s2 = S2;
59 thx = NOT_YET_COMPUTED;
60 thy = NOT_YET_COMPUTED;
61 x0 = NOT_YET_COMPUTED;
62 y0 = NOT_YET_COMPUTED;
63 energy = NOT_YET_COMPUTED;
64 virtuality = NOT_YET_COMPUTED;
65// matrp1 = new TMatrix(MDIM,MDIM);
66// matrp2 = new TMatrix(MDIM,MDIM);
67 thebeam = new H_AbstractBeamLine(beamline);
68 H_RomanPot * rp1 = new H_RomanPot("rp1",s1,0);
69 thebeam->add(rp1);
70 H_RomanPot * rp2 = new H_RomanPot("rp2",s2,0);
71 thebeam->add(rp2);
72 matrp1 = new TMatrix(*(thebeam->getPartialMatrix("rp1",0.,MP,QP)));
73 matrp2 = new TMatrix(*(thebeam->getPartialMatrix("rp2",0.,MP,QP)));
74
75 corr1_TM = getECorrectionFactor(0,TM);
76 corr2_TM = getECorrectionFactor(1,TM);
77 corr1_AM = getECorrectionFactor(0,AM);
78 corr2_AM = getECorrectionFactor(1,AM);
79// cout << corr1_TM << " " << corr2_TM << endl;
80// cout << corr1_AM << " " << corr2_AM << endl;
81}
82
83H_RecRPObject::H_RecRPObject(const H_RecRPObject& r) {
84 x1 = r.x1;
85 x2 = r.x2;
86 y1 = r.y1;
87 y2 = r.y2;
88 s1 = r.s1;
89 s2 = r.s2;
90 x0 = r.x0;
91 y0 = r.y0;
92 thx = r.thx;
93 thy = r.thy;
94 energy = r.energy;
95 virtuality = r.virtuality;
96 matrp1 = new TMatrix(*(r.matrp1));
97 matrp2 = new TMatrix(*(r.matrp2));
98 corr1_TM = r.corr1_TM;
99 corr2_TM = r.corr2_TM;
100 corr1_AM = r.corr1_AM;
101 corr2_AM = r.corr2_AM;
102 thebeam = new H_AbstractBeamLine(*(r.thebeam));
103}
104
105H_RecRPObject& H_RecRPObject::operator=(const H_RecRPObject& r) {
106 if(this==&r) return *this;
107 x1 = r.x1;
108 x2 = r.x2;
109 y1 = r.y1;
110 y2 = r.y2;
111 s1 = r.s1;
112 s2 = r.s2;
113 x0 = r.x0;
114 y0 = r.y0;
115 thx = r.thx;
116 thy = r.thy;
117 energy = r.energy;
118 virtuality = r.virtuality;
119 matrp1 = new TMatrix(*(r.matrp1));
120 matrp2 = new TMatrix(*(r.matrp2));
121 corr1_TM = r.corr1_TM;
122 corr2_TM = r.corr2_TM;
123 corr1_AM = r.corr1_AM;
124 corr2_AM = r.corr2_AM;
125 thebeam = new H_AbstractBeamLine(*(r.thebeam));
126 return *this;
127}
128
129float H_RecRPObject::getECorrectionFactor(const unsigned int facn, const unsigned int method) {
130/*
131 * commented out because CMSSW does not want any TGraph/TCanvas/TFit
132 *
133 * to be fixed !
134 * X.R. 07/05/2009
135 *
136 float beta1 = ((thebeam->getPartialMatrix("rp1",0,MP,QP))->GetMatrixArray())[1*MDIM];
137 float beta2 = ((thebeam->getPartialMatrix("rp2",0,MP,QP))->GetMatrixArray())[1*MDIM];
138 float disp1 = ((thebeam->getPartialMatrix("rp1",0,MP,QP))->GetMatrixArray())[4*MDIM]*URAD;
139 float disp2 = ((thebeam->getPartialMatrix("rp2",0,MP,QP))->GetMatrixArray())[4*MDIM]*URAD;
140 const int n = 20; //using 20 points to get a good quadratic fit
141 float ee[n], rece[n];
142
143 for(int i = 0; i < n; i++) {
144 ee[i] = 10 + i*200./((float)n-1);
145 H_BeamParticle p1;
146 p1.emitGamma(ee[i],0.);
147 p1.computePath(thebeam,1);
148 p1.propagate(s1);
149 float x1 = p1.getX();
150 p1.propagate(s2);
151 float x2 = p1.getX();
152 switch (method) {
153 case TM: { rece[i] = -x1/disp1; }; break;
154 case AM: { rece[i] = -(beta2*x1-beta1*x2)/(beta2*disp1-beta1*disp2); }; break;
155 case PM: { rece[i] = -x1/disp1; cout<<"this method has not been implemented, using trivial reconstruction"<<endl; } break;
156 default: { rece[i] = -(beta2*x1-beta1*x2)/(beta2*disp1-beta1*disp2); }; break;
157 }
158 ee[i] = ee[i] - rece[i];
159 }
160 char mytitle[50];
161 sprintf(mytitle,"c_%d",method);
162// TCanvas*c = new TCanvas();
163// c->SetTitle(mytitle);
164 TGraph* g1 = new TGraph(n,rece,ee);
165 TF1* fit1 = new TF1("fit1","[0]*x + [1]*x*x",0,100);
166 g1->Fit("fit1","Q");
167 float xfact = fit1->GetParameter(facn);
168// g1->Draw("APL");
169 delete g1;
170 delete fit1;
171
172 return xfact;
173*/
174 return 1.;
175}
176
177
178void H_RecRPObject::setPositions(const float X1, const float Y1, const float X2, const float Y2) {
179 thx = NOT_YET_COMPUTED;
180 thy = NOT_YET_COMPUTED;
181 x0 = NOT_YET_COMPUTED;
182 y0 = NOT_YET_COMPUTED;
183 energy = NOT_YET_COMPUTED;
184 virtuality = NOT_YET_COMPUTED;
185 x1 = X1;
186 x2 = X2;
187 y1 = Y1;
188 y2 = Y2;
189
190 return;
191}
192
193void H_RecRPObject::printProperties() const {
194 cout << "Roman pot variables :" << endl;
195 cout << "\t pot 1 : (x,y,s) = (" << x1 << " , " << y1 << " , " << s1 << " )" << endl;
196 cout << "\t pot 2 : (x,y,s) = (" << x2 << " , " << y2 << " , " << s2 << " )" << endl;
197 cout << endl << "Reconstructed variables :" << endl;
198 cout << "\t IP : (x,y) = (" << x0 << " , " << y0 << ") and (theta_x, theta_y) = (" << thx << " , " << thy << " )" << endl;
199 if (energy==NOT_YET_COMPUTED) cout << "\t Energy not yet computed" << endl;
200 else cout << "\t Energy = " << energy << " GeV" << endl;
201 if (virtuality==NOT_YET_COMPUTED) cout << "\t Virtuality not yet computed" << endl;
202 else cout << "\t Virtuality = " << virtuality << " GeV^2" << endl;
203 cout << endl;
204}
205
206float H_RecRPObject::getE() {
207 if(energy==NOT_YET_COMPUTED) {
208 cout<<"Please first compute energy using your favourite method"<<endl;
209 return NOT_YET_COMPUTED;
210 }
211 return energy;
212}
213
214float H_RecRPObject::getE(const unsigned int method) {
215 switch (method) {
216 case TM: {energy = computeE_TM();} break;
217 case AM: {energy = computeE_AM();} break;
218 case PM: {energy = computeE_PM();} break;
219 default: {energy = computeE_AM();} break;
220 };
221 return energy;
222}
223
224float H_RecRPObject::computeX0() {
225 if(energy==NOT_YET_COMPUTED) {
226 cout<<"Please first compute energy using your favourite method"<<endl;
227 return NOT_YET_COMPUTED;
228 }
229 float alpha1 = (matrp1->GetMatrixArray())[0];
230 float alpha2 = (matrp2->GetMatrixArray())[0];
231 float disp1 = (matrp1->GetMatrixArray())[4*MDIM]*URAD;
232 float disp2 = (matrp2->GetMatrixArray())[4*MDIM]*URAD;
233 x0 = (disp2*x1-disp1*x2)/(disp2*alpha1-disp1*alpha2);
234 return x0;
235}
236
237float H_RecRPObject::computeY0() {
238 if(energy==NOT_YET_COMPUTED) {
239 cout<<"Please first compute energy using your favourite method"<<endl;
240 return NOT_YET_COMPUTED;
241 }
242 float gamma1 = (matrp1->GetMatrixArray())[2*MDIM+2];
243 float gamma2 = (matrp2->GetMatrixArray())[2*MDIM+2];
244 float delta1 = (matrp1->GetMatrixArray())[3*MDIM+2];
245 float delta2 = (matrp2->GetMatrixArray())[3*MDIM+2];
246 y0 = (delta2*y1-delta1*y2)/(delta2*gamma1-delta1*gamma2);
247 return y0;
248}
249
250float H_RecRPObject::computeTX() {
251 if(energy==NOT_YET_COMPUTED) {
252 cout<<"Please first compute energy using your favourite method"<<endl;
253 return NOT_YET_COMPUTED;
254 }
255 float beta1 = (matrp1->GetMatrixArray())[1*MDIM];
256 float beta2 = (matrp2->GetMatrixArray())[1*MDIM];
257 float disp1 = (matrp1->GetMatrixArray())[4*MDIM]*URAD;
258 float disp2 = (matrp2->GetMatrixArray())[4*MDIM]*URAD;
259 // computes thx (murad)
260 thx = (x1*disp2-x2*disp1)/(beta1*disp2-beta2*disp1);
261 return thx;
262}
263
264float H_RecRPObject::computeTY() {
265 if(energy==NOT_YET_COMPUTED) {
266 cout<<"Please first compute energy using your favourite method"<<endl;
267 return NOT_YET_COMPUTED;
268 }
269 float gamma1 = (matrp1->GetMatrixArray())[2*MDIM+2];
270 float gamma2 = (matrp2->GetMatrixArray())[2*MDIM+2];
271 float delta1 = (matrp1->GetMatrixArray())[3*MDIM+2];
272 float delta2 = (matrp2->GetMatrixArray())[3*MDIM+2];
273 // computes thy (murad)
274 thy = (y1*gamma2-y2*gamma1)/(delta1*gamma2-delta2*gamma1);
275 return thy;
276}
277
278float H_RecRPObject::computeE_TM() {
279 // computes the emitted particle energy, from the trivial method
280 float disp = ((thebeam->getPartialMatrix("rp1",0,MP,QP))->GetMatrixArray())[4*MDIM]*URAD;
281 energy = -x1/disp;
282 // corrects for nonlinear effects
283 energy = (1+corr1_TM)*energy + corr2_TM*energy*energy;
284 // sets the rp matrices at obtained energy
285 delete matrp1;
286 delete matrp2;
287 matrp1 = new TMatrix(*(thebeam->getPartialMatrix("rp1",energy,MP,QP)));
288 matrp2 = new TMatrix(*(thebeam->getPartialMatrix("rp2",energy,MP,QP)));
289 // returns ...
290 return energy;
291}
292
293float H_RecRPObject::computeE_AM() {
294 // computes the emitted particle energy, from the angle compensation method, iterative way
295 const int N = 10;
296 delete matrp1;
297 delete matrp2;
298 matrp1 = new TMatrix(*(thebeam->getPartialMatrix("rp1",0,MP,QP)));
299 float disp = (matrp1->GetMatrixArray())[4*MDIM]*URAD;
300 delete matrp1;
301 energy = -x1/disp;
302 matrp1 = new TMatrix(*(thebeam->getPartialMatrix("rp1",energy,MP,QP)));
303 matrp2 = new TMatrix(*(thebeam->getPartialMatrix("rp2",energy,MP,QP)));
304 float beta1 = (matrp1->GetMatrixArray())[1*MDIM];
305 float beta2 = (matrp2->GetMatrixArray())[1*MDIM];
306 float disp1 = (matrp1->GetMatrixArray())[4*MDIM]*URAD;
307 float disp2 = (matrp2->GetMatrixArray())[4*MDIM]*URAD;
308 for(int i = 0; i < N; i++) {
309 energy = -(beta2*x1-beta1*x2)/(beta2*disp1-beta1*disp2);
310 delete matrp1;
311 delete matrp2;
312 matrp1 = new TMatrix(*(thebeam->getPartialMatrix("rp1",energy,MP,QP)));
313 matrp2 = new TMatrix(*(thebeam->getPartialMatrix("rp2",energy,MP,QP)));
314 beta1 = (matrp1->GetMatrixArray())[1*MDIM];
315 beta2 = (matrp2->GetMatrixArray())[1*MDIM];
316 disp1 = (matrp1->GetMatrixArray())[4*MDIM]*URAD;
317 disp2 = (matrp2->GetMatrixArray())[4*MDIM]*URAD;
318 }
319 // returns ...
320 return energy;
321}
322
323float H_RecRPObject::computeE_PM() {
324 cout<<"Not yet implemented, nothing done"<<endl;
325 energy = NOT_YET_COMPUTED;
326 return energy;
327}
328
329float H_RecRPObject::computeQ2() {
330 // computes the emitted particle virtuality
331 // energy should be teconstructed first
332 if(energy==NOT_YET_COMPUTED) {
333 cout<<"Please first compute energy using your favourite method"<<endl;
334 return NOT_YET_COMPUTED;
335 }
336 // getting parameters for reconstructed particle energy
337 float beta1 = (matrp1->GetMatrixArray())[1*MDIM];
338 float beta2 = (matrp2->GetMatrixArray())[1*MDIM];
339 float gamma1 = (matrp1->GetMatrixArray())[2*MDIM+2];
340 float gamma2 = (matrp2->GetMatrixArray())[2*MDIM+2];
341 float delta1 = (matrp1->GetMatrixArray())[3*MDIM+2];
342 float delta2 = (matrp2->GetMatrixArray())[3*MDIM+2];
343 float disp1 = (matrp1->GetMatrixArray())[4*MDIM]*URAD;
344 float disp2 = (matrp2->GetMatrixArray())[4*MDIM]*URAD;
345 // angles reconstruction
346 float rec_thx = (x1*disp2-x2*disp1)/(beta1*disp2-beta2*disp1)/URAD;
347 float rec_thy = (y1*gamma2-y2*gamma1)/(delta1*gamma2-delta2*gamma1)/URAD;
348 // Q² reconstruction
349 virtuality = BE*(BE-energy)*(rec_thx*rec_thx+rec_thy*rec_thy);
350 // returns ...
351 return virtuality;
352}
353
Note: See TracBrowser for help on using the repository browser.