#include "Performance.h"
	
#include "../Algo/L1Algo/L1Algo.h"
#include "../Algo/L1AlgoInter.h"

#include "../Algo/L1Algo/L1Branch.h" // contain L1Track

#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;


	  /// 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, L1AlgoInter *algoint_)
{
  cout << " -I- run Performance::Init()" << endl;	

	algo = &(algoint_->algo);
	
	maxNEvent = algoint_->maxNEvent;
	
	iVerbose = algoint_->iVerbose;
	
	strcpy (work_dir,_work_dir);
	
	nEvent = 1;
	
// 	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;
      // collect N events into 1.
//   const int NMergeEvents = 2;
//   const int MCTrStep = 10000;
	ReadStAPPerfoData();
//   Performance perf_tmp = *this;
//   for ( int iEv = 1; iEv < NMergeEvents; iEv++ ){
//     ReadStAPPerfoData();
//     for (int i = 0; i < vMCPoints.size(); i++){
//       CbmL1MCPoint p = vMCPoints[i];
//       p.ID += MCTrStep*iEv;
//       p.mother_ID += MCTrStep*iEv;
//       perf_tmp.vMCPoints.push_back(p);
//     }
//     for (int i = 0; i < vMCTracks.size(); i++){
//       CbmL1MCTrack t = vMCTracks[i];
//       t.ID += MCTrStep*iEv;
//       t.mother_ID += MCTrStep*iEv;
//       perf_tmp.vMCTracks.push_back(t);
//     }
//     for (int i = 0; i < vHitMCRef.size(); i++){
//       int mc = vHitMCRef[i] + MCTrStep*iEv;
//       perf_tmp.vHitMCRef.push_back(mc);
//     }
//     for (int i = 0; i < vHitStore.size(); i++){
//       perf_tmp.vHitStore.push_back(vHitStore[i]);
//     }
//   }
//   *this = perf_tmp;
	GetAlgoResults();
	ExecPerf();

	nEvent++;
};

		/// read data for performance(vMCPoints, vMCTracks, vHitMCRef, vHitStore) of one event from file
void Performance::ReadStAPPerfoData()
{

	static int iEvent = 1; 
	static ifstream fadata;
	static char fname[100];

	if (iEvent <= maxNEvent){
	
		if ( iEvent == 1 ){  			// file open on begin of all work class and close at end
			strcpy(fname, work_dir);
			strcat(fname, "data_perfo.txt");
			fadata.open(fname);	
		};
		vMCPoints.clear();
 		vMCTracks.clear();
		vHitStore.clear();
		vHitMCRef.clear();
			// check if it is right position in file
		char s[] = "Event:  ";  // buffer
		int nEv;  							// event number
		fadata >> s;
		fadata >> nEv;
		if (nEv != iEvent)	cout << " -E- Performance: can't read event number " << iEvent << " from file " << fname << endl;
			// vMCPoints
		int n;  // number of elements
		fadata >> n;
		for (int i = 0; i < n; i++){
			CbmL1MCPoint element;
	
			fadata >> element.x; 
			fadata >> element.y;
			fadata >> element.z;
			fadata >> element.px;
			fadata >> element.py;
			fadata >> element.pz;
			fadata >> element.p;
			fadata >> element.q;
			fadata >> element.mass;
			
			fadata >> element.pdg;
			fadata >> element.ID;
			fadata >> element.mother_ID;
			fadata >> element.iStation; 
			vMCPoints.push_back(element);
		};
		
			// vMCTracks . without Points
		fadata >> n;	
		for (int i = 0; i < n; i++){
			CbmL1MCTrack element;
	
			fadata >> element.x; 				
			fadata >> element.y; 				
			fadata >> element.z; 				
			fadata >> element.px; 			
			fadata >> element.py; 			
			fadata >> element.pz; 			
			fadata >> element.p; 				
			fadata >> element.q; 				
			fadata >> element.mass; 		
			
			fadata >> element.pdg; 			
			fadata >> element.ID; 			
			fadata >> element.mother_ID;
			
			int nhits;
			fadata >> nhits; 		
			for (int k = 0; k < nhits; k++){
				int helement;
				fadata >> helement;
				element.StsHits.push_back(helement);
			};
			
			fadata >> element.nMCContStations; 
			fadata >> element.nHitContStations; 
			vMCTracks.push_back(element);
		};
				
				// find the corrspondece between MCPoint MCTrack
		
		for (unsigned k = 0; k < vMCPoints.size(); k++){
					// find track with ID
				int itr = 0;
				for (;vMCTracks[itr].ID != vMCPoints[k].ID;itr++); 
				vMCTracks[itr].Points.push_back(&vMCPoints[k]);
		};
		
				// vHitMCRef
		fadata >> n;
		for (int i = 0; i < n; i++){
			int element;
			fadata >> element; 
			vHitMCRef.push_back(element);
		};

			// vHitStore
		fadata >> n;									
		for (int i = 0; i < n; i++){
			CbmL1HitStore element;
			fadata >> element.ExtIndex; 
			fadata >> element.iStation; 
			
			fadata >> element.x;  
			fadata >> element.y; 
			vHitStore.push_back(element);
		};
			
		if (iEvent == maxNEvent) { // file open on begin of all work class and close at end
      fadata.close(); 
      cout << " -I- Performance: data read from file " << fname << " successfully"<< endl;
    }
	} // if (iEvent <= maxNEvent)
	iEvent++;
};  // Performance::ReadOneEvData()

		/// fill vRTracks from algo members	
