/[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.35 by pam-fi, Thu May 24 14:32:14 2007 UTC revision 1.40 by pam-fi, Fri Aug 31 14:56:51 2007 UTC
# Line 52  TrkTrack::TrkTrack(){ Line 52  TrkTrack::TrkTrack(){
52          ayv[ip]    = 0;          ayv[ip]    = 0;
53          dedx_x[ip] = 0;          dedx_x[ip] = 0;
54          dedx_y[ip] = 0;          dedx_y[ip] = 0;
55            multmaxx[ip] = 0;
56            multmaxy[ip] = 0;
57            seedx[ip]  = 0;  
58            seedy[ip]  = 0;
59            xpu[ip]    = 0;  
60            ypu[ip]    = 0;  
61    
62      };      };
 //    clx = 0;  
 //    cly = 0;  
 //    clx = new TRefArray(6,0); //forse causa memory leak???  
 //    cly = new TRefArray(6,0); //forse causa memory leak???  
 //    clx = TRefArray(6,0);  
 //    cly = TRefArray(6,0);  
63    
64      TrkParams::SetTrackingMode();      TrkParams::SetTrackingMode();
65      TrkParams::SetPrecisionFactor();      TrkParams::SetPrecisionFactor();
# Line 96  TrkTrack::TrkTrack(const TrkTrack& t){ Line 97  TrkTrack::TrkTrack(const TrkTrack& t){
97          ayv[ip]    = t.ayv[ip];          ayv[ip]    = t.ayv[ip];
98          dedx_x[ip] = t.dedx_x[ip];          dedx_x[ip] = t.dedx_x[ip];
99          dedx_y[ip] = t.dedx_y[ip];          dedx_y[ip] = t.dedx_y[ip];
100            multmaxx[ip] = t.multmaxx[ip];
101            multmaxy[ip] = t.multmaxy[ip];
102            seedx[ip]    = t.seedx[ip];  
103            seedy[ip]    = t.seedy[ip];
104            xpu[ip]      = t.xpu[ip];  
105            ypu[ip]      = t.ypu[ip];  
106      };      };
 //    clx = 0;  
 //    cly = 0;  
 //    if(t.clx)clx = new TRefArray(*(t.clx));  
 //    if(t.cly)cly = new TRefArray(*(t.cly));  
 //    clx = TRefArray(t.clx);  
 //    cly = TRefArray(t.cly);  
107    
108      TrkParams::SetTrackingMode();      TrkParams::SetTrackingMode();
109      TrkParams::SetPrecisionFactor();      TrkParams::SetPrecisionFactor();
# Line 141  void TrkTrack::Copy(TrkTrack& t){ Line 142  void TrkTrack::Copy(TrkTrack& t){
142          t.ayv[ip]    = ayv[ip];          t.ayv[ip]    = ayv[ip];
143          t.dedx_x[ip] = dedx_x[ip];          t.dedx_x[ip] = dedx_x[ip];
144          t.dedx_y[ip] = dedx_y[ip];          t.dedx_y[ip] = dedx_y[ip];
145            t.multmaxx[ip] = multmaxx[ip];
146            t.multmaxy[ip] = multmaxy[ip];
147            t.seedx[ip]    = seedx[ip];  
148            t.seedy[ip]    = seedy[ip];
149            t.xpu[ip]      = xpu[ip];  
150            t.ypu[ip]      = ypu[ip];  
151                            
152      };      };
153    
 //    t.clx = TRefArray(clx);  
 //    t.cly = TRefArray(cly);  
