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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.31 - (hide annotations) (download)
Wed Mar 28 09:22:28 2007 UTC (17 years, 8 months ago) by pam-fi
Branch: MAIN
Changes since 1.30: +17 -0 lines
methods to set FACT and ISTEPMIN for tracking precision

1 mocchiut 1.1 /**
2     * \file TrkLevel2.cpp
3     * \author Elena Vannuccini
4     */
5     #include <TrkLevel2.h>
6     #include <iostream>
7 pam-fi 1.14 #include <math.h>
8 mocchiut 1.1 using namespace std;
9     //......................................
10     // F77 routines
11     //......................................
12     extern "C" {
13     void dotrack_(int*, double*, double*, double*, double*, int*);
14 pam-fi 1.2 void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);
15 pam-fi 1.16 // int readb_(const char*);
16 pam-fi 1.26 // int readb_();
17     void mini2_(int*,int*,int*);
18     void guess_();
19     void gufld_(float*, float*);
20 mocchiut 1.1 }
21 pam-fi 1.26
22 mocchiut 1.1 //--------------------------------------
23     //
24     //
25     //--------------------------------------
26     TrkTrack::TrkTrack(){
27 pam-fi 1.21 // cout << "TrkTrack::TrkTrack()" << endl;
28 pam-fi 1.3 seqno = -1;
29     image = -1;
30 mocchiut 1.1 chi2 = 0;
31 pam-fi 1.14 nstep = 0;
32     for(int it1=0;it1<5;it1++){
33     al[it1] = 0;
34     for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;
35 mocchiut 1.1 };
36     for(int ip=0;ip<6;ip++){
37 pam-fi 1.21 xgood[ip] = 0;
38     ygood[ip] = 0;
39     xm[ip] = 0;
40     ym[ip] = 0;
41     zm[ip] = 0;
42     resx[ip] = 0;
43     resy[ip] = 0;
44     xv[ip] = 0;
45     yv[ip] = 0;
46     zv[ip] = 0;
47     axv[ip] = 0;
48     ayv[ip] = 0;
49     dedx_x[ip] = 0;
50     dedx_y[ip] = 0;
51     };
52     clx = 0;
53     cly = 0;
54 pam-fi 1.29 // clx = new TRefArray(6,0); //forse causa memory leak???
55     // cly = new TRefArray(6,0); //forse causa memory leak???
56 mocchiut 1.1 };
57     //--------------------------------------
58     //
59     //
60     //--------------------------------------
61     TrkTrack::TrkTrack(const TrkTrack& t){
62 pam-fi 1.3 seqno = t.seqno;
63 mocchiut 1.1 image = t.image;
64     chi2 = t.chi2;
65 pam-fi 1.14 nstep = t.nstep;
66     for(int it1=0;it1<5;it1++){
67     al[it1] = t.al[it1];
68     for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
69 mocchiut 1.1 };
70     for(int ip=0;ip<6;ip++){
71 pam-fi 1.21 xgood[ip] = t.xgood[ip];
72     ygood[ip] = t.ygood[ip];
73     xm[ip] = t.xm[ip];
74     ym[ip] = t.ym[ip];
75     zm[ip] = t.zm[ip];
76     resx[ip] = t.resx[ip];
77     resy[ip] = t.resy[ip];
78     xv[ip] = t.xv[ip];
79     yv[ip] = t.yv[ip];
80     zv[ip] = t.zv[ip];
81     axv[ip] = t.axv[ip];
82     ayv[ip] = t.ayv[ip];
83     dedx_x[ip] = t.dedx_x[ip];
84     dedx_y[ip] = t.dedx_y[ip];
85     };
86 pam-fi 1.22 clx = 0;
87     cly = 0;
88     if(t.clx)clx = new TRefArray(*(t.clx));
89     if(t.cly)cly = new TRefArray(*(t.cly));
90 pam-fi 1.9
91 mocchiut 1.1 };
92     //--------------------------------------
93     //
94     //
95     //--------------------------------------
96 pam-fi 1.21 void TrkTrack::Copy(TrkTrack& t){
97    
98     t.seqno = seqno;
99     t.image = image;
100     t.chi2 = chi2;
101     t.nstep = nstep;
102     for(int it1=0;it1<5;it1++){
103     t.al[it1] = al[it1];
104     for(int it2=0;it2<5;it2++)t.coval[it1][it2] = coval[it1][it2];
105     };
106     for(int ip=0;ip<6;ip++){
107     t.xgood[ip] = xgood[ip];
108     t.ygood[ip] = ygood[ip];
109     t.xm[ip] = xm[ip];
110     t.ym[ip] = ym[ip];
111     t.zm[ip] = zm[ip];
112     t.resx[ip] = resx[ip];
113     t.resy[ip] = resy[ip];
114     t.xv[ip] = xv[ip];
115     t.yv[ip] = yv[ip];
116     t.zv[ip] = zv[ip];
117     t.axv[ip] = axv[ip];
118     t.ayv[ip] = ayv[ip];
119     t.dedx_x[ip] = dedx_x[ip];
120     t.dedx_y[ip] = dedx_y[ip];
121    
122     };
123    
124     };
125     //--------------------------------------
126     //
127     //
128     //--------------------------------------
129 mocchiut 1.1 /**
130     * Evaluates the trajectory in the apparatus associated to the track.
131     * 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.
132     * @param t pointer to an object of the class Trajectory,
133     * which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ).
134     * @return error flag.
135     */
136     int TrkTrack::DoTrack(Trajectory* t){
137    
138     double *dxout = new double[t->npoint];
139     double *dyout = new double[t->npoint];
140     double *dzin = new double[t->npoint];
141     double dal[5];
142    
143     int ifail = 0;
144    
145     for (int i=0; i<5; i++) dal[i] = (double)al[i];
146     for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i];
147    
148 pam-fi 1.26 TrkParams::Load(1);
149     if( !TrkParams::IsLoaded(1) ){
150     cout << "int TrkTrack::DoTrack(Trajectory* t) --- ERROR --- m.field not loaded"<<endl;
151     return 0;
152     }
153 mocchiut 1.1 dotrack_(&(t->npoint),dzin,dxout,dyout,dal,&ifail);
154    
155     for (int i=0; i<t->npoint; i++){
156     t->x[i] = (float)*dxout++;
157     t->y[i] = (float)*dyout++;
158     }
159    
160     // delete [] dxout;
161     // delete [] dyout;
162     // delete [] dzin;
163    
164     return ifail;
165     };
166     //--------------------------------------
167     //
168     //
169     //--------------------------------------
170 pam-fi 1.2 /**
171     * Evaluates the trajectory in the apparatus associated to the track.
172     * 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.
173     * @param t pointer to an object of the class Trajectory,
174     * which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ).
175     * @return error flag.
176     */
177     int TrkTrack::DoTrack2(Trajectory* t){
178    
179     double *dxout = new double[t->npoint];
180     double *dyout = new double[t->npoint];
181     double *dthxout = new double[t->npoint];
182     double *dthyout = new double[t->npoint];
183     double *dtlout = new double[t->npoint];
184     double *dzin = new double[t->npoint];
185     double dal[5];
186    
187     int ifail = 0;
188    
189     for (int i=0; i<5; i++) dal[i] = (double)al[i];
190     for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i];
191    
192 pam-fi 1.26 TrkParams::Load(1);
193     if( !TrkParams::IsLoaded(1) ){
194     cout << "int TrkTrack::DoTrack2(Trajectory* t) --- ERROR --- m.field not loaded"<<endl;
195     return 0;
196     }
197 pam-fi 1.2 dotrack2_(&(t->npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail);
198    
199     for (int i=0; i<t->npoint; i++){
200     t->x[i] = (float)*dxout++;
201     t->y[i] = (float)*dyout++;
202     t->thx[i] = (float)*dthxout++;
203     t->thy[i] = (float)*dthyout++;
204     t->tl[i] = (float)*dtlout++;
205     }
206    
207     // delete [] dxout;
208     // delete [] dyout;
209     // delete [] dzin;
210    
211     return ifail;
212     };
213     //--------------------------------------
214     //
215     //
216     //--------------------------------------
217 mocchiut 1.1 //float TrkTrack::BdL(){
218     //};
219     //--------------------------------------
220     //
221     //
222     //--------------------------------------
223     Float_t TrkTrack::GetRigidity(){
224     Float_t rig=0;
225     if(chi2>0)rig=1./al[4];
226     if(rig<0) rig=-rig;
227     return rig;
228     };
229     //
230     Float_t TrkTrack::GetDeflection(){
231     Float_t def=0;
232     if(chi2>0)def=al[4];
233     return def;
234     };
235     //
236     Float_t TrkTrack::GetDEDX(){
237     Float_t dedx=0;
238 pam-fi 1.30 // for(Int_t i=0; i<6; i++)dedx+=dedx_x[i]*xgood[i]+dedx_y[i]*ygood[i];
239     for(Int_t i=0; i<6; i++)dedx+=dedx_x[i]*XGood(i)+dedx_y[i]*YGood(i);
240 mocchiut 1.1 dedx = dedx/(this->GetNX()+this->GetNY());
241     return dedx;
242     };
243    
244     //--------------------------------------
245     //
246     //
247     //--------------------------------------
248     void TrkTrack::Dump(){
249     cout << endl << "========== Track " ;
250 pam-fi 1.13 cout << endl << "seq. n. : "<< seqno;
251     cout << endl << "image n. : "<< image;
252     cout << endl << "al : "; for(int i=0; i<5; i++)cout << al[i] << " ";
253 mocchiut 1.1 cout << endl << "chi^2 : "<< chi2;
254 pam-fi 1.13 cout << endl << "n.step : "<< nstep;
255 pam-fi 1.30 cout << endl << "xgood : "; for(int i=0; i<6; i++)cout << XGood(i) ;
256     cout << endl << "ygood : "; for(int i=0; i<6; i++)cout << YGood(i) ;
257 mocchiut 1.1 cout << endl << "xm : "; for(int i=0; i<6; i++)cout << xm[i] << " ";
258     cout << endl << "ym : "; for(int i=0; i<6; i++)cout << ym[i] << " ";
259     cout << endl << "zm : "; for(int i=0; i<6; i++)cout << zm[i] << " ";
260 pam-fi 1.13 cout << endl << "xv : "; for(int i=0; i<6; i++)cout << xv[i] << " ";
261     cout << endl << "yv : "; for(int i=0; i<6; i++)cout << yv[i] << " ";
262     cout << endl << "zv : "; for(int i=0; i<6; i++)cout << zv[i] << " ";
263     cout << endl << "resx : "; for(int i=0; i<6; i++)cout << resx[i] << " ";
264     cout << endl << "resy : "; for(int i=0; i<6; i++)cout << resy[i] << " ";
265     cout << endl << "coval : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" ";
266     cout << endl << " "; for(int i=0; i<5; i++)cout << coval[1][i]<<" ";
267     cout << endl << " "; for(int i=0; i<5; i++)cout << coval[2][i]<<" ";
268     cout << endl << " "; for(int i=0; i<5; i++)cout << coval[3][i]<<" ";
269     cout << endl << " "; for(int i=0; i<5; i++)cout << coval[4][i]<<" ";
270 mocchiut 1.1 cout << endl << "dedx_x : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " ";
271     cout << endl << "dedx_y : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " ";
272 pam-fi 1.14 cout << endl;
273 mocchiut 1.1 }
274 pam-fi 1.13 /**
275     * Set the TrkTrack position measurements
276     */
277     void TrkTrack::SetMeasure(double *xmeas, double *ymeas, double *zmeas){
278     for(int i=0; i<6; i++) xm[i]=*xmeas++;
279     for(int i=0; i<6; i++) ym[i]=*ymeas++;
280     for(int i=0; i<6; i++) zm[i]=*zmeas++;
281     }
282     /**
283     * Set the TrkTrack position resolution
284     */
285     void TrkTrack::SetResolution(double *rx, double *ry){
286     for(int i=0; i<6; i++) resx[i]=*rx++;
287     for(int i=0; i<6; i++) resy[i]=*ry++;
288     }
289     /**
290     * Set the TrkTrack good measurement
291     */
292     void TrkTrack::SetGood(int *xg, int *yg){
293 pam-fi 1.30 // NB! si perdera` l'informazione sul numero del cluster
294 pam-fi 1.13 for(int i=0; i<6; i++) xgood[i]=*xg++;
295     for(int i=0; i<6; i++) ygood[i]=*yg++;
296     }
297    
298     /**
299     * Load the magnetic field
300     */
301 pam-fi 1.16 void TrkTrack::LoadField(TString path){
302    
303 pam-fi 1.26 // strcpy(path_.path,path.Data());
304     // path_.pathlen = path.Length();
305     // path_.error = 0;
306     // readb_();
307    
308     TrkParams::Set(path,1);
309 pam-fi 1.28 TrkParams::Load(1);
310 pam-fi 1.16
311 pam-fi 1.13 };
312 pam-fi 1.19
313 pam-fi 1.25
314 pam-fi 1.19 /**
315     * Method to fill minimization-routine common
316     */
317     void TrkTrack::FillMiniStruct(cMini2track& track){
318    
319     for(int i=0; i<6; i++){
320    
321 pam-fi 1.30 // track.xgood[i]=xgood[i];
322     // track.ygood[i]=ygood[i];
323     track.xgood[i]=XGood(i);
324     track.ygood[i]=YGood(i);
325 pam-fi 1.19
326     track.xm[i]=xm[i];
327     track.ym[i]=ym[i];
328     track.zm[i]=zm[i];
329    
330     // --- temporaneo ----------------------------
331     // andrebbe inserita la dimensione del sensore
332     float segment = 100.;
333     track.xm_a[i]=xm[i];
334     track.xm_b[i]=xm[i];
335     track.ym_a[i]=ym[i];
336     track.ym_b[i]=ym[i];
337 pam-fi 1.30 if( XGood(i) && !YGood(i) ){
338 pam-fi 1.19 track.ym_a[i] = track.ym_a[i]+segment;
339     track.ym_b[i] = track.ym_b[i]-segment;
340 pam-fi 1.30 }else if( !XGood(i) && YGood(i)){
341 pam-fi 1.19 track.xm_a[i] = track.xm_a[i]+segment;
342     track.xm_b[i] = track.xm_b[i]-segment;
343     }
344     // --- temporaneo ----------------------------
345    
346     track.resx[i]=resx[i];
347     track.resy[i]=resy[i];
348     }
349    
350     for(int i=0; i<5; i++) track.al[i]=al[i];
351     track.zini = 23.5;
352     // ZINI = 23.5 !!! it should be the same parameter in all codes
353    
354     }
355     /**
356     * Method to set values from minimization-routine common
357     */
358     void TrkTrack::SetFromMiniStruct(cMini2track *track){
359    
360     for(int i=0; i<5; i++) {
361     al[i]=track->al[i];
362     for(int j=0; j<5; j++) coval[i][j]=track->cov[i][j];
363     }
364     chi2 = track->chi2;
365     nstep = track->nstep;
366     for(int i=0; i<6; i++){
367     xv[i] = track->xv[i];
368     yv[i] = track->yv[i];
369     zv[i] = track->zv[i];
370     xm[i] = track->xm[i];
371     ym[i] = track->ym[i];
372     zm[i] = track->zm[i];
373     axv[i] = track->axv[i];
374     ayv[i] = track->ayv[i];
375     }
376    
377     }
378 pam-fi 1.13 /**
379     * Tracking method. It calls F77 mini routine.
380     */
381     void TrkTrack::Fit(double pfixed, int& fail, int iprint){
382 pam-fi 1.15
383     float al_ini[] = {0.,0.,0.,0.,0.};
384    
385 pam-fi 1.13 extern cMini2track track_;
386     fail = 0;
387 pam-fi 1.19 FillMiniStruct(track_);
388 pam-fi 1.15
389 pam-fi 1.19 // if fit variables have been reset, evaluate the initial guess
390 pam-fi 1.15 if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.)guess_();
391    
392 pam-fi 1.16 // --------------------- free momentum
393 pam-fi 1.14 if(pfixed==0.) {
394 pam-fi 1.19 track_.pfixed=0.;
395 pam-fi 1.14 }
396 pam-fi 1.16 // --------------------- fixed momentum
397 pam-fi 1.13 if(pfixed!=0.) {
398 pam-fi 1.19 al[4]=1./pfixed;
399     track_.pfixed=pfixed;
400 pam-fi 1.13 }
401 pam-fi 1.15
402 pam-fi 1.19 // store temporarily the initial guess
403 pam-fi 1.15 for(int i=0; i<5; i++) al_ini[i]=track_.al[i];
404    
405 pam-fi 1.19 // ------------------------------------------
406     // call mini routine
407 pam-fi 1.26 TrkParams::Load(1);
408     if( !TrkParams::IsLoaded(1) ){
409     cout << "void TrkTrack::Fit(double pfixed, int& fail, int iprint) --- ERROR --- m.field not loaded"<<endl;
410     return;
411     }
412 pam-fi 1.13 int istep=0;
413     int ifail=0;
414     mini2_(&istep,&ifail, &iprint);
415     if(ifail!=0) {
416 pam-fi 1.19 if(iprint)cout << "ERROR: ifail= " << ifail << endl;
417 pam-fi 1.13 fail = 1;
418     }
419 pam-fi 1.19 // ------------------------------------------
420 pam-fi 1.15
421 pam-fi 1.19 SetFromMiniStruct(&track_);
422 pam-fi 1.13 // cout << endl << "eta ===> " << track_.al[4] << endl;
423    
424 pam-fi 1.19 // for(int i=0; i<5; i++) al[i]=track_.al[i];
425     // chi2=track_.chi2;
426     // nstep=track_.nstep;
427     // for(int i=0; i<6; i++) xv[i]=track_.xv[i];
428     // for(int i=0; i<6; i++) yv[i]=track_.yv[i];
429     // for(int i=0; i<6; i++) zv[i]=track_.zv[i];
430     // for(int i=0; i<6; i++) axv[i]=track_.axv[i];
431     // for(int i=0; i<6; i++) ayv[i]=track_.ayv[i];
432     // for(int i=0; i<5; i++) {
433     // for(int j=0; j<5; j++) coval[i][j]=track_.cov[i][j];
434     // }
435 pam-fi 1.15
436 pam-fi 1.20 if(fail){
437     if(iprint)cout << " >>>> fit failed >>>> drawing initial par"<<endl;
438     for(int i=0; i<5; i++) al[i]=al_ini[i];
439     }
440 pam-fi 1.15
441 pam-fi 1.13 };
442     /*
443 pam-fi 1.15 * Reset the fit parameters
444 pam-fi 1.13 */
445     void TrkTrack::FitReset(){
446     for(int i=0; i<5; i++) al[i]=-9999.;
447     chi2=0.;
448     nstep=0;
449     for(int i=0; i<6; i++) xv[i]=0.;
450     for(int i=0; i<6; i++) yv[i]=0.;
451     for(int i=0; i<6; i++) zv[i]=0.;
452     for(int i=0; i<6; i++) axv[i]=0.;
453     for(int i=0; i<6; i++) ayv[i]=0.;
454     for(int i=0; i<5; i++) {
455     for(int j=0; j<5; j++) coval[i][j]=0.;
456     }
457     }
458 pam-fi 1.31 /*
459     * Set the tracking mode
460     */
461 pam-fi 1.27 void TrkTrack::SetTrackingMode(int trackmode){
462     extern cMini2track track_;
463     track_.trackmode = trackmode;
464     }
465 pam-fi 1.31 /*
466     * Set the factor scale for tracking precision
467     */
468     void TrkTrack::SetPrecisionFactor(double fact){
469     extern cMini2track track_;
470     track_.fact = fact;
471     }
472     /*
473     * Set the factor scale for tracking precision
474     */
475     void TrkTrack::SetStepMin(int istepmin){
476     extern cMini2track track_;
477     track_.istepmin = istepmin;
478     }
479 pam-fi 1.13
480 pam-fi 1.29
481     /*
482     * Method to retrieve the X-view clusters associated to the track.
483     * @param ip Tracker plane (0-5)
484     */
485     TrkCluster *TrkTrack::GetClusterX(int ip){
486     cout << " TrkCluster *TrkTrack::GetClusterX(int ip) -- momentaneamente fuori servizio --"<< endl;
487     if(!clx)return NULL;
488     TrkCluster *pt = (TrkCluster*)(clx->At(ip));
489     return pt;
490     };
491     /*
492     * Method to retrieve the Y-view clusters associated to the track.
493     * @param ip Tracker plane (0-5)
494     */
495     TrkCluster *TrkTrack::GetClusterY(int ip){
496     cout << " TrkCluster *TrkTrack::GetClusterY(int ip) -- momentaneamente fuori servizio --"<< endl;
497     if(!cly)return NULL;
498     TrkCluster *pt = (TrkCluster*)(cly->At(ip));
499     return pt;
500     };
501    
502    
503 mocchiut 1.1 //--------------------------------------
504     //
505     //
506     //--------------------------------------
507 pam-fi 1.10 void TrkTrack::Clear(){
508 pam-fi 1.21 // cout << "TrkTrack::Clear()"<<endl;
509     seqno = -1;
510     image = -1;
511     chi2 = 0;
512     nstep = 0;
513     for(int it1=0;it1<5;it1++){
514     al[it1] = 0;
515     for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;
516     };
517     for(int ip=0;ip<6;ip++){
518     xgood[ip] = 0;
519     ygood[ip] = 0;
520     xm[ip] = 0;
521     ym[ip] = 0;
522     zm[ip] = 0;
523     resx[ip] = 0;
524     resy[ip] = 0;
525     xv[ip] = 0;
526     yv[ip] = 0;
527     zv[ip] = 0;
528     axv[ip] = 0;
529     ayv[ip] = 0;
530     dedx_x[ip] = 0;
531     dedx_y[ip] = 0;
532    
533     };
534     if(clx)clx->Clear();
535     if(cly)cly->Clear();
536 pam-fi 1.10 };
537     //--------------------------------------
538     //
539     //
540     //--------------------------------------
541 pam-fi 1.11 void TrkTrack::Delete(){
542 pam-fi 1.21 // cout << "TrkTrack::Delete()"<<endl;
543     // Clear();
544     if(clx)delete clx;
545     if(cly)delete cly;
546 pam-fi 1.11 };
547 pam-fi 1.21 //--------------------------------------
548 pam-fi 1.11 //
549     //
550     //--------------------------------------
551 pam-fi 1.10
552     //--------------------------------------
553     //
554     //
555     //--------------------------------------
556 mocchiut 1.1 TrkSinglet::TrkSinglet(){
557 pam-fi 1.21 // cout << "TrkSinglet::TrkSinglet() " << GetUniqueID()<<endl;
558 mocchiut 1.1 plane = 0;
559     coord[0] = 0;
560     coord[1] = 0;
561     sgnl = 0;
562 pam-fi 1.21 // cls = 0;
563 mocchiut 1.1 };
564     //--------------------------------------
565     //
566     //
567     //--------------------------------------
568     TrkSinglet::TrkSinglet(const TrkSinglet& s){
569 pam-fi 1.21 // cout << "TrkSinglet::TrkSinglet(const TrkSinglet& s) " << GetUniqueID()<<endl;
570 mocchiut 1.1 plane = s.plane;
571     coord[0] = s.coord[0];
572     coord[1] = s.coord[1];
573     sgnl = s.sgnl;
574 pam-fi 1.9 // cls = 0;//<<<<pointer
575 pam-fi 1.21 cls = TRef(s.cls);
576 mocchiut 1.1 };
577     //--------------------------------------
578     //
579     //
580     //--------------------------------------
581     void TrkSinglet::Dump(){
582     int i=0;
583     cout << endl << "========== Singlet " ;
584     cout << endl << "plane : " << plane;
585     cout << endl << "coord[2] : "; while( i<2 && cout << coord[i] << " ") i++;
586     cout << endl << "sgnl : " << sgnl;
587     }
588     //--------------------------------------
589     //
590     //
591     //--------------------------------------
592 pam-fi 1.21 void TrkSinglet::Clear(){
593     // cout << "TrkSinglet::Clear() " << GetUniqueID()<<endl;
594     // cls=0;
595     plane=-1;
596     coord[0]=-999;
597     coord[1]=-999;
598     sgnl=0;
599    
600     }
601     //--------------------------------------
602     //
603     //
604     //--------------------------------------
605 mocchiut 1.1 TrkLevel2::TrkLevel2(){
606 mocchiut 1.24 // cout <<"TrkLevel2::TrkLevel2()"<<endl;
607 mocchiut 1.1 for(Int_t i=0; i<12 ; i++){
608 pam-fi 1.10 good[i] = -1;
609     };
610 pam-fi 1.23 // okkio!! memory-leak
611     // Track = new TClonesArray("TrkTrack");
612     // SingletX = new TClonesArray("TrkSinglet");
613     // SingletY = new TClonesArray("TrkSinglet");
614 pam-fi 1.21 Track = 0;
615     SingletX = 0;
616     SingletY = 0;
617 pam-fi 1.6
618 mocchiut 1.1 }
619     //--------------------------------------
620     //
621     //
622     //--------------------------------------
623 pam-fi 1.23 void TrkLevel2::Set(){
624     if(!Track)Track = new TClonesArray("TrkTrack");
625     if(!SingletX)SingletX = new TClonesArray("TrkSinglet");
626     if(!SingletY)SingletY = new TClonesArray("TrkSinglet");
627     }
628     //--------------------------------------
629     //
630     //
631     //--------------------------------------
632 mocchiut 1.1 void TrkLevel2::Dump(){
633 pam-fi 1.10
634     //
635 mocchiut 1.1 cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-";
636 pam-fi 1.10 cout << endl << "good : "; for(int i=0; i<12; i++) cout << good[i]<<" ";
637 mocchiut 1.1 cout << endl << "ntrk() : " << this->ntrk() ;
638     cout << endl << "nclsx() : " << this->nclsx();
639     cout << endl << "nclsy() : " << this->nclsy();
640 pam-fi 1.21 // TClonesArray &t = *Track;
641     // TClonesArray &sx = *SingletX;
642     // TClonesArray &sy = *SingletY;
643     // for(int i=0; i<ntrk(); i++) ((TrkTrack *)t[i])->Dump();
644     // for(int i=0; i<nclsx(); i++) ((TrkSinglet *)sx[i])->Dump();
645     // for(int i=0; i<nclsy(); i++) ((TrkSinglet *)sy[i])->Dump();
646     if(Track){
647     TClonesArray &t = *Track;
648     for(int i=0; i<ntrk(); i++) ((TrkTrack *)t[i])->Dump();
649     }
650     if(SingletX){
651     TClonesArray &sx = *SingletX;
652     for(int i=0; i<nclsx(); i++) ((TrkSinglet *)sx[i])->Dump();
653     }
654     if(SingletY){
655     TClonesArray &sy = *SingletY;
656     for(int i=0; i<nclsy(); i++) ((TrkSinglet *)sy[i])->Dump();
657     }
658 mocchiut 1.1 }
659     //--------------------------------------
660     //
661     //
662     //--------------------------------------
663     /**
664     * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common).
665     */
666 pam-fi 1.26 // void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2){
667 pam-fi 1.10
668 pam-fi 1.26 // // temporary objects:
669     // TrkSinglet* t_singlet = new TrkSinglet();
670     // TrkTrack* t_track = new TrkTrack();
671    
672     // // **** general variables ****
673     // // good2 = l2->good2;
674     // for(Int_t i=0; i<12 ; i++){
675     // // crc[i] = l2->crc[i];
676     // good[i] = l2->good[i];
677     // };
678     // // *** TRACKS ***
679     // if(!Track) Track = new TClonesArray("TrkTrack");
680     // TClonesArray &t = *Track;
681     // for(int i=0; i<l2->ntrk; i++){
682     // t_track->seqno = i;// NBNBNBNB deve sempre essere = i
683     // t_track->image = l2->image[i]-1;
684     // // cout << "track "<<i<<t_track->seqno << t_track->image<<endl;
685     // t_track->chi2 = l2->chi2_nt[i];
686     // t_track->nstep = l2->nstep_nt[i];
687     // for(int it1=0;it1<5;it1++){
688     // t_track->al[it1] = l2->al_nt[i][it1];
689     // for(int it2=0;it2<5;it2++)
690     // t_track->coval[it1][it2] = l2->coval[i][it2][it1];
691     // };
692     // for(int ip=0;ip<6;ip++){
693     // t_track->xgood[ip] = l2->xgood_nt[i][ip];
694     // t_track->ygood[ip] = l2->ygood_nt[i][ip];
695     // t_track->xm[ip] = l2->xm_nt[i][ip];
696     // t_track->ym[ip] = l2->ym_nt[i][ip];
697     // t_track->zm[ip] = l2->zm_nt[i][ip];
698     // t_track->resx[ip] = l2->resx_nt[i][ip];
699     // t_track->resy[ip] = l2->resy_nt[i][ip];
700     // t_track->xv[ip] = l2->xv_nt[i][ip];
701     // t_track->yv[ip] = l2->yv_nt[i][ip];
702     // t_track->zv[ip] = l2->zv_nt[i][ip];
703     // t_track->axv[ip] = l2->axv_nt[i][ip];
704     // t_track->ayv[ip] = l2->ayv_nt[i][ip];
705     // t_track->dedx_x[ip] = l2->dedx_x[i][ip];
706     // t_track->dedx_y[ip] = l2->dedx_y[i][ip];
707     // // t_track->clx[ip] = 0;
708     // // t_track->cly[ip] = 0;
709     // };
710     // new(t[i]) TrkTrack(*t_track);
711     // t_track->Clear();
712     // };
713     // // *** SINGLETS ***
714     // if(!SingletX)SingletX = new TClonesArray("TrkSinglet");
715     // TClonesArray &sx = *SingletX;
716     // for(int i=0; i<l2->nclsx; i++){
717     // t_singlet->plane = l2->planex[i];
718     // t_singlet->coord[0] = l2->xs[i][0];
719     // t_singlet->coord[1] = l2->xs[i][1];
720     // t_singlet->sgnl = l2->signlxs[i];
721     // // t_singlet->cls = 0;
722     // new(sx[i]) TrkSinglet(*t_singlet);
723     // t_singlet->Clear();
724     // }
725     // if(!SingletY)SingletY = new TClonesArray("TrkSinglet");
726     // TClonesArray &sy = *SingletY;
727     // for(int i=0; i<l2->nclsy; i++){
728     // t_singlet->plane = l2->planey[i];
729     // t_singlet->coord[0] = l2->ys[i][0];
730     // t_singlet->coord[1] = l2->ys[i][1];
731     // t_singlet->sgnl = l2->signlys[i];
732     // // t_singlet->cls = 0;
733     // new(sy[i]) TrkSinglet(*t_singlet);
734     // t_singlet->Clear();
735     // };
736 pam-fi 1.9
737 pam-fi 1.26 // delete t_track;
738     // delete t_singlet;
739     // }
740 pam-fi 1.9 //--------------------------------------
741     //
742     //
743     //--------------------------------------
744     /**
745     * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common).
746 pam-fi 1.29 * Ref to Level1 data (clusters) is also set. If l1==NULL no references are set.
747     * (NB It make sense to set references only if events are stored in a tree that contains also the Level1 branch)
748 pam-fi 1.9 */
749     void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1){
750    
751 pam-fi 1.29 // cout << "void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1)"<<endl;
752     Clear();
753 pam-fi 1.9 // temporary objects:
754 pam-fi 1.29 TrkSinglet* t_singlet = new TrkSinglet();
755     TrkTrack* t_track = new TrkTrack();
756    
757     // -----------------
758     // general variables
759     // -----------------
760     for(Int_t i=0; i<12 ; i++){
761     good[i] = l2->good[i];
762     };
763     // --------------
764     // *** TRACKS ***
765     // --------------
766     if(!Track) Track = new TClonesArray("TrkTrack");
767     TClonesArray &t = *Track;
768     //-----------------------------------------------------
769     if( l1 && !t_track->clx )t_track->clx = new TRefArray(6,0);
770     if( l1 && !t_track->cly )t_track->cly = new TRefArray(6,0);
771     //-----------------------------------------------------
772     for(int i=0; i<l2->ntrk; i++){
773     // cout <<" TRACK "<<i<<" ------------------ "<<endl;
774     t_track->seqno = i;// NBNBNBNB deve sempre essere = i
775     t_track->image = l2->image[i]-1;
776     t_track->chi2 = l2->chi2_nt[i];
777     t_track->nstep = l2->nstep_nt[i];
778     for(int it1=0;it1<5;it1++){
779     t_track->al[it1] = l2->al_nt[i][it1];
780     for(int it2=0;it2<5;it2++)
781     t_track->coval[it1][it2] = l2->coval[i][it2][it1];
782 pam-fi 1.9 };
783 pam-fi 1.29 for(int ip=0;ip<6;ip++){
784     t_track->xgood[ip] = l2->cltrx[i][ip];//l2->xgood_nt[i][ip];
785     t_track->ygood[ip] = l2->cltry[i][ip];//l2->ygood_nt[i][ip];
786     t_track->xm[ip] = l2->xm_nt[i][ip];
787     t_track->ym[ip] = l2->ym_nt[i][ip];
788     t_track->zm[ip] = l2->zm_nt[i][ip];
789     t_track->resx[ip] = l2->resx_nt[i][ip];
790     t_track->resy[ip] = l2->resy_nt[i][ip];
791     t_track->xv[ip] = l2->xv_nt[i][ip];
792     t_track->yv[ip] = l2->yv_nt[i][ip];
793     t_track->zv[ip] = l2->zv_nt[i][ip];
794     t_track->axv[ip] = l2->axv_nt[i][ip];
795     t_track->ayv[ip] = l2->ayv_nt[i][ip];
796     t_track->dedx_x[ip] = l2->dedx_x[i][ip];
797     t_track->dedx_y[ip] = l2->dedx_y[i][ip];
798     // cout << " ip "<<ip<<" l2->cltrx[i][ip] "<< l2->cltrx[i][ip]<<endl;
799     // cout << " ip "<<ip<<" l2->cltry[i][ip] "<< l2->cltry[i][ip]<<endl;
800     //-----------------------------------------------------
801     //-----------------------------------------------------
802     // if(l1 && t_track->xgood[ip])t_track->clx->AddAt(l1->GetCluster(l2->cltrx[i][ip]-1),ip);
803     // if(l1 && t_track->ygood[ip])t_track->cly->AddAt(l1->GetCluster(l2->cltry[i][ip]-1),ip);
804     if(l2->xgood_nt[i][ip]){
805     // cout << " ip "<<ip<<" l2->cltrx[i][ip] "<< l2->cltrx[i][ip]<<endl;
806     // cout << " ip "<<ip<<" l2->cltrx[i][ip] "<< l2->cltrx[i][ip]<<" ";
807     // if( l1->GetCluster(l2->cltrx[i][ip]-1)->TestBit(l1->GetCluster(l2->cltrx[i][ip]-1)->kIsReferenced) )cout << ">> is referenced ";
808     if(l1)t_track->clx->AddAt(l1->GetCluster(l2->cltrx[i][ip]-1),ip);
809     // cout << " --- "<<l1->GetCluster(l2->cltrx[i][ip]-1)->GetUniqueID()<<endl;
810     // t_track->xgood[ip] = l2->cltrx[i][ip]; // WORK-AROUND *****
811     }else{
812     if(l1)t_track->clx->RemoveAt(ip);
813     }
814     if(l2->ygood_nt[i][ip]){
815     // cout << " ip "<<ip<<" l2->cltry[i][ip] "<< l2->cltry[i][ip]<<endl;
816     // cout << " ip "<<ip<<" l2->cltry[i][ip] "<< l2->cltry[i][ip]<<" ";
817     // if( l1->GetCluster(l2->cltry[i][ip]-1)->TestBit(l1->GetCluster(l2->cltry[i][ip]-1)->kIsReferenced) )cout << ">> is referenced ";
818     if(l1)t_track->cly->AddAt(l1->GetCluster(l2->cltry[i][ip]-1),ip);
819     // cout << " --- "<<l1->GetCluster(l2->cltry[i][ip]-1)->GetUniqueID()<<endl;
820     // t_track->ygood[ip] = l2->cltry[i][ip]; // WORK-AROUND *****
821     }else{
822     if(l1)t_track->cly->RemoveAt(ip);
823     }
824     //-----------------------------------------------------
825     //-----------------------------------------------------
826     };
827     new(t[i]) TrkTrack(*t_track);
828     t_track->Clear();
829     };
830     // ----------------
831     // *** SINGLETS ***
832     // ----------------
833     if(!SingletX)SingletX = new TClonesArray("TrkSinglet");
834     TClonesArray &sx = *SingletX;
835     for(int i=0; i<l2->nclsx; i++){
836     t_singlet->plane = l2->planex[i];
837     t_singlet->coord[0] = l2->xs[i][0];
838     t_singlet->coord[1] = l2->xs[i][1];
839     t_singlet->sgnl = l2->signlxs[i];
840     //-----------------------------------------------------
841     if(l1) t_singlet->cls = l1->GetCluster(l2->clsx[i]-1);
842     //-----------------------------------------------------
843     new(sx[i]) TrkSinglet(*t_singlet);
844     t_singlet->Clear();
845     }
846     if(!SingletY)SingletY = new TClonesArray("TrkSinglet");
847     TClonesArray &sy = *SingletY;
848     for(int i=0; i<l2->nclsy; i++){
849     t_singlet->plane = l2->planey[i];
850     t_singlet->coord[0] = l2->ys[i][0];
851     t_singlet->coord[1] = l2->ys[i][1];
852     t_singlet->sgnl = l2->signlys[i];
853 pam-fi 1.26 //-----------------------------------------------------
854 pam-fi 1.29 if(l1) t_singlet->cls = l1->GetCluster(l2->clsy[i]-1);
855 pam-fi 1.26 //-----------------------------------------------------
856 pam-fi 1.29 new(sy[i]) TrkSinglet(*t_singlet);
857     t_singlet->Clear();
858     };
859 pam-fi 1.5
860 pam-fi 1.29 delete t_track;
861     delete t_singlet;
862 mocchiut 1.1 }
863 pam-fi 1.7 /**
864     * Fills a struct cTrkLevel2 with values from a TrkLevel2 object (to put data into a F77 common).
865     */
866    
867     void TrkLevel2::GetLevel2Struct(cTrkLevel2 *l2) const {
868    
869     // general variables
870 pam-fi 1.10 // l2->good2 = good2 ;
871 pam-fi 1.7 for(Int_t i=0; i<12 ; i++){
872 pam-fi 1.10 // l2->crc[i] = crc[i];
873     l2->good[i] = good[i];
874 pam-fi 1.7 };
875     // *** TRACKS ***
876    
877 pam-fi 1.21 if(Track){
878     l2->ntrk = Track->GetEntries();
879     for(Int_t i=0;i<l2->ntrk;i++){
880     l2->image[i] = 1 + ((TrkTrack *)Track->At(i))->image;
881     l2->chi2_nt[i] = ((TrkTrack *)Track->At(i))->chi2;
882     l2->nstep_nt[i] = ((TrkTrack *)Track->At(i))->nstep;
883     for(int it1=0;it1<5;it1++){
884     l2->al_nt[i][it1] = ((TrkTrack *)Track->At(i))->al[it1];
885     for(int it2=0;it2<5;it2++)
886     l2->coval[i][it2][it1] = ((TrkTrack *)Track->At(i))->coval[it1][it2];
887     };
888     for(int ip=0;ip<6;ip++){
889 pam-fi 1.30 l2->xgood_nt[i][ip] = ((TrkTrack *)Track->At(i))->XGood(ip);
890     l2->ygood_nt[i][ip] = ((TrkTrack *)Track->At(i))->YGood(ip);
891 pam-fi 1.21 l2->xm_nt[i][ip] = ((TrkTrack *)Track->At(i))->xm[ip];
892     l2->ym_nt[i][ip] = ((TrkTrack *)Track->At(i))->ym[ip];
893     l2->zm_nt[i][ip] = ((TrkTrack *)Track->At(i))->zm[ip];
894     l2->resx_nt[i][ip] = ((TrkTrack *)Track->At(i))->resx[ip];
895     l2->resy_nt[i][ip] = ((TrkTrack *)Track->At(i))->resy[ip];
896     l2->xv_nt[i][ip] = ((TrkTrack *)Track->At(i))->xv[ip];
897     l2->yv_nt[i][ip] = ((TrkTrack *)Track->At(i))->yv[ip];
898     l2->zv_nt[i][ip] = ((TrkTrack *)Track->At(i))->zv[ip];
899     l2->axv_nt[i][ip] = ((TrkTrack *)Track->At(i))->axv[ip];
900     l2->ayv_nt[i][ip] = ((TrkTrack *)Track->At(i))->ayv[ip];
901     l2->dedx_x[i][ip] = ((TrkTrack *)Track->At(i))->dedx_x[ip];
902     l2->dedx_y[i][ip] = ((TrkTrack *)Track->At(i))->dedx_y[ip];
903     };
904     }
905     }
906     // *** SINGLETS ***
907     if(SingletX){
908     l2->nclsx = SingletX->GetEntries();
909     for(Int_t i=0;i<l2->nclsx;i++){
910     l2->planex[i] = ((TrkSinglet *)SingletX->At(i))->plane;
911     l2->xs[i][0] = ((TrkSinglet *)SingletX->At(i))->coord[0];
912     l2->xs[i][1] = ((TrkSinglet *)SingletX->At(i))->coord[1];
913     l2->signlxs[i] = ((TrkSinglet *)SingletX->At(i))->sgnl;
914     }
915 pam-fi 1.7 }
916    
917 pam-fi 1.21 if(SingletY){
918     l2->nclsy = SingletY->GetEntries();
919     for(Int_t i=0;i<l2->nclsy;i++){
920     l2->planey[i] = ((TrkSinglet *)SingletY->At(i))->plane;
921     l2->ys[i][0] = ((TrkSinglet *)SingletY->At(i))->coord[0];
922     l2->ys[i][1] = ((TrkSinglet *)SingletY->At(i))->coord[1];
923     l2->signlys[i] = ((TrkSinglet *)SingletY->At(i))->sgnl;
924     }
925 pam-fi 1.7 }
926     }
927 mocchiut 1.1 //--------------------------------------
928     //
929     //
930     //--------------------------------------
931     void TrkLevel2::Clear(){
932     for(Int_t i=0; i<12 ; i++){
933 pam-fi 1.21 good[i] = -1;
934     };
935     // if(Track)Track->Clear("C");
936     // if(SingletX)SingletX->Clear("C");
937     // if(SingletY)SingletY->Clear("C");
938     if(Track)Track->Delete();
939     if(SingletX)SingletX->Delete();
940     if(SingletY)SingletY->Delete();
941     }
942     // //--------------------------------------
943     // //
944     // //
945     // //--------------------------------------
946 pam-fi 1.11 void TrkLevel2::Delete(){
947    
948 pam-fi 1.21 // cout << "void TrkLevel2::Delete()"<<endl;
949     Clear();
950     if(Track)delete Track;
951     if(SingletX)delete SingletX;
952     if(SingletY)delete SingletY;
953    
954 pam-fi 1.11 }
955     //--------------------------------------
956     //
957     //
958     //--------------------------------------
959 mocchiut 1.1 /**
960     * 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).
961     * This method is overridden by PamLevel2::GetTracks(), where calorimeter and TOF information is used.
962     */
963 pam-fi 1.8 TRefArray *TrkLevel2::GetTracks_NFitSorted(){
964 pam-fi 1.3
965 pam-fi 1.21 if(!Track)return 0;
966    
967     TRefArray *sorted = new TRefArray();
968 pam-fi 1.8
969 pam-fi 1.21 TClonesArray &t = *Track;
970 pam-fi 1.8 // TClonesArray &ts = *PhysicalTrack;
971 pam-fi 1.21 int N = ntrk();
972     vector<int> m(N); for(int i=0; i<N; i++)m[i]=1;
973 pam-fi 1.8 // int m[50]; for(int i=0; i<N; i++)m[i]=1;
974    
975 pam-fi 1.21 int indo=0;
976     int indi=0;
977 pam-fi 1.26 while(N > 0){
978     // while(N != 0){
979 pam-fi 1.21 int nfit =0;
980     float chi2ref = numeric_limits<float>::max();
981 pam-fi 1.9
982 pam-fi 1.21 // first loop to search maximum num. of fit points
983     for(int i=0; i < ntrk(); i++){
984     if( ((TrkTrack *)t[i])->GetNtot() >= nfit && m[i]==1){
985     nfit = ((TrkTrack *)t[i])->GetNtot();
986     }
987     }
988     //second loop to search minimum chi2 among selected
989 pam-fi 1.26 for(int i=0; i<ntrk(); i++){
990 pam-fi 1.21 Float_t chi2 = ((TrkTrack *)t[i])->chi2;
991 pam-fi 1.26 if(chi2 < 0) chi2 = -chi2*1000;
992 pam-fi 1.21 if( chi2 < chi2ref
993     && ((TrkTrack *)t[i])->GetNtot() == nfit
994     && m[i]==1){
995     chi2ref = ((TrkTrack *)t[i])->chi2;
996     indi = i;
997     };
998     };
999     if( ((TrkTrack *)t[indi])->HasImage() ){
1000     m[((TrkTrack *)t[indi])->image] = 0;
1001     N--;
1002 pam-fi 1.8
1003 pam-fi 1.26 // cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl;
1004 pam-fi 1.21 };
1005     sorted->Add( (TrkTrack*)t[indi] );
1006 pam-fi 1.3
1007 pam-fi 1.21 m[indi] = 0;
1008 pam-fi 1.26 // cout << "SORTED "<< indo << " "<< indi << " "<< N << " "<<((TrkTrack *)t[indi])->image<<" "<<chi2ref<<endl;
1009 pam-fi 1.21 N--;
1010     indo++;
1011     }
1012     m.clear();
1013 pam-fi 1.26 // cout << "GetTracks_NFitSorted(it): Done"<< endl;
1014 pam-fi 1.8
1015 pam-fi 1.21 return sorted;
1016 pam-fi 1.6 // return PhysicalTrack;
1017 pam-fi 1.3 }
1018 mocchiut 1.1 //--------------------------------------
1019     //
1020     //
1021     //--------------------------------------
1022     /**
1023     * Retrieves the is-th stored track.
1024     * @param it Track number, ranging from 0 to ntrk().
1025     * Fitted tracks ( images included ) are stored in a TObjectArray ( TrkLevel2::Track ) in the same order they are returned by the F77 fitting routine.
1026     */
1027     TrkTrack *TrkLevel2::GetStoredTrack(int is){
1028    
1029     if(is >= this->ntrk()){
1030     cout << "** TrkLevel2 ** Track "<< is << "doen not exits! " << endl;
1031     cout << " Stored tracks ntrk() = "<< this->ntrk() << endl;
1032     return 0;
1033     }
1034 pam-fi 1.21 if(!Track){
1035     cout << "TrkTrack *TrkLevel2::GetStoredTrack(int is) >> (TClonesArray*) Track ==0 "<<endl;
1036     };
1037 mocchiut 1.1 TClonesArray &t = *(Track);
1038     TrkTrack *track = (TrkTrack*)t[is];
1039     return track;
1040     }
1041     //--------------------------------------
1042     //
1043     //
1044     //--------------------------------------
1045     /**
1046 pam-fi 1.6 * Retrieves the is-th stored X singlet.
1047     * @param it Singlet number, ranging from 0 to nclsx().
1048     */
1049     TrkSinglet *TrkLevel2::GetSingletX(int is){
1050    
1051     if(is >= this->nclsx()){
1052     cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl;
1053     cout << " Stored x-singlets nclsx() = "<< this->nclsx() << endl;
1054     return 0;
1055     }
1056 pam-fi 1.21 if(!SingletX)return 0;
1057 pam-fi 1.6 TClonesArray &t = *(SingletX);
1058     TrkSinglet *singlet = (TrkSinglet*)t[is];
1059     return singlet;
1060     }
1061     //--------------------------------------
1062     //
1063     //
1064     //--------------------------------------
1065     /**
1066     * Retrieves the is-th stored Y singlet.
1067     * @param it Singlet number, ranging from 0 to nclsx().
1068     */
1069     TrkSinglet *TrkLevel2::GetSingletY(int is){
1070    
1071     if(is >= this->nclsy()){
1072     cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl;
1073     cout << " Stored y-singlets nclsy() = "<< this->nclsx() << endl;
1074     return 0;
1075     }
1076 pam-fi 1.21 if(!SingletY)return 0;
1077 pam-fi 1.6 TClonesArray &t = *(SingletY);
1078     TrkSinglet *singlet = (TrkSinglet*)t[is];
1079     return singlet;
1080     }
1081     //--------------------------------------
1082     //
1083     //
1084     //--------------------------------------
1085     /**
1086 mocchiut 1.1 * Retrieves the it-th "physical" track, sorted by the method GetNTracks().
1087     * @param it Track number, ranging from 0 to GetNTracks().
1088     */
1089 pam-fi 1.10
1090 pam-fi 1.8 TrkTrack *TrkLevel2::GetTrack(int it){
1091    
1092     if(it >= this->GetNTracks()){
1093     cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl;
1094     cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl;
1095     return 0;
1096     }
1097    
1098     TRefArray *sorted = GetTracks(); //TEMPORANEO
1099 pam-fi 1.21 if(!sorted)return 0;
1100 pam-fi 1.8 TrkTrack *track = (TrkTrack*)sorted->At(it);
1101 pam-fi 1.21 sorted->Clear();
1102     delete sorted;
1103 pam-fi 1.8 return track;
1104 mocchiut 1.1 }
1105 pam-fi 1.6 /**
1106     * Give the number of "physical" tracks, sorted by the method GetTracks().
1107     */
1108 pam-fi 1.5 Int_t TrkLevel2::GetNTracks(){
1109 pam-fi 1.8
1110     Float_t ntot=0;
1111 pam-fi 1.21 if(!Track)return 0;
1112 pam-fi 1.8 TClonesArray &t = *Track;
1113 mocchiut 1.12 for(int i=0; i<ntrk(); i++) {
1114 pam-fi 1.8 if( ((TrkTrack *)t[i])->GetImageSeqNo() == -1 ) ntot+=1.;
1115     else ntot+=0.5;
1116     }
1117     return (Int_t)ntot;
1118    
1119 pam-fi 1.5 };
1120 mocchiut 1.1 //--------------------------------------
1121     //
1122     //
1123     //--------------------------------------
1124     /**
1125     * Retrieves (if present) the image of the it-th "physical" track, sorted by the method GetNTracks().
1126     * @param it Track number, ranging from 0 to GetNTracks().
1127     */
1128 pam-fi 1.8 TrkTrack *TrkLevel2::GetTrackImage(int it){
1129    
1130 pam-fi 1.21 if(it >= this->GetNTracks()){
1131     cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl;
1132     cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl;
1133     return 0;
1134     }
1135 pam-fi 1.8
1136 pam-fi 1.21 TRefArray* sorted = GetTracks(); //TEMPORANEO
1137     if(!sorted)return 0;
1138     TrkTrack *track = (TrkTrack*)sorted->At(it);
1139 pam-fi 1.8
1140 pam-fi 1.21 if(!track->HasImage()){
1141     cout << "** TrkLevel2 ** Track "<< it << "does not have image! " << endl;
1142     return 0;
1143     }
1144     if(!Track)return 0;
1145     TrkTrack *image = (TrkTrack*)(*Track)[track->image];
1146    
1147     sorted->Delete();
1148     delete sorted;
1149 pam-fi 1.8
1150 pam-fi 1.21 return image;
1151 pam-fi 1.8
1152 mocchiut 1.1 }
1153     //--------------------------------------
1154     //
1155     //
1156     //--------------------------------------
1157     /**
1158     * Loads the magnetic field.
1159     * @param s Path of the magnetic-field files.
1160     */
1161 pam-fi 1.16 void TrkLevel2::LoadField(TString path){
1162     //
1163 pam-fi 1.26 // strcpy(path_.path,path.Data());
1164     // path_.pathlen = path.Length();
1165     // path_.error = 0;
1166     // readb_();
1167    
1168     TrkParams::Set(path,1);
1169 pam-fi 1.28 TrkParams::Load(1);
1170 pam-fi 1.26
1171 pam-fi 1.16 //
1172 mocchiut 1.1 };
1173 pam-fi 1.25 /**
1174     * Get BY (kGauss)
1175     * @param v (x,y,z) coordinates in cm
1176     */
1177     float TrkLevel2::GetBX(float* v){
1178     float b[3];
1179     gufld_(v,b);
1180     return b[0]/10.;
1181     }
1182     /**
1183     * Get BY (kGauss)
1184     * @param v (x,y,z) coordinates in cm
1185     */
1186     float TrkLevel2::GetBY(float* v){
1187     float b[3];
1188     gufld_(v,b);
1189     return b[1]/10.;
1190     }
1191     /**
1192     * Get BY (kGauss)
1193     * @param v (x,y,z) coordinates in cm
1194     */
1195     float TrkLevel2::GetBZ(float* v){
1196     float b[3];
1197     gufld_(v,b);
1198     return b[2]/10.;
1199     }
1200 mocchiut 1.1 //--------------------------------------
1201     //
1202     //
1203     //--------------------------------------
1204     /**
1205 pam-fi 1.6 * Get tracker-plane (mechanical) z-coordinate
1206     * @param plane_id plane index (1=TOP,2,3,4,5,6=BOTTOM)
1207     */
1208     Float_t TrkLevel2::GetZTrk(Int_t plane_id){
1209     switch(plane_id){
1210     case 1: return ZTRK1;
1211     case 2: return ZTRK2;
1212     case 3: return ZTRK3;
1213     case 4: return ZTRK4;
1214     case 5: return ZTRK5;
1215     case 6: return ZTRK6;
1216     default: return 0.;
1217     };
1218     };
1219     //--------------------------------------
1220     //
1221     //
1222     //--------------------------------------
1223     /**
1224 pam-fi 1.2 * Trajectory default constructor.
1225     * (By default is created with z-coordinates inside the tracking volume)
1226     */
1227     Trajectory::Trajectory(){
1228     npoint = 10;
1229     x = new float[npoint];
1230     y = new float[npoint];
1231     z = new float[npoint];
1232     thx = new float[npoint];
1233     thy = new float[npoint];
1234     tl = new float[npoint];
1235 pam-fi 1.6 float dz = ((ZTRK1)-(ZTRK6))/(npoint-1);
1236 pam-fi 1.2 for(int i=0; i<npoint; i++){
1237     x[i] = 0;
1238     y[i] = 0;
1239 pam-fi 1.6 z[i] = (ZTRK1) - i*dz;
1240 pam-fi 1.2 thx[i] = 0;
1241     thy[i] = 0;
1242     tl[i] = 0;
1243     }
1244     }
1245     //--------------------------------------
1246     //
1247     //
1248     //--------------------------------------
1249     /**
1250 mocchiut 1.1 * Trajectory constructor.
1251 pam-fi 1.2 * (By default is created with z-coordinates inside the tracking volume)
1252 mocchiut 1.1 * \param n Number of points
1253     */
1254     Trajectory::Trajectory(int n){
1255 pam-fi 1.2 if(n<=0){
1256     cout << "NB! Trajectory must have at least 1 point >>> created with 10 points" << endl;
1257     n=10;
1258     }
1259 mocchiut 1.1 npoint = n;
1260     x = new float[npoint];
1261     y = new float[npoint];
1262     z = new float[npoint];
1263 pam-fi 1.2 thx = new float[npoint];
1264     thy = new float[npoint];
1265     tl = new float[npoint];
1266 pam-fi 1.6 float dz = ((ZTRK1)-(ZTRK6))/(npoint-1);
1267 mocchiut 1.1 for(int i=0; i<npoint; i++){
1268 pam-fi 1.2 x[i] = 0;
1269 mocchiut 1.1 y[i] = 0;
1270 pam-fi 1.6 z[i] = (ZTRK1) - i*dz;
1271 pam-fi 1.2 thx[i] = 0;
1272     thy[i] = 0;
1273     tl[i] = 0;
1274 mocchiut 1.1 }
1275     }
1276     //--------------------------------------
1277     //
1278     //
1279     //--------------------------------------
1280     /**
1281     * Trajectory constructor.
1282     * \param n Number of points
1283     * \param pz Pointer to float array, defining z coordinates
1284     */
1285     Trajectory::Trajectory(int n, float* zin){
1286 pam-fi 1.2 npoint = 10;
1287     if(n>0)npoint = n;
1288 mocchiut 1.1 x = new float[npoint];
1289     y = new float[npoint];
1290     z = new float[npoint];
1291 pam-fi 1.2 thx = new float[npoint];
1292     thy = new float[npoint];
1293     tl = new float[npoint];
1294     int i=0;
1295     do{
1296 pam-fi 1.21 x[i] = 0;
1297     y[i] = 0;
1298     z[i] = zin[i];
1299     thx[i] = 0;
1300     thy[i] = 0;
1301     tl[i] = 0;
1302     i++;
1303 pam-fi 1.2 }while(zin[i-1] > zin[i] && i < npoint);
1304     npoint=i;
1305     if(npoint != n)cout << "NB! Trajectory created with "<<npoint<<" points"<<endl;
1306 mocchiut 1.1 }
1307 pam-fi 1.21 void Trajectory::Delete(){
1308    
1309     if(x) delete [] x;
1310     if(y) delete [] y;
1311     if(z) delete [] z;
1312     if(thx) delete [] thx;
1313     if(thy) delete [] thy;
1314     if(tl) delete [] tl;
1315    
1316     }
1317 mocchiut 1.1 //--------------------------------------
1318     //
1319     //
1320     //--------------------------------------
1321     /**
1322     * Dump the trajectory coordinates.
1323     */
1324     void Trajectory::Dump(){
1325     cout <<endl<< "Trajectory ========== "<<endl;
1326     for (int i=0; i<npoint; i++){
1327 pam-fi 1.2 cout << i <<" >> " << x[i] <<" "<< y[i] <<" "<< z[i] ;
1328     cout <<" -- " << thx[i] <<" "<< thy[i] ;
1329     cout <<" -- " << tl[i] << endl;
1330 mocchiut 1.1 };
1331     }
1332 pam-fi 1.2 //--------------------------------------
1333     //
1334     //
1335     //--------------------------------------
1336     /**
1337     * Get trajectory length between two points
1338     * @param ifirst first point (default 0)
1339     * @param ilast last point (default npoint)
1340     */
1341     float Trajectory::GetLength(int ifirst, int ilast){
1342     if( ifirst<0 ) ifirst = 0;
1343     if( ilast>=npoint) ilast = npoint-1;
1344     float l=0;
1345     for(int i=ifirst;i<=ilast;i++){
1346     l=l+tl[i];
1347     };
1348     if(z[ilast] > ZINI)l=l-tl[ilast];
1349     if(z[ifirst] < ZINI) l=l-tl[ifirst];
1350    
1351     return l;
1352 mocchiut 1.1
1353 pam-fi 1.2 }
1354 pam-fi 1.6
1355 pam-fi 1.19 /**
1356     * Evaluates the trajectory in the apparatus associated to the track.
1357     * 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.
1358     * @param t pointer to an object of the class Trajectory,
1359     * which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ).
1360     * @return error flag.
1361     */
1362     int Trajectory::DoTrack2(float* al){
1363    
1364     double *dxout = new double[npoint];
1365     double *dyout = new double[npoint];
1366     double *dthxout = new double[npoint];
1367     double *dthyout = new double[npoint];
1368     double *dtlout = new double[npoint];
1369     double *dzin = new double[npoint];
1370     double dal[5];
1371    
1372     int ifail = 0;
1373    
1374     for (int i=0; i<5; i++) dal[i] = (double)al[i];
1375     for (int i=0; i<npoint; i++) dzin[i] = (double)z[i];
1376    
1377 pam-fi 1.26 TrkParams::Load(1);
1378     if( !TrkParams::IsLoaded(1) ){
1379     cout << "int Trajectory::DoTrack2(float* al) --- ERROR --- m.field not loaded"<<endl;
1380     return 0;
1381     }
1382 pam-fi 1.19 dotrack2_(&(npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail);
1383    
1384     for (int i=0; i<npoint; i++){
1385     x[i] = (float)*dxout++;
1386     y[i] = (float)*dyout++;
1387     thx[i] = (float)*dthxout++;
1388     thy[i] = (float)*dthyout++;
1389     tl[i] = (float)*dtlout++;
1390     }
1391    
1392     return ifail;
1393     };
1394 pam-fi 1.10
1395 mocchiut 1.1 ClassImp(TrkLevel2);
1396     ClassImp(TrkSinglet);
1397     ClassImp(TrkTrack);
1398     ClassImp(Trajectory);

  ViewVC Help
Powered by ViewVC 1.1.23