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

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

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

revision 1.12 by mocchiut, Tue Oct 24 07:28:38 2006 UTC revision 1.13 by pam-fi, Thu Oct 26 16:22:37 2006 UTC
# Line 12  extern "C" {     Line 12  extern "C" {    
12      void dotrack_(int*, double*, double*, double*, double*, int*);      void dotrack_(int*, double*, double*, double*, double*, int*);
13      void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);      void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);
14      int  readb_(const char*);      int  readb_(const char*);
15        void mini2_(int*,int*,int*);
16  }  }
17  //--------------------------------------  //--------------------------------------
18  //  //
# Line 196  Float_t TrkTrack::GetDEDX(){ Line 197  Float_t TrkTrack::GetDEDX(){
197  //--------------------------------------  //--------------------------------------
198  void TrkTrack::Dump(){  void TrkTrack::Dump(){
199      cout << endl << "========== Track " ;      cout << endl << "========== Track " ;
200          cout << endl << "seq.  n. : "<< seqno;      cout << endl << "seq.  n. : "<< seqno;
201          cout << endl << "image n. : "<< image;      cout << endl << "image n. : "<< image;
202          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] << " ";
203      cout << endl << "chi^2    : "<< chi2;      cout << endl << "chi^2    : "<< chi2;
204          cout << endl << "n.step   : "<< nstep;      cout << endl << "n.step   : "<< nstep;
205          cout << endl << "xgood    : "; for(int i=0; i<6; i++)cout << xgood[i] ;      cout << endl << "xgood    : "; for(int i=0; i<6; i++)cout << xgood[i] ;
206      cout << endl << "ygood    : "; for(int i=0; i<6; i++)cout << ygood[i] ;      cout << endl << "ygood    : "; for(int i=0; i<6; i++)cout << ygood[i] ;
207      cout << endl << "xm       : "; for(int i=0; i<6; i++)cout << xm[i] << " ";      cout << endl << "xm       : "; for(int i=0; i<6; i++)cout << xm[i] << " ";
208      cout << endl << "ym       : "; for(int i=0; i<6; i++)cout << ym[i] << " ";      cout << endl << "ym       : "; for(int i=0; i<6; i++)cout << ym[i] << " ";
209      cout << endl << "zm       : "; for(int i=0; i<6; i++)cout << zm[i] << " ";      cout << endl << "zm       : "; for(int i=0; i<6; i++)cout << zm[i] << " ";
210        cout << endl << "xv       : "; for(int i=0; i<6; i++)cout << xv[i] << " ";
211        cout << endl << "yv       : "; for(int i=0; i<6; i++)cout << yv[i] << " ";
212        cout << endl << "zv       : "; for(int i=0; i<6; i++)cout << zv[i] << " ";
213        cout << endl << "resx     : "; for(int i=0; i<6; i++)cout << resx[i] << " ";
214        cout << endl << "resy     : "; for(int i=0; i<6; i++)cout << resy[i] << " ";
215        cout << endl << "coval    : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" ";
216        cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[1][i]<<" ";
217        cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[2][i]<<" ";
218        cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[3][i]<<" ";
219        cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[4][i]<<" ";
220      cout << endl << "dedx_x   : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " ";      cout << endl << "dedx_x   : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " ";
221      cout << endl << "dedx_y   : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " ";      cout << endl << "dedx_y   : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " ";
222  }  }
223    /**
224     * Set the TrkTrack position measurements
225     */
226    void TrkTrack::SetMeasure(double *xmeas, double *ymeas, double *zmeas){
227        for(int i=0; i<6; i++) xm[i]=*xmeas++;
228        for(int i=0; i<6; i++) ym[i]=*ymeas++;
229        for(int i=0; i<6; i++) zm[i]=*zmeas++;
230    }
231    /**
232     * Set the TrkTrack position resolution
233     */
234    void TrkTrack::SetResolution(double *rx, double *ry){
235        for(int i=0; i<6; i++) resx[i]=*rx++;
236        for(int i=0; i<6; i++) resy[i]=*ry++;
237    }
238    /**
239     * Set the TrkTrack good measurement
240     */
241    void TrkTrack::SetGood(int *xg, int *yg){
242        for(int i=0; i<6; i++) xgood[i]=*xg++;
243        for(int i=0; i<6; i++) ygood[i]=*yg++;
244    }
245    
246    /**
247     * Load the magnetic field
248     */
249    void TrkTrack::LoadField(TString s){
250        readb_(s.Data());
251    };
252    /**
253     * Tracking method. It calls F77 mini routine.
254     */
255    void TrkTrack::Fit(double pfixed, int& fail, int iprint){
256        extern cMini2track track_;
257        fail = 0;
258    //    extern cMini2fitinfo fit_info_;
259    //    extern void mini_2_(int*,int*);
260    //    track_.xm[0]=1.0;
261        for(int i=0; i<6; i++) track_.xm[i]=xm[i];
262        for(int i=0; i<6; i++) track_.ym[i]=ym[i];
263        for(int i=0; i<6; i++) track_.zm[i]=zm[i];
264        for(int i=0; i<6; i++) track_.resx[i]=resx[i];
265        for(int i=0; i<6; i++) track_.resy[i]=resy[i];
266        for(int i=0; i<6; i++) track_.xgood[i]=xgood[i];
267        for(int i=0; i<6; i++) track_.ygood[i]=ygood[i];
268    // initial guess of "al" with linear fit
269        if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.){
270            double szz=0., szx=0., szy=0., ssx=0., ssy=0., sz=0., s1=0.;
271            double det, ax, ay, bx, by;
272            for(int i=0; i<NPLANE; i++) {
273                szz=szz+zm[i]*zm[i];
274                szx=szx+zm[i]*xm[i];
275                szy=szy+zm[i]*ym[i];
276                ssx=ssx+xm[i];
277                ssy=ssy+ym[i];
278                sz=sz+zm[i];
279                s1=s1+1.;
280            }
281            det=szz*s1-sz*sz;
282            ax=(szx*s1-sz*ssx)/det;
283            bx=(szz*ssx-szx*sz)/det;
284            ay=(szy*s1-sz*ssy)/det;
285            by=(szz*ssy-szy*sz)/det;
286            al[0]=ax*23.5+bx; // ZINI = 23.5 !!! it should be the same parameter in all codes
287            al[1]=ay*23.5+by; //  "  
288            al[2]= sqrt(pow(ax,2)+pow(ay,2))/ sqrt(pow(ax,2)+pow(ay,2)+1.);
289            al[3]=0.;
290            if( (ax!=0.)||(ay!=0.) ) {
291                al[3]= asin(ay/ sqrt(pow(ax,2)+pow(ay,2)));
292                if(ax<0.) al[3]=acos(-1.)-al[3];
293            }
294            al[4]=0.;
295        }
296    // end guess
297        for(int i=0; i<5; i++) track_.al[i]=al[i];
298        if(pfixed!=0.) {
299            al[4]=1./pfixed;      // to fix the momentum
300            track_.pfixed=pfixed; //         "
301        }
302        int istep=0;
303        int ifail=0;
304        mini2_(&istep,&ifail, &iprint);
305        if(ifail!=0) {
306            if(iprint==1)cout << "ERROR: ifail= " << ifail << endl;
307            fail = 1;
308    //      return;
309        }
310    
311    //    cout << endl << "eta ===> " << track_.al[4] << endl;
312    
313        for(int i=0; i<5; i++) al[i]=track_.al[i];
314        chi2=track_.chi2;
315        nstep=track_.nstep;
316        for(int i=0; i<6; i++) xv[i]=track_.xv[i];
317        for(int i=0; i<6; i++) yv[i]=track_.yv[i];
318        for(int i=0; i<6; i++) zv[i]=track_.zv[i];
319        for(int i=0; i<6; i++) axv[i]=track_.axv[i];
320        for(int i=0; i<6; i++) ayv[i]=track_.ayv[i];
321        for(int i=0; i<5; i++) {
322            for(int j=0; j<5; j++) coval[i][j]=track_.cov[i][j];
323        }
324    };
325    /*
326     * Reset the parameters fit
327     */
328    void TrkTrack::FitReset(){
329        for(int i=0; i<5; i++) al[i]=-9999.;
330        chi2=0.;
331        nstep=0;
332        for(int i=0; i<6; i++) xv[i]=0.;
333        for(int i=0; i<6; i++) yv[i]=0.;
334        for(int i=0; i<6; i++) zv[i]=0.;
335        for(int i=0; i<6; i++) axv[i]=0.;
336        for(int i=0; i<6; i++) ayv[i]=0.;
337        for(int i=0; i<5; i++) {
338            for(int j=0; j<5; j++) coval[i][j]=0.;
339        }
340    }
341    
342  //--------------------------------------  //--------------------------------------
343  //  //
344  //  //
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.23