/[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.22 - (hide annotations) (download)
Mon Jan 15 12:49:05 2007 UTC (17 years, 11 months ago) by pam-fi
Branch: MAIN
Changes since 1.21: +4 -2 lines
bug fixed

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

  ViewVC Help
Powered by ViewVC 1.1.23