#include "Performance.h"
#include "PerformanceConstants.h"

#include "L1Algo.h"

#include "L1Branch.h" // contain L1Track

#include "CbmL1MCTrack.h"
#include "L1Extrapolation.h"
#include "CbmKFVertex.h"
#include "TDHelper.h"
#include "TMath.h"


#include <cstring>
#include <iostream>
#include <fstream>
#include <vector>
#include <map>

  // set draw oerformance info for each one seperate event
#define XXX 

using std::cout;
using std::endl;
using std::vector;
using std::string;

Performance *Performance::fInstance = 0;

Performance::Performance():iVerbose(1), nEvent(0), 
//                MinRecoMom(0.1), MinRefMom(1.), MinPurity(0.7), MinNStations(4.),
                L1_NEVENTS(0), L1_CATIME(0.)
{
  p_eff_all_vs_mom.clear(); p_eff_prim_vs_mom.clear();
  p_eff_sec_vs_mom.clear(); p_eff_d0_vs_mom.clear();   
  for( int i = 0; i < 5; i++ )  
    for ( int k = 0; k < 3; k++)
      hist[i][k].clear();

  if( !fInstance ) fInstance = this;
}

Performance *Performance::Instance(){ return fInstance; }

    /// init members
    /// @param algoint_   - interface for L1Algo, witch contain algo
    /// @param _work_dir  - first part of path in witch read and write files (data_perfo.txt, data_algo.txt, geo_algo.txt)
void Performance::Init(const char*  _work_dir, L1Algo* algoint_, vector<CbmKFParticle> *kfpart , int iVerbose_)
{
//  if (iVerbose >= 1) cout << " -I- run Performance::Init()" << endl;  

  algo = algoint_;
  fKFParticles = kfpart;
  iVerbose = iVerbose_;
  strcpy (work_dir,_work_dir);
  nEvent = 1;

  remove( "../Output/track_param.dat" ); // clean
  remove( "../Output/primary_vertex.dat" ); // clean
  remove( "../Output/masses.dat" ); // clean
  remove( "../Output/Prob_reco.dat" ); // clean
  remove( "../Output/Prob_ghost.dat" ); // clean
  remove( "../Output/Prob_reco_3.dat" ); // clean
  remove( "../Output/Prob_reco_4.dat" ); // clean
  remove( "../Output/Prob_ghost_3.dat" ); // clean
  remove( "../Output/Prob_ghost_4.dat" ); // clean

//   if (iVerbose >= 1) cout << " void Performance::Init(const char*  _work_dir, L1AlgoInter *algoint_)   done " << endl;
}

    /// read event data from file, run performance. Before call it you must full the vRTracks
void Performance::RunEv()
{
//  if (iVerbose >= 1) cout << " Perfomance for Event number " << nEvent << endl;
//  ReadStAPPerfoData();
  ExecPerf();

  nEvent++;
}


void Performance::SetData( const vector< CbmL1HitStore > & HitStore_,
                           const vector< CbmL1StsHit >   & StsHits_,
                           const vector< CbmL1MCPoint >  & MCPoints_,
                           const vector< CbmL1MCTrack >  & MCTracks_,
                           const vector< int >           & HitMCRef_ )
{
  vHitStore.resize(HitStore_.size());
  vStsHits.resize(StsHits_.size());
  vMCPoints.resize(MCPoints_.size());
  vMCTracks.resize(MCTracks_.size());
  vHitMCRef.resize(HitMCRef_.size());

  for(unsigned int i=0; i<HitStore_.size(); ++i ) vHitStore[i] = HitStore_[i];
  for(unsigned int i=0; i<StsHits_.size(); ++i )  vStsHits[i] = StsHits_[i];
  for(unsigned int i=0; i<MCPoints_.size(); ++i ) vMCPoints[i] = MCPoints_[i];
  for(unsigned int i=0; i<MCTracks_.size(); ++i ) vMCTracks[i] = MCTracks_[i];
  for(unsigned int i=0; i<HitMCRef_.size(); ++i ) vHitMCRef[i] = HitMCRef_[i];
}