154            
155  };  };
156  //--------------------------------------  //--------------------------------------
# Line 365  Int_t TrkTrack::GetLeverArmY(){ Line 370  Int_t TrkTrack::GetLeverArmY(){
370      }      }
371      return (last_plane-first_plane+1);      return (last_plane-first_plane+1);
372  }  }
373    /**
374     * Returns the reduced chi-square of track x-projection
375     */
376    Float_t  TrkTrack::GetChi2X(){
377        float chiq=0;
378        for(int ip=0; ip<6; ip++)if(XGood(ip))chiq+= pow((xv[ip]-xm[ip])/resx[ip],2.);
379        if(GetNX()>3)chiq=chiq/(GetNX()-3);
380        else chiq=0;
381        if(chiq==0)cout << " Float_t  TrkTrack::GetChi2X() -- WARNING -- value not defined "<<chiq<<endl;
382        return chiq;
383    }
384    /**
385     * Returns the reduced chi-square of track y-projection
386     */
387    Float_t  TrkTrack::GetChi2Y(){
388        float chiq=0;
389        for(int ip=0; ip<6; ip++)if(YGood(ip))chiq+= pow((yv[ip]-ym[ip])/resy[ip],2.);
390        if(GetNY()>2)chiq=chiq/(GetNY()-2);
391        else chiq=0;
392        if(chiq==0)cout << " Float_t  TrkTrack::GetChi2Y() -- WARNING -- value not defined "<<chiq<<endl;
393        return chiq;
394    }
395    /**
396     * Returns the logarythm of the likeliwood-function of  track x-projection
397     */
398    Float_t TrkTrack::GetLnLX(){
399        float lnl=0;
400        for(int ip=0; ip<6; ip++)
401            if( XGood(ip) && tailx[ip]!=0 )
402                lnl += (tailx[ip]+1.) * log( (tailx[ip]*pow(resx[ip],2.) + pow(xv[ip]-xm[ip],2.)) / (tailx[ip]*pow(resx[ip],2)) );
403        if(GetNX()>3)lnl=lnl/(GetNX()-3);
404        else lnl=0;
405        if(lnl==0){
406            cout << " Float_t  TrkTrack::GetLnLX() -- WARNING -- value not defined "<<lnl<<endl;
407            Dump();
408        }
409        return lnl;
410        
411    }
412    /**
413     * Returns the logarythm of the likeliwood-function of  track y-projection
414     */
415    Float_t TrkTrack::GetLnLY(){
416        float lnl=0;
417        for(int ip=0; ip<6; ip++)
418            if( YGood(ip) && taily[ip]!=0 )
419                lnl += (taily[ip]+1.) * log( (taily[ip]*pow(resy[ip],2.) + pow(yv[ip]-ym[ip],2.)) / (taily[ip]*pow(resy[ip],2)) );
420        if(GetNY()>2)lnl=lnl/(GetNY()-2);
421        else lnl=0;
422        if(lnl==0){
423            cout << " Float_t  TrkTrack::GetLnLY() -- WARNING -- value not defined "<<lnl<<endl;
424            Dump();
425        }
426        return lnl;
427        
428    }
429    /**
430     * Returns the effective angle, relative to the sensor, on each plane.
431     * @param ip plane (0-5)
432     * @param iv view (0=x 1=y)
433     */
434    Float_t TrkTrack::GetEffectiveAngle(int ip, int iv){
435    
436        if(ip<0 || ip>5){
437            cout << "Float_t TrkTrack::GetEffectiveAngle(int "<<ip<<", int "<<iv<<") ==> wrong input"<<endl;
438            return 0.;
439        }
440    
441        float v[3]={xv[ip],yv[ip],zv[ip]};
442        //-----------------------------------------
443        // effective angle (relative to the sensor)
444        //-----------------------------------------
445        float axv_geo  = axv[ip];
446        float muhall_h = 297.61; //cm**2/Vs
447        float BY = TrkParams::GetBY(v);
448        float axv_eff = 0;
449        if(ip==5) axv_geo = -1*axv_geo;
450        if(ip==5) BY      = -1*BY;
451        axv_eff = 180.*atan( tan(axv_geo*acos(-1.)/180.) + muhall_h * BY * 0.0001)/acos(-1.);
452        //-----------------------------------------
453        // effective angle (relative to the sensor)
454        //-----------------------------------------
455        float ayv_geo = ayv[ip];
456        float muhall_e = 1258.18; //cm**2/Vs
457        float BX = TrkParams::GetBX(v);
458        float ayv_eff = 0;
459        ayv_eff = 180.*atan( tan(ayv_geo*acos(-1.)/180.) + muhall_e * BX * 0.0001)/acos(-1.);
460      
461        if     (iv==0)return axv_eff;
462        else if(iv==1)return ayv_eff;
463        else{
464            cout << "Float_t TrkTrack::GetEffectiveAngle(int "<<ip<<", int "<<iv<<") ==> wrong input"<<endl;
465            return 0.;
466        }
467      
468    };
469    
470  //--------------------------------------  //--------------------------------------
471  //  //
472  //  //
# Line 395  void TrkTrack::Dump(){ Line 497  void TrkTrack::Dump(){
497      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]<<" ";
498      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] << " ";
499      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] << " ";
500        cout << endl << "maxs x   : "; for(int i=0; i<6; i++)cout << GetClusterX_MaxStrip(i) << " ";
501        cout << endl << "maxs y   : "; for(int i=0; i<6; i++)cout << GetClusterY_MaxStrip(i) << " ";
502        cout << endl << "mult x   : "; for(int i=0; i<6; i++)cout << GetClusterX_Multiplicity(i) << " ";
503        cout << endl << "mult y   : "; for(int i=0; i<6; i++)cout << GetClusterY_Multiplicity(i) << " ";
504        cout << endl << "seed x   : "; for(int i=0; i<6; i++)cout << GetClusterX_Seed(i) << " ";
505        cout << endl << "seed y   : "; for(int i=0; i<6; i++)cout << GetClusterY_Seed(i) << " ";
506        cout << endl << "xpu      : "; for(int i=0; i<6; i++)cout << xpu[i] << " ";
507        cout << endl << "ypu      : "; for(int i=0; i<6; i++)cout << ypu[i] << " ";
508    
509      cout << endl;      cout << endl;
510  }  }
511  /**  /**
# Line 611  void TrkTrack::SetFromMiniStruct(cMini2t Line 722  void TrkTrack::SetFromMiniStruct(cMini2t
722   *   *
723   * @see TrkParams::SetPFA(int)   * @see TrkParams::SetPFA(int)
724   */   */
725  void TrkTrack::EvaluateClusterPositions(){  Bool_t TrkTrack::EvaluateClusterPositions(){
726            
727  //     cout << "void TrkTrack::GetClusterPositions() "<<endl;  //     cout << "void TrkTrack::GetClusterPositions() "<<endl;
728    
729      TrkParams::Load( );      TrkParams::Load( );
730      if( !TrkParams::IsLoaded() )return;      if( !TrkParams::IsLoaded() )return false;
731            
732      for(int ip=0; ip<6; ip++){      for(int ip=0; ip<6; ip++){
733  //      cout << ip<<" ** "<<xm[ip]<<" / "<<ym[ip]<<endl;;  //      cout << ip<<" ** "<<xm[ip]<<" / "<<ym[ip]<<endl;;
# Line 635  void TrkTrack::EvaluateClusterPositions( Line 746  void TrkTrack::EvaluateClusterPositions(
746          float bfy = 10*TrkParams::GetBY(v);//Tesla          float bfy = 10*TrkParams::GetBY(v);//Tesla
747          int ipp=ip+1;          int ipp=ip+1;
748          xyzpam_(&ipp,&icx,&icy,&ladder,&sensor,&ax,&ay,&bfx,&bfy);          xyzpam_(&ipp,&icx,&icy,&ladder,&sensor,&ax,&ay,&bfx,&bfy);
749            if(icx<0 || icy<0)return false;
750      }      }
751        return true;
752  }  }
753  /**  /**
754   * \brief Tracking method. It calls F77 mini routine.   * \brief Tracking method. It calls F77 mini routine.
# Line 673  void TrkTrack::Fit(double pfixed, int& f Line 786  void TrkTrack::Fit(double pfixed, int& f
786    
787      extern cMini2track track_;      extern cMini2track track_;
788      fail = 0;      fail = 0;
789    
790      FillMiniStruct(track_);      FillMiniStruct(track_);
791                
792      if(froml1!=0)EvaluateClusterPositions();      if(froml1!=0){
793            if( !EvaluateClusterPositions() ){
794                cout << "void TrkTrack::Fit("<<pfixed<<","<<fail<<","<<iprint<<","<<froml1<<") --- ERROR evaluating cluster positions "<<endl;
795                FillMiniStruct(track_) ;
796                fail = 1;
797                return;
798            }
799        }else{
800            FillMiniStruct(track_);
801        }
802            
803      // if fit variables have been reset, evaluate the initial guess      // if fit variables have been reset, evaluate the initial guess
804      if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.)guess_();      if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.)guess_();
# Line 786  Int_t TrkTrack::GetClusterX_ID(int ip){ Line 909  Int_t TrkTrack::GetClusterX_ID(int ip){
909  Int_t TrkTrack::GetClusterY_ID(int ip){  Int_t TrkTrack::GetClusterY_ID(int ip){
910      return ((Int_t)fabs(ygood[ip]))%10000000-1;      return ((Int_t)fabs(ygood[ip]))%10000000-1;
911  };  };
912    
913  /**  /**
914   * Method to retrieve the ladder (0-4, increasing x) traversed by the track on this plane.   * Method to retrieve the ladder (0-4, increasing x) traversed by the track on this plane.
915   * If no ladder is traversed (dead area) the metod retuns -1.   * If no ladder is traversed (dead area) the metod retuns -1.
# Line 969  void TrkLevel2::Dump(){ Line 1093  void TrkLevel2::Dump(){
1093                    
1094          //          //
1095      cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-";      cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-";
1096      cout << endl << "good     : "; for(int i=0; i<12; i++) cout << good[i]<<" ";      cout << endl << "good     : "; for(int i=0; i<12; i++) cout << hex <<" 0x"<< good[i]<<dec;
1097      cout << endl << "ntrk()   : " << this->ntrk() ;      cout << endl << "ntrk()   : " << ntrk() ;
1098      cout << endl << "nclsx()  : " << this->nclsx();      cout << endl << "nclsx()  : " << nclsx();
1099      cout << endl << "nclsy()  : " << this->nclsy();      cout << endl << "nclsy()  : " << nclsy();
1100      if(Track){      if(Track){
1101          TClonesArray &t  = *Track;          TClonesArray &t  = *Track;
1102          for(int i=0; i<ntrk(); i++)     ((TrkTrack *)t[i])->Dump();          for(int i=0; i<ntrk(); i++)     ((TrkTrack *)t[i])->Dump();
1103      }            }      
1104      if(SingletX){  //     if(SingletX){
1105          TClonesArray &sx = *SingletX;  //      TClonesArray &sx = *SingletX;
1106          for(int i=0; i<nclsx(); i++) ((TrkSinglet *)sx[i])->Dump();  //      for(int i=0; i<nclsx(); i++) ((TrkSinglet *)sx[i])->Dump();
1107      }  //     }
1108      if(SingletY){  //     if(SingletY){
1109          TClonesArray &sy = *SingletY;  //      TClonesArray &sy = *SingletY;
1110          for(int i=0; i<nclsy(); i++) ((TrkSinglet *)sy[i])->Dump();  //      for(int i=0; i<nclsy(); i++) ((TrkSinglet *)sy[i])->Dump();
1111      }  //     }
1112        cout << endl;
1113  }  }
1114    /**
1115     * \brief Dump processing status
1116     */
1117    void TrkLevel2::StatusDump(int view){
1118        cout << "DSP n. "<<view+1<<" status: "<<hex<<good[view]<<endl;    
1119    };
1120    /**
1121     * \brief Check event status
1122     *
1123     * Check the event status, according to a flag-mask given as input.
1124     * Return true if the view passes the check.
1125     *
1126     * @param view View number (0-11)
1127     * @param flagmask Mask of flags to check (eg. flagmask=0x111 no missing packet,
1128     *  no crc error, no software alarm)
1129     *
1130     * @see TrkLevel2 class definition to know how the status flag is defined
1131     *
1132     */
1133    Bool_t TrkLevel2::StatusCheck(int view, int flagmask){
1134    
1135        if( view<0 || view >= 12)return false;
1136        return !(good[view]&flagmask);
1137    
1138    };
1139    
1140    
1141  //--------------------------------------  //--------------------------------------
1142  //  //
1143  //  //
# Line 1101  void TrkLevel2::SetFromLevel2Struct(cTrk Line 1253  void TrkLevel2::SetFromLevel2Struct(cTrk
1253              t_track->ayv[ip]    = l2->ayv_nt[i][ip];              t_track->ayv[ip]    = l2->ayv_nt[i][ip];
1254              t_track->dedx_x[ip] = l2->dedx_x[i][ip];              t_track->dedx_x[ip] = l2->dedx_x[i][ip];
1255              t_track->dedx_y[ip] = l2->dedx_y[i][ip];              t_track->dedx_y[ip] = l2->dedx_y[i][ip];
1256                t_track->multmaxx[ip] = l2->multmaxx[i][ip];
1257                t_track->multmaxy[ip] = l2->multmaxy[i][ip];
1258                t_track->seedx[ip]  = l2->seedx[i][ip];  
1259                t_track->seedy[ip]  = l2->seedy[i][ip];
1260                t_track->xpu[ip]    = l2->xpu[i][ip];  
1261                t_track->ypu[ip]    = l2->ypu[i][ip];  
1262              //-----------------------------------------------------              //-----------------------------------------------------
1263              //-----------------------------------------------------              //-----------------------------------------------------
1264              //-----------------------------------------------------              //-----------------------------------------------------
Legend:
Removed from v.1.35  
changed lines
  Added in v.1.40
  ViewVC Help
Powered by ViewVC 1.1.23