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

#include "L1Algo.h"

#include "L1Branch.h" // contain L1Track

#include "CbmL1MCTrack.h"
#include "L1Extrapolation.h"
#include "KFPVertex.h"
#include "TDHelper.h"
#include "TMath.h"
#include "KFMCTrack.h"

namespace CbmL1Constants {
   /// Performance constants
   const double MinRecoMom = 0.1;   // Extra set of tracks = (MinRecoMom, MinRefMom)
   const double MinRefMom = 1.;     // Primary set of tracks = (MinRefMom, +inf)  //All reco tracks = (MinRecoMom, +inf)
   const double MinPurity = 0.7;    // min NStationInRecoTrack/NStationInMCTrack
   const double MinNStations = 4.;  // min number of stations in track to be reconstructable
}



#include <cstring>
#include <iostream>
#include <fstream>
#include <vector>
#include <map>
//ROOT includes 
#ifdef ROOT_IS_USED
#include "TFile.h"
#include "TDirectory.h"
#include "TH1F.h"
#endif
  // 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.),
                           algo(0), fKFParticles(0),
                           L1_NEVENTS(0), L1_CATIME(0.)
{
#ifdef ROOT_IS_USED
  fPerfFile = new TFile("FLESQA.root","RECREATE");
  
  TDirectory *curdir = gDirectory;
  
  fPerfFile->cd();
  fHistoDir = fPerfFile->mkdir( "FLESQA" );
  fHistoDir->cd();
  
  gDirectory->mkdir("TimeHisto");
  gDirectory->cd("TimeHisto");
  {
    TH1FParameters timeHisto[fNTimeHistos] = 
    {
      {"NTracksVsTime","NTracksVsTime", 12100, -100.f, 12000.f},
      {"TracksTimeResidual","TracksTimeResidual", 300,-30.f,30.f},
      {"TracksTimePull","TracksTimePull",100,-10.f,10.f},
      {"NHitsVsTime","NHitsVsTime", 12100, -100.f, 12000.f},
      {"HitsTimeResidual","HitsTimeResidual", 300,-30.f,30.f},
      {"HitsTimePull","HitsTimePull",100,-10.f,10.f},
      {"NTracksInEventsVsTime","NTracksInEventsVsTime", 12100, -100.f, 12000.f},
      {"NHitsInEventsVsTime","NHitsInEventsVsTime", 12100, -100.f, 12000.f},
      {"NTracksVsMCTime","NTracksVsMCTime", 12100, -100.f, 12000.f},
      {"NHitsVsMCTime","NHitsVsMCTime", 12100, -100.f, 12000.f},
      {"dtDistribution","dtDistribution", 100, 0.f, 1000.f},
      {"NTracksInEventsVsTime1","NTracksInEventsVsTime1", 12100, -100.f, 12000.f},
      {"NTracksInEventsVsTime2","NTracksInEventsVsTime2", 12100, -100.f, 12000.f},
      {"NTracksInEventsVsTime3","NTracksInEventsVsTime3", 12100, -100.f, 12000.f},
      {"NTracksInEventsVsTime4","NTracksInEventsVsTime4", 12100, -100.f, 12000.f},
      {"NTracksInEventsVsTime5","NTracksInEventsVsTime5", 12100, -100.f, 12000.f},
        
      {"NHitsInEventsVsTime1","NHitsInEventsVsTime1", 12100, -100.f, 12000.f},
      {"NHitsInEventsVsTime2","NHitsInEventsVsTime2", 12100, -100.f, 12000.f},
      {"NHitsInEventsVsTime3","NHitsInEventsVsTime3", 12100, -100.f, 12000.f},
      {"NHitsInEventsVsTime4","NHitsInEventsVsTime4", 12100, -100.f, 12000.f},
      {"NHitsInEventsVsTime5","NHitsInEventsVsTime5", 12100, -100.f, 12000.f},
        
      {"Number of merged events","Number of merged events", 6, -0.5f, 5.5f}
    };
    for(int iH=0; iH<fNTimeHistos; iH++)
      fTimeHisto[iH] = new TH1F(timeHisto[iH].name.Data(),
                                timeHisto[iH].title.Data(),
                                timeHisto[iH].nbins,
                                timeHisto[iH].xMin,
                                timeHisto[iH].xMax );

    fTimeHisto[0]->GetXaxis()->SetTitle("Time, ns");
    fTimeHisto[0]->GetYaxis()->SetTitle("Number of tracks");
    fTimeHisto[16]->GetYaxis()->SetTitle("Number of Events");
    fTimeHisto[16]->GetXaxis()->SetTitle("Number of MCEvents in Event");
  
    fTimeHisto[11]->SetLineColor(5);
    fTimeHisto[11]->SetFillColor(5);
    fTimeHisto[12]->SetLineColor(2);
    fTimeHisto[12]->SetFillColor(2);
    fTimeHisto[13]->SetLineColor(3);
    fTimeHisto[13]->SetFillColor(3);
    fTimeHisto[14]->SetLineColor(7);
    fTimeHisto[14]->SetFillColor(7);
    fTimeHisto[15]->SetLineColor(6);
    fTimeHisto[15]->SetFillColor(6);
      
    fTimeHisto[16]->SetLineColor(5);
    fTimeHisto[16]->SetFillColor(5);
    fTimeHisto[17]->SetLineColor(2);
    fTimeHisto[17]->SetFillColor(2);
    fTimeHisto[18]->SetLineColor(3);
    fTimeHisto[18]->SetFillColor(3);
    fTimeHisto[19]->SetLineColor(7);
    fTimeHisto[19]->SetFillColor(7);
    fTimeHisto[20]->SetLineColor(6);
    fTimeHisto[20]->SetFillColor(6);
  }
  gDirectory->cd(".."); //FLESQA
  
  gDirectory->mkdir("Ghost");
  gDirectory->cd("Ghost");
  {
    TH1FParameters histoParam[fNGhostHistos] = 
    {
      {"h_ghost_mom", "Momentum of ghost tracks", 500, 0.0, 5.0},
      {"h_ghost_nhits", "Number of hits in ghost track", 51, -0.1, 10.1},
      {"h_ghost_fstation", "First station of ghost track", 50, -0.5, 10.0},
      {"h_ghost_chi2", "Chi2/NDF of ghost track", 50, -0.5, 10.0},
      {"h_ghost_prob", "Prob of ghost track", 505, 0., 1.01},
      {"h_ghost_r", "R of ghost track at the first hit", 50, 0.0, 15.0},
      {"h_ghost_tx", "tx of ghost track at the first hit", 50, -5.0, 5.0},
      {"h_ghost_ty", "ty of ghost track at the first hit", 50, -1.0, 1.0},
      {"h_ghost_Rmom", "Ghost tracks", 100, 0.0, 5.0},
    };
    
    for(int iH=0; iH<fNGhostHistos; iH++)
      fGhostHisto[iH] = new TH1F(histoParam[iH].name.Data(),
                                 histoParam[iH].title.Data(),
                                 histoParam[iH].nbins,
                                 histoParam[iH].xMin,
                                 histoParam[iH].xMax );
  }
  gDirectory->cd(".."); //FLESQA
  
  curdir->cd();  
#endif
  
  fTopoPerformance.CreateHistos("TopoPerformance",fPerfFile);

  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()
{
#ifdef ROOT_IS_USED

  TDirectory *curr = gDirectory;
  TFile *currentFile = gFile;
  // Open output file and write histograms

  fPerfFile->cd();
  WriteHistosCurFile(fHistoDir);
  WriteHistosCurFile(fTopoPerformance.GetHistosDirectory());
  fPerfFile->Close();
  fPerfFile->Delete();

  gFile = currentFile;
  gDirectory = curr;
#endif
}

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_, int iVerbose_)
{
//  if (iVerbose >= 1) cout << " -I- run Performance::Init()" << endl;  

  algo = algoint_;
  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::SetRecoEvents (const vector< L1Event >  *Event_ )
{
 // int nEvent = Event_->size();
 // vEvent->resize(nEvent);
 // vEvent.resize(Event_.size());
  //for(unsigned int i=0; i<Event_->size(); ++i )  vEvent[i] = Event_[i];
  vEvent = Event_;
}

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 KFPVertex& 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  = vRTracks[iR].StsHits.size(); // number of hits in track
   
    // count how many hits from each mcTrack belong to current recoTrack
    map<int, int > &hitmap = vRTracks[iR].hitMap; // how many hits from each mcTrack belong to current recoTrack
    for (vector<int>::iterator ih  = (vRTracks[iR].StsHits).begin(); ih != (vRTracks[iR].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);

//         std::cout<<" prtraMC "<<prtra->GetMCTracks().size()<<std::endl;
       
      } 
      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();

  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() && ! mtra.IsAdditional() ) 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;
    double mc_length_hits = mtra.NStations();
    int mc_length = mtra.NMCStations();
    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;
