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

#include "CbmL1.h"
#include "CbmKF.h"
#include "CbmKFMath.h"
#include "L1Algo/L1Algo.h"
#include "L1Algo/L1StsHit.h"

#include "iostream.h"


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;

static int L1_NMTRA_D0_RECO         = 0;
static int L1_NMTRA_TOTAL_RECO      = 0;
static int L1_NMTRA_FAST_RECO       = 0;
static int L1_NMTRA_FAST_PRIM_RECO  = 0;
static int L1_NMTRA_FAST_SEC_RECO   = 0;
static int L1_NMTRA_SLOW_RECO       = 0;
static int L1_NMTRA_SLOW_PRIM_RECO  = 0;
static int L1_NMTRA_SLOW_SEC_RECO   = 0;
static int L1_NTRACKS               = 0;
static int L1_NCLONES               = 0;
static int L1_NGHOSTS               = 0;

static int L1_NEVENTS               = 0;

static map<CbmL1Track*, CbmL1MCTrack* > RTRAtoMTRAmap;
static int RTRAtoMTRAmap_clear = 1;


void CbmL1::Performance()
{
  typedef  map<CbmL1MCTrack*, CbmL1Track*>  Mtrartramap;
  typedef  map<CbmL1MCTrack*, CbmL1Track*>::iterator MtrartramapIt;
  Mtrartramap   mtrartramap;
  typedef  map<CbmL1Track*, CbmL1MCTrack*>           Rtramtramap;
  typedef  map<CbmL1Track*, CbmL1MCTrack*>::iterator RtramtramapIt;
  Rtramtramap   rtramtramap;

  //CbmKF &KF = *CbmKF::Instance();
  // histogramms 

  static TProfile *p_eff_all_vs_mom, *p_eff_prim_vs_mom, *p_eff_sec_vs_mom, *p_eff_d0_vs_mom; 
  
  static TH1F *h_mom_track;
  static TH1F 
    *h_ghost_mom, *h_ghost_nhits, *h_ghost_station, 
    *h_ghost_chi2, *h_ghost_prob, *h_ghost_tx, *h_ghost_ty;
  static TH1F *h_reco_mom,  *h_reco_d0_mom, *h_reco_nhits, *h_reco_station, 
    *h_reco_chi2, *h_reco_prob, *h_reco_clean, *h_reco_time;
  static TProfile *h_reco_timeNtr, *h_reco_timeNhit ;
  static TProfile *h_reco_fakeNtr, *h_reco_fakeNhit ;
  static TH1F *h_tx, *h_ty, *h_sec_r, *h_ghost_r;

  static TH1F *h_notfound_mom, *h_notfound_nhits, *h_notfound_station, *h_notfound_r, *h_notfound_tx, *h_notfound_ty;
  
     
  static TH1F *MC_P, *nhits, *h_chi2, *h_prob, *MC_vx, *MC_vy, *MC_vz, *VtxChiPrim, *VtxChiSec;
    

  static TH2F *P_vs_P ;

  static bool first_call = 1;

  if ( first_call )
    {
      first_call = 0;

      //TDirectory *maindir = gDirectory;
      //histodir->cd();

      TDirectory *curdir = gDirectory;
      //histodir = gROOT->mkdir("L1");
      histodir->cd();

 
      p_eff_all_vs_mom = new TProfile("p_eff_all_vs_mom", "AllSet Efficiency vs Momentum",
				     100, 0.0, 5.0, 0.0, 100.0);
      p_eff_prim_vs_mom = new TProfile("p_eff_prim_vs_mom", "All Primary Set Efficiency vs Momentum",
				     100, 0.0, 5.0, 0.0, 100.0);
      p_eff_sec_vs_mom = new TProfile("p_eff_sec_vs_mom", "All Secondary Set Efficiency vs Momentum",
				     100, 0.0, 5.0, 0.0, 100.0);
      p_eff_d0_vs_mom = new TProfile("p_eff_d0_vs_mom", "D0 Secondary Tracks Efficiency vs Momentum",
				     150, 0.0, 15.0, 0.0, 100.0);

      h_mom_track   = new TH1F("h_mom_track", "Momentum of accepted tracks", 100, 0.0, 20.0);

      h_ghost_mom = new TH1F("h_ghost_mom", "Momentum of ghost track", 50, 0.0, 5.0);
      h_ghost_nhits = new TH1F("h_ghost_nhits", "Number of hits in ghost track", 50, 0.0, 10.0);
      h_ghost_station = new TH1F("h_ghost_station", "First station of ghost track", 50, -0.5, 10.0);
      h_ghost_chi2 = new TH1F("h_ghost_chi2", "Chi2/NDF of ghost track", 50, -0.5, 10.0);
      h_ghost_prob = new TH1F("h_ghost_prob", "Prob of ghost track", 50, 0., 1.01);
      h_ghost_r = new TH1F("h_ghost_r", "R of ghost track at the first hit", 50, 0.0, 15.0);
      h_ghost_tx = new TH1F("h_ghost_tx", "tx of ghost track at the first hit", 50, -5.0, 5.0);
      h_ghost_ty = new TH1F("h_ghost_ty", "ty of ghost track at the first hit", 50, -1.0, 1.0);

      h_reco_mom = new TH1F("h_reco_mom", "Momentum of reco track", 50, 0.0, 5.0);
      h_reco_nhits = new TH1F("h_reco_nhits", "Number of hits in reco track", 50, 0.0, 10.0);
      h_reco_station = new TH1F("h_reco_station", "First station of reco track", 50, -0.5, 10.0);
      h_reco_chi2 = new TH1F("h_reco_chi2", "Chi2/NDF of reco track", 50, -0.5, 10.0);
      h_reco_prob = new TH1F("h_reco_prob", "Prob of reco track", 50, 0., 1.01);
      h_reco_clean = new TH1F("h_reco_clean", "Percentage of correct hits", 100, -0.5, 100.5);
      h_reco_time = new TH1F("h_reco_time", "CA Track Finder Time (s/ev)", 20, 0.0, 20.0);
      h_reco_timeNtr = new TProfile("h_reco_timeNtr", 
				    "CA Track Finder Time (s/ev) vs N Tracks", 
				    200, 0.0, 1000.0);
      h_reco_timeNhit = new TProfile("h_reco_timeNhit", 
				    "CA Track Finder Time (s/ev) vs N Hits", 
				    200, 0.0, 30000.0);
      h_reco_fakeNtr = new TProfile("h_reco_fakeNtr", 
				    "N Fake hits vs N Tracks", 
				    200, 0.0, 1000.0);
      h_reco_fakeNhit = new TProfile("h_reco_fakeNhit", 
				    "N Fake hits vs N Hits", 
				    200, 0.0, 30000.0);

      h_reco_d0_mom = new TH1F("h_reco_d0_mom", "Momentum of reco D0 track", 150, 0.0, 15.0);

      h_tx = new TH1F("h_tx", "tx of MC track at the vertex", 50, -0.5, 0.5);
      h_ty = new TH1F("h_ty", "ty of MC track at the vertex", 50, -0.5, 0.5);

      h_sec_r = new TH1F("h_sec_r", "R of sec MC track at the first hit", 50, 0.0, 15.0);

      h_notfound_mom = new TH1F("h_notfound_mom", "Momentum of not found track", 50, 0.0, 5.0);
      h_notfound_nhits = new TH1F("h_notfound_nhits", "Num hits of not found track", 50, 0.0, 10.0);
      h_notfound_station = new TH1F("h_notfound_station", "First station of not found track", 50, 0.0, 10.0);
      h_notfound_r = new TH1F("h_notfound_r", "R at first hit of not found track", 50, 0.0, 15.0);
      h_notfound_tx = new TH1F("h_notfound_tx", "tx of not found track at the first hit", 50, -5.0, 5.0);
      h_notfound_ty = new TH1F("h_notfound_ty", "ty of not found track at the first hit", 50, -1.0, 1.0);

      h_chi2 = new TH1F("chi2", "Chi^2", 100, 0.0, 10.);
      h_prob = new TH1F("prob", "Prob", 100, 0.0, 1.01);
      VtxChiPrim = new TH1F("vtxChiPrim", "Chi to primary vtx for primary tracks", 100, 0.0, 10.);
      VtxChiSec = new TH1F("vtxChiSec", "Chi to primary vtx for secondary tracks", 100, 0.0, 10.);
      MC_P = new TH1F("MC_P", "MC P/Q ", 1000, 0., 10.);
      MC_vx = new TH1F("MC_vx", "MC Vertex X", 100, -.05, .05);
      MC_vy = new TH1F("MC_vy", "MC Vertex Y", 100, -.05, .05);
      MC_vz = new TH1F("MC_vz", "MC Vertex Z", 100, -.05, .05);
      nhits = new TH1F("nhits", "N STS hits on MC track", 20, 0., 20.);
 
      P_vs_P = new TH2F("P_vs_P", "Resolution P/Q vs P", 20, 0., 20.,100, -.05, .05);
      
      //maindir->cd(); 

      // ----- Create list of all histogram pointers
      
      curdir->cd();      
      
    }// first_call 

  
  //track finding part 

  typedef  map<int, int >           Hitmap;   // <trackID, nhits>
  typedef  map<int, int >::iterator HitmapIt;
  Hitmap   hitmap;
    
  map <int, int> recostamap, mcstamap;

  map <CbmL1MCTrack*, int> pmtramap; // reco track + clones  

  map <int, CbmL1MCTrack*> pmctrackmap; pmctrackmap.clear();


  //track finding performance

  pmtramap.clear();
    
  mtrartramap.clear();
  rtramtramap.clear();
    
  if (RTRAtoMTRAmap_clear){
    RTRAtoMTRAmap.clear();
    RTRAtoMTRAmap_clear = 0;
  }
    
  int ntracks = 0;
  int nghosts = 0;
  //cout<<"reco.."<<endl;
  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;
    }
  }
  
  //reconstructed tracks
  for (vector<CbmL1Track>::iterator rtraIt = vRTracks.begin(); rtraIt != vRTracks.end(); ++rtraIt){
      
    ntracks++;
    
    int hitsum  = 0;
    hitmap.clear();
    
    for (vector<int>::iterator ih = (rtraIt->StsHits).begin();
	 ih != (rtraIt->StsHits).end(); ++ih){
      //L1StsHit &h = algo->vStsHits[*ih];
      int ID = -1;
      int iMC = vHitMCRef[*ih];
      if (iMC >= 0)
	ID = vMCPoints[iMC].ID;
      if(hitmap.find(ID) == hitmap.end())
	hitmap.insert(pair<int, int>(ID, 1));
      else{
	hitmap[ID] += 1;
      }
      hitsum++;
    }

    {      
      CbmL1Track* prtra = &(*rtraIt);
      if( fabs(prtra->T[4])>1.e-10 ) h_reco_mom->Fill(fabs(1.0/prtra->T[4]));
      h_reco_nhits->Fill((prtra->StsHits).size());
      CbmL1HitStore &mh = vHitStore[prtra->StsHits[0]];
      h_reco_station->Fill(mh.iStation);
      if (prtra->NDF > 0){
	h_reco_chi2->Fill(prtra->chi2/prtra->NDF);
	h_reco_prob->Fill(TMath::Erf(sqrt(fabs(prtra->chi2/prtra->NDF))));
      }
    }
    
    // RTRA <-> MTRA identification 
    int ghost = 1;
    double momentum = 0.0;
    double max_percent = 0.0;
    for( HitmapIt posIt = hitmap.begin(); posIt != hitmap.end(); posIt++ ){ // for Hitmap
      
      if (posIt->first < 0) continue; //not a MC track - based on fake hits

      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 (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;
	  
	if(pmtramap.find(pmtra) == pmtramap.end())
	  pmtramap.insert(pair<CbmL1MCTrack*, int>(pmtra, 1));
	else{
	  pmtramap[pmtra]++;
	}
	
	RTRAtoMTRAmap.insert(pair<CbmL1Track*, CbmL1MCTrack*>(prtra, pmtra));
	
	ghost = 0;

      }
    }

    h_reco_clean->Fill(max_percent);
    
    nghosts += ghost;
    
    if (ghost == 1){ //ghost
      CbmL1Track* prtra = &(*rtraIt);
      if( fabs(prtra->T[4])>1.e-10) h_ghost_mom->Fill(fabs(1.0/prtra->T[4]));
      h_ghost_nhits->Fill((prtra->StsHits).size());
      CbmL1HitStore &h1 = vHitStore[prtra->StsHits[0]];
      CbmL1HitStore &h2 = vHitStore[prtra->StsHits[1]];
      h_ghost_station->Fill(h1.iStation);
      if (prtra->NDF > 0){
	h_ghost_chi2->Fill(prtra->chi2/prtra->NDF);
	h_ghost_prob->Fill(TMath::Erf(sqrt(fabs(prtra->chi2/prtra->NDF))));
      }
      h_ghost_r->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 ){
	h_ghost_tx->Fill((h2.x-h1.x)/(z2-z1));
	h_ghost_ty->Fill((h2.y-h1.y)/(z2-z1));
      }
    }
    
  }//reco tracks

  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;


  for ( map<int, CbmL1MCTrack*>::iterator mtraIt = pmctrackmap.begin();
	mtraIt != pmctrackmap.end(); mtraIt++ ) {
      
      
    // No Sts hits? 
    int nmchits = (mtraIt->second)->StsHits.size();
    if (nmchits == 0) continue;
      
    double momentum = (mtraIt->second)->p;
      
    if (nmchits < 4) continue;// allset  3 or 4 or 5 stations ?
 
    if (momentum < 0.2) continue;// reject very slow tracks from analysis
   
    int sta_count = 1;
    
    vector<int>::iterator previous_hit = (mtraIt->second)->StsHits.begin();
    vector<int>::iterator next_hit = previous_hit; next_hit++;
    
    /*
      {for (vector<L1StsHit*>::reverse_iterator imh = previous_hit;
      imh != (mtraIt->second)->StsHits.rend(); ++imh){
      cout << " " << (*imh)->strip_f->iStation;
      }}
    */

    for (vector<int>::iterator imh = next_hit;
	 imh != (mtraIt->second)->StsHits.end(); ++imh){
      CbmL1HitStore &mh = vHitStore[*imh];
      CbmL1HitStore &ph = vHitStore[*previous_hit];
      if (mh.iStation - ph.iStation == 1) sta_count++;      
      else sta_count = 1;      
      if (sta_count == 4) break;
      previous_hit = imh;
    }

    //cout << " => " << sta_count << endl;
    if (sta_count < 4) continue;//less than 4 continuous hits
    
    
    //continuous stations
    int iph = (mtraIt->second)->StsHits[0];
    CbmL1HitStore &ph = vHitStore[iph];
    
    h_sec_r->Fill(sqrt(fabs(ph.x*ph.x+ph.y*ph.y)));
    if( fabs( (mtraIt->second)->pz )>1.e-8 ){
      h_tx->Fill((mtraIt->second)->px/(mtraIt->second)->pz);
      h_ty->Fill((mtraIt->second)->py/(mtraIt->second)->pz);
    }

    h_mom_track->Fill(momentum);

    int reco = 0;
    if (pmtramap.find(mtraIt->second) != pmtramap.end()){
      reco = 1;
    }
    if(!reco){
      h_notfound_mom->Fill(momentum);
    }
    if (reco){
      p_eff_all_vs_mom->Fill(momentum, 100.0);
      if ((mtraIt->second)->mother_ID < 0){ // primary
	p_eff_prim_vs_mom->Fill(momentum, 100.0);
      }else{
	p_eff_sec_vs_mom->Fill(momentum, 100.0);
      }
    }else{
      p_eff_all_vs_mom->Fill(momentum, 0.0);
      if ((mtraIt->second)->mother_ID < 0){ // primary
	p_eff_prim_vs_mom->Fill(momentum, 0.0);
      }else{
	p_eff_sec_vs_mom->Fill(momentum, 0.0);
      }
    }


    if (pmtramap.find(mtraIt->second) != pmtramap.end()){
      nclones += pmtramap[mtraIt->second]-1;
    }
    
    if(!reco){
      h_notfound_nhits->Fill(nmchits);
      h_notfound_station->Fill(ph.iStation);
      h_notfound_r->Fill(sqrt(fabs(ph.x*ph.x+ph.y*ph.y)));
      CbmL1HitStore &ph21 = vHitStore[mtraIt->second->StsHits[0]];
      CbmL1HitStore &ph22 = vHitStore[mtraIt->second->StsHits[1]];
      
      double z21 = algo->vStations[ph21.iStation].z[0];
      double z22 = algo->vStations[ph22.iStation].z[0];
      if( fabs(z22-z21)>1.e-4 ){
	h_notfound_tx->Fill((ph22.x-ph21.x)/(z22-z21));
	h_notfound_ty->Fill((ph22.y-ph21.y)/(z22-z21));
      }
    }
    
    if (((mtraIt->second)->mother_ID < 0)&&((mtraIt->second)->z > 0)){ // D0
      nmtra_d0++;
      if (reco) nmtra_d0_reco++;
      if (reco) p_eff_d0_vs_mom->Fill(momentum, 100.0);
      else      p_eff_d0_vs_mom->Fill(momentum,   0.0);
      
      h_reco_d0_mom->Fill(momentum);
    }

    nmtra_total++;
    
    if (reco) nmtra_total_reco++;
          
    if (momentum>1.0){ // reference tracks
      nmtra_fast++;
      if (reco) nmtra_fast_reco++;
	
      if ((mtraIt->second)->mother_ID < 0){ // primary
	nmtra_fast_prim++;
	if (reco) nmtra_fast_prim_reco++;
      }
      else{
	nmtra_fast_sec++;
	if (reco) nmtra_fast_sec_reco++;
      }
    }
    else{   // extra set of tracks
      nmtra_slow++;
      if (reco){ 
	nmtra_slow_reco++;
      }
      
      if ((mtraIt->second)->mother_ID < 0){ // primary
	nmtra_slow_prim++;
	if (reco) nmtra_slow_prim_reco++;
      }
      else{
	nmtra_slow_sec++;
	if (reco) nmtra_slow_sec_reco++;
      }
      
    }
    
  }

  h_reco_time->Fill(algo->CATime);
  h_reco_timeNtr->Fill(nmtra_total,algo->CATime);
  h_reco_timeNhit->Fill(algo->vStsHits.size(),algo->CATime);

  int NFakes = 0;
  for( unsigned ih=0; ih<algo->vStsHits.size(); ih++){
    int iMC = vHitMCRef[ih];
    if (iMC < 0) NFakes++;
  }

  h_reco_fakeNtr->Fill(nmtra_total,NFakes);
  h_reco_fakeNhit->Fill(algo->vStsHits.size()-NFakes,NFakes);
  

  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_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;
  double P_CLONE=0.0,P_GHOST=0.0;
  if (L1_NMTRA_D0!=0)
    P_D0 = double(L1_NMTRA_D0_RECO)/double(L1_NMTRA_D0);
  if (L1_NMTRA_TOTAL!=0)
    P_TOTAL = double(L1_NMTRA_TOTAL_RECO)/double(L1_NMTRA_TOTAL);
  if (L1_NMTRA_FAST!=0)
    P_FAST  = double(L1_NMTRA_FAST_RECO)/double(L1_NMTRA_FAST);
  if (L1_NMTRA_FAST_PRIM!=0)
    P_FAST_PRIM  = double(L1_NMTRA_FAST_PRIM_RECO)/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);
  if (L1_NMTRA_SLOW!=0)
    P_SLOW  = double(L1_NMTRA_SLOW_RECO)/double(L1_NMTRA_SLOW);
  if (L1_NMTRA_SLOW_PRIM!=0)
    P_SLOW_PRIM  = double(L1_NMTRA_SLOW_PRIM_RECO)/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);
  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);
    
    
  cout.setf(ios::fixed);
  cout.setf(ios::showpoint);
  cout.precision(4);
