TrackerLevel2/src/ExtTrack.cpp

/**
 * \file ExtTrack.cpp
 * \author Elena Vannuccini
 */
#include <ExtTrack.h>
#include <iostream>
#include <math.h>
using namespace std;

ExtTrack::ExtTrack(int dim){

  //cout << " ExtTrack::ExtTrack("<<dim<<")"<<endl;
  xgood   = NULL;
  ygood   = NULL;
  multmaxx   = NULL;
  multmaxy   = NULL;
  xm      = NULL;
  ym      = NULL;
  zm      = NULL;
  xma     = NULL;
  yma     = NULL;
  zma     = NULL;
  xmb     = NULL;
  ymb     = NULL;
  zmb     = NULL;
  resx    = NULL;
  resy    = NULL;
  xv      = NULL;
  yv      = NULL;
  zv      = NULL;
  axv     = NULL;
  ayv     = NULL;
  dedx_x  = NULL;
  dedx_y  = NULL;

  // xGF = new float[TrkParams::nGF];
  // yGF = new float[TrkParams::nGF];
  // init vectors  
  SetDimension(dim);//allocate vectors
  // init fit variables
  Reset();
  // set fit parameters (default)
  SetZ0(23.5);        
  SetMiniDefault();  
};
//--------------------------------------
//
//
//--------------------------------------
ExtTrack::ExtTrack( const ExtTrack& t ){
  //cout << " ExtTrack::ExtTrack( const ExtTrack& t ) "<<this<<endl;

  // --------------------
  // initialize pointers
  // --------------------

  xgood   = NULL;
  ygood   = NULL;
  multmaxx   = NULL;
  multmaxy   = NULL;
  xm      = NULL;
  ym      = NULL;
  zm      = NULL;
  xma     = NULL;
  yma     = NULL;
  zma     = NULL;
  xmb     = NULL;
  ymb     = NULL;
  zmb     = NULL;
  resx    = NULL;
  resy    = NULL;
  xv      = NULL;
  yv      = NULL;
  zv      = NULL;
  axv     = NULL;
  ayv     = NULL;
  dedx_x  = NULL;
  dedx_y  = NULL;

  // xGF = new float[TrkParams::nGF];
  // yGF = new float[TrkParams::nGF];
  
  SetDimension(t.nplanes);//allocate vectors

  // --------------------
  // copy values
  // --------------------
  zini  = t.zini;
  chi2  = t.chi2;
  nstep = t.nstep;
  for(int it1=0;it1<5;it1++){
    al[it1] = t.al[it1];
    for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
  };


  for(int ip=0;ip<nplanes;ip++){
    xgood[ip]  = t.xgood[ip];
    ygood[ip]  = t.ygood[ip];
    multmaxx[ip]  = t.multmaxx[ip];
    multmaxy[ip]  = t.multmaxy[ip];
    xm[ip]     = t.xm[ip];
    ym[ip]     = t.ym[ip];
    zm[ip]     = t.zm[ip];
    xma[ip]     = t.xma[ip];
    yma[ip]     = t.yma[ip];
    zma[ip]     = t.zma[ip];
    xmb[ip]     = t.xmb[ip];
    ymb[ip]     = t.ymb[ip];
    zmb[ip]     = t.zmb[ip];
    resx[ip]   = t.resx[ip];
    resy[ip]   = t.resy[ip];
    xv[ip]     = t.xv[ip];
    yv[ip]     = t.yv[ip];
    zv[ip]     = t.zv[ip];
    axv[ip]    = t.axv[ip];
    ayv[ip]    = t.ayv[ip];
    dedx_x[ip] = t.dedx_x[ip];
    dedx_y[ip] = t.dedx_y[ip];
  };
  
  int ngf = TrkParams::nGF;
  for(int i=0; i<ngf; i++){
    xGF[i] = t.xGF[i];
    yGF[i] = t.yGF[i];
  }
  
};
//--------------------------------------
//
//
//--------------------------------------
void ExtTrack::Copy( ExtTrack& t ){

  //  cout << " ExtTrack::Copy( ExtTrack& t ) "<<endl;

  t.Reset();//reset variables
  t.Clear("C");//deallocate vectors
  t.SetDimension(nplanes);//allocate vectors

  t.zini  = zini;
  t.chi2  = chi2;
  t.nstep = nstep;
  for(int it1=0;it1<5;it1++){
    t.al[it1] = al[it1];
    for(int it2=0;it2<5;it2++)t.coval[it1][it2] = coval[it1][it2];
  };    
//  cout << "ExtTrack::Copy()........nplanes "<<nplanes<<endl;
  for(int ip=0;ip<nplanes;ip++){
    t.xgood[ip]  = xgood[ip];
    t.ygood[ip]  = ygood[ip];
    t.multmaxx[ip]  = multmaxx[ip];
    t.multmaxy[ip]  = multmaxy[ip];
    t.xm[ip]     = xm[ip];
    t.ym[ip]     = ym[ip];
    t.zm[ip]     = zm[ip];
    t.xma[ip]     = xma[ip];
    t.yma[ip]     = yma[ip];
    t.zma[ip]     = zma[ip];
    t.xmb[ip]     = xmb[ip];
    t.ymb[ip]     = ymb[ip];
    t.zmb[ip]     = zmb[ip];
    t.resx[ip]   = resx[ip];
    t.resy[ip]   = resy[ip];
    t.xv[ip]     = xv[ip];
    t.yv[ip]     = yv[ip];
    t.zv[ip]     = zv[ip];
    t.axv[ip]    = axv[ip];
    t.ayv[ip]    = ayv[ip];
    t.dedx_x[ip] = dedx_x[ip];
    t.dedx_y[ip] = dedx_y[ip];
  };
  int ngf = TrkParams::nGF;
  for(int i=0; i<ngf; i++){
    t.xGF[i] = xGF[i];
    t.yGF[i] = yGF[i];
  }
  
};
//--------------------------------------
//
//
//--------------------------------------
void ExtTrack::Set( TrkTrack& t , int index){

  //  cout << "Set("<< &t <<", "<<index<<")"<<endl;
  if(index < 0 ){ //NON FUNZIONA
    cout << "Clear(\"C\") "<<endl; 
    Reset();//reset variables
    Clear("C");//deallocate vectors
    cout << "SetDimension("<<t.GetNhit()<<")"<<endl; 
    SetDimension(t.GetNhit()); //allocate vectors
//   }else{
//       SetDimension(6+index);
  }
  // default fit settings (the same as TrkTrack)
  SetMiniDefault();
  SetZ0(23.5); 
  //
  int ipp;
  if(index<0) ipp = 0;
  else        ipp = index;
  //
//  cout << "Assign... nplanes = "<<nplanes<<endl;
  for(int ip=0;ip<6;ip++){
    
    if(ipp>=nplanes){
      cout << "void ExtTrack::Set( TrkTrack& , "<<index<<") -- ipp="<<ipp<<" exceed vector dimention "<<nplanes<<endl;
      continue;
    }

    if( t.XGood(ip) || t.YGood(ip) || !(index<0) ){

//      cout << "ipp "<<ipp<<" ip "<<ip<<endl;


      xgood[ipp]  = t.XGood(ip); //NB! 
      ygood[ipp]  = t.YGood(ip); //NB! 
      multmaxx[ipp]  = t.multmaxx[ip]; 
      multmaxy[ipp]  = t.multmaxy[ip]; 
      xm[ipp]     = t.xm[ip];
      ym[ipp]     = t.ym[ip];
      zm[ipp]     = t.zm[ip];

      float dummy = -100.;
      xma[ipp]     = dummy;
      yma[ipp]     = dummy;
      zma[ipp]     = zm[ipp];
      xmb[ipp]     = dummy;
      ymb[ipp]     = dummy;
      zmb[ipp]     = zm[ipp];

      if( t.XGood(ip) && t.YGood(ip) ){ //double hit
        xma[ipp]     = t.xm[ip];
        yma[ipp]     = t.ym[ip];
        zma[ipp]     = t.zm[ip];
        xmb[ipp]     = t.xm[ip];
        ymb[ipp]     = t.ym[ip];
        zmb[ipp]     = t.zm[ip];
      } else if (t.XGood(ip) || t.YGood(ip) ){                         //single hit
        TrkParams::Load(5); 
        if( !TrkParams::IsLoaded(5) ){
          cout << "void ExtTrack::FillMiniStruct(TrkTrack&) ---ERROR--- trk align.param. not loaded "<<endl;
          return;
        }
        double segment = 7.;// 2.;//cm //Elena 10th
        int is = (int)t.GetSensor(ip); if(ip==5)is=1-is;
        int ippp = 5-ip;
        int il = (int)t.GetLadder(ip);      
        // NB: i parametri di allineamento hanno una notazione particolare!!!
        // sensor = 0 (hybrid side), 1 
        // ladder = 0-2 (increasing x)
        // plane  = 0-5 (from bottom to top!!!) 
        double omega   = 0.; 
        double beta   = 0.; 
        double gamma   = 0.; 
        if( 
           (is < 0 || is > 1 || ip < 0 || ippp > 5 || il < 0 || il > 2) &&
           true){
          cout << " void ExtTrack::FillMiniStruct(TrkTrack&) --- WARNING ---trk sensor not defined, cannot read alignment parameters "<<endl;
          cout << " is ip il = "<<is<<" "<<ippp<<" "<<il<<endl;
        }else{
          omega   = alignparameters_.omega[is][il][ippp];
          beta    = alignparameters_.beta[is][il][ippp];
          gamma   = alignparameters_.gamma[is][il][ippp];
        }         
        if(       t.XGood(ip) && !t.YGood(ip) ){
          xma[ipp] = t.xm[ip] - omega * segment;
          yma[ipp] = t.ym[ip] + segment;
          zma[ipp] = t.zm[ip] + beta * segment;//not used yet
          xmb[ipp] = t.xm[ip] + omega * segment;
          ymb[ipp] = t.ym[ip] - segment;
          zmb[ipp] = t.zm[ip] - beta * segment;//not used yet
        }else if( !t.XGood(ip) && t.YGood(ip) ){
          xma[ipp] = t.xm[ip] + segment;
          yma[ipp] = t.ym[ip] + omega * segment;
          zma[ipp] = t.zm[ip] - gamma * segment;//not used yet
          xmb[ipp] = t.xm[ip] - segment;
          ymb[ipp] = t.ym[ip] - omega * segment;
          zmb[ipp] = t.zm[ip] + gamma * segment;//not used yet
        }
      }
        
      // tailx[ipp]  = t.tailx[ip];
      // taily[ipp]  = t.taily[ip];
      resx[ipp]   = t.resx[ip];
      resy[ipp]   = t.resy[ip];
      xv[ipp]     = t.xv[ip];
      yv[ipp]     = t.yv[ip];
      zv[ipp]     = t.zv[ip];
      axv[ipp]    = t.axv[ip];
      ayv[ipp]    = t.ayv[ip];
      dedx_x[ipp] = t.dedx_x[ip];
      dedx_y[ipp] = t.dedx_y[ip];

      ipp++;

    }
  };//endl loop on planes

  chi2  = t.chi2;
  nstep = t.nstep;
  for(int it1=0;it1<5;it1++){
    al[it1] = t.al[it1];
    for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
  }; 
 
};
//--------------------------------------
//
//
//--------------------------------------
void ExtTrack::SetMiniDefault(){
  //  cout << " ExtTrack::SetMiniDefault()"<<endl;
  exttrack_.trackmode = TrkParams::init__mini_trackmode;
  exttrack_.fact      = TrkParams::init__mini_fact;   
  exttrack_.istepmin  = TrkParams::init__mini_istepmin; 
  TrkParams::SetDeltaB();
  TrkParams::SetDLT();
}
//--------------------------------------
//
//
//--------------------------------------
void ExtTrack::SetDimension(int dim){

  //cout << " ExtTrack::SetDimension("<<dim<<")"<<endl;;
  if(dim<=0){
    nplanes = 0;
    if(dim<0)cout << "void ExtTrack::SetDimension("<<dim<<") --- invalid dimension "<<endl;
    return;
  }

  //  Clear();

  nplanes = dim;
  
  xgood   = new int[nplanes];
  ygood   = new int[nplanes];
  multmaxx   = new int[nplanes];
  multmaxy   = new int[nplanes];
  xm      = new float[nplanes];
  ym      = new float[nplanes];
  zm      = new float[nplanes];
  xma     = new float[nplanes];
  yma     = new float[nplanes];
  zma     = new float[nplanes];
  xmb     = new float[nplanes];
  ymb     = new float[nplanes];
  zmb     = new float[nplanes];
  resx    = new float[nplanes];
  resy    = new float[nplanes];
  xv      = new float[nplanes];
  yv      = new float[nplanes];
  zv      = new float[nplanes];
  axv     = new float[nplanes];
  ayv     = new float[nplanes];
  dedx_x  = new float[nplanes];
  dedx_y  = new float[nplanes];
  

}
//--------------------------------------
//
//
//--------------------------------------
bool ExtTrack::SetZ(int ip,float zmeas){
  if(ip<0 || ip>=nplanes){
    cout << " ExtTrack::SetZ("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
    return false;
  }
  zm[ip] = zmeas;
  zma[ip] = zmeas;
  zmb[ip] = zmeas;
  return true;

};

