#define pmcsana_cxx
#include "pmcsana.h"
#include <TH2.h>
#include <TStyle.h>
#include <TFile.h>
#include <TTree.h>
#include "PCalTools.hpp"
#include "WZ_Utils.hpp"

#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

static double PDFWeightFactor[41];

double pmcsana::computePDF(int iset, double x, double q, int flav)
{
  double f[13];
  for (int i=0;i<13;i++) f[i]=0.;

  evolvepdfm_(iset,x,q,f);
  double pdf = f[flav+6];

  return pdf;
}


void pmcsana::Begin(TTree *tree) {

   TString option = GetOption();
   cout<<"Start pmcsana::Begin"<<endl;

   cout<<"Initialize some global parameters:"<<endl;
   _globalparameters = new GlobalParameters(_parameter_filename);

   // some parameters for photon merging
   // for photon merging probability
   _simulate_merge_fsr_photons_in_cone_prob = _globalparameters->simulate_merge_fsr_photons_in_cone_prob();
   _merge_fsr_photons_in_cone_prob_p0 = _globalparameters->merge_fsr_photons_in_cone_prob_p0();
   _merge_fsr_photons_in_cone_prob_p1 = _globalparameters->merge_fsr_photons_in_cone_prob_p1();
   _merge_fsr_photons_in_cone_prob_p2 = _globalparameters->merge_fsr_photons_in_cone_prob_p2();
   _merge_fsr_photons_in_cone_prob_p3 = _globalparameters->merge_fsr_photons_in_cone_prob_p3();
   _simulate_merge_fsr_photons_out_cone_prob = _globalparameters->simulate_merge_fsr_photons_out_cone_prob();
   _merge_fsr_photons_out_cone_prob_p0 = _globalparameters->merge_fsr_photons_out_cone_prob_p0();
   _merge_fsr_photons_out_cone_prob_p1 = _globalparameters->merge_fsr_photons_out_cone_prob_p1();
   _merge_fsr_photons_out_cone_prob_p2 = _globalparameters->merge_fsr_photons_out_cone_prob_p2();
   _merge_fsr_photons_out_cone_prob_p3 = _globalparameters->merge_fsr_photons_out_cone_prob_p3();
   _merge_fsr_photons_out_cone_prob_p4 = _globalparameters->merge_fsr_photons_out_cone_prob_p4();
   // for merged photon response
   _simulate_merge_fsr_photons_in_cone_response = _globalparameters->simulate_merge_fsr_photons_in_cone_response();
   _merge_fsr_photons_in_cone_response_p0 = _globalparameters->merge_fsr_photons_in_cone_response_p0();
   _merge_fsr_photons_in_cone_response_p1 = _globalparameters->merge_fsr_photons_in_cone_response_p1();
   _merge_fsr_photons_in_cone_response_p2 = _globalparameters->merge_fsr_photons_in_cone_response_p2();
   _merge_fsr_photons_in_cone_response_p3 = _globalparameters->merge_fsr_photons_in_cone_response_p3();
   _merge_fsr_photons_in_cone_response_p4 = _globalparameters->merge_fsr_photons_in_cone_response_p4();
   _merge_fsr_photons_in_cone_response_p5 = _globalparameters->merge_fsr_photons_in_cone_response_p5();
   _simulate_merge_fsr_photons_out_cone_response = _globalparameters->simulate_merge_fsr_photons_out_cone_response();
   _merge_fsr_photons_out_cone_response_p0 = _globalparameters->merge_fsr_photons_out_cone_response_p0();
   _merge_fsr_photons_out_cone_response_p1 = _globalparameters->merge_fsr_photons_out_cone_response_p1();
   _merge_fsr_photons_out_cone_response_p2 = _globalparameters->merge_fsr_photons_out_cone_response_p2();
   _merge_fsr_photons_out_cone_response_p3 = _globalparameters->merge_fsr_photons_out_cone_response_p3();
   _merge_fsr_photons_out_cone_response_p4 = _globalparameters->merge_fsr_photons_out_cone_response_p4();
   _merge_fsr_photons_out_cone_response_p5 = _globalparameters->merge_fsr_photons_out_cone_response_p5();

   // if we want to do PDF reweight
   // will calculate weight for all 41 PDFs, hardcoded
   _doPdfReweight = _globalparameters->pdfreweight();
   if(_doPdfReweight) {
     
     // pdf path
     string pdfpath = "/d0usr/products/lhapdfsets/NULL/v3_0/";
     
     char _nameini[100];
     string tmp = pdfpath;
     tmp += _globalparameters->pdfori();
     strcpy(_nameini, tmp.c_str());
     cout<<"Files are generated with PDF: "<<_nameini<<endl;
     
     int set1 = 1;
     int pdfori_subset = _globalparameters->pdfori_subset();
     vector<string> vpdfnew;
     vector<int> vpdfnew_subset;
     for (int i=0; i<41; ++i) {
       vpdfnew.push_back("cteq61.LHgrid");
       vpdfnew_subset.push_back(i);
     }   
     
     initpdfsetm_(set1,*_nameini);
     initpdfm_(set1, pdfori_subset);
     
     // Initialize the new PDF sets :
     cout << "loading PDFs for reweighting" << endl;
     for (int i=0; i<vpdfnew.size();++i) {     
       char namenew[160];
       string tmp2 = pdfpath;
       tmp2 += vpdfnew[i];
       strcpy(namenew, tmp2.c_str());
       
       int iset = i+2;     
       initpdfsetm_(iset, *namenew);
       initpdfm_(iset, vpdfnew_subset[i]);
       cout << "Reading "<<namenew << " subset " << vpdfnew_subset[i] << endl;     
     }
   } // initialize all PDFs
   
   // initialize 
   if(_runoption == 0) {
     cout<<"Initialize ZAnalysis"<<endl;
     _zanalysis = new ZAnalysis(_parameter_filename);
   } else if(_runoption == 1) {
     cout<<"Initialize WAnalysis"<<endl;
     _wanalysis = new WAnalysis(_parameter_filename);
   } else if(_runoption == 2) {
     cout<<"Initialize ZMuAnalysis"<<endl;
     _zmuanalysis = new ZMuAnalysis(_parameter_filename);
   } else if(_runoption == 3) {
     cout<<"Initialize JPsiMuAnalysis"<<endl;
     _jpsimuanalysis = new JPsiMuAnalysis(_parameter_filename);
   } else if(_runoption == 4) {
     cout<<"Initialize ZNuNuAnalysis"<<endl;
     _znunuanalysis = new ZNuNuAnalysis(_parameter_filename);
     // also need to read parameters for ZNuNuAnalysis to count Loos and Tight neutrinos
     PParameterReader parm(_parameter_filename.Data());  
     _NuLoose_EtaMax = parm.GetDouble("ZNuNu_NuLoose_EtaMax", 1000.0);  // no cut
     _NuLoose_PtMin  = parm.GetDouble("ZNuNu_NuLoose_PtMin",  0.0);     // no cut
     _NuTight_EtaMax = parm.GetDouble("ZNuNu_NuTight_EtaMax", 1000.0);  // no cut
     _NuTight_PtMin  = parm.GetDouble("ZNuNu_NuTight_PtMin",  0.0);     // no cut
     if(( _NuTight_PtMin < _NuLoose_PtMin ) || 
	( _NuTight_EtaMax > _NuLoose_EtaMax )) {
       cerr << "pmcsana: Wrong specification of Loose and Tight cuts!!!" << endl
	    << "         pT(Loose nu)>=" << _NuLoose_PtMin 
	    << " GeV, |eta(Loose nu)|<" << _NuLoose_EtaMax << endl
	    << "         pT(Tight nu)>=" << _NuTight_PtMin 
	    << " GeV, |eta(Tight nu)|<" << _NuTight_EtaMax << endl;
       throw;
     }
   } 

   // initialize random number generator
   if (_randomseed==0) {
     _dummy = new TRandom3(_globalparameters->rand_seed());
     cout << "pmcsana: Using random seed = " << _globalparameters->rand_seed() << endl;
   }
   else {
     _dummy = new TRandom3(_randomseed);
     cout << "pmcsana: Using random seed = " << _randomseed << endl;
   }

   // also set gRandom points to the instance of TRandom3 
   gRandom = _dummy;
 
   // initialize beam spot 
   if(_globalparameters->VtxSmear_Option()==1) {
     cout<<"Initializing beamspot dependence on the run ranges and instantaneous luminosity now ..."<<endl;
     _beam = new BeamWeight((string)(_globalparameters->BeamSpotShape_file()));
     cout<<"Finished beamspot"<<endl;
   }

   // read instantaneous luminosity profile
   _instlumi_profile = _globalparameters->instlumi_profile();

   cout<<"End pmcsana::Begin"<<endl;
}