void Performance::SetRecoTracks( const vector< CbmL1Track >  & RTracks_ )
{
  vRTracks.resize(RTracks_.size());
  for(unsigned int i=0; i<RTracks_.size(); ++i )  vRTracks[i] = RTracks_[i];
}
void Performance::SetPrimaryVertex(const CbmKFVertex& PV_) { fPV = &PV_; }

/// read data for performance(vMCPoints, vMCTracks, vHitMCRef, vHitStore) of one event from file

    /// L1 Fit performance
void Performance::ExecPerf(){
  TrackMatch();
  EfficienciesPerformance();
  WriteTrackParamInFile();
  WritePV();
  WriteMasses();
  WriteProb();
}


#include "CbmL1Counters.h"

#include <iostream>
#include <vector>
#include <list>
#include <map>

using std::cout;
using std::endl;
using std::ios;
using std::vector;
using std::pair;
using std::map;

void Performance::TrackMatch()
{
  map <int, CbmL1MCTrack*> pMCTrackMap; pMCTrackMap.clear();

    // fill pMCTrackMap
  for( vector<CbmL1MCTrack>::iterator i = vMCTracks.begin(); i != vMCTracks.end(); ++i){
    CbmL1MCTrack &MC = *i;
    
    if (pMCTrackMap.find(MC.ID) == pMCTrackMap.end()){
      pMCTrackMap.insert(pair<int, CbmL1MCTrack*>(MC.ID, &MC));
    }
    else{
      cout << "*** L1: Track ID " << MC.ID << " encountered twice! ***" << endl;
    }
  }
    // -- prepare information about reconstructed tracks
  const int nRTracks = vRTracks.size();
  for (int iR = 0; iR < nRTracks; iR++){
    CbmL1Track* prtra = &(vRTracks[iR]);
    
    int hitsum  = prtra->StsHits.size(); // number of hits in track

      // count how many hits from each mcTrack belong to current recoTrack
    map<int, int > &hitmap = prtra->hitMap; // how many hits from each mcTrack belong to current recoTrack
    for (vector<int>::iterator ih  = (prtra->StsHits).begin(); ih != (prtra->StsHits).end(); ++ih){

      const int nMCPoints = vStsHits[*ih].mcPointIds.size();
      for (int iP = 0; iP < nMCPoints; iP++){
        int iMC = vStsHits[*ih].mcPointIds[iP];
        int ID = -1;
        if (iMC >= 0) ID = vMCPoints[iMC].ID;
        if(hitmap.find(ID) == hitmap.end())
          hitmap[ID] = 1;
        else{
          hitmap[ID] += 1;
        }
      } // for iPoint
    } // for iHit

      // RTrack <-> MCTrack identification
    double max_percent = 0.0;  // [%]. maximum persent of hits, which belong to one mcTrack.
    for( map<int, int >::iterator posIt = hitmap.begin(); posIt != hitmap.end(); posIt++ ){ // loop over all touched MCTracks
      
      if (posIt->first < 0) continue; // not a MC track - based on fake hits
      
        // count max-purity
      if (double(posIt->second) > max_percent*double(hitsum))
        max_percent = double(posIt->second)/double(hitsum);
      
        // set relation to the mcTrack
      if ( double(posIt->second) > PerformanceConstants::MinPurity * double(hitsum) ){ // found correspondent MCTrack
        if (pMCTrackMap.find(posIt->first) == pMCTrackMap.end()) continue;
        CbmL1MCTrack* pmtra = pMCTrackMap[posIt->first];

        pmtra->AddRecoTrack(prtra);
        prtra->AddMCTrack(pmtra);
      } 
      else{
        if (pMCTrackMap.find(posIt->first) == pMCTrackMap.end()) continue;
        CbmL1MCTrack* pmtra = pMCTrackMap[posIt->first];

        pmtra->AddTouchTrack(prtra);
      }
    } // for hitmap
    prtra->SetMaxPurity(max_percent);
  } // for reco tracks
} // void Performance::TrackMatch() 