//     }
    
    if ( mtra.IsAdditional() ){ // short
      ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "short");
      switch ( mtra.pdg ) {
        case 11:
        case -11:
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortE");
          break;
        case 211:
        case -211:
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortPion");
          break;
        case 321:
        case -321:
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortKaon");
          break;
        case 2212:
        case -2212:
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortProton");
          break;
        default:
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "shortRest");
      }
    }
    else { // separate all efficiecies from short eff
    
      ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "total");

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

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

      if ( mtra.pdg == 11 || mtra.pdg == -11 ) {
        ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "total_e");

        if ( mtra.p > CbmL1Constants::MinRefMom ){                        // reference tracks
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_e");
      
          if ( mtra.IsPrimary() ){                         // reference primary
          }
          else{                                             // reference secondary
            ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "fast_sec_e");
          }
        }
        else{                                               // extra set of tracks
          ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_e");
       
          if ( mtra.IsPrimary() ){             // extra primary
          }
          else{
            ntra.Inc(reco, killed, ratio_length, ratio_fakes, nclones, mc_length, mc_length_hits, "slow_sec_e");
          }
        } // 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();
  }
  
  std::cout << "Event reconstruction efficiency" << std::endl;
  fEventEfficiency.PrintEff();
  cout << endl;
  
  cout.precision(7);
  cout.setf(ios::scientific);
}

    /// 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);  

  float covPV[6];
  fPV->GetCovarianceMatrix(covPV);
  
  fout << 0 << " " << 0 << " " << 0 << endl;
  fout << fPV->GetX() << " " << fPV->GetY() << " " << fPV->GetZ() << endl;
  for( int k = 0; k < 6; k++ ) fout << covPV[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);  

  if(fKFParticles)
  {
    const unsigned int NParticles = fKFParticles->size();
    for( unsigned int iT = 0; iT < NParticles; ++iT ) {
      float 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();
}

#ifdef ROOT_IS_USED


void Performance::TimeHist() 
{  
   for(unsigned int i=0;i<vStsHits.size();i++)
   {
     int hitId = vStsHits[i].hitId;
     const L1StsHit& recoHit = algo->vStsHits->at(hitId);
    // int nEvent = recoHit.n;
   //  vStsHits[i].n = nEvent;
   }
  
  for(unsigned int i=0;i<vStsHits.size();i++)
  {
    int hitId = vStsHits[i].hitId;
    const L1StsHit& recoHit = algo->vStsHits->at(hitId);
    
    fTimeHisto[3]->Fill(recoHit.t_reco);
      
    int iCol = 0;// (recoHit.n)%5;
    
    fTimeHisto[16+iCol]->Fill(recoHit.t_reco);
    
    if(vStsHits[i].mcPointIds.size() == 0) continue;
    int mcPointId = vStsHits[i].mcPointIds[0];
    double mcTime = vMCPoints[mcPointId].time;
    double residual = recoHit.t_reco - mcTime;
    double pull =residual/5.;
    fTimeHisto[4]->Fill(residual);
    fTimeHisto[5]->Fill(pull);
  }
  
  float tEvent = 0.f;
  int eventId = -1;
  if(vMCTracks.size()>0)
    tEvent = vMCTracks[0].time;
  for(unsigned int i=0;i<vMCTracks.size();i++)
  {
    fTimeHisto[8]->Fill(vMCTracks[i].time);
      
    if(vMCTracks[i].mother_ID<eventId)
    {
      fTimeHisto[10]->Fill( vMCTracks[i].time - tEvent );
      tEvent = vMCTracks[i].time;
      eventId = vMCTracks[i].mother_ID;
    }
  }
  for(unsigned int i=0;i<vMCPoints.size();i++) 
    fTimeHisto[9]->Fill(vMCPoints[i].time);
  
  for(unsigned int i=0;i<vRTracks.size();i++)
  {
    fTimeHisto[0]->Fill(vRTracks[i].fTrackTime);
    
    if(vRTracks[i].IsGhost())
    {
      if( fabs(vRTracks[i].T[4])>1.e-10) 
      {
        fGhostHisto[0]->Fill(fabs(1.0/vRTracks[i].T[4]));
        fGhostHisto[8]->Fill(fabs(1.0/vRTracks[i].T[4]));
      }
      fGhostHisto[1]->Fill((vRTracks[i].StsHits).size());
      CbmL1HitStore &h1 = vHitStore[vRTracks[i].StsHits[0]];
      CbmL1HitStore &h2 = vHitStore[vRTracks[i].StsHits[1]];
      fGhostHisto[2]->Fill(h1.iStation);
      fGhostHisto[5]->Fill(sqrt(fabs(h1.x*h1.x+h1.y*h1.y)));
      double z1 = algo->vStations[h1.iStation].z[0];
      double z2 = algo->vStations[h2.iStation].z[0];
      if( fabs(z2-z1)>1.e-4 )
      {
        fGhostHisto[6]->Fill((h2.x-h1.x)/(z2-z1));
        fGhostHisto[7]->Fill((h2.y-h1.y)/(z2-z1));
      }
      
      if (vRTracks[i].NDF > 0)
      {
        fGhostHisto[3]->Fill(vRTracks[i].chi2/vRTracks[i].NDF);
        fGhostHisto[4]->Fill(TMath::Prob( vRTracks[i].chi2, vRTracks[i].NDF ));   
      }
      
      continue;
    }
    
    if(vRTracks[i].GetMaxPurity()<0.7) continue;
    CbmL1MCTrack* mcTrack = vRTracks[i].GetMCTracks()[0];
    double mcTime = mcTrack->time;
    double residual = vRTracks[i].fTrackTime - mcTime;
    double pull =residual/5.;
    fTimeHisto[1]->Fill(residual);
    fTimeHisto[2]->Fill(pull);    
  }
 
  for(unsigned int k=0;k<vEvent->size();k++)
  {
    const KFPTrackVector& trackId1 = vEvent->at(0).getTracks();
    int iCol = k%5;

    for(unsigned int i=0;i<vEvent->at(k).getHits().size();i++)  
    { 
      int hitId = vEvent->at(k).getHits()[i];
      const L1StsHit& recoHit = algo->vStsHits->at(hitId);
      fTimeHisto[7]->Fill(recoHit.t_reco);
      fTimeHisto[11+iCol]->Fill(recoHit.t_reco);
    }
      
    const KFPTrackVector& trackId = vEvent->at(k).getTracks();
      
    vector<int>id;
    for(unsigned int i=0;i<trackId.Id().size();i++)
    {  
      CbmL1Track& recoTrack = vRTracks.at(trackId.Id()[i]);
      fTimeHisto[6]->Fill(recoTrack.fTrackTime);
    }
  }
    
  vector<KFMCTrack> mcTracks(vMCTracks.size());
  for(unsigned int iMC=0; iMC<vMCTracks.size(); iMC++)
  {
    mcTracks[iMC].SetX ( vMCTracks[iMC].x );
    mcTracks[iMC].SetY ( vMCTracks[iMC].y );
    mcTracks[iMC].SetZ ( vMCTracks[iMC].z );
    mcTracks[iMC].SetPx( vMCTracks[iMC].px );
    mcTracks[iMC].SetPy( vMCTracks[iMC].py );
    mcTracks[iMC].SetPz( vMCTracks[iMC].pz );
      
    double p = vMCTracks[iMC].p;
    double q = vMCTracks[iMC].q;
      
    mcTracks[iMC].SetMotherId(vMCTracks[iMC].mother_ID);
    mcTracks[iMC].SetQP(q/p);
    mcTracks[iMC].SetPDG(vMCTracks[iMC].pdg);  
  }
    
  vector<int> trackMatch(vRTracks.size());
    
  for(unsigned int i=0;i<vRTracks.size();i++)
  { 
    CbmL1Track& recoTrack = vRTracks[i];
      
    if(recoTrack.GetMCTracks().size()>0)
    {
      CbmL1MCTrack* mcTrack = recoTrack.GetMCTracks()[0];
      trackMatch[i] = mcTrack->ID;
      mcTracks[mcTrack->ID].SetReconstructed();
    }
    else trackMatch[i] = -1;
  }
  KFParticleTopoReconstructor topoReconstructor;
  vector<int> PVTrackIndexArray;
  int indexAdd = 0;
  
  for(unsigned int iEvent=0; iEvent<vEvent->size(); iEvent++)
  {
    const KFParticleTopoReconstructor* eventTR = vEvent->at(iEvent).getTopoReconstructor();
    for( int iPV=0;iPV<eventTR->NPrimaryVertices();iPV++)
    {
      PVTrackIndexArray = eventTR->GetPVTrackIndexArray(iPV);
      for(unsigned int iTr=0; iTr<PVTrackIndexArray.size(); iTr++)
        PVTrackIndexArray[iTr]=PVTrackIndexArray[iTr] + indexAdd;
            
      topoReconstructor.AddPV(eventTR->GetPrimKFVertex(iPV),PVTrackIndexArray);
      PVTrackIndexArray.clear();
    }
    
    for(unsigned int iP=0;iP<eventTR->GetParticles().size();iP++)
    {
      KFParticle particle;
      const KFParticle& particleEvent =eventTR->GetParticles()[iP];
      particle = particleEvent;
      particle.CleanDaughtersId();
      int idP = particleEvent.Id()+indexAdd;
      int idDaughter=0;
      for(int nD=0;nD<particleEvent.NDaughters();nD++)
      {
//         cout<< " idDaughter "<<particleEvent.DaughterIds()[nD] <<endl;
        if(particleEvent.NDaughters()==1) idDaughter = particleEvent.DaughterIds()[nD];
        if(particleEvent.NDaughters()>1) idDaughter = particleEvent.DaughterIds()[nD]+indexAdd;
//          cout<< " idDaughter "<< idDaughter<<" particleEvent.DaughterIds()[nD] "<<particleEvent.DaughterIds()[nD]<<" nD "<<nD<<endl;
        particle.AddDaughterId(idDaughter);

      }
      particle.SetId(idP);
      topoReconstructor.AddParticle(particle);
    }
    indexAdd +=eventTR->GetParticles().size();
  }

  fTopoPerformance.SetMCTracks(mcTracks);
  fTopoPerformance.SetTrackMatch(trackMatch);
  fTopoPerformance.SetTopoReconstructor(&topoReconstructor);
  
  fTopoPerformance.CheckMCTracks();
  fTopoPerformance.MatchTracks();
  fTopoPerformance.FillHistos();
  

  int motherId=0;
  int mcId=0;
  int iCurrentEvent = -1;
  
  fMCEvents.clear();
  for(unsigned int iMC=0;iMC<vMCTracks.size(); iMC++ )
  {
    motherId = vMCTracks[iMC].mother_ID;
    mcId = vMCTracks[iMC].ID;
    while(motherId>=0)
      motherId = vMCTracks[motherId].mother_ID;

    int iEvent = -motherId-1;
    if(iEvent > iCurrentEvent)
    {
      L1MCEvent mcEvent;
      mcEvent.SetId(iEvent);
      fMCEvents.push_back(mcEvent);
      iCurrentEvent = iEvent;
    }
    if(iEvent != fMCEvents.size() - 1)
    {
      std::cout << "MC tracks are not sorted according to the event number." << std::endl;
      std::cout << " motherId " << motherId << "   iEvent " << iEvent << " " << "  fMCEvents.size() - 1: "<< (fMCEvents.size() - 1)  << std::endl;
      exit(0);
    }
      
    vector<int> & track = fMCEvents[iEvent].GetMCTrackIds();
    track.push_back(mcId);
    
    if(vMCTracks[iMC].IsReconstructable()){fMCEvents[iEvent].SetReconstructable(1);}
  }
  
  
  int idMcTrack=0;
  L1EventMatch EventMatch;
  vector<int> tracksId;
  for(unsigned int iEv=0;iEv<vEvent->size();iEv++)
  {
    EventMatch.SetNEventTracks(vEvent->at(iEv).getTracks().Size());
//       tracksId.resize(vEvent->at(iEv).getTracks().Size());
    for(int iTr=0;iTr<vEvent->at(iEv).getTracks().Size();iTr++)
    {
      CbmL1Track& recoTrack = vRTracks.at(vEvent->at(iEv).getTracks().Id()[iTr]);
      if(recoTrack.GetMCTracks().size()>0)
      {
        CbmL1MCTrack* mcTrack = recoTrack.GetMCTracks()[0];
        idMcTrack = mcTrack->ID;
        motherId = vMCTracks[idMcTrack].mother_ID;
        while(motherId>=0)
          motherId = vMCTracks[motherId].mother_ID;
        int iMCEvent = -motherId-1;
        EventMatch.AddTrack(iMCEvent);
      } 
      tracksId.push_back(vEvent->at(iEv).getTracks().Id()[iTr]);
    }
    EventMatch.SetTracks(tracksId);
    fEventMatch.push_back(EventMatch);
    
    for(std::map<int,int>::iterator it=EventMatch.GetMCEvents().begin(); it!=EventMatch.GetMCEvents().end(); ++it)
      fMCEvents[it->first].AddRecoEvent(iEv);

    EventMatch.Clear(); 
  }
  
  static int NEventsWith2MC = 0;
  static int NEventsWith2Reco = 0;
  
  bool is2MCReconstructable = true;
  bool is2MCReconstructed = true;
  
  for(unsigned int iMCPV=0; iMCPV<fTopoPerformance.GetMCtoRPVId().size(); iMCPV++)
  {
    is2MCReconstructable &= fTopoPerformance.GetPrimVertices()[iMCPV].IsReconstructable();
    is2MCReconstructed &= fTopoPerformance.GetPrimVertices()[iMCPV].IsReconstructed();
  }
  
  if(is2MCReconstructable && fTopoPerformance.GetMCtoRPVId().size() == 2) 
  {
    NEventsWith2MC++;
    if(is2MCReconstructed) 
      NEventsWith2Reco++;
  }
  
  std::cout << "NEventsWith2MC " << NEventsWith2MC << " NEventsWith2Reco " << NEventsWith2Reco << " efficiency " << float(NEventsWith2Reco)/float(NEventsWith2MC) << std::endl;
    
//     float r[2][3];  //distance between MCPV
//     
//     for(unsigned int iMCPV=0; iMCPV<fTopoPerformance.GetMCtoRPVId().size(); iMCPV++)
//     {
//       std::cout<<"MCPV num: "<< iMCPV <<"  ";
//       std::cout<<"N Tracks in MCVertex: "<<fTopoPerformance.GetPrimVertices()[iMCPV].NDaughterTracks()<<"  ";
//       std::cout<<"N RecTracks in MCVertex: "<<fTopoPerformance.GetPrimVertices()[iMCPV].NReconstructedDaughterTracks()<<"  ";
//       std::cout<<"MC PV x: "<< fTopoPerformance.GetPrimVertices()[iMCPV].X()<<" y "<<fTopoPerformance.GetPrimVertices()[iMCPV].Y()<<" z "<<fTopoPerformance.GetPrimVertices()[iMCPV].Z()<<"  ";                
//       std::cout<<std::endl;
//       
//       r[iMCPV][0] = fTopoPerformance.GetPrimVertices()[iMCPV].X();
//       r[iMCPV][1] = fTopoPerformance.GetPrimVertices()[iMCPV].Y();
//       r[iMCPV][2] = fTopoPerformance.GetPrimVertices()[iMCPV].Z();
//       
//       if(fTopoPerformance.GetMCtoRPVId()[iMCPV].IsMatched())  
//       {  
//         std::cout<<std::endl<<std::endl;
//         for(unsigned int iRecoPV=0; iRecoPV<fTopoPerformance.GetMCtoRPVId()[iMCPV].ids.size(); iRecoPV++)
//         {
//           KFParticle pvTmp = topoReconstructor.GetPrimVertex( fTopoPerformance.GetMCtoRPVId()[iMCPV].ids[iRecoPV] );
//           std::cout<<"     PV Number in Event "<<fTopoPerformance.GetMCtoRPVId()[iMCPV].ids[iRecoPV]<<"  ";
//           std::cout<<"Quantity of RecTracks in PV: "<<topoReconstructor.GetPVTrackIndexArray(fTopoPerformance.GetMCtoRPVId()[iMCPV].ids[iRecoPV]).size()<<"  ";
//           std::cout<<"PV: x: "<<pvTmp.X()<<" y: "<<pvTmp.Y()<<" z: "<<pvTmp.Z()<<std::endl;
//           if (pvTmp.Chi2()/pvTmp.NDF()>=4) 
//             std::cout<<" Chi2/NDF "<< pvTmp.Chi2()/pvTmp.NDF() << std::endl;
//           std::cout<<std::endl;
//         }
//       } 
//      
//       std::cout << std::endl;
//       if ( fTopoPerformance.GetPrimVertices()[iMCPV].IsReconstructable() && !fTopoPerformance.GetPrimVertices()[iMCPV].IsReconstructed()) std::cin.get();
//     }
//      std::cout << " distance " << sqrt( (r[0][0] - r[1][0])*(r[0][0] - r[1][0]) + (r[0][1] - r[1][1])*(r[0][1] - r[1][1]) + (r[0][2] - r[1][2])*(r[0][2] - r[1][2]) ) << " cm " << std::endl << std::endl; 
//   for( int iPV=0; iPV<topoReconstructor.NPrimaryVertices(); iPV++)
//   {
//     std::cout << "PV " << iPV << std::endl;
//     for(unsigned int iTrack=0; iTrack<topoReconstructor.GetPVTrackIndexArray(iPV).size(); iTrack++)
//     {
//       int particleIndex = topoReconstructor.GetPVTrackIndexArray(iPV)[iTrack];
//       int trackIndex = topoReconstructor.GetParticles()[particleIndex].DaughterIds()[0];
//       const KFMCTrack& mcTrack = fTopoPerformance.GetMCTrack(particleIndex);
//       std::cout << "iPVTrack " << iTrack << " iParticle " << particleIndex <<  " iTrack " << trackIndex << " iMCPV " << mcTrack.MotherId() <<" PDG " << mcTrack.PDG() << std::endl;
//     }
//   }
//   std::cin.get();
}    

void Performance::EventEfficiency() 
{ 
  L1EventEfficiencies eventEfficiency; // efficiencies for current event
  for(unsigned int iEvent=0; iEvent<fEventMatch.size(); iEvent++)
    eventEfficiency.AddGhost( fEventMatch[iEvent].IsGhost() );

  for(unsigned int iMCEvent=0; iMCEvent<fMCEvents.size(); iMCEvent++) {
    if( fMCEvents[iMCEvent].IsReconstructable() )
    {
      const vector< int >& recoEvents = fMCEvents[iMCEvent].GetRecoEvents(); // for length calculations
      bool reco = 0;
      if(recoEvents.size() ==1)
        reco = 1;
      // number of cloned events
      int nclones = 0;
      if( fMCEvents[iMCEvent].NClones()>0 )
        nclones = 1;
      eventEfficiency.Inc(reco,nclones, "reconstructable");
    }
    if( fMCEvents[iMCEvent].IsReconstructed() )
    {
      const vector< int >& recoEvents = fMCEvents[iMCEvent].GetRecoEvents(); // for length calculations
      bool reco = 0;
      if(recoEvents.size() ==1)
        reco = 1;
      // number of cloned events
      int nclones = 0;
      if( fMCEvents[iMCEvent].NClones()>0 )
        nclones = 1;
      eventEfficiency.Inc(reco,nclones, "reconstructed");
    }
  } // for mcTracks
  eventEfficiency.IncNEvents();
  fEventEfficiency += eventEfficiency;
  eventEfficiency.CalcEff();
  fEventEfficiency.CalcEff();
  
  for(int iM=0;iM<fEventMatch.size();iM++)
    fTimeHisto[fNTimeHistos-1]->Fill(fEventMatch[iM].NMCEvents());  
  
//   output Reco vs MC PV
//   int iDMCTrack;
//   for(int iM=0;iM<fEventMatch.size();iM++)
//   {
//     fTimeHisto[fNTimeHistos-1]->Fill(fEventMatch[iM].NMCEvents());
//     if(fEventMatch[iM].NMCEvents()>=2)
//     {
//       int ntracks;
//       ntracks = fEventMatch[iM].NRecoTracks();
//       for(std::map<int,int>::iterator it=fEventMatch[iM].GetMCEvents().begin(); it!=fEventMatch[iM].GetMCEvents().end(); ++it)
//       {
//         
//         for(int q=0;q<fMCEvents.size();q++){
//           if(fMCEvents[q].GetId()==it->first)
//           {
//             std::cout<<"Event number  "<< iM<<"  MCEvent  "<<fMCEvents[q].GetId()<<std::endl;
//             int index = fMCEvents[q].GetMCTrackIds()[0];
//             std::cout<<"MCVertex coord: "<< vMCTracks[index].x<<" "<<vMCTracks[index].y<<" "<<vMCTracks[index].z<<std::endl;
//             int n=0;
//             for(int p=0;p<fEventMatch[iM].GetTracksId().size();p++)
//             {
//               CbmL1Track& recoTrack = vRTracks.at(fEventMatch[iM].GetTracksId()[p]);
//               if(recoTrack.GetMCTracks().size()>0)
//               {
//                 CbmL1MCTrack* mcTrack = recoTrack.GetMCTracks()[0];
//                 iDMCTrack = mcTrack->ID;
//                 for(int w=0;w<fMCEvents[q].GetMCTrackIds().size();w++)
//                 {
//                   if(fMCEvents[q].GetMCTrackIds()[w]==iDMCTrack){n++;}
//                 }
//               }
//             }
//             std::cout<<"Quantity of MC Tracks in Event   "<<n<< std::endl<<std::endl;
//           }
//         }  
//       }
//       std::cout<<"Quantity of RecVertices in Event   "<<(*vEvent)[iM].getTopoReconstructor()->NPrimaryVertices()<<std::endl;
//       std::cout<<"RecVertex coord:"<<std::endl;
//         for(int iPV=0; iPV<(*vEvent)[iM].getTopoReconstructor()->NPrimaryVertices(); iPV++)
//         {
//           KFParticle pvTmp = (*vEvent)[iM].getTopoReconstructor()->GetPrimVertex(iPV);
//           std::cout<<"  "<<pvTmp.X()<<"  "<<pvTmp.Y()<<"  "<<pvTmp.Z()<<std::endl;
//           std::cout<<"Quantity of RecTracks in PV: "<<(*vEvent)[iM].getTopoReconstructor()->GetPVTrackIndexArray(iPV).size()<< std::endl;
//         }
//       std::cout<<std::endl;
//     }
//   }

//   for(unsigned int iM=0;iM<fMCEvents.size();iM++)
//   {
//     if(fMCEvents[iM].IsReconstructable()) {
//       
//       if(!fMCEvents[iM].IsReconstructed()) {
//         for(unsigned int g=0;g<fMCEvents[iM].GetMCTrackIds().size();g++){
//         CbmL1MCTrack& mcTr = vMCTracks.at(fMCEvents[iM].GetMCTrackIds()[g]);
// 
//         if(mcTr.IsReconstructable()) std::cout<<"TrId "<<mcTr.ID <<"  mcTr   "<<mcTr.IsReconstructable()<<"      " << mcTr.IsReconstructed()<<"           fMCEvents[iM]  "<<fMCEvents[iM].GetId()<<std::endl;}}
//       
//     }
//   }
}

void Performance::FillHist()
{
  TimeHist();
  EventEfficiency();
}

void Performance::WriteHistosCurFile( TObject *obj )
{
  if( !obj->IsFolder() )
    obj->Write();
  else
  {
    TDirectory *cur = gDirectory;
    TFile *currentFile = gFile;

    TDirectory *sub = cur->GetDirectory(obj->GetName());
    sub->cd();
    TList *listSub = (static_cast<TDirectory*>(obj))->GetList();
    TIter it(listSub);
    while( TObject *obj1=it() ) WriteHistosCurFile(obj1);
    cur->cd();
    gFile = currentFile;
    gDirectory = cur;
  }
}
#endif