#define anal_sim_cxx
#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"
//#include <TPostScript.h>
//#include <iostream.h>
#include <iostream>
#include <math.h>
#include <fstream>
#include <string.h> 
//#include <fstream.h>
#include <stdlib.h>
//#include <string.h>
//#include <stdlib>
#include <string> 
#include "anal_sim.h"
#include "sim_hist.h"
Float_t p_old=0.;
Int_t ntot_mu=0;
Int_t calo_type;

void anal_sim::set_psfile_name(string file_in){
  sprintf(psfile_name,"%s",file_in.c_str());
  cout << " plot are written on  file=" << psfile_name  << endl;
  file_in.erase (file_in.length() - 2, 2);
  sprintf(rootfile_name,"hist_%s%s",file_in.c_str(),"root");
  cout <<" rootfile name "<<rootfile_name<<endl;

}
void anal_sim::Loop(){
  //   In a ROOT session, you can do:
  //      Root > .L anal_sim.C
  //      Root > anal_sim t
  //      Root > t.GetEntry(12); // Fill t data members with entry number 12
  //      Root > t.Show();       // Show values of entry 12
  //      Root > t.Show(16);     // Read and show values of entry 16
  //      Root > t.Loop();       // Loop on all entries
  //

  //     This is the loop skeleton
  //       To read only selected branches, Insert statements like:
  // METHOD1:
  //    fChain->SetBranchStatus("*",0);  // disable all branches
  //    fChain->SetBranchStatus("branchname",1);  // activate branchname
  // METHOD2: replace line
  //    fChain->GetEntry(i);  // read all branches
  //by  b_branchname->GetEntry(i); //read only this branchq

  Int_t n_to_read=100000000;
  //Int_t n_to_read=100;
  Int_t ntot=0;

  const Int_t nbp=10;
  const Int_t  nbe=20;
  const Int_t nkas=2;
  Bool_t mu_eff=1;
  Int_t n1=0;
  Int_t n2=0;
  Int_t n3=0;
  Int_t n2_a=0;
  Float_t pmin[]={1.6,3.0,3.0};
  Float_t pmax[]={3.6,8.0,8.0};
  Float_t y_limit[]={280.,418.,551.};
  // sim_hist * my_histos =new sim_hist(psfile_name,nbp,nbe,nkas,pmin,pmax);
  sim_hist * my_histos =new sim_hist(psfile_name,nbp,nbe,nkas);
  my_histos->sim_hist_init(pmin,pmax);
  // const Int_t nbp=my_histos->nbp;
  //const Int_t nbe=my_histos->nbe;
  //const Int_t nkas=my_histos->nkas;
  Float_t bin_nb=static_cast <Float_t>(nbp);
  Float_t ptmin=my_histos->ptmin;
  Float_t ptmax=my_histos->ptmax;
  Float_t etamin=my_histos->etamin;
  Float_t etamax=my_histos->etamax;
  Float_t pt_bin=my_histos->pt_bin;
  Float_t eta_bin=my_histos-> eta_bin;
  Float_t pbin_A=(pmax[0]-pmin[0])/bin_nb;
  Float_t pbin_B=(pmax[1]-pmin[1])/bin_nb;
  Float_t pbin_C=(pmax[2]-pmin[2])/bin_nb;

  Int_t bad_index=0;
  Float_t phimin=3.93;
  Float_t phimax=5.50;
  Float_t phi_hole_inf=4.42;
  Float_t phi_hole_sup=5.02;
  cout <<" avant test sur fchain"<<endl;

  if (fChain == 0) return;
  cout <<" apres test sur fchain"<<endl;
  ofstream out("my_output");
   cout<<" j'ecris A sur   my_output "<<endl;
   out <<"A "<<endl;
  Int_t nentries = Int_t(fChain->GetEntriesFast());
  Int_t nmax=(nentries > n_to_read)?n_to_read:nentries;
  cout <<" Nmax / Nentries: "<< nmax << " / " << nentries << endl;
  Int_t ifreq=nmax/10;
  Int_t freq_accum=1;
  if(nmax >100)freq_accum=100;
  if(ifreq <=0)ifreq=1;
  ifreq=1;
  cout <<" Processing" << endl;
  Bool_t doprint,debug=0;
  Int_t nbytes = 0, nb = 0;
 
  Int_t Ngener[nbp][nbe][nkas];
  Int_t Ngener_A[nbp][nbe][nkas];
  Int_t Ngener_B[nbp][nbe][nkas];
  Int_t Ngener_C[nbp][nbe][nkas];

  Int_t Naccep_A[nbp][nbe][nkas], Naccep_B[nbp][nbe][nkas], Naccep_C[nbp][nbe][nkas];
  Float_t muoH = 292.059; // distance from the beam to "A" layer

  for (Int_t i=0;i<nbp;i++){
    // cout <<" i "<<i<<endl;
    for (Int_t j=0;j<nbe;j++){
      // cout <<" j" <<j<<endl;
      for (Int_t k=0;k<nkas;k++){
	Ngener[i][j][k]=0;
	Ngener_A[i][j][k]=0;
	Ngener_B[i][j][k]=0;
	Ngener_C[i][j][k]=0;
	Naccep_A[i][j][k]=0;
	Naccep_B[i][j][k]=0;
	Naccep_C[i][j][k]=0;
      }}} 
    Float_t etot_cal;

  for (Int_t jentry=0; jentry<nmax;jentry++) {
    Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the entry number in the current file
    ntot++;
    //    cout <<" j'ecris "<<ientry<<" sur my_file"<<endl;
    //out <<jentry<<endl;
    doprint=ntot<10;
    debug=ntot<-1;
    if(ntot%ifreq == 0 || doprint) cout <<" ientry "<<ientry<<endl;
    if(doprint) cout <<" ientry "<<ientry<<endl;
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    // if (Cut(ientry) < 0) continue;
    float x,y,z;
    float steplength,g3code;
    Float_t edep=0.;
    Float_t ptot=0; 
    //Float_t dist,tof;
    if(doprint){
      cout <<"ncpsh,ncfth,nfpsh,nmuonh,nsmth "<< CPSH_ncpsh<<" "<<CFTH_ncfth<<" "<<FPSH_nfpsh<<" "<<MUONH_nmuonh<<" "<<SMTH_nsmth<<endl;
    }
    ntot_mu=MUONH_nmuonh;

    if(debug){
      cout <<" avant premier fill"<<endl;
      //  Int_t nbinx=my_histos->HSMT_SIM->GetNbinsX();  
      //cout <<" nbinx "<<nbinx<<endl;  
    }
    my_histos->HSMT_SIM ->Fill(SMTH_nsmth);
    my_histos->HCFT_SIM ->Fill(CFTH_ncfth);
    my_histos->HCPS_SIM ->Fill(CPSH_ncpsh);
    my_histos->HFPS_SIM ->Fill(FPSH_nfpsh);
    my_histos->HMUON_SIM->Fill(MUONH_nmuonh);
    if(debug)cout<<" after premier fill"<<endl;
    //
    // Loop on Kine
    //
    Int_t nstab=0;
    Int_t ncent=0;
    Int_t nprim=0;
    Float_t* eta_det= new Float_t[KINE_mpart];
    Int_t* used_track_A = new Int_t[KINE_mpart];
    Int_t* used_track_B = new Int_t[KINE_mpart];
    Int_t* used_track_C = new Int_t[KINE_mpart];
    Float_t* t_lengthA  = new Float_t[KINE_mpart];
    Float_t* t_lengthB  = new Float_t[KINE_mpart];
    Float_t* t_lengthC  = new Float_t[KINE_mpart];
    Int_t* dir_part     = new Int_t[KINE_mpart];
    for (Int_t i = 0;i<   KINE_mpart;i++){
       eta_det[i]=0.;
      used_track_A[i]=0;
      used_track_B[i]=0;
      used_track_C[i]=0;
      t_lengthA[i]=0.;
      t_lengthB[i]=0.;
      t_lengthC[i]=0.;
      dir_part[i]=1;  //forward

      if(KINE_e_kin[i] < 0.1) continue;
      my_histos->H_1011->Fill (float(KINE_id_kin[i]));
      my_histos->H_1012->Fill (KINE_pt_kin[i]);
      my_histos->H_1013->Fill (KINE_eta_kin[i],KINE_phi_kin[i]);
      my_histos->H_1015->Fill (KINE_eta_kin[i]);
      my_histos->H_1016->Fill (KINE_phi_kin[i]);

      if(KINE_prim_kin[i]==1) nprim++;
      if(KINE_stab_kin[i]==1) nstab++;
      if(KINE_cent_kin[i]==1) ncent++;
  
      // accumulation for muon efficiency
      if( !mu_eff)continue;
      n1++;
      if(abs(KINE_id_kin[i])!=14 )continue;
      n2++;
      if(KINE_e_kin[i]<pmin[0]){
	//cout <<" pmuon "<<KINE_e_kin[i]<<endl;
	continue;
      } n2_a++;
      if(doprint )cout <<" pos 1 muon "<<i<<" id "<<KINE_id_kin[i]<<" pt "<<KINE_pt_kin[i]<<" prim "<< KINE_prim_kin[i]<<" pdg_id "<< KINE_pdg_id[i]<<endl;	
      Int_t k=(KINE_phi_kin[i]>phi_hole_inf && KINE_phi_kin[i]< phi_hole_sup)?1:0; 
      // Calculate detector eta for muon "A" layer:
      
      Float_t Zv= VERGEN_zver_gen[KINE_vtx_kin[i]-1];
      Float_t z0      = Zv + muoH*sinh(KINE_eta_kin[i]);
      Float_t detaMuo = log(z0/muoH + sqrt(pow((z0/muoH),2) + 1.));
      eta_det[i]=fabs(detaMuo);
      // cout <<" id_kine "<<KINE_id_kin[i]<<" abs "<<abs(KINE_id_kin[i])<<endl;
 
      if(fabs(eta_det[i]) <etamin || fabs(eta_det[i])>etamax){
	cout <<" pmuon "<<KINE_e_kin[i]<<" eta "<<KINE_eta_kin[i]<<" eta det "<<detaMuo <<" zv z0 "<<Zv <<" , "<<z0<<" etamin etamax "<<etamin<<" "<<etamax<<endl ;

	continue;
      }
  
      n3++;
      //      if(jentry>300)    cout <<" pos 1"<<endl;
      //if(jentry>300)      out <<" pos 1"<<endl;
      Int_t ipb=static_cast <int> ((KINE_e_kin[i]- ptmin)/ pt_bin) ;
      Int_t ipe  =static_cast <int>((eta_det[i]-etamin)/eta_bin);
      if(ipb>=nbp)ipb=nbp-1;
      //if(jentry>300)    cout <<" ipb 1"<<ipb<<endl;
      //if(jentry>300)      out <<" ipb 1"<<endl;
      Ngener[ipb][ipe][k]++;
      ipb =static_cast <int>( (KINE_e_kin[i]- pmin[0])/pbin_A);
      if(ipb <0)continue;
      if(ipb>=nbp)ipb=nbp-1;
      //if(jentry>300)    cout <<" ipb 2"<<ipb<<endl;
      //if(jentry>300)      out <<" ipb 2"<<endl;
      Ngener_A[ipb][ipe][k]++;

      ipb =static_cast <int>( (KINE_e_kin[i]- pmin[1])/pbin_B);
      if(ipb <0)continue;
      if(ipb>=nbp)ipb=nbp-1;
      Ngener_B[ipb][ipe][k]++;

      ipb =static_cast <int>( (KINE_e_kin[i]- pmin[2])/pbin_C);
      if(ipb>=nbp)ipb=nbp-1;
      if(ipb <0)continue;
      //if(jentry>300)    cout <<" ipb 4"<<endl;
      //if(jentry>300)      out <<" ipb 4"<<endl;
      Ngener_C[ipb][ipe][k]++;
      if(Ngener_B[ipb][ipe][k] != Ngener_C[ipb][ipe][k]){
	cout <<" Ngener_B "<<Ngener_B[ipb][ipe][k]<<"Ngener_C "<<  Ngener_C[ipb][ipe][k]<<" ipb,k "<<ipb <<" "<<k<<endl;
}
      //      if(jentry>300)cout <<"dsn kine fin de boucle"<<endl;
      //if(jentry>300)out <<"dsn kine fin de boucle"<<endl;
      if(doprint && KINE_e_kin[i]<2.)
	cout <<" pmuon "<<KINE_e_kin[i]<<" eta "<<KINE_eta_kin[i]<<" eta det "<<detaMuo <<" zv z0 "<<Zv <<" , "<<z0<<" Bin pt "<<ipb<<" bin eta "<<ipe<<endl ;
    } //end of loop on KINE
    //if(jentry>300)cout <<" apres kine"<<endl;
    //if(jentry>300)    out <<" apres kine "<<endl;

    //cout <<" ntot "<<ntot <<" n1 n2 n2a n3 "<<n1<<" "<<n2<<" "<<n2_a<<" "<<n3<<endl;

    if(debug)cout <<" apres loop sur kine"<<endl;
    //
    // Loop on SMT
    //
    Int_t freq=(SMTH_nsmth<500)?1:SMTH_nsmth/500;




    for (Int_t i = 0;i< SMTH_nsmth;i++){

      g3code       = SMTH_g3_smth[i];   
      ptot         = SMTH_ptot_smth[i];   
      steplength   = SMTH_st_smth[i];
      edep         = SMTH_ed_smth[i];   
      x            = SMTH_xg_smth[i];   
      y            = SMTH_yg_smth[i];   
      z            = SMTH_zg_smth[i];   
      Int_t indice = SMTH_ind_smth[i];  

      if(ptot >0.1){
	my_histos-> H_2121 -> Fill (steplength);
	my_histos->H_2131 -> Fill (edep);
	if(steplength>0)  my_histos->H_2141 -> Fill (edep/steplength);
      }else{
	my_histos->H_2122 -> Fill (steplength);
	my_histos->H_2132 -> Fill (edep);
	if(steplength>0)  my_histos->H_2142 -> Fill (edep/steplength);
      }


      if(ptot>0.1){
	if(indice==1) {
      
	  my_histos->H_2111 -> Fill (x,y);
	  my_histos->H_2114 -> Fill (z,sqrt(x*x+y*y));
	  if( z > 0. ) my_histos->  H_2112 -> Fill (x,y);
	  if( z < 0. ) my_histos->  H_2113 -> Fill (x,y);
	}


	if(indice==2) {
	  my_histos->H_2211 -> Fill (x,y);
	  my_histos->H_2214 -> Fill (z,sqrt(x*x+y*y));
	  if( z > 0. )  my_histos->H_2212 -> Fill (x,y);
	  if( z < 0. )  my_histos->H_2213 -> Fill (x,y);
	}

	if(indice==3) {
	  my_histos->H_2311 -> Fill (x,y);
	  my_histos->H_2314 -> Fill (z,sqrt(x*x+y*y));
	  if( z > 0. )  my_histos->H_2312 -> Fill (x,y);
	  if( z < 0. )  my_histos->H_2313 -> Fill (x,y);
	}

      }//end of ptot


    }//end of loop on smt hits} 	

    //end of loop on SMT
    if(debug)cout <<" apres loop sur SMT"<<endl;

    //
    // Loop on CFT
    //
    freq=(CFTH_ncfth<1000)?1:CFTH_ncfth/1000;

    for (Int_t i =0;i<CFTH_ncfth;i++){

      g3code     = CFTH_g3_cfth[i];
      ptot       = CFTH_ptot_cfth[i];   
      steplength = CFTH_st_cfth[i];
      edep       = CFTH_ed_cfth[i];
      x          = CFTH_xg_cfth[i];  
      y          = CFTH_yg_cfth[i];   
      z          = CFTH_zg_cfth[i];   

      if(ptot>0.1){
	my_histos->H_3111 -> Fill (x,y);
	my_histos->H_3112 -> Fill (z,sqrt(x*x+y*y));
	my_histos->H_3121 -> Fill (steplength);
	my_histos->H_3131 -> Fill (edep);
	if(steplength>0.)my_histos->H_3141 -> Fill (edep/steplength);
      }else {	  
 	my_histos->H_3122 -> Fill (steplength);
	my_histos->H_3132 -> Fill (edep);
	if(steplength>0.)my_histos->H_3142 -> Fill (edep/steplength);
      }

    }//end of loop on CFT
    if(debug)cout <<" apres loop sur cft"<<endl;

    //
    // Loop on CPS
    //
    freq=(CFTH_ncfth<1000)?1:CFTH_ncfth/1000;


    for (Int_t i =0;i<CPSH_ncpsh;i++){
      g3code     = CPSH_g3_cpsh[i];  
      ptot       = CPSH_ptot_cpsh[i];  
      steplength = CPSH_st_cpsh[i]; 
      x          = CPSH_xg_cpsh[i]; 
      y          = CPSH_yg_cpsh[i]; 
      z          = CPSH_zg_cpsh[i]; 
      edep       = CPSH_ed_cpsh[i];
      if(ptot > 0.1 ){
	my_histos->H_4111 -> Fill (x,y);
	my_histos->H_4112 -> Fill (z,sqrt(x*x+y*y));

	my_histos->H_4121 -> Fill (steplength);
	my_histos->H_4131 -> Fill (edep);
	if(steplength>0.) my_histos->H_4141 -> Fill (edep/steplength);
      }else{
	my_histos->H_4122 -> Fill (steplength);
	my_histos->H_4132 -> Fill (edep);
	if(steplength>0.) my_histos->H_4142 -> Fill (edep/steplength);
      }
    }//end of loop on CPS
    if(debug)cout <<" apres loop sur CPS"<<endl;
    //


    // Loop on FPS
    //
    freq=(FPSH_nfpsh<1000)?1:FPSH_nfpsh/1000;


    for (Int_t i =0;i< FPSH_nfpsh;i++){
      g3code     = FPSH_g3_fpsh[i];   
      ptot       = FPSH_ptot_fpsh[i];  
      steplength = FPSH_st_fpsh[i];
      edep       = FPSH_ed_fpsh[i];  
      x          = FPSH_xg_fpsh[i];
      y          = FPSH_yg_fpsh[i]; 
      z          = FPSH_zg_fpsh[i]; 
      if( ptot > 0.1){ 
	if(i%freq == 0){
	  my_histos->H_5111 -> Fill (x,y);
	  my_histos->H_5112 -> Fill (z,sqrt(x*x+y*y));
	  my_histos->H_5113 -> Fill (z,sqrt(x*x+y*y));
	}
	my_histos->H_5121 -> Fill (steplength);
	my_histos->H_5131 -> Fill (edep);
	if(steplength>0.) my_histos->H_5141 -> Fill (edep/steplength);
      }else{
	my_histos->H_5122 -> Fill (steplength);
	my_histos->H_5132 -> Fill (edep);
	if(steplength>0.) my_histos->H_5142 -> Fill (edep/steplength);
      }
    }//end of loop on FPS
    if(debug)cout <<" apres loop sur FPS"<<endl;


    //===============================================================
    //-------+
    //| calo |
    //-------+
    if(debug)	cout <<" posit 10 ncal"<<CAL_ncal<<endl;

    etot_cal=0.;
    Float_t etot_cal_type[3];
    etot_cal_type[0]=etot_cal_type[1]=etot_cal_type[2]=0;
    Int_t ntower[3];
    Int_t ntwr[3];// index=0 EM, =1 Hadr, 3=other
    ntower[0]=ntower[1]=ntower[2]=0;
    Float_t edig_layer[17];
    Int_t ndig_layer[17];
    const Int_t ntmax=3000;
    Int_t phit[ntmax][3],etat[ntmax][3];
    Float_t et[ntmax][3];
    // layer_name[]={"EM1" ,"EM2","EM3_1","EM3_2","EM3_3","EM3_4","EM4",
    //                "CCGM","ICD","ECMG" ,"FH1"  ,"FH2"  ,"FH3"  ,"FH4",
    //			"CH1","CH2","CH3"};
    static int layToLyr[] = {0, 1, 2, 2, 2, 2, 3, 4, 5, 6, 7, 7, 7, 7, 8, 8, 8};

    for( Int_t i=0;i< 17 ;i++){
      ndig_layer[i]=0;
      edig_layer[i]=0.;
    }
    freq=(CAL_ncal<3000)?1:CAL_ncal/3000;
    
    for( Int_t i=0;i< CAL_ncal ;i++){
      if(CAL_lyr_cal[i] >21)break;
      etot_cal+=CAL_e_cal[i];
      Int_t layer;
      if(CAL_lyr_cal[i]<20){ 
	layer = layToLyr[CAL_lyr_cal[i]-1];
	if(CAL_lyr_cal[i]<=7)		                calo_type=0;
	else if(CAL_lyr_cal[i]>=11)                     calo_type=1;
	else if(CAL_lyr_cal[i]>=8 &&CAL_lyr_cal[i]<=10) calo_type=2;
      }else{
	layer=( CAL_lyr_cal[i]==20)?1:11;
	calo_type=CAL_lyr_cal[i]-20;
      }
      etot_cal_type[calo_type]+=CAL_e_cal[i];

      Int_t lyr=CAL_lyr_cal[i];
	
      if(CAL_lyr_cal[i]<20){
	//EM1 EM2 (lyr = 1, 2)
	if      (CAL_lyr_cal[i]>= 3 && lyr <= 6)        lyr = 3;     //EM3
	else if (CAL_lyr_cal[i] == 7)                   lyr = 4;     //EM4
	else if (CAL_lyr_cal[i] >= 11 && lyr <=14 )     lyr = lyr-6; //Fine Hadronic
	else if (CAL_lyr_cal[i] >= 15 && CAL_lyr_cal[i]<20) lyr = 9; //Coarse hadronic
	else if (CAL_lyr_cal[i] == 8)                   lyr = 10;    //CCMG
	else if (CAL_lyr_cal[i] == 9)                   lyr = 11;    //ICD
	else if (CAL_lyr_cal[i] == 10)                  lyr = 12;    //ECMG
      }
      else lyr =1+(CAL_lyr_cal[i]-20)*4;

      lyr=lyr-1; 

      ndig_layer[lyr]++;  
      edig_layer[lyr]+= CAL_e_cal[i];
     
      if(lyr <=8){
	//cout <<"lyr "<<lyr<<endl;     
	my_histos-> H_El_dig[lyr]->Fill(CAL_e_cal[i]);
	my_histos-> H_Eh_dig[lyr]->Fill(CAL_e_cal[i]);
      }
      if(i%freq_accum == 0){
	my_histos->H_etaxphi[calo_type*3]->Fill(CAL_ieta_cal[i],CAL_iphi_cal[i]);
	my_histos->H_etaxphi[calo_type*3+1]->Fill(CAL_ieta_cal[i],CAL_iphi_cal[i],CAL_pt_cal[i]);
      }
      Int_t newt=1;
      if( ntower[calo_type] <ntmax){
  
    
	for(Int_t tow=0;tow<ntower[calo_type];tow++){
	  if(CAL_iphi_cal[i] != phit[tow][calo_type] || CAL_ieta_cal[i] != etat[tow][calo_type])continue;

	  et[tow][calo_type]+= CAL_pt_cal[i];
	  newt=0;
	  break;  
	}//end of loop on towers
	if(newt==1){
	  Int_t ntt=ntower[calo_type];
	  phit[ntt][calo_type]=CAL_iphi_cal[i];
	  etat[ntt][calo_type]=CAL_ieta_cal[i];
	  et[ntt][calo_type]= CAL_pt_cal[i]; 
	  ntwr[calo_type]++;
	  ntower[calo_type]++;

	}//end of test on new tower
      }//end of test on nb. of towers

    }//end of loop on calo

    my_histos->cal_ndig->Fill(CAL_ncal);
    my_histos-> cal_edig->Fill(etot_cal); 

    //Energy and nb. of hits per layer
    for( Int_t i=0;i< 12 ;i++){
      //	  if(ndig_layer[i]<=0)continue;
      my_histos-> cal_ndig_layer->Fill(i+1,ndig_layer[i]);
      my_histos-> cal_edig_layer->Fill(i+1,edig_layer[i]);
    }
	
    for(Int_t it=0;it<3;it++){
      my_histos->  cal_etot[it]->Fill(etot_cal_type[it]);
      my_histos-> cal_ntower[it]->Fill(ntower[it]);
      for( Int_t i=0;i< ntower[it] ;i++){
	my_histos-> cal_etower[it]->Fill(et[i][it]);
	my_histos-> eta_phi_tow[it*3]->Fill( etat[i][it], phit[i][it]);
	my_histos-> eta_phi_tow[it*3+1]->Fill( etat[i][it], phit[i][it],et[i][it]);  
      }//end of loop on towers
    }//end of loop on calo_type


    //------------+
    //|end of calo|
    //+-----------+
    if(debug)cout <<" apres loop sur calo"<<endl;
    //===============================================================

    //
    // Loop on muon chambers
    //
    Bool_t first_hit_A[2];
    first_hit_A[0] =true;
    first_hit_A[1] =true;
    Bool_t first_hit_B[2];
    first_hit_B[0] =true;
    first_hit_B[1] =true;

    Bool_t first_hit_C[2];
    first_hit_C[0] =true;
    first_hit_C[1] =true;
    if(doprint)cout <<" MUONH_nmuonh "<<MUONH_nmuonh<<endl;
    //freq=(MUONH_nmuonh<100)?1:MUONH_nmuonh/100;
    freq=1;
    //    cout <<"avant les muons "<<endl;
    //out <<"avant les muons "<<endl;
    Float_t xa,ya,za;
    Float_t xb,yb,zb;
    Float_t xc,yc,zc;

    for (Int_t i =0;i<MUONH_nmuonh;i++){
      //if(i>=60)cout <<" retour en debut de boucle mu "<<endl;
      g3code       = MUONH_g3_muh[i];  
      Float_t px   = MUONH_pxg_muh[i];
      Float_t py   = MUONH_pyg_muh[i];
      Float_t pz   = MUONH_pzg_muh[i];
      ptot         = sqrt(px*px+py*py+pz*pz);  
      steplength   =  MUONH_step_muh[i];  
      x            =  MUONH_xg_muh[i]; 
      y            =  MUONH_yg_muh[i]; 
      z            =  MUONH_zg_muh[i]; 
      edep         =  MUONH_edep_muh[i];
      Int_t indice =  MUONH_ind_muh[i]%10-1;

      Float_t dx=steplength*fabs(px)/ptot;
      Float_t dy=steplength*fabs(py)/ptot;
      Float_t dz=steplength*fabs(pz)/ptot;

if(indice<0){
	bad_index++;
	if(bad_index!=0 && bad_index%100 ==0)cout <<"  MUONH_ind_muh"<<MUONH_ind_muh[i]<<" indice "<<indice<<" file " <<psfile_name<<" nb. of evts with bad index "<<bad_index<<endl;
	continue;
      }  
      Int_t j= MUONH_itr_muh[i]-1;
      if(ptot > 0.1){
        my_histos-> H_mu_xy[indice] -> Fill (x,y);
	my_histos-> H_mu_zr[indice] -> Fill (z,sqrt(x*x+y*y));
	my_histos-> H_mu_step[indice] -> Fill (steplength);
	my_histos-> H_mu_e[indice]    -> Fill (edep);
	if (steplength>0.)my_histos-> H_mu_enorm[indice] -> Fill (edep/steplength);
	//Int_t j=MUONH_itr_muh[i]-1;
      }else{	
	my_histos->H_mu_step_d[indice] -> Fill (steplength);
	my_histos->H_mu_e_d[indice]    -> Fill (edep);
	if (steplength>0.)my_histos->H_mu_enorm_d[indice]->Fill (edep/steplength);
      }
      if(!mu_eff)continue;

if(MUONH_g3_muh[i]!=5 && MUONH_g3_muh[i]!=6)continue;
      if(abs(KINE_id_kin[j]) !=14)continue;
      if(ptot <0.1)continue;
      p_old=ptot;
      if(KINE_e_kin[j]< ptmin )continue;
      if(eta_det[j] <etamin || eta_det[j]>etamax)continue;
      if(indice ==1 || indice==3)continue;//select chambers (exclude scintillators)

      Int_t k=(KINE_phi_kin[j]> phi_hole_inf && KINE_phi_kin[j]<phi_hole_sup)?1:0; 
      Int_t ipb =static_cast <int>( (KINE_e_kin[j]- pmin[0])/ pbin_A);
      if(ipb>=nbp)ipb=nbp-1;
      if(ipb <0)continue;
      Int_t ipe =static_cast <int>((eta_det[j]-etamin)/ eta_bin);
      Float_t diff =0.;
      if(doprint)cout <<" MUONH_ind_muh "<<MUONH_ind_muh[i]<<" first_hit "<<first_hit_A[k]<<" used_track "<<used_track_A[j]<<" "<<used_track_B[j]<<" "<<used_track_C[j]<<endl;
      diff =  KINE_e_kin[j] -ptot;
      if(diff<0.){
	cout <<" pmuon "<<ptot<<" e_kine "<<KINE_e_kin[j]  <<" diff p "<<  KINE_e_kin[j] -ptot<<" eta part "<< KINE_eta_kin[j]<<" eta_det "<<eta_det[j]<<" index "<<MUONH_ind_muh[i]<<" kine mpart "<<KINE_mpart<<" x y z "<<x<<" , "<<y<<" ,"<<z<<endl;
	for (Int_t ik = 0;ik<   KINE_mpart;ik++){
	  cout <<" pos 2 muon "<<ik<<" id "<<KINE_id_kin[ik]<<" pt "<<KINE_e_kin[ik]<<" pdg_id "<< KINE_pdg_id[ik]<<endl;	
	}
	continue;
      }//end of test on diff
      Float_t dst=0.;

      if(MUONH_ind_muh[i]==11 || MUONH_ind_muh[i]==13){
        
	if( MUONH_ind_muh[i]==11){
          dir_part[j]=0;  //central
	  dst=(fabs(y) < y_limit[0])?dx:dy;
	}else  dst=dz;
	t_lengthA[j]+=fabs(dst);
        if(doprint)cout <<" layer A projected length "<<dst<<" total length "<<t_lengthA[j]<<endl;
	//	cout <<" initial muon  nb: "<<j<<" pmuon "<<ptot<<" e_kine "<<KINE_e_kin[j]  <<" eta part "<< KINE_eta_kin[j]<<" index "<<MUONH_ind_muh[i]<<" x y z "<<x<<" , "<<y<<" ,"<<z<<" used_track_A"<<used_track_A[j]<<endl;
	if(used_track_A[j]==1)        continue;

	xa=x;
        ya=y;
        za=z;
        used_track_A[j]=1;
	first_hit_A[k]=false;
	if(k == 1 && MUONH_ind_muh[i]==11 ){
	  cout <<" in hole layer A:  MUONH_ind_mu "<< MUONH_ind_muh[i]<<" indice "<<indice<<" ptot "<<ptot<<" xg yg zg "<<x<<" , "<<y<<" ,"<<z<<endl;
	  cout <<" pmuon "<<ptot<<" eta part "<< KINE_eta_kin[j]<<" eta_det "<<eta_det[j]<<" phi "<<KINE_phi_kin[j]<<endl;
	  Float_t phip=KINE_phi_kin[j];
	  if(phip >=phimin && phip<4.5)phimin=phip;
	  if(phip <=phimax && phip>5.0)phimax=phip;
	  cout <<" phimin "<<phimin <<" phimax "<<phimax <<" phi mu "<<atan2(-y,-x)+3.14<<endl;
	}
	if(Naccep_A[ipb][ipe][k] > Ngener_A[ipb][ipe][k]){
	  cout <<" ipb ipe k "<< ipb <<","<<ipe <<","<<k<<" naccep "<<Naccep_A[ipb][ipe][k]<<" ngener "<<Ngener_A[ipb][ipe][k]<<endl;
	}//end of test on naccep
	if(doprint)cout <<" avant accum de naccep "<<endl;
	Naccep_A[ipb][ipe][k]++;
        my_histos-> mu_E_loss[ipb][0]->Fill(eta_det[i],KINE_e_kin[j]-ptot);
	if(doprint){
	  cout <<" accum de  H_acc: laye 0: ipb , k "<<ipb<<" "<<k<<endl;
	  //  out <<" accum de  H_acc: ipb , k ,layer"<<ipb<<" "<<k<<"  0 "<<endl;
	}//end of test on doprint
	if(doprint)cout <<" after A dir_part "<<dir_part[j]<<endl;
      }else if(MUONH_ind_muh[i]==21 || MUONH_ind_muh[i]==23){
	if( MUONH_ind_muh[i]==21){
          dir_part[j]=0;  //central
	  dst=(fabs(y) < y_limit[1])?dx:dy;
	}else  dst=dz;

	t_lengthB[j]+=fabs(dst);
        if(doprint)cout <<" layer B projected length "<<dst<<" total length "<<t_lengthB[j]<<endl;

	if(used_track_B[j]==1)continue;
	used_track_B[j]=1;
	first_hit_B[k]=false;
	ipb =static_cast <int>( (KINE_e_kin[j]- pmin[1])/ pbin_B);
	if(ipb <0)continue;
	if(ipb>=nbp)ipb=nbp-1;
	Naccep_B[ipb][ipe][k]++;
        my_histos-> mu_E_loss[ipb][1]->Fill(eta_det[i],KINE_e_kin[j]-ptot);

	if(doprint)cout <<" accum de  H_acc: laye 1: ipb , k "<<ipb<<" "<<k<<endl;
		if(doprint)cout <<" after B dir_part "<<dir_part[j]<<endl;


      }else if(MUONH_ind_muh[i]==31 || MUONH_ind_muh[i]==33){
	if( MUONH_ind_muh[i]== 31){
          dir_part[j]=0;  //central
	  dst=(fabs(y) < y_limit[2])?dx:dy;
	}else  dst=dz;
	t_lengthC[j]+=fabs(dst);
        if(doprint)cout <<" layer C projected length "<<dst<<" total length "<<t_lengthC[j]<<endl;
	if(used_track_C[j]==1)continue;
	used_track_C[j]=1;
	first_hit_C[k]=false;
	Int_t ipb =static_cast <int>( (KINE_e_kin[j]- pmin[2])/ pbin_C);
	if(ipb <0)continue;
	if(ipb>=nbp)ipb=nbp-1;
	Naccep_C[ipb][ipe][k]++;  
        my_histos-> mu_E_loss[ipb][2]->Fill(eta_det[i],KINE_e_kin[j]-ptot);

	if(doprint)cout <<" accum de  H_acc: laye 2: ipb , k "<<ipb<<" "<<k<<endl;

      }//end of test on  C Layer
	if(doprint)	cout <<" after C dir_part "<<dir_part[j]<<endl;

      //cout <<" after test on ptot "<<endl;
    }//end of loop on muons 
    if(debug) cout <<"after loop on muons"<<endl;
    if( mu_eff){
      for (Int_t i = 0;i<   KINE_mpart;i++){
      if(KINE_e_kin[i] < 0.1) continue;
      if(abs(KINE_id_kin[i])!=14 )continue;
      if(fabs(eta_det[i]) <etamin || fabs(eta_det[i])>etamax)continue;
Int_t k=(KINE_phi_kin[i]>phi_hole_inf && KINE_phi_kin[i]< phi_hole_sup)?1:0; 
      Int_t d=dir_part[i];
     Int_t ipb =static_cast <int>( (KINE_e_kin[i]- pmin[0])/pbin_A);
      if(ipb >-1){
      if(ipb>=nbp)ipb=nbp-1;
      my_histos-> E_dep_cal[ipb]->Fill(eta_det[i],etot_cal);
        my_histos->range_A[ipb][d][k]->Fill(t_lengthA[i]);
      my_histos-> E_dep_cal[ipb]->Fill(eta_det[i],etot_cal);
      }
       ipb =static_cast <int>( (KINE_e_kin[i]- pmin[1])/pbin_B);
      if(ipb <0)continue;
      if(ipb>=nbp)ipb=nbp-1;
      my_histos->range_B[ipb][d]->Fill(t_lengthB[i]);
      my_histos->range_C[ipb][d]->Fill(t_lengthC[i]);
       if(doprint)cout<<" part. direction "<<dir_part[i]<<" total lengths in muon chambers "<< t_lengthA[i] <<" , "<<  t_lengthB[i] <<" , "<<     t_lengthC[i] <<" k hole "<<k<<" momentum "<<KINE_e_kin[i]<<endl;
   
 }//end of loop on kine 
    }//end of test on mu_eff

    if(eta_det != NULL)delete [] eta_det;
    if (used_track_A != NULL)delete [] used_track_A;
    if (used_track_B != NULL)delete [] used_track_B;
    if (used_track_C != NULL)delete [] used_track_C;
    if (t_lengthA!= NULL)delete [] t_lengthA;
    if (t_lengthB!= NULL)delete [] t_lengthB;
    if (t_lengthC!= NULL)delete [] t_lengthC;
    if (dir_part!= NULL)delete [] dir_part;
    //+------------------------------+
  }// |end of loop on events         |
   // +------------------------------+
  if(bad_index!=0 )cout<<"  nb of bax index MUONH_ind_muh "<<bad_index<<" file " <<psfile_name<<endl;
  if(debug) cout <<"after loop on events"<<endl;

  //cout <<"j'ecris C sur my_output"<<endl;
  //out <<"C"<<endl;
  //
  //cout <<" avant definition de psfile "<<endl;
  //TPostScript* psfile= new  TPostScript("file_out.ps",112);
  //cout <<" apres definition de psfile "<<endl;
  //cout<<" pt vs eta  genere"<<endl;
  if( mu_eff){
    for  (Int_t k=0;k<nkas;k++){
      Int_t sgener=0;
      Int_t sgener_A=0;
      Int_t sgener_B=0;
      Int_t sgener_C=0;
      Int_t saccepA=0; 
      Int_t saccepB=0;
      Int_t saccepC=0;
   
      for (Int_t i=0;i<nbp;i++){
	for (Int_t j=0;j<nbe;j++){
	  sgener+= Ngener[i][j][k];
	  sgener_A+= Ngener_A[i][j][k];
	  sgener_B+= Ngener_B[i][j][k];
	  sgener_C+= Ngener_C[i][j][k];
	  saccepA+= Naccep_A[i][j][k];
	  saccepB+= Naccep_B[i][j][k];
	  saccepC+= Naccep_C[i][j][k];

          my_histos->H_gen[i][k][0]-> SetBinContent(j, Ngener_A[i][j][k]);
          
          my_histos->H_acc[i][k][0]->SetBinContent(j,Naccep_A[i][j][k]) ;
          my_histos->H_gen[i][k][1]->SetBinContent(j,Ngener_B[i][j][k]);
          my_histos-> H_acc[i][k][1]->SetBinContent(j,Naccep_B[i][j][k]) ;
          my_histos->H_gen[i][k][2]->SetBinContent(j,Ngener_C[i][j][k]);
          my_histos-> H_acc[i][k][2]->SetBinContent(j,Naccep_C[i][j][k]) ;
	  // cout <<" i,j "<<i<<","<<j<<" bin content h_gen_a "<<my_histos->H_gen[i][k][0]
	  //  -> GetBinContent(j)<<" ngener "<<Ngener_A[i][j][k]<<endl;
	}
                
      }
      cout <<"sgener "<<sgener_A<<","<<sgener_B<<","<<sgener_C<<endl;
      cout <<"saccept "<<saccepA<<","<<saccepB<<","<<saccepC<<endl;

      //  cout <<" case "<<k <<" Ngener_A  "<<sgener_A<<endl;
      //      out <<" case "<<k <<" Ngener_A  "<<sgener_A<<endl;
 
      for (Int_t i=0;i<nbp;i++){
	//for (Int_t j=0;j<nbe;j++)cout <<"|"<<setw(4)<<Ngener_A[i][j][k];
	for (Int_t j=0;j<nbe;j++) out <<"|"<<setw(4)<<Ngener_A[i][j][k];
	//cout <<"|"<<endl;
	out <<"|"<<endl;
      }//end of loop on Pt bins
      cout<<" pt vs eta  accepte"<<endl;
      out<<" pt vs eta  accepte"<<endl;
      cout <<"accept layer A "<<saccepA<<endl;
      out <<"accept layer A "<<saccepA<<endl;

      for (Int_t i=0;i<nbp;i++){
	for (Int_t j=0;j<nbe;j++){
	  //cout <<"|"<<setw(4)<< Naccep_A[i][j][k];
	  out <<"|"<<setw(4)<<  Naccep_A[i][j][k];
	}
	//cout <<"|"<<endl;
	out <<"|"<<endl;
      }//end of loop on PT bins
      cout <<" case "<<k <<" Ngener_B  "<<sgener_B<<endl;
      out <<" case "<<k <<" Ngener_B  "<<sgener_B<<endl;
 
      for (Int_t i=0;i<nbp;i++){
	//for (Int_t j=0;j<nbe;j++)cout <<"|"<<setw(4)<<      Ngener_B[i][j][k];
        for (Int_t j=0;j<nbe;j++) out <<"|"<<setw(4)<<      Ngener_B[i][j][k];
	//	cout <<"|"<<endl;
	out <<"|"<<endl;
      }//end of loop on Pt bins

      cout <<"accept layer B"<<saccepB<<endl;
      out <<"accept layer B "<<saccepB<<endl;

      for (Int_t i=0;i<nbp;i++){
	//	for (Int_t j=0;j<nbe;j++) cout <<"|"<<setw(4)<<      Naccep_B[i][j][k];
	for (Int_t j=0;j<nbe;j++)  out <<"|"<<setw(4)<<      Naccep_B[i][j][k];
	//cout <<"|"<<endl;
	out <<"|"<<endl;
      }//end of loop on PT bins
      cout <<" case "<<k <<" Ngener_C  "<<sgener_C<<endl;
      out <<" case "<<k <<" Ngener_C  "<<sgener_C<<endl;
 
      for (Int_t i=0;i<nbp;i++){
	//	for (Int_t j=0;j<nbe;j++) cout <<"|"<<setw(4)<<      Ngener_C[i][j][k];
	for (Int_t j=0;j<nbe;j++)  out <<"|"<<setw(4)<<      Ngener_C[i][j][k];
	//cout <<"|"<<endl;
	out <<"|"<<endl;
      }//end of loop on Pt bins

      cout <<"accept layer C"<<saccepC<<endl;
      out <<"accept layer C "<<saccepC<<endl;

      for (Int_t i=0;i<nbp;i++){
	//for (Int_t j=0;j<nbe;j++) cout <<"|"<<setw(4)<<      Naccep_C[i][j][k];
        for (Int_t j=0;j<nbe;j++)  out <<"|"<<setw(4)<<      Naccep_C[i][j][k];
	//cout <<"|"<<endl;
	out <<"|"<<endl;
      }//end of loop on PT bins

    }//end of loop on cases 
    string prob_name []={"int probA_no_hole[]={",
			 "int probA_hole[]={",
			 "int probB_no_hole[]={",
			 "int probB_hole[]={",
			 "int probC_no_hole[]={",
			 "int probC_hole[]={"};
    out <<" //----------cut and paste from this line-----------------" <<endl;

    out <<"void muon_prob(float ptrack,float phi,float eta, float Zv, float & probA ,float & probB, float & probC){"<<endl;
   out <<"/*"<<endl;
out <<"----------------------------------------------------------------------"<<endl;

   out <<"eta min "<<etamin <<" etamax "<<etamax<<" width eta "<<eta_bin<<endl;
   out << " line = (Ptrack - pmin[0])/binw[0]  for layer A"<<endl;
   out << " line = (Ptrack - pmin[1])/binw[1]  for layer B"<<endl;
   out << " line = (Ptrack - pmin[2])/binw[2]  for layer C"<<endl;
   out << "    Column = (Eta_det _ "<<etamin<<")/"<<eta_bin<<endl;
    out <<" if (line >=nbp) line=nbp;"<<endl;
     out <<" bin defined as column + line*"<<nbe<<endl;
    out<<"  hole defined as: "<<phi_hole_inf<<"<phi_part<"<<phi_hole_sup<<endl;
    out <<"- phi (at DCA) is used to compute what table to be used ( for instance"<<endl;
out <<"A_hole or A_no_hole) "<<endl;
     out <<"-The probability to be seen in layer A is"<<endl;
     out <<" float prob =probA_no_hole[column+20*line]  (for phi outside the hole)"<<endl;
     out <<"or float prob =probA_hole[column+20*line]     (for phi in the hole)"<<endl;


out <<" - calculation of eta_det:(same for 3 layers)"<<endl;

out <<"      // Calculate detector eta for muon <A> layer:"<<endl;
out <<"      Float_t Zv= VERGEN_zver_gen[KINE_vtx_kin[i]-1];"<<endl;
out <<"      Float_t z0      = Zv + muoH*sinh(KINE_eta_kin[i]);"<<endl;
out <<"      Float_t detaMuo = log(z0/muoH + sqrt(pow((z0/muoH),2) + 1.));"<<endl;
out <<"      eta_det[i]=fabs(detaMuo);"<<endl;
out <<"       with  muoH = 292.059; // distance from the beam to <A> layer"<<endl;
 out <<"      -----------------note the absolute value ---------------------"<<endl;
 out <<" "<<endl;
out <<"   -The tables represent the percentage of time, integrated over phi distinguishing the region"<<endl;
out <<" in and out the bottom hole that  a muon generated with p, eta at Zv has at least one hit in"<<endl;
out <<" layer A, B and C."<<endl;
 

    out <<"-eta_detector is computed and, if |eta|> etamax ( ="<<etamax<<" in our case) then the muon is ;"<<endl;
   out <<"considered as outside the acceptance."<<endl;
   out <<"-The tables extend from "<<pmin[0] <<"GeV/c  to "<<pmax[0]<<" for layer A"<<endl;
   out <<"-The tables extend from "<<pmin[1] <<"GeV/c  to "<<pmax[1]<<" for layer B"<<endl;
   out <<"-The tables extend from "<<pmin[2] <<"GeV/c  to "<<pmax[2]<<" for layer C"<<endl;
   out <<"particles with P above P max are considered to have the same probability tables as"<<endl; 
   out <<"particles with p= pmax"<<endl;
   out <<" -the probabilities ARE NOT independant"<<endl;
   out <<" -assumptions (not correct when part. in phi hole)"<<endl;
   out <<"     if a particle reaches layer C  it has hits in A and B"<<endl;
   out <<"     if      --     does not reach layer C but reaches layer B then it has hits in A"<<endl;

 out <<" Input : ptrack, phi (between 0 and 2*pi) , eta  at DCA Z vertex"<<endl;
 out <<" Output: 0 or 1 for the 3 layers"<<endl;
 out <<" "<<endl; 
 out <<" Written by A. Zylberstejn July 2003"<<endl;
 out <<" Latest version based on MC single muons generated by Lei (Oct 2003) "<<endl;
out <<"----------------------------------------------------------------------"<<endl;
   out <<"*/"<<endl; 
   out <<" "<<endl;
    out <<"int nbp ="<<nbp<<";"<<endl;
    out <<"float binw[] = {"<< pbin_A<<","<< pbin_B<<","<<pbin_C<<"};"<<endl;
    out <<"float pmin[]={"<<pmin[0]<<","<<pmin[1]<<","<<pmin[2]<<"};"<<endl;
    out <<"float pmax[]={"<<pmax[0]<<","<<pmax[1]<<","<<pmax[2]<<"};"<<endl;
    out <<"float eta_min ="<<etamin <<";"<<endl;
    out <<"float eta_max ="<<etamax<<";"<<endl;
    out <<"float eta_bin= "<<eta_bin<<";"<<endl;
    out <<"float muoH = "<<muoH<<";"<<endl;
    out<<" float phi_min= "<<phi_hole_inf<<";"<<endl;
    out<<" float phi_max= "<<phi_hole_sup<<";"<<endl;

    for (Int_t layer=0;layer<3;layer++){
      for (Int_t k=0;k<2;k++){
	out <<prob_name[k+layer*2] <<endl;
	cout <<" layer "<<layer <<" cas "<<k<<prob_name[k+layer*2] <<endl;
	cout<<" nbp, nbe "<<nbp <<" , "<<nbe<<endl;
	for (Int_t i=0;i<nbp;i++){
	  Float_t eta_ref= etamin+0.5*eta_bin;
	  for (Int_t j=0;j<nbe;j++){
	    Float_t num=0.;
	    Int_t prob=0;
	    //cout <<"i j k "<<i<<" , "<<j<<" , "<<k<< " negener "<< Ngener[i][j][k] <<" layer "<<layer <<endl;
	    if( Ngener[i][j][k] >=1){
	      Float_t den=static_cast <float>( Ngener[i][j][k]);
	      switch (layer){
	      case 0 :
		den=static_cast <float>( Ngener_A[i][j][k]);
		num=static_cast <float>( Naccep_A[i][j][k]);
		break;
	      case 1 :
		den=static_cast <float>( Ngener_B[i][j][k]);
		num=static_cast <float>( Naccep_B[i][j][k]);
		break;
	      case 2 :
		den=static_cast <float>( Ngener_C[i][j][k]);
		num=static_cast <float>( Naccep_C[i][j][k]);
		break;

	      }
	      // cout <<" den,num "<<den<<" , "<<num<<endl;
                
	      Float_t rap =( den>=1)?num/den:0.;
	      //	cout <<" avant accum de h_prob eta_ref "<<eta_ref<<" rap "<<rap<<" i,k,layer "<<i<<","<<k<<","<<layer<<endl;
	      //my_histos->H_prob[i ][k][layer]->Fill(eta_ref,rap);
	      eta_ref+=eta_bin;
	      //	cout <<" apres accum "<<endl;
	      prob =static_cast <Int_t>(100*rap);
	      //cout <<" prob "<<prob<<endl;	
	    }//end of test on number of ngener
	    //	    cout<<"|"<<setw(3)<<prob;
	    //cout <<"|"<<  std::setw(5)<<std::setprecision(2)<<    prob[i][j];
	    //	cout <<","<< setw(5)<<setprecision(2)<<    prob[i][j];
	    out<< setw(3)<<    prob ;
	    if(j!=nbe-1 || i!=nbp-1) out <<",";
	    if(prob>100){
	      cout <<" prob " <<prob<<" ngener "<<  Ngener[i][j][k]<<" naccep "<<Naccep_A[i][j][k]<<" i,j,k "<<i<<","<<j<<","<<k<<endl;
	      out <<" prob " <<prob<<" ngener "<<  Ngener[i][j][k]<<" naccep "<<Naccep_A[i][j][k]<<" i,j,k "<<i<<","<<j<<","<<k<<endl;
	    }
	  }//end of loop on eta
	  out<<endl;
	}//end of loop on nbp
  
	//out<<"};"<<endl;
	cout <<" end of loop on cases "<<endl;
      }//end of loop on cases
    }//end of loop on layers
    out <<" //----------cut and paste to this line-----------------" <<endl;

    cout <<" end of loop on layers "<<endl;
    my_histos->sim_hist_draw_effic(nentries, ntot,psfile_name);

    //  outputfile->Close();
  }//end of test on mu_eff
  //my_histos->sim_hist_draw(nentries, ntot,psfile_name);

  //my_histos->sim_hist_write(ntot,rootfile_name);
  //my_histos->_ps_file->Close();
  cout <<" end of anal_sim"<<endl;

}