void Performance::GetAlgoResults()
{
    // filling vRTracks
    vRTracks.clear();
    int start_hit = 0; // for interation in algo->vRecoHits[]
    static int ngood = 0, nall = 0;
    for (unsigned i = 0; i < algo->vTracks.size(); i++){
      L1Track it = algo->vTracks[i];
      CbmL1Track t;
      // check for nans
      bool ok = 1;
      nall++;
      for( int k=0; k<6; k++) ok = ok && finite(it.TFirst[k]) && finite(it.TLast[k]);
      for( int k=0; k<15; k++) ok = ok && finite(it.CFirst[k]) && finite(it.CLast[k]);
      ok = ok && finite(it.chi2) && finite(it.NDF);
      ok = ok && it.CFirst[0]>=0 && it.CFirst[2]>=0 && it.CFirst[5]>=0 && it.CFirst[9]>=0 && it.CFirst[14]>=0 &&
	it.CLast[0]>=0 && it.CLast[2]>=0 && it.CLast[5]>=0 && it.CLast[9]>=0 && it.CLast[14]>=0 ;
      if( !ok ) continue;
      ngood++;

      for( int k=0; k<6; k++) t.T[k] = it.TFirst[k];
      for( int k=0; k<15; k++) t.C[k] = it.CFirst[k];
      for( int k=0; k<6; k++) t.TLast[k] = it.TLast[k];
      for( int k=0; k<15; k++) t.CLast[k] = it.CLast[k];
      t.chi2 = it.chi2;
      t.NDF = it.NDF;
      for( int k = 0; k < it.NHits; k++ ){ 
	t.StsHits.push_back( algo->vRecoHits[start_hit++]);
      }
      t.mass = 0.1395679; // pion mass
      t.is_electron = 0;      



      vRTracks.push_back(t);
    } // for vTracks
} // GetAlgoResults()

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

	ofstream faout;   // file with histograms data

	static char fname[100];
	
	strcpy(fname, work_dir);
	strcat(fname, "results_perfo.txt");
	faout.open(fname);	

	
		// Profiles 
	int Nall = p_eff_all_vs_mom.size();
	faout << Nall << endl;
	for (int i = 0; i < Nall; i++){
		faout << p_eff_all_vs_mom[i].x << " " << p_eff_all_vs_mom[i].y << endl; 
	};
	
	int Nprim = p_eff_prim_vs_mom.size();
	faout << Nprim << endl;
	for (int i = 0; i < Nprim; i++){
		faout << p_eff_prim_vs_mom[i].x << " " << p_eff_prim_vs_mom[i].y << endl; 
	};

	int Nsec = p_eff_sec_vs_mom.size();
	faout << Nsec << endl;
	for (int i = 0; i < Nsec; i++){
		faout << p_eff_sec_vs_mom[i].x << " " << p_eff_sec_vs_mom[i].y << endl; 
	};
	
	int Nd0 = p_eff_d0_vs_mom.size();
	faout << Nd0 << endl;
	for (int i = 0; i < Nd0; i++){
	  faout << p_eff_d0_vs_mom[i].x << " " << p_eff_d0_vs_mom[i].y << endl; 
	};
		// Histograms
	int N;
	
	for( int i=0; i<5; i++ ){	
		for ( int k = 0; k < 3; k++){
		
		  N = hist[i][k].size();
			faout << N << endl;
			for (int j = 0; j < N; j++){
				faout << hist[i][k][j].x << endl;
			};	
			
		}		
	}	
	
	faout.close();
	
	cout << " -I- histograms data was written in " << fname << endl;
}





/*
 *====================================================================
 *
 *  CBM Level 1 Reconstruction 
 *  
 *  Authors: I.Kisel,  S.Gorbunov
 *
 *  e-mail : ikisel@kip.uni-heidelberg.de 
 *
 *====================================================================
 *
 *  L1 Fit performance 
 *
 *====================================================================
 */
 
 // iklm
 // remove histogram
 // add output in file some profile histogram

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

					// number of all MCtrack  (see reco)
static int L1_NMTRA_D0              = 0;
static int L1_NMTRA_TOTAL           = 0;
static int L1_NMTRA_FAST            = 0;
static int L1_NMTRA_FAST_PRIM       = 0;  
static int L1_NMTRA_FAST_SEC        = 0;
static int L1_NMTRA_SLOW            = 0; 
static int L1_NMTRA_SLOW_PRIM       = 0;
static int L1_NMTRA_SLOW_SEC        = 0;

					// number of reconstracted
static int L1_NMTRA_D0_RECO         = 0;
static int L1_NMTRA_TOTAL_RECO      = 0;  // MC tracks found 
static int L1_NMTRA_FAST_RECO       = 0;  // all reference tracks:    momentum>1.0 (outputing as Refset)
static int L1_NMTRA_FAST_PRIM_RECO  = 0;  // Reference primary track 			( RefPrim )
static int L1_NMTRA_FAST_SEC_RECO   = 0;  // Reference secondary 			( RefSec )
static int L1_NMTRA_SLOW_RECO       = 0;  // Extra     (extra set of tracks  -  0.1<momentum<=1.0)
static int L1_NMTRA_SLOW_PRIM_RECO  = 0;  // Extra primary 				( ExtraPrim )
static int L1_NMTRA_SLOW_SEC_RECO   = 0;  // Extra secondary				( ExtraSec )
static int L1_NTRACKS               = 0;  // finding track (vRTrack)
static int L1_NCLONES               = 0;  // Clone track
static int L1_NGHOSTS               = 0;  // Ghost track

					// number of killed  (see reco)
