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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.12 - (show annotations) (download)
Tue Oct 24 07:28:38 2006 UTC (18 years, 1 month ago) by mocchiut
Branch: MAIN
Changes since 1.11: +1 -1 lines
Compilation bug  + bug in LoadRunInfoTree method fixed

1 /**
2 * \file TrkLevel2.cpp
3 * \author Elena Vannuccini
4 */
5 #include <TrkLevel2.h>
6 #include <iostream>
7 using namespace std;
8 //......................................
9 // F77 routines
10 //......................................
11 extern "C" {
12 void dotrack_(int*, double*, double*, double*, double*, int*);
13 void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);
14 int readb_(const char*);
15 }
16 //--------------------------------------
17 //
18 //
19 //--------------------------------------
20 TrkTrack::TrkTrack(){
21 seqno = -1;
22 image = -1;
23 chi2 = 0;
24 nstep = 0;
25 for(int it1=0;it1<5;it1++){
26 al[it1] = 0;
27 for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;
28 };
29 for(int ip=0;ip<6;ip++){
30 xgood[ip] = 0;
31 ygood[ip] = 0;
32 xm[ip] = 0;
33 ym[ip] = 0;
34 zm[ip] = 0;
35 resx[ip] = 0;
36 resy[ip] = 0;
37 xv[ip] = 0;
38 yv[ip] = 0;
39 zv[ip] = 0;
40 axv[ip] = 0;
41 ayv[ip] = 0;
42 dedx_x[ip] = 0;
43 dedx_y[ip] = 0;
44 // clx[ip] = 0;
45 // cly[ip] = 0;
46 };
47 clx = new TRefArray(6,0);
48 cly = new TRefArray(6,0);
49 };
50 //--------------------------------------
51 //
52 //
53 //--------------------------------------
54 TrkTrack::TrkTrack(const TrkTrack& t){
55 seqno = t.seqno;
56 image = t.image;
57 chi2 = t.chi2;
58 nstep = t.nstep;
59 for(int it1=0;it1<5;it1++){
60 al[it1] = t.al[it1];
61 for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
62 };
63 for(int ip=0;ip<6;ip++){
64 xgood[ip] = t.xgood[ip];
65 ygood[ip] = t.ygood[ip];
66 xm[ip] = t.xm[ip];
67 ym[ip] = t.ym[ip];
68 zm[ip] = t.zm[ip];
69 resx[ip] = t.resx[ip];
70 resy[ip] = t.resy[ip];
71 xv[ip] = t.xv[ip];
72 yv[ip] = t.yv[ip];
73 zv[ip] = t.zv[ip];
74 axv[ip] = t.axv[ip];
75 ayv[ip] = t.ayv[ip];
76 dedx_x[ip] = t.dedx_x[ip];
77 dedx_y[ip] = t.dedx_y[ip];
78 // clx[ip] = 0;//<<<<pointer
79 // cly[ip] = 0;//<<<<pointer
80 };
81 clx = new TRefArray(*(t.clx));
82 cly = new TRefArray(*(t.cly));
83
84 };
85 //--------------------------------------
86 //
87 //
88 //--------------------------------------
89 /**
90 * Evaluates the trajectory in the apparatus associated to the track.
91 * It integrates the equations of motion in the magnetic field. The magnetic field should be previously loaded ( by calling TrkLevel2::LoadField() ), otherwise an error message is returned.
92 * @param t pointer to an object of the class Trajectory,
93 * which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ).
94 * @return error flag.
95 */
96 int TrkTrack::DoTrack(Trajectory* t){
97
98 double *dxout = new double[t->npoint];
99 double *dyout = new double[t->npoint];
100 double *dzin = new double[t->npoint];
101 double dal[5];
102
103 int ifail = 0;
104
105 for (int i=0; i<5; i++) dal[i] = (double)al[i];
106 for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i];
107
108 dotrack_(&(t->npoint),dzin,dxout,dyout,dal,&ifail);
109
110 for (int i=0; i<t->npoint; i++){
111 t->x[i] = (float)*dxout++;
112 t->y[i] = (float)*dyout++;
113 }
114
115 // delete [] dxout;
116 // delete [] dyout;
117 // delete [] dzin;
118
119 return ifail;
120 };
121 //--------------------------------------
122 //
123 //
124 //--------------------------------------
125 /**
126 * Evaluates the trajectory in the apparatus associated to the track.
127 * It integrates the equations of motion in the magnetic field. The magnetic field should be previously loaded ( by calling TrkLevel2::LoadField() ), otherwise an error message is returned.
128 * @param t pointer to an object of the class Trajectory,
129 * which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ).
130 * @return error flag.
131 */
132 int TrkTrack::DoTrack2(Trajectory* t){
133
134 double *dxout = new double[t->npoint];
135 double *dyout = new double[t->npoint];
136 double *dthxout = new double[t->npoint];
137 double *dthyout = new double[t->npoint];
138 double *dtlout = new double[t->npoint];
139 double *dzin = new double[t->npoint];
140 double dal[5];
141
142 int ifail = 0;
143
144 for (int i=0; i<5; i++) dal[i] = (double)al[i];
145 for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i];
146
147 dotrack2_(&(t->npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail);
148
149 for (int i=0; i<t->npoint; i++){
150 t->x[i] = (float)*dxout++;
151 t->y[i] = (float)*dyout++;
152 t->thx[i] = (float)*dthxout++;
153 t->thy[i] = (float)*dthyout++;
154 t->tl[i] = (float)*dtlout++;
155 }
156
157 // delete [] dxout;
158 // delete [] dyout;
159 // delete [] dzin;
160
161 return ifail;
162 };
163 //--------------------------------------
164 //
165 //
166 //--------------------------------------
167 //float TrkTrack::BdL(){
168 //};
169 //--------------------------------------
170 //
171 //
172 //--------------------------------------
173 Float_t TrkTrack::GetRigidity(){
174 Float_t rig=0;
175 if(chi2>0)rig=1./al[4];
176 if(rig<0) rig=-rig;
177 return rig;
178 };
179 //
180 Float_t TrkTrack::GetDeflection(){
181 Float_t def=0;
182 if(chi2>0)def=al[4];
183 return def;
184 };
185 //
186 Float_t TrkTrack::GetDEDX(){
187 Float_t dedx=0;
188 for(Int_t i=0; i<6; i++)dedx+=dedx_x[i]*xgood[i]+dedx_y[i]*ygood[i];
189 dedx = dedx/(this->GetNX()+this->GetNY());
190 return dedx;
191 };
192
193 //--------------------------------------
194 //
195 //
196 //--------------------------------------
197 void TrkTrack::Dump(){
198 cout << endl << "========== Track " ;
199 cout << endl << "seq. n. : "<< seqno;
200 cout << endl << "image n. : "<< image;
201 cout << endl << "al : "; for(int i=0; i<5; i++)cout << al[i] << " ";
202 cout << endl << "chi^2 : "<< chi2;
203 cout << endl << "n.step : "<< nstep;
204 cout << endl << "xgood : "; for(int i=0; i<6; i++)cout << xgood[i] ;
205 cout << endl << "ygood : "; for(int i=0; i<6; i++)cout << ygood[i] ;
206 cout << endl << "xm : "; for(int i=0; i<6; i++)cout << xm[i] << " ";
207 cout << endl << "ym : "; for(int i=0; i<6; i++)cout << ym[i] << " ";
208 cout << endl << "zm : "; for(int i=0; i<6; i++)cout << zm[i] << " ";
209 cout << endl << "dedx_x : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " ";
210 cout << endl << "dedx_y : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " ";
211 }
212 //--------------------------------------
213 //
214 //
215 //--------------------------------------
216 void TrkTrack::Clear(){
217 seqno = -1;
218 image = -1;
219 chi2 = 0;
220 nstep = 0;
221 for(int it1=0;it1<5;it1++){
222 al[it1] = 0;
223 for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;
224 };
225 for(int ip=0;ip<6;ip++){
226 xgood[ip] = 0;
227 ygood[ip] = 0;
228 xm[ip] = 0;
229 ym[ip] = 0;
230 zm[ip] = 0;
231 resx[ip] = 0;
232 resy[ip] = 0;
233 xv[ip] = 0;
234 yv[ip] = 0;
235 zv[ip] = 0;
236 axv[ip] = 0;
237 ayv[ip] = 0;
238 dedx_x[ip] = 0;
239 dedx_y[ip] = 0;
240 // clx[ip] = 0;
241 // cly[ip] = 0;
242 };
243 clx->Clear();
244 cly->Clear();
245 };
246 //--------------------------------------
247 //
248 //
249 //--------------------------------------
250 void TrkTrack::Delete(){
251 Clear();
252 clx->Delete();
253 cly->Delete();
254 };
255 //--------------------------------------
256 //
257 //
258 //--------------------------------------
259
260 //--------------------------------------
261 //
262 //
263 //--------------------------------------
264 TrkSinglet::TrkSinglet(){
265 plane = 0;
266 coord[0] = 0;
267 coord[1] = 0;
268 sgnl = 0;
269 cls = 0;
270 };
271 //--------------------------------------
272 //
273 //
274 //--------------------------------------
275 TrkSinglet::TrkSinglet(const TrkSinglet& s){
276 plane = s.plane;
277 coord[0] = s.coord[0];
278 coord[1] = s.coord[1];
279 sgnl = s.sgnl;
280 // cls = 0;//<<<<pointer
281 cls = TRef(s.cls);
282 };
283 //--------------------------------------
284 //
285 //
286 //--------------------------------------
287 void TrkSinglet::Dump(){
288 int i=0;
289 cout << endl << "========== Singlet " ;
290 cout << endl << "plane : " << plane;
291 cout << endl << "coord[2] : "; while( i<2 && cout << coord[i] << " ") i++;
292 cout << endl << "sgnl : " << sgnl;
293 }
294 //--------------------------------------
295 //
296 //
297 //--------------------------------------
298 TrkLevel2::TrkLevel2(){
299 // good2 = -1;
300 for(Int_t i=0; i<12 ; i++){
301 // crc[i] = -1;
302 good[i] = -1;
303 };
304 Track = new TClonesArray("TrkTrack");
305 SingletX = new TClonesArray("TrkSinglet");
306 SingletY = new TClonesArray("TrkSinglet");
307
308 // PhysicalTrack = new TClonesArray("TrkTrack");
309 //sostituire con TRefArray... appena ho capito come si usa
310 }
311 //--------------------------------------
312 //
313 //
314 //--------------------------------------
315 void TrkLevel2::Dump(){
316
317 TClonesArray &t = *Track;
318 TClonesArray &sx = *SingletX;
319 TClonesArray &sy = *SingletY;
320 //
321 cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-";
322 cout << endl << "good : "; for(int i=0; i<12; i++) cout << good[i]<<" ";
323 cout << endl << "ntrk() : " << this->ntrk() ;
324 cout << endl << "nclsx() : " << this->nclsx();
325 cout << endl << "nclsy() : " << this->nclsy();
326 for(int i=0; i<this->ntrk(); i++) ((TrkTrack *)t[i])->Dump();
327 for(int i=0; i<this->nclsx(); i++) ((TrkSinglet *)sx[i])->Dump();
328 for(int i=0; i<this->nclsy(); i++) ((TrkSinglet *)sy[i])->Dump();
329 }
330 //--------------------------------------
331 //
332 //
333 //--------------------------------------
334 /**
335 * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common).
336 */
337 void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2){
338
339 // temporary objects:
340 TrkSinglet* t_singlet = new TrkSinglet();
341 TrkTrack* t_track = new TrkTrack();
342
343 // **** general variables ****
344 // good2 = l2->good2;
345 for(Int_t i=0; i<12 ; i++){
346 // crc[i] = l2->crc[i];
347 good[i] = l2->good[i];
348 };
349 // *** TRACKS ***
350 TClonesArray &t = *Track;
351 for(int i=0; i<l2->ntrk; i++){
352 t_track->seqno = i;// NBNBNBNB deve sempre essere = i
353 t_track->image = l2->image[i]-1;
354 // cout << "track "<<i<<t_track->seqno << t_track->image<<endl;
355 t_track->chi2 = l2->chi2_nt[i];
356 t_track->nstep = l2->nstep_nt[i];
357 for(int it1=0;it1<5;it1++){
358 t_track->al[it1] = l2->al_nt[i][it1];
359 for(int it2=0;it2<5;it2++)
360 t_track->coval[it1][it2] = l2->coval[i][it2][it1];
361 };
362 for(int ip=0;ip<6;ip++){
363 t_track->xgood[ip] = l2->xgood_nt[i][ip];
364 t_track->ygood[ip] = l2->ygood_nt[i][ip];
365 t_track->xm[ip] = l2->xm_nt[i][ip];
366 t_track->ym[ip] = l2->ym_nt[i][ip];
367 t_track->zm[ip] = l2->zm_nt[i][ip];
368 t_track->resx[ip] = l2->resx_nt[i][ip];
369 t_track->resy[ip] = l2->resy_nt[i][ip];
370 t_track->xv[ip] = l2->xv_nt[i][ip];
371 t_track->yv[ip] = l2->yv_nt[i][ip];
372 t_track->zv[ip] = l2->zv_nt[i][ip];
373 t_track->axv[ip] = l2->axv_nt[i][ip];
374 t_track->ayv[ip] = l2->ayv_nt[i][ip];
375 t_track->dedx_x[ip] = l2->dedx_x[i][ip];
376 t_track->dedx_y[ip] = l2->dedx_y[i][ip];
377 // t_track->clx[ip] = 0;
378 // t_track->cly[ip] = 0;
379 };
380 new(t[i]) TrkTrack(*t_track);
381 t_track->Clear();
382 };
383 // *** SINGLETS ***
384 TClonesArray &sx = *SingletX;
385 for(int i=0; i<l2->nclsx; i++){
386 t_singlet->plane = l2->planex[i];
387 t_singlet->coord[0] = l2->xs[i][0];
388 t_singlet->coord[1] = l2->xs[i][1];
389 t_singlet->sgnl = l2->signlxs[i];
390 // t_singlet->cls = 0;
391 new(sx[i]) TrkSinglet(*t_singlet);
392 t_singlet->Clear();
393 }
394 TClonesArray &sy = *SingletY;
395 for(int i=0; i<l2->nclsy; i++){
396 t_singlet->plane = l2->planey[i];
397 t_singlet->coord[0] = l2->ys[i][0];
398 t_singlet->coord[1] = l2->ys[i][1];
399 t_singlet->sgnl = l2->signlys[i];
400 // t_singlet->cls = 0;
401 new(sy[i]) TrkSinglet(*t_singlet);
402 t_singlet->Clear();
403 };
404
405 delete t_track;
406 delete t_singlet;
407 }
408 //--------------------------------------
409 //
410 //
411 //--------------------------------------
412 /**
413 * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common).
414 * Ref to Level1 data (clusters) is also set.
415 */
416 void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1){
417
418 // temporary objects:
419 TrkSinglet* t_singlet = new TrkSinglet();
420 TrkTrack* t_track = new TrkTrack();
421 // general variables
422 // good2 = l2->good2;
423 for(Int_t i=0; i<12 ; i++){
424 // crc[i] = l2->crc[i];
425 good[i] = l2->good[i];
426 };
427 // *** TRACKS ***
428 TClonesArray &t = *Track;
429 for(int i=0; i<l2->ntrk; i++){
430 t_track->seqno = i;// NBNBNBNB deve sempre essere = i
431 t_track->image = l2->image[i]-1;
432 // cout << "track "<<i<<t_track->seqno << t_track->image<<endl;
433 t_track->chi2 = l2->chi2_nt[i];
434 t_track->nstep = l2->nstep_nt[i];
435 for(int it1=0;it1<5;it1++){
436 t_track->al[it1] = l2->al_nt[i][it1];
437 for(int it2=0;it2<5;it2++)
438 t_track->coval[it1][it2] = l2->coval[i][it2][it1];
439 };
440 for(int ip=0;ip<6;ip++){
441 t_track->xgood[ip] = l2->xgood_nt[i][ip];
442 t_track->ygood[ip] = l2->ygood_nt[i][ip];
443 t_track->xm[ip] = l2->xm_nt[i][ip];
444 t_track->ym[ip] = l2->ym_nt[i][ip];
445 t_track->zm[ip] = l2->zm_nt[i][ip];
446 t_track->resx[ip] = l2->resx_nt[i][ip];
447 t_track->resy[ip] = l2->resy_nt[i][ip];
448 t_track->xv[ip] = l2->xv_nt[i][ip];
449 t_track->yv[ip] = l2->yv_nt[i][ip];
450 t_track->zv[ip] = l2->zv_nt[i][ip];
451 t_track->axv[ip] = l2->axv_nt[i][ip];
452 t_track->ayv[ip] = l2->ayv_nt[i][ip];
453 t_track->dedx_x[ip] = l2->dedx_x[i][ip];
454 t_track->dedx_y[ip] = l2->dedx_y[i][ip];
455 // cout << "traccia "<<i<<" -- "<< ip << " "<< l2->cltrx[i][ip] <<" "<< l2->cltry[i][ip] <<" "<< t_track->xgood[ip] << t_track->ygood[ip]<<endl;
456 //-----------------------------------------------------
457 // t_track->clx[ip] = l1->GetCluster(l2->cltrx[i][ip]-1);
458 // t_track->cly[ip] = l1->GetCluster(l2->cltry[i][ip]-1);
459 if(t_track->xgood[ip])t_track->clx->AddAt(l1->GetCluster(l2->cltrx[i][ip]-1),ip);
460 if(t_track->ygood[ip])t_track->cly->AddAt(l1->GetCluster(l2->cltry[i][ip]-1),ip);
461 // if(t_track->ygood[ip])cout<<" i "<<i<<" ip "<<ip<<" l2->cltry[i][ip] "<<l2->cltry[i][ip]<< " l1->GetCluster(l2->cltry[i][ip]-1) "<<l1->GetCluster(l2->cltry[i][ip]-1)<<endl;
462 // if(t_track->xgood[ip])cout<<" i "<<i<<" ip "<<ip<<" l2->cltrx[i][ip] "<<l2->cltrx[i][ip]<< " l1->GetCluster(l2->cltrx[i][ip]-1) "<<l1->GetCluster(l2->cltrx[i][ip]-1)<<endl;
463 //-----------------------------------------------------
464 };
465 new(t[i]) TrkTrack(*t_track);
466 t_track->Clear();
467 };
468 // *** SINGLETS ***
469 TClonesArray &sx = *SingletX;
470 for(int i=0; i<l2->nclsx; i++){
471 t_singlet->plane = l2->planex[i];
472 t_singlet->coord[0] = l2->xs[i][0];
473 t_singlet->coord[1] = l2->xs[i][1];
474 t_singlet->sgnl = l2->signlxs[i];
475 //-----------------------------------------------------
476 // cout << "singolo x "<<i<<" -- "<< l2->clsx[i] <<endl;
477 t_singlet->cls = l1->GetCluster(l2->clsx[i]-1);
478 // cout<<" i "<<i<<" l2->clsx[i] "<<l2->clsx[i]<< " l1->GetCluster(l2->clsx[i]-1) "<<l1->GetCluster(l2->clsx[i]-1)<<endl;
479 //-----------------------------------------------------
480 new(sx[i]) TrkSinglet(*t_singlet);
481 t_singlet->Clear();
482 }
483 TClonesArray &sy = *SingletY;
484 for(int i=0; i<l2->nclsy; i++){
485 t_singlet->plane = l2->planey[i];
486 t_singlet->coord[0] = l2->ys[i][0];
487 t_singlet->coord[1] = l2->ys[i][1];
488 t_singlet->sgnl = l2->signlys[i];
489 //-----------------------------------------------------
490 // cout << "singolo y "<<i<<" -- "<< l2->clsy[i] <<endl;
491 t_singlet->cls = l1->GetCluster(l2->clsy[i]-1);
492 // cout<<" i "<<i<<" l2->clsy[i] "<<l2->clsy[i]<< " l1->GetCluster(l2->clsy[i]-1) "<<l1->GetCluster(l2->clsy[i]-1)<<endl;
493 //-----------------------------------------------------
494 new(sy[i]) TrkSinglet(*t_singlet);
495 t_singlet->Clear();
496 };
497
498 delete t_track;
499 delete t_singlet;
500 }
501 /**
502 * Fills a struct cTrkLevel2 with values from a TrkLevel2 object (to put data into a F77 common).
503 */
504
505 void TrkLevel2::GetLevel2Struct(cTrkLevel2 *l2) const {
506
507 // general variables
508 // l2->good2 = good2 ;
509 for(Int_t i=0; i<12 ; i++){
510 // l2->crc[i] = crc[i];
511 l2->good[i] = good[i];
512 };
513 // *** TRACKS ***
514
515 l2->ntrk = Track->GetEntries();
516 for(Int_t i=0;i<l2->ntrk;i++){
517 l2->image[i] = 1 + ((TrkTrack *)Track->At(i))->image;
518 l2->chi2_nt[i] = ((TrkTrack *)Track->At(i))->chi2;
519 l2->nstep_nt[i] = ((TrkTrack *)Track->At(i))->nstep;
520 for(int it1=0;it1<5;it1++){
521 l2->al_nt[i][it1] = ((TrkTrack *)Track->At(i))->al[it1];
522 for(int it2=0;it2<5;it2++)
523 l2->coval[i][it2][it1] = ((TrkTrack *)Track->At(i))->coval[it1][it2];
524 };
525 for(int ip=0;ip<6;ip++){
526 l2->xgood_nt[i][ip] = ((TrkTrack *)Track->At(i))->xgood[ip];
527 l2->ygood_nt[i][ip] = ((TrkTrack *)Track->At(i))->ygood[ip];
528 l2->xm_nt[i][ip] = ((TrkTrack *)Track->At(i))->xm[ip];
529 l2->ym_nt[i][ip] = ((TrkTrack *)Track->At(i))->ym[ip];
530 l2->zm_nt[i][ip] = ((TrkTrack *)Track->At(i))->zm[ip];
531 l2->resx_nt[i][ip] = ((TrkTrack *)Track->At(i))->resx[ip];
532 l2->resy_nt[i][ip] = ((TrkTrack *)Track->At(i))->resy[ip];
533 l2->xv_nt[i][ip] = ((TrkTrack *)Track->At(i))->xv[ip];
534 l2->yv_nt[i][ip] = ((TrkTrack *)Track->At(i))->yv[ip];
535 l2->zv_nt[i][ip] = ((TrkTrack *)Track->At(i))->zv[ip];
536 l2->axv_nt[i][ip] = ((TrkTrack *)Track->At(i))->axv[ip];
537 l2->ayv_nt[i][ip] = ((TrkTrack *)Track->At(i))->ayv[ip];
538 l2->dedx_x[i][ip] = ((TrkTrack *)Track->At(i))->dedx_x[ip];
539 l2->dedx_y[i][ip] = ((TrkTrack *)Track->At(i))->dedx_y[ip];
540 };
541 }
542
543 // *** SINGLETS ***
544 l2->nclsx = SingletX->GetEntries();
545 for(Int_t i=0;i<l2->nclsx;i++){
546 l2->planex[i] = ((TrkSinglet *)SingletX->At(i))->plane;
547 l2->xs[i][0] = ((TrkSinglet *)SingletX->At(i))->coord[0];
548 l2->xs[i][1] = ((TrkSinglet *)SingletX->At(i))->coord[1];
549 l2->signlxs[i] = ((TrkSinglet *)SingletX->At(i))->sgnl;
550 }
551 l2->nclsy = SingletY->GetEntries();
552 for(Int_t i=0;i<l2->nclsy;i++){
553 l2->planey[i] = ((TrkSinglet *)SingletY->At(i))->plane;
554 l2->ys[i][0] = ((TrkSinglet *)SingletY->At(i))->coord[0];
555 l2->ys[i][1] = ((TrkSinglet *)SingletY->At(i))->coord[1];
556 l2->signlys[i] = ((TrkSinglet *)SingletY->At(i))->sgnl;
557 }
558 }
559 //--------------------------------------
560 //
561 //
562 //--------------------------------------
563 void TrkLevel2::Clear(){
564 // good2 = -1;
565 for(Int_t i=0; i<12 ; i++){
566 // crc[i] = -1;
567 good[i] = -1;
568 };
569 /* Track->RemoveAll();
570 SingletX->RemoveAll();
571 SingletY->RemoveAll();*/
572 // modify to avoid memory leakage
573 Track->Clear();
574 SingletX->Clear();
575 SingletY->Clear();
576 }
577 //--------------------------------------
578 //
579 //
580 //--------------------------------------
581 void TrkLevel2::Delete(){
582
583 Clear();
584 Track->Delete();
585 SingletX->Delete();
586 SingletY->Delete();
587 }
588 //--------------------------------------
589 //
590 //
591 //--------------------------------------
592 /**
593 * Sort physical tracks and stores them in a TObjectArray, ordering by increasing chi**2 value (in case of track image, it selects the one with lower chi**2). The total number of physical tracks is given by GetNTracks() and the it-th physical track can be retrieved by means of the method GetTrack(int it).
594 * This method is overridden by PamLevel2::GetTracks(), where calorimeter and TOF information is used.
595 */
596 TRefArray *TrkLevel2::GetTracks_NFitSorted(){
597
598 TRefArray *sorted = new TRefArray();
599
600 TClonesArray &t = *Track;
601 // TClonesArray &ts = *PhysicalTrack;
602 int N = ntrk();
603 vector<int> m(N); for(int i=0; i<N; i++)m[i]=1;
604 // int m[50]; for(int i=0; i<N; i++)m[i]=1;
605
606 int indo=0;
607 int indi=0;
608 while(N != 0){
609 int nfit =0;
610 float chi2ref = numeric_limits<float>::max();
611
612 // first loop to search maximum num. of fit points
613 for(int i=0; i < ntrk(); i++){
614 if( ((TrkTrack *)t[i])->GetNtot() >= nfit && m[i]==1){
615 nfit = ((TrkTrack *)t[i])->GetNtot();
616 }
617 }
618 //second loop to search minimum chi2 among selected
619 for(int i=0; i<this->ntrk(); i++){
620 Float_t chi2 = ((TrkTrack *)t[i])->chi2;
621 if(chi2 < 0) chi2 = chi2*1000;
622 if( chi2 < chi2ref
623 && ((TrkTrack *)t[i])->GetNtot() == nfit
624 && m[i]==1){
625 chi2ref = ((TrkTrack *)t[i])->chi2;
626 indi = i;
627 };
628 };
629 if( ((TrkTrack *)t[indi])->HasImage() ){
630 m[((TrkTrack *)t[indi])->image] = 0;
631 N--;
632
633 // cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl;
634 };
635 sorted->Add( (TrkTrack*)t[indi] );
636
637 m[indi] = 0;
638 // cout << "SORTED "<< indo << " "<< indi << " "<< N << endl;
639 N--;
640 indo++;
641 }
642 m.clear();
643 // cout << "GetTracks_NFitSorted(it): Done"<< endl;
644
645 return sorted;
646 // return PhysicalTrack;
647 }
648 //--------------------------------------
649 //
650 //
651 //--------------------------------------
652 /**
653 * Retrieves the is-th stored track.
654 * @param it Track number, ranging from 0 to ntrk().
655 * Fitted tracks ( images included ) are stored in a TObjectArray ( TrkLevel2::Track ) in the same order they are returned by the F77 fitting routine.
656 */
657 TrkTrack *TrkLevel2::GetStoredTrack(int is){
658
659 if(is >= this->ntrk()){
660 cout << "** TrkLevel2 ** Track "<< is << "doen not exits! " << endl;
661 cout << " Stored tracks ntrk() = "<< this->ntrk() << endl;
662 return 0;
663 }
664 TClonesArray &t = *(Track);
665 TrkTrack *track = (TrkTrack*)t[is];
666 return track;
667 }
668 //--------------------------------------
669 //
670 //
671 //--------------------------------------
672 /**
673 * Retrieves the is-th stored X singlet.
674 * @param it Singlet number, ranging from 0 to nclsx().
675 */
676 TrkSinglet *TrkLevel2::GetSingletX(int is){
677
678 if(is >= this->nclsx()){
679 cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl;
680 cout << " Stored x-singlets nclsx() = "<< this->nclsx() << endl;
681 return 0;
682 }
683 TClonesArray &t = *(SingletX);
684 TrkSinglet *singlet = (TrkSinglet*)t[is];
685 return singlet;
686 }
687 //--------------------------------------
688 //
689 //
690 //--------------------------------------
691 /**
692 * Retrieves the is-th stored Y singlet.
693 * @param it Singlet number, ranging from 0 to nclsx().
694 */
695 TrkSinglet *TrkLevel2::GetSingletY(int is){
696
697 if(is >= this->nclsy()){
698 cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl;
699 cout << " Stored y-singlets nclsy() = "<< this->nclsx() << endl;
700 return 0;
701 }
702 TClonesArray &t = *(SingletY);
703 TrkSinglet *singlet = (TrkSinglet*)t[is];
704 return singlet;
705 }
706 //--------------------------------------
707 //
708 //
709 //--------------------------------------
710 /**
711 * Retrieves the it-th "physical" track, sorted by the method GetNTracks().
712 * @param it Track number, ranging from 0 to GetNTracks().
713 */
714
715 TrkTrack *TrkLevel2::GetTrack(int it){
716
717 if(it >= this->GetNTracks()){
718 cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl;
719 cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl;
720 return 0;
721 }
722
723 TRefArray *sorted = GetTracks(); //TEMPORANEO
724 TrkTrack *track = (TrkTrack*)sorted->At(it);
725 sorted->Delete();
726 return track;
727 }
728 /**
729 * Give the number of "physical" tracks, sorted by the method GetTracks().
730 */
731 Int_t TrkLevel2::GetNTracks(){
732
733 Float_t ntot=0;
734 TClonesArray &t = *Track;
735 for(int i=0; i<ntrk(); i++) {
736 if( ((TrkTrack *)t[i])->GetImageSeqNo() == -1 ) ntot+=1.;
737 else ntot+=0.5;
738 }
739 return (Int_t)ntot;
740
741 };
742 //--------------------------------------
743 //
744 //
745 //--------------------------------------
746 /**
747 * Retrieves (if present) the image of the it-th "physical" track, sorted by the method GetNTracks().
748 * @param it Track number, ranging from 0 to GetNTracks().
749 */
750 TrkTrack *TrkLevel2::GetTrackImage(int it){
751
752 if(it >= this->GetNTracks()){
753 cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl;
754 cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl;
755 return 0;
756 }
757
758 TRefArray* sorted = GetTracks(); //TEMPORANEO
759 TrkTrack *track = (TrkTrack*)sorted->At(it);
760
761 if(!track->HasImage()){
762 cout << "** TrkLevel2 ** Track "<< it << "does not have image! " << endl;
763 return 0;
764 }
765 TrkTrack *image = (TrkTrack*)(*Track)[track->image];
766
767 sorted->Delete();
768
769 return image;
770
771 }
772 //--------------------------------------
773 //
774 //
775 //--------------------------------------
776 /**
777 * Loads the magnetic field.
778 * @param s Path of the magnetic-field files.
779 */
780 void TrkLevel2::LoadField(TString s){
781 readb_(s.Data());
782 };
783 //--------------------------------------
784 //
785 //
786 //--------------------------------------
787 /**
788 * Get tracker-plane (mechanical) z-coordinate
789 * @param plane_id plane index (1=TOP,2,3,4,5,6=BOTTOM)
790 */
791 Float_t TrkLevel2::GetZTrk(Int_t plane_id){
792 switch(plane_id){
793 case 1: return ZTRK1;
794 case 2: return ZTRK2;
795 case 3: return ZTRK3;
796 case 4: return ZTRK4;
797 case 5: return ZTRK5;
798 case 6: return ZTRK6;
799 default: return 0.;
800 };
801 };
802 //--------------------------------------
803 //
804 //
805 //--------------------------------------
806 /**
807 * Trajectory default constructor.
808 * (By default is created with z-coordinates inside the tracking volume)
809 */
810 Trajectory::Trajectory(){
811 npoint = 10;
812 x = new float[npoint];
813 y = new float[npoint];
814 z = new float[npoint];
815 thx = new float[npoint];
816 thy = new float[npoint];
817 tl = new float[npoint];
818 float dz = ((ZTRK1)-(ZTRK6))/(npoint-1);
819 for(int i=0; i<npoint; i++){
820 x[i] = 0;
821 y[i] = 0;
822 z[i] = (ZTRK1) - i*dz;
823 thx[i] = 0;
824 thy[i] = 0;
825 tl[i] = 0;
826 }
827 }
828 //--------------------------------------
829 //
830 //
831 //--------------------------------------
832 /**
833 * Trajectory constructor.
834 * (By default is created with z-coordinates inside the tracking volume)
835 * \param n Number of points
836 */
837 Trajectory::Trajectory(int n){
838 if(n<=0){
839 cout << "NB! Trajectory must have at least 1 point >>> created with 10 points" << endl;
840 n=10;
841 }
842 npoint = n;
843 x = new float[npoint];
844 y = new float[npoint];
845 z = new float[npoint];
846 thx = new float[npoint];
847 thy = new float[npoint];
848 tl = new float[npoint];
849 float dz = ((ZTRK1)-(ZTRK6))/(npoint-1);
850 for(int i=0; i<npoint; i++){
851 x[i] = 0;
852 y[i] = 0;
853 z[i] = (ZTRK1) - i*dz;
854 thx[i] = 0;
855 thy[i] = 0;
856 tl[i] = 0;
857 }
858 }
859 //--------------------------------------
860 //
861 //
862 //--------------------------------------
863 /**
864 * Trajectory constructor.
865 * \param n Number of points
866 * \param pz Pointer to float array, defining z coordinates
867 */
868 Trajectory::Trajectory(int n, float* zin){
869 npoint = 10;
870 if(n>0)npoint = n;
871 x = new float[npoint];
872 y = new float[npoint];
873 z = new float[npoint];
874 thx = new float[npoint];
875 thy = new float[npoint];
876 tl = new float[npoint];
877 int i=0;
878 do{
879 x[i] = 0;
880 y[i] = 0;
881 z[i] = zin[i];
882 thx[i] = 0;
883 thy[i] = 0;
884 tl[i] = 0;
885 i++;
886 }while(zin[i-1] > zin[i] && i < npoint);
887 npoint=i;
888 if(npoint != n)cout << "NB! Trajectory created with "<<npoint<<" points"<<endl;
889 }
890 //--------------------------------------
891 //
892 //
893 //--------------------------------------
894 /**
895 * Dump the trajectory coordinates.
896 */
897 void Trajectory::Dump(){
898 cout <<endl<< "Trajectory ========== "<<endl;
899 for (int i=0; i<npoint; i++){
900 cout << i <<" >> " << x[i] <<" "<< y[i] <<" "<< z[i] ;
901 cout <<" -- " << thx[i] <<" "<< thy[i] ;
902 cout <<" -- " << tl[i] << endl;
903 };
904 }
905 //--------------------------------------
906 //
907 //
908 //--------------------------------------
909 /**
910 * Get trajectory length between two points
911 * @param ifirst first point (default 0)
912 * @param ilast last point (default npoint)
913 */
914 float Trajectory::GetLength(int ifirst, int ilast){
915 if( ifirst<0 ) ifirst = 0;
916 if( ilast>=npoint) ilast = npoint-1;
917 float l=0;
918 for(int i=ifirst;i<=ilast;i++){
919 l=l+tl[i];
920 };
921 if(z[ilast] > ZINI)l=l-tl[ilast];
922 if(z[ifirst] < ZINI) l=l-tl[ifirst];
923
924 return l;
925
926 }
927
928
929 ClassImp(TrkLevel2);
930 ClassImp(TrkSinglet);
931 ClassImp(TrkTrack);
932 ClassImp(Trajectory);

  ViewVC Help
Powered by ViewVC 1.1.23