bool ExtTrack::ResetXY(int ip){

  // cout << " ExtTrack::ResetXY("<<ip<<") "<<endl;

  if(ip<0 || ip>=nplanes){
    cout << " ExtTrack::ResetXY("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
    return false;
  } 
  xm[ip] = -100.;
  ym[ip] = -100.;
  xma[ip] = -100.;
  yma[ip] = -100.;
  xmb[ip] = -100.;
  ymb[ip] = -100.;
  resx[ip] = 1000.;
  resy[ip] = 1000.;
  dedx_x[ip] = 0.;
  dedx_y[ip] = 0.;
  xgood[ip] = 0;
  ygood[ip] = 0;
  multmaxx[ip] = 0;
  multmaxy[ip] = 0;

  return true;

};

bool ExtTrack::SetXY(int ip,float xmeas, float ymeas, float rx, float ry){

  if(ip<0 || ip>=nplanes){
    cout << " ExtTrack::SetXY("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
    return false;
  }
  
  xm[ip] = xmeas;
  ym[ip] = ymeas;
  xma[ip] = xmeas;
  yma[ip] = ymeas;
  xmb[ip] = xmeas;
  ymb[ip] = ymeas;
  resx[ip] = rx;
  resy[ip] = ry;
  // dedx_x[ip] = dedxx;
  // dedx_y[ip] = dedxy;
  xgood[ip] = 1;
  ygood[ip] = 1;

  return true;

}; 
bool ExtTrack::SetX(int ip,float xa, float xb, float ya, float yb,float res ){

  // cout << " ExtTrack::SetX("<<ip<<",...)  -- nplanes"<<nplanes<<endl;

  if(ip<0 || ip>=nplanes){
    cout << " ExtTrack::SetX("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
    return false;
  }
  xm[ip] = -100.;
  ym[ip] = -100.;
  xma[ip] = xa;
  yma[ip] = ya;
  xmb[ip] = xb;
  ymb[ip] = yb;
  resx[ip] = res;
  resy[ip] = 1000.;
  xgood[ip] = 1;
  ygood[ip] = 0;
  // dedx_x[ip] = dedx;
  // dedx_y[ip] = 0.;
  return true;
}; 
bool ExtTrack::SetY(int ip,float xa, float xb, float ya, float yb,float res ){

  // cout << " ExtTrack::SetY("<<ip<<",...)  -- nplanes"<<nplanes<<endl;

  if(ip<0 || ip>=nplanes){
    cout << " ExtTrack::SetX("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
    return false;
  }
  xm[ip] = -100.;
  ym[ip] = -100.;
  xma[ip] = xa;
  yma[ip] = ya;
  xmb[ip] = xb;
  ymb[ip] = yb;
  resx[ip] = 1000.;
  resy[ip] = res;
  xgood[ip] = 0;
  ygood[ip] = 1;
  // dedx_x[ip] = 0.;
  // dedx_y[ip] = dedx;
  return true;
}; 