void Performance::EfficienciesPerformance()
{
  TL1PerfEfficiencies ntra; // efficiencies for current event

  for (vector<CbmL1Track>::iterator rtraIt = vRTracks.begin(); rtraIt != vRTracks.end(); ++rtraIt){
    ntra.ghosts += rtraIt->IsGhost();
//     if(rtraIt->IsGhost()){ // Debug.
//       cout << " " << rtraIt->GetNOfHits() << " " << 1./rtraIt->T[5] << " " << rtraIt->GetMaxPurity() << " | ";
//       for( map<int, int>::iterator posIt = rtraIt->hitMap.begin(); posIt != rtraIt->hitMap.end(); posIt++ ){
//         cout << " (" << posIt->second << " " << posIt->first << ") ";
//       }
//       cout << endl;
//     }
  }


  int sta_nhits[algo->MaxNStations], sta_nfakes[algo->MaxNStations];

  for (int i = 0; i < algo->NStations; i++){
    sta_nhits[i]  = 0;
    sta_nfakes[i] = 0;
  }

  for ( vector<CbmL1MCTrack>::iterator mtraIt = vMCTracks.begin(); mtraIt != vMCTracks.end(); mtraIt++ ) {
    CbmL1MCTrack &mtra = *(mtraIt);
    
//    if( !( mtra.pdg == -11 && mtra.mother_ID == -1 ) ) continue; // electrons only
    if( ! mtra.IsReconstructable() ) continue;

      // -- find used constans --
      // is track reconstructed
    const bool reco = (mtra.IsReconstructed());
      // is track killed. At least one hit of it belong to some recoTrack
    const bool killed = !reco && mtra.IsDisturbed();
      // ration length for current mcTrack
    vector< CbmL1Track* >& rTracks = mtra.GetRecoTracks(); // for length calculations
    double ratio_length = 0;
    double ratio_fakes  = 0;
    if (reco){
      double mc_length = mtra.NStations();
      for (unsigned int irt = 0; irt < rTracks.size(); irt++) {
        ratio_length += static_cast<double>( rTracks[irt]->GetNOfHits() )*rTracks[irt]->GetMaxPurity() / mc_length;
        ratio_fakes += 1 - rTracks[irt]->GetMaxPurity();
      }
    }
      // number of clones
    int nclones = 0;
    if (reco) nclones = mtra.GetNClones();
//     if (nclones){ // Debug. Look at clones
//       for (int irt = 0; irt < rTracks.size(); irt++){
//         const int ista = vHitStore[rTracks[irt]->StsHits[0]].iStation;
//         cout << rTracks[irt]->GetNOfHits() << "(" << ista << ") ";
//       }
//       cout << mtra.NStations() << endl;
//     }

    ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, "total");
    
    if (( mtra.IsPrimary() )&&(mtra.z > 0)){ // D0
      ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, "d0");
    }

    if ( mtra.p > PerformanceConstants::MinRefMom ){                        // reference tracks
      ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, "fast");
      
      if ( mtra.IsPrimary() ){                         // reference primary
        if ( mtra.NStations() == algo->NStations ){ // long reference primary
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, "long_fast_prim");
        }
        ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, "fast_prim");
      }
      else{                                             // reference secondary
        ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, "fast_sec");
      }
    }
    else{                                               // extra set of tracks
      ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, "slow");
       
      if ( mtra.IsPrimary() ){             // extra primary
        ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, "slow_prim");
      }
      else{
        ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, "slow_sec");
      }
    } // if extra

  } // for mcTracks

  L1_CATIME += algo->CATime;
  L1_NEVENTS++;
  ntra.IncNEvents();
  
  L1_NTRA += ntra;

  ntra.CalcEff();
  L1_NTRA.CalcEff();

} // void CbmL1::Performance()

