/[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.33 by pam-fi, Mon May 14 11:03:05 2007 UTC revision 1.43 by pam-fi, Tue Jan 22 08:55:07 2008 UTC
# Line 12  using namespace std; Line 12  using namespace std;
12  extern "C" {      extern "C" {    
13      void dotrack_(int*, double*, double*, double*, double*, int*);      void dotrack_(int*, double*, double*, double*, double*, int*);
14      void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);      void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);
15       void mini2_(int*,int*,int*);      void mini2_(int*,int*,int*);
16       void guess_();      void guess_();
17        void gufld_(float*, float*);
18        float risxeta2_(float *);
19        float risxeta3_(float *);
20        float risxeta4_(float *);
21        float risyeta2_(float *);
22  }  }
23    
24  //--------------------------------------  //--------------------------------------
# Line 47  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      };      };
63  //    clx = 0;  
64  //    cly = 0;  //     TrkParams::SetTrackingMode();
65  //    clx = new TRefArray(6,0); //forse causa memory leak???  //     TrkParams::SetPrecisionFactor();
66  //    cly = new TRefArray(6,0); //forse causa memory leak???  //     TrkParams::SetStepMin();
67  //    clx = TRefArray(6,0);      TrkParams::SetMiniDefault();
 //    cly = TRefArray(6,0);  
   
     TrkParams::SetTrackingMode();  
     TrkParams::SetPrecisionFactor();  
     TrkParams::SetStepMin();  
68      TrkParams::SetPFA();      TrkParams::SetPFA();
69    
70  };  };
# Line 91  TrkTrack::TrkTrack(const TrkTrack& t){ Line 98  TrkTrack::TrkTrack(const TrkTrack& t){
98          ayv[ip]    = t.ayv[ip];          ayv[ip]    = t.ayv[ip];
99          dedx_x[ip] = t.dedx_x[ip];          dedx_x[ip] = t.dedx_x[ip];
100          dedx_y[ip] = t.dedx_y[ip];          dedx_y[ip] = t.dedx_y[ip];
101            multmaxx[ip] = t.multmaxx[ip];
102            multmaxy[ip] = t.multmaxy[ip];
103            seedx[ip]    = t.seedx[ip];  
104            seedy[ip]    = t.seedy[ip];
105            xpu[ip]      = t.xpu[ip];  
106            ypu[ip]      = t.ypu[ip];  
107      };      };
108  //    clx = 0;  
109  //    cly = 0;  //     TrkParams::SetTrackingMode();
110  //    if(t.clx)clx = new TRefArray(*(t.clx));  //     TrkParams::SetPrecisionFactor();
111  //    if(t.cly)cly = new TRefArray(*(t.cly));  //     TrkParams::SetStepMin();  
112  //    clx = TRefArray(t.clx);      TrkParams::SetMiniDefault();
 //    cly = TRefArray(t.cly);  
   
     TrkParams::SetTrackingMode();  
     TrkParams::SetPrecisionFactor();  
     TrkParams::SetStepMin();    
113      TrkParams::SetPFA();      TrkParams::SetPFA();
114    
115  };  };
# Line 136  void TrkTrack::Copy(TrkTrack& t){ Line 144  void TrkTrack::Copy(TrkTrack& t){
144          t.ayv[ip]    = ayv[ip];          t.ayv[ip]    = ayv[ip];
145          t.dedx_x[ip] = dedx_x[ip];          t.dedx_x[ip] = dedx_x[ip];
146          t.dedx_y[ip] = dedx_y[ip];          t.dedx_y[ip] = dedx_y[ip];
147            t.multmaxx[ip] = multmaxx[ip];
148            t.multmaxy[ip] = multmaxy[ip];
149            t.seedx[ip]    = seedx[ip];  
150            t.seedy[ip]    = seedy[ip];
151            t.xpu[ip]      = xpu[ip];  
152            t.ypu[ip]      = ypu[ip];  
153                            
154      };      };
155    
 //    t.clx = TRefArray(clx);  
 //    t.cly = TRefArray(cly);  