bool ExtTrack::SetXGood(int ip,int icl_piu_uno, int il, int is ){

  // cout << " ExtTrack::SetXGood("<<ip<<",...)  -- nplanes"<<nplanes<<endl;

  if(ip<0 || ip>=nplanes){
    cout << " ExtTrack::SetXGood("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
    return false;
  }
  xgood[ip]=(il+1)*100000000+(is+1)*10000000+icl_piu_uno; 
  return true;
}; 
bool ExtTrack::SetYGood(int ip,int icl_piu_uno, int il, int is ){ 
  if(ip<0 || ip>=nplanes){
    cout << " ExtTrack::SetYGood("<<ip<<",...)  -- nplanes"<<nplanes<<endl;
    cout << " ExtTrack::SetYGood("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
    return false;
  }
  ygood[ip]=(il+1)*100000000+(is+1)*10000000+icl_piu_uno; 
  return true;
};

//--------------------------------------
//
//
//--------------------------------------
/**
 * Clear the track. If option C is specified, deallocates the hit vectors.
 */
void ExtTrack::Clear(Option_t* option){
    
//    cout << " ExtTrack::Clear( "<<option<<" )"<<this<<endl;
        
    Reset();

   if (option && option[0] == 'C') {

//       cout << " uccidi!! "<<endl;

       if(xgood) delete []      xgood;
       if(ygood) delete []      ygood;
       if(multmaxx) delete [] multmaxx;
       if(multmaxy) delete [] multmaxy;
       if(xm)    delete []      xm;
       if(ym)    delete []      ym;
       if(zm)    delete []      zm;
       if(xma)   delete []      xma;
       if(yma)   delete []      yma;
       if(zma)   delete []      zma;
       if(xmb)   delete []      xmb;
       if(ymb)   delete []      ymb;
       if(zmb)   delete []      zmb;
       if(resx)  delete []      resx;
       if(resy)  delete []      resy;
       if(xv)    delete []      xv;
       if(yv)    delete []      yv;
       if(zv)    delete []      zv;
       if(axv)   delete []      axv;
       if(ayv)   delete []      ayv;
       if(dedx_x)delete []      dedx_x;
       if(dedx_y)delete []      dedx_y;
       //}
//        delete []     xgood;
//        delete []     ygood;
//        delete [] multmaxx;
//        delete [] multmaxy;
//        delete []     xm;
//        delete []     ym;
//        delete []     zm;
//        delete []     xma;
//        delete []     yma;
//        delete []     zma;
//        delete []     xmb;
//        delete []     ymb;
//        delete []     zmb;
//        delete []     resx;
//        delete []     resy;
//        delete []     xv;
//        delete []     yv;
//        delete []     zv;
//        delete []     axv;
//        delete []     ayv;
//        delete []     dedx_x;
//        delete []     dedx_y;

       xgood   = NULL;
       ygood   = NULL;
       multmaxx   = NULL;
       multmaxy   = NULL;
       xm      = NULL;
       ym      = NULL;
       zm      = NULL;
       xma     = NULL;
       yma     = NULL;
       zma     = NULL;
       xmb     = NULL;
       ymb     = NULL;
       zmb     = NULL;
       resx    = NULL;
       resy    = NULL;
       xv      = NULL;
       yv      = NULL;
       zv      = NULL;
       axv     = NULL;
       ayv     = NULL;
       dedx_x  = NULL;
       dedx_y  = NULL;

       nplanes = 0;

   }

};

