/*
 *====================================================================
 *
 *  CBM Level 1 Reconstruction 
 *  
 *  Authors: I.Kisel,  S.Gorbunov
 *
 *  e-mail : ikisel@kip.uni-heidelberg.de 
 *
 *====================================================================
 *
 *  L1 Monte Carlo information
 *
 *====================================================================
 */
#include "../TaskManager.h"
#include "../L1Algo/L1Algo.h"
#include "CbmL1MCTrack.h"
#include "Performance.h"
#include "PerformanceConstants.h"

void CbmL1MCTrack::CalculateMCCont(const vector<CbmL1MCPoint>& vMCPoints)
{
  
  int nPoints = Points.size();
  nMCContStations = 0;
  int istaold = -1, ncont=0;
  for( int ih=0; ih<nPoints; ih++ ){
    const CbmL1MCPoint &h = vMCPoints[Points[ih]];
    int ista = h.iStation;
    if (ista - istaold == 1) ncont++;
    else if(ista - istaold > 1){
      if( nMCContStations < ncont ) nMCContStations = ncont;
      ncont = 1;
    }
    if (ista <= istaold ) continue; // backward direction
    istaold = ista;
  }
  if( nMCContStations < ncont ) nMCContStations = ncont;
}; // void CbmL1MCTrack::CalculateMCCont()

void CbmL1MCTrack::CalculateHitCont(const vector<CbmL1MCPoint>& vMCPoints)
{ 
  
  TaskManager* fTracker = TaskManager::Instance();
  const L1Algo& algo = fTracker->GetAlgo();

  int nhits = StsHits.size();
  nHitContStations = 0;
  int istaold = -1, ncont=0;
  {
    for( int ih=0; ih<nhits; ih++ ){
      int jh = StsHits[ih];

      const L1StsHit &h = (*(algo.vStsHits))[jh];
      int ista = (*(algo.vSFlag))[h.f]/4;
      if (ista - istaold == 1) ncont++;
      else if(ista - istaold > 1){
        if( nHitContStations < ncont ) nHitContStations = ncont;
        ncont = 1;
      }

      if ( !( ista >= istaold ) ) { // tracks going in backward direction are not reconstructable
        nHitContStations = 0;
        return;
      }
      if (ista == istaold ) continue; // backward direction
      istaold = ista;
    }
  }
  if( nHitContStations<ncont ) nHitContStations=ncont;
}; // void CbmL1MCTrack::CalculateHitCont()

void CbmL1MCTrack::CalculateMaxNStaHits(const vector< CbmL1HitStore >& vHitStore)
{

  maxNStaHits = 0;
  nStations = 0;
  int lastSta = -1;
  int cur_maxNStaHits = 0;
  for(unsigned int iH = 0; iH < StsHits.size(); iH++){
    
    
    const CbmL1HitStore& sh = vHitStore[StsHits[iH]];
    
    if (sh.iStation == lastSta){ 
      cur_maxNStaHits++;
    }
    else{ // new station
      if ( !(sh.iStation > lastSta) ) {// tracks going in backward direction are not reconstructable
        maxNStaHits = 0;
        return;
      }
      if (cur_maxNStaHits > maxNStaHits) maxNStaHits = cur_maxNStaHits;
      cur_maxNStaHits = 1;
      lastSta = sh.iStation;
      nStations++;
    }
  };
  if (cur_maxNStaHits > maxNStaHits) maxNStaHits = cur_maxNStaHits;
//   cout << pdg << " " << p << " " << StsHits.size() << " > " << maxNStaHits << endl;
}; // void CbmL1MCTrack::CalculateHitCont()

void CbmL1MCTrack::CalculateIsReconstructable( vector<CbmL1MCPoint>& vMCPoints,  vector< L1StsHit >& vStsHits,  vector< unsigned char >& vSFlag,  vector< CbmL1HitStore >& vHitStore)
{

  CalculateMaxNStaMC(vMCPoints);

  CalculateMCCont(vMCPoints);

  //CalculateMaxNStaHits(vHitStore);
   
  //CalculateHitCont(vMCPoints, vStsHits, vSFlag);
    
  
  bool f = 1;
  
    // reject very slow tracks from analysis
  f &= (p > PerformanceConstants::MinRecoMom);
    // detected at least in 4 stations
//   f &= (nMCContStations >= 4);
//   f &= (nStations        >= L1->MinNStations);     // QA definition
    // f &= (nHitContStations >= Performance::MinNStations); // L1 definition
    // maximul 4 layers for a station.
//   f &= (maxNStaHits <= 4);
  f &= (maxNStaMC <= 4);
//   f &= (maxNSensorMC <= 1);
  isReconstructable = f & (nMCContStations >= PerformanceConstants::MinNStations);
  if(Points.size() > 0)
  {
    isAdditional = f &
      (nHitContStations == nStations) & (nMCContStations == nStations) & (nMCStations == nStations) &
      (nHitContStations >= 3) &
      (vMCPoints[Points[0]].iStation == 0);
    isAdditional &= !isReconstructable;
  }
  else
    isAdditional = 0;
  
}; // bool CbmL1MCTrack::IsReconstructable()

void CbmL1MCTrack::CalculateMaxNStaMC(const vector<CbmL1MCPoint>& vMCPoints)
{
  maxNStaMC = 0;
  maxNSensorMC = 0;
  nMCStations = 0;
  int lastSta = -1;
  float lastz = -1;
  int cur_maxNStaMC = 0 , cur_maxNSensorMC = 0;
  for(unsigned int iH = 0; iH < Points.size(); iH++){
    const CbmL1MCPoint& mcP = vMCPoints[Points[iH]];
    if (mcP.iStation == lastSta)
      cur_maxNStaMC++;
    else{ // new station
      if (cur_maxNStaMC > maxNStaMC) maxNStaMC = cur_maxNStaMC;
      cur_maxNStaMC = 1;
      lastSta = mcP.iStation;
      nMCStations++;
    }
    
    if (mcP.z == lastz)   // TODO: works incorrect because of different z
      cur_maxNSensorMC++;
    else{ // new z
      if (cur_maxNSensorMC > maxNSensorMC) maxNSensorMC = cur_maxNSensorMC;
      cur_maxNSensorMC = 1;
      lastz = mcP.z;
    }
  };
  if (cur_maxNStaMC > maxNStaMC) maxNStaMC = cur_maxNStaMC;
  if (cur_maxNSensorMC > maxNSensorMC) maxNSensorMC = cur_maxNSensorMC;
//   cout << pdg << " " << p << " " << Points.size() << " > " << maxNStaMC << " " << maxNSensorMC << endl;
}; // void CbmL1MCTrack::CalculateMaxNStaMC()