/[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.9 - (show annotations) (download)
Fri Feb 26 12:03:36 2016 UTC (8 years, 9 months ago) by pam-fi
Branch: MAIN
CVS Tags: HEAD
Changes since 1.8: +8 -0 lines
patch to assign (locally) the xGF and yGF coordinates

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("C");//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(\"C\") "<<endl;
189 Reset();//reset variables
190 Clear("C");//deallocate vectors
191 cout << "SetDimension("<<t.GetNhit()<<")"<<endl;
192 SetDimension(t.GetNhit()); //allocate vectors
193 // }else{
194 // SetDimension(6+index);
195 }
196 // default fit settings (the same as TrkTrack)
197 SetMiniDefault();
198 SetZ0(23.5);
199 //
200 int ipp;
201 if(index<0) ipp = 0;
202 else ipp = index;
203 //
204 // cout << "Assign... nplanes = "<<nplanes<<endl;
205 for(int ip=0;ip<6;ip++){
206
207 if(ipp>=nplanes){
208 cout << "void ExtTrack::Set( TrkTrack& , "<<index<<") -- ipp="<<ipp<<" exceed vector dimention "<<nplanes<<endl;
209 continue;
210 }
211
212 if( t.XGood(ip) || t.YGood(ip) || !(index<0) ){
213
214 // cout << "ipp "<<ipp<<" ip "<<ip<<endl;
215
216
217 xgood[ipp] = t.XGood(ip); //NB!
218 ygood[ipp] = t.YGood(ip); //NB!
219 multmaxx[ipp] = t.multmaxx[ip];
220 multmaxy[ipp] = t.multmaxy[ip];
221 xm[ipp] = t.xm[ip];
222 ym[ipp] = t.ym[ip];
223 zm[ipp] = t.zm[ip];
224
225 float dummy = -100.;
226 xma[ipp] = dummy;
227 yma[ipp] = dummy;
228 zma[ipp] = zm[ipp];
229 xmb[ipp] = dummy;
230 ymb[ipp] = dummy;
231 zmb[ipp] = zm[ipp];
232
233 if( t.XGood(ip) && t.YGood(ip) ){ //double hit
234 xma[ipp] = t.xm[ip];
235 yma[ipp] = t.ym[ip];
236 zma[ipp] = t.zm[ip];
237 xmb[ipp] = t.xm[ip];
238 ymb[ipp] = t.ym[ip];
239 zmb[ipp] = t.zm[ip];
240 } else if (t.XGood(ip) || t.YGood(ip) ){ //single hit
241 TrkParams::Load(5);
242 if( !TrkParams::IsLoaded(5) ){
243 cout << "void ExtTrack::FillMiniStruct(TrkTrack&) ---ERROR--- trk align.param. not loaded "<<endl;
244 return;
245 }
246 double segment = 7.;// 2.;//cm //Elena 10th
247 int is = (int)t.GetSensor(ip); if(ip==5)is=1-is;
248 int ippp = 5-ip;
249 int il = (int)t.GetLadder(ip);
250 // NB: i parametri di allineamento hanno una notazione particolare!!!
251 // sensor = 0 (hybrid side), 1
252 // ladder = 0-2 (increasing x)
253 // plane = 0-5 (from bottom to top!!!)
254 double omega = 0.;
255 double beta = 0.;
256 double gamma = 0.;
257 if(
258 (is < 0 || is > 1 || ip < 0 || ippp > 5 || il < 0 || il > 2) &&
259 true){
260 cout << " void ExtTrack::FillMiniStruct(TrkTrack&) --- WARNING ---trk sensor not defined, cannot read alignment parameters "<<endl;
261 cout << " is ip il = "<<is<<" "<<ippp<<" "<<il<<endl;
262 }else{
263 omega = alignparameters_.omega[is][il][ippp];
264 beta = alignparameters_.beta[is][il][ippp];
265 gamma = alignparameters_.gamma[is][il][ippp];
266 }
267 if( t.XGood(ip) && !t.YGood(ip) ){
268 xma[ipp] = t.xm[ip] - omega * segment;
269 yma[ipp] = t.ym[ip] + segment;
270 zma[ipp] = t.zm[ip] + beta * segment;//not used yet
271 xmb[ipp] = t.xm[ip] + omega * segment;
272 ymb[ipp] = t.ym[ip] - segment;
273 zmb[ipp] = t.zm[ip] - beta * segment;//not used yet
274 }else if( !t.XGood(ip) && t.YGood(ip) ){
275 xma[ipp] = t.xm[ip] + segment;
276 yma[ipp] = t.ym[ip] + omega * segment;
277 zma[ipp] = t.zm[ip] - gamma * segment;//not used yet
278 xmb[ipp] = t.xm[ip] - segment;
279 ymb[ipp] = t.ym[ip] - omega * segment;
280 zmb[ipp] = t.zm[ip] + gamma * segment;//not used yet
281 }
282 }
283
284 // tailx[ipp] = t.tailx[ip];
285 // taily[ipp] = t.taily[ip];
286 resx[ipp] = t.resx[ip];
287 resy[ipp] = t.resy[ip];
288 xv[ipp] = t.xv[ip];
289 yv[ipp] = t.yv[ip];
290 zv[ipp] = t.zv[ip];
291 axv[ipp] = t.axv[ip];
292 ayv[ipp] = t.ayv[ip];
293 dedx_x[ipp] = t.dedx_x[ip];
294 dedx_y[ipp] = t.dedx_y[ip];
295
296 ipp++;
297
298 }
299 };//endl loop on planes
300
301 chi2 = t.chi2;
302 nstep = t.nstep;
303 for(int it1=0;it1<5;it1++){
304 al[it1] = t.al[it1];
305 for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
306 };
307
308 };
309 //--------------------------------------
310 //
311 //
312 //--------------------------------------
313 void ExtTrack::SetMiniDefault(){
314 // cout << " ExtTrack::SetMiniDefault()"<<endl;
315 exttrack_.trackmode = TrkParams::init__mini_trackmode;
316 exttrack_.fact = TrkParams::init__mini_fact;
317 exttrack_.istepmin = TrkParams::init__mini_istepmin;
318 TrkParams::SetDeltaB();
319 TrkParams::SetDLT();
320 }
321 //--------------------------------------
322 //
323 //
324 //--------------------------------------
325 void ExtTrack::SetDimension(int dim){
326
327 //cout << " ExtTrack::SetDimension("<<dim<<")"<<endl;;
328 if(dim<=0){
329 nplanes = 0;
330 if(dim<0)cout << "void ExtTrack::SetDimension("<<dim<<") --- invalid dimension "<<endl;
331 return;
332 }
333
334 // Clear();
335
336 nplanes = dim;
337
338 xgood = new int[nplanes];
339 ygood = new int[nplanes];
340 multmaxx = new int[nplanes];
341 multmaxy = new int[nplanes];
342 xm = new float[nplanes];
343 ym = new float[nplanes];
344 zm = new float[nplanes];
345 xma = new float[nplanes];
346 yma = new float[nplanes];
347 zma = new float[nplanes];
348 xmb = new float[nplanes];
349 ymb = new float[nplanes];
350 zmb = new float[nplanes];
351 resx = new float[nplanes];
352 resy = new float[nplanes];
353 xv = new float[nplanes];
354 yv = new float[nplanes];
355 zv = new float[nplanes];
356 axv = new float[nplanes];
357 ayv = new float[nplanes];
358 dedx_x = new float[nplanes];
359 dedx_y = new float[nplanes];
360
361
362 }
363 //--------------------------------------
364 //
365 //
366 //--------------------------------------
367 bool ExtTrack::SetZ(int ip,float zmeas){
368 if(ip<0 || ip>=nplanes){
369 cout << " ExtTrack::SetZ("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
370 return false;
371 }
372 zm[ip] = zmeas;
373 zma[ip] = zmeas;
374 zmb[ip] = zmeas;
375 return true;
376
377 };
378
379 bool ExtTrack::ResetXY(int ip){
380
381 // cout << " ExtTrack::ResetXY("<<ip<<") "<<endl;
382
383 if(ip<0 || ip>=nplanes){
384 cout << " ExtTrack::ResetXY("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
385 return false;
386 }
387 xm[ip] = -100.;
388 ym[ip] = -100.;
389 xma[ip] = -100.;
390 yma[ip] = -100.;
391 xmb[ip] = -100.;
392 ymb[ip] = -100.;
393 resx[ip] = 1000.;
394 resy[ip] = 1000.;
395 dedx_x[ip] = 0.;
396 dedx_y[ip] = 0.;
397 xgood[ip] = 0;
398 ygood[ip] = 0;
399 multmaxx[ip] = 0;
400 multmaxy[ip] = 0;
401
402 return true;
403
404 };
405
406 bool ExtTrack::SetXY(int ip,float xmeas, float ymeas, float rx, float ry){
407
408 if(ip<0 || ip>=nplanes){
409 cout << " ExtTrack::SetXY("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
410 return false;
411 }
412
413 xm[ip] = xmeas;
414 ym[ip] = ymeas;
415 xma[ip] = xmeas;
416 yma[ip] = ymeas;
417 xmb[ip] = xmeas;
418 ymb[ip] = ymeas;
419 resx[ip] = rx;
420 resy[ip] = ry;
421 // dedx_x[ip] = dedxx;
422 // dedx_y[ip] = dedxy;
423 xgood[ip] = 1;
424 ygood[ip] = 1;
425
426 return true;
427
428 };
429 bool ExtTrack::SetX(int ip,float xa, float xb, float ya, float yb,float res ){
430
431 // cout << " ExtTrack::SetX("<<ip<<",...) -- nplanes"<<nplanes<<endl;
432
433 if(ip<0 || ip>=nplanes){
434 cout << " ExtTrack::SetX("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
435 return false;
436 }
437 xm[ip] = -100.;
438 ym[ip] = -100.;
439 xma[ip] = xa;
440 yma[ip] = ya;
441 xmb[ip] = xb;
442 ymb[ip] = yb;
443 resx[ip] = res;
444 resy[ip] = 1000.;
445 xgood[ip] = 1;
446 ygood[ip] = 0;
447 // dedx_x[ip] = dedx;
448 // dedx_y[ip] = 0.;
449 return true;
450 };
451 bool ExtTrack::SetY(int ip,float xa, float xb, float ya, float yb,float res ){
452
453 // cout << " ExtTrack::SetY("<<ip<<",...) -- nplanes"<<nplanes<<endl;
454
455 if(ip<0 || ip>=nplanes){
456 cout << " ExtTrack::SetX("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
457 return false;
458 }
459 xm[ip] = -100.;
460 ym[ip] = -100.;
461 xma[ip] = xa;
462 yma[ip] = ya;
463 xmb[ip] = xb;
464 ymb[ip] = yb;
465 resx[ip] = 1000.;
466 resy[ip] = res;
467 xgood[ip] = 0;
468 ygood[ip] = 1;
469 // dedx_x[ip] = 0.;
470 // dedx_y[ip] = dedx;
471 return true;
472 };
473
474 bool ExtTrack::SetXGood(int ip,int icl_piu_uno, int il, int is ){
475
476 // cout << " ExtTrack::SetXGood("<<ip<<",...) -- nplanes"<<nplanes<<endl;
477
478 if(ip<0 || ip>=nplanes){
479 cout << " ExtTrack::SetXGood("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
480 return false;
481 }
482 xgood[ip]=(il+1)*100000000+(is+1)*10000000+icl_piu_uno;
483 return true;
484 };
485 bool ExtTrack::SetYGood(int ip,int icl_piu_uno, int il, int is ){
486 if(ip<0 || ip>=nplanes){
487 cout << " ExtTrack::SetYGood("<<ip<<",...) -- nplanes"<<nplanes<<endl;
488 cout << " ExtTrack::SetYGood("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
489 return false;
490 }
491 ygood[ip]=(il+1)*100000000+(is+1)*10000000+icl_piu_uno;
492 return true;
493 };
494
495 //--------------------------------------
496 //
497 //
498 //--------------------------------------
499 /**
500 * Clear the track. If option C is specified, deallocates the hit vectors.
501 */
502 void ExtTrack::Clear(Option_t* option){
503
504 // cout << " ExtTrack::Clear( "<<option<<" )"<<this<<endl;
505
506 Reset();
507
508 if (option && option[0] == 'C') {
509
510 // cout << " uccidi!! "<<endl;
511
512 if(xgood) delete [] xgood;
513 if(ygood) delete [] ygood;
514 if(multmaxx) delete [] multmaxx;
515 if(multmaxy) delete [] multmaxy;
516 if(xm) delete [] xm;
517 if(ym) delete [] ym;
518 if(zm) delete [] zm;
519 if(xma) delete [] xma;
520 if(yma) delete [] yma;
521 if(zma) delete [] zma;
522 if(xmb) delete [] xmb;
523 if(ymb) delete [] ymb;
524 if(zmb) delete [] zmb;
525 if(resx) delete [] resx;
526 if(resy) delete [] resy;
527 if(xv) delete [] xv;
528 if(yv) delete [] yv;
529 if(zv) delete [] zv;
530 if(axv) delete [] axv;
531 if(ayv) delete [] ayv;
532 if(dedx_x)delete [] dedx_x;
533 if(dedx_y)delete [] dedx_y;
534 //}
535 // delete [] xgood;
536 // delete [] ygood;
537 // delete [] multmaxx;
538 // delete [] multmaxy;
539 // delete [] xm;
540 // delete [] ym;
541 // delete [] zm;
542 // delete [] xma;
543 // delete [] yma;
544 // delete [] zma;
545 // delete [] xmb;
546 // delete [] ymb;
547 // delete [] zmb;
548 // delete [] resx;
549 // delete [] resy;
550 // delete [] xv;
551 // delete [] yv;
552 // delete [] zv;
553 // delete [] axv;
554 // delete [] ayv;
555 // delete [] dedx_x;
556 // delete [] dedx_y;
557
558 xgood = NULL;
559 ygood = NULL;
560 multmaxx = NULL;
561 multmaxy = NULL;
562 xm = NULL;
563 ym = NULL;
564 zm = NULL;
565 xma = NULL;
566 yma = NULL;
567 zma = NULL;
568 xmb = NULL;
569 ymb = NULL;
570 zmb = NULL;
571 resx = NULL;
572 resy = NULL;
573 xv = NULL;
574 yv = NULL;
575 zv = NULL;
576 axv = NULL;
577 ayv = NULL;
578 dedx_x = NULL;
579 dedx_y = NULL;
580
581 nplanes = 0;
582
583 }
584
585 };
586
587 void ExtTrack::Delete(){
588 Clear("C");
589 // delete [] xGF;
590 // delete [] yGF;
591 }
592
593 //--------------------------------------
594 //
595 //
596 //--------------------------------------
597
598 // /**
599 // * Set the position measurements
600 // */
601 // void ExtTrack::SetMeasure(Double_t *xmeas, Double_t *ymeas, Double_t *zmeas){
602 // for(int i=0; i<nplanes; i++) xm[i]= (xmeas? *xmeas++ : 0.);
603 // for(int i=0; i<nplanes; i++) ym[i]= (ymeas? *ymeas++ : 0.);
604 // for(int i=0; i<nplanes; i++) zm[i]= (ymeas? *zmeas++ : 0.);
605 // }
606 // /**
607 // * Set the TrkTrack position resolution
608 // */
609 // void ExtTrack::SetResolution(Double_t *rx, Double_t *ry){
610 // for(int i=0; i<nplanes; i++) resx[i]= (rx? *rx++ : 0.);
611 // for(int i=0; i<nplanes; i++) resy[i]= (ry? *ry++ : 0.);
612 // }
613 // /**
614 // * Set the TrkTrack good measurement
615 // */
616 // void ExtTrack::SetGood(int *xg, int *yg){
617
618 // for(int i=0; i<nplanes; i++) xgood[i]=(xg? *xg++ : 0 );
619 // for(int i=0; i<nplanes; i++) ygood[i]=(yg? *yg++ : 0 );
620 // }
621
622 //--------------------------------------
623 //
624 //
625 //--------------------------------------
626 void ExtTrack::Dump(){
627 cout << endl << "========== Track " ;
628 cout << endl << "zini : "<< zini;
629 cout << endl << "al : "; for(int i=0; i<5; i++)cout << al[i] << " ";
630 cout << endl << "chi^2 : "<< chi2;
631 cout << endl << "n.step : "<< nstep;
632 cout << endl << "xgood : "; for(int i=0; i<nplanes; i++)cout << XGood(i) ;
633 cout << endl << "ygood : "; for(int i=0; i<nplanes; i++)cout << YGood(i) ;
634 cout << endl << "xm : "; for(int i=0; i<nplanes; i++)cout << xm[i] << " ";
635 cout << endl << "ym : "; for(int i=0; i<nplanes; i++)cout << ym[i] << " ";
636 cout << endl << "zm : "; for(int i=0; i<nplanes; i++)cout << zm[i] << " ";
637 cout << endl << "xma : "; for(int i=0; i<nplanes; i++)cout << xma[i] << " ";
638 cout << endl << "yma : "; for(int i=0; i<nplanes; i++)cout << yma[i] << " ";
639 cout << endl << "zma : "; for(int i=0; i<nplanes; i++)cout << zma[i] << " ";
640 cout << endl << "xmb : "; for(int i=0; i<nplanes; i++)cout << xmb[i] << " ";
641 cout << endl << "ymb : "; for(int i=0; i<nplanes; i++)cout << ymb[i] << " ";
642 cout << endl << "zmb : "; for(int i=0; i<nplanes; i++)cout << zmb[i] << " ";
643 cout << endl << "xv : "; for(int i=0; i<nplanes; i++)cout << xv[i] << " ";
644 cout << endl << "yv : "; for(int i=0; i<nplanes; i++)cout << yv[i] << " ";
645 cout << endl << "zv : "; for(int i=0; i<nplanes; i++)cout << zv[i] << " ";
646 cout << endl << "resx : "; for(int i=0; i<nplanes; i++)cout << resx[i] << " ";
647 cout << endl << "resy : "; for(int i=0; i<nplanes; i++)cout << resy[i] << " ";
648 cout << endl << "coval : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" ";
649 cout << endl << " "; for(int i=0; i<5; i++)cout << coval[1][i]<<" ";
650 cout << endl << " "; for(int i=0; i<5; i++)cout << coval[2][i]<<" ";
651 cout << endl << " "; for(int i=0; i<5; i++)cout << coval[3][i]<<" ";
652 cout << endl << " "; for(int i=0; i<5; i++)cout << coval[4][i]<<" ";
653 cout << endl << "MDR : "<< (coval[4][4]>0. ? 1./sqrt(coval[4][4]) : 0. ) ;
654 cout << endl << "dedx_x : "; for(int i=0; i<nplanes; i++)cout << dedx_x[i] << " ";
655 cout << endl << "dedx_y : "; for(int i=0; i<nplanes; i++)cout << dedx_y[i] << " ";
656
657 cout << endl << "=================";
658 cout << endl;
659 }
660 /**
661 * Reset the fit parameters
662 */
663 void ExtTrack::ResetFit(){
664
665 // cout << " ExtTrack::ResetFit() "<<endl;
666 for(int i=0; i<5; i++) al[i]=-9999.;
667 chi2=0.;
668 nstep=0;
669
670 for(int i=0; i<5; i++) {
671 for(int j=0; j<5; j++) coval[i][j]=0.;
672 }
673
674 for(int i=0; i<nplanes; i++){
675 xv[i] = 0.;
676 yv[i] = 0.;
677 zv[i] = 0.;
678 axv[i] = 0.;
679 ayv[i] = 0.;
680 }
681
682 int ngf = TrkParams::nGF;
683 for(int i=0; i<ngf; i++){
684 xGF[i] = 0.;
685 yGF[i] = 0.;
686 }
687
688
689 }
690 /**
691 * Reset x/y measured coordinates (NB it does not reset z coordinates)
692 */
693 void ExtTrack::ResetXY(){
694
695 for(int i=0; i<nplanes; i++)ResetXY(i);
696
697
698 }
699 /**
700 * \brief Tracking method. It calls F77 mini routine.
701 *
702 * @see void TrkTrack::Fit(double, int&, int, int)
703 */
704 void ExtTrack::Fit(double pfixed, int& fail, int iprint){
705
706 TrkParams::Load(1);
707 if( !TrkParams::IsLoaded(1) ){
708 cout << "void ExtTrack::Fit(double,int&,int,int) ---ERROR--- m.field not loaded "<<endl;
709 return;
710 }
711
712 float al_ini[] = {0.,0.,0.,0.,0.};
713
714 fail = 0;
715
716 FillMiniStruct(); //fill F77 common
717
718
719 // if fit variables have been reset, evaluate the initial guess
720 if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.){
721 if(iprint==2)cout << "calling guessext_()"<<endl;
722 guessext_();
723 }
724 // --------------------- free momentum
725 if(pfixed==0.) {
726 exttrack_.pfixed=0.;
727 }
728 // --------------------- fixed momentum
729 if(pfixed!=0.) {
730 al[4]=1./pfixed;
731 exttrack_.pfixed=pfixed;
732 }
733
734 // store temporarily the initial guess
735 for(int i=0; i<5; i++) al_ini[i]=exttrack_.al[i];
736
737 if(iprint==2){
738 cout<<endl<<" init. guess: ";
739 for(int i=0; i<5; i++)cout << al_ini[i]<<" ";
740 cout<<endl;
741 }
742 // ------------------------------------------
743 // call mini routine
744 // ------------------------------------------
745 int istep=0;
746 int ifail=0;
747 if(iprint==2)cout << "calling miniext_("<<istep<<","<<ifail<<", "<<iprint<<")"<<endl;
748 miniext_(&istep,&ifail, &iprint);
749 if(ifail!=0) {
750 if(iprint)cout << "ERROR: ifail= " << ifail << endl;
751 fail = 1;
752 }
753 if(chi2!=chi2){
754 if(iprint)cout << "ERROR: chi2= " << chi2 << endl;
755 ResetFit();
756 fail = 1;
757 }
758 // ------------------------------------------
759
760 SetFromMiniStruct(); //copy content of F77 common into the ExtTrack object
761
762 if(fail){
763 if(iprint)cout << " >>>> fit failed "<<endl;
764 for(int i=0; i<5; i++) al[i]=al_ini[i];
765 }
766
767 int ngf = TrkParams::nGF;
768 float* zgf = TrkParams::zGF;
769 Trajectory tj = Trajectory(ngf,zgf);
770 tj.DoTrack(al,zini);
771 for(int i=0; i<14; i++){
772 xGF[i] = tj.x[i];
773 yGF[i] = tj.y[i];
774 }
775
776
777 };
778
779 /**
780 * Method to fill minimization-routine common
781 * NB (zini=23.5 is implicitelly set)
782 */
783 void ExtTrack::FillMiniStruct(cMiniExtTrack& track){
784
785 track.nplanes = nplanes;
786 for(int i=0; i<nplanes; i++){
787
788 track.xgood[i] = (double) XGood(i);
789 track.ygood[i] = (double) YGood(i);
790
791 track.xm[i] = (double) xm[i];
792 track.ym[i] = (double) ym[i];
793 track.zm[i] = (double) zm[i];
794
795
796 track.xm_a[i] = (double) xma[i];
797 track.xm_b[i] = (double) xmb[i];
798 track.ym_a[i] = (double) yma[i];
799 track.ym_b[i] = (double) ymb[i];
800 track.zm_a[i] = (double) zma[i];
801 track.zm_b[i] = (double) zmb[i];
802
803 track.resx[i] = (double) resx[i];
804 track.resy[i] = (double) resy[i];
805 }
806
807 for(int i=0; i<5; i++) track.al[i] = (double) al[i];
808
809 track.zini = (double) zini;
810
811
812 }
813 /**
814 * Method to set values from minimization-routine common
815 */
816 void ExtTrack::SetFromMiniStruct(cMiniExtTrack *track){
817
818 for(int i=0; i<5; i++) {
819 al[i] = (float) (track->al[i]);
820 for(int j=0; j<5; j++) coval[i][j] = (float) (track->cov[i][j]);
821 }
822 chi2 = (float) (track->chi2);
823 nstep = (float) (track->nstep);
824 for(int i=0; i<nplanes; i++){
825 xv[i] = (float) (track->xv[i]);
826 yv[i] = (float) (track->yv[i]);
827 zv[i] = (float) (track->zv[i]);
828 xm[i] = (float) (track->xm[i]);
829 ym[i] = (float) (track->ym[i]);
830 zm[i] = (float) (track->zm[i]);
831 xma[i] = (float) (track->xm_a[i]);
832 yma[i] = (float) (track->ym_a[i]);
833 zma[i] = (float) (track->zm_a[i]);
834 xmb[i] = (float) (track->xm_b[i]);
835 ymb[i] = (float) (track->ym_b[i]);
836 zmb[i] = (float) (track->zm_b[i]);
837 axv[i] = (float) (track->axv[i]);
838 ayv[i] = (float) (track->ayv[i]);
839 }
840 zini = (float) (track->zini);
841 }
842 /**
843 * Method to retrieve ID (0,1,...) of x-cluster (if any) associated to this track.
844 * If no cluster is associated, ID=-1.
845 *
846 */
847 int ExtTrack::GetClusterX_ID(int ip){
848 if(ip<0 || ip>=nplanes)return -1;
849 return ((int)fabs(xgood[ip]))%10000000-1;
850 };
851 /**
852 * Method to retrieve ID (0-xxx) of y-cluster (if any) associated to this track.
853 *
854 */
855 int ExtTrack::GetClusterY_ID(int ip){
856 if(ip<0 || ip>=nplanes)return -1;
857 return ((int)fabs(ygood[ip]))%10000000-1;
858 };
859
860 /**
861 * Method to retrieve the dE/dx measured on a tracker view.
862 * @param ip plane (0-5)
863 * @param iv view (0=x 1=y)
864 */
865 Float_t ExtTrack::GetDEDX(int ip, int iv){
866 if(iv==0 && ip>=0 && ip<nplanes)return fabs(dedx_x[ip]);
867 else if(iv==1 && ip>=0 && ip<nplanes)return fabs(dedx_y[ip]);
868 else {
869 cout << "TrkTrack::GetDEDX(int ip, int iv) -- wrong input parameters "<<ip<<iv<<endl;
870 return 0.;
871 }
872 }
873 /**
874 * Method to evaluate the dE/dx measured on a tracker plane.
875 * The two measurements on x- and y-view are averaged.
876 * @param ip plane (0-5)
877 */
878 Float_t ExtTrack::GetDEDX(int ip){
879 if( (Int_t)XGood(ip)+(Int_t)YGood(ip) == 0 ) return 0;
880 return (GetDEDX(ip,0)+GetDEDX(ip,1))/((Int_t)XGood(ip)+(Int_t)YGood(ip));
881 };
882
883 /**
884 * Method to evaluate the dE/dx averaged over all planes.
885 */
886 Float_t ExtTrack::GetDEDX(){
887 Float_t dedx=0;
888 for(Int_t ip=0; ip<nplanes; ip++)dedx+=GetDEDX(ip,0)*XGood(ip)+GetDEDX(ip,1)*YGood(ip);
889 dedx = dedx/(GetNX()+GetNY());
890 return dedx;
891 };
892
893
894 /**
895 * Method to evaluate the dE/dx averaged over all X views.
896 */
897 Float_t ExtTrack::GetDEDXX(bool cutSat){
898 Float_t dedx=0;
899 Int_t np=0;
900 for(Int_t ip=0; ip<nplanes; ip++){
901 dedx += GetDEDX(ip,0)*XGood(ip)*(cutSat ? !IsSaturated(ip,0) : 1);
902 np += XGood(ip)*(cutSat ? !IsSaturated(ip,0) : 1);
903 }
904 if(np>0)return dedx/np;
905 return -1;
906 };
907 /**
908 * Method to evaluate the dE/dx averaged over all Y views.
909 */
910 Float_t ExtTrack::GetDEDXY(bool cutSat){
911 Float_t dedx=0;
912 Int_t np=0;
913 for(Int_t ip=0; ip<nplanes; ip++){
914 dedx += GetDEDX(ip,1)*YGood(ip)*(cutSat ? !IsSaturated(ip,1) : 1);
915 np += YGood(ip)*(cutSat ? !IsSaturated(ip,1) : 1);
916 }
917 if(np>0)return dedx/np;
918 return -1;
919 };
920
921 /**
922 * Returns 1 if the cluster on a tracker view includes bad strips
923 * (at least one bad strip among the four strip used by p.f.a.)
924 * @param ip plane (0-5)
925 * @param iv view (0=x 1=y)
926 */
927 Bool_t ExtTrack::IsBad(int ip,int iv){
928 if(ip>=6&&ip<8)return false;
929 if(iv==0 && ip>=0 && ip<6)return (xgood[ip]<0) ;
930 else if(iv==1 && ip>=0 && ip<6)return (ygood[ip]<0) ;
931 else {
932 cout << "ExtTrack::IsBad(int ip, int iv) -- wrong input parameters "<<ip<<iv<<endl;
933 return 0.;
934 }
935 };
936
937 /**
938 * Method to retrieve the ladder, if assigned
939 */
940 int ExtTrack::GetLadder(int ip){
941 if(XGood(ip))return (int)fabs(xgood[ip]/100000000)-1;
942 if(YGood(ip))return (int)fabs(ygood[ip]/100000000)-1;
943 return -1;
944 };
945 /**
946 * Method to retrieve the sensor, if assigned
947 */
948 int ExtTrack::GetSensor(int ip){
949 if(XGood(ip))return (int)((int)fabs(xgood[ip]/10000000)%10)-1;
950 if(YGood(ip))return (int)((int)fabs(ygood[ip]/10000000)%10)-1;
951 return -1;
952 };
953 /**
954 *
955 */
956 Int_t ExtTrack::GetClusterX_Multiplicity(int ip){
957 if(ip<0 || ip>=nplanes){
958 cout << " ExtTrack::GetClusterX_Multiplicity("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
959 return -1;
960 }
961 return (Int_t)(multmaxx[ip]/10000);
962 }
963 /**
964 *
965 */
966
967 Int_t ExtTrack::GetClusterY_Multiplicity(int ip){
968 if(ip<0 || ip>=nplanes){
969 cout << " ExtTrack::GetClusterY_Multiplicity("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
970 return -1;
971 }
972 return (Int_t)(multmaxy[ip]/10000);
973 }
974 /**
975 *
976 */
977
978 Int_t ExtTrack::GetClusterX_MaxStrip(int ip){
979 if(ip<0 || ip>=nplanes){
980 cout << " ExtTrack::GetClusterX_MaxStrip("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
981 return -1;
982 }
983 return (Int_t)(multmaxx[ip]%10000);
984 }
985 /**
986 *
987 */
988
989 Int_t ExtTrack::GetClusterY_MaxStrip(int ip){
990 if(ip<0 || ip>=nplanes){
991 cout << " ExtTrack::GetClusterY_MaxStrip("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
992 return -1;
993 }
994 return (Int_t)(multmaxy[ip]%10000);
995 }
996 Float_t ExtTrack::GetRigidity(){
997 Float_t rig=0;
998 if(chi2>0)rig=1./al[4];
999 if(rig<0) rig=-rig;
1000 return rig;
1001 };
1002 //
1003 Float_t ExtTrack::GetDeflection(){
1004 Float_t def=0;
1005 if(chi2>0)def=al[4];
1006 return def;
1007 };
1008
1009
1010 //
1011 // all that follows: EM porting from TrkLevel2
1012 //
1013 Bool_t ExtTrack::IsInsideAcceptance(float toll){
1014 int ngf = TrkParams::nGF;
1015 float* zgf = TrkParams::zGF;
1016
1017 //this is a patch, to assign correctly xGF and yGF cordinates.
1018
1019 Trajectory tj = Trajectory(ngf,zgf);
1020 tj.DoTrack(al,zini);
1021
1022 //NB the folowing assignment is NOT permanent!!!
1023 for(int i=0; i<14; i++){
1024 xGF[i] = tj.x[i];
1025 yGF[i] = tj.y[i];
1026 }
1027
1028 for(int i=0; i<ngf; i++){
1029 //
1030 // cout << endl << TrkParams::GF_element[i];
1031 // if(
1032 // TrkParams::GF_element[i].CompareTo("S11") &&
1033 // TrkParams::GF_element[i].CompareTo("S12") &&
1034 // TrkParams::GF_element[i].CompareTo("S21") &&
1035 // TrkParams::GF_element[i].CompareTo("S22") &&
1036 // TrkParams::GF_element[i].CompareTo("T1") &&
1037 // TrkParams::GF_element[i].CompareTo("CUF") &&
1038 // TrkParams::GF_element[i].CompareTo("T2") &&
1039 // TrkParams::GF_element[i].CompareTo("T3") &&
1040 // TrkParams::GF_element[i].CompareTo("T4") &&
1041 // TrkParams::GF_element[i].CompareTo("T5") &&
1042 // TrkParams::GF_element[i].CompareTo("CLF") &&
1043 // TrkParams::GF_element[i].CompareTo("T6") &&
1044 // TrkParams::GF_element[i].CompareTo("S31") &&
1045 // TrkParams::GF_element[i].CompareTo("S32") &&
1046 // true)continue;
1047 // apply condition only within the cavity
1048 // cout << " -- "<<xGF[i]<<" "<<yGF[i];
1049 if(
1050 // xGF[i] <= TrkParams::xGF_min[i] + toll ||
1051 // xGF[i] >= TrkParams::xGF_max[i] - toll ||
1052 // yGF[i] <= TrkParams::yGF_min[i] + toll ||
1053 // yGF[i] >= TrkParams::yGF_max[i] - toll ||
1054 tj.x[i] <= TrkParams::xGF_min[i] + toll ||
1055 tj.x[i] >= TrkParams::xGF_max[i] - toll ||
1056 tj.y[i] <= TrkParams::yGF_min[i] + toll ||
1057 tj.y[i] >= TrkParams::yGF_max[i] - toll ||
1058 false){
1059
1060 return false;
1061 }
1062 }
1063 return true;
1064 }
1065
1066 /**
1067 * Returns the reduced chi-square of track x-projection
1068 */
1069 Float_t ExtTrack::GetChi2X(){
1070 float chiq=0;
1071 for(int ip=0; ip<nplanes; ip++)if(XGood(ip))chiq+= pow((xv[ip]-xm[ip])/resx[ip],2.);
1072 if(GetNX()>3)chiq=chiq/(GetNX()-3);
1073 else chiq=0;
1074 // if(chiq==0)cout << " Float_t ExtTrack::GetChi2X() -- WARNING -- value not defined "<<chiq<<endl;
1075 return chiq;
1076 }
1077 /**
1078 * Returns the reduced chi-square of track y-projection
1079 */
1080 Float_t ExtTrack::GetChi2Y(){
1081 float chiq=0;
1082 for(int ip=0; ip<nplanes; ip++)if(YGood(ip))chiq+= pow((yv[ip]-ym[ip])/resy[ip],2.);
1083 if(GetNY()>2)chiq=chiq/(GetNY()-2);
1084 else chiq=0;
1085 // if(chiq==0)cout << " Float_t ExtTrack::GetChi2Y() -- WARNING -- value not defined "<<chiq<<endl;
1086 return chiq;
1087 }
1088
1089 /**
1090 * Returns the track "lever-arm" on the x view, defined as the distance (in planes) between
1091 * the upper and lower x measurements (the maximum value of lever-arm is 6).
1092 */
1093 Int_t ExtTrack::GetLeverArmX(){
1094 int first_plane = -1;
1095 int last_plane = -1;
1096 for(Int_t ip=0; ip<nplanes; ip++){
1097 if( XGood(ip) && first_plane == -1 )first_plane = ip;
1098 if( XGood(ip) && first_plane != -1 )last_plane = ip;
1099 }
1100 if( first_plane == -1 || last_plane == -1){
1101 cout<< "Int_t ExtTrack::GetLeverArmX() -- XGood(ip) always false ??? "<<endl;
1102 return 0;
1103 }
1104 return (last_plane-first_plane+1);
1105 }
1106 /**
1107 * Returns the track "lever-arm" on the y view, defined as the distance (in planes) between
1108 * the upper and lower y measurements (the maximum value of lever-arm is 6).
1109 */
1110 Int_t ExtTrack::GetLeverArmY(){
1111 int first_plane = -1;
1112 int last_plane = -1;
1113 for(Int_t ip=0; ip<nplanes; ip++){
1114 if( YGood(ip) && first_plane == -1 )first_plane = ip;
1115 if( YGood(ip) && first_plane != -1 )last_plane = ip;
1116 }
1117 if( first_plane == -1 || last_plane == -1){
1118 cout<< "Int_t ExtTrack::GetLeverArmY() -- YGood(ip) always false ??? "<<endl;
1119 return 0;
1120 }
1121 return (last_plane-first_plane+1);
1122 }
1123 /**
1124 * Returns the track "lever-arm" on the x+y view, defined as the distance (in planes) between
1125 * the upper and lower x,y (couple) measurements (the maximum value of lever-arm is 6).
1126 */
1127 Int_t ExtTrack::GetLeverArmXY(){
1128 int first_plane = -1;
1129 int last_plane = -1;
1130 for(Int_t ip=0; ip<nplanes; ip++){
1131 if( XGood(ip) && YGood(ip) && first_plane == -1 )first_plane = ip;
1132 if( XGood(ip) && YGood(ip) && first_plane != -1 )last_plane = ip;
1133 }
1134 if( first_plane == -1 || last_plane == -1){
1135 // cout<< "Int_t ExtTrack::GetLeverArmXY() -- XGood(ip)*YGood(ip) always false ??? "<<endl;
1136 return 0;
1137 }
1138 return (last_plane-first_plane+1);
1139 }
1140
1141
1142 /**
1143 * Returns 1 if the signal on a tracker view is saturated.
1144 * @param ip plane (0-5)
1145 * @param iv view (0=x 1=y)
1146 */
1147 Bool_t ExtTrack::IsSaturated(int ip,int iv){
1148 if(ip>=6 && ip<8)return false;
1149 if(iv==0 && ip>=0 && ip<6)return (dedx_x[ip]<0) ;
1150 else if(iv==1 && ip>=0 && ip<6)return (dedx_y[ip]<0) ;
1151 else {
1152 cout << "ExtTrack::IsSaturated(int ip, int iv) -- wrong input parameters "<<ip<<iv<<endl;
1153 return 0.;
1154 }
1155 };
1156 /**
1157 * Returns 1 if either the x or the y signal on a tracker plane is saturated.
1158 * @param ip plane (0-5)
1159 */
1160 Bool_t ExtTrack::IsSaturated(int ip){
1161 return (IsSaturated(ip,0)||IsSaturated(ip,1));
1162 };
1163 /**
1164 * Returns 1 if there is at least a saturated signal along the track.
1165 */
1166 Bool_t ExtTrack::IsSaturated(){
1167 for(int ip=0; ip<nplanes; ip++)for(int iv=0; iv<2; iv++)if(IsSaturated(ip,iv))return true;
1168 return false;
1169 }
1170
1171
1172
1173
1174
1175 ClassImp(ExtTrack);

  ViewVC Help
Powered by ViewVC 1.1.23