void ExtTrack::Delete(){ 
  Clear("C");
  // delete [] xGF; 
  // delete [] yGF; 
}        

//--------------------------------------
//
//
//--------------------------------------

// /**
//  * Set the  position measurements 
//  */
// void ExtTrack::SetMeasure(Double_t *xmeas, Double_t *ymeas, Double_t *zmeas){
//     for(int i=0; i<nplanes; i++) xm[i]= (xmeas? *xmeas++ : 0.);
//     for(int i=0; i<nplanes; i++) ym[i]= (ymeas? *ymeas++ : 0.);
//     for(int 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 i=0; i<nplanes; i++) resx[i]= (rx? *rx++ : 0.);
//     for(int i=0; i<nplanes; i++) resy[i]= (ry? *ry++ : 0.);
// }
// /**
//  * Set the TrkTrack good measurement 
//  */
// void ExtTrack::SetGood(int *xg, int *yg){

//     for(int i=0; i<nplanes; i++) xgood[i]=(xg? *xg++ : 0 );
//     for(int i=0; i<nplanes; i++) ygood[i]=(yg? *yg++ : 0 );
// }

//--------------------------------------
//
//
//--------------------------------------
void ExtTrack::Dump(){
    cout << endl << "========== Track " ;
    cout << endl << "zini     : "<< zini;
    cout << endl << "al       : "; for(int i=0; i<5; i++)cout << al[i] << " ";
    cout << endl << "chi^2    : "<< chi2;
    cout << endl << "n.step   : "<< nstep;
    cout << endl << "xgood    : "; for(int i=0; i<nplanes; i++)cout << XGood(i) ;
    cout << endl << "ygood    : "; for(int i=0; i<nplanes; i++)cout << YGood(i) ;
    cout << endl << "xm       : "; for(int i=0; i<nplanes; i++)cout << xm[i] << " ";
    cout << endl << "ym       : "; for(int i=0; i<nplanes; i++)cout << ym[i] << " ";
    cout << endl << "zm       : "; for(int i=0; i<nplanes; i++)cout << zm[i] << " ";
    cout << endl << "xma       : "; for(int i=0; i<nplanes; i++)cout << xma[i] << " ";
    cout << endl << "yma       : "; for(int i=0; i<nplanes; i++)cout << yma[i] << " ";
    cout << endl << "zma       : "; for(int i=0; i<nplanes; i++)cout << zma[i] << " ";
    cout << endl << "xmb       : "; for(int i=0; i<nplanes; i++)cout << xmb[i] << " ";
    cout << endl << "ymb       : "; for(int i=0; i<nplanes; i++)cout << ymb[i] << " ";
    cout << endl << "zmb       : "; for(int i=0; i<nplanes; i++)cout << zmb[i] << " ";
    cout << endl << "xv       : "; for(int i=0; i<nplanes; i++)cout << xv[i] << " ";
    cout << endl << "yv       : "; for(int i=0; i<nplanes; i++)cout << yv[i] << " ";
    cout << endl << "zv       : "; for(int i=0; i<nplanes; i++)cout << zv[i] << " ";
    cout << endl << "resx     : "; for(int i=0; i<nplanes; i++)cout << resx[i] << " ";
    cout << endl << "resy     : "; for(int i=0; i<nplanes; i++)cout << resy[i] << " ";
    cout << endl << "coval    : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" ";
    cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[1][i]<<" ";
    cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[2][i]<<" ";
    cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[3][i]<<" ";
    cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[4][i]<<" ";
    cout << endl << "MDR      : "<< (coval[4][4]>0. ?  1./sqrt(coval[4][4]) : 0. ) ;
    cout << endl << "dedx_x   : "; for(int i=0; i<nplanes; i++)cout << dedx_x[i] << " ";
    cout << endl << "dedx_y   : "; for(int i=0; i<nplanes; i++)cout << dedx_y[i] << " ";

    cout << endl << "=================";
    cout << endl;
}
/**
 * Reset the fit parameters 
 */
