CbmL1TrackFitter.cxx

Go to the documentation of this file.

/*
 *=====================================================
 *
 *  CBM Level 1 Reconstruction 
 *  
 *  Authors: I.Kisel,  S.Gorbunov
 *
 *  e-mail : ikisel@kip.uni-heidelberg.de 
 *
 *=====================================================
 *
 *  Fit reconstructed tracks and find primary vertex
 *
 */
 
#include "CbmL1.h"
 
#include "L1Algo/L1Algo.h"
#include "L1Algo/L1TrackPar.h"
#include "L1Algo/L1Extrapolation.h"
#include "L1Algo/L1Filtration.h"
#include "L1Algo/L1AddMaterial.h"
#include "L1Algo/L1StsHit.h"
 
#include "CbmKF.h"
#include "CbmKFMath.h"
#include "CbmKFPrimaryVertexFinder.h"
 
#include "TStopwatch.h"
 
void CbmL1::TrackFitter( vector<CbmL1Track> &Tracks, CbmL1Vtx *V )
{
  TStopwatch timer;
  timer.Start();
  /*
  //int nt=0;
  for( vector<CbmL1Track>::iterator i = Tracks.begin(); i!=Tracks.end(); ++i)
    {
      CbmL1Track &T = *i;
      
      int stmin = 1000, stmax = -1000;
      for ( vector<L1StsHit*>::iterator j = T.StsHits.begin(); j != T.StsHits.end(); ++j)
  {      
    if ((*j)->iStation < stmin ) stmin = (*j)->iStation;
    if ((*j)->iStation > stmax ) stmax = (*j)->iStation;
  }
      
      T.nStations = stmax - stmin + 1;
      T.is_long = ( T.StsHits.size() >= 4 ) && ( T.nStations >=4 );
 
      //T.is_long = T.is_long &&( vMCPoints[T.MapsHits.back()->MC_Point].mother_ID==-1 );
      //if ( !T.is_long ) continue;
      
      static L1FieldRegion fld[20];
      {
  vector<L1StsHit*>::iterator ih2 = T.StsHits.begin(), ih1, ih0;
  ih1 = ih2;
  ih1++;
  ih0 = ih1; 
  ih0++;
  int ist0 = (*ih0)->iStation;
  int ist1 = (*ih1)->iStation;
  int ist2 = (*ih2)->iStation;
  L1Station *st0 = &(algo->vStations[ist0]);
  L1Station *st1 = &(algo->vStations[ist1]);
  L1Station *st2 = &(algo->vStations[ist2]);
  L1FieldValue B0, B1, B2;
  fvec Z0= st0->z, Z1= st1->z, Z2= st2->z ; 
  st0->fieldSlice.GetFieldValue( (*ih0)->x, (*ih0)->y, B0 );
  st1->fieldSlice.GetFieldValue( (*ih1)->x, (*ih1)->y, B1 );
  st2->fieldSlice.GetFieldValue( (*ih2)->x, (*ih2)->y, B2 );
  fld[ist0].Set( B0, Z0, B1, Z1, B2, Z2 );
  fld[ist1] = fld[ist0];
  fld[ist2] = fld[ist0];
  if( ih0!=T.StsHits.end() ){
    while(1){ 
      ih2 = ih1; B2 = B1; Z2 = Z1; st2 = st1;
      ih1 = ih0; B1 = B0; Z1 = Z0; st1 = st0;
      ih0++;
      if( ih0==T.StsHits.end() ) break;
      ist0 = (*ih0)->iStation;
      st0 = &(algo->vStations[ist0]);
      Z0 = st0->z;
      st0->fieldSlice.GetFieldValue( (*ih0)->x, (*ih0)->y, B0 );
      fld[ist0].Set( B0, Z0, B1, Z1, B2, Z2 );
    }
  }
      }
      fvec vINF = 100.;
      L1TrackPar tp;
      for( int i=0; i<5; i++ ) *((&tp.x)+i) = T.T[i];
      {
  fvec qp0  = tp.qp;
  tp.NDF  = 2;
  tp.chi2 = 0;
  vector<L1StsHit*>::iterator ih = T.StsHits.begin();
  int ist = (*ih)->iStation;
  L1Station *st = &(algo->vStations[ist]);  
  tp.x  = (*ih)->x;
  tp.y  = (*ih)->y;
  tp.z  = st->z;
  tp.C00 = st->XYInfo.C00;
  tp.C10 = st->XYInfo.C10;
  tp.C11 = st->XYInfo.C11;
  for( int i=3; i<15; i++ ) *((&tp.C00)+i) = 0;
  tp.C22 = tp.C33 = tp.C44 = vINF;
  for( ; ih!= T.StsHits.end(); ih++ ){
    ist = (*ih)->iStation;
    st = &(algo->vStations[ist]);
    L1Extrapolate( tp, st->z, qp0, fld[ist] );
    fvec u_f = (*ih)->u_front, u_b = (*ih)->u_back;
    L1Filter( tp, st->frontInfo, u_f  );
    L1Filter( tp, st->backInfo, u_b  );
    L1AddMaterial( tp, st->materialInfo, qp0 );
  }
      } 
      for( int i=0; i<6;  i++ ) T.TLast[i] = (*((&tp.x)+i))[0] ;
      for( int i=0; i<15; i++ ) T.CLast[i] = (*((&tp.C00)+i))[0];
 
      {
  fvec qp0  = tp.qp;
  tp.NDF  = 2;
  tp.chi2 = 0;
  vector<L1StsHit*>::reverse_iterator ih = T.StsHits.rbegin();
  int ist = (*ih)->iStation;
  L1Station *st = &(algo->vStations[ist]);  
  tp.x  = (*ih)->x;
  tp.y  = (*ih)->y;
  tp.z  = st->z;
  tp.C00 = st->XYInfo.C00;
  tp.C10 = st->XYInfo.C10;
  tp.C11 = st->XYInfo.C11;
  for( int i=3; i<15; i++ ) *((&tp.C00)+i) = 0;
  tp.C22 = tp.C33 = tp.C44 = vINF;
  for( ; ih!= T.StsHits.rend(); ih++ ){
    ist = (*ih)->iStation;
    st = &(algo->vStations[ist]);
    L1Extrapolate( tp, st->z, qp0, fld[ist] );
    fvec u_f = (*ih)->u_front, u_b = (*ih)->u_back;
    L1Filter( tp, st->frontInfo, u_f  );
    L1Filter( tp, st->backInfo, u_b  );
    L1AddMaterial( tp, st->materialInfo, qp0 );
  }
      } 
      for( int i=0; i<6;  i++ ) T.T[i] = (*((&tp.x)+i))[0] ;
      for( int i=0; i<15; i++ ) T.C[i] = (*((&tp.C00)+i))[0];
      T.chi2 = tp.chi2[0];
      T.NDF = tp.NDF[0]-5;    
    }
  timer.Stop();
 
  stat_vtx_ntracks = 0; 
 
  if(V){ // primary vertex fit  
    
    vector<CbmKFTrackInterface*> vTracks;
    
   { for( vector<CbmL1Track>::iterator i = Tracks.begin(); i!=Tracks.end(); ++i)
      {
  if( ( i->StsHits.size() < 5 ) || ( i->nStations <5 ) ) continue;
  vTracks.push_back( &*i );    
  stat_vtx_ntracks ++;
      }
    }  
   cout<<"Fit vertex"<<endl;
    if ( vTracks.size()>=3 )
      {
  CbmKFPrimaryVertexFinder PVFinder;
  PVFinder.SetTracks( vTracks );
  PVFinder.Fit(*V);
      }
  }// prim. vtx 
  */
  //c_time.Stop();    
  //stat_fit_time += double(c_time.Time());
  //timer.Stop();
  stat_fit_time += timer.CpuTime();
  stat_fit_time_ntracks += Tracks.size();
}