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

Diff of /DarthVader/TrackerLevel2/src/ExtTrack.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by pam-fi, Thu Feb 27 11:24:43 2014 UTC revision 1.3 by mocchiut, Fri Jun 6 11:24:00 2014 UTC
# Line 6  Line 6 
6  #include <iostream>  #include <iostream>
7  #include <math.h>  #include <math.h>
8  using namespace std;  using namespace std;
 //......................................  
 // F77 routines  
 //......................................  
 extern "C" {      
 //     void dotrack_(int*, double*, double*, double*, double*, int*);  
 //     void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);  
 //     void dotrack3_(int*, double*, double*, double*, double*,double*, double*, double*,double*,int*);  
 //     void mini2_(int*,int*,int*);  
 //     void guess_();  
 //     void gufld_(float*, float*);  
 //     float risxeta2_(float *);  
 //     float risxeta3_(float *);  
 //     float risxeta4_(float *);  
 //     float risyeta2_(float *);  
 }  
9    
10    ExtTrack::ExtTrack(int dim){
11    
12      //cout << " ExtTrack::ExtTrack("<<dim<<")"<<endl;
13      xgood   = NULL;
14      ygood   = NULL;
15      multmaxx   = NULL;
16      multmaxy   = NULL;
17      xm      = NULL;
18      ym      = NULL;
19      zm      = NULL;
20      xma     = NULL;
21      yma     = NULL;
22      zma     = NULL;
23      xmb     = NULL;
24      ymb     = NULL;
25      zmb     = NULL;
26      resx    = NULL;
27      resy    = NULL;
28      xv      = NULL;
29      yv      = NULL;
30      zv      = NULL;
31      axv     = NULL;
32      ayv     = NULL;
33      dedx_x  = NULL;
34      dedx_y  = NULL;
35    
36      // xGF = new float[TrkParams::nGF];
37      // yGF = new float[TrkParams::nGF];
38      // init vectors  
39      SetDimension(dim);//allocate vectors
40      // init fit variables
41      Reset();
42      // set fit parameters (default)
43      SetZ0(23.5);        
44      SetMiniDefault();  
45    };
46  //--------------------------------------  //--------------------------------------
47  //  //
48  //  //
49  //--------------------------------------  //--------------------------------------
50  ExtTrack::ExtTrack(){  ExtTrack::ExtTrack( const ExtTrack& t ){
51      chi2  = 0;    //cout << " ExtTrack::ExtTrack( const ExtTrack& t ) "<<this<<endl;
     nstep = 0;  
     for(int it1=0;it1<5;it1++){  
         al[it1] = 0;  
         for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;  
     };  
     nplanes=0;  
   
     xgood   = 0;  
     ygood   = 0;  
     xm      = 0;  
     ym      = 0;  
     zm      = 0;  
     resx    = 0;  
     resy    = 0;  
     xv      = 0;  
     yv      = 0;  
     zv      = 0;  
     axv     = 0;  
     ayv     = 0;  
     dedx_x  = 0;  
     dedx_y  = 0;  
   
     TrkParams::SetMiniDefault();  
52    
53      // --------------------
54      // initialize pointers
55      // --------------------
56    
57      xgood   = NULL;
58      ygood   = NULL;
59      multmaxx   = NULL;
60      multmaxy   = NULL;
61      xm      = NULL;
62      ym      = NULL;
63      zm      = NULL;
64      xma     = NULL;
65      yma     = NULL;
66      zma     = NULL;
67      xmb     = NULL;
68      ymb     = NULL;
69      zmb     = NULL;
70      resx    = NULL;
71      resy    = NULL;
72      xv      = NULL;
73      yv      = NULL;
74      zv      = NULL;
75      axv     = NULL;
76      ayv     = NULL;
77      dedx_x  = NULL;
78      dedx_y  = NULL;
79    
80      // xGF = new float[TrkParams::nGF];
81      // yGF = new float[TrkParams::nGF];
82      
83      SetDimension(t.nplanes);//allocate vectors
84    
85      // --------------------
86      // copy values
87      // --------------------
88      zini  = t.zini;
89      chi2  = t.chi2;
90      nstep = t.nstep;
91      for(int it1=0;it1<5;it1++){
92        al[it1] = t.al[it1];
93        for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
94      };
95    
96    
97      for(int ip=0;ip<nplanes;ip++){
98        xgood[ip]  = t.xgood[ip];
99        ygood[ip]  = t.ygood[ip];
100        multmaxx[ip]  = t.multmaxx[ip];
101        multmaxy[ip]  = t.multmaxy[ip];
102        xm[ip]     = t.xm[ip];
103        ym[ip]     = t.ym[ip];
104        zm[ip]     = t.zm[ip];
105        xma[ip]     = t.xma[ip];
106        yma[ip]     = t.yma[ip];
107        zma[ip]     = t.zma[ip];
108        xmb[ip]     = t.xmb[ip];
109        ymb[ip]     = t.ymb[ip];
110        zmb[ip]     = t.zmb[ip];
111        resx[ip]   = t.resx[ip];
112        resy[ip]   = t.resy[ip];
113        xv[ip]     = t.xv[ip];
114        yv[ip]     = t.yv[ip];
115        zv[ip]     = t.zv[ip];
116        axv[ip]    = t.axv[ip];
117        ayv[ip]    = t.ayv[ip];
118        dedx_x[ip] = t.dedx_x[ip];
119        dedx_y[ip] = t.dedx_y[ip];
120      };
121      
122      int ngf = TrkParams::nGF;
123      for(int i=0; i<ngf; i++){
124        xGF[i] = t.xGF[i];
125        yGF[i] = t.yGF[i];
126      }
127      
128  };  };
129  //--------------------------------------  //--------------------------------------
130  //  //
131  //  //
132  //--------------------------------------  //--------------------------------------
133  ExtTrack::ExtTrack(Int_t dim){  void ExtTrack::Copy( ExtTrack& t ){
134      chi2  = 0;  
135      nstep = 0;    //  cout << " ExtTrack::Copy( ExtTrack& t ) "<<endl;
     for(int it1=0;it1<5;it1++){  
         al[it1] = 0;  
         for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;  
     };  
       
     SetDimension(dim);  
136    
137      TrkParams::SetMiniDefault();    t.Reset();//reset variables
138      t.Clear();//deallocate vectors
139      t.SetDimension(nplanes);//allocate vectors
140    
141      t.zini  = zini;
142      t.chi2  = chi2;
143      t.nstep = nstep;
144      for(int it1=0;it1<5;it1++){
145        t.al[it1] = al[it1];
146        for(int it2=0;it2<5;it2++)t.coval[it1][it2] = coval[it1][it2];
147      };    
148      //    cout << "ExtTrack::Copy()........nplanes "<<nplanes<<endl;
149      for(int ip=0;ip<nplanes;ip++){
150        t.xgood[ip]  = xgood[ip];
151        t.ygood[ip]  = ygood[ip];
152        t.multmaxx[ip]  = multmaxx[ip];
153        t.multmaxy[ip]  = multmaxy[ip];
154        t.xm[ip]     = xm[ip];
155        t.ym[ip]     = ym[ip];
156        t.zm[ip]     = zm[ip];
157        t.xma[ip]     = xma[ip];
158        t.yma[ip]     = yma[ip];
159        t.zma[ip]     = zma[ip];
160        t.xmb[ip]     = xmb[ip];
161        t.ymb[ip]     = ymb[ip];
162        t.zmb[ip]     = zmb[ip];
163        t.resx[ip]   = resx[ip];
164        t.resy[ip]   = resy[ip];
165        t.xv[ip]     = xv[ip];
166        t.yv[ip]     = yv[ip];
167        t.zv[ip]     = zv[ip];
168        t.axv[ip]    = axv[ip];
169        t.ayv[ip]    = ayv[ip];
170        t.dedx_x[ip] = dedx_x[ip];
171        t.dedx_y[ip] = dedx_y[ip];
172      };
173      int ngf = TrkParams::nGF;
174      for(int i=0; i<ngf; i++){
175        t.xGF[i] = xGF[i];
176        t.yGF[i] = yGF[i];
177      }
178      
179  };  };
180  //--------------------------------------  //--------------------------------------
181  //  //
182  //  //
183  //--------------------------------------  //--------------------------------------
184  ExtTrack::ExtTrack(const ExtTrack& t){  void ExtTrack::Set( TrkTrack& t , int index){
185      chi2  = t.chi2;  
186      nstep = t.nstep;    //  cout << "Set("<< &t <<", "<<index<<")"<<endl;
187      for(int it1=0;it1<5;it1++){    if(index < 0 ){ //NON FUNZIONA
188          al[it1] = t.al[it1];      cout << "Clear() "<<endl;
189          for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];      Reset();//reset variables
190      };      Clear();//deallocate vectors
191      SetDimension(t.nplanes);      cout << "SetDimension("<<t.GetNhit()<<")"<<endl;
192      for(int ip=0;ip<nplanes;ip++){      SetDimension(t.GetNhit()); //allocate vectors
193          xgood[ip]  = t.xgood[ip];    }
194          ygood[ip]  = t.ygood[ip];    // default fit settings (the same as TrkTrack)
195          xm[ip]     = t.xm[ip];    SetMiniDefault();
196          ym[ip]     = t.ym[ip];    SetZ0(23.5);
197          zm[ip]     = t.zm[ip];    //
198          resx[ip]   = t.resx[ip];    int ipp;
199          resy[ip]   = t.resy[ip];    if(index<0) ipp = 0;
200          xv[ip]     = t.xv[ip];    else        ipp = index;
201          yv[ip]     = t.yv[ip];    //
202          zv[ip]     = t.zv[ip];  //  cout << "Assign... nplanes = "<<nplanes<<endl;
203          axv[ip]    = t.axv[ip];    for(int ip=0;ip<6;ip++){
204          ayv[ip]    = t.ayv[ip];      
205          dedx_x[ip] = t.dedx_x[ip];      if(ipp>=nplanes){
206          dedx_y[ip] = t.dedx_y[ip];        cout << "void ExtTrack::Set( TrkTrack& , "<<index<<") -- ipp="<<ipp<<" exceed vector dimention "<<nplanes<<endl;
207      };        continue;
208        }
209      TrkParams::SetMiniDefault();  
210        if( t.XGood(ip) || t.YGood(ip) || !(index<0) ){
211    
212    //      cout << "ipp "<<ipp<<" ip "<<ip<<endl;
213    
214    
215          xgood[ipp]  = t.XGood(ip); //NB!
216          ygood[ipp]  = t.YGood(ip); //NB!
217          multmaxx[ipp]  = t.multmaxx[ip];
218          multmaxy[ipp]  = t.multmaxy[ip];
219          xm[ipp]     = t.xm[ip];
220          ym[ipp]     = t.ym[ip];
221          zm[ipp]     = t.zm[ip];
222    
223          float dummy = -100.;
224          xma[ipp]     = dummy;
225          yma[ipp]     = dummy;
226          zma[ipp]     = zm[ipp];
227          xmb[ipp]     = dummy;
228          ymb[ipp]     = dummy;
229          zmb[ipp]     = zm[ipp];
230    
231          if( t.XGood(ip)*t.YGood(ip) ){ //double hit
232            xma[ipp]     = t.xm[ip];
233            yma[ipp]     = t.ym[ip];
234            zma[ipp]     = t.zm[ip];
235            xmb[ipp]     = t.xm[ip];
236            ymb[ipp]     = t.ym[ip];
237            zmb[ipp]     = t.zm[ip];
238          } else if (t.XGood(ip) || t.YGood(ip) ){                         //single hit
239            TrkParams::Load(5);
240            if( !TrkParams::IsLoaded(5) ){
241              cout << "void ExtTrack::FillMiniStruct(TrkTrack&) ---ERROR--- trk align.param. not loaded "<<endl;
242              return;
243            }
244            double segment = 7.;// 2.;//cm //Elena 10th
245            int is = (int)t.GetSensor(ip); if(ip==5)is=1-is;
246            int ippp = 5-ip;
247            int il = (int)t.GetLadder(ip);      
248            // NB: i parametri di allineamento hanno una notazione particolare!!!
249            // sensor = 0 (hybrid side), 1
250            // ladder = 0-2 (increasing x)
251            // plane  = 0-5 (from bottom to top!!!)
252            double omega   = 0.;
253            double beta   = 0.;
254            double gamma   = 0.;
255            if(
256               (is < 0 || is > 1 || ip < 0 || ippp > 5 || il < 0 || il > 2) &&
257               true){
258              cout << " void ExtTrack::FillMiniStruct(TrkTrack&) --- WARNING ---trk sensor not defined, cannot read alignment parameters "<<endl;
259              cout << " is ip il = "<<is<<" "<<ippp<<" "<<il<<endl;
260            }else{
261              omega   = alignparameters_.omega[is][il][ippp];
262              beta    = alignparameters_.beta[is][il][ippp];
263              gamma   = alignparameters_.gamma[is][il][ippp];
264            }        
265            if(       t.XGood(ip) && !t.YGood(ip) ){
266              xma[ipp] = t.xm[ip] - omega * segment;
267              yma[ipp] = t.ym[ip] + segment;
268              zma[ipp] = t.zm[ip] + beta * segment;//not used yet
269              xmb[ipp] = t.xm[ip] + omega * segment;
270              ymb[ipp] = t.ym[ip] - segment;
271              zmb[ipp] = t.zm[ip] - beta * segment;//not used yet
272            }else if( !t.XGood(ip) && t.YGood(ip) ){
273              xma[ipp] = t.xm[ip] + segment;
274              yma[ipp] = t.ym[ip] + omega * segment;
275              zma[ipp] = t.zm[ip] - gamma * segment;//not used yet
276              xmb[ipp] = t.xm[ip] - segment;
277              ymb[ipp] = t.ym[ip] - omega * segment;
278              zmb[ipp] = t.zm[ip] + gamma * segment;//not used yet
279            }
280          }
281            
282          // tailx[ipp]  = t.tailx[ip];
283          // taily[ipp]  = t.taily[ip];
284          resx[ipp]   = t.resx[ip];
285          resy[ipp]   = t.resy[ip];
286          xv[ipp]     = t.xv[ip];
287          yv[ipp]     = t.yv[ip];
288          zv[ipp]     = t.zv[ip];
289          axv[ipp]    = t.axv[ip];
290          ayv[ipp]    = t.ayv[ip];
291          dedx_x[ipp] = t.dedx_x[ip];
292          dedx_y[ipp] = t.dedx_y[ip];
293    
294          ipp++;
295    
296        }
297      };//endl loop on planes
298    
299      chi2  = t.chi2;
300      nstep = t.nstep;
301      for(int it1=0;it1<5;it1++){
302        al[it1] = t.al[it1];
303        for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
304      };
305    
306  };  };
307  //--------------------------------------  //--------------------------------------
308  //  //
309  //  //
310  //--------------------------------------  //--------------------------------------
311  void ExtTrack::SetDimension(Int_t dim){  void ExtTrack::SetMiniDefault(){
312      //  cout << " ExtTrack::SetMiniDefault()"<<endl;
313      nplanes = dim;    exttrack_.trackmode = TrkParams::init__mini_trackmode;
314      exttrack_.fact      = TrkParams::init__mini_fact;  
315      xgood   = new Int_t[nplanes];    exttrack_.istepmin  = TrkParams::init__mini_istepmin;
316      ygood   = new Int_t[nplanes];    TrkParams::SetDeltaB();
317      xm      = new Float_t[nplanes];    TrkParams::SetDLT();
     ym      = new Float_t[nplanes];  
     zm      = new Float_t[nplanes];  
     resx    = new Float_t[nplanes];  
     resy    = new Float_t[nplanes];  
     xv      = new Float_t[nplanes];  
     yv      = new Float_t[nplanes];  
     zv      = new Float_t[nplanes];  
     axv     = new Float_t[nplanes];  
     ayv     = new Float_t[nplanes];  
     dedx_x  = new Float_t[nplanes];  
     dedx_y  = new Float_t[nplanes];  
       
     SetMeasure(0,0,0);  
     SetResolution(0,0);  
     SetGood(0,0);  
     FitReset();  
318  }  }
319  //--------------------------------------  //--------------------------------------
320  //  //
321  //  //
322  //--------------------------------------  //--------------------------------------
323  void ExtTrack::Clear(Option_t* option){  void ExtTrack::SetDimension(int dim){
   
     FitReset();  
     if(nplanes){  
         delete []       xgood;  
         delete []       ygood;  
         delete []       xm;  
         delete []       ym;  
         delete []       zm;  
         delete []       resx;  
         delete []       resy;  
         delete []       xv;  
         delete []       yv;  
         delete []       zv;  
         delete []       axv;  
         delete []       ayv;  
         delete []       dedx_x;  
         delete []       dedx_y;  
     }  
324    
325      //cout << " ExtTrack::SetDimension("<<dim<<")"<<endl;;
326      if(dim<=0){
327      nplanes = 0;      nplanes = 0;
328      xgood   = 0;      if(dim<0)cout << "void ExtTrack::SetDimension("<<dim<<") --- invalid dimension "<<endl;
329      ygood   = 0;      return;
330      xm      = 0;    }
331      ym      = 0;  
332      zm      = 0;    //  Clear();
333      resx    = 0;  
334      resy    = 0;    nplanes = dim;
335      xv      = 0;    
336      yv      = 0;    xgood   = new int[nplanes];
337      zv      = 0;    ygood   = new int[nplanes];
338      axv     = 0;    multmaxx   = new int[nplanes];
339      ayv     = 0;    multmaxy   = new int[nplanes];
340      dedx_x  = 0;    xm      = new float[nplanes];
341      dedx_y  = 0;    ym      = new float[nplanes];
342      zm      = new float[nplanes];
343      xma     = new float[nplanes];
344      yma     = new float[nplanes];
345      zma     = new float[nplanes];
346      xmb     = new float[nplanes];
347      ymb     = new float[nplanes];
348      zmb     = new float[nplanes];
349      resx    = new float[nplanes];
350      resy    = new float[nplanes];
351      xv      = new float[nplanes];
352      yv      = new float[nplanes];
353      zv      = new float[nplanes];
354      axv     = new float[nplanes];
355      ayv     = new float[nplanes];
356      dedx_x  = new float[nplanes];
357      dedx_y  = new float[nplanes];
358      
359    
360    }
361    //--------------------------------------
362    //
363    //
364    //--------------------------------------
365    bool ExtTrack::SetZ(int ip,float zmeas){
366      if(ip<0 || ip>=nplanes){
367        cout << " ExtTrack::SetZ("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
368        return false;
369      }
370      zm[ip] = zmeas;
371      zma[ip] = zmeas;
372      zmb[ip] = zmeas;
373      return true;
374    
375    };
376    
377    bool ExtTrack::ResetXY(int ip){
378    
379      // cout << " ExtTrack::ResetXY("<<ip<<") "<<endl;
380    
381      if(ip<0 || ip>=nplanes){
382        cout << " ExtTrack::ResetXY("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
383        return false;
384      }
385      xm[ip] = -100.;
386      ym[ip] = -100.;
387      xma[ip] = -100.;
388      yma[ip] = -100.;
389      xmb[ip] = -100.;
390      ymb[ip] = -100.;
391      resx[ip] = 1000.;
392      resy[ip] = 1000.;
393      dedx_x[ip] = 0.;
394      dedx_y[ip] = 0.;
395      xgood[ip] = 0;
396      ygood[ip] = 0;
397      multmaxx[ip] = 0;
398      multmaxy[ip] = 0;
399    
400      return true;
401    
402  };  };
403    
404    bool ExtTrack::SetXY(int ip,float xmeas, float ymeas, float rx, float ry){
405    
406      if(ip<0 || ip>=nplanes){
407        cout << " ExtTrack::SetXY("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
408        return false;
409      }
410      
411      xm[ip] = xmeas;
412      ym[ip] = ymeas;
413      xma[ip] = xmeas;
414      yma[ip] = ymeas;
415      xmb[ip] = xmeas;
416      ymb[ip] = ymeas;
417      resx[ip] = rx;
418      resy[ip] = ry;
419      // dedx_x[ip] = dedxx;
420      // dedx_y[ip] = dedxy;
421      xgood[ip] = 1;
422      ygood[ip] = 1;
423    
424      return true;
425    
426    };
427    bool ExtTrack::SetX(int ip,float xa, float xb, float ya, float yb,float res ){
428    
429      // cout << " ExtTrack::SetX("<<ip<<",...)  -- nplanes"<<nplanes<<endl;
430    
431      if(ip<0 || ip>=nplanes){
432        cout << " ExtTrack::SetX("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
433        return false;
434      }
435      xm[ip] = -100.;
436      ym[ip] = -100.;
437      xma[ip] = xa;
438      yma[ip] = ya;
439      xmb[ip] = xb;
440      ymb[ip] = yb;
441      resx[ip] = res;
442      resy[ip] = 1000.;
443      xgood[ip] = 1;
444      ygood[ip] = 0;
445      // dedx_x[ip] = dedx;
446      // dedx_y[ip] = 0.;
447      return true;
448    };
449    bool ExtTrack::SetY(int ip,float xa, float xb, float ya, float yb,float res ){
450    
451      // cout << " ExtTrack::SetY("<<ip<<",...)  -- nplanes"<<nplanes<<endl;
452    
453      if(ip<0 || ip>=nplanes){
454        cout << " ExtTrack::SetX("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
455        return false;
456      }
457      xm[ip] = -100.;
458      ym[ip] = -100.;
459      xma[ip] = xa;
460      yma[ip] = ya;
461      xmb[ip] = xb;
462      ymb[ip] = yb;
463      resx[ip] = 1000.;
464      resy[ip] = res;
465      xgood[ip] = 0;
466      ygood[ip] = 1;
467      // dedx_x[ip] = 0.;
468      // dedx_y[ip] = dedx;
469      return true;
470    };
471    
472    bool ExtTrack::SetXGood(int ip,int icl_piu_uno, int il, int is ){
473    
474      // cout << " ExtTrack::SetXGood("<<ip<<",...)  -- nplanes"<<nplanes<<endl;
475    
476      if(ip<0 || ip>=nplanes){
477        cout << " ExtTrack::SetXGood("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
478        return false;
479      }
480      xgood[ip]=(il+1)*100000000+(is+1)*10000000+icl_piu_uno;
481      return true;
482    };
483    bool ExtTrack::SetYGood(int ip,int icl_piu_uno, int il, int is ){
484      if(ip<0 || ip>=nplanes){
485        cout << " ExtTrack::SetYGood("<<ip<<",...)  -- nplanes"<<nplanes<<endl;
486        cout << " ExtTrack::SetYGood("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
487        return false;
488      }
489      ygood[ip]=(il+1)*100000000+(is+1)*10000000+icl_piu_uno;
490      return true;
491    };
492    
493  //--------------------------------------  //--------------------------------------
494  //  //
495  //  //
496  //--------------------------------------  //--------------------------------------
497  void ExtTrack::Delete(){  void ExtTrack::Clear(Option_t* option){
498      Clear();  
499      //  cout << " ExtTrack::Clear("<<option<<")  "<<this<<endl;
500    
501    
502      //  cout << " xgood "<<xgood<<endl;
503      
504      //  if(nplanes>0){
505      if(xgood) delete []   xgood;
506      if(ygood) delete []   ygood;
507      if(multmaxx) delete [] multmaxx;
508      if(multmaxy) delete [] multmaxy;
509      if(xm)    delete []   xm;
510      if(ym)    delete []   ym;
511      if(zm)    delete []   zm;
512      if(xma)   delete []   xma;
513      if(yma)   delete []   yma;
514      if(zma)   delete []   zma;
515      if(xmb)   delete []   xmb;
516      if(ymb)   delete []   ymb;
517      if(zmb)   delete []   zmb;
518      if(resx)  delete []   resx;
519      if(resy)  delete []   resy;
520      if(xv)    delete []   xv;
521      if(yv)    delete []   yv;
522      if(zv)    delete []   zv;
523      if(axv)   delete []   axv;
524      if(ayv)   delete []   ayv;
525      if(dedx_x)delete []   dedx_x;
526      if(dedx_y)delete []   dedx_y;
527      //}
528    
529      xgood   = NULL;
530      ygood   = NULL;
531      multmaxx   = NULL;
532      multmaxy   = NULL;
533      xm      = NULL;
534      ym      = NULL;
535      zm      = NULL;
536      xma     = NULL;
537      yma     = NULL;
538      zma     = NULL;
539      xmb     = NULL;
540      ymb     = NULL;
541      zmb     = NULL;
542      resx    = NULL;
543      resy    = NULL;
544      xv      = NULL;
545      yv      = NULL;
546      zv      = NULL;
547      axv     = NULL;
548      ayv     = NULL;
549      dedx_x  = NULL;
550      dedx_y  = NULL;
551    
552      nplanes = 0;
553    
554      //  Reset();
555    
556  };  };
557    
558  /**  void ExtTrack::Delete(){
559   * Set the  position measurements    Clear();
560   */    // delete [] xGF;
561  void ExtTrack::SetMeasure(Double_t *xmeas, Double_t *ymeas, Double_t *zmeas){    // delete [] yGF;
562      for(Int_t i=0; i<nplanes; i++) xm[i]= (xmeas? *xmeas++ : 0.);  }        
     for(Int_t i=0; i<nplanes; i++) ym[i]= (ymeas? *ymeas++ : 0.);  
     for(Int_t i=0; i<nplanes; i++) zm[i]= (ymeas? *zmeas++ : 0.);  
 }  
 /**  
  * Set the TrkTrack position resolution  
  */  
 void ExtTrack::SetResolution(Double_t *rx, Double_t *ry){  
     for(Int_t i=0; i<nplanes; i++) resx[i]= (rx? *rx++ : 0.);  
     for(Int_t i=0; i<nplanes; i++) resy[i]= (ry? *ry++ : 0.);  
 }  
 /**  
  * Set the TrkTrack good measurement  
  */  
 void ExtTrack::SetGood(Int_t *xg, Int_t *yg){  
563    
564      for(Int_t i=0; i<nplanes; i++) xgood[i]=(xg? *xg++ : 0 );  //--------------------------------------
565      for(Int_t i=0; i<nplanes; i++) ygood[i]=(yg? *yg++ : 0 );  //
566  }  //
567    //--------------------------------------
568    
569    // /**
570    //  * Set the  position measurements
571    //  */
572    // void ExtTrack::SetMeasure(Double_t *xmeas, Double_t *ymeas, Double_t *zmeas){
573    //     for(int i=0; i<nplanes; i++) xm[i]= (xmeas? *xmeas++ : 0.);
574    //     for(int i=0; i<nplanes; i++) ym[i]= (ymeas? *ymeas++ : 0.);
575    //     for(int i=0; i<nplanes; i++) zm[i]= (ymeas? *zmeas++ : 0.);
576    // }
577    // /**
578    //  * Set the TrkTrack position resolution
579    //  */
580    // void ExtTrack::SetResolution(Double_t *rx, Double_t *ry){
581    //     for(int i=0; i<nplanes; i++) resx[i]= (rx? *rx++ : 0.);
582    //     for(int i=0; i<nplanes; i++) resy[i]= (ry? *ry++ : 0.);
583    // }
584    // /**
585    //  * Set the TrkTrack good measurement
586    //  */
587    // void ExtTrack::SetGood(int *xg, int *yg){
588    
589    //     for(int i=0; i<nplanes; i++) xgood[i]=(xg? *xg++ : 0 );
590    //     for(int i=0; i<nplanes; i++) ygood[i]=(yg? *yg++ : 0 );
591    // }
592    
593  //--------------------------------------  //--------------------------------------
594  //  //
# Line 207  void ExtTrack::SetGood(Int_t *xg, Int_t Line 596  void ExtTrack::SetGood(Int_t *xg, Int_t
596  //--------------------------------------  //--------------------------------------
597  void ExtTrack::Dump(){  void ExtTrack::Dump(){
598      cout << endl << "========== Track " ;      cout << endl << "========== Track " ;
599        cout << endl << "zini     : "<< zini;
600      cout << endl << "al       : "; for(int i=0; i<5; i++)cout << al[i] << " ";      cout << endl << "al       : "; for(int i=0; i<5; i++)cout << al[i] << " ";
601      cout << endl << "chi^2    : "<< chi2;      cout << endl << "chi^2    : "<< chi2;
602      cout << endl << "n.step   : "<< nstep;      cout << endl << "n.step   : "<< nstep;
# Line 215  void ExtTrack::Dump(){ Line 605  void ExtTrack::Dump(){
605      cout << endl << "xm       : "; for(int i=0; i<nplanes; i++)cout << xm[i] << " ";      cout << endl << "xm       : "; for(int i=0; i<nplanes; i++)cout << xm[i] << " ";
606      cout << endl << "ym       : "; for(int i=0; i<nplanes; i++)cout << ym[i] << " ";      cout << endl << "ym       : "; for(int i=0; i<nplanes; i++)cout << ym[i] << " ";
607      cout << endl << "zm       : "; for(int i=0; i<nplanes; i++)cout << zm[i] << " ";      cout << endl << "zm       : "; for(int i=0; i<nplanes; i++)cout << zm[i] << " ";
608        cout << endl << "xma       : "; for(int i=0; i<nplanes; i++)cout << xma[i] << " ";
609        cout << endl << "yma       : "; for(int i=0; i<nplanes; i++)cout << yma[i] << " ";
610        cout << endl << "zma       : "; for(int i=0; i<nplanes; i++)cout << zma[i] << " ";
611        cout << endl << "xmb       : "; for(int i=0; i<nplanes; i++)cout << xmb[i] << " ";
612        cout << endl << "ymb       : "; for(int i=0; i<nplanes; i++)cout << ymb[i] << " ";
613        cout << endl << "zmb       : "; for(int i=0; i<nplanes; i++)cout << zmb[i] << " ";
614      cout << endl << "xv       : "; for(int i=0; i<nplanes; i++)cout << xv[i] << " ";      cout << endl << "xv       : "; for(int i=0; i<nplanes; i++)cout << xv[i] << " ";
615      cout << endl << "yv       : "; for(int i=0; i<nplanes; i++)cout << yv[i] << " ";      cout << endl << "yv       : "; for(int i=0; i<nplanes; i++)cout << yv[i] << " ";
616      cout << endl << "zv       : "; for(int i=0; i<nplanes; i++)cout << zv[i] << " ";      cout << endl << "zv       : "; for(int i=0; i<nplanes; i++)cout << zv[i] << " ";
# Line 225  void ExtTrack::Dump(){ Line 621  void ExtTrack::Dump(){
621      cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[2][i]<<" ";      cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[2][i]<<" ";
622      cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[3][i]<<" ";      cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[3][i]<<" ";
623      cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[4][i]<<" ";      cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[4][i]<<" ";
624        cout << endl << "MDR      : "<< (coval[4][4]>0. ?  1./sqrt(coval[4][4]) : 0. ) ;
625      cout << endl << "dedx_x   : "; for(int i=0; i<nplanes; i++)cout << dedx_x[i] << " ";      cout << endl << "dedx_x   : "; for(int i=0; i<nplanes; i++)cout << dedx_x[i] << " ";
626      cout << endl << "dedx_y   : "; for(int i=0; i<nplanes; i++)cout << dedx_y[i] << " ";      cout << endl << "dedx_y   : "; for(int i=0; i<nplanes; i++)cout << dedx_y[i] << " ";
627    
628        cout << endl << "=================";
629      cout << endl;      cout << endl;
630  }  }
631  /**  /**
632   * Reset the fit parameters   * Reset the fit parameters
633   */   */
634  void ExtTrack::FitReset(){  void ExtTrack::ResetFit(){
635    
636      // cout << "  ExtTrack::ResetFit() "<<endl;
637      for(int i=0; i<5; i++) al[i]=-9999.;      for(int i=0; i<5; i++) al[i]=-9999.;
638      chi2=0.;      chi2=0.;
639      nstep=0;      nstep=0;
640      for(int i=0; i<nplanes; i++) xv[i]=0.;  
     for(int i=0; i<nplanes; i++) yv[i]=0.;  
     for(int i=0; i<nplanes; i++) zv[i]=0.;  
     for(int i=0; i<nplanes; i++) axv[i]=0.;  
     for(int i=0; i<nplanes; i++) ayv[i]=0.;  
641      for(int i=0; i<5; i++) {      for(int i=0; i<5; i++) {
642          for(int j=0; j<5; j++) coval[i][j]=0.;          for(int j=0; j<5; j++) coval[i][j]=0.;
643      }      }
644    
645        for(int i=0; i<nplanes; i++){
646          xv[i] = 0.;
647          yv[i] = 0.;
648          zv[i] = 0.;
649          axv[i] = 0.;
650          ayv[i] = 0.;
651        }
652    
653        int ngf = TrkParams::nGF;
654        for(int i=0; i<ngf; i++){
655          xGF[i] = 0.;
656          yGF[i] = 0.;
657        }
658    
659    
660  }  }
661  /**  /**
662   * Method to fill minimization-routine common   * Reset x/y measured coordinates (NB it does not reset z coordinates)
  * NB (zini=23.5 is implicitelly set)  
663   */   */
664  void ExtTrack::FillMiniStruct(cMini2track& track){  void ExtTrack::ResetXY(){
665    
666      for(int i=0; i<nplanes; i++){    for(int i=0; i<nplanes; i++)ResetXY(i);
667          
668          track.xgood[i]=XGood(i);      
669          track.ygood[i]=YGood(i);  }
670            /**
671          track.xm[i]=xm[i];   * \brief Tracking method. It calls F77 mini routine.
672          track.ym[i]=ym[i];   *
673          track.zm[i]=zm[i];   * @see void TrkTrack::Fit(double, int&, int, int)
674             */
675          float segment = 100.; //temporaneo  void ExtTrack::Fit(double pfixed, int& fail, int iprint){
         track.xm_a[i]=xm[i];  
         track.xm_b[i]=xm[i];  
         track.ym_a[i]=ym[i];  
         track.ym_b[i]=ym[i];  
         if(       XGood(i) && !YGood(i) ){  
             track.ym_a[i] = track.ym_a[i]+segment;  
             track.ym_b[i] = track.ym_b[i]-segment;  
         }else if( !XGood(i) && YGood(i)){  
             track.xm_a[i] = track.xm_a[i]+segment;  
             track.xm_b[i] = track.xm_b[i]-segment;  
         }  
676    
677              TrkParams::Load(1);
678          track.resx[i]=resx[i];    if( !TrkParams::IsLoaded(1) ){
679          track.resy[i]=resy[i];      cout << "void ExtTrack::Fit(double,int&,int,int) ---ERROR--- m.field not loaded "<<endl;
680      }      return;
681      }
682    
683      float al_ini[] = {0.,0.,0.,0.,0.};
684      
685      fail = 0;
686      
687      FillMiniStruct();   //fill F77 common
688        
689    
690      for(int i=0; i<5; i++) track.al[i]=al[i];    // if fit variables have been reset, evaluate the initial guess
691      if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.){
692        if(iprint==2)cout << "calling guessext_()"<<endl;
693        guessext_();
694      }
695      // --------------------- free momentum
696      if(pfixed==0.) {
697        exttrack_.pfixed=0.;
698      }
699      // --------------------- fixed momentum
700      if(pfixed!=0.) {
701        al[4]=1./pfixed;    
702        exttrack_.pfixed=pfixed;
703      }
704      
705      //  store temporarily the initial guess
706      for(int i=0; i<5; i++) al_ini[i]=exttrack_.al[i];
707    
708      if(iprint==2){
709        cout<<endl<<" init. guess: ";
710        for(int i=0; i<5; i++)cout << al_ini[i]<<" ";
711        cout<<endl;
712      }
713      //  ------------------------------------------
714      //  call mini routine
715      //  ------------------------------------------
716      int istep=0;
717      int ifail=0;
718      if(iprint==2)cout << "calling miniext_("<<istep<<","<<ifail<<", "<<iprint<<")"<<endl;
719      miniext_(&istep,&ifail, &iprint);
720      if(ifail!=0) {
721        if(iprint)cout << "ERROR: ifail= " << ifail << endl;
722        fail = 1;
723      }
724      if(chi2!=chi2){
725        if(iprint)cout << "ERROR: chi2= " << chi2 << endl;  
726        ResetFit();
727        fail = 1;  
728      }
729      //  ------------------------------------------
730      
731      SetFromMiniStruct(); //copy content of F77 common into the ExtTrack object
732      
733      if(fail){
734        if(iprint)cout << " >>>> fit failed "<<endl;
735        for(int i=0; i<5; i++) al[i]=al_ini[i];
736      }
737    };
738    
739      track.zini = 55.;  /**
740     * Method to fill minimization-routine common
741     * NB (zini=23.5 is implicitelly set)
742     */
743    void ExtTrack::FillMiniStruct(cMiniExtTrack& track){
744    
745  //    track.zini = 23.5;    track.nplanes = nplanes;
746  // ZINI = 23.5 !!! it should be the same parameter in all codes    for(int i=0; i<nplanes; i++){
747        
748        track.xgood[i] = (double) XGood(i);
749        track.ygood[i] = (double) YGood(i);
750        
751        track.xm[i] = (double) xm[i];
752        track.ym[i] = (double) ym[i];
753        track.zm[i] = (double) zm[i];
754        
755        
756        track.xm_a[i] = (double) xma[i];
757        track.xm_b[i] = (double) xmb[i];
758        track.ym_a[i] = (double) yma[i];
759        track.ym_b[i] = (double) ymb[i];
760        track.zm_a[i] = (double) zma[i];
761        track.zm_b[i] = (double) zmb[i];
762        
763        track.resx[i] = (double) resx[i];
764        track.resy[i] = (double) resy[i];
765      }
766      
767      for(int i=0; i<5; i++) track.al[i] = (double) al[i];
768      
769      track.zini = (double) zini;
770      
771            
772  }  }
773  /**  /**
774   * Method to set values from  minimization-routine common   * Method to set values from  minimization-routine common
775   */   */
776  void ExtTrack::SetFromMiniStruct(cMini2track *track){  void ExtTrack::SetFromMiniStruct(cMiniExtTrack *track){
777    
778      for(int i=0; i<5; i++) {      for(int i=0; i<5; i++) {
779          al[i]=track->al[i];        al[i] = (float) (track->al[i]);
780          for(int j=0; j<5; j++) coval[i][j]=track->cov[i][j];        for(int j=0; j<5; j++) coval[i][j] = (float) (track->cov[i][j]);
781      }      }
782      chi2  = track->chi2;      chi2  = (float) (track->chi2);
783      nstep = track->nstep;      nstep = (float) (track->nstep);
784      for(int i=0; i<6; i++){      for(int i=0; i<nplanes; i++){
785          xv[i]  = track->xv[i];          xv[i]  = (float) (track->xv[i]);
786          yv[i]  = track->yv[i];          yv[i]  = (float) (track->yv[i]);
787          zv[i]  = track->zv[i];          zv[i]  = (float) (track->zv[i]);
788          xm[i]  = track->xm[i];          xm[i]  = (float) (track->xm[i]);
789          ym[i]  = track->ym[i];          ym[i]  = (float) (track->ym[i]);
790          zm[i]  = track->zm[i];          zm[i]  = (float) (track->zm[i]);
791          axv[i] = track->axv[i];          xma[i]  = (float) (track->xm_a[i]);
792          ayv[i] = track->ayv[i];          yma[i]  = (float) (track->ym_a[i]);
793            zma[i]  = (float) (track->zm_a[i]);
794            xmb[i]  = (float) (track->xm_b[i]);
795            ymb[i]  = (float) (track->ym_b[i]);
796            zmb[i]  = (float) (track->zm_b[i]);
797            axv[i] = (float) (track->axv[i]);
798            ayv[i] = (float) (track->ayv[i]);      
799      }      }
800            zini = (float) (track->zini);
801    }
802    /**
803     * Method to retrieve ID (0,1,...) of x-cluster (if any) associated to this track.
804     * If no cluster is associated, ID=-1.
805     *
806     */
807    int ExtTrack::GetClusterX_ID(int ip){
808        return ((int)fabs(xgood[ip]))%10000000-1;
809    };
810    /**
811     * Method to retrieve ID (0-xxx) of y-cluster (if any) associated to this track.
812     *
813     */
814    int ExtTrack::GetClusterY_ID(int ip){
815        return ((int)fabs(ygood[ip]))%10000000-1;
816    };
817    
818    /**
819     * Method to retrieve the dE/dx measured on a tracker view.
820     * @param ip plane (0-5)
821     * @param iv view (0=x 1=y)
822     */
823    Float_t ExtTrack::GetDEDX(int ip, int iv){
824        if(iv==0 && ip>=0 && ip<6)return fabs(dedx_x[ip]);
825        else if(iv==1 && ip>=0 && ip<6)return fabs(dedx_y[ip]);
826        else {
827            cout << "TrkTrack::GetDEDX(int ip, int iv) -- wrong input parameters "<<ip<<iv<<endl;
828            return 0.;
829        }
830    }
831    /**
832     * Method to evaluate the dE/dx measured on a tracker plane.
833     * The two measurements on x- and y-view are averaged.
834     * @param ip plane (0-5)
835     */
836    Float_t ExtTrack::GetDEDX(int ip){
837        if( (Int_t)XGood(ip)+(Int_t)YGood(ip) == 0 ) return 0;
838        return (GetDEDX(ip,0)+GetDEDX(ip,1))/((Int_t)XGood(ip)+(Int_t)YGood(ip));
839    };
840    
841    /**
842     * Method to evaluate the dE/dx averaged over all planes.
843     */
844    Float_t ExtTrack::GetDEDX(){
845        Float_t dedx=0;
846        for(Int_t ip=0; ip<6; ip++)dedx+=GetDEDX(ip,0)*XGood(ip)+GetDEDX(ip,1)*YGood(ip);
847        dedx = dedx/(GetNX()+GetNY());
848        return dedx;
849    };
850    
851    /**
852     * Method to retrieve the ladder, if assigned
853     */
854    int ExtTrack::GetLadder(int ip){
855        if(XGood(ip))return (int)fabs(xgood[ip]/100000000)-1;
856        if(YGood(ip))return (int)fabs(ygood[ip]/100000000)-1;
857        return -1;
858    };
859    /**
860     * Method to retrieve the sensor, if assigned
861     */
862    int ExtTrack::GetSensor(int ip){
863        if(XGood(ip))return (int)((int)fabs(xgood[ip]/10000000)%10)-1;
864        if(YGood(ip))return (int)((int)fabs(ygood[ip]/10000000)%10)-1;
865        return -1;
866    };
867    /**
868     *
869     */
870    Int_t ExtTrack::GetClusterX_Multiplicity(int ip){
871      if(ip<0 || ip>=nplanes){
872        cout << " ExtTrack::GetClusterX_Multiplicity("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
873        return -1;
874      }
875      return (Int_t)(multmaxx[ip]/10000);
876  }  }
877    /**
878     *
879     */
880    
881    Int_t ExtTrack::GetClusterY_Multiplicity(int ip){
882      if(ip<0 || ip>=nplanes){
883        cout << " ExtTrack::GetClusterY_Multiplicity("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
884        return -1;
885      }
886      return (Int_t)(multmaxy[ip]/10000);
887    }
888    /**
889     *
890     */
891    
892    Int_t ExtTrack::GetClusterX_MaxStrip(int ip){
893      if(ip<0 || ip>=nplanes){
894        cout << " ExtTrack::GetClusterX_MaxStrip("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
895        return -1;
896      }
897      return (Int_t)(multmaxx[ip]%10000);
898    }
899    /**
900     *
901     */
902    
903    Int_t ExtTrack::GetClusterY_MaxStrip(int ip){
904      if(ip<0 || ip>=nplanes){
905        cout << " ExtTrack::GetClusterY_MaxStrip("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
906        return -1;
907      }
908      return (Int_t)(multmaxy[ip]%10000);
909    }
910    Float_t ExtTrack::GetRigidity(){
911            Float_t rig=0;
912            if(chi2>0)rig=1./al[4];
913            if(rig<0) rig=-rig;
914            return rig;
915    };
916    //
917    Float_t ExtTrack::GetDeflection(){
918            Float_t def=0;
919            if(chi2>0)def=al[4];
920            return def;
921    };
922  ClassImp(ExtTrack);  ClassImp(ExtTrack);
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23