/[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.2 - (show annotations) (download)
Thu Jun 1 09:03:09 2006 UTC (18 years, 6 months ago) by pam-fi
Branch: MAIN
Changes since 1.1: +121 -6 lines
Tracking routine updated in order to get track projected angles and track length

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 image = 0;
22 chi2 = 0;
23 for(int it1=0;it1<5;it1++){
24 al[it1] = 0;
25 for(int it2=0;it2<5;it2++)
26 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 };
44 };
45 //--------------------------------------
46 //
47 //
48 //--------------------------------------
49 TrkTrack::TrkTrack(const TrkTrack& t){
50 image = t.image;
51 chi2 = t.chi2;
52 for(int it1=0;it1<5;it1++){
53 al[it1] = t.al[it1];
54 for(int it2=0;it2<5;it2++)
55 coval[it1][it2] = t.coval[it1][it2];
56 };
57 for(int ip=0;ip<6;ip++){
58 xgood[ip] = t.xgood[ip];
59 ygood[ip] = t.ygood[ip];
60 xm[ip] = t.xm[ip];
61 ym[ip] = t.ym[ip];
62 zm[ip] = t.zm[ip];
63 resx[ip] = t.resx[ip];
64 resy[ip] = t.resy[ip];
65 xv[ip] = t.xv[ip];
66 yv[ip] = t.yv[ip];
67 zv[ip] = t.zv[ip];
68 axv[ip] = t.axv[ip];
69 ayv[ip] = t.ayv[ip];
70 dedx_x[ip] = t.dedx_x[ip];
71 dedx_y[ip] = t.dedx_y[ip];
72 };
73 };
74 //--------------------------------------
75 //
76 //
77 //--------------------------------------
78 /**
79 * Evaluates the trajectory in the apparatus associated to the track.
80 * 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.
81 * @param t pointer to an object of the class Trajectory,
82 * which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ).
83 * @return error flag.
84 */
85 int TrkTrack::DoTrack(Trajectory* t){
86
87 double *dxout = new double[t->npoint];
88 double *dyout = new double[t->npoint];
89 double *dzin = new double[t->npoint];
90 double dal[5];
91
92 int ifail = 0;
93
94 for (int i=0; i<5; i++) dal[i] = (double)al[i];
95 for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i];
96
97 dotrack_(&(t->npoint),dzin,dxout,dyout,dal,&ifail);
98
99 for (int i=0; i<t->npoint; i++){
100 t->x[i] = (float)*dxout++;
101 t->y[i] = (float)*dyout++;
102 }
103
104 // delete [] dxout;
105 // delete [] dyout;
106 // delete [] dzin;
107
108 return ifail;
109 };
110 //--------------------------------------
111 //
112 //
113 //--------------------------------------
114 /**
115 * Evaluates the trajectory in the apparatus associated to the track.
116 * 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.
117 * @param t pointer to an object of the class Trajectory,
118 * which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ).
119 * @return error flag.
120 */
121 int TrkTrack::DoTrack2(Trajectory* t){
122
123 double *dxout = new double[t->npoint];
124 double *dyout = new double[t->npoint];
125 double *dthxout = new double[t->npoint];
126 double *dthyout = new double[t->npoint];
127 double *dtlout = new double[t->npoint];
128 double *dzin = new double[t->npoint];
129 double dal[5];
130
131 int ifail = 0;
132
133 for (int i=0; i<5; i++) dal[i] = (double)al[i];
134 for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i];
135
136 dotrack2_(&(t->npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail);
137
138 for (int i=0; i<t->npoint; i++){
139 t->x[i] = (float)*dxout++;
140 t->y[i] = (float)*dyout++;
141 t->thx[i] = (float)*dthxout++;
142 t->thy[i] = (float)*dthyout++;
143 t->tl[i] = (float)*dtlout++;
144 }
145
146 // delete [] dxout;
147 // delete [] dyout;
148 // delete [] dzin;
149
150 return ifail;
151 };
152 //--------------------------------------
153 //
154 //
155 //--------------------------------------
156 //float TrkTrack::BdL(){
157 //};
158 //--------------------------------------
159 //
160 //
161 //--------------------------------------
162 Float_t TrkTrack::GetRigidity(){
163 Float_t rig=0;
164 if(chi2>0)rig=1./al[4];
165 if(rig<0) rig=-rig;
166 return rig;
167 };
168 //
169 Float_t TrkTrack::GetDeflection(){
170 Float_t def=0;
171 if(chi2>0)def=al[4];
172 return def;
173 };
174 //
175 Float_t TrkTrack::GetDEDX(){
176 Float_t dedx=0;
177 for(Int_t i=0; i<6; i++)dedx+=dedx_x[i]*xgood[i]+dedx_y[i]*ygood[i];
178 dedx = dedx/(this->GetNX()+this->GetNY());
179 return dedx;
180 };
181
182 //--------------------------------------
183 //
184 //
185 //--------------------------------------
186 void TrkTrack::Dump(){
187 cout << endl << "========== Track " ;
188 cout << endl << "al : "; for(int i=0; i<5; i++)cout << al[i] << " ";
189 cout << endl << "chi^2 : "<< chi2;
190 cout << endl << "xgood : "; for(int i=0; i<6; i++)cout << xgood[i] ;
191 cout << endl << "ygood : "; for(int i=0; i<6; i++)cout << ygood[i] ;
192 cout << endl << "xm : "; for(int i=0; i<6; i++)cout << xm[i] << " ";
193 cout << endl << "ym : "; for(int i=0; i<6; i++)cout << ym[i] << " ";
194 cout << endl << "zm : "; for(int i=0; i<6; i++)cout << zm[i] << " ";
195 cout << endl << "dedx_x : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " ";
196 cout << endl << "dedx_y : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " ";
197 }
198 //--------------------------------------
199 //
200 //
201 //--------------------------------------
202 TrkSinglet::TrkSinglet(){
203 plane = 0;
204 coord[0] = 0;
205 coord[1] = 0;
206 sgnl = 0;
207 };
208 //--------------------------------------
209 //
210 //
211 //--------------------------------------
212 TrkSinglet::TrkSinglet(const TrkSinglet& s){
213 plane = s.plane;
214 coord[0] = s.coord[0];
215 coord[1] = s.coord[1];
216 sgnl = s.sgnl;
217 };
218 //--------------------------------------
219 //
220 //
221 //--------------------------------------
222 void TrkSinglet::Dump(){
223 int i=0;
224 cout << endl << "========== Singlet " ;
225 cout << endl << "plane : " << plane;
226 cout << endl << "coord[2] : "; while( i<2 && cout << coord[i] << " ") i++;
227 cout << endl << "sgnl : " << sgnl;
228 }
229 //--------------------------------------
230 //
231 //
232 //--------------------------------------
233 TrkLevel2::TrkLevel2(){
234 good2 = -1;
235 for(Int_t i=0; i<12 ; i++){
236 crc[i] = -1;
237 };
238 Track = new TClonesArray("TrkTrack");
239 SingletX = new TClonesArray("TrkSinglet");
240 SingletY = new TClonesArray("TrkSinglet");
241 }
242 //--------------------------------------
243 //
244 //
245 //--------------------------------------
246 void TrkLevel2::Dump(){
247 TClonesArray &t = *Track;
248 TClonesArray &sx = *SingletX;
249 TClonesArray &sy = *SingletY;
250
251 cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-";
252 cout << endl << "good2 : " << good2;
253 cout << endl << "crc : "; for(int i=0; i<12; i++) cout << crc[i];
254 cout << endl << "ntrk() : " << this->ntrk() ;
255 cout << endl << "nclsx() : " << this->nclsx();
256 cout << endl << "nclsy() : " << this->nclsy();
257 for(int i=0; i<this->ntrk(); i++) ((TrkTrack *)t[i])->Dump();
258 for(int i=0; i<this->nclsx(); i++) ((TrkSinglet *)sx[i])->Dump();
259 for(int i=0; i<this->nclsy(); i++) ((TrkSinglet *)sy[i])->Dump();
260 }
261 //--------------------------------------
262 //
263 //
264 //--------------------------------------
265 /**
266 * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common).
267 */
268 void TrkLevel2::FillCommonVar(cTrkLevel2 *l2){
269 // temporary objects:
270 TrkSinglet* t_singlet = new TrkSinglet();
271 TrkTrack* t_track = new TrkTrack();
272 // general variables
273 good2 = l2->good2;
274 for(Int_t i=0; i<12 ; i++){
275 crc[i] = l2->crc[i];
276 };
277 // *** TRACKS ***
278 TClonesArray &t = *Track;
279 for(int i=0; i<l2->ntrk; i++){
280 t_track->image = l2->image[i]-1;
281 t_track->chi2 = l2->chi2_nt[i];
282 for(int it1=0;it1<5;it1++){
283 t_track->al[it1] = l2->al_nt[i][it1];
284 for(int it2=0;it2<5;it2++)
285 t_track->coval[it1][it2] = l2->coval[i][it2][it1];
286 };
287 for(int ip=0;ip<6;ip++){
288 t_track->xgood[ip] = l2->xgood_nt[i][ip];
289 t_track->ygood[ip] = l2->ygood_nt[i][ip];
290 t_track->xm[ip] = l2->xm_nt[i][ip];
291 t_track->ym[ip] = l2->ym_nt[i][ip];
292 t_track->zm[ip] = l2->zm_nt[i][ip];
293 t_track->resx[ip] = l2->resx_nt[i][ip];
294 t_track->resy[ip] = l2->resy_nt[i][ip];
295 t_track->xv[ip] = l2->xv_nt[i][ip];
296 t_track->yv[ip] = l2->yv_nt[i][ip];
297 t_track->zv[ip] = l2->zv_nt[i][ip];
298 t_track->axv[ip] = l2->axv_nt[i][ip];
299 t_track->ayv[ip] = l2->ayv_nt[i][ip];
300 t_track->dedx_x[ip] = l2->dedx_x[i][ip];
301 t_track->dedx_y[ip] = l2->dedx_y[i][ip];
302 };
303 new(t[i]) TrkTrack(*t_track);
304 t_track->Clear();
305 };
306 // *** SINGLETS ***
307 TClonesArray &sx = *SingletX;
308 for(int i=0; i<l2->nclsx; i++){
309 t_singlet->plane = l2->planex[i];
310 t_singlet->coord[0] = l2->xs[i][0];
311 t_singlet->coord[1] = l2->xs[i][1];
312 t_singlet->sgnl = l2->signlxs[i];
313 new(sx[i]) TrkSinglet(*t_singlet);
314 t_singlet->Clear();
315 }
316 TClonesArray &sy = *SingletY;
317 for(int i=0; i<l2->nclsy; i++){
318 t_singlet->plane = l2->planey[i];
319 t_singlet->coord[0] = l2->ys[i][0];
320 t_singlet->coord[1] = l2->ys[i][1];
321 t_singlet->sgnl = l2->signlys[i];
322 new(sy[i]) TrkSinglet(*t_singlet);
323 t_singlet->Clear();
324 };
325 }
326 //--------------------------------------
327 //
328 //
329 //--------------------------------------
330 void TrkLevel2::Clear(){
331 good2 = -1;
332 for(Int_t i=0; i<12 ; i++){
333 crc[i] = -1;
334 };
335 Track->RemoveAll();
336 SingletX->RemoveAll();
337 SingletY->RemoveAll();
338 }
339 //--------------------------------------
340 //
341 //
342 //--------------------------------------
343 /**
344 * 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).
345 * This method is overridden by PamLevel2::GetTracks(), where calorimeter and TOF information is used.
346 */
347 TClonesArray *TrkLevel2::GetTracks(){
348
349 TClonesArray *sorted = new TClonesArray("TrkTrack");
350 TClonesArray &t = *Track;
351 TClonesArray &ts = *sorted;
352 int N=this->ntrk();
353 vector<int> m(N); for(int i=0; i<N; i++)m[i]=1;
354
355 int indo=0;
356 int indi=0;
357 while(N != 0){
358 float chi2ref=1000000;
359 for(int i=0; i<this->ntrk(); i++){
360 if(((TrkTrack *)t[i])->chi2 < chi2ref && m[i]==1){
361 chi2ref = ((TrkTrack *)t[i])->chi2;
362 indi = i;
363 }
364 }
365 if( ((TrkTrack *)t[indi])->image != -1 ){
366 m[((TrkTrack *)t[indi])->image] = 0;
367 N--;
368 }
369 new(ts[indo]) TrkTrack(*(TrkTrack*)t[indi]);
370 m[indi] = 0;
371 N--;
372 indo++;
373 }
374 return sorted;
375 }
376 //--------------------------------------
377 //
378 //
379 //--------------------------------------
380 /**
381 * Retrieves the is-th stored track.
382 * @param it Track number, ranging from 0 to ntrk().
383 * Fitted tracks ( images included ) are stored in a TObjectArray ( TrkLevel2::Track ) in the same order they are returned by the F77 fitting routine.
384 */
385 TrkTrack *TrkLevel2::GetStoredTrack(int is){
386
387 if(is >= this->ntrk()){
388 cout << "** TrkLevel2 ** Track "<< is << "doen not exits! " << endl;
389 cout << " Stored tracks ntrk() = "<< this->ntrk() << endl;
390 return 0;
391 }
392 TClonesArray &t = *(Track);
393 TrkTrack *track = (TrkTrack*)t[is];
394 return track;
395 }
396 //--------------------------------------
397 //
398 //
399 //--------------------------------------
400 /**
401 * Retrieves the it-th "physical" track, sorted by the method GetNTracks().
402 * @param it Track number, ranging from 0 to GetNTracks().
403 */
404 TrkTrack *TrkLevel2::GetTrack(int it){
405
406 if(it >= this->GetNTracks()){
407 cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl;
408 cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl;
409 return 0;
410 }
411 TrkTrack *track = (TrkTrack*)(*(this->GetTracks()))[it];
412 return track;
413 }
414 //--------------------------------------
415 //
416 //
417 //--------------------------------------
418 /**
419 * Retrieves (if present) the image of the it-th "physical" track, sorted by the method GetNTracks().
420 * @param it Track number, ranging from 0 to GetNTracks().
421 */
422 TrkTrack *TrkLevel2::GetTrackImage(int it){
423
424 if(it >= this->GetNTracks()){
425 cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl;
426 cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl;
427 return 0;
428 }
429 TrkTrack *track = (TrkTrack*)(*(this->GetTracks()))[it];
430 if(!track->HasImage()){
431 cout << "** TrkLevel2 ** Track "<< it << "does not have image! " << endl;
432 return 0;
433 }
434 TrkTrack *image = (TrkTrack*)(*Track)[track->image];
435 return image;
436
437 }
438 //--------------------------------------
439 //
440 //
441 //--------------------------------------
442 /**
443 * Loads the magnetic field.
444 * @param s Path of the magnetic-field files.
445 */
446 void TrkLevel2::LoadField(TString s){
447 readb_(s.Data());
448 };
449 //--------------------------------------
450 //
451 //
452 //--------------------------------------
453 /**
454 * Trajectory default constructor.
455 * (By default is created with z-coordinates inside the tracking volume)
456 */
457 Trajectory::Trajectory(){
458 npoint = 10;
459 x = new float[npoint];
460 y = new float[npoint];
461 z = new float[npoint];
462 thx = new float[npoint];
463 thy = new float[npoint];
464 tl = new float[npoint];
465 float dz = ((ZTRKUP)-(ZTRKDW))/(npoint-1);
466 for(int i=0; i<npoint; i++){
467 x[i] = 0;
468 y[i] = 0;
469 z[i] = (ZTRKUP) - i*dz;
470 thx[i] = 0;
471 thy[i] = 0;
472 tl[i] = 0;
473 }
474 }
475 //--------------------------------------
476 //
477 //
478 //--------------------------------------
479 /**
480 * Trajectory constructor.
481 * (By default is created with z-coordinates inside the tracking volume)
482 * \param n Number of points
483 */
484 Trajectory::Trajectory(int n){
485 if(n<=0){
486 cout << "NB! Trajectory must have at least 1 point >>> created with 10 points" << endl;
487 n=10;
488 }
489 npoint = n;
490 x = new float[npoint];
491 y = new float[npoint];
492 z = new float[npoint];
493 thx = new float[npoint];
494 thy = new float[npoint];
495 tl = new float[npoint];
496 float dz = ((ZTRKUP)-(ZTRKDW))/(npoint-1);
497 for(int i=0; i<npoint; i++){
498 x[i] = 0;
499 y[i] = 0;
500 z[i] = (ZTRKUP) - i*dz;
501 thx[i] = 0;
502 thy[i] = 0;
503 tl[i] = 0;
504 }
505 }
506 //--------------------------------------
507 //
508 //
509 //--------------------------------------
510 /**
511 * Trajectory constructor.
512 * \param n Number of points
513 * \param pz Pointer to float array, defining z coordinates
514 */
515 Trajectory::Trajectory(int n, float* zin){
516 npoint = 10;
517 if(n>0)npoint = n;
518 x = new float[npoint];
519 y = new float[npoint];
520 z = new float[npoint];
521 thx = new float[npoint];
522 thy = new float[npoint];
523 tl = new float[npoint];
524 int i=0;
525 do{
526 x[i] = 0;
527 y[i] = 0;
528 z[i] = zin[i];
529 thx[i] = 0;
530 thy[i] = 0;
531 tl[i] = 0;
532 i++;
533 }while(zin[i-1] > zin[i] && i < npoint);
534 npoint=i;
535 if(npoint != n)cout << "NB! Trajectory created with "<<npoint<<" points"<<endl;
536 }
537 //--------------------------------------
538 //
539 //
540 //--------------------------------------
541 /**
542 * Dump the trajectory coordinates.
543 */
544 void Trajectory::Dump(){
545 cout <<endl<< "Trajectory ========== "<<endl;
546 for (int i=0; i<npoint; i++){
547 cout << i <<" >> " << x[i] <<" "<< y[i] <<" "<< z[i] ;
548 cout <<" -- " << thx[i] <<" "<< thy[i] ;
549 cout <<" -- " << tl[i] << endl;
550 };
551 }
552 //--------------------------------------
553 //
554 //
555 //--------------------------------------
556 /**
557 * Get trajectory length between two points
558 * @param ifirst first point (default 0)
559 * @param ilast last point (default npoint)
560 */
561 float Trajectory::GetLength(int ifirst, int ilast){
562 if( ifirst<0 ) ifirst = 0;
563 if( ilast>=npoint) ilast = npoint-1;
564 float l=0;
565 for(int i=ifirst;i<=ilast;i++){
566 l=l+tl[i];
567 };
568 if(z[ilast] > ZINI)l=l-tl[ilast];
569 if(z[ifirst] < ZINI) l=l-tl[ifirst];
570
571 return l;
572
573 }
574 ClassImp(TrkLevel2);
575 ClassImp(TrkSinglet);
576 ClassImp(TrkTrack);
577 ClassImp(Trajectory);

  ViewVC Help
Powered by ViewVC 1.1.23