void ExtTrack::ResetFit(){

  // cout << "  ExtTrack::ResetFit() "<<endl;
    for(int i=0; i<5; i++) al[i]=-9999.;
    chi2=0.;
    nstep=0;

    for(int i=0; i<5; i++) {
        for(int j=0; j<5; j++) coval[i][j]=0.;
    }

    for(int i=0; i<nplanes; i++){
      xv[i] = 0.;
      yv[i] = 0.;
      zv[i] = 0.;
      axv[i] = 0.;
      ayv[i] = 0.;
    }

    int ngf = TrkParams::nGF;
    for(int i=0; i<ngf; i++){
      xGF[i] = 0.;
      yGF[i] = 0.;
    }


}
/**
 * Reset x/y measured coordinates (NB it does not reset z coordinates) 
 */
void ExtTrack::ResetXY(){

  for(int i=0; i<nplanes; i++)ResetXY(i);
      
    
}
/**
 * \brief Tracking method. It calls F77 mini routine.
 *
 * @see void TrkTrack::Fit(double, int&, int, int)
 */
void ExtTrack::Fit(double pfixed, int& fail, int iprint){

  TrkParams::Load(1); 
  if( !TrkParams::IsLoaded(1) ){
    cout << "void ExtTrack::Fit(double,int&,int,int) ---ERROR--- m.field not loaded "<<endl;
    return;
  }
 
  float al_ini[] = {0.,0.,0.,0.,0.};
  
  fail = 0;
  
  FillMiniStruct();   //fill F77 common
    

  // if fit variables have been reset, evaluate the initial guess
  if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.){
    if(iprint==2)cout << "calling guessext_()"<<endl;
    guessext_();
  }
  // --------------------- free momentum 
  if(pfixed==0.) {
    exttrack_.pfixed=0.; 
  }
  // --------------------- fixed momentum 
  if(pfixed!=0.) {
    al[4]=1./pfixed;     
    exttrack_.pfixed=pfixed;
  }
  
  //  store temporarily the initial guess
  for(int i=0; i<5; i++) al_ini[i]=exttrack_.al[i]; 

  if(iprint==2){
    cout<<endl<<" init. guess: ";
    for(int i=0; i<5; i++)cout << al_ini[i]<<" ";
    cout<<endl;
  }
  //  ------------------------------------------
  //  call mini routine
  //  ------------------------------------------
  int istep=0;
  int ifail=0;
  if(iprint==2)cout << "calling miniext_("<<istep<<","<<ifail<<", "<<iprint<<")"<<endl;
  miniext_(&istep,&ifail, &iprint);
  if(ifail!=0) {
    if(iprint)cout << "ERROR: ifail= " << ifail << endl;
    fail = 1;
  }
  if(chi2!=chi2){
    if(iprint)cout << "ERROR: chi2= " << chi2 << endl;  
    ResetFit();
    fail = 1;   
  }
  //  ------------------------------------------
  
  SetFromMiniStruct(); //copy content of F77 common into the ExtTrack object
  
  if(fail){
    if(iprint)cout << " >>>> fit failed "<<endl;
    for(int i=0; i<5; i++) al[i]=al_ini[i]; 
  }

  int ngf = TrkParams::nGF;
  float* zgf =  TrkParams::zGF;    
  Trajectory tj = Trajectory(ngf,zgf);
  tj.DoTrack(al,zini);
  for(int i=0; i<14; i++){
      xGF[i] = tj.x[i];
      yGF[i] = tj.y[i];
  }


};

