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

CbmL1MCTrack::CbmL1MCTrack(double mass_, double q_, TVector3 vr, TLorentzVector vp, int _ID, int _mother_ID, int _pdg):
  mass(mass_),q(q_),p(0),x(0),y(0),z(0),px(0),py(0),pz(0),ID(_ID), mother_ID(_mother_ID), pdg(_pdg),Points(),StsHits(),
  nMCContStations(0),nHitContStations(0),maxNStaMC(0),maxNSensorMC(0),maxNStaHits(0),nStations(0),nMCStations(0),isReconstructable(0),isAdditional(),
     rTracks(),tTracks()
{
  x = vr.X();
  y = vr.Y();
  z = vr.Z();
  px = vp.Px();
  py = vp.Py();
  pz = vp.Pz();
  p = sqrt( fabs(px*px + py*py + pz*pz ));
};

// CbmL1MCTrack::CbmL1MCTrack(TmpMCPoints &mcPoint, TVector3 vr, TLorentzVector vp, int ID, int mother_ID):
//      ID(_ID), mother_ID(_mother_ID)
// {
//   mass = mcPoint->mass;
//   q = mcPoint->q;
//   pdg  = mcPoint->pdg;
//   
//   x = vr.X();
//   y = vr.Y();
//   z = vr.Z();
//   px = vp.Px();
//   py = vp.Py();
//   pz = vp.Pz();
//   p = sqrt( fabs(px*px + py*py + pz*pz ));
// };

void CbmL1MCTrack::Init()
{
  CbmL1* L1 = CbmL1::Instance();
     // get stsHits
  StsHits.clear();
  for (unsigned int iP = 0; iP < Points.size(); iP++){
    CbmL1MCPoint* point = &(L1->vMCPoints[Points[iP]]);
    for (unsigned int iH = 0; iH < point->hitIds.size(); iH++){
      const int iih = point->hitIds[iH];
      if( std::find(StsHits.begin(), StsHits.end(), iih) == StsHits.end() )
        StsHits.push_back(iih);
    }
  }
  
  CalculateMCCont();
  CalculateHitCont();
  CalculateMaxNStaMC();
  CalculateMaxNStaHits();
  CalculateIsReconstructable();
} // void CbmL1MCTrack::Init()

void CbmL1MCTrack::CalculateMCCont()
{
  CbmL1* L1 = CbmL1::Instance();
  
  int nPoints = Points.size();
  nMCContStations = 0;
  int istaold = -1, ncont=0;
  for( int ih=0; ih<nPoints; ih++ ){
    CbmL1MCPoint &h = L1->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()
{
  CbmL1* L1 = CbmL1::Instance();
  L1Algo* algo = L1->algo;

  int nhits = StsHits.size();
  nHitContStations = 0;
  int istaold = -1, ncont=0;
  for( int ih=0; ih<nhits; ih++ ){
    int jh = StsHits[ih];
    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()
{
  CbmL1* L1 = CbmL1::Instance();

  maxNStaHits = 0;
  nStations = 0;
  int lastSta = -1;
  int cur_maxNStaHits = 0;
  for(unsigned int iH = 0; iH < StsHits.size(); iH++){
    CbmL1HitStore& sh = L1->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::CalculateMaxNStaMC()
{
  CbmL1* L1 = CbmL1::Instance();

  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++){
    CbmL1MCPoint& mcP = L1->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()

void CbmL1MCTrack::CalculateIsReconstructable()
{
  CbmL1* L1 = CbmL1::Instance();
  
  bool f = 1;
  
    // reject very slow tracks from analysis
  f &= (p > CbmL1Constants::MinRecoMom);
    // detected at least in 4 stations
    //   f &= (nMCContStations >= 4);
  
    // maximul 4 layers for a station.
    //   f &= (maxNStaHits <= 4);
  f &= (maxNStaMC <= 4);
    //   f &= (maxNSensorMC <= 1);
  
  if (L1->fPerformance == 3) isReconstructable = f & (nMCContStations  >= CbmL1Constants::MinNStations); // L1-MC
  if (L1->fPerformance == 2) isReconstructable = f & (nStations        >= CbmL1Constants::MinNStations); // QA definition
  if (L1->fPerformance == 1) isReconstructable = f & (nHitContStations >= CbmL1Constants::MinNStations); // L1 definition

  isAdditional = f &
    (nHitContStations == nStations) & (nMCContStations == nStations) & (nMCStations == nStations) &
    (nHitContStations >= 3) &
    (L1->vMCPoints[Points[0]].iStation == 0);
  isAdditional &= !isReconstructable;
}; // bool CbmL1MCTrack::IsReconstructable()