void Performance::PrintEff() {
  cout.precision(3);
  if( iVerbose ){
    // if( iVerbose > 1 ){
    //   ntra.PrintEff();
    //   cout << "Number of true and fake hits in stations: " << endl;
    //   for (int i = 0; i < algo->NStations; i++){
    //     cout << sta_nhits[i]-sta_nfakes[i] << "+" << sta_nfakes[i] << "   ";
    //   }
    //   cout << endl;
    // } // fVerbose > 1
    cout << endl;
    cout << "L1 ACCUMULATED STAT    : " << L1_NEVENTS << " EVENTS "               << endl << endl;
    L1_NTRA.PrintEff();
    cout << "MC tracks/event found  : " << int(double(L1_NTRA.reco.counters[L1_NTRA.indices["total"]])/double(L1_NEVENTS)) << endl;
    cout<<endl;
    cout<<"CA Track Finder: " << L1_CATIME/L1_NEVENTS << " s/ev" << endl << endl;
	  fstream filestr;
      filestr.open("speedUp.log", fstream::out | fstream::app);
      filestr << (double(L1_NEVENTS)/L1_CATIME) << std::endl;
      filestr.close();
  }
}

    /// write histogram data to file
void Performance::WriteTrackParamInFile(){

  ofstream fout;   // file with histograms data

  static char fname[100];
  
  strcpy(fname, work_dir);
  strcat(fname, "track_param.dat");
  fout.open(fname, ofstream::app);
  
  fout << setprecision(8);
  for( unsigned int iT = 0; iT < vRTracks.size(); ++iT ) {
    CbmL1Track &t = vRTracks[iT];
    if( t.IsGhost() ) continue;
    
    const int iMCF = vHitMCRef[t.StsHits.front()]; // TODO: mcTrack t.GetMCTracks()[0];
    const int iMCL = vHitMCRef[t.StsHits.back()];
    if (iMCF < 0 || iMCL < 0) continue;

    //if ( (t.StsHits.size() < 4) || (t.StsHits.size() >= 8) ) continue;
    //if ( (t.StsHits.size() < 8) ) continue;
    
    fout << t.StsHits.size() << endl;
    {
      CbmL1MCPoint &mc = vMCPoints[iMCF];
        // mc
      const float mcz = mc.zIn;
      const float mcx = mc.xIn;
      const float mcy = mc.yIn;

      fout << mcx << " " << mcy << " " << mcz << " " << mc.pxIn << " " << mc.pyIn << " " << mc.pzIn << " " << mc.q << endl;
      
        // reco
      L1TrackPar T;
      for( int k=0; k<6; k++) (&(T.x))[k] = t.T[k];
      for( int k=0; k<15; k++) (&(T.C00))[k] = t.C[k];
      L1ExtrapolateLine( T, mcz );
      
      for( int k = 0; k < 6; k++ )
        fout << (&(T.x))[k][0] << " ";
      fout << endl;
      for( int k = 0; k < 15; k++ )
        fout << (&(T.C00))[k][0] << " ";
      fout << endl;
      
    }
    {
      CbmL1MCPoint &mc = vMCPoints[iMCL];

      const float mcz = mc.zOut;
      const float mcx = mc.xOut;
      const float mcy = mc.yOut;
      
      fout << mcx << " " << mcy << " " << mcz << " " << mc.pxOut << " " << mc.pyOut << " " << mc.pzOut << " " << mc.q << endl;
      
        // reco
      L1TrackPar T;
      for( int k=0; k<6; k++) (&(T.x))[k] = t.TLast[k];
      for( int k=0; k<15; k++) (&(T.C00))[k] = t.CLast[k];
      L1ExtrapolateLine( T, mcz );
      
      for( int k = 0; k < 6; k++ )
        fout << (&(T.x))[k][0] << " ";
      fout << endl;
      for( int k = 0; k < 15; k++ )
        fout << (&(T.C00))[k][0] << " ";
      fout << endl;
    }
  }
    
  fout.close();
  if ( iVerbose >= 2 ) {
    static int iEvent = 1;
    cout << " -I- Event " << iEvent << ". " << vRTracks.size() << " tracks has been written in " << fname << endl;
    iEvent++;
  }
}

