#define anal_rec_cxx
#include "anal_rec.h"
#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"
#include <iostream>
#include <math.h>
#include <fstream>
//#include <stdlib.h>
#include <string>
#include "rec_hist.h"
using namespace std;
void anal_rec::set_psfile_name(std::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);
    cout <<" file_in apres reduction "<<file_in<<endl;
    sprintf(rootfile_name,"hist_%s%s",file_in.c_str(),"root");
    cout <<" rootfile name "<<rootfile_name<<endl;

}

void anal_rec::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L anal_rec.C
//      Root > anal_rec 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 where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       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(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
   if (fChain == 0) return;

   Int_t nentries = Int_t(fChain->GetEntriesFast());
    rec_hist * my_histos =new rec_hist(psfile_name);
    my_histos->rec_hist_init();
    Int_t ntot=0;
    Int_t n_to_read=10000;         // 0 to read all evts.
    Int_t nmax=(nentries > n_to_read) ? n_to_read:nentries;
    cout <<" Nmax / Nentries: "<< nmax << " / " << nentries << endl;
    Int_t ifreq=nmax/10;
    if(ifreq <=0)ifreq=1;
    cout <<" Processing" << endl;
    Bool_t doprint,debug=0;
    ntot++;
    doprint=ntot<10 ;
    debug=ntot>-1;
   Int_t nbytes = 0, nb = 0;
   for (Int_t jentry=0; jentry<nmax;jentry++) {
      Int_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;
      //	 cout <<" ientry "<<jentry<<" nb of clusters: in CFT"<<	     CFTC_ncftc<<" in CPS "<<CPSC_ncpsc<<" njets "<<JCCB_JCCBnjet<<" nmuons "<<Muonid_Nmuonid<<" ntracks "<<TRACK_ntra<<endl;
 	if (ientry < 0) break;
        ntot++;
        doprint=ntot<10 ;
        debug=ntot<-1;
	if(jentry%ifreq == 0 || doprint) cout <<" ientry "<<jentry<<" nb of clusters: in CFT "<< CFTC_ncftc<<" in CPS "<<CPSC_ncpsc<<" njets "<<JCCB_JCCBnjet<<" nmuons "<<Muonid_Nmuonid<<" ntracks "<<TRACK_ntra<<endl;
	nb = fChain->GetEntry(jentry);   nbytes += nb;
	// if (Cut(ientry) < 0) continue;
	//+-----+
	//| CFT | 
	//+-----+

	//	const int nfib[16] ={1280,1600,1920,2240,2560,2880,3200,3520,
	//	     1280,1920,2560,3200,1600,2240,2880,3520};
	Int_t nclus [16];
	Int_t nstrip [16];
	for(Int_t i=0;i<16;i++){
	  nclus [i]=0;
          nstrip[i]=0;
    }
	if(debug)cout <<" avant histos CFT"<<endl;
   my_histos->cft_nclus->Fill(CFTC_ncftc);
	for (Int_t i =0;i< CFTC_ncftc;i++){
	    Int_t layer=CFTC_lay_cftc[i];
	    nclus [layer]++;
	    Float_t phi=(CFTC_phi_cftc[i]>0)?CFTC_phi_cftc[i]:6.28318+CFTC_phi_cftc[i];
            Float_t x_cft=CFTC_r_cftc[i]*cos(phi);
            Float_t y_cft=CFTC_r_cftc[i]*sin(phi);
	    nstrip[layer]+= CFTC_nst_cftc[i];
            my_histos->cft_xy_clust->Fill(x_cft,y_cft);
	    //           my_histos->cft_rz_clust->Fill(CFTC_r_cftc[i]); //z deos not exist
	}//end of loop on cft clusters

	for(Int_t i=0;i<16;i++){
	    if(doprint) cout <<" layer "<<i<<" nclus "<<nclus[i]<<" nstrip "<<nstrip[i]<<endl;
            //my_histos->H_cft_Nclus[i]->Fill(nclus[i]);
            my_histos->cft_nclus_l->Fill(i,nclus[i]);
            my_histos->cft_nstrip_l->Fill(i,nstrip[i]);	
}
	//CPS
	if(debug)cout <<" avant histos CPS"<<endl;

       	Float_t r;
        my_histos->	cps_nclus->Fill(CPSC_ncpsc);
        Float_t cps_etot=0.;
	for (Int_t i =0;i< CPSC_ncpsc;i++){
            cps_etot+=CPSC_e_cpsc[i];
            my_histos->H_cps_Eclus->Fill(CPSC_e_cpsc[i]);   
            my_histos->H_cps_xy_clust->Fill(CPSC_x_cpsc[i], CPSC_y_cpsc[i]);
	    r = sqrt(CPSC_x_cpsc[i]*CPSC_x_cpsc[i]+CPSC_y_cpsc[i]*CPSC_y_cpsc[i]);
            my_histos->cps_z_clust->Fill(CPSC_z_cpsc[i]);
	}//end of loop on CFT layers	
        my_histos->H_cps_etot->Fill(cps_etot);

	//FPS

	Int_t fpsn=0;
	Int_t fpss=0;
	if(debug)cout <<" avant histos FPS"<<endl;
	for (Int_t i =0;i< FPSC_nfpsc;i++){
	    if(FPSC_z_fpsc[i] >0.)fpsn++;
	    else  fpss++;
	}
	my_histos->H_fps_Nclus[0]->Fill(fpsn);
	my_histos->H_fps_Nclus[1]->Fill(fpss);
   
	//SMT
	if(debug)cout <<" before smt "<<endl;
	
	Int_t nst1=0;
	Int_t nst2=0;
	for (Int_t i =0;i< SMTB_nsmtb;i++){
	    nst1++;
	    //nst2+= SMTB_st2_smtb[i];
	    my_histos->H_smtb_xy_clust->Fill(SMTB_xg_smtb[i],SMTB_yg_smtb[i]);
	    r=sqrt(SMTB_xg_smtb[i]*SMTB_xg_smtb[i]+SMTB_yg_smtb[i]*SMTB_yg_smtb[i]);
	    my_histos->H_smtb_rz_clust->Fill( SMTB_zg_smtb[i],r);
	}//end of loop on SMT barrel
	my_histos->H_smtb_nclus[0]->Fill(nst1);
// my_histos->H_smtb_nclus[1]->Fill(nst2);
	if(debug)cout <<" before disks "<<endl;

	//disks
	nst1=0;
	nst2=0;
	for (Int_t i =0;i< SMTC_nsmtc;i++){
	    // cout <<" in loop on disks i "<<i<<" imax "<<SMTC_nsmtc<<endl;
	    r=sqrt(SMTC_xg_smtc[i]*SMTC_xg_smtc[i]+SMTC_yg_smtc[i]*SMTC_yg_smtc[i]);
            if(fabs( SMTC_zg_smtc[i])<55.){
		//   cout <<"f disk "<<endl;
		nst1++;
		my_histos->H_xy_clust_F->Fill(SMTC_xg_smtc[i],SMTC_yg_smtc[i]);
		my_histos->H_rz_clust_F->Fill(SMTC_zg_smtc[i],r);

            }else{
		//               cout <<"h disk "<<endl;
		nst2++;
		my_histos->H_xy_clust_H->Fill(SMTC_xg_smtc[i],SMTC_yg_smtc[i]);
		my_histos->H_rz_clust_H->Fill(SMTC_zg_smtc[i],r);
	    }
	}//end of loop on disks clusters
		my_histos->H_nclus_H->Fill(nst2);
	my_histos->H_nclus_F->Fill(nst1);
        if(doprint) cout <<" nb of disks "<<SMTC_nsmtc<<" dont F "<<nst1<<" dont H "<<nst2<<endl;
	if(debug)cout <<" before tracks "<<endl;
	//Tracks 
	my_histos->H_tra_nb->Fill(TRACK_ntra); //number of track per evt
	for (Int_t i =0;i< TRACK_ntra;i++){
	    my_histos->H_pt_tra->Fill(TRACK_pt_tra[i]);//Pt recons track
	    r = sqrt(TRACK_x_tra[i]*TRACK_x_tra[i]+TRACK_y_tra[i]*TRACK_y_tra[i]);
	    my_histos->H_rz_tra ->Fill(TRACK_z_tra[i],r);
	    my_histos->H_etaphi_tra->Fill( TRACK_eta_tra[i], TRACK_phi_tra[i]);//" theta phi recons. track"
	    my_histos->H_id_tra ->Fill(TRACK_id_tra[i]);//," id recons. track"
	    my_histos->H_cft_smt_tra->Fill( TRACK_smt_tra[i], TRACK_cft_tra[i]);// Nb. of CFt vs Nb of SMt on track"
 	}//end of loop on tracks

	//Vertices

        my_histos->H_vert_nb->Fill(VERTEX_nvert);
        if(doprint)cout <<" nvertex "<<VERTEX_nvert<<endl;
	for (Int_t i =0;i< VERTEX_nvert;i++){
	    //Int_t itr=VERTEX_gtr_ver[i];
	    //Int_t ntr=0;
	    //  Int_t tr [1000];
	    Int_t k=0;
	    //Int_t no_axial=1;
	    // if(doprint)cout <<" accum vertex "<<endl;
             Float_t r= sqrt( VERTEX_x_ver[i]*VERTEX_x_ver[i]+VERTEX_y_ver[i]*VERTEX_y_ver[i]);
             my_histos-> H_vert_nb_tra->Fill(VERTEX_ntrk_ver[i]);
             my_histos->H_vert_z->Fill(VERTEX_z_ver[i]);
             my_histos->H_vert_r->Fill(r);
             my_histos->H_vert_typ->Fill(VERTEX_ityp_ver[i]);
	     /*
	    while(itr!=0){
		if(itr%2 !=0){ 
		    //  tr[ntr]=k;
		    ntr++;
		    if(TRACK_cft_tra[k]<=8 && TRACK_smt_tra[k]==0)no_axial=0;
		    if(VERTEX_z_ver[i]>-.1 && VERTEX_z_ver[i]<.1 ){
		    }// end of test on Z vertex
		} 
		itr=itr/2;
		k++;
	    }//end of loop on tracks
	    if(no_axial)my_histos->Zv_no_axial->Fill(VERTEX_z_ver[i]);
	    else        my_histos->Zv_axial->Fill(VERTEX_z_ver[i]);
	     */
	}// end of loop on vertices

        // Muons
        if(debug)cout <<" nmuons "<<Muonid_Nmuonid<<endl;
        my_histos->H_nmu->Fill(Muonid_Nmuonid);
	if(Muonid_Nmuonid>=2){
	    Float_t spx = Muonid_Idpx[0]+Muonid_Idpx[1];	 
	    Float_t spy = Muonid_Idpy[0]+Muonid_Idpy[1];	 
	    Float_t spz = Muonid_Idpz[0]+Muonid_Idpz[1];	 
	    Float_t se  = Muonid_IdE[0]+Muonid_IdE[1];
	    my_histos->H_mass_dimul->Fill(sqrt(se*se-(spx*spx+spy*spy+spz*spz)));
	    my_histos->H_mass_dimuh->Fill(sqrt(se*se-(spx*spx+spy*spy+spz*spz)));
	}  

	for (Int_t i =0;i< Muonid_Nmuonid;i++){
	    Int_t nwa=Muonid_Idndeck[i]%10;
	    Int_t nwbc=(Muonid_Idndeck[i]/10)%10;
	    Int_t nsa=(Muonid_Idndeck[i]/1000)%10;
	    Int_t nsbc=Muonid_Idndeck[i]/10000;
	    Int_t is_nothing=1;
	    Int_t is_tight=0;
	    Int_t is_medium=0;
	    if(doprint)cout <<" accum muons "<<endl;
   my_histos-> H_nseg_mu ->Fill(Muonid_Idnseg[i]);
     my_histos->H_central_match_mu ->Fill(Muonid_Idncentralmatch[i]);
     my_histos->H_central_rank_mu->Fill(Muonid_Idcentralrank[i]);
     my_histos->H_theta_phi_mu->Fill(Muonid_Ideta[i],Muonid_Idphi[i]);
     my_histos->H_tight_mu->Fill(Muonid_tightmuo[i]);
     my_histos->H_pt_mu->Fill(Muonid_Idpt[i]);
     my_histos->e33_mu ->Fill(Muonid_e33[i]);
     my_histos-> e55_mu ->Fill(Muonid_e55[i]);

	    //Int_t is_loose=0; 
	    if(nwa >1 && nsa>0 && nwbc>1 && nsbc >0){
		is_nothing=0;
		if( nwbc>2 &&   Muonid_Idchisq[i] >0.)   is_tight=1;
		else    is_medium=1;                                 
	    } 
   
	    for (Int_t cas=0;cas<2;cas++){
		Int_t ik=cas;
		//		if(ik==1 && Muonid_tightmuo[i]==0) continue;
		//if(ik>=2 && Muonid_Nmuonid<2)      continue;
		//if(ik==3 &&  Muonid_tightmuo[i]==0)continue; 
		my_histos->PTH[cas]->Fill(Muonid_Idpt[i]);
		my_histos->ETA[cas]->Fill(Muonid_IdetaA[i]);
		my_histos->PHI[cas]->Fill(Muonid_Idphi[i]);
		my_histos->SEG[cas] ->Fill(Float_t(Muonid_Idnseg[i]));
		my_histos->TIGHT[cas] ->Fill(float(Muonid_tightmuo[i]));
		my_histos->CHI[cas] ->Fill(Muonid_Idchisqloc[i]);
		my_histos->TA[cas]  ->Fill(Muonid_IdsctimeA[i]);
		my_histos->TB[cas]  ->Fill(Muonid_IdsctimeC[i]);
		my_histos->DTBA[cas]->Fill(Muonid_IdsctimeC[i]-Muonid_IdsctimeA[i]);
		my_histos->E33[cas] ->Fill(Muonid_e33[i]);
		my_histos->E55[cas] ->Fill(Muonid_e55[i]);
		//ISOCA[cas] ->Fill(Muonid_IddrJet7[i]);
		//my_histos->ISOCB[cas] ->Fill(Muonid_IddrJet5[i]);
		//my_histos->DRJET5[cas] ->Fill(Muonid_IddrJet5[i]); 
		//my_histos->DRJET7[cas] ->Fill(Muonid_IddrJet7[i]);
		//my_histos->DRTRK[cas] ->Fill(Muonid_IddrTrk[i]);
		my_histos->		EIC1[cas] ->Fill(Muonid_EInCone1[i]); 
		my_histos->		EIC15[cas] ->Fill(Muonid_EInCone15[i]);
		//my_histos->EIC2[cas] ->Fill(Muonid_IdeInCone2[i]);
		//my_histos->EIC4[cas] ->Fill(Muonid_IdeInCone4[i]);
		//my_histos->EIC6[cas] ->Fill(Muonid_IdeInCone6[i]);
		//my_histos->EMTC[cas] ->Fill(Muonid_IdeMTC[i]);
		//my_histos->ETHALO[cas]->Fill(Muonid_IdetHalo[i]);
		//my_histos->ETKC5[cas]->Fill(Muonid_IdetTrkCone5[i]);
		//my_histos->ETKC7[cas]->Fill(Muonid_IdetTrkCone7[i]);
		//my_histos->EH4M1[cas]->Fill(Muonid_IdeInCone4[i]-Muonid_IdeInCone1[i]);
		my_histos->		SASBC[cas] ->Fill(nsa,nsbc);
		my_histos->		WAWBC[cas] ->Fill(nwa,nwbc);
		//my_histos->CALS[ik]  ->Fill(Muonid_IdcaleSig[i]);
		//  if(debug==1)cout<< " in accum position 6_n"<<endl;

	    }
	}//end of loop on muons

	//electrons
       my_histos->em_nb->Fill(EMPART_S_S_n_EM);
       for (Int_t i=0;i<EMPART_S_S_n_EM;i++){
         my_histos->em_id->Fill(EMPART_S_S_id[i]);
        if(EMPART_S_S_isol[i]<0.5)continue;
          my_histos->em_pt->Fill(EMPART_S_S_pT[i]);
          my_histos->em_frac->Fill(EMPART_S_S_EMfrac[i]);
          my_histos->em_iso->Fill(EMPART_S_S_iso[i]);
          my_histos->em_hmatx8->Fill(EMPART_S_S_HMx8[i]);
          my_histos->em_id_cut->Fill(EMPART_S_S_id[i]);
          my_histos->em_eta_phi->Fill(EMPART_S_S_eta[i],EMPART_S_S_phi[i]);
     }//end of loop on electrons
 	if(debug)cout <<"before jets"<<endl;

     //jets
	my_histos->	jet_nb->Fill(JCCB_JCCBnjet);
	for (Int_t i =0;i<  JCCB_JCCBnjet;i++){
	   my_histos->jet_et->Fill(JCCB_JCCBPT[i]);
	   my_histos->jet_phi_eta->Fill(JCCB_JCCBEta[i],JCCB_JCCBPhi[i]);
	   my_histos->jet_phi_eta_w->Fill(JCCB_JCCBEta[i],JCCB_JCCBPhi[i],JCCB_JCCBPT[i]);
	   Float_t f90=float(JCCB_JCCBn90[i])/float(JCCB_JCCBNitm[i]);
           my_histos->	    jet_f90->Fill(f90);
	}
	//Et miss
 	if(debug)cout <<"before Et miss"<<endl;

        my_histos-> etmis_x->Fill(ETMIS_etmis_x[0]);
        my_histos-> etmis_y->Fill(ETMIS_etmis_y[0]);
        my_histos-> etmis -> Fill(sqrt(ETMIS_etmis_x[0]*ETMIS_etmis_x[0]+ETMIS_etmis_y[0]*ETMIS_etmis_y[0]));
       my_histos-> scalar_et -> Fill(ETMIS_scal_et[0]);

	if(debug)cout <<"all detectors done"<<endl;
   }//end of loop on events
    //    my_histos->rec_hist_draw(nentries, nmax,psfile_name);
      my_histos->rec_hist_write(ntot,rootfile_name);

}