[3] | 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"
|
---|
| 22 | using 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 |
|
---|
| 30 | H_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 |
|
---|
| 47 | H_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 |
|
---|
| 65 | void 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 |
|
---|
| 166 | void 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 |
|
---|
| 174 | void 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 |
|
---|
| 183 | void H_RecRPObject::setPosition_det1(const float xx1, const float yy1) {
|
---|
| 184 | energy = NOT_YET_COMPUTED;
|
---|
| 185 | x1 = xx1;
|
---|
| 186 | y1 = yy1;
|
---|
| 187 | }
|
---|
| 188 |
|
---|
| 189 | void H_RecRPObject::setPosition_det2(const float xx2, const float yy2) {
|
---|
| 190 | energy = NOT_YET_COMPUTED;
|
---|
| 191 | x2 = xx2;
|
---|
| 192 | y2 = yy2;
|
---|
| 193 | }
|
---|
| 194 |
|
---|
| 195 | void 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 |
|
---|
| 213 | void 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 |
|
---|
| 239 | void 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 |
|
---|
| 279 | float 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 |
|
---|
| 285 | float H_RecRPObject::getE() {
|
---|
| 286 | if(energy==NOT_YET_COMPUTED) { computeAll(); };
|
---|
| 287 | return energy;
|
---|
| 288 | }
|
---|
| 289 |
|
---|
| 290 | float H_RecRPObject::getTX() {
|
---|
| 291 | if(energy==NOT_YET_COMPUTED) { computeAll(); };
|
---|
| 292 | return txip;
|
---|
| 293 | }
|
---|
| 294 |
|
---|
| 295 | float H_RecRPObject::getTY() {
|
---|
| 296 | if(energy==NOT_YET_COMPUTED) { computeAll(); };
|
---|
| 297 | return tyip;
|
---|
| 298 | }
|
---|
| 299 |
|
---|
| 300 | float H_RecRPObject::getQ2() {
|
---|
| 301 | if(energy==NOT_YET_COMPUTED) { computeAll(); };
|
---|
| 302 | cout<<"<H_RecRPObject::getQ2> Not implemented yet"<<endl;
|
---|
| 303 | return 0;
|
---|
| 304 | }
|
---|
| 305 |
|
---|
| 306 | float H_RecRPObject::getPt() {
|
---|
| 307 | if(energy==NOT_YET_COMPUTED) { computeAll(); };
|
---|
| 308 | return pt;
|
---|
| 309 | }
|
---|