/*
 *====================================================================
 *
 *  CBM Level 1 Reconstruction 
 *  
 *  Authors: I.Kisel,  S.Gorbunov
 *
 *  e-mail : ikisel@kip.uni-heidelberg.de 
 *
 *====================================================================
 *
 *  Read Event information to L1
 *
 *====================================================================
 */

#include "CbmL1.h"

struct TmpHitSort{
  int UID;
  double z;
  CbmL1StsHit *StsHit;
  CbmL1MCPoint *MCPoint;

  static bool compareIDz( const TmpHitSort &a, const TmpHitSort &b )
  {
    return ( a.UID < b.UID ) || ( ( a.UID == b.UID ) && (a.z < b.z) );
  }
};

void CbmL1::ReadEvent()
{
  vMCPoints.clear();
  vMCTracks.clear();
  vRTracks.clear();

  vStsHits.clear();
  vStsStrips.clear();
  vTrdHits.clear();

  // produce Sts hits from space points

  for(int i = 0; i < NStation; i++){
    StsHitsStartIndex[i] = -1;
    StsHitsStopIndex[i]  = -2;
  }
  
  // get STS hits

  int nhits=0;
  if( listStsHits ){ 
    map<int, int > Event2PileUpmap;
    Event2PileUpmap.insert(pair<int, int>(0, 0));
    int nPileUp = 0;
    int nHits = listStsHits->GetEntries();
    for (int j=0; j <  nHits; j++ )
      {
	CbmStsHit *sh = (CbmStsHit*) listStsHits->At(j);	
	CbmL1StsHit h;
	ReadStsHit( &h, j );	
	int iMC = sh->GetRefIndex();
	if( listStsPts && iMC>=0 ){
	  CbmL1MCPoint MC;
	  if( ! ReadMCPoint( &MC, iMC, 0 ) ){
	    vMCPoints.push_back( MC );
	    h.MC_Point = vMCPoints.size()-1;	    
	  }
	}
	//if(  h.MC_Point >=0 ) 
	{
	  vStsHits.push_back(h);
	  nhits++;
	}
      }
  }
  cout<<" L1 Sts hits :"<< nhits<<endl;

  sort( vStsHits.begin(), vStsHits.end(), CbmL1StsHit::compareStation );      
   
  int index = 0;
  for ( vector<CbmL1StsHit>::iterator im = vStsHits.begin(); im != vStsHits.end(); ++im){
    CbmL1StsHit* pm = &(*im);
    int sta = pm->iStation;
    if (StsHitsStartIndex[sta] < 0) StsHitsStartIndex[sta] = index;    
    StsHitsStopIndex[sta] = index;	    
    index++;
  }      

  // create pointers strip -> hit

  cout<<"Making strips"<<endl;

  for ( vector<CbmL1StsHit>::iterator ih = vStsHits.begin(); ih != vStsHits.end(); ++ih){
    if( !ih->isStrip ) continue;
    int sf=-1, sb=-1, ii=0;
    for ( vector<CbmL1StsStrip>::iterator is = vStsStrips.begin(); is != vStsStrips.end(); ++is, ++ii)
      {
	if( is->iStation!=ih->iStation ||  is->iSector!=ih->iSector ) continue;
	if( is->iStrip==ih->iStripF ) sf = ii;
	if( is->iStrip==ih->iStripB ) sb = ii;
      }
    if( sf<0 && ih->iStripF>=0 )
      {
	CbmL1StsStrip tmp;
	tmp.iStation = ih->iStation;
	tmp.iSector = ih->iSector;
	tmp.iStrip = ih->iStripF;
	tmp.isFront = 1;
	vStsStrips.push_back(tmp);
	sf = vStsStrips.size()-1;
	ih->indStripF = sf;
	vStsStrips[sf].vHits.push_back(&*ih);
     }
    if( sb<0 && ih->iStripB>=0)
      {
	CbmL1StsStrip tmp;
	tmp.iStation = ih->iStation;
	tmp.iSector = ih->iSector;
	tmp.iStrip = ih->iStripB;
	tmp.isFront = 0;
	vStsStrips.push_back(tmp);
	sb = vStsStrips.size()-1;
	ih->indStripB = sb;
	vStsStrips[sb].vHits.push_back(&*ih);
      }
  }


  // produce Trd hits from space points

  if (0&& listTrdPts )
    {
      double const Sigma = 500.E-4;
      /*
	for (int j=0; j <  listTrdPts->GetEntries(); j++ )
	{
	CbmTrdPoint *pt = (CbmTrdPoint*) listTrdPts->At(j);
	
	  if ( !pt ) continue;      
	  
	  CbmL1TrdHit h;
	  
	  if ( pt->GetDetectorID() < FirstTrdStationID )
	    {
	      FirstTrdStationID = pt->GetDetectorID();
	    }
	  
	  h.iStation = ( pt->GetDetectorID() - FirstTrdStationID);

	  CbmL1MCPoint MC;
	  MC.ID = pt->GetTrackID();	  
	  CbmMCTrack *MCTrack = (CbmMCTrack*) listMCTracks->At( MC.ID );
	  if ( !MCTrack ) continue;	  
	  MC.pdg  = MCTrack->GetPdgCode();	  
	  MC.mother_ID = MCTrack->GetMotherID(); 
	  if ( ( MC.pdg == 10010020) || // what are these particles?????
	       ( MC.pdg == 10010030) ||
	       ( MC.pdg == 50000050) ||
	       ( MC.pdg == 50010051) ||
	       ( MC.pdg == 10020040)     ) continue;	  
	  TVector3 xyz,P;
	  pt->Position(xyz);
	  pt->Momentum(P);
	  MC.x  = xyz.X();
	  MC.y  = xyz.Y();
	  MC.z  = xyz.Z();
	  MC.px = P.X();
	  MC.py = P.Y();
	  MC.pz = P.Z();
	  
	  h.z()      = MC.z;
	  h.sigma2() = Sigma*Sigma;
	  
	  h.phi()     = 0;
	  if ( 406.<=h.z() && h.z()<412. || 
	       606.<=h.z() && h.z()<612. || 
	       806.<=h.z() && h.z()<812.    ) h.phi() = 90./180.*3.1415927; 
	  //if(h.iStation%2) h.phi() = 90./180.*3.1415927;
	  h.phi_s()   = sin(h.phi());
	  h.phi_c()   = cos(h.phi());
	  h.phi_ss()  = h.phi_s() * h.phi_s();
	  h.phi_cc()  = h.phi_c() * h.phi_c();
	  h.phi_2sc() = 2.* h.phi_s() * h.phi_c();
	  
	  h.u() = h.phi_c() * MC.x  +  h.phi_s() * MC.y + gRandom->Gaus(0,Sigma);
	  
	  vMCPoints.push_back( MC );
	  h.MC_Point = vMCPoints.size()-1;
	  vTrdHits.push_back(h);	  
	}

      sort( vTrdHits.begin(), vTrdHits.end(), CbmL1TrdHit::compareStation );
      */
    }

  cout<<" sorting MC points.."<<endl;
  vector<TmpHitSort> vTmp;

  {for ( vector<CbmL1StsHit>::iterator i = vStsHits.begin(); i != vStsHits.end(); ++i)
    {
      int iMC = i->MC_Point;
      if( iMC < 0 ) continue;
      CbmL1MCPoint &MC = vMCPoints[iMC];
      TmpHitSort tmp;
      tmp.UID = MC.UID;
      tmp.z = MC.z;
      tmp.StsHit = &(*i);
      tmp.MCPoint = &MC;
      vTmp.push_back(tmp);
    }
  }
  sort( vTmp.begin(), vTmp.end(), TmpHitSort::compareIDz );
  cout<<"Fill MC tracks.."<<endl;

  // fill MC tracks

  PrimVtx.MC_ID = 999;

  {
    CbmL1Vtx Vtxcurr;
    int nvtracks=0, nvtrackscurr=0;
    int UID = -10000;
    for ( vector<TmpHitSort>::iterator i= vTmp.begin(); i!=vTmp.end(); ++i)
      {
	if( vMCTracks.empty() || i->UID!= UID )
	  {
	    //cout<<endl;
	    UID = i->UID;
	    CbmL1MCTrack T;
	    T.ID = UID/100;
	    T.iPileUp = UID%100;
	    T.UID = UID;
	    if( !T.iPileUp && listMCTracks)
	      {
		T.q = i->MCPoint->q;
		T.mass = i->MCPoint->mass;
		T.pdg  = i->MCPoint->pdg;
		CbmMCTrack *MCTrack = (CbmMCTrack*) listMCTracks->At( T.ID );
		T.mother_ID = MCTrack->GetMotherID();
		TVector3 vr = MCTrack->GetStartVertex();
		TLorentzVector vp = MCTrack->Get4Momentum();
		T.px = vp.Px();
		T.py = vp.Py();
		T.pz = vp.Pz();
		T.x = vr.X();
		T.y = vr.Y();
		T.z = vr.Z();
		T.p = sqrt( T.px*T.px + T.py*T.py + T.pz*T.pz );
		T.T[0] = T.x;
		T.T[1] = T.y;
		T.T[2] = T.px/T.pz;
		T.T[3] = T.py/T.pz;
		T.T[4] = T.q/T.p;
		T.T[5] = T.z;
	      }
	    vMCTracks.push_back( T );	      
	    if ( T.mother_ID ==-1 && !T.iPileUp ){ // vertex track

	      if (  PrimVtx.MC_ID == 999 || fabs(T.z-Vtxcurr.MC_z)>1.e-7 ){// new vertex
		if( nvtrackscurr > nvtracks ){
		  PrimVtx = Vtxcurr;
		  nvtracks = nvtrackscurr;
		}
		Vtxcurr.MC_x  = T.x;
		Vtxcurr.MC_y  = T.y;
		Vtxcurr.MC_z  = T.z;
		Vtxcurr.MC_ID = T.mother_ID;
		nvtrackscurr = 1;
	      }else nvtrackscurr++;
	    }
	  }
	vMCTracks.back().Points.push_back(i->MCPoint);
	vMCTracks.back().StsHits.push_back(i->StsHit);
 	//cout<<i->UID/100<<" ";
     }
    if( nvtrackscurr > nvtracks ) PrimVtx = Vtxcurr;

    //cout<<endl;
  } // fill MC tracks
  
  if (  PrimVtx.MC_ID == 999 ) cout<<"No primary vertex !!!"<<endl;

}