/[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.9 - (show annotations) (download)
Tue Sep 5 12:52:20 2006 UTC (18 years, 3 months ago) by pam-fi
Branch: MAIN
CVS Tags: v2r00BETA
Changes since 1.8: +192 -90 lines
implemented class TrkLevel1

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

  ViewVC Help
Powered by ViewVC 1.1.23