/**
 * Method to fill minimization-routine common
 * NB (zini=23.5 is implicitelly set)
 */
void ExtTrack::FillMiniStruct(cMiniExtTrack& track){

  track.nplanes = nplanes;
  for(int i=0; i<nplanes; i++){ 
    
    track.xgood[i] = (double) XGood(i);
    track.ygood[i] = (double) YGood(i);
    
    track.xm[i] = (double) xm[i];
    track.ym[i] = (double) ym[i];
    track.zm[i] = (double) zm[i];
    
    
    track.xm_a[i] = (double) xma[i];
    track.xm_b[i] = (double) xmb[i];
    track.ym_a[i] = (double) yma[i];
    track.ym_b[i] = (double) ymb[i];
    track.zm_a[i] = (double) zma[i];
    track.zm_b[i] = (double) zmb[i];
    
    track.resx[i] = (double) resx[i];
    track.resy[i] = (double) resy[i];
  }
  
  for(int i=0; i<5; i++) track.al[i] = (double) al[i]; 
  
  track.zini = (double) zini;
  
    
}
/**
 * Method to set values from  minimization-routine common
 */
void ExtTrack::SetFromMiniStruct(cMiniExtTrack *track){

    for(int i=0; i<5; i++) {
      al[i] = (float) (track->al[i]);
      for(int j=0; j<5; j++) coval[i][j] = (float) (track->cov[i][j]);
    }
    chi2  = (float) (track->chi2);
    nstep = (float) (track->nstep);
    for(int i=0; i<nplanes; i++){ 
        xv[i]  = (float) (track->xv[i]);
        yv[i]  = (float) (track->yv[i]);
        zv[i]  = (float) (track->zv[i]);
        xm[i]  = (float) (track->xm[i]);
        ym[i]  = (float) (track->ym[i]);
        zm[i]  = (float) (track->zm[i]);
        xma[i]  = (float) (track->xm_a[i]);
        yma[i]  = (float) (track->ym_a[i]);
        zma[i]  = (float) (track->zm_a[i]);
        xmb[i]  = (float) (track->xm_b[i]);
        ymb[i]  = (float) (track->ym_b[i]);
        zmb[i]  = (float) (track->zm_b[i]);
        axv[i] = (float) (track->axv[i]);
        ayv[i] = (float) (track->ayv[i]);       
    }
    zini = (float) (track->zini);
}
/**
 * Method to retrieve ID (0,1,...) of x-cluster (if any) associated to this track.
 * If no cluster is associated, ID=-1.
 *
 */
int ExtTrack::GetClusterX_ID(int ip){
    if(ip<0 || ip>=nplanes)return -1;
    return ((int)fabs(xgood[ip]))%10000000-1;
};
/**
 * Method to retrieve ID (0-xxx) of y-cluster (if any) associated to this track.
 * 
 */
int ExtTrack::GetClusterY_ID(int ip){
    if(ip<0 || ip>=nplanes)return -1;
    return ((int)fabs(ygood[ip]))%10000000-1;
};

/**
 * Method to retrieve the dE/dx measured on a tracker view.
 * @param ip plane (0-5)
 * @param iv view (0=x 1=y)
 */
Float_t ExtTrack::GetDEDX(int ip, int iv){
    if(iv==0 && ip>=0 && ip<nplanes)return fabs(dedx_x[ip]);
    else if(iv==1 && ip>=0 && ip<nplanes)return fabs(dedx_y[ip]);
    else {
        cout << "TrkTrack::GetDEDX(int ip, int iv) -- wrong input parameters "<<ip<<iv<<endl;
        return 0.;
    }
}
/**
 * Method to evaluate the dE/dx measured on a tracker plane. 
 * The two measurements on x- and y-view are averaged.
 * @param ip plane (0-5)
 */
Float_t ExtTrack::GetDEDX(int ip){
    if( (Int_t)XGood(ip)+(Int_t)YGood(ip) == 0 ) return 0; 
    return (GetDEDX(ip,0)+GetDEDX(ip,1))/((Int_t)XGood(ip)+(Int_t)YGood(ip));
};

/**
 * Method to evaluate the dE/dx averaged over all planes. 
 */
Float_t ExtTrack::GetDEDX(){
    Float_t dedx=0; 
    for(Int_t ip=0; ip<nplanes; ip++)dedx+=GetDEDX(ip,0)*XGood(ip)+GetDEDX(ip,1)*YGood(ip); 
    dedx = dedx/(GetNX()+GetNY());
    return dedx;
};


/**
 * Method to evaluate the dE/dx averaged over all X views. 
 */
