/[PAMELA software]/DarthVader/TrackerLevel2/src/ExtTrack.cpp
ViewVC logotype

Contents of /DarthVader/TrackerLevel2/src/ExtTrack.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (show annotations) (download)
Wed Jun 4 07:57:03 2014 UTC (10 years, 6 months ago) by pam-ts
Branch: MAIN
Changes since 1.1: +778 -209 lines
New tracking algorythm implementation (extended to up to 2 calorimeter planes and with level1 cleaning for nuclei)

1 /**
2 * \file ExtTrack.cpp
3 * \author Elena Vannuccini
4 */
5 #include <ExtTrack.h>
6 #include <iostream>
7 #include <math.h>
8 using namespace std;
9
10 ExtTrack::ExtTrack(int dim){
11
12 //cout << " ExtTrack::ExtTrack("<<dim<<")"<<endl;
13 xgood = NULL;
14 ygood = NULL;
15 multmaxx = NULL;
16 multmaxy = NULL;
17 xm = NULL;
18 ym = NULL;
19 zm = NULL;
20 xma = NULL;
21 yma = NULL;
22 zma = NULL;
23 xmb = NULL;
24 ymb = NULL;
25 zmb = NULL;
26 resx = NULL;
27 resy = NULL;
28 xv = NULL;
29 yv = NULL;
30 zv = NULL;
31 axv = NULL;
32 ayv = NULL;
33 dedx_x = NULL;
34 dedx_y = NULL;
35
36 // xGF = new float[TrkParams::nGF];
37 // yGF = new float[TrkParams::nGF];
38 // init vectors
39 SetDimension(dim);//allocate vectors
40 // init fit variables
41 Reset();
42 // set fit parameters (default)
43 SetZ0(23.5);
44 SetMiniDefault();
45 };
46 //--------------------------------------
47 //
48 //
49 //--------------------------------------
50 ExtTrack::ExtTrack( const ExtTrack& t ){
51 //cout << " ExtTrack::ExtTrack( const ExtTrack& t ) "<<this<<endl;
52
53 // --------------------
54 // initialize pointers
55 // --------------------
56
57 xgood = NULL;
58 ygood = NULL;
59 multmaxx = NULL;
60 multmaxy = NULL;
61 xm = NULL;
62 ym = NULL;
63 zm = NULL;
64 xma = NULL;
65 yma = NULL;
66 zma = NULL;
67 xmb = NULL;
68 ymb = NULL;
69 zmb = NULL;
70 resx = NULL;
71 resy = NULL;
72 xv = NULL;
73 yv = NULL;
74 zv = NULL;
75 axv = NULL;
76 ayv = NULL;
77 dedx_x = NULL;
78 dedx_y = NULL;
79
80 // xGF = new float[TrkParams::nGF];
81 // yGF = new float[TrkParams::nGF];
82
83 SetDimension(t.nplanes);//allocate vectors
84
85 // --------------------
86 // copy values
87 // --------------------
88 zini = t.zini;
89 chi2 = t.chi2;
90 nstep = t.nstep;
91 for(int it1=0;it1<5;it1++){
92 al[it1] = t.al[it1];
93 for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
94 };
95
96
97 for(int ip=0;ip<nplanes;ip++){
98 xgood[ip] = t.xgood[ip];
99 ygood[ip] = t.ygood[ip];
100 multmaxx[ip] = t.multmaxx[ip];
101 multmaxy[ip] = t.multmaxy[ip];
102 xm[ip] = t.xm[ip];
103 ym[ip] = t.ym[ip];
104 zm[ip] = t.zm[ip];
105 xma[ip] = t.xma[ip];
106 yma[ip] = t.yma[ip];
107 zma[ip] = t.zma[ip];
108 xmb[ip] = t.xmb[ip];
109 ymb[ip] = t.ymb[ip];
110 zmb[ip] = t.zmb[ip];
111 resx[ip] = t.resx[ip];
112 resy[ip] = t.resy[ip];
113 xv[ip] = t.xv[ip];
114 yv[ip] = t.yv[ip];
115 zv[ip] = t.zv[ip];
116 axv[ip] = t.axv[ip];
117 ayv[ip] = t.ayv[ip];
118 dedx_x[ip] = t.dedx_x[ip];
119 dedx_y[ip] = t.dedx_y[ip];
120 };
121
122 int ngf = TrkParams::nGF;
123 for(int i=0; i<ngf; i++){
124 xGF[i] = t.xGF[i];
125 yGF[i] = t.yGF[i];
126 }
127
128 };
129 //--------------------------------------
130 //
131 //
132 //--------------------------------------
133 void ExtTrack::Copy( ExtTrack& t ){
134
135 // cout << " ExtTrack::Copy( ExtTrack& t ) "<<endl;
136
137 t.Reset();//reset variables
138 t.Clear();//deallocate vectors
139 t.SetDimension(nplanes);//allocate vectors
140
141 t.zini = zini;
142 t.chi2 = chi2;
143 t.nstep = nstep;
144 for(int it1=0;it1<5;it1++){
145 t.al[it1] = al[it1];
146 for(int it2=0;it2<5;it2++)t.coval[it1][it2] = coval[it1][it2];
147 };
148 // cout << "ExtTrack::Copy()........nplanes "<<nplanes<<endl;
149 for(int ip=0;ip<nplanes;ip++){
150 t.xgood[ip] = xgood[ip];
151 t.ygood[ip] = ygood[ip];
152 t.multmaxx[ip] = multmaxx[ip];
153 t.multmaxy[ip] = multmaxy[ip];
154 t.xm[ip] = xm[ip];
155 t.ym[ip] = ym[ip];
156 t.zm[ip] = zm[ip];
157 t.xma[ip] = xma[ip];
158 t.yma[ip] = yma[ip];
159 t.zma[ip] = zma[ip];
160 t.xmb[ip] = xmb[ip];
161 t.ymb[ip] = ymb[ip];
162 t.zmb[ip] = zmb[ip];
163 t.resx[ip] = resx[ip];
164 t.resy[ip] = resy[ip];
165 t.xv[ip] = xv[ip];
166 t.yv[ip] = yv[ip];
167 t.zv[ip] = zv[ip];
168 t.axv[ip] = axv[ip];
169 t.ayv[ip] = ayv[ip];
170 t.dedx_x[ip] = dedx_x[ip];
171 t.dedx_y[ip] = dedx_y[ip];
172 };
173 int ngf = TrkParams::nGF;
174 for(int i=0; i<ngf; i++){
175 t.xGF[i] = xGF[i];
176 t.yGF[i] = yGF[i];
177 }
178
179 };
180 //--------------------------------------
181 //
182 //
183 //--------------------------------------
184 void ExtTrack::Set( TrkTrack& t , int index){
185
186 // cout << "Set("<< &t <<", "<<index<<")"<<endl;
187 if(index < 0 ){ //NON FUNZIONA
188 cout << "Clear() "<<endl;
189 Reset();//reset variables
190 Clear();//deallocate vectors
191 cout << "SetDimension("<<t.GetNhit()<<")"<<endl;
192 SetDimension(t.GetNhit()); //allocate vectors
193 }
194 // default fit settings (the same as TrkTrack)
195 SetMiniDefault();
196 SetZ0(23.5);
197 //
198 int ipp;
199 if(index<0) ipp = 0;
200 else ipp = index;
201 //
202 // cout << "Assign... nplanes = "<<nplanes<<endl;
203 for(int ip=0;ip<6;ip++){
204
205 if(ipp>=nplanes){
206 cout << "void ExtTrack::Set( TrkTrack& , "<<index<<") -- ipp="<<ipp<<" exceed vector dimention "<<nplanes<<endl;
207 continue;
208 }
209
210 if( t.XGood(ip) || t.YGood(ip) || !(index<0) ){
211
212 // cout << "ipp "<<ipp<<" ip "<<ip<<endl;
213
214
215 xgood[ipp] = t.XGood(ip); //NB!
216 ygood[ipp] = t.YGood(ip); //NB!
217 multmaxx[ipp] = t.multmaxx[ip];
218 multmaxy[ipp] = t.multmaxy[ip];
219 xm[ipp] = t.xm[ip];
220 ym[ipp] = t.ym[ip];
221 zm[ipp] = t.zm[ip];
222
223 float dummy = -100.;
224 xma[ipp] = dummy;
225 yma[ipp] = dummy;
226 zma[ipp] = zm[ipp];
227 xmb[ipp] = dummy;
228 ymb[ipp] = dummy;
229 zmb[ipp] = zm[ipp];
230
231 if( t.XGood(ip)*t.YGood(ip) ){ //double hit
232 xma[ipp] = t.xm[ip];
233 yma[ipp] = t.ym[ip];
234 zma[ipp] = t.zm[ip];
235 xmb[ipp] = t.xm[ip];
236 ymb[ipp] = t.ym[ip];
237 zmb[ipp] = t.zm[ip];
238 } else if (t.XGood(ip) || t.YGood(ip) ){ //single hit
239 TrkParams::Load(5);
240 if( !TrkParams::IsLoaded(5) ){
241 cout << "void ExtTrack::FillMiniStruct(TrkTrack&) ---ERROR--- trk align.param. not loaded "<<endl;
242 return;
243 }
244 double segment = 7.;// 2.;//cm //Elena 10th
245 int is = (int)t.GetSensor(ip); if(ip==5)is=1-is;
246 int ippp = 5-ip;
247 int il = (int)t.GetLadder(ip);
248 // NB: i parametri di allineamento hanno una notazione particolare!!!
249 // sensor = 0 (hybrid side), 1
250 // ladder = 0-2 (increasing x)
251 // plane = 0-5 (from bottom to top!!!)
252 double omega = 0.;
253 double beta = 0.;
254 double gamma = 0.;
255 if(
256 (is < 0 || is > 1 || ip < 0 || ippp > 5 || il < 0 || il > 2) &&
257 true){
258 cout << " void ExtTrack::FillMiniStruct(TrkTrack&) --- WARNING ---trk sensor not defined, cannot read alignment parameters "<<endl;
259 cout << " is ip il = "<<is<<" "<<ippp<<" "<<il<<endl;
260 }else{
261 omega = alignparameters_.omega[is][il][ippp];
262 beta = alignparameters_.beta[is][il][ippp];
263 gamma = alignparameters_.gamma[is][il][ippp];
264 }
265 if( t.XGood(ip) && !t.YGood(ip) ){
266 xma[ipp] = t.xm[ip] - omega * segment;
267 yma[ipp] = t.ym[ip] + segment;
268 zma[ipp] = t.zm[ip] + beta * segment;//not used yet
269 xmb[ipp] = t.xm[ip] + omega * segment;
270 ymb[ipp] = t.ym[ip] - segment;
271 zmb[ipp] = t.zm[ip] - beta * segment;//not used yet
272 }else if( !t.XGood(ip) && t.YGood(ip) ){
273 xma[ipp] = t.xm[ip] + segment;
274 yma[ipp] = t.ym[ip] + omega * segment;
275 zma[ipp] = t.zm[ip] - gamma * segment;//not used yet
276 xmb[ipp] = t.xm[ip] - segment;
277 ymb[ipp] = t.ym[ip] - omega * segment;
278 zmb[ipp] = t.zm[ip] + gamma * segment;//not used yet
279 }
280 }
281
282 // tailx[ipp] = t.tailx[ip];
283 // taily[ipp] = t.taily[ip];
284 resx[ipp] = t.resx[ip];
285 resy[ipp] = t.resy[ip];
286 xv[ipp] = t.xv[ip];
287 yv[ipp] = t.yv[ip];
288 zv[ipp] = t.zv[ip];
289 axv[ipp] = t.axv[ip];
290 ayv[ipp] = t.ayv[ip];
291 dedx_x[ipp] = t.dedx_x[ip];
292 dedx_y[ipp] = t.dedx_y[ip];
293
294 ipp++;
295
296 }
297 };//endl loop on planes
298
299 chi2 = t.chi2;
300 nstep = t.nstep;
301 for(int it1=0;it1<5;it1++){
302 al[it1] = t.al[it1];
303 for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
304 };
305
306 };
307 //--------------------------------------
308 //
309 //
310 //--------------------------------------
311 void ExtTrack::SetMiniDefault(){
312 // cout << " ExtTrack::SetMiniDefault()"<<endl;
313 exttrack_.trackmode = TrkParams::init__mini_trackmode;
314 exttrack_.fact = TrkParams::init__mini_fact;
315 exttrack_.istepmin = TrkParams::init__mini_istepmin;
316 TrkParams::SetDeltaB();
317 TrkParams::SetDLT();
318 }
319 //--------------------------------------
320 //
321 //
322 //--------------------------------------
323 void ExtTrack::SetDimension(int dim){
324
325 //cout << " ExtTrack::SetDimension("<<dim<<")"<<endl;;
326 if(dim<=0){
327 nplanes = 0;
328 if(dim<0)cout << "void ExtTrack::SetDimension("<<dim<<") --- invalid dimension "<<endl;
329 return;
330 }
331
332 // Clear();
333
334 nplanes = dim;
335
336 xgood = new int[nplanes];
337 ygood = new int[nplanes];
338 multmaxx = new int[nplanes];
339 multmaxy = new int[nplanes];
340 xm = new float[nplanes];
341 ym = new float[nplanes];
342 zm = new float[nplanes];
343 xma = new float[nplanes];
344 yma = new float[nplanes];
345 zma = new float[nplanes];
346 xmb = new float[nplanes];
347 ymb = new float[nplanes];
348 zmb = new float[nplanes];
349 resx = new float[nplanes];
350 resy = new float[nplanes];
351 xv = new float[nplanes];
352 yv = new float[nplanes];
353 zv = new float[nplanes];
354 axv = new float[nplanes];
355 ayv = new float[nplanes];
356 dedx_x = new float[nplanes];
357 dedx_y = new float[nplanes];
358
359
360 }
361 //--------------------------------------
362 //
363 //
364 //--------------------------------------
365 bool ExtTrack::SetZ(int ip,float zmeas){
366 if(ip<0 || ip>=nplanes){
367 cout << " ExtTrack::SetZ("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
368 return false;
369 }
370 zm[ip] = zmeas;
371 zma[ip] = zmeas;
372 zmb[ip] = zmeas;
373 return true;
374
375 };
376
377 bool ExtTrack::ResetXY(int ip){
378
379 // cout << " ExtTrack::ResetXY("<<ip<<") "<<endl;
380
381 if(ip<0 || ip>=nplanes){
382 cout << " ExtTrack::ResetXY("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
383 return false;
384 }
385 xm[ip] = -100.;
386 ym[ip] = -100.;
387 xma[ip] = -100.;
388 yma[ip] = -100.;
389 xmb[ip] = -100.;
390 ymb[ip] = -100.;
391 resx[ip] = 1000.;
392 resy[ip] = 1000.;
393 dedx_x[ip] = 0.;
394 dedx_y[ip] = 0.;
395 xgood[ip] = 0;
396 ygood[ip] = 0;
397 multmaxx[ip] = 0;
398 multmaxy[ip] = 0;
399
400 return true;
401
402 };
403
404 bool ExtTrack::SetXY(int ip,float xmeas, float ymeas, float rx, float ry){
405
406 if(ip<0 || ip>=nplanes){
407 cout << " ExtTrack::SetXY("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
408 return false;
409 }
410
411 xm[ip] = xmeas;
412 ym[ip] = ymeas;
413 xma[ip] = xmeas;
414 yma[ip] = ymeas;
415 xmb[ip] = xmeas;
416 ymb[ip] = ymeas;
417 resx[ip] = rx;
418 resy[ip] = ry;
419 // dedx_x[ip] = dedxx;
420 // dedx_y[ip] = dedxy;
421 xgood[ip] = 1;
422 ygood[ip] = 1;
423
424 return true;
425
426 };
427 bool ExtTrack::SetX(int ip,float xa, float xb, float ya, float yb,float res ){
428
429 // cout << " ExtTrack::SetX("<<ip<<",...) -- nplanes"<<nplanes<<endl;
430
431 if(ip<0 || ip>=nplanes){
432 cout << " ExtTrack::SetX("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
433 return false;
434 }
435 xm[ip] = -100.;
436 ym[ip] = -100.;
437 xma[ip] = xa;
438 yma[ip] = ya;
439 xmb[ip] = xb;
440 ymb[ip] = yb;
441 resx[ip] = res;
442 resy[ip] = 1000.;
443 xgood[ip] = 1;
444 ygood[ip] = 0;
445 // dedx_x[ip] = dedx;
446 // dedx_y[ip] = 0.;
447 return true;
448 };
449 bool ExtTrack::SetY(int ip,float xa, float xb, float ya, float yb,float res ){
450
451 // cout << " ExtTrack::SetY("<<ip<<",...) -- nplanes"<<nplanes<<endl;
452
453 if(ip<0 || ip>=nplanes){
454 cout << " ExtTrack::SetX("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
455 return false;
456 }
457 xm[ip] = -100.;
458 ym[ip] = -100.;
459 xma[ip] = xa;
460 yma[ip] = ya;
461 xmb[ip] = xb;
462 ymb[ip] = yb;
463 resx[ip] = 1000.;
464 resy[ip] = res;
465 xgood[ip] = 0;
466 ygood[ip] = 1;
467 // dedx_x[ip] = 0.;
468 // dedx_y[ip] = dedx;
469 return true;
470 };
471
472 bool ExtTrack::SetXGood(int ip,int icl_piu_uno, int il, int is ){
473
474 // cout << " ExtTrack::SetXGood("<<ip<<",...) -- nplanes"<<nplanes<<endl;
475
476 if(ip<0 || ip>=nplanes){
477 cout << " ExtTrack::SetXGood("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
478 return false;
479 }
480 xgood[ip]=(il+1)*100000000+(is+1)*10000000+icl_piu_uno;
481 return true;
482 };
483 bool ExtTrack::SetYGood(int ip,int icl_piu_uno, int il, int is ){
484 if(ip<0 || ip>=nplanes){
485 cout << " ExtTrack::SetYGood("<<ip<<",...) -- nplanes"<<nplanes<<endl;
486 cout << " ExtTrack::SetYGood("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
487 return false;
488 }
489 ygood[ip]=(il+1)*100000000+(is+1)*10000000+icl_piu_uno;
490 return true;
491 };
492
493 //--------------------------------------
494 //
495 //
496 //--------------------------------------
497 void ExtTrack::Clear(Option_t* option){
498
499 // cout << " ExtTrack::Clear("<<option<<") "<<this<<endl;
500
501
502 // cout << " xgood "<<xgood<<endl;
503
504 // if(nplanes>0){
505 if(xgood) delete [] xgood;
506 if(ygood) delete [] ygood;
507 if(multmaxx) delete [] multmaxx;
508 if(multmaxy) delete [] multmaxy;
509 if(xm) delete [] xm;
510 if(ym) delete [] ym;
511 if(zm) delete [] zm;
512 if(xma) delete [] xma;
513 if(yma) delete [] yma;
514 if(zma) delete [] zma;
515 if(xmb) delete [] xmb;
516 if(ymb) delete [] ymb;
517 if(zmb) delete [] zmb;
518 if(resx) delete [] resx;
519 if(resy) delete [] resy;
520 if(xv) delete [] xv;
521 if(yv) delete [] yv;
522 if(zv) delete [] zv;
523 if(axv) delete [] axv;
524 if(ayv) delete [] ayv;
525 if(dedx_x)delete [] dedx_x;
526 if(dedx_y)delete [] dedx_y;
527 //}
528
529 xgood = NULL;
530 ygood = NULL;
531 multmaxx = NULL;
532 multmaxy = NULL;
533 xm = NULL;
534 ym = NULL;
535 zm = NULL;
536 xma = NULL;
537 yma = NULL;
538 zma = NULL;
539 xmb = NULL;
540 ymb = NULL;
541 zmb = NULL;
542 resx = NULL;
543 resy = NULL;
544 xv = NULL;
545 yv = NULL;
546 zv = NULL;
547 axv = NULL;
548 ayv = NULL;
549 dedx_x = NULL;
550 dedx_y = NULL;
551
552 nplanes = 0;
553
554 // Reset();
555
556 };
557
558 void ExtTrack::Delete(){
559 Clear();
560 // delete [] xGF;
561 // delete [] yGF;
562 }
563
564 //--------------------------------------
565 //
566 //
567 //--------------------------------------
568
569 // /**
570 // * Set the position measurements
571 // */
572 // void ExtTrack::SetMeasure(Double_t *xmeas, Double_t *ymeas, Double_t *zmeas){
573 // for(int i=0; i<nplanes; i++) xm[i]= (xmeas? *xmeas++ : 0.);
574 // for(int i=0; i<nplanes; i++) ym[i]= (ymeas? *ymeas++ : 0.);
575 // for(int i=0; i<nplanes; i++) zm[i]= (ymeas? *zmeas++ : 0.);
576 // }
577 // /**
578 // * Set the TrkTrack position resolution
579 // */
580 // void ExtTrack::SetResolution(Double_t *rx, Double_t *ry){
581 // for(int i=0; i<nplanes; i++) resx[i]= (rx? *rx++ : 0.);
582 // for(int i=0; i<nplanes; i++) resy[i]= (ry? *ry++ : 0.);
583 // }
584 // /**
585 // * Set the TrkTrack good measurement
586 // */
587 // void ExtTrack::SetGood(int *xg, int *yg){
588
589 // for(int i=0; i<nplanes; i++) xgood[i]=(xg? *xg++ : 0 );
590 // for(int i=0; i<nplanes; i++) ygood[i]=(yg? *yg++ : 0 );
591 // }
592
593 //--------------------------------------
594 //
595 //
596 //--------------------------------------
597 void ExtTrack::Dump(){
598 cout << endl << "========== Track " ;
599 cout << endl << "zini : "<< zini;
600 cout << endl << "al : "; for(int i=0; i<5; i++)cout << al[i] << " ";
601 cout << endl << "chi^2 : "<< chi2;
602 cout << endl << "n.step : "<< nstep;
603 cout << endl << "xgood : "; for(int i=0; i<nplanes; i++)cout << XGood(i) ;
604 cout << endl << "ygood : "; for(int i=0; i<nplanes; i++)cout << YGood(i) ;
605 cout << endl << "xm : "; for(int i=0; i<nplanes; i++)cout << xm[i] << " ";
606 cout << endl << "ym : "; for(int i=0; i<nplanes; i++)cout << ym[i] << " ";
607 cout << endl << "zm : "; for(int i=0; i<nplanes; i++)cout << zm[i] << " ";
608 cout << endl << "xma : "; for(int i=0; i<nplanes; i++)cout << xma[i] << " ";
609 cout << endl << "yma : "; for(int i=0; i<nplanes; i++)cout << yma[i] << " ";
610 cout << endl << "zma : "; for(int i=0; i<nplanes; i++)cout << zma[i] << " ";
611 cout << endl << "xmb : "; for(int i=0; i<nplanes; i++)cout << xmb[i] << " ";
612 cout << endl << "ymb : "; for(int i=0; i<nplanes; i++)cout << ymb[i] << " ";
613 cout << endl << "zmb : "; for(int i=0; i<nplanes; i++)cout << zmb[i] << " ";
614 cout << endl << "xv : "; for(int i=0; i<nplanes; i++)cout << xv[i] << " ";
615 cout << endl << "yv : "; for(int i=0; i<nplanes; i++)cout << yv[i] << " ";
616 cout << endl << "zv : "; for(int i=0; i<nplanes; i++)cout << zv[i] << " ";
617 cout << endl << "resx : "; for(int i=0; i<nplanes; i++)cout << resx[i] << " ";
618 cout << endl << "resy : "; for(int i=0; i<nplanes; i++)cout << resy[i] << " ";
619 cout << endl << "coval : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" ";
620 cout << endl << " "; for(int i=0; i<5; i++)cout << coval[1][i]<<" ";
621 cout << endl << " "; for(int i=0; i<5; i++)cout << coval[2][i]<<" ";
622 cout << endl << " "; for(int i=0; i<5; i++)cout << coval[3][i]<<" ";
623 cout << endl << " "; for(int i=0; i<5; i++)cout << coval[4][i]<<" ";
624 cout << endl << "MDR : "<< (coval[4][4]>0. ? 1./sqrt(coval[4][4]) : 0. ) ;
625 cout << endl << "dedx_x : "; for(int i=0; i<nplanes; i++)cout << dedx_x[i] << " ";
626 cout << endl << "dedx_y : "; for(int i=0; i<nplanes; i++)cout << dedx_y[i] << " ";
627
628 cout << endl << "=================";
629 cout << endl;
630 }
631 /**
632 * Reset the fit parameters
633 */
634 void ExtTrack::ResetFit(){
635
636 // cout << " ExtTrack::ResetFit() "<<endl;
637 for(int i=0; i<5; i++) al[i]=-9999.;
638 chi2=0.;
639 nstep=0;
640
641 for(int i=0; i<5; i++) {
642 for(int j=0; j<5; j++) coval[i][j]=0.;
643 }
644
645 for(int i=0; i<nplanes; i++){
646 xv[i] = 0.;
647 yv[i] = 0.;
648 zv[i] = 0.;
649 axv[i] = 0.;
650 ayv[i] = 0.;
651 }
652
653 int ngf = TrkParams::nGF;
654 for(int i=0; i<ngf; i++){
655 xGF[i] = 0.;
656 yGF[i] = 0.;
657 }
658
659
660 }
661 /**
662 * Reset x/y measured coordinates (NB it does not reset z coordinates)
663 */
664 void ExtTrack::ResetXY(){
665
666 for(int i=0; i<nplanes; i++)ResetXY(i);
667
668
669 }
670 /**
671 * \brief Tracking method. It calls F77 mini routine.
672 *
673 * @see void TrkTrack::Fit(double, int&, int, int)
674 */
675 void ExtTrack::Fit(double pfixed, int& fail, int iprint){
676
677 TrkParams::Load(1);
678 if( !TrkParams::IsLoaded(1) ){
679 cout << "void ExtTrack::Fit(double,int&,int,int) ---ERROR--- m.field not loaded "<<endl;
680 return;
681 }
682
683 float al_ini[] = {0.,0.,0.,0.,0.};
684
685 fail = 0;
686
687 FillMiniStruct(); //fill F77 common
688
689
690 // if fit variables have been reset, evaluate the initial guess
691 if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.){
692 if(iprint==2)cout << "calling guessext_()"<<endl;
693 guessext_();
694 }
695 // --------------------- free momentum
696 if(pfixed==0.) {
697 exttrack_.pfixed=0.;
698 }
699 // --------------------- fixed momentum
700 if(pfixed!=0.) {
701 al[4]=1./pfixed;
702 exttrack_.pfixed=pfixed;
703 }
704
705 // store temporarily the initial guess
706 for(int i=0; i<5; i++) al_ini[i]=exttrack_.al[i];
707
708 if(iprint==2){
709 cout<<endl<<" init. guess: ";
710 for(int i=0; i<5; i++)cout << al_ini[i]<<" ";
711 cout<<endl;
712 }
713 // ------------------------------------------
714 // call mini routine
715 // ------------------------------------------
716 int istep=0;
717 int ifail=0;
718 if(iprint==2)cout << "calling miniext_("<<istep<<","<<ifail<<", "<<iprint<<")"<<endl;
719 miniext_(&istep,&ifail, &iprint);
720 if(ifail!=0) {
721 if(iprint)cout << "ERROR: ifail= " << ifail << endl;
722 fail = 1;
723 }
724 if(chi2!=chi2){
725 if(iprint)cout << "ERROR: chi2= " << chi2 << endl;
726 ResetFit();
727 fail = 1;
728 }
729 // ------------------------------------------
730
731 SetFromMiniStruct(); //copy content of F77 common into the ExtTrack object
732
733 if(fail){
734 if(iprint)cout << " >>>> fit failed "<<endl;
735 for(int i=0; i<5; i++) al[i]=al_ini[i];
736 }
737 };
738
739 /**
740 * Method to fill minimization-routine common
741 * NB (zini=23.5 is implicitelly set)
742 */
743 void ExtTrack::FillMiniStruct(cMiniExtTrack& track){
744
745 track.nplanes = nplanes;
746 for(int i=0; i<nplanes; i++){
747
748 track.xgood[i] = (double) XGood(i);
749 track.ygood[i] = (double) YGood(i);
750
751 track.xm[i] = (double) xm[i];
752 track.ym[i] = (double) ym[i];
753 track.zm[i] = (double) zm[i];
754
755
756 track.xm_a[i] = (double) xma[i];
757 track.xm_b[i] = (double) xmb[i];
758 track.ym_a[i] = (double) yma[i];
759 track.ym_b[i] = (double) ymb[i];
760 track.zm_a[i] = (double) zma[i];
761 track.zm_b[i] = (double) zmb[i];
762
763 track.resx[i] = (double) resx[i];
764 track.resy[i] = (double) resy[i];
765 }
766
767 for(int i=0; i<5; i++) track.al[i] = (double) al[i];
768
769 track.zini = (double) zini;
770
771
772 }
773 /**
774 * Method to set values from minimization-routine common
775 */
776 void ExtTrack::SetFromMiniStruct(cMiniExtTrack *track){
777
778 for(int i=0; i<5; i++) {
779 al[i] = (float) (track->al[i]);
780 for(int j=0; j<5; j++) coval[i][j] = (float) (track->cov[i][j]);
781 }
782 chi2 = (float) (track->chi2);
783 nstep = (float) (track->nstep);
784 for(int i=0; i<nplanes; i++){
785 xv[i] = (float) (track->xv[i]);
786 yv[i] = (float) (track->yv[i]);
787 zv[i] = (float) (track->zv[i]);
788 xm[i] = (float) (track->xm[i]);
789 ym[i] = (float) (track->ym[i]);
790 zm[i] = (float) (track->zm[i]);
791 xma[i] = (float) (track->xm_a[i]);
792 yma[i] = (float) (track->ym_a[i]);
793 zma[i] = (float) (track->zm_a[i]);
794 xmb[i] = (float) (track->xm_b[i]);
795 ymb[i] = (float) (track->ym_b[i]);
796 zmb[i] = (float) (track->zm_b[i]);
797 axv[i] = (float) (track->axv[i]);
798 ayv[i] = (float) (track->ayv[i]);
799 }
800 zini = (float) (track->zini);
801 }
802 /**
803 * Method to retrieve ID (0,1,...) of x-cluster (if any) associated to this track.
804 * If no cluster is associated, ID=-1.
805 *
806 */
807 int ExtTrack::GetClusterX_ID(int ip){
808 return ((int)fabs(xgood[ip]))%10000000-1;
809 };
810 /**
811 * Method to retrieve ID (0-xxx) of y-cluster (if any) associated to this track.
812 *
813 */
814 int ExtTrack::GetClusterY_ID(int ip){
815 return ((int)fabs(ygood[ip]))%10000000-1;
816 };
817
818 /**
819 * Method to retrieve the ladder, if assigned
820 */
821 int ExtTrack::GetLadder(int ip){
822 if(XGood(ip))return (int)fabs(xgood[ip]/100000000)-1;
823 if(YGood(ip))return (int)fabs(ygood[ip]/100000000)-1;
824 return -1;
825 };
826 /**
827 * Method to retrieve the sensor, if assigned
828 */
829 int ExtTrack::GetSensor(int ip){
830 if(XGood(ip))return (int)((int)fabs(xgood[ip]/10000000)%10)-1;
831 if(YGood(ip))return (int)((int)fabs(ygood[ip]/10000000)%10)-1;
832 return -1;
833 };
834 /**
835 *
836 */
837 Int_t ExtTrack::GetClusterX_Multiplicity(int ip){
838 if(ip<0 || ip>=nplanes){
839 cout << " ExtTrack::GetClusterX_Multiplicity("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
840 return -1;
841 }
842 return (Int_t)(multmaxx[ip]/10000);
843 }
844 /**
845 *
846 */
847
848 Int_t ExtTrack::GetClusterY_Multiplicity(int ip){
849 if(ip<0 || ip>=nplanes){
850 cout << " ExtTrack::GetClusterY_Multiplicity("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
851 return -1;
852 }
853 return (Int_t)(multmaxy[ip]/10000);
854 }
855 /**
856 *
857 */
858
859 Int_t ExtTrack::GetClusterX_MaxStrip(int ip){
860 if(ip<0 || ip>=nplanes){
861 cout << " ExtTrack::GetClusterX_MaxStrip("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
862 return -1;
863 }
864 return (Int_t)(multmaxx[ip]%10000);
865 }
866 /**
867 *
868 */
869
870 Int_t ExtTrack::GetClusterY_MaxStrip(int ip){
871 if(ip<0 || ip>=nplanes){
872 cout << " ExtTrack::GetClusterY_MaxStrip("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
873 return -1;
874 }
875 return (Int_t)(multmaxy[ip]%10000);
876 }
877 Float_t ExtTrack::GetRigidity(){
878 Float_t rig=0;
879 if(chi2>0)rig=1./al[4];
880 if(rig<0) rig=-rig;
881 return rig;
882 };
883 ClassImp(ExtTrack);

  ViewVC Help
Powered by ViewVC 1.1.23