#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 << "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;

#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   : " << P_FAST_PRIM <<" | "<<L1_NMTRA_FAST_PRIM_RECO   << endl;
  cout << "RefSec    efficiency   : " << P_FAST_SEC  <<" | "<<L1_NMTRA_FAST_SEC_RECO    << endl;
  cout << "Refset    efficiency   : " << P_FAST      <<" | "<<L1_NMTRA_FAST_RECO        << endl;
  cout << "Allset    efficiency   : " << P_TOTAL     <<" | "<<L1_NMTRA_TOTAL_RECO       << endl;
  cout << "ExtraPrim efficiency   : " << P_SLOW_PRIM <<" | "<<L1_NMTRA_SLOW_PRIM_RECO   << endl;
  cout << "ExtraSec  efficiency   : " << P_SLOW_SEC  <<" | "<<L1_NMTRA_SLOW_SEC_RECO    << endl;
  cout << "Extra     efficiency   : " << P_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;

  //#endif//XXX


  bool write = false;

  if (L1_NEVENTS < 10){//write at every event
    write = true;
  }

  if (L1_NEVENTS < 100){//write at every 10th event
    if ( L1_NEVENTS % 10 == 0) 
      write = true;
  }

  if (L1_NEVENTS >= 100){//write at every 100th event
    if ( L1_NEVENTS % 100 == 0) 
      write = true;
  }
  //write = 0;
  if (write){
    TDirectory *curr = gDirectory;
    // Open output file and write histograms
    TFile* outfile = new TFile("L1_histo.root","RECREATE");
    outfile->cd();
    writedir2current(histodir);
    outfile->Close();
    curr->cd();
  }

}