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

#include "CbmL1.h"
#include "CbmKF.h"
#include "CbmRunAna.h"

ClassImp(CbmL1)

CbmL1 *CbmL1::fInstance = 0;


CbmL1::CbmL1():ref(0)
{
  if( !fInstance ) fInstance = this;
}

CbmL1::~CbmL1()
{
  CbmRootManager::Instance()->Write();
  fInstance = 0;
}


CbmL1::CbmL1(const char *name, const char *title):CbmTask(name)
{
  ref= new TRefArray();
}

void CbmL1::SetParContainers()
{
  CbmRunAna* ana = CbmRunAna::Instance();
  CbmRuntimeDb* rtdb=ana->GetRuntimeDb();
  rtdb->getContainer("CbmGeoPassivePar");
  rtdb->getContainer("CbmGeoStsPar");
  rtdb->getContainer("CbmGeoTrdPar");
}

InitStatus CbmL1::ReInit()
{
  SetParContainers();
  return Init();
}

InitStatus CbmL1::Init()
{
  CbmRootManager *fManger = CbmRootManager::Instance();

  histodir = gROOT->mkdir("L1histo");

  listStsPts = (TClonesArray *)  fManger->GetObject("STSPoint");
  listTrdPts = 0;//(TClonesArray *)  fManger->GetObject("TRDPoint");

  listMCTracks = (TClonesArray*) fManger->GetObject("MCTrack");
  
  listStsHits = (TClonesArray *)  fManger->GetObject("STSHit");

  //Kalman.Init();

  NStation = CbmKF::Instance()->StsStationIDMap.size();
  stat_fit_time = 0;
  stat_fit_time_ntracks = 0;

  return kSUCCESS;
}


void CbmL1::Exec(Option_t * option)
{
  static int nevent=1;
  cout << endl << "CbmL1::Exec event " << nevent++ << " ..." << endl << endl;
        
  ReadEvent();  

  cout<<"L1 Track finder..."<<endl;
  CATrackFinder();
  //MCTrackFinder();

  cout<<"L1 Track fitter..."<<endl;

  TrackFitter( vRTracks, &PrimVtx );

  cout<<"L1 Performance..."<<endl;
  Performance();

  cout<<"End of L1"<<endl;

  static bool ask = 0; 
  char symbol;
  if (ask){
    std::cout << std::endl << "L1 run (any/r/q) > ";
    do{
      std::cin.get(symbol);
      if (symbol == 'r')
	ask = false;
      if ((symbol == 'e')||(symbol == 'q'))
	exit(0);
    } while (symbol != '\n');
  }
}


void CbmL1::ReadStsHit( CbmL1StsHit *h, int index )
{
  CbmStsHit *mh = (CbmStsHit*) listStsHits->At(index);

  h->ExtIndex = index;
  CbmKF &KF = *CbmKF::Instance();
  int id = 1000+mh->GetStationNr();

  if( KF.StsStationIDMap.find(id) != KF.StsStationIDMap.end())
    h->iStation = KF.StsStationIDMap[id];
  else  h->iStation = 0;

  h->iSector  = mh->GetSectorNr(); 
  h->isStrip  = 0;
  h->iStripF = mh->GetDigi(0);
  h->iStripB = mh->GetDigi(1);
  if( h->iStripF>=0 && h->iStripB>=0 ) h->isStrip  = 1;
  h->indStripF = h->indStripB = -1;
  h->first_hit     = -1; h->last_hit      = -1;
  h->first_triplet = -1; h->last_triplet  = -1;
  h->level = -2;
  h->used = false;
  
  h->KFStsHit.Create(mh);
  
  TVector3 pos, err;
  mh->Position(pos);
  mh->PositionError(err);

  h->X = pos.X();
  h->Y = pos.Y();
  h->Z = pos.Z();
  h->MC_Point = -1;
  h->iPileUp = 0;
}


bool CbmL1::ReadMCPoint( CbmL1MCPoint *MC, int iPoint, int iPileUp )
{
  CbmStsPoint *pt = (CbmStsPoint*) listStsPts->At(iPoint);
  if ( !pt ) return 1;
  TVector3 xyz,P;
  pt->Position(xyz);
  pt->Momentum(P);
  MC->iPileUp = iPileUp;
  MC->x  = xyz.X();
  MC->y  = xyz.Y();
  MC->z  = xyz.Z();
  MC->px = P.X();
  MC->py = P.Y();
  MC->pz = P.Z();  
  MC->p = sqrt( MC->px*MC->px + MC->py*MC->py + MC->pz*MC->pz );
  MC->ID = pt->GetTrackID();
  MC->UID = MC->ID*100+MC->iPileUp;

  CbmMCTrack *MCTrack = (CbmMCTrack*) listMCTracks->At( MC->ID );
  if ( !MCTrack ) return 1;	  
  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)     ) return 1;
  MC->q = TDatabasePDG::Instance()->GetParticle(MC->pdg)->Charge()/3.0;
  MC->mass = TDatabasePDG::Instance()->GetParticle(MC->pdg)->Mass();      
  MC->T[0] = MC->x;
  MC->T[1] = MC->y;
  MC->T[2] = MC->px/MC->pz;
  MC->T[3] = MC->py/MC->pz;
  MC->T[4] = MC->q/MC->p;
  MC->T[5] = MC->z;

  pt->PositionOut(xyz);
  pt->MomentumOut(P);

  MC->TOut[0] = xyz.X();
  MC->TOut[1] = xyz.Y();
  MC->TOut[2] = P.X()/P.Z();
  MC->TOut[3] = P.Y()/P.Z();
  MC->TOut[4] = MC->q/P.Mag();
  MC->TOut[5] = xyz.Z();

  return 0;
}