void pmcsana::SlaveBegin(TTree *tree) {
  Init(tree);
  TString option = GetOption();
}

pmcsana::~pmcsana() {
  if(_zanalysis!=NULL) delete _zanalysis;
  if(_wanalysis!=NULL) delete _wanalysis;
  if(_zmuanalysis!=NULL) delete _zmuanalysis;
  if(_jpsimuanalysis!=NULL) delete _jpsimuanalysis;
  if(_znunuanalysis!=NULL) delete _znunuanalysis;
  if(_dummy!=NULL) delete _dummy;
}

// process event
Bool_t pmcsana::Process(Long64_t entry) {

  //skip this part if doing ZAnalysis in a MINUIT mode with event selection disabled
  if(_runoption==0 && _zanalysis->MINUIT_skipPreselection()) return false;
 
  //how often do we print out the number of events processed?
  static int update_frequency = _globalparameters->update_frequency();
  static int eventsProcessed = 0; 
  if(eventsProcessed%update_frequency==0) cout<<"Processing event number: "<< eventsProcessed 
				    << " / " << _numEvents << endl;
  eventsProcessed++;
  
  //  
  // construct PMCS event information and put all information into PMCSEvent class
  //
  double evtweight;
  PMCSParticle wzboson_gen;
  vector<PMCSEMObj> emobjs_gen;
  vector<PMCSMuon> muons_gen;
  PMCSMet met_gen;
  PMCSRecoil recoil_gen;
  PMCSVtx vtx_gen;

  // run number and instantaneous luminosity for this zerobias event
  // also keep the index so that we will use the same event for overlay
  int runNo = 1, ZBEvtIndex = -1, MBEvtIndex = -1;
  double instlumi = 0.;
  
  //
  // get run number and instantaneous luminosity from zerobias overlay
  // also keep the indices for the overlaid zerobias and minbias events
  // only need for ZAnalysis and WAnalysis, not relevant for ZMuAnalysis and JPsiMuAnalysis
  // if we read recoil file, then we do not need to get the indices to the overlaid Minbias and Zerobias events
  //
  if(_runoption == 0) {
    if(!_zanalysis->getRecoilSmear()->read_recoil_file()){
      _zanalysis->getRecoilSmear()->setRunNoInstLumiMBZBEvtIndex(_dummy); 
      _zanalysis->getRecoilSmear()->getRunNoInstLumiMBZBEvtIndex(runNo, instlumi, MBEvtIndex, ZBEvtIndex);
    }
    
    // if we do not read the mb and zb library files, will get the luminosity from a histogram
    if(_zanalysis->getRecoilSmear()->METSmear_Option()!=3 || _zanalysis->getRecoilSmear()->read_recoil_file())  
      instlumi =_instlumi_profile->GetRandom();
  } else if(_runoption == 1) {
    if(!_wanalysis->getRecoilSmear()->read_recoil_file()){
      _wanalysis->getRecoilSmear()->setRunNoInstLumiMBZBEvtIndex(_dummy);  
      _wanalysis->getRecoilSmear()->getRunNoInstLumiMBZBEvtIndex(runNo, instlumi, MBEvtIndex, ZBEvtIndex); 
    }

    // if we do not read the mb and zb library files, will get the luminosity from a histogram
    if(_wanalysis->getRecoilSmear()->METSmear_Option()!=3 || _wanalysis->getRecoilSmear()->read_recoil_file())  
      instlumi =_instlumi_profile->GetRandom();
  }

  // read all generator level information from pmcs_ana branch
  constructPMCSEvent(entry, runNo, instlumi, evtweight, wzboson_gen, emobjs_gen, 
		     met_gen, recoil_gen, vtx_gen, muons_gen, _dummy);    

  // construct pmcsEvent for future use
  PMCSEvent pmcsevent(runNo, MBEvtIndex, ZBEvtIndex, evtweight, instlumi, wzboson_gen, emobjs_gen, 
		      met_gen, recoil_gen, vtx_gen, muons_gen);

  // whether we want to do vertex cut
  if(_globalparameters->do_zvtx_cut() && 
     fabs(vtx_gen.vtxz()) > _globalparameters->zvtx_cut()) {  
    return kFALSE;               // skip this event
  } else{                        // analyze this event
    if(_runoption == 0) _zanalysis->analyzeEvent(pmcsevent, _dummy);  
    if(_runoption == 1) _wanalysis->analyzeEvent(pmcsevent, _dummy);    
    if(_runoption == 2) _zmuanalysis->analyzeEvent(pmcsevent, _dummy);    
    if(_runoption == 3) _jpsimuanalysis->analyzeEvent(pmcsevent, _dummy);        
    if(_runoption == 4) {
      // used for Z->nunu analysis, need to apply cuts on generator level neutrinos
      b_pmcs_ana->GetEntry(entry);
      int nLooseNeutrinos = 0, nTightNeutrinos = 0;
      for(int ipart=0; ipart<pmcs_ana_npart; ipart++) {
        double pid = pmcs_ana_ppid[ipart];
        if(abs(pid)==12 || abs(pid)==14 || abs(pid)==16) {
          double eta = pmcs_ana_peta[ipart];
	  double pt  = pmcs_ana_ppt[ipart];
	  if( fabs(eta)<_NuLoose_EtaMax && pt>_NuLoose_PtMin) nLooseNeutrinos++;
	  if( fabs(eta)<_NuTight_EtaMax && pt>_NuTight_PtMin) nTightNeutrinos++;
        }
      } // loop over all particles to find neutrinos
      _znunuanalysis->analyzeEvent(pmcsevent, nLooseNeutrinos, nTightNeutrinos, _dummy);        
    }
    return kTRUE;
  }

}