156            
157  };  };
158  //--------------------------------------  //--------------------------------------
# Line 360  Int_t TrkTrack::GetLeverArmY(){ Line 372  Int_t TrkTrack::GetLeverArmY(){
372      }      }
373      return (last_plane-first_plane+1);      return (last_plane-first_plane+1);
374  }  }
375    /**
376     * Returns the reduced chi-square of track x-projection
377     */
378    Float_t  TrkTrack::GetChi2X(){
379        float chiq=0;
380        for(int ip=0; ip<6; ip++)if(XGood(ip))chiq+= pow((xv[ip]-xm[ip])/resx[ip],2.);
381        if(GetNX()>3)chiq=chiq/(GetNX()-3);
382        else chiq=0;
383        if(chiq==0)cout << " Float_t  TrkTrack::GetChi2X() -- WARNING -- value not defined "<<chiq<<endl;
384        return chiq;
385    }
386    /**
387     * Returns the reduced chi-square of track y-projection
388     */
389    Float_t  TrkTrack::GetChi2Y(){
390        float chiq=0;
391        for(int ip=0; ip<6; ip++)if(YGood(ip))chiq+= pow((yv[ip]-ym[ip])/resy[ip],2.);
392        if(GetNY()>2)chiq=chiq/(GetNY()-2);
393        else chiq=0;
394        if(chiq==0)cout << " Float_t  TrkTrack::GetChi2Y() -- WARNING -- value not defined "<<chiq<<endl;
395        return chiq;
396    }
397    /**
398     * Returns the logarythm of the likeliwood-function of  track x-projection
399     */
400    Float_t TrkTrack::GetLnLX(){
401        float lnl=0;
402        for(int ip=0; ip<6; ip++)
403            if( XGood(ip) && tailx[ip]!=0 )
404                lnl += (tailx[ip]+1.) * log( (tailx[ip]*pow(resx[ip],2.) + pow(xv[ip]-xm[ip],2.)) / (tailx[ip]*pow(resx[ip],2)) );
405        if(GetNX()>3)lnl=lnl/(GetNX()-3);
406        else lnl=0;
407        if(lnl==0){
408            cout << " Float_t  TrkTrack::GetLnLX() -- WARNING -- value not defined "<<lnl<<endl;
409            Dump();
410        }
411        return lnl;
412        
413    }
414    /**
415     * Returns the logarythm of the likeliwood-function of  track y-projection
416     */
417    Float_t TrkTrack::GetLnLY(){
418        float lnl=0;
419        for(int ip=0; ip<6; ip++)
420            if( YGood(ip) && taily[ip]!=0 )
421                lnl += (taily[ip]+1.) * log( (taily[ip]*pow(resy[ip],2.) + pow(yv[ip]-ym[ip],2.)) / (taily[ip]*pow(resy[ip],2)) );
422        if(GetNY()>2)lnl=lnl/(GetNY()-2);
423        else lnl=0;
424        if(lnl==0){
425            cout << " Float_t  TrkTrack::GetLnLY() -- WARNING -- value not defined "<<lnl<<endl;
426            Dump();
427        }
428        return lnl;
429        
430    }
431    /**
432     * Returns the effective angle, relative to the sensor, on each plane.
433     * @param ip plane (0-5)
434     * @param iv view (0=x 1=y)
435     */
436    Float_t TrkTrack::GetEffectiveAngle(int ip, int iv){
437    
438        if(ip<0 || ip>5){
439            cout << "Float_t TrkTrack::GetEffectiveAngle(int "<<ip<<", int "<<iv<<") ==> wrong input"<<endl;
440            return 0.;
441        }
442    
443        float v[3]={xv[ip],yv[ip],zv[ip]};
444        //-----------------------------------------
445        // effective angle (relative to the sensor)
446        //-----------------------------------------
447        float axv_geo  = axv[ip];
448        float muhall_h = 297.61; //cm**2/Vs
449        float BY = TrkParams::GetBY(v);
450        float axv_eff = 0;
451        if(ip==5) axv_geo = -1*axv_geo;
452        if(ip==5) BY      = -1*BY;
453        axv_eff = 180.*atan( tan(axv_geo*acos(-1.)/180.) + muhall_h * BY * 0.0001)/acos(-1.);
454        //-----------------------------------------
455        // effective angle (relative to the sensor)
456        //-----------------------------------------
457        float ayv_geo = ayv[ip];
458        float muhall_e = 1258.18; //cm**2/Vs
459        float BX = TrkParams::GetBX(v);
460        float ayv_eff = 0;
461        ayv_eff = 180.*atan( tan(ayv_geo*acos(-1.)/180.) + muhall_e * BX * 0.0001)/acos(-1.);
462      
463        if     (iv==0)return axv_eff;
464        else if(iv==1)return ayv_eff;
465        else{
466            cout << "Float_t TrkTrack::GetEffectiveAngle(int "<<ip<<", int "<<iv<<") ==> wrong input"<<endl;
467            return 0.;
468        }
469      
470    };
471    
472  //--------------------------------------  //--------------------------------------
473  //  //
474  //  //
# Line 381  void TrkTrack::Dump(){ Line 490  void TrkTrack::Dump(){
490      cout << endl << "zv       : "; for(int i=0; i<6; i++)cout << zv[i] << " ";      cout << endl << "zv       : "; for(int i=0; i<6; i++)cout << zv[i] << " ";
491      cout << endl << "resx     : "; for(int i=0; i<6; i++)cout << resx[i] << " ";      cout << endl << "resx     : "; for(int i=0; i<6; i++)cout << resx[i] << " ";
492      cout << endl << "resy     : "; for(int i=0; i<6; i++)cout << resy[i] << " ";      cout << endl << "resy     : "; for(int i=0; i<6; i++)cout << resy[i] << " ";
493        cout << endl << "tailx    : "; for(int i=0; i<6; i++)cout << tailx[i] << " ";
494        cout << endl << "taily    : "; for(int i=0; i<6; i++)cout << taily[i] << " ";
495      cout << endl << "coval    : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" ";      cout << endl << "coval    : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" ";
496      cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[1][i]<<" ";      cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[1][i]<<" ";
497      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]<<" ";
# Line 388  void TrkTrack::Dump(){ Line 499  void TrkTrack::Dump(){
499      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]<<" ";
500      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] << " ";
501      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] << " ";
502        cout << endl << "maxs x   : "; for(int i=0; i<6; i++)cout << GetClusterX_MaxStrip(i) << " ";
503        cout << endl << "maxs y   : "; for(int i=0; i<6; i++)cout << GetClusterY_MaxStrip(i) << " ";
504        cout << endl << "mult x   : "; for(int i=0; i<6; i++)cout << GetClusterX_Multiplicity(i) << " ";
505        cout << endl << "mult y   : "; for(int i=0; i<6; i++)cout << GetClusterY_Multiplicity(i) << " ";
506        cout << endl << "seed x   : "; for(int i=0; i<6; i++)cout << GetClusterX_Seed(i) << " ";
507        cout << endl << "seed y   : "; for(int i=0; i<6; i++)cout << GetClusterY_Seed(i) << " ";
508        cout << endl << "xpu      : "; for(int i=0; i<6; i++)cout << xpu[i] << " ";
509        cout << endl << "ypu      : "; for(int i=0; i<6; i++)cout << ypu[i] << " ";
510    
511      cout << endl;      cout << endl;
512  }  }
513  /**  /**
# Line 406  void TrkTrack::SetResolution(double *rx, Line 526  void TrkTrack::SetResolution(double *rx,
526      for(int i=0; i<6; i++) resy[i]=*ry++;      for(int i=0; i<6; i++) resy[i]=*ry++;
527  }  }
528  /**  /**
529   * Set the TrkTrack good measurement.   * Set the TrkTrack tails position resolution
530     */
531    void TrkTrack::SetTail(double *tx, double *ty, double factor){
532        for(int i=0; i<6; i++) tailx[i]=factor*(*tx++);
533        for(int i=0; i<6; i++) taily[i]=factor*(*ty++);
534    }
535    /**
536     * Set the TrkTrack Student parameter (resx,resy,tailx,taily)
537     * from previous gausian fit
538     *@param flag =0 standard, =1 with noise correction
539     */
540    void TrkTrack::SetStudentParam(int flag){
541        float sx[11]={0.000128242,
542                       0.000136942,
543                       0.000162718,
544                       0.000202644,
545                       0.00025597,
546                       0.000317456,
547                       0.000349048,
548                       0.000384638,
549                       0.000457295,
550                       0.000512319,
551                       0.000538573};
552        float tx[11]={1.79402,
553                       2.04876,
554                       2.88376,
555                       3.3,
556                       3.14084,
557                       4.07686,
558                       4.44736,
559                       3.5179,
560                       3.38697,
561                       3.45739,
562                       3.18627};
563        float sy[11]={0.000483075,
564                       0.000466925,
565                       0.000431658,
566                       0.000428317,
567                       0.000433854,
568                       0.000444044,
569                       0.000482098,
570                       0.000537579,
571                       0.000636279,
572                       0.000741998,
573                       0.000864261};
574        float ty[11]={0.997032,
575                       1.11147,
576                       1.18526,
577                       1.61404,
578                       2.21908,
579                       3.08959,
580                       4.48833,
581                       4.42687,
582                       4.65253,
583                       4.52043,
584                       4.29926};
585        int index;
586        float fact;
587        for(int i=0; i<6; i++) {
588            index = int((fabs(axv[i])+1.)/2.);
589            if(index>10) index=10;
590            tailx[i]=tx[index];
591            if(flag==1) {
592                if(fabs(axv[i])<=10.) fact = resx[i]/risxeta2_(&(axv[i]));
593                if(fabs(axv[i])>10.&&fabs(axv[i])<=15.) fact = resx[i]/risxeta3_(&(axv[i]));
594                if(fabs(axv[i])>15.) fact = resx[i]/risxeta4_(&(axv[i]));
595            } else fact = 1.;
596            resx[i] = sx[index]*fact;
597        }
598        for(int i=0; i<6; i++) {
599            index = int((fabs(ayv[i])+1.)/2.);
600            if(index>10) index=10;
601            taily[i]=ty[index];
602            if(flag==1) fact = resy[i]/risyeta2_(&(ayv[i]));
603            else fact = 1.;
604            resy[i] = sy[index]*fact;
605        }
606    }
607    /**
608     * Set the TrkTrack good measurement
609   */   */
610  void TrkTrack::SetGood(int *xg, int *yg){  void TrkTrack::SetGood(int *xg, int *yg){
611    
# Line 424  void TrkTrack::LoadField(TString path){ Line 623  void TrkTrack::LoadField(TString path){
623  //     path_.error   = 0;  //     path_.error   = 0;
624  //     readb_();  //     readb_();
625    
626      TrkParams::SetTrackingMode();  //     TrkParams::SetTrackingMode();
627      TrkParams::SetPrecisionFactor();  //     TrkParams::SetPrecisionFactor();
628      TrkParams::SetStepMin();  //     TrkParams::SetStepMin();
629        TrkParams::SetMiniDefault();
630    
631      TrkParams::Set(path,1);      TrkParams::Set(path,1);
632      TrkParams::Load(1);      TrkParams::Load(1);
# Line 468  void TrkTrack::FillMiniStruct(cMini2trac Line 668  void TrkTrack::FillMiniStruct(cMini2trac
668                    
669          track.resx[i]=resx[i];          track.resx[i]=resx[i];
670          track.resy[i]=resy[i];          track.resy[i]=resy[i];
671            track.tailx[i]=tailx[i];
672            track.taily[i]=taily[i];
673      }      }
674    
675      for(int i=0; i<5; i++) track.al[i]=al[i];      for(int i=0; i<5; i++) track.al[i]=al[i];
# Line 523  void TrkTrack::SetFromMiniStruct(cMini2t Line 725  void TrkTrack::SetFromMiniStruct(cMini2t
725   *   *
726   * @see TrkParams::SetPFA(int)   * @see TrkParams::SetPFA(int)
727   */   */
728  void TrkTrack::EvaluateClusterPositions(){  Bool_t TrkTrack::EvaluateClusterPositions(){
729            
730  //     cout << "void TrkTrack::GetClusterPositions() "<<endl;  //     cout << "void TrkTrack::GetClusterositions() "<<endl;
731    
732      TrkParams::Load( );      TrkParams::Load( );
733      if( !TrkParams::IsLoaded() )return;      if( !TrkParams::IsLoaded() )return false;
734            
735      for(int ip=0; ip<6; ip++){      for(int ip=0; ip<6; ip++){
736  //      cout << ip<<" ** "<<xm[ip]<<" / "<<ym[ip]<<endl;;  //      cout << ip<<" ** "<<xm[ip]<<" / "<<ym[ip]<<endl;;
# Line 547  void TrkTrack::EvaluateClusterPositions( Line 749  void TrkTrack::EvaluateClusterPositions(
749          float bfy = 10*TrkParams::GetBY(v);//Tesla          float bfy = 10*TrkParams::GetBY(v);//Tesla
750          int ipp=ip+1;          int ipp=ip+1;
751          xyzpam_(&ipp,&icx,&icy,&ladder,&sensor,&ax,&ay,&bfx,&bfy);          xyzpam_(&ipp,&icx,&icy,&ladder,&sensor,&ax,&ay,&bfx,&bfy);
752            if(icx<0 || icy<0)return false;
753      }      }
754        return true;
755  }  }
756  /**  /**
757   * \brief Tracking method. It calls F77 mini routine.   * \brief Tracking method. It calls F77 mini routine.
# Line 571  void TrkTrack::EvaluateClusterPositions( Line 775  void TrkTrack::EvaluateClusterPositions(
775   * @see EvaluateClusterPositions()   * @see EvaluateClusterPositions()
776   *   *
777   * The fitting procedure can be varied by changing the tracking mode,   * The fitting procedure can be varied by changing the tracking mode,
778   * the fit-precision factor and the minimum number of step.   * the fit-precision factor, the minimum number of step, etc.
779   * @see SetTrackingMode(int)   * @see SetTrackingMode(int)
780   * @see SetPrecisionFactor(double)   * @see SetPrecisionFactor(double)
781   * @see SetStepMin(int)   * @see SetStepMin(int)
782     * @see SetDeltaB(int,double)
783   */   */
784  void TrkTrack::Fit(double pfixed, int& fail, int iprint, int froml1){  void TrkTrack::Fit(double pfixed, int& fail, int iprint, int froml1){
785    
# Line 585  void TrkTrack::Fit(double pfixed, int& f Line 790  void TrkTrack::Fit(double pfixed, int& f
790    
791      extern cMini2track track_;      extern cMini2track track_;
792      fail = 0;      fail = 0;
793    
794      FillMiniStruct(track_);      FillMiniStruct(track_);
795                
796      if(froml1!=0)EvaluateClusterPositions();      if(froml1!=0){
797            if( !EvaluateClusterPositions() ){
798                cout << "void TrkTrack::Fit("<<pfixed<<","<<fail<<","<<iprint<<","<<froml1<<") --- ERROR evaluating cluster positions "<<endl;
799                FillMiniStruct(track_) ;
800                fail = 1;
801                return;
802            }
803        }else{
804            FillMiniStruct(track_);
805        }
806            
807      // if fit variables have been reset, evaluate the initial guess      // if fit variables have been reset, evaluate the initial guess
808      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 660  void TrkTrack::SetPrecisionFactor(double Line 875  void TrkTrack::SetPrecisionFactor(double
875      track_.fact = fact;      track_.fact = fact;
876  }  }
877  /**  /**
878   * Set the factor scale for tracking precision   * Set the minimum number of steps for tracking precision
879   */   */
880  void TrkTrack::SetStepMin(int istepmin){  void TrkTrack::SetStepMin(int istepmin){
881      extern cMini2track track_;      extern cMini2track track_;
882      track_.istepmin = istepmin;      track_.istepmin = istepmin;
883  }  }
884    /**
885     * Set deltaB parameters (id=0,1). By default they are set to zero.
886     */
887    void TrkTrack::SetDeltaB(int id, double db){
888        if(id!=0 && id!=1)cout << "void TrkTrack::SetDeltaB(int id,double db) -- wrong input parameters: "<<id<<" "<<db<<endl;
889        TrkParams::SetDeltaB(id,db);
890    }
891    
892  /**  /**
893     * Returns 1 if the track is inside the magnet cavity
894     * Set the minimum number of steps for tracking precision
895     */
896    Bool_t TrkTrack::IsInsideCavity(){
897        float xmagntop, ymagntop, xmagnbottom, ymagnbottom;
898        xmagntop = xv[0] + (ZMAGNHIGH-zv[0])*tan(acos(-1.0)*axv[0]/180.);
899        ymagntop = yv[0] + (ZMAGNHIGH-zv[0])*tan(acos(-1.0)*ayv[0]/180.);
900        xmagnbottom = xv[5] + (ZMAGNLOW-zv[5])*tan(acos(-1.0)*axv[5]/180.);
901        ymagnbottom = yv[5] + (ZMAGNLOW-zv[5])*tan(acos(-1.0)*ayv[5]/180.);
902        if( xmagntop>XMAGNLOW && xmagntop<XMAGNHIGH &&
903            ymagntop>YMAGNLOW && ymagntop<YMAGNHIGH &&
904            xmagnbottom>XMAGNLOW && xmagnbottom<XMAGNHIGH &&
905            ymagnbottom>YMAGNLOW && ymagnbottom<YMAGNHIGH ) return(true);
906        else return(false);
907    }
908    /**
909   * Method to retrieve ID (0,1,...) of x-cluster (if any) associated to this track.   * Method to retrieve ID (0,1,...) of x-cluster (if any) associated to this track.
910   * If no cluster is associated, ID=-1.   * If no cluster is associated, ID=-1.
911   * @param ip Tracker plane (0-5)   * @param ip Tracker plane (0-5)
# Line 684  Int_t TrkTrack::GetClusterX_ID(int ip){ Line 921  Int_t TrkTrack::GetClusterX_ID(int ip){
921  Int_t TrkTrack::GetClusterY_ID(int ip){  Int_t TrkTrack::GetClusterY_ID(int ip){
922      return ((Int_t)fabs(ygood[ip]))%10000000-1;      return ((Int_t)fabs(ygood[ip]))%10000000-1;
923  };  };
924    
925  /**  /**
926   * 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.
927   * If no ladder is traversed (dead area) the metod retuns -1.   * If no ladder is traversed (dead area) the metod retuns -1.
# Line 732  void TrkTrack::SetYGood(int ip, int clid Line 970  void TrkTrack::SetYGood(int ip, int clid
970      if(bad)ygood[ip]=-ygood[ip];      if(bad)ygood[ip]=-ygood[ip];
971  };  };
972    
973    /**
974     * \brief Average X
975     * Average value of <xv>, evaluated from the first to the last hit x view.
976     */
977    Float_t TrkTrack::GetXav(){
978    
979        int first_plane = -1;
980        int last_plane  = -1;
981        for(Int_t ip=0; ip<6; ip++){
982            if( XGood(ip) && first_plane == -1 )first_plane = ip;
983            if( XGood(ip) && first_plane != -1 )last_plane = ip;
984        }
985        if( first_plane == -1 || last_plane == -1){
986            return -100;
987        }
988        if( last_plane-first_plane+1 ==0 )return -100;
989    
990        Float_t av = 0;    
991        for(int ip=first_plane; ip<=last_plane; ip++)av+=xv[ip];
992        
993        return (av/(last_plane-first_plane+1));    
994    }
995    /**
996     * \brief Average Y
997     * Average value of <yv>, evaluated from the first to the last hit x view.
998     */
999    Float_t TrkTrack::GetYav(){
1000    
1001        int first_plane = -1;
1002        int last_plane  = -1;
1003        for(Int_t ip=0; ip<6; ip++){
1004            if( XGood(ip) && first_plane == -1 )first_plane = ip;
1005            if( XGood(ip) && first_plane != -1 )last_plane = ip;
1006        }
1007        if( first_plane == -1 || last_plane == -1){
1008            return -100;
1009        }
1010        if( last_plane-first_plane+1 ==0 )return -100;
1011    
1012        Float_t av = 0;    
1013        for(int ip=first_plane; ip<=last_plane; ip++)av+=yv[ip];
1014    
1015        return (av/(last_plane-first_plane+1));    
1016    }
1017    /**
1018     * \brief Average Z
1019     * Average value of <zv>, evaluated from the first to the last hit x view.
1020     */
1021    Float_t TrkTrack::GetZav(){
1022    
1023        int first_plane = -1;
1024        int last_plane  = -1;
1025        for(Int_t ip=0; ip<6; ip++){
1026            if( XGood(ip) && first_plane == -1 )first_plane = ip;
1027            if( XGood(ip) && first_plane != -1 )last_plane = ip;
1028        }
1029        if( first_plane == -1 || last_plane == -1){
1030            return -100;
1031        }
1032        if( last_plane-first_plane+1 ==0 )return -100;
1033    
1034        Float_t av = 0;    
1035        for(int ip=first_plane; ip<=last_plane; ip++)av+=zv[ip];
1036        
1037        return (av/(last_plane-first_plane+1));    
1038    }
1039    
1040    /**
1041     * \brief Number of column traversed
1042     */
1043    Int_t TrkTrack::GetNColumns(){
1044        int sensors[] = {0,0,0,0,0,0};
1045        for(int ip=0; ip<6; ip++){
1046            int sensorid = GetLadder(ip)+3*GetSensor(ip);    
1047            if(XGood(ip)||YGood(ip))
1048                if(sensorid>=0 && sensorid<6)sensors[sensorid]=1;
1049        }
1050        int nsensors=0;
1051        for(int is=0; is<6; is++)nsensors += sensors[is];
1052        return nsensors;
1053    };
1054    /**
1055     * \brief Give the maximum energy release
1056     */
1057    Float_t TrkTrack::GetDEDX_max(int ip, int iv){
1058        Float_t max=0;
1059        int pfrom = 0;
1060        int pto   = 6;
1061        int vfrom = 0;
1062        int vto   = 2;
1063        if(ip>=0&&ip<6){
1064            pfrom = ip;
1065            pto   = ip+1;
1066        }
1067        if(iv>=0&&iv<2){
1068            vfrom = iv;
1069            vto   = iv+1;
1070        }
1071        for(int i=pfrom; i<pto; i++)
1072            for(int j=0; j<vto; j++)
1073                if(GetDEDX(i,j)>max)max=GetDEDX(i,j);
1074    
1075        return max;
1076    
1077    };
1078    
1079    /**
1080     * \brief Give the minimum energy release
1081     */
1082    Float_t TrkTrack::GetDEDX_min(int ip, int iv){
1083        Float_t min=100000000;
1084        int pfrom = 0;
1085        int pto   = 6;
1086        int vfrom = 0;
1087        int vto   = 2;
1088        if(ip>=0&&ip<6){
1089            pfrom = ip;
1090            pto   = ip+1;
1091        }
1092        if(iv>=0&&iv<2){
1093            vfrom = iv;
1094            vto   = iv+1;
1095        }
1096        for(int i=pfrom; i<pto; i++)
1097            for(int j=0; j<vto; j++)
1098                if(GetDEDX(i,j)<min)min=GetDEDX(i,j);
1099    
1100        return min;
1101    
1102    };
1103    
1104    /**
1105     * \brief Give the maximum spatial residual release
1106     */
1107    Float_t TrkTrack::GetResidual_max(int ip, int iv){
1108        Float_t max=0;
1109        int pfrom = 0;
1110        int pto   = 6;
1111        int vfrom = 0;
1112        int vto   = 2;
1113        if(ip>=0&&ip<6){
1114            pfrom = ip;
1115            pto   = ip+1;
1116        }
1117        if(iv>=0&&iv<2){
1118            vfrom = iv;
1119            vto   = iv+1;
1120        }
1121        for(int i=pfrom; i<pto; i++){
1122            for(int j=0; j<vto; j++){
1123                if(j==0 && XGood(i) && fabs(xm[i]-xv[i])>fabs(max))max=xv[i]-xm[i];
1124                if(j==1 && YGood(i) && fabs(ym[i]-yv[i])>fabs(max))max=yv[i]-ym[i];
1125            }
1126        }
1127        return max;
1128    
1129    };
1130    
1131    
1132    /**
1133     * \brief Give the maximum multiplicity on the x view
1134     */
1135    Int_t TrkTrack::GetClusterX_Multiplicity_max(){
1136        int max=0;
1137        for(int ip=0; ip<6; ip++)
1138            if(GetClusterX_Multiplicity(ip)>max)max=GetClusterX_Multiplicity(ip);
1139        return max;
1140    };
1141    /**
1142     * \brief Give the minimum multiplicity on the x view
1143     */
1144    Int_t TrkTrack::GetClusterX_Multiplicity_min(){
1145        int min=50;
1146        for(int ip=0; ip<6; ip++)
1147            if(GetClusterX_Multiplicity(ip)<min)min=GetClusterX_Multiplicity(ip);
1148        return min;
1149    };
1150    /**
1151     * \brief Give the maximum multiplicity on the x view
1152     */
1153    Int_t TrkTrack::GetClusterY_Multiplicity_max(){
1154        int max=0;
1155        for(int ip=0; ip<6; ip++)
1156            if(GetClusterY_Multiplicity(ip)>max)max=GetClusterY_Multiplicity(ip);
1157        return max;
1158    };
1159    /**
1160     * \brief Give the minimum multiplicity on the x view
1161     */
1162    Int_t TrkTrack::GetClusterY_Multiplicity_min(){
1163        int min=50;
1164        for(int ip=0; ip<6; ip++)
1165            if(GetClusterY_Multiplicity(ip)<min)min=GetClusterY_Multiplicity(ip);
1166        return min;
1167    };
1168    
1169    /**
1170     * \brief Give the minimum seed on the x view
1171     */
1172    Float_t TrkTrack::GetClusterX_Seed_min(){
1173        Float_t min=100000;
1174        for(int ip=0; ip<6; ip++)
1175            if(XGood(ip) && GetClusterX_Seed(ip)<min)min=GetClusterX_Seed(ip);
1176        return min;
1177    };
1178    /**
1179     * \brief Give the minimum seed on the x view
1180     */
1181    Float_t TrkTrack::GetClusterY_Seed_min(){
1182        Float_t min=100000;
1183        for(int ip=0; ip<6; ip++)
1184            if(YGood(ip) && GetClusterY_Seed(ip)<min)min=GetClusterY_Seed(ip);
1185        return min;
1186    };
1187    
1188    
1189  //--------------------------------------  //--------------------------------------
1190  //  //
1191  //  //
# Line 867  void TrkLevel2::Dump(){ Line 1321  void TrkLevel2::Dump(){
1321                    
1322          //          //
1323      cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-";      cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-";
1324      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;
1325      cout << endl << "ntrk()   : " << this->ntrk() ;      cout << endl << "ntrk()   : " << ntrk() ;
1326      cout << endl << "nclsx()  : " << this->nclsx();      cout << endl << "nclsx()  : " << nclsx();
1327      cout << endl << "nclsy()  : " << this->nclsy();      cout << endl << "nclsy()  : " << nclsy();
1328      if(Track){      if(Track){
1329          TClonesArray &t  = *Track;          TClonesArray &t  = *Track;
1330          for(int i=0; i<ntrk(); i++)     ((TrkTrack *)t[i])->Dump();          for(int i=0; i<ntrk(); i++)     ((TrkTrack *)t[i])->Dump();
1331      }            }      
1332      if(SingletX){  //     if(SingletX){
1333          TClonesArray &sx = *SingletX;  //      TClonesArray &sx = *SingletX;
1334          for(int i=0; i<nclsx(); i++) ((TrkSinglet *)sx[i])->Dump();  //      for(int i=0; i<nclsx(); i++) ((TrkSinglet *)sx[i])->Dump();
1335      }  //     }
1336      if(SingletY){  //     if(SingletY){
1337          TClonesArray &sy = *SingletY;  //      TClonesArray &sy = *SingletY;
1338          for(int i=0; i<nclsy(); i++) ((TrkSinglet *)sy[i])->Dump();  //      for(int i=0; i<nclsy(); i++) ((TrkSinglet *)sy[i])->Dump();
1339      }  //     }
1340        cout << endl;
1341  }  }
1342    /**
1343     * \brief Dump processing status
1344     */
1345    void TrkLevel2::StatusDump(int view){
1346        cout << "DSP n. "<<view+1<<" status: "<<hex<<good[view]<<endl;    
1347    };
1348    /**
1349     * \brief Check event status
1350     *
1351     * Check the event status, according to a flag-mask given as input.
1352     * Return true if the view passes the check.
1353     *
1354     * @param view View number (0-11)
1355     * @param flagmask Mask of flags to check (eg. flagmask=0x111 no missing packet,
1356     *  no crc error, no software alarm)
1357     *
1358     * @see TrkLevel2 class definition to know how the status flag is defined
1359     *
1360     */
1361    Bool_t TrkLevel2::StatusCheck(int view, int flagmask){
1362    
1363        if( view<0 || view >= 12)return false;
1364        return !(good[view]&flagmask);
1365    
1366    };
1367    
1368    
1369  //--------------------------------------  //--------------------------------------
1370  //  //
1371  //  //
# Line 999  void TrkLevel2::SetFromLevel2Struct(cTrk Line 1481  void TrkLevel2::SetFromLevel2Struct(cTrk
1481              t_track->ayv[ip]    = l2->ayv_nt[i][ip];              t_track->ayv[ip]    = l2->ayv_nt[i][ip];
1482              t_track->dedx_x[ip] = l2->dedx_x[i][ip];              t_track->dedx_x[ip] = l2->dedx_x[i][ip];
1483              t_track->dedx_y[ip] = l2->dedx_y[i][ip];              t_track->dedx_y[ip] = l2->dedx_y[i][ip];
1484                t_track->multmaxx[ip] = l2->multmaxx[i][ip];
1485                t_track->multmaxy[ip] = l2->multmaxy[i][ip];
1486                t_track->seedx[ip]  = l2->seedx[i][ip];  
1487                t_track->seedy[ip]  = l2->seedy[i][ip];
1488                t_track->xpu[ip]    = l2->xpu[i][ip];  
1489                t_track->ypu[ip]    = l2->ypu[i][ip];  
1490              //-----------------------------------------------------              //-----------------------------------------------------
1491              //-----------------------------------------------------              //-----------------------------------------------------
1492              //-----------------------------------------------------              //-----------------------------------------------------
# Line 1213  TRefArray *TrkLevel2::GetTracks_NFitSort Line 1701  TRefArray *TrkLevel2::GetTracks_NFitSort
1701  TrkTrack *TrkLevel2::GetStoredTrack(int is){  TrkTrack *TrkLevel2::GetStoredTrack(int is){
1702    
1703      if(is >= this->ntrk()){      if(is >= this->ntrk()){
1704          cout << "** TrkLevel2 ** Track "<< is << "doen not exits! " << endl;          cout << "TrkTrack *TrkLevel2::GetStoredTrack(int) >> Track "<< is << "doen not exits! " << endl;
1705          cout << "                Stored tracks ntrk() = "<< this->ntrk() << endl;          cout << "Stored tracks ntrk() = "<< this->ntrk() << endl;
1706          return 0;          return 0;
1707      }      }
1708      if(!Track){      if(!Track){
# Line 1235  TrkTrack *TrkLevel2::GetStoredTrack(int Line 1723  TrkTrack *TrkLevel2::GetStoredTrack(int
1723  TrkSinglet *TrkLevel2::GetSingletX(int is){  TrkSinglet *TrkLevel2::GetSingletX(int is){
1724    
1725          if(is >= this->nclsx()){          if(is >= this->nclsx()){
1726                  cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl;                  cout << "TrkSinglet *TrkLevel2::GetSingletX(int) >> Singlet "<< is << "doen not exits! " << endl;
1727                  cout << "                Stored x-singlets nclsx() = "<< this->nclsx() << endl;                  cout << "Stored x-singlets nclsx() = "<< this->nclsx() << endl;
1728                  return 0;                  return 0;
1729          }          }
1730          if(!SingletX)return 0;          if(!SingletX)return 0;
# Line 1255  TrkSinglet *TrkLevel2::GetSingletX(int i Line 1743  TrkSinglet *TrkLevel2::GetSingletX(int i
1743  TrkSinglet *TrkLevel2::GetSingletY(int is){  TrkSinglet *TrkLevel2::GetSingletY(int is){
1744    
1745          if(is >= this->nclsy()){          if(is >= this->nclsy()){
1746                  cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl;                  cout << "TrkSinglet *TrkLevel2::GetSingletY(int) >> Singlet "<< is << "doen not exits! " << endl;
1747                  cout << "                Stored y-singlets nclsy() = "<< this->nclsx() << endl;                  cout << "Stored y-singlets nclsx() = "<< this->nclsx() << endl;
1748                  return 0;                  return 0;
1749          }          }
1750          if(!SingletY)return 0;          if(!SingletY)return 0;
# Line 1276  TrkSinglet *TrkLevel2::GetSingletY(int i Line 1764  TrkSinglet *TrkLevel2::GetSingletY(int i
1764  TrkTrack *TrkLevel2::GetTrack(int it){  TrkTrack *TrkLevel2::GetTrack(int it){
1765            
1766          if(it >= this->GetNTracks()){          if(it >= this->GetNTracks()){
1767                  cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl;                  cout << "TrkTrack *TrkLevel2::GetTrack(int) >> Track "<< it << "does not exits! " << endl;
1768                  cout << "                Physical tracks GetNTracks() = "<< this->ntrk() << endl;                  cout << "Physical tracks GetNTracks() = "<< this->ntrk() << endl;
1769                  return 0;                  return 0;
1770          }          }
1771                    
# Line 1314  Int_t TrkLevel2::GetNTracks(){ Line 1802  Int_t TrkLevel2::GetNTracks(){
1802  TrkTrack *TrkLevel2::GetTrackImage(int it){  TrkTrack *TrkLevel2::GetTrackImage(int it){
1803    
1804      if(it >= this->GetNTracks()){      if(it >= this->GetNTracks()){
1805          cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl;          cout << "TrkTrack *TrkLevel2::GetTrackImage(int) >> Track "<< it << "does not exits! " << endl;
1806          cout << "                Physical tracks GetNTracks() = "<< this->ntrk() << endl;          cout << "Physical tracks GetNTracks() = "<< this->ntrk() << endl;
1807          return 0;          return 0;
1808      }      }
1809                    
# Line 1324  TrkTrack *TrkLevel2::GetTrackImage(int i Line 1812  TrkTrack *TrkLevel2::GetTrackImage(int i
1812      TrkTrack *track = (TrkTrack*)sorted->At(it);      TrkTrack *track = (TrkTrack*)sorted->At(it);
1813                    
1814      if(!track->HasImage()){      if(!track->HasImage()){
1815          cout << "** TrkLevel2 ** Track "<< it << "does not have image! " << endl;          cout << "TrkTrack *TrkLevel2::GetTrackImage(int) >> Track "<< it << "does not have image! " << endl;
1816          return 0;          return 0;
1817      }      }
1818      if(!Track)return 0;      if(!Track)return 0;
# Line 1351  void TrkLevel2::LoadField(TString path){ Line 1839  void TrkLevel2::LoadField(TString path){
1839  //     path_.error   = 0;  //     path_.error   = 0;
1840  //     readb_();  //     readb_();
1841    
1842      TrkParams::SetTrackingMode();  //     TrkParams::SetTrackingMode();
1843      TrkParams::SetPrecisionFactor();  //     TrkParams::SetPrecisionFactor();
1844      TrkParams::SetStepMin();  //     TrkParams::SetStepMin();
1845        TrkParams::SetMiniDefault();
1846    
1847      TrkParams::Set(path,1);      TrkParams::Set(path,1);
1848      TrkParams::Load(1);      TrkParams::Load(1);
Legend:
Removed from v.1.33  
changed lines
  Added in v.1.43
  ViewVC Help
Powered by ViewVC 1.1.23