void Performance::WritePV(){

  ofstream fout;   // file with histograms data

  static char fname[100];
  
  strcpy(fname, work_dir);
  strcat(fname, "primary_vertex.dat");
  fout.open(fname,ofstream::app);  

  fout << 0 << " " << 0 << " " << 0 << endl;
  fout << fPV->GetRefX() << " " << fPV->GetRefY() << " " << fPV->GetRefZ() << endl;
  for( int k = 0; k < 6; k++ ) fout << fPV->GetCovMatrix()[k] << " ";
  fout << endl;

  fout.close();
}

void Performance::WriteMasses(){

  ofstream fout;   // file with histograms data

  static char fname[100];
  
  strcpy(fname, work_dir);
  strcat(fname, "masses.dat");
  fout.open(fname,ofstream::app);  

  const unsigned int NParticles = fKFParticles->size();
  for( unsigned int iT = 0; iT < NParticles; ++iT ) {
    Double_t mass, errMass;
    int pdg = fKFParticles->at(iT).GetPDG();
    fKFParticles->at(iT).GetMass(mass,errMass);
    fout << iT << " " << mass << " " << pdg << std::endl;
  }

  fout.close();
}

void Performance::WriteProb(){
  
   const int n_histo_prob = 6;
    ofstream fout[n_histo_prob];   // file with histograms data
    
   string title[n_histo_prob]= {"Prob_reco.dat","Prob_ghost.dat","Prob_reco_3.dat","Prob_ghost_3.dat", "Prob_reco_4.dat",
                                          "Prob_ghost_4.dat"};
    for(int i=0; i<n_histo_prob; i++){
      
      static char fname[100];
      strcpy(fname, work_dir);
      strcat(fname, title[i].data());
      fout[i].open(fname,ofstream::app);  
      
    }

  
   for (vector<CbmL1Track>::iterator rtraIt = vRTracks.begin(); rtraIt != vRTracks.end(); ++rtraIt){
    CbmL1Track* prtra = &(*rtraIt);
    if((prtra->StsHits).size() < 1)continue;

    if ( prtra->IsGhost() ) {
      
      if (prtra->NDF > 0){
 //       h_ghost_chi2->Fill(prtra->chi2/prtra->NDF);
  //      h_ghost_prob->Fill(TMath::Prob( prtra->chi2, prtra->NDF ));
	//float prob_ghost = (TDHelper<float>::Chi2IProbability(  prtra->NDF, prtra->chi2 ));
          float prob_ghost =(TMath::Prob( prtra->chi2, prtra->NDF ));
	fout[1] << prob_ghost <<" "<<(prtra->StsHits).size()<<std::endl;
	
	
	if (((prtra->StsHits).size()<5)&&((prtra->StsHits).size()>3)){
	  
	  fout[5] << prob_ghost <<std::endl;
	}
	  
	if (((prtra->StsHits).size()<4)){
	  
	  fout[3] << prob_ghost <<std::endl;
	}    
      }
    }
    else
    {  // fill histos
      if (prtra->NDF > 0){

      //  float prob_ghost = (TDHelper<float>::Chi2IProbability(  prtra->NDF, prtra->chi2 ));
          float prob_ghost = (TMath::Prob( prtra->chi2, prtra->NDF ));
        fout[0] << prob_ghost <<" "<<(prtra->StsHits).size()<< std::endl;
	
	
        if (((prtra->StsHits).size()<5)&&((prtra->StsHits).size()>3)){
	  
	  fout[4] << prob_ghost <<std::endl;
	}
	  
	if (((prtra->StsHits).size()<4)){
	  
	  fout[2] << prob_ghost <<std::endl;
	
        }
      }
    }
  } // for reco tracks
  
  

  
   for(int i=0; i<n_histo_prob; i++){
      
     fout[i].close();
      
    }
}