// SlaveTerminate
void pmcsana::SlaveTerminate() {}

// Terminate
void pmcsana::Terminate() {
  if(_runoption == 0) _zanalysis -> jobSummary(_dummy, _numEvents);  
  if(_runoption == 1) _wanalysis -> jobSummary();
  if(_runoption == 2) _zmuanalysis -> jobSummary();  
  if(_runoption == 3) _jpsimuanalysis -> jobSummary();
  if(_runoption == 4) _znunuanalysis -> jobSummary();
}

// construct PMCSEvent
void pmcsana::constructPMCSEvent(Int_t entry, int runNo, double instlumi, double& evtweight, PMCSParticle& wzboson, vector<PMCSEMObj>& emobjs,
				 PMCSMet& met, PMCSRecoil& recoil, PMCSVtx& vtx, vector<PMCSMuon>& muons, TRandom3 *dummy) {
  b_pmcs_ana->GetEntry(entry);
  b_pmcs_vtx->GetEntry(entry);

  // Event Weight, if we do not want to read event weight, set all of them to 1
  evtweight = pmcs_ana_evwt[0];
  if(!(_globalparameters->reweight())) evtweight = 1.;

/*
  // calculate PDF weight factor
  if(_doPdfReweight) {
    double x1 = pmcs_ana_evx1[0];
    double x2 = pmcs_ana_evx2[0];
    int flav1 = (int)(pmcs_ana_evflav1[0]);
    int flav2 = (int)(pmcs_ana_evflav2[0]);

    // need to convert flav2 to -flav2 since flav2 is coming from antiproton
    // while computerPDF gets values from proton
    flav2 = -flav2; 
    double q = sqrt(pmcs_ana_evqsq[0]);
    
    int iset_pdf = 1;

    double rpdfini1 = computePDF(iset_pdf, x1, q, flav1);
    double rpdfini2 = computePDF(iset_pdf, x2, q, flav2);

    for(int i=0; i<41; ++i) {
      iset_pdf = i+2;

      double rpdfnew1 = computePDF(iset_pdf, x1, q, flav1);
      double rpdfnew2 = computePDF(iset_pdf, x2, q, flav2);

      double rpdfweight = 1.;
      if(rpdfini1!=0 && rpdfnew1!=0) rpdfweight *= rpdfnew1/rpdfini1;
      if(rpdfini2!=0 && rpdfnew2!=0) rpdfweight *= rpdfnew2/rpdfini2;

      // fill the static variables
      PDFWeightFactor[i] = rpdfweight;

    } // 41 PDF sets
    cout<<endl;
  } // PDF reweight
*/

  // primary vertex
  vtx = PMCSVtx(0., 0., pmcs_vtx_vtxzg[0]);

  // need to smear vertex here because I need to 
  // calculate electron detector eta and phi
  bool smearVertex = _globalparameters->smearVertex();
  double vtx_res = _globalparameters->VtxResolution();
  int vtxSmearOption = _globalparameters->VtxSmear_Option();
  if(smearVertex) {
    if(vtxSmearOption == 0) vtx.SmearVtx(vtx_res, dummy);
    else {
     // first smear primary vertex using a 40 cm gaussian function
     vtx.SmearVtx(40., dummy);
     // then calculate the event weight that is relative to this 40 cm gaussian
     double z = vtx.vtxz();
     double beamshape = _beam->beamshape(runNo, instlumi, z);
     double gaus = _beam->Gaussian(40, 0., z);

     // now we modify event weight
     evtweight *= beamshape/gaus;
    } 
  }

  // Merge electron and photons
  if( _globalparameters->merge_fsr_photons() ) 
    fixEMBlock(entry, vtx, _globalparameters->merge_fsr_radius_CC(), _globalparameters->merge_fsr_radius_EC());
  else 
    fixEMBlock(entry, vtx, -1.0, -1.0);

  // clear emobjs vector first  
  emobjs.clear();

  double metx_raw = 0., mety_raw = 0., metz_raw = 0., met_E = 0.;

  // loop over all generator level particles
  for(int ipart=0; ipart<pmcs_ana_npart; ipart++) {
    double pid = pmcs_ana_ppid[ipart];

    // electrons and photons are treated differently since
    // we have to do electron-photon merging
    double eta = pmcs_ana_peta[ipart];
    double phi = pmcs_ana_pphi[ipart];
    double energy = pmcs_ana_pE[ipart];
    double pt = pmcs_ana_ppt[ipart];
    double ppx = pmcs_ana_ppx[ipart];
    double ppy = pmcs_ana_ppy[ipart];
    double ppz = pmcs_ana_ppz[ipart];
     
    if(abs(pid)==24 || abs(pid)==23) 
      wzboson = PMCSParticle(pid, energy, eta, phi, ppx, ppy, ppz);

    if(abs(pid)==13) {
      PMCSMuon muon(pid, energy, eta, phi);
      muons.push_back(muon);
    } //muons

    // calculate raw missing Et using neutrino information
    if(abs(pid)==12 || abs(pid)==14 || abs(pid)==16) {
      metx_raw += ppx;
      mety_raw += ppy;
      metz_raw += ppz;
      met_E    += energy;
     
    }
  } // loop over all particles

  // Missing energy
  met = PMCSMet(metx_raw, mety_raw);
  
  // for electrons, track pT and real pT may be different since 
  // electron and photons may merged (pmcs_em_elfid, pmcs_em_elpts are two unused variables in pmcs_em, 
  // I just use them to pass electron ID and track pT information)
  for(int j=0;j<pmcs_em_nelg;j++){ 
    if(pmcs_em_elptg[j]>_globalparameters->pT_Cut()) {
      PMCSEMObj emobj(pmcs_em_elfid[j], pmcs_em_eleg[j], pmcs_em_eletag[j], 
		      pmcs_em_elphig[j], pmcs_em_elpts[j]);
      emobjs.push_back(emobj);
    }//pt cut
  }//electrons

  // electrons and photons are already merged here since I use pmcs_em branch information
  // for photons, track pT is the same as photon pTset to 0
  // so that later we will not calculate E/p for EMobjs that are really photons
  // all photons are treated as EMObjects
  for(int j=0;j<pmcs_em_nphg;j++){ 
    //    if(pmcs_em_phptg[j]>_globalparameters->pT_Cut()) {
    PMCSEMObj emobj(22, pmcs_em_pheg[j], pmcs_em_phetag[j], pmcs_em_phphig[j], 0.);
    emobjs.push_back(emobj);
    //    }//pt cut
  }//photons
 
  // mainly used for Sahal's ratio method to pass generator level neutrino information
  if(_runoption == 1){
    TLorentzVector e(metx_raw, mety_raw,metz_raw, met_E);
    PMCSEMObj Nugen(12,e.E(),e.Eta(),e.Phi(),e.Pt());
    emobjs.push_back(Nugen);
  }
 
  // have to loop over all EMObjs and calculate detector eta, phi in global and local coordinate systems  
  for(int iem=0; iem<emobjs.size(); iem++) {
    emobjs[iem].SetDetEtaDetPhi(vtx);
    emobjs[iem].SetDetEtaDetPhiGlobal(vtx);
  }

  // sort generator level electrons according to pT
  sort(emobjs.begin(), emobjs.end(), ComparePt);

  // recoil jet, same as WZ boson, except opposite direction
  double recoilx_system = -wzboson.ppx();
  double recoily_system = -wzboson.ppy();

  // the following code is commented out on 04/23/2007
  // we decided to treat fsr photons differently than the raw recoil
  // photons will always be assumed as an EM cluster, the energy response will be different
  // no UE simulation below electron window, and the photons will be added back to the smeared recoil system later

  // for photons that have very low pT and well separated from electrons
  // they should be treated as part of the recoil system 
  //  if(_globalparameters->merge_fsr_photons_recoil()) {
  //    for(int iph=0; iph<pmcs_em_nphg; iph++) {
  //      if(pmcs_em_phptg[iph] < _globalparameters->pT_Cut()) {
  //	recoilx_system += pmcs_em_phptg[iph] * cos(pmcs_em_phphig[iph]); 
  //	recoily_system += pmcs_em_phptg[iph] * sin(pmcs_em_phphig[iph]); 
  //      }
  //    } // loop over all photons
  //  } 
  
  // create vector for the recoil system
  recoil = PMCSRecoil(recoilx_system, recoily_system);

  //
  // print out information
  //
  if(_globalparameters->debug()) {
    cout<<"raw event weight: "<<evtweight<<endl;
    cout<<"    event vertex: ";
    vtx.Print();
    cout<<"    boson:        ";
    wzboson.Print();
    cout<<"    electrons:    ";
    for(int iem=0; iem<emobjs.size(); iem++) emobjs[iem].Print();
    cout<<"    muons:    ";
    for(int imu=0; imu<muons.size(); imu++) muons[imu].Print();
    cout<<"    event recoil: ";
    recoil.Print();
    cout<<"    event met:    ";
    met.Print();
  }
}