static int L1_NMTRA_D0_KILLED         = 0;
static int L1_NMTRA_TOTAL_KILLED      = 0;
static int L1_NMTRA_FAST_KILLED       = 0;
static int L1_NMTRA_FAST_PRIM_KILLED  = 0;
static int L1_NMTRA_FAST_SEC_KILLED   = 0;
static int L1_NMTRA_SLOW_KILLED       = 0;
static int L1_NMTRA_SLOW_PRIM_KILLED  = 0;
static int L1_NMTRA_SLOW_SEC_KILLED   = 0;

static int L1_NHITS[20], L1_NFAKES[20];  // number of hits(fakes) on station

static int L1_NEVENTS               = 0; 		// number of Performanced events
static double L1_CATIME             = 0.0;  // time spent for performance


void Performance::ExecPerf()
{

	ProfileData elem; // will use for full profiles vectors
	
	
  typedef  map<CbmL1MCTrack*, CbmL1Track*>  Mtrartramap;  // Monte Carlo tracks(MCtrack) to finded tracks map
  typedef  map<CbmL1MCTrack*, CbmL1Track*>::iterator MtrartramapIt;
  Mtrartramap   mtrartramap, mtrartradistortmap; // mtrartramap - for reconstracted track, // mtrartradistortmap - for distorted(ghost) track
  typedef  map<CbmL1Track*, CbmL1MCTrack*>           Rtramtramap;  // finded track to MC_track map
  typedef  map<CbmL1Track*, CbmL1MCTrack*>::iterator RtramtramapIt;
  Rtramtramap   rtramtramap, rtramtradistortmap; // rtramtramap - for reconstracted(unghost) track, // rtramtradistortmap - for ghost track
	
	static bool first_call = 1;
  if ( first_call )
    {
      first_call = 0;
     
      for (int i = 0; i < algo->NStations; i++){        // go through all algorithm stations
				L1_NHITS[i]  = 0;
				L1_NFAKES[i] = 0;
      }
    }// first_call 
  
	  //----------      TRACK FINDING PART         -----------------------

  map <int, int> recostamap, mcstamap;  //  DON USES!

  map <CbmL1MCTrack*, int> pmtramap; 		// pmtramap 	- number of finding track(Rtrack) used hits from definite Mtrack | reco track + clones    | for unghost Rtrack
  map <CbmL1Track*, int> rtradistortmap; 	// rtradistortmap - analog of pmtramap for ghost Rtrack
  
	// pmctrackmap - map of pointers on MCtrack  indexed by ID
  map <int, CbmL1MCTrack*> pmctrackmap; 	
  pmctrackmap.clear(); 
		// filling of pmctrackmap by pointers on vMCTracks 
  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{  // track repeated in array vMCTracks
      cout << "*** L1: Track ID " << MC.ID << " encountered twice! ***" << endl;
    }
  }

  		//  ---    track finding performance  ---

  pmtramap.clear();
  rtradistortmap.clear();
    
  mtrartramap.clear();
  rtramtramap.clear();
    
  mtrartradistortmap.clear();
  rtramtradistortmap.clear();
    
  bool vHitUsed[vHitStore.size()];      // vHitUsed - for mark all used in Rtracks hits from vHitStore
  for( unsigned int ih = 0; ih < vHitStore.size(); ih++) vHitUsed[ih]=0;

  int ntracks = 0;  // number of considered(from vRTracks)  tracks 
  int nghosts = 0;	// number of considered(from vRTracks) ghost tracks 


  	//  --- reconstructed tracks --- 
		
  for (vector<CbmL1Track>::iterator 	rtraIt = vRTracks.begin(); 
  																		rtraIt != vRTracks.end(); ++rtraIt){ // rtraIt - interator for vRTracks
      
    ntracks++;
    
    int hitsum  = 0;  // will count number of hit in current Rtrack
    
    typedef  map<int, int >           Hitmap;   // Hitmap<trackID, nhits>
    typedef  map<int, int >::iterator HitmapIt;
    Hitmap   hitmap;  // hitmap - number of hits belong each MCtrack(with define trackID)
    hitmap.clear();
     	// add hits from current Rtrack to hitmap,  and fill vHitUsed
    for (vector<int>::iterator 	ih  = (rtraIt->StsHits).begin();
	 															ih != (rtraIt->StsHits).end(); ++ih){
      //L1StsHit &h = algo->vStsHits[*ih];
      vHitUsed[*ih]=1;  // if hit belong some Rtrack - mark it
      int ID = -1;
      int iMC = vHitMCRef[*ih];  // get index of MC point
      if (iMC >= 0)  // all hits must have corespondend point, which marked in vHitMCRef. else ID = -1
				ID = vMCPoints[iMC].ID;
			// registred hit in hitmap
      if(hitmap.find(ID) == hitmap.end())  
				hitmap.insert(pair<int, int>(ID, 1));  
      else{
				hitmap[ID] += 1;
      }
      
      hitsum++;
    }

// iklm		
// for (unsigned	ih2  = 0;	ih2 < algo->vStsHits.size(); ++ih2){		
// 	int iMC = vHitMCRef[ih2];
// 	if (iMC == -1) continue;
// 	int ID = vMCPoints[iMC].ID;
// 	if (ID == 96) cout << ih2 << " ";
// };
// iklm
// cout << hitsum << endl;   
   
    // RTRA(real(finded) track?) <-> MTRA(MC track?) identification 
    int ghost = 1;
    double momentum = 0.0;	// momentum of particle. NOT USED
    double max_percent = 0.0; 	// max._number_of_hits_with_one_ID_in_hitmap/number_of_all_hits*100%. Only for histo??
    for( HitmapIt 	posIt = hitmap.begin(); 
    								posIt != hitmap.end(); posIt++ ){ // for Hitmap
      
      if (posIt->first < 0) continue; //if there are hits with ID = -1 it is fault
	// find the ID of most hits in Rtrack
				// count max_percent
      if (100.0*double(posIt->second) > max_percent*double(hitsum)) 
			max_percent = 100.0*double(posIt->second)/double(hitsum);
      
      if ( double(posIt->second) > 0.7*double(hitsum) ){ //MTRA reconstructed       // if more then 70% hits from current Rtrack have same ID(belong to one MCtrack) it is count as reconstructed(not ghost)!!!!
// iklm
// cout << posIt->first << " " << hitsum << endl;				
				// link current Rtrack to MCtrack with ID correspondent to most of hits in current Rtrack, and vice versa  in maps
				if (pmctrackmap.find(posIt->first) == pmctrackmap.end()) continue;
				CbmL1MCTrack* pmtra = pmctrackmap[posIt->first]; // pointer to MC track
				CbmL1Track* prtra = &(*rtraIt);
				mtrartramap.insert(pair<CbmL1MCTrack*, CbmL1Track*>(pmtra, prtra));
				rtramtramap.insert(pair<CbmL1Track*, CbmL1MCTrack*>(prtra, pmtra));
				
				momentum = pmtra->p;  // get momentum from correspondent to current Rtrack MCtrack
					
					// registred finding correspondet to define MCtrack(=pmtra) Rtrack(=current track) | reconstracted(unghost)
				if(pmtramap.find(pmtra) == pmtramap.end())
					pmtramap.insert(pair<CbmL1MCTrack*, int>(pmtra, 1));
				else{
					pmtramap[pmtra]++;
				}
				
				ghost = 0; 
			}
			else{ // MTRA destorted  - ghost Rtrack
      		// link current Rtrack to MCtrack with ID correspondent to current hit, and vice versa in maps(I don't see the reason!??)
				if (pmctrackmap.find(posIt->first) == pmctrackmap.end()) continue;
				CbmL1MCTrack* pmtra = pmctrackmap[posIt->first]; //pointer to MC track
				CbmL1Track* prtra = &(*rtraIt);
				mtrartradistortmap.insert(pair<CbmL1MCTrack*, CbmL1Track*>(pmtra, prtra));
				rtramtradistortmap.insert(pair<CbmL1Track*, CbmL1MCTrack*>(prtra, pmtra));
				
					// registred finding correspondet to define MCtrack(=pmtra) Rtrack(=current track) | ghost
				if(rtradistortmap.find(prtra) == rtradistortmap.end())
					rtradistortmap.insert(pair<CbmL1Track*, int>(prtra, 1));
				else{
					rtradistortmap[prtra]++;
				}
      }
    }
    
    nghosts += ghost;

  }//reco tracks

  // clone 
  /*
  for (vector<CbmL1Track>::iterator rtraIt = vRTracks.begin(); rtraIt != vRTracks.end(); ++rtraIt){
         	
    CbmL1Track* prtra = &(*rtraIt);
    
    if( rtramtramap.find(prtra) == rtramtramap.end() ) continue;
    CbmL1MCTrack *pmtra = rtramtramap[prtra];
    
    if(pmtramap.find(pmtra) == pmtramap.end()) continue;
    if( pmtramap[pmtra]<=1 ) continue;
    CbmL1HitStore &fh = vHitStore[prtra->StsHits[0]];
    CbmL1HitStore &lh = vHitStore[prtra->StsHits[prtra->StsHits.size()-1]];
    cout<<"clone "<<pmtra->ID<<" "<<fh.iStation<<" "<<lh.iStation<<endl;	
  }
  */
	
		 
	if(0){	 // ---- store data for histogram  ---
  
	 HistoData e;		 
	 for (unsigned i = 0; i < vRTracks.size(); i++){
		  CbmL1Track *rtr = &vRTracks[i];
		  CbmL1MCTrack *mctr = rtramtramap[rtr]; // get MCTrack corresponding to recoTrack
		  if (!mctr) continue;
		  
		  // find the MCpoint, witch correspond first hit of recoTrack
		  double z = rtr->T[5];
		  const double dmaxz = 0.1;
		  unsigned ip = 0;
  
		  for (; ip < mctr->Points.size(); ip++){
			 if ((mctr->Points[ip]->z <= z+dmaxz) && (mctr->Points[ip]->z >= z-dmaxz)) break;
		  }
		  if (ip == mctr->Points.size()){
			 //cout << " Performance - warning - dont find the mcPoint z=" << z << endl;
			 //cout << " mctr->Points.size() " << mctr->Points.size() << " mctr->ID " << mctr->ID << endl;	
			 continue;
		  };
  // 		cout << " dz " << mctr->Points[ip]->z - z << " MCz " << mctr->Points[ip]->z << endl ;
	    double mc[5], c[5];
		  double pz = mctr->Points[ip]->pz;
      if( pz <1.e-4 ) continue;
		  mc[0]=mctr->Points[ip]->x;
		  mc[1]=mctr->Points[ip]->y;
		  mc[2]=mctr->Points[ip]->px/pz;
		  mc[3]=mctr->Points[ip]->py/pz;
		  mc[4]= mctr->Points[ip]->q/mctr->Points[ip]->p;
	    c[0]=	rtr->C[0];
	    c[1]=	rtr->C[2];
	    c[2]=	rtr->C[5];
	    c[3]=	rtr->C[9];
	    c[4]=	rtr->C[14];
	   
		  bool ok = 1;
		  for( int i=0; i<5; i++ ){
			 ok = ok && c[i]>=0 && finite(c[i]) && finite(rtr->T[i]);						
		  }
		  if( !ok ) continue;
			 
		  
		  
		  if (fabs(rtr->T[4]) < 1.0/0.5)
		  for( int i=0; i<5; i++ ){		
		    //if (i == 4) e.x = rtr->T[4]/mc[4]-1;
		    e.x = rtr->T[i] - mc[i];
		    hist[i][0].push_back(e);
		    e.x = sqrt(c[i]);
		    hist[i][1].push_back(e);
		    e.x = (rtr->T[i] - mc[i])/sqrt(c[i]);
		    //if (i == 2) e.x = mctr->Points[ip]->q*(rtr->T[i] - mc[i])/sqrt(c[i]);
		    hist[i][2].push_back(e);
		  }		
	 } // for vRTracks rtramtramap
	}
		 
	
			//  ---- calculation of efficiency  ----
			
	// local counters. after will be added to L1_NMTRA_D0 and else correspondently (see begin of file)
  int nmtra_d0         = 0;
  int nmtra_total      = 0;
  int nmtra_fast       = 0;
  int nmtra_fast_prim  = 0;
  int nmtra_fast_sec   = 0;
  int nmtra_slow       = 0;
  int nmtra_slow_prim  = 0;
  int nmtra_slow_sec   = 0;
  
  int nmtra_d0_reco    = 0;
  int nmtra_total_reco = 0;
  int nmtra_fast_reco  = 0;
  int nmtra_fast_prim_reco  = 0;
  int nmtra_fast_sec_reco   = 0;
  int nmtra_slow_reco  = 0;
  int nmtra_slow_prim_reco  = 0;
  int nmtra_slow_sec_reco   = 0;
  int nclones          = 0;

  int nmtra_d0_killed    = 0;
  int nmtra_total_killed = 0;
  int nmtra_fast_killed  = 0;
  int nmtra_fast_prim_killed  = 0;
  int nmtra_fast_sec_killed   = 0;
  int nmtra_slow_killed  = 0;
  int nmtra_slow_prim_killed  = 0;
  int nmtra_slow_sec_killed   = 0;



  int sta_nhits[algo->NStations], sta_nfakes[algo->NStations];  // sta_nhits - number of hits on each station // sta_nfakes - number of fakes on each station ?? 
  for (int i = 0; i < algo->NStations; i++){
    sta_nhits[i]  = 0;
    sta_nfakes[i] = 0;
  }

  for ( map<int, CbmL1MCTrack*>::iterator 	mtraIt = pmctrackmap.begin();
						mtraIt != pmctrackmap.end(); mtraIt++ ) {   // mtraIt  - loop for MCtrack
      
    // // No Sts hits? 
    int nmchits = (mtraIt->second)->StsHits.size();  // nmchits - count of hits. For fillter short track.

  
    double momentum = (mtraIt->second)->p;  // momentum - momentum of particle. Need for define if track is fast
    
    if (nmchits < 2) continue; // allset  3 or 4 or 5 stations ?  !!!!
    
    int ista = -1; 	// index of station
    
    	// counts the stations
    int sta_count = 0; 	// count the hits(=stations)
    for( vector<int>::iterator 	ih = (mtraIt->second)->StsHits.begin();
	 															ih!= (mtraIt->second)->StsHits.end(); ++ih){  // loop for hits in MCtrack
      CbmL1HitStore &mh = vHitStore[*ih];
      if (mh.iStation <= ista ) continue;
      sta_count++;
    }
/*
    	// fill the gisto. for registred
    CbmL1HitStore &fh = vHitStore[*((mtraIt->second)->StsHits.begin())];
    CbmL1HitStore &lh = vHitStore[*((mtraIt->second)->StsHits.rbegin())];
    if ((mtraIt->second)->mother_ID < 0){ 	// primary
      h_reg_mom_prim->Fill(momentum);
      h_reg_nhits_prim->Fill(sta_count);
      h2_reg_nhits_vs_mom_prim->Fill(momentum, sta_count);
      h2_reg_fstation_vs_mom_prim->Fill(momentum, fh.iStation+1);
      if (lh.iStation >= fh.iStation)
	h2_reg_lstation_vs_fstation_prim->Fill(fh.iStation+1, lh.iStation+1);
    }else{ 					// secondary
      h_reg_mom_sec->Fill(momentum);
      h_reg_nhits_sec->Fill(sta_count);
      if (momentum > 0.01){
	h2_reg_nhits_vs_mom_sec->Fill(momentum, sta_count);
	h2_reg_fstation_vs_mom_sec->Fill(momentum, fh.iStation+1);
	if (lh.iStation >= fh.iStation)
	  h2_reg_lstation_vs_fstation_sec->Fill(fh.iStation+1, lh.iStation+1);
      }
    }*/
      

    //cout << " => " << sta_count << endl;
    //if (sta_count < 4) continue;      //  detected in less than 4 stations !!!!
		
    if( mtraIt->second->nMCContStations < 4 ) continue;
      
    if (momentum < 0.1) continue;	// reject very slow tracks from analysis !!!!

    bool reco = (pmtramap.find(mtraIt->second) != pmtramap.end());  // reconstracted = not ghost!!!!

			// fill histo data
		elem.x = momentum;
		if (reco){
			elem.y = 100.0;
		}else{
			elem.y = 0.0;
		}
 		p_eff_all_vs_mom.push_back(elem);
		if ((mtraIt->second)->mother_ID < 0){ // primary
			p_eff_prim_vs_mom.push_back(elem);
		}else{
			p_eff_sec_vs_mom.push_back(elem);
		}	
		
 	// count the number of clones
    if (pmtramap.find(mtraIt->second) != pmtramap.end()){
      nclones += pmtramap[mtraIt->second]-1;
    }
 
    bool killed = 0;  // killed - same unreco track will be killed      
    if(!reco){
      int nusedhits = 0; // number of used from current MCtrack hits
      for( vector<int>::iterator 	ih = (mtraIt->second)->StsHits.begin();
	   															ih != (mtraIt->second)->StsHits.end(); ih++ ){
				if( vHitUsed[*ih] ) nusedhits++;
      }
      if( nusedhits > 0 ) killed = 1;      // if some hits from current(non reconstracted!) MCtrack used in some Rtrack then it was killed
    }
    
			// fill histo
		elem.x = momentum;
    if (((mtraIt->second)->mother_ID < 0)&&((mtraIt->second)->z > 0)){ // D0
      nmtra_d0++;
      if (reco) nmtra_d0_reco++;
			if (reco) elem.y = 100.0;
      else      elem.y = 0.0;
			
      if (killed) nmtra_d0_killed++;
    }
		p_eff_d0_vs_mom.push_back(elem);
    
		
		nmtra_total++;
    if (reco) nmtra_total_reco++;
    if (killed) nmtra_total_killed++;

    if (momentum>1.0 /*&&  mtraIt->second->nMCContStations >=6*/){ // reference tracks !!!!
      nmtra_fast++;
      if (reco) nmtra_fast_reco++;
      if (killed) nmtra_fast_killed++;
	
      if ((mtraIt->second)->mother_ID < 0){ 	// primary
				nmtra_fast_prim++;
				if (reco) nmtra_fast_prim_reco++;
				if (killed) nmtra_fast_prim_killed++;
     }
      else{					// secondary
				nmtra_fast_sec++;
				if (reco) nmtra_fast_sec_reco++;
				if (killed) nmtra_fast_sec_killed++;
      }
    }
    else{   // extra set of tracks  -  0.1<momentum<=1.0
      nmtra_slow++;
      if (reco) nmtra_slow_reco++;      
      if (killed) nmtra_slow_killed++;
       
      if ((mtraIt->second)->mother_ID < 0){ 	// primary
				nmtra_slow_prim++;
				if (reco) nmtra_slow_prim_reco++;
				if (killed) nmtra_slow_prim_killed++;
      }
      else{					// secondary
				nmtra_slow_sec++;
				if (reco) nmtra_slow_sec_reco++;
				if (killed) nmtra_slow_sec_killed++;
      }
    } // else
  } // for ( map<int, CbmL1MCTrack*>::iterator 	mtraIt

  L1_CATIME += algo->CATime;

  	// --- count a Fakes and hits on each stations ---  only for histo
  int NFakes = 0;
  for( unsigned ih=0; ih < algo->vStsHits.size(); ih++){
    int iMC = vHitMCRef[ih];
    if (iMC < 0) NFakes++;  // if hit not have corresponding MCpoint it is fake
    //station
    int ista = vHitStore[ih].iStation;
    sta_nhits[ista]++;
    L1_NHITS[ista]++;
    if (iMC < 0){   // hit is Fake
      sta_nfakes[ista]++;
      L1_NFAKES[ista]++;
    }
  }

  L1_NMTRA_D0              += nmtra_d0;
  L1_NMTRA_TOTAL           += nmtra_total;
  L1_NMTRA_FAST            += nmtra_fast;
  L1_NMTRA_FAST_PRIM       += nmtra_fast_prim;
  L1_NMTRA_FAST_SEC        += nmtra_fast_sec;
  L1_NMTRA_SLOW            += nmtra_slow;
  L1_NMTRA_SLOW_PRIM       += nmtra_slow_prim;
  L1_NMTRA_SLOW_SEC        += nmtra_slow_sec;
  
  L1_NMTRA_D0_RECO         += nmtra_d0_reco;
  L1_NMTRA_TOTAL_RECO      += nmtra_total_reco;
  L1_NMTRA_FAST_RECO       += nmtra_fast_reco;
  L1_NMTRA_FAST_PRIM_RECO  += nmtra_fast_prim_reco;
  L1_NMTRA_FAST_SEC_RECO   += nmtra_fast_sec_reco;
  L1_NMTRA_SLOW_RECO       += nmtra_slow_reco;
  L1_NMTRA_SLOW_PRIM_RECO  += nmtra_slow_prim_reco;
  L1_NMTRA_SLOW_SEC_RECO   += nmtra_slow_sec_reco;
  L1_NCLONES               += nclones;
  L1_NGHOSTS               += nghosts;
  L1_NTRACKS               += ntracks;

  L1_NMTRA_D0_KILLED         += nmtra_d0_killed;
  L1_NMTRA_TOTAL_KILLED      += nmtra_total_killed;
  L1_NMTRA_FAST_KILLED       += nmtra_fast_killed;
  L1_NMTRA_FAST_PRIM_KILLED  += nmtra_fast_prim_killed;
  L1_NMTRA_FAST_SEC_KILLED   += nmtra_fast_sec_killed;
  L1_NMTRA_SLOW_KILLED       += nmtra_slow_killed;
  L1_NMTRA_SLOW_PRIM_KILLED  += nmtra_slow_prim_killed;
  L1_NMTRA_SLOW_SEC_KILLED   += nmtra_slow_sec_killed;
  
  L1_NEVENTS++;
  
  double p_total=0.0, p_fast=0.0, p_fast_prim=0.0, p_fast_sec=0.0, p_slow=0.0, p_slow_prim=0.0, p_slow_sec=0.0;
  double p_clone=0.0, p_ghost=0.0;
  if (nmtra_total    !=0) p_total = double(nmtra_total_reco)/double(nmtra_total);
  if (nmtra_fast     !=0) p_fast  = double(nmtra_fast_reco)/double(nmtra_fast);
  if (nmtra_fast_prim!=0) p_fast_prim  = double(nmtra_fast_prim_reco)/double(nmtra_fast_prim);
  if (nmtra_fast_sec !=0) p_fast_sec  = double(nmtra_fast_sec_reco)/double(nmtra_fast_sec);
  if (nmtra_slow     !=0) p_slow  = double(nmtra_slow_reco)/double(nmtra_slow);
  if (nmtra_slow_prim!=0) p_slow_prim  = double(nmtra_slow_prim_reco)/double(nmtra_slow_prim);
  if (nmtra_slow_sec !=0) p_slow_sec  = double(nmtra_slow_sec_reco)/double(nmtra_slow_sec);
  if (ntracks        !=0) p_clone = double(nclones)/double(ntracks);
  if (ntracks        !=0) p_ghost = double(nghosts)/double(ntracks);
  
  double P_D0=0.0, P_TOTAL=0.0, P_FAST=0.0, P_FAST_PRIM=0.0, P_FAST_SEC=0.0, P_SLOW=0.0, P_SLOW_PRIM=0.0, P_SLOW_SEC=0.0;  // ratio between reconstracted number of track and initial //e.g. P_FAST_SEC = L1_NMTRA_FAST_SEC_RECO/L1_NMTRA_FAST_SEC
  double K_D0=0.0, K_TOTAL=0.0, K_FAST=0.0, K_FAST_PRIM=0.0, K_FAST_SEC=0.0, K_SLOW=0.0, K_SLOW_PRIM=0.0, K_SLOW_SEC=0.0;  // ratio between killed number of track and initial //e.g.  K_FAST_SEC = L1_NMTRA_FAST_SEC_KILLED/L1_NMTRA_FAST_SEC
  double P_CLONE=0.0, P_GHOST=0.0; // ratio between clone(ghost) number of track and reconstracted //e.g. _GHOST = L1_NGHOSTS/L1_NMTRA_TOTAL_RECO
  if (L1_NMTRA_D0!=0){
    P_D0 = double(L1_NMTRA_D0_RECO)/double(L1_NMTRA_D0);
    K_D0 = double(L1_NMTRA_D0_KILLED)/double(L1_NMTRA_D0);
  }
  if (L1_NMTRA_TOTAL!=0){
    P_TOTAL = double(L1_NMTRA_TOTAL_RECO)/double(L1_NMTRA_TOTAL);
    K_TOTAL = double(L1_NMTRA_TOTAL_KILLED)/double(L1_NMTRA_TOTAL);
  }
  if (L1_NMTRA_FAST!=0){
    P_FAST  = double(L1_NMTRA_FAST_RECO)/double(L1_NMTRA_FAST);
    K_FAST  = double(L1_NMTRA_FAST_KILLED)/double(L1_NMTRA_FAST);
  }
  if (L1_NMTRA_FAST_PRIM!=0){
    P_FAST_PRIM  = double(L1_NMTRA_FAST_PRIM_RECO)/double(L1_NMTRA_FAST_PRIM);
    K_FAST_PRIM  = double(L1_NMTRA_FAST_PRIM_KILLED)/double(L1_NMTRA_FAST_PRIM);
  }
  if (L1_NMTRA_FAST_SEC!=0){
    P_FAST_SEC  = double(L1_NMTRA_FAST_SEC_RECO)/double(L1_NMTRA_FAST_SEC);
    K_FAST_SEC  = double(L1_NMTRA_FAST_SEC_KILLED)/double(L1_NMTRA_FAST_SEC);
  }
  if (L1_NMTRA_SLOW!=0){
    P_SLOW  = double(L1_NMTRA_SLOW_RECO)/double(L1_NMTRA_SLOW);
    K_SLOW  = double(L1_NMTRA_SLOW_KILLED)/double(L1_NMTRA_SLOW);
  }
  if (L1_NMTRA_SLOW_PRIM!=0){
    P_SLOW_PRIM  = double(L1_NMTRA_SLOW_PRIM_RECO)/double(L1_NMTRA_SLOW_PRIM);
    K_SLOW_PRIM  = double(L1_NMTRA_SLOW_PRIM_KILLED)/double(L1_NMTRA_SLOW_PRIM);
  }
  if (L1_NMTRA_SLOW_SEC!=0){
    P_SLOW_SEC  = double(L1_NMTRA_SLOW_SEC_RECO)/double(L1_NMTRA_SLOW_SEC);
    K_SLOW_SEC  = double(L1_NMTRA_SLOW_SEC_KILLED)/double(L1_NMTRA_SLOW_SEC);
  }
  if (L1_NMTRA_TOTAL_RECO!=0)
    //P_CLONE = double(L1_NCLONES)/double(L1_NTRACKS);
    P_CLONE = double(L1_NCLONES)/double(L1_NMTRA_TOTAL);
  if (L1_NMTRA_TOTAL_RECO!=0)
    //P_GHOST = double(L1_NGHOSTS)/double(L1_NTRACKS);
    P_GHOST = double(L1_NGHOSTS)/double(L1_NMTRA_TOTAL_RECO);
    
  
  if( iVerbose){
    cout.setf(ios::fixed);
    cout.setf(ios::showpoint);
    cout.precision(3);
#ifdef XXX
  cout                                                   << endl;
  cout << "L1 PER EVENT STAT      : "                    << endl;
  cout << "MC Refset              : " << nmtra_fast      << endl;
  cout << "MC Extras              : " << nmtra_slow      << endl;
  cout << "ALL SIMULATED          : " << nmtra_total     << endl;
  cout                                                   << endl;
  cout << "RC Refset              : " << nmtra_fast_reco << endl;
  cout << "RC Extras              : " << nmtra_slow_reco << endl;
  cout << "ghosts                 : " << nghosts         << endl;
  cout << "clones                 : " << nclones         << endl;
  cout << "ALL RECONSTRUCTED      : " << ntracks         << endl;
  cout                                                   << endl;
  cout << "RefPrim efficiency     : " << p_fast_prim     << endl;
  cout << "RefSec efficiency      : " << p_fast_sec      << endl;
  cout << "Refset efficiency      : " << p_fast          << endl;
  cout << "Allset efficiency      : " << p_total         << endl;
  cout << "Extra  efficiency      : " << p_slow          << endl;
  cout << "clone  probability     : " << p_clone         << endl;
  cout << "ghost  probability     : " << p_ghost         << endl;
  cout                                                   << endl;

  	
  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;

#endif//XXX
  cout                                                   << endl;
  cout << "L1 ACCUMULATED STAT    : " << L1_NEVENTS << " EVENTS "               << endl << endl;
  
  if (L1_NMTRA_D0!=0)
    cout << "D0        efficiency   : " << P_D0        <<" | "<<L1_NMTRA_D0_RECO   << endl;

// 	  cout << "RefPrim   efficiency   : " << L1_NMTRA_FAST_PRIM_RECO <<" / "<< L1_NMTRA_FAST_PRIM <<" = "<<P_FAST_PRIM   << endl;
	
  cout << "RefPrim   efficiency   : " << P_FAST_PRIM <<" / "<< K_FAST_PRIM <<" | "<<L1_NMTRA_FAST_PRIM_RECO   << endl;
  cout << "RefSec    efficiency   : " << P_FAST_SEC  <<" / "<< K_FAST_SEC  <<" | "<<L1_NMTRA_FAST_SEC_RECO    << endl;
  cout << "Refset    efficiency   : " << P_FAST      <<" / "<< K_FAST      <<" | "<<L1_NMTRA_FAST_RECO        << endl;
  cout << "Allset    efficiency   : " << P_TOTAL     <<" / "<< K_TOTAL     <<" | "<<L1_NMTRA_TOTAL_RECO       << endl;
  cout << "ExtraPrim efficiency   : " << P_SLOW_PRIM <<" / "<< K_SLOW_PRIM <<" | "<<L1_NMTRA_SLOW_PRIM_RECO   << endl;
  cout << "ExtraSec  efficiency   : " << P_SLOW_SEC  <<" / "<< K_SLOW_SEC  <<" | "<<L1_NMTRA_SLOW_SEC_RECO    << endl;
  cout << "Extra     efficiency   : " << P_SLOW      <<" / "<< K_SLOW      <<" | "<<L1_NMTRA_SLOW_RECO        << endl;
  cout << "Clone     probability  : " << P_CLONE     <<" | "<<L1_NCLONES                << endl;
  cout << "Ghost     probability  : " << P_GHOST     <<" | "<<L1_NGHOSTS                << endl;
  cout << "MC tracks/event found  : " << int(double(L1_NMTRA_TOTAL_RECO)/double(L1_NEVENTS)) << endl;
  
  cout<<endl;

  cout<<"CA Track Finder: " << L1_CATIME/L1_NEVENTS << " s/ev" << endl << endl;
  /*
  cout << "Mean number of all (true+fake) hits in stations: " << endl;
  for (int i = 0; i < algo->NStations; i++){
    cout << (int)((L1_NHITS[i]-L1_NFAKES[i])/L1_NEVENTS) + (int)(L1_NFAKES[i]/L1_NEVENTS)  << "(" << (int)((L1_NHITS[i]-L1_NFAKES[i])/L1_NEVENTS) << "+" << (int)(L1_NFAKES[i]/L1_NEVENTS) << ")   ";
  }
  cout << endl;

  cout.precision(1);
  cout << "Mean percentage (%) of true+fake hits in stations: " << endl;
  for (int i = 0; i < algo->NStations; i++){
    cout << 100.0-100.0*L1_NFAKES[i]/L1_NHITS[i] << "+" << 100.0*L1_NFAKES[i]/L1_NHITS[i] << "   ";
  }
  cout << endl;
  */
  cout.precision(3);
  
  //#endif//XXX
  }
}