Float_t ExtTrack::GetDEDXX(bool cutSat){
    Float_t dedx=0;
    Int_t   np=0;
    for(Int_t ip=0; ip<nplanes; ip++){
        dedx += GetDEDX(ip,0)*XGood(ip)*(cutSat ? !IsSaturated(ip,0) : 1); 
        np   +=               XGood(ip)*(cutSat ? !IsSaturated(ip,0) : 1);
    }
    if(np>0)return dedx/np;
    return -1;
};
/**
 * Method to evaluate the dE/dx averaged over all Y views. 
 */
Float_t ExtTrack::GetDEDXY(bool cutSat){
    Float_t dedx=0;
    Int_t   np=0;
    for(Int_t ip=0; ip<nplanes; ip++){
        dedx += GetDEDX(ip,1)*YGood(ip)*(cutSat ? !IsSaturated(ip,1) : 1); 
        np   +=               YGood(ip)*(cutSat ? !IsSaturated(ip,1) : 1);
    }
    if(np>0)return dedx/np;
    return -1;
};

/**
 * Returns 1 if the cluster on a tracker view includes bad strips 
 * (at least one bad strip among the four strip used by p.f.a.)
 * @param ip plane (0-5)
 * @param iv view (0=x 1=y)
 */
Bool_t ExtTrack::IsBad(int ip,int iv){
    if(ip>=6&&ip<8)return false;
    if(iv==0 && ip>=0 && ip<6)return (xgood[ip]<0) ;
    else if(iv==1 && ip>=0 && ip<6)return (ygood[ip]<0) ;
    else {
        cout << "ExtTrack::IsBad(int ip, int iv) -- wrong input parameters "<<ip<<iv<<endl;
        return 0.;
    }
};

/**
 * Method to retrieve the ladder, if assigned
 */
int ExtTrack::GetLadder(int ip){
    if(XGood(ip))return (int)fabs(xgood[ip]/100000000)-1;
    if(YGood(ip))return (int)fabs(ygood[ip]/100000000)-1;
    return -1;
};
/**
 * Method to retrieve the sensor, if assigned 
 */
int ExtTrack::GetSensor(int ip){
    if(XGood(ip))return (int)((int)fabs(xgood[ip]/10000000)%10)-1;
    if(YGood(ip))return (int)((int)fabs(ygood[ip]/10000000)%10)-1;
    return -1;
};
/**
 * 
 */
Int_t ExtTrack::GetClusterX_Multiplicity(int ip){ 
  if(ip<0 || ip>=nplanes){
    cout << " ExtTrack::GetClusterX_Multiplicity("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
    return -1;
  }
  return (Int_t)(multmaxx[ip]/10000); 
}
/**
 * 
 */

Int_t ExtTrack::GetClusterY_Multiplicity(int ip){ 
  if(ip<0 || ip>=nplanes){
    cout << " ExtTrack::GetClusterY_Multiplicity("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
    return -1;
  }
  return (Int_t)(multmaxy[ip]/10000); 
}
/**
 * 
 */

Int_t ExtTrack::GetClusterX_MaxStrip(int ip){ 
  if(ip<0 || ip>=nplanes){
    cout << " ExtTrack::GetClusterX_MaxStrip("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
    return -1;
  }
  return (Int_t)(multmaxx[ip]%10000); 
}
/**
 * 
 */

Int_t ExtTrack::GetClusterY_MaxStrip(int ip){ 
  if(ip<0 || ip>=nplanes){
    cout << " ExtTrack::GetClusterY_MaxStrip("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
    return -1;
  }
  return (Int_t)(multmaxy[ip]%10000); 
}
Float_t ExtTrack::GetRigidity(){ 
        Float_t rig=0;
        if(chi2>0)rig=1./al[4];
        if(rig<0) rig=-rig;
        return rig;
};
//
Float_t ExtTrack::GetDeflection(){
        Float_t def=0;
        if(chi2>0)def=al[4];
        return def;
};


//
// all that follows: EM porting from TrkLevel2
//
Bool_t ExtTrack::IsInsideAcceptance(float toll){
    int ngf = TrkParams::nGF;
    float* zgf =  TrkParams::zGF;
    
    Trajectory tj = Trajectory(ngf,zgf);
    tj.DoTrack(al,zini);

    for(int i=0; i<ngf; i++){
        //
//      cout << endl << TrkParams::GF_element[i];
//      if( 
//          TrkParams::GF_element[i].CompareTo("S11") &&
//          TrkParams::GF_element[i].CompareTo("S12") &&
//          TrkParams::GF_element[i].CompareTo("S21") &&
//          TrkParams::GF_element[i].CompareTo("S22") &&
//          TrkParams::GF_element[i].CompareTo("T1")  &&
//          TrkParams::GF_element[i].CompareTo("CUF") &&
//          TrkParams::GF_element[i].CompareTo("T2")  &&
//          TrkParams::GF_element[i].CompareTo("T3")  &&
//          TrkParams::GF_element[i].CompareTo("T4")  &&
//          TrkParams::GF_element[i].CompareTo("T5")  &&
//          TrkParams::GF_element[i].CompareTo("CLF") &&
//          TrkParams::GF_element[i].CompareTo("T6")  &&
//          TrkParams::GF_element[i].CompareTo("S31") &&
//          TrkParams::GF_element[i].CompareTo("S32") &&
//          true)continue;
        // apply condition only within the cavity
//      cout << " -- "<<xGF[i]<<" "<<yGF[i];
        if( 
//          xGF[i] <= TrkParams::xGF_min[i] + toll ||
//          xGF[i] >= TrkParams::xGF_max[i] - toll ||
//          yGF[i] <= TrkParams::yGF_min[i] + toll ||
//          yGF[i] >= TrkParams::yGF_max[i] - toll ||
            tj.x[i] <= TrkParams::xGF_min[i] + toll ||
            tj.x[i] >= TrkParams::xGF_max[i] - toll ||
            tj.y[i] <= TrkParams::yGF_min[i] + toll ||
            tj.y[i] >= TrkParams::yGF_max[i] - toll ||
            false){
            
            return false;
        }
    }
    return true;
}

