/[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.54 by pam-fi, Wed Mar 11 14:19:10 2009 UTC revision 1.60 by pam-fi, Thu Jan 28 14:38:24 2016 UTC
# Line 381  Int_t TrkTrack::GetLeverArmXY(){ Line 381  Int_t TrkTrack::GetLeverArmXY(){
381      int first_plane = -1;      int first_plane = -1;
382      int last_plane  = -1;      int last_plane  = -1;
383      for(Int_t ip=0; ip<6; ip++){      for(Int_t ip=0; ip<6; ip++){
384          if( XGood(ip)*YGood(ip) && first_plane == -1 )first_plane = ip;          if( XGood(ip) && YGood(ip) && first_plane == -1 )first_plane = ip;
385          if( XGood(ip)*YGood(ip) && first_plane != -1 )last_plane = ip;          if( XGood(ip) && YGood(ip) && first_plane != -1 )last_plane = ip;
386      }      }
387      if( first_plane == -1 || last_plane == -1){      if( first_plane == -1 || last_plane == -1){
388          cout<< "Int_t TrkTrack::GetLeverArmXY() -- XGood(ip)*YGood(ip) always false ??? "<<endl;          cout<< "Int_t TrkTrack::GetLeverArmXY() -- XGood(ip)*YGood(ip) always false ??? "<<endl;
# Line 391  Int_t TrkTrack::GetLeverArmXY(){ Line 391  Int_t TrkTrack::GetLeverArmXY(){
391      return (last_plane-first_plane+1);      return (last_plane-first_plane+1);
392  }  }
393  /**  /**
394     * Returns the number of hit planes
395     */
396    Int_t TrkTrack::GetNhit()  {
397      int np=0;
398      for(Int_t ip=0; ip<6; ip++) np += (XGood(ip)||YGood(ip)) ;
399      return np;
400    };
401    /**
402   * Returns the reduced chi-square of track x-projection   * Returns the reduced chi-square of track x-projection
403   */   */
404  Float_t  TrkTrack::GetChi2X(){  Float_t  TrkTrack::GetChi2X(){
# Line 398  Float_t  TrkTrack::GetChi2X(){ Line 406  Float_t  TrkTrack::GetChi2X(){
406      for(int ip=0; ip<6; ip++)if(XGood(ip))chiq+= pow((xv[ip]-xm[ip])/resx[ip],2.);      for(int ip=0; ip<6; ip++)if(XGood(ip))chiq+= pow((xv[ip]-xm[ip])/resx[ip],2.);
407      if(GetNX()>3)chiq=chiq/(GetNX()-3);      if(GetNX()>3)chiq=chiq/(GetNX()-3);
408      else chiq=0;      else chiq=0;
409      if(chiq==0)cout << " Float_t  TrkTrack::GetChi2X() -- WARNING -- value not defined "<<chiq<<endl;      //    if(chiq==0)cout << " Float_t  TrkTrack::GetChi2X() -- WARNING -- value not defined "<<chiq<<endl;
410      return chiq;      return chiq;
411  }  }
412  /**  /**
# Line 409  Float_t  TrkTrack::GetChi2Y(){ Line 417  Float_t  TrkTrack::GetChi2Y(){
417      for(int ip=0; ip<6; ip++)if(YGood(ip))chiq+= pow((yv[ip]-ym[ip])/resy[ip],2.);      for(int ip=0; ip<6; ip++)if(YGood(ip))chiq+= pow((yv[ip]-ym[ip])/resy[ip],2.);
418      if(GetNY()>2)chiq=chiq/(GetNY()-2);      if(GetNY()>2)chiq=chiq/(GetNY()-2);
419      else chiq=0;      else chiq=0;
420      if(chiq==0)cout << " Float_t  TrkTrack::GetChi2Y() -- WARNING -- value not defined "<<chiq<<endl;      //    if(chiq==0)cout << " Float_t  TrkTrack::GetChi2Y() -- WARNING -- value not defined "<<chiq<<endl;
421      return chiq;      return chiq;
422  }  }
423  /**  /**
# Line 664  void TrkTrack::FillMiniStruct(cMini2trac Line 672  void TrkTrack::FillMiniStruct(cMini2trac
672    
673  //      cout << i<<" - "<<xgood[i]<<" "<<XGood(i)<<endl;  //      cout << i<<" - "<<xgood[i]<<" "<<XGood(i)<<endl;
674  //      cout << i<<" - "<<ygood[i]<<" "<<YGood(i)<<endl;  //      cout << i<<" - "<<ygood[i]<<" "<<YGood(i)<<endl;
675          track.xgood[i]=XGood(i);        track.xgood[i] = (double) XGood(i);
676          track.ygood[i]=YGood(i);          track.ygood[i] = (double) YGood(i);
677                    
678          track.xm[i]=xm[i];          track.xm[i]= (double) xm[i];
679          track.ym[i]=ym[i];          track.ym[i]= (double) ym[i];
680          track.zm[i]=zm[i];          track.zm[i]= (double) zm[i];
681                    
682  //      --- temporaneo ----------------------------  //      --- temporaneo ----------------------------
683  //      float segment = 100.;  //      float segment = 100.;
# Line 687  void TrkTrack::FillMiniStruct(cMini2trac Line 695  void TrkTrack::FillMiniStruct(cMini2trac
695  //      --- temporaneo ----------------------------  //      --- temporaneo ----------------------------
696    
697          if( XGood(i) || YGood(i) ){          if( XGood(i) || YGood(i) ){
698              double segment = 2.;//cm              //NB!! the length of the sensor is not exactely taken into account    
699                double segment = 7.;// 2.;//cm //Elena 10th
700              // NB: i parametri di allineamento hanno una notazione particolare!!!              // NB: i parametri di allineamento hanno una notazione particolare!!!
701              // sensor = 0 (hybrid side), 1              // sensor = 0 (hybrid side), 1
702              // ladder = 0-2 (increasing x)              // ladder = 0-2 (increasing x)
# Line 697  void TrkTrack::FillMiniStruct(cMini2trac Line 706  void TrkTrack::FillMiniStruct(cMini2trac
706              int il = (int)GetLadder(i);              int il = (int)GetLadder(i);
707                            
708              double omega   = 0.;              double omega   = 0.;
709              double beta    = 0.;              //      double beta    = 0.;// EM GCC 4.7
710              double gamma   = 0.;              //      double gamma   = 0.;
711              if(              if(
712                  (is < 0 || is > 1 || ip < 0 || ip > 5 || il < 0 || il > 2) &&                  (is < 0 || is > 1 || ip < 0 || ip > 5 || il < 0 || il > 2) &&
713                  true){                  true){
# Line 708  void TrkTrack::FillMiniStruct(cMini2trac Line 717  void TrkTrack::FillMiniStruct(cMini2trac
717                  cout << " is ip il = "<<is<<" "<<ip<<" "<<il<<endl;                  cout << " is ip il = "<<is<<" "<<ip<<" "<<il<<endl;
718              }else{              }else{
719                  omega   = alignparameters_.omega[is][il][ip];                  omega   = alignparameters_.omega[is][il][ip];
720                  beta    = alignparameters_.beta[is][il][ip];                  //              beta    = alignparameters_.beta[is][il][ip];// EM GCC 4.7 unused
721                  gamma   = alignparameters_.gamma[is][il][ip];                  //              gamma   = alignparameters_.gamma[is][il][ip];// EM GCC 4.7 unused
722              }              }
723                            
724              if(       XGood(i) && !YGood(i) ){              if(       XGood(i) && !YGood(i) ){
725                  track.xm_a[i] = xm[i] - omega * segment;                  track.xm_a[i] =  (double) xm[i] - omega * segment;
726                  track.ym_a[i] = ym[i] + segment;                  track.ym_a[i] =  (double) ym[i] + segment;
727  //          track.zm_a[i] = zm[i] + beta * segment;//not used yet  //          track.zm_a[i] = zm[i] + beta * segment;//not used yet
728                  track.xm_b[i] = xm[i] + omega * segment;                  track.xm_b[i] =  (double) xm[i] + omega * segment;
729                  track.ym_b[i] = ym[i] - segment;                  track.ym_b[i] =  (double) ym[i] - segment;
730  //          track.zm_b[i] = zm[i] - beta * segment;//not used yet  //          track.zm_b[i] = zm[i] - beta * segment;//not used yet
731              }else if( !XGood(i) && YGood(i) ){              }else if( !XGood(i) && YGood(i) ){
732                  track.xm_a[i] = xm[i] + segment;                  track.xm_a[i] =  (double) xm[i] + segment;
733                  track.ym_a[i] = ym[i] + omega * segment;                  track.ym_a[i] =  (double) ym[i] + omega * segment;
734  //          track.zm_a[i] = zm[i] - gamma * segment;//not used yet  //          track.zm_a[i] = zm[i] - gamma * segment;//not used yet
735                  track.xm_b[i] = xm[i] - segment;                  track.xm_b[i] =  (double) xm[i] - segment;
736                  track.ym_b[i] = ym[i] - omega * segment;                  track.ym_b[i] =  (double) ym[i] - omega * segment;
737  //          track.zm_b[i] = zm[i] + gamma * segment;//not used yet  //          track.zm_b[i] = zm[i] + gamma * segment;//not used yet
738              }              }
739          }          }
740                    
741          track.resx[i]=resx[i];          track.resx[i]= (double) resx[i];
742          track.resy[i]=resy[i];          track.resy[i]= (double) resy[i];
743          track.tailx[i]=tailx[i];          track.tailx[i]= (double) tailx[i];
744          track.taily[i]=taily[i];          track.taily[i]= (double) taily[i];
745      }      }
746    
747      for(int i=0; i<5; i++) track.al[i]=al[i];      for(int i=0; i<5; i++) track.al[i]= (double) al[i];
748      track.zini = 23.5;      track.zini = 23.5;
749  // ZINI = 23.5 !!! it should be the same parameter in all codes  // ZINI = 23.5 !!! it should be the same parameter in all codes
750            
# Line 746  void TrkTrack::FillMiniStruct(cMini2trac Line 755  void TrkTrack::FillMiniStruct(cMini2trac
755  void TrkTrack::SetFromMiniStruct(cMini2track *track){  void TrkTrack::SetFromMiniStruct(cMini2track *track){
756    
757      for(int i=0; i<5; i++) {      for(int i=0; i<5; i++) {
758          al[i]=track->al[i];        al[i]= (float) (track->al[i]);
759          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]);
760      }      }
761      chi2  = track->chi2;      chi2  =  (float) (track->chi2);
762      nstep = track->nstep;      nstep =  (float) (track->nstep);
763      for(int i=0; i<6; i++){      for(int i=0; i<6; i++){
764          xv[i]  = track->xv[i];          xv[i]  =  (float) (track->xv[i]);
765          yv[i]  = track->yv[i];          yv[i]  =  (float) (track->yv[i]);
766          zv[i]  = track->zv[i];          zv[i]  =  (float) (track->zv[i]);
767          xm[i]  = track->xm[i];          xm[i]  = (float)  (track->xm[i]);
768          ym[i]  = track->ym[i];          ym[i]  =  (float) (track->ym[i]);
769          zm[i]  = track->zm[i];          zm[i]  =  (float) (track->zm[i]);
770          axv[i] = track->axv[i];          axv[i] =  (float) (track->axv[i]);
771          ayv[i] = track->ayv[i];          ayv[i] =  (float) (track->ayv[i]);      
772            resx[i] = (float)  (track->resx[i]); //Elena 10th
773            resy[i] =  (float) (track->resy[i]);
774      }      }
775            
776  }  }
# Line 810  Bool_t TrkTrack::EvaluateClusterPosition Line 821  Bool_t TrkTrack::EvaluateClusterPosition
821    
822      for(int ip=0; ip<6; ip++){      for(int ip=0; ip<6; ip++){
823  //      cout << ip<<" ** "<<xm[ip]<<" / "<<ym[ip]<<endl;;  //      cout << ip<<" ** "<<xm[ip]<<" / "<<ym[ip]<<endl;;
824          int icx = GetClusterX_ID(ip)+1;          int icx = GetClusterX_ID(ip)+1;//0=no-cluster,1-N
825          int icy = GetClusterY_ID(ip)+1;          int icy = GetClusterY_ID(ip)+1;//0=no-cluster,1-N
826          int sensor = GetSensor(ip)+1;//<< convenzione "Paolo"          int sensor = GetSensor(ip)+1;//<< convenzione "Paolo"
827          if(ip==5 && sensor!=0)sensor=3-sensor;//<< convenzione "Elena"          if(ip==5 && sensor!=0)sensor=3-sensor;//<< convenzione "Elena"
828          int ladder = GetLadder(ip)+1;          int ladder = GetLadder(ip)+1;
# Line 825  Bool_t TrkTrack::EvaluateClusterPosition Line 836  Bool_t TrkTrack::EvaluateClusterPosition
836          float bfy = 10*TrkParams::GetBY(v);//Tesla          float bfy = 10*TrkParams::GetBY(v);//Tesla
837          int ipp=ip+1;          int ipp=ip+1;
838          xyzpam_(&ipp,&icx,&icy,&ladder,&sensor,&ax,&ay,&bfx,&bfy);          xyzpam_(&ipp,&icx,&icy,&ladder,&sensor,&ax,&ay,&bfx,&bfy);
839          if(icx<0 || icy<0)return false;          //      if(icx<0 || icy<0)return false;
840      }      }
841      return true;      return true;
842  }  }
# Line 875  void TrkTrack::Fit(double pfixed, int& f Line 886  void TrkTrack::Fit(double pfixed, int& f
886      extern cMini2track track_;      extern cMini2track track_;
887      fail = 0;      fail = 0;
888    
889      FillMiniStruct(track_);      //    FillMiniStruct(track_);
890                    
891      if(froml1!=0){      if(froml1!=0){
892          if( !EvaluateClusterPositions() ){          if( !EvaluateClusterPositions() ){
# Line 914  void TrkTrack::Fit(double pfixed, int& f Line 925  void TrkTrack::Fit(double pfixed, int& f
925          if(iprint)cout << "ERROR: ifail= " << ifail << endl;          if(iprint)cout << "ERROR: ifail= " << ifail << endl;
926          fail = 1;          fail = 1;
927      }      }
928        if(chi2!=chi2){
929            if(iprint)cout << "ERROR: chi2= " << chi2 << endl;      
930            FitReset();
931            fail = 1;      
932        }
933      //  ------------------------------------------      //  ------------------------------------------
934            
935      SetFromMiniStruct(&track_);      SetFromMiniStruct(&track_);
# Line 923  void TrkTrack::Fit(double pfixed, int& f Line 939  void TrkTrack::Fit(double pfixed, int& f
939          for(int i=0; i<5; i++) al[i]=al_ini[i];          for(int i=0; i<5; i++) al[i]=al_ini[i];
940      }      }
941    
942    
943        //ELENA 2015
944        int ngf = TrkParams::nGF;
945        float* zgf =  TrkParams::zGF;    
946        Trajectory tj = Trajectory(ngf,zgf);
947        tj.DoTrack(al,ZINI);
948        for(int i=0; i<14; i++){
949            xGF[i] = tj.x[i];
950            yGF[i] = tj.y[i];
951        }
952        if(ngf!=14){
953            cout << "TrkTrack::Fit("<<pfixed<<","<<fail<<","<<iprint<<","<<froml1<<"<<) "<<endl;;
954            cout << "TrkParams::nGF = "<<TrkParams::nGF<<" != 14 "<<endl;
955            cout << "back-incompatibility !?!?! acceptance check not reliable"<<endl;
956        }
957        
958        
959    
960  };  };
961  /**  /**
962   * Reset the fit parameters   * Reset the fit parameters
# Line 1152  Int_t TrkTrack::GetSensor(int ip){ Line 1186  Int_t TrkTrack::GetSensor(int ip){
1186  void TrkTrack::SetXGood(int ip, int clid, int il, int is, bool bad){  void TrkTrack::SetXGood(int ip, int clid, int il, int is, bool bad){
1187  //    int il=0;       //ladder (temporary)  //    int il=0;       //ladder (temporary)
1188  //    bool bad=false; //ladder (temporary)  //    bool bad=false; //ladder (temporary)
1189      if(ip<0||ip>5||clid<0||il<-1||il>2||is<-1||is>1)      if(ip<0||ip>5||clid<1||il<-1||il>2||is<-1||is>1)
1190          cout << " void TrkTrack::SetXGood(int,int,int,int,bool) --> MA SEI DI COCCIO?!?!"<<endl;          cout << " void TrkTrack::SetXGood(int,int,int,int,bool) --> MA SEI DI COCCIO?!?!"<<endl;
1191      xgood[ip]=(il+1)*100000000+(is+1)*10000000+clid;      xgood[ip]=(il+1)*100000000+(is+1)*10000000+clid;
1192      if(bad)xgood[ip]=-xgood[ip];      if(bad)xgood[ip]=-xgood[ip];
# Line 1169  void TrkTrack::SetXGood(int ip, int clid Line 1203  void TrkTrack::SetXGood(int ip, int clid
1203  void TrkTrack::SetYGood(int ip, int clid, int il, int is, bool bad){  void TrkTrack::SetYGood(int ip, int clid, int il, int is, bool bad){
1204  //    int il=0;       //ladder (temporary)  //    int il=0;       //ladder (temporary)
1205  //    bool bad=false; //ladder (temporary)  //    bool bad=false; //ladder (temporary)
1206      if(ip<0||ip>5||clid<0||il<-1||il>2||is<-1||is>1)      if(ip<0||ip>5||clid<1||il<-1||il>2||is<-1||is>1)
1207          cout << " void TrkTrack::SetYGood(int,int,int,int,bool) --> MA SEI DI COCCIO?!?!"<<endl;          cout << " void TrkTrack::SetYGood(int,int,int,int,bool) --> MA SEI DI COCCIO?!?!"<<endl;
1208      ygood[ip]=(il+1)*100000000+(is+1)*10000000+clid;      ygood[ip]=(il+1)*100000000+(is+1)*10000000+clid;
1209      if(bad)ygood[ip]=-ygood[ip];      if(bad)ygood[ip]=-ygood[ip];
# Line 1568  void TrkLevel2::Set(){ Line 1602  void TrkLevel2::Set(){
1602  //  //
1603  //  //
1604  //--------------------------------------  //--------------------------------------
1605    void TrkLevel2::SetTrackArray(TClonesArray *track){
1606        if(track && strcmp(track->GetClass()->GetName(),"TrkTrack")==0){
1607            if(Track)Track->Clear("C");    
1608            Track = track;
1609        }
1610    }
1611    //--------------------------------------
1612    //
1613    //
1614    //--------------------------------------
1615  void TrkLevel2::Dump(){  void TrkLevel2::Dump(){
1616                    
1617          //          //
# Line 2178  Float_t TrkLevel2::GetZTrk(Int_t plane_i Line 2222  Float_t TrkLevel2::GetZTrk(Int_t plane_i
2222   * (By default is created with z-coordinates inside the tracking volume)   * (By default is created with z-coordinates inside the tracking volume)
2223    */    */
2224  Trajectory::Trajectory(){  Trajectory::Trajectory(){
2225      npoint = 10;      npoint = 6;
2226      x = new float[npoint];      x = new float[npoint];
2227      y = new float[npoint];      y = new float[npoint];
2228      z = new float[npoint];      z = new float[npoint];
# Line 2245  Trajectory::Trajectory(int n, float* zin Line 2289  Trajectory::Trajectory(int n, float* zin
2289      thy = new float[npoint];      thy = new float[npoint];
2290      tl = new float[npoint];      tl = new float[npoint];
2291      int i=0;      int i=0;
2292      do{      do{      
2293          x[i] = 0;          x[i] = 0.;
2294          y[i] = 0;          y[i] = 0.;
2295          z[i] = zin[i];          z[i] = zin[i];
2296          thx[i] = 0;          thx[i] = 0.;
2297          thy[i] = 0;          thy[i] = 0.;
2298          tl[i] = 0;          tl[i] = 0.;
2299          i++;                      i++;            
2300      }while(zin[i-1] > zin[i] && i < npoint);      }while(zin[i-1] > zin[i] && i < npoint);
2301      npoint=i;      npoint=i;
2302      if(npoint != n)cout << "NB! Trajectory created with "<<npoint<<" points"<<endl;      if(npoint != n)cout << "NB! Trajectory created with "<<npoint<<" points instean of "<<n<<endl;
2303        //    Dump();
2304  }  }
2305  void Trajectory::Delete(){  void Trajectory::Delete(){
2306            
Legend:
Removed from v.1.54  
changed lines
  Added in v.1.60
  ViewVC Help
Powered by ViewVC 1.1.23