// merge electron and photon if they are close to each other
void pmcsana::fixEMBlock(Int_t entry, PMCSVtx& vtx, double cut_dR_CC, double cut_dR_EC) {
  //Recreate pmcs_em from the pmcs_ana block.
  //Include FSR photon merging.
  
  // Shortcomings:
  //   * the electron-photon matching is too simple. it is possible
  //     that a photon is not assigned to the closest electron. i
  //     think the effect is incredibly tiny, though.
  //   * photon-electron merging is by four-vector addition. This should
  //     probably be changed to something resembling the EM clustering
  //     algorithm in the data.
  // 
  b_pmcs_ana->GetEntry(entry);

  // get primary vertex
  double vtxZ = vtx.vtxz();
  
  const Bool_t FDBG = kFALSE;

  Int_t nphi=0, neli=0;
  Int_t ph_indices[50];
  Int_t el_indices[50];
  for(int i=0; i<50; i++) {
    ph_indices[i]=-1;
    el_indices[i]=-1;
  }

  if(FDBG) printf("\n");
  
  //index all photons and electrons in the pmcs_ana block
  //if there are many particles in this event, this will save some time
  for(int ip = 0; ip < pmcs_ana_npart; ip++) {    
    if(TMath::Abs(pmcs_ana_ppid[ip])==22 && pmcs_ana_pE[ip]>0.001){ //photon      
      TLorentzVector ph(pmcs_ana_ppx[ip],pmcs_ana_ppy[ip], pmcs_ana_ppz[ip],pmcs_ana_pE[ip]);      
      if(FDBG) printf("Original photon: E=%.3f pT=%.3f eta=%0.3f phi=%0.3f ",
                      ph.E(),ph.Pt(),ph.Eta(),ph.Phi());
      
      ph_indices[nphi]=ip;
      nphi++;
      if(FDBG) cout<<endl;
    } else if(TMath::Abs(pmcs_ana_ppid[ip])==11){ //electron
      el_indices[neli]=ip;
      neli++;
    } 
  }//loop over all generator-level particles
  
  pmcs_em_nelg=0;
  
  //loop over electrons
  for(int ke=0;ke<neli;ke++){
    Int_t ie = el_indices[ke];
    
    TLorentzVector e(pmcs_ana_ppx[ie],pmcs_ana_ppy[ie], pmcs_ana_ppz[ie],pmcs_ana_pE[ie]);
    TLorentzVector e_before(e);
    if(FDBG) printf("Electron: E=%.2f pT=%.2f eta=%0.2f phi=%0.2f\n",
                    e.E(),e.Pt(),e.Eta(),e.Phi());
    
    // track pT is the raw electron pT before it is merged with photons
    double trkpT = e.Pt();
    
    // merge closest photons
    for(int j=0;j<nphi;j++){ 
      if(ph_indices[j] < 0) continue; //unused entry in ph_indices
      Int_t ip = ph_indices[j];
      
      // use a simple cone (deltaR in CC and fixed radius in EC) to merge electrons and photons
      TLorentzVector ph(pmcs_ana_ppx[ip],pmcs_ana_ppy[ip], pmcs_ana_ppz[ip],pmcs_ana_pE[ip]);
      TVector3 calp_e=PCalTools::floorXYZ(vtxZ, e.Eta(), e.Phi());
      TVector3 calp_ph=PCalTools::floorXYZ(vtxZ, ph.Eta(), ph.Phi());
      bool isCC = false;
      if(fabs(calp_e.Eta())<1.3) isCC = true;
      
      if(isCC && FDBG) printf("    ph: eta=%.3f phi=%.3f E=%.3f R_CC=%.3f",
                              ph.Eta(), ph.Phi(), ph.E(), ph.DeltaR(e));
      
      if(!isCC && FDBG) printf("    ph: eta=%.3f phi=%.3f E=%.3f R_EC=%.3f",
                               ph.Eta(), ph.Phi(), ph.E(), PCalTools::Deltar(calp_e, calp_ph));

      // decide whether to merge photons with the electron or not      
      bool merge = false;
      if(isCC) merge = (e.DeltaR(ph) < cut_dR_CC);
      else merge = (PCalTools::Deltar(calp_e, calp_ph) < cut_dR_EC);
     
      //
      // assume CC and EC photons have the same probability to reach the calorimetr
      // and the same response
      //
      if(merge) {
	// probability for photons to reach the calorimeter
	double prob_photon = 1.;
	if(_simulate_merge_fsr_photons_in_cone_prob)      
	  prob_photon = wz_utils::photon_prob_in_cone(ph.Pt(), _merge_fsr_photons_in_cone_prob_p0, 
						      _merge_fsr_photons_in_cone_prob_p1, _merge_fsr_photons_in_cone_prob_p2, 
						      _merge_fsr_photons_in_cone_prob_p3);
	
	// remove the photons that are not supposed to reach the calorimeter permanently
	double prob = _dummy->Uniform(0., 1.);
	if(prob < prob_photon) {	  

	  // calculate the fraction of photon energy merged with the EM cluster
	  double fraction = 1.;
	  if(_simulate_merge_fsr_photons_in_cone_response) 
	    fraction = wz_utils::photon_response_in_cone(ph.Eta(), ph.Pt(), _merge_fsr_photons_in_cone_response_p0,
							 _merge_fsr_photons_in_cone_response_p1, _merge_fsr_photons_in_cone_response_p2,
							 _merge_fsr_photons_in_cone_response_p3, _merge_fsr_photons_in_cone_response_p4,
							 _merge_fsr_photons_in_cone_response_p5);

	  if(FDBG) printf(" ==> MERGED energy fraction=%.3f ", fraction);
	  ph_indices[j]=-1; //mark as found, remove this partially merged photons, is this the best way ???
	  e += fraction * ph;  //Merge. Is this the best way?
	} // only a fraction of photons can reach the calorimeter
      } // merged

      if(FDBG) printf("\n");
    }
    
    //Create pmcs_em object (generated quantities only)
    if(FDBG){
      printf("    pmcs_em entry: %.2f->%.2f  Pt: %.2f->%.2f  "
             "Eta: %.2f->%.2f  Phi: %.2f->%.2f\n",
             e_before.E(),e.E(),
             e_before.Pt(),e.Pt(),
             e_before.Eta(),e.Eta(),
             e_before.Phi()>0.0?e_before.Phi():e_before.Phi()+2*TMath::Pi(),
             e.Phi()>0.0?e.Phi():e.Phi()+2*TMath::Pi());
    }
    
    pmcs_em_eleg[ke]=e.E();
    pmcs_em_elptg[ke]=e.Pt();
    pmcs_em_eletag[ke]=e.Eta();
    pmcs_em_elphig[ke]=e.Phi();
    if(pmcs_em_elphig[ke]<0.0) pmcs_em_elphig[ke] += 2*TMath::Pi();
    
    // need to keep track pT and particle ID
    // currently choose to use two unused variables pmcs_em_elpts, pmcs_em_elfid
    // which means these two variables can not be used for other purposes
    
    pmcs_em_elpts[ke] = trkpT;                // track pT before electron/photon merging
    pmcs_em_elfid[ke]  = pmcs_ana_ppid[ie];   //use this as pid
    
    pmcs_em_nelg++;
  }
  
  //Recreate the list of photons, excluding those that are collinear
  //to electrons. 
  int NPhotonAfterMerging = 0;
  
  for(int j=0;j<nphi;j++) {
    if(ph_indices[j] < 0) continue;    
    Int_t iph=ph_indices[j];
    
    TLorentzVector photon(pmcs_ana_ppx[iph],pmcs_ana_ppy[iph],
                          pmcs_ana_ppz[iph],pmcs_ana_pE[iph]);
    
    // require photon with pT at least 1 MeV 
    if(ph_indices[j]!=-1 && photon.Pt()>0.001){

      // these photons will be included in the recoil system
      // and I will apply the fraction of photons that can reach the calorimeter
      // and the fraction of energy deposited in the calorimeter
      double prob_photon = 1.;
      if(_simulate_merge_fsr_photons_out_cone_prob)      // assume CC and EC photons have the same probability
	prob_photon = wz_utils::photon_prob_out_cone(photon.Pt(), _merge_fsr_photons_out_cone_prob_p0, 
						     _merge_fsr_photons_out_cone_prob_p1, _merge_fsr_photons_out_cone_prob_p2, 
						     _merge_fsr_photons_out_cone_prob_p3, _merge_fsr_photons_out_cone_prob_p4);
      
      // remove the photons that are not supposed to reach the calorimeter permanently
      double prob = _dummy->Uniform(0., 1.);
      if(prob < prob_photon) {	  
	// calculate the fraction of photon energy merged with the EM cluster
	double fraction = 1.;
	if(_simulate_merge_fsr_photons_out_cone_response) 
	  fraction = wz_utils::photon_response_out_cone(photon.Eta(), photon.Pt(), _merge_fsr_photons_out_cone_response_p0,
						       _merge_fsr_photons_out_cone_response_p1, _merge_fsr_photons_out_cone_response_p2,
						       _merge_fsr_photons_out_cone_response_p3, _merge_fsr_photons_out_cone_response_p4,
						       _merge_fsr_photons_out_cone_response_p5);
		
	if(FDBG) printf(" ==> MERGED energy fraction=%.3f ", fraction);
	ph_indices[j]=-1; //mark as found, remove this partially merged photons, is this the best way ???
	photon = fraction * photon;  // change photon energy that will be included in the recoil system
	      
	pmcs_em_phptg[NPhotonAfterMerging]  = photon.Pt();
	pmcs_em_phetag[NPhotonAfterMerging] = photon.Eta();
	pmcs_em_phphig[NPhotonAfterMerging] = photon.Phi();
	pmcs_em_pheg[NPhotonAfterMerging]   = photon.E();
	pmcs_em_phfid[NPhotonAfterMerging]  = 0; // this variable does not matter
	
	NPhotonAfterMerging += 1;	
      } // only a fraction of photon can reach the calorimeter
      
      if(NPhotonAfterMerging>=120){
        printf("Warning: Too many photons (%d)!\n",nphi);
        break;
      }
    }//get rid of merged photons and photon with pT at least 1 MeV
  } // loop over all photons
  
  pmcs_em_nphg = NPhotonAfterMerging;
  
  return;
}