/**
 * Returns the reduced chi-square of track x-projection
 */
Float_t  ExtTrack::GetChi2X(){
    float chiq=0;
    for(int ip=0; ip<nplanes; ip++)if(XGood(ip))chiq+= pow((xv[ip]-xm[ip])/resx[ip],2.); 
    if(GetNX()>3)chiq=chiq/(GetNX()-3);
    else chiq=0;
    //    if(chiq==0)cout << " Float_t  ExtTrack::GetChi2X() -- WARNING -- value not defined "<<chiq<<endl;
    return chiq;
}
/**
 * Returns the reduced chi-square of track y-projection
 */
Float_t  ExtTrack::GetChi2Y(){
    float chiq=0;
    for(int ip=0; ip<nplanes; ip++)if(YGood(ip))chiq+= pow((yv[ip]-ym[ip])/resy[ip],2.); 
    if(GetNY()>2)chiq=chiq/(GetNY()-2);
    else chiq=0;
    //    if(chiq==0)cout << " Float_t  ExtTrack::GetChi2Y() -- WARNING -- value not defined "<<chiq<<endl;
    return chiq;
}

/**
 * Returns the track "lever-arm" on the x view, defined as the distance (in planes) between 
 * the upper and lower x measurements (the maximum value of lever-arm is 6). 
 */
Int_t ExtTrack::GetLeverArmX(){
    int first_plane = -1;
    int last_plane  = -1;
    for(Int_t ip=0; ip<nplanes; ip++){
        if( XGood(ip) && first_plane == -1 )first_plane = ip;
        if( XGood(ip) && first_plane != -1 )last_plane = ip;
    }
    if( first_plane == -1 || last_plane == -1){
        cout<< "Int_t ExtTrack::GetLeverArmX() -- XGood(ip) always false ??? "<<endl; 
        return 0;
    }
    return (last_plane-first_plane+1);
}
/**
 * Returns the track "lever-arm" on the y view, defined as the distance (in planes) between 
 * the upper and lower y measurements (the maximum value of lever-arm is 6). 
 */
Int_t ExtTrack::GetLeverArmY(){
    int first_plane = -1;
    int last_plane  = -1;
    for(Int_t ip=0; ip<nplanes; ip++){
        if( YGood(ip) && first_plane == -1 )first_plane = ip;
        if( YGood(ip) && first_plane != -1 )last_plane = ip;
    }
    if( first_plane == -1 || last_plane == -1){
        cout<< "Int_t ExtTrack::GetLeverArmY() -- YGood(ip) always false ??? "<<endl; 
        return 0;
    }
    return (last_plane-first_plane+1);
}
/**
 * Returns the track "lever-arm" on the x+y view, defined as the distance (in planes) between 
 * the upper and lower x,y (couple) measurements (the maximum value of lever-arm is 6). 
 */
Int_t ExtTrack::GetLeverArmXY(){
    int first_plane = -1;
    int last_plane  = -1;
    for(Int_t ip=0; ip<nplanes; ip++){
        if( XGood(ip) && YGood(ip) && first_plane == -1 )first_plane = ip;
        if( XGood(ip) && YGood(ip) && first_plane != -1 )last_plane = ip;
    }
    if( first_plane == -1 || last_plane == -1){
//      cout<< "Int_t ExtTrack::GetLeverArmXY() -- XGood(ip)*YGood(ip) always false ??? "<<endl; 
        return 0;
    }
    return (last_plane-first_plane+1);
}


/**
 * Returns 1 if the signal on a tracker view is saturated.
 * @param ip plane (0-5)
 * @param iv view (0=x 1=y)
 */
Bool_t ExtTrack::IsSaturated(int ip,int iv){
    if(ip>=6 && ip<8)return false;
    if(iv==0 && ip>=0 && ip<6)return (dedx_x[ip]<0) ;
    else if(iv==1 && ip>=0 && ip<6)return (dedx_y[ip]<0) ;
    else {
        cout << "ExtTrack::IsSaturated(int ip, int iv) -- wrong input parameters "<<ip<<iv<<endl;
        return 0.;
    }
};
/**
 * Returns 1 if either the x or the y signal on a tracker plane is saturated.
 * @param ip plane (0-5)
 */
Bool_t ExtTrack::IsSaturated(int ip){
    return (IsSaturated(ip,0)||IsSaturated(ip,1));
};
/**
 * Returns 1 if there is at least a saturated signal along the track.
 */
Bool_t ExtTrack::IsSaturated(){
    for(int ip=0; ip<nplanes; ip++)for(int iv=0; iv<2; iv++)if(IsSaturated(ip,iv))return true;
    return false;
}


ClassImp(ExtTrack);