#include "caf_util/ElectronSelector.hpp"
#include "caf_util/TRefFinder.hpp"
#include "cafe/Config.hpp"
#include "cafe/DefaultFileFinder.hpp"
#include "caf_util/WZ_Utils.hpp"
#include "caf_util/CaloGeometryManager.hpp"

#include <iostream>
#include <sstream>

struct xyz{double x; double y; double z;};
std::map<long, xyz> CELLtable;
const int minEta=-37;
const int maxEta=37;
const int minPhi=1;
const int maxPhi=64;
const int minLayer=1;
const int maxLayer=17;

using namespace cafe ;
//using namespace edm ;
using namespace std ;

namespace {
  bool moreThan(const TMBEMCluster& a, const TMBEMCluster& b) {
    return a.Pt() > b.Pt();    
  }
}

namespace caf_util {

  ElectronSelector::ElectronSelector(const char *name) : 
    cafe::SelectUserObjects<TMBEMCluster>(name), _nselected(0), 
    _event(0), _finder(0), _emidfilepath(""), _doNewFiducial(false),
    _treehandler(0)
  {
    cafe::Config config(name);
    
    _quality = config.get("Quality", "Tight");
    _ver = config.get("Version", 1);

    _fiducial = config.get("FiducialCutSet","LEGACY");
    string fiducialCutsPath = "";
    string fiducialDefsPath = "";
    if( _fiducial != "LEGACY") {
      _doNewFiducial = true;
      _fiducialVersion = config.get("FiducialCutSetVersion",-1);
      fiducialCutsPath = config.get("FiducialCutsFile","fiducial_cuts.txt");
      fiducialDefsPath = config.get("FiducialDefsFile","fiducial_defs.txt");
      _fidCutsToApply = config.getVString("FiducialCutsToApply"," ");
      if( _fidCutsToApply.size() == 0 )
	_fidCutsToApply.push_back("all");
    }

    _pt       = config.get("pT",-1.0);
    _ptmax = config.get("pTmax",-1.0);
    _etaCC    = config.get("etaCC",-1.0);
    _etaECmin = config.get("etaECmin",-1.0);
    _etaECmax = config.get("etaECmax",-1.0);

    std::string toBranch = config.get("To","");
    _toTreeName = config.get("Tree","");
    _leadingBranchName = config.get("LeadingBranchName",("Leading"+toBranch).c_str());
    _createLeadingBranch = config.get("CreateLeadingBranch",false);

    _nelectrons = config.get("nElectrons", 0);
    _nelectronsmax = config.get("nElectronsMax", -1);
    _sort = config.get("Sort", 0);
    string emidfilepath = config.get("SelectionFile","emid_cuts/support/emid_cuts.txt");
    
    DefaultFileFinder ffinder ;
    ffinder.push_front("./emid_cuts/support") ;

    if (emidfilepath != "" )
      _emidfilepath = ffinder.findFile(emidfilepath) ;

    if ( _emidfilepath == "") {
      err() << "ElectronSelector(" << name << ")\nError: can not find the file emid_cuts with name [" << emidfilepath << "]. " 
	    << "Please specify the correct path in your config file:\n   "
	    << name << ".SelectionFile: path/emid_cuts.txt\n";
      throw runtime_error("Can not find emid_cuts file");
    }
  
    if ( _createLeadingBranch ) {
	_treehandler = new TreeHandler();
	_sort = true;
    }

    out() << "\n======== ElectronSelector("<<name<<") ==================\n\n";
    out() << "Using ElectronID criteria from ["<<_emidfilepath<<"]\n\n";
    _id = new emID(_emidfilepath, _quality, _ver, 0);

    if( _doNewFiducial ) {
      string cutsfile="";
      string defsfile="";
      if( fiducialCutsPath != "")
	cutsfile = ffinder.findFile(fiducialCutsPath);
      if( fiducialDefsPath != "")
	defsfile = ffinder.findFile(fiducialDefsPath);
      if( cutsfile == "") {
	err() << "ElectronSelector(" << name << ")\nError: can not find the file fiducial_cuts with name [" << fiducialCutsPath << "]. " 
	    << "Please specify the correct path in your config file:\n   "
	    << name << ".FiducialCutsFile: path/fiducial_cuts.txt\n";
      throw runtime_error("Can not find fiducial_cuts file");
      }
      if( defsfile == "") {
	err() << "ElectronSelector(" << name << ")\nError: can not find the file fiducial_defs with name [" << fiducialDefsPath << "]. " 
	    << "Please specify the correct path in your config file:\n   "
	    << name << ".FiducialCutsFile: path/fiducial_defs.txt\n";
	throw runtime_error("Can not find fiducial_defs file");
      }
      _fid = new emIDFiducial(cutsfile,defsfile,_fiducial,_fiducialVersion);
      cout << "Fiducial Cuts *enabled*" << endl;
      cout << "Using cut set: " << _fiducial << " Version " << _fiducialVersion << endl;
      cout << "from cut file [" << cutsfile << "] " << endl 
	   << "from def file [" << defsfile << "]" << endl ;
      cout << "Applying cuts: ";

      for( int i =0,n=_fidCutsToApply.size();i<n;++i)
	cout << _fidCutsToApply[i] << " ";
      cout << endl;
      

      _fid->PtCut( _pt );
      if (_pt>0) cout << "   pT cut = " << _pt << endl;
      else cout << "  NO pT cut" << endl ;
    }else {
      _fid = 0;
      if ( _pt>0 || _ptmax || _etaCC>0 || (_etaECmin>0 && _etaECmax >0) ) {
	printf("\nAdditional object selection:\n");
	
	if ( _pt > 0 ) {
	  char sbuf[256];
	  sprintf(sbuf,"Electron pT > %4.1f GeV",_pt);
	  _ptCutString=sbuf;
	  out() <<"  "<<_ptCutString<<"\n";
	}
	
	if ( _ptmax > 0 ) {
	  char sbuf[256];
	  sprintf(sbuf,"Electron pT < %4.1f GeV",_ptmax);
	  _ptmaxCutString=sbuf;
	  out() <<"  "<<_ptmaxCutString<<"\n";
	}
	_etaCutString="|CalDetectorEta|";
	char sbuf[256];
	if ( _etaCC > 0 ) { sprintf(sbuf," < %3.1f",_etaCC); _etaCutString+=sbuf; }
	if ( _etaCC > 0 && _etaECmin>0 && _etaECmax >0 ) _etaCutString+=" OR";
	if ( _etaECmin>0 && _etaECmax >0 ) { 
	  sprintf(sbuf," in (%3.1f, %3.1f)",_etaECmin,_etaECmax);
	  _etaCutString+=sbuf;
	}
	
	if ( _etaCC>0 || (_etaECmin>0 && _etaECmax >0) ) out() <<"  "<<_etaCutString<<"\n";
	else out() <<"  No cuts on CalDetectorEta\n";
      }
    else out() <<"\nNo additional pT or CalDetectorEta selection\n";
    }

    if ( _createLeadingBranch )
	printf("\n  The branch %s will be created (copy of the leading object)\n",
	       _leadingBranchName.c_str());

    if ( _nelectrons > 0 || _nelectronsmax > 0 ) {
      printf("\nEvent selection:\n");
      if ( _nelectrons == _nelectronsmax ) printf("  N selected electrons == %d\n",_nelectrons);
      else {
	if ( _nelectrons    > 0 ) printf("  N selected electrons >= %d\n",_nelectrons);
	if ( _nelectronsmax > 0 ) printf("  N selected electrons <= %d\n",_nelectronsmax);
      }
    }

    out() << "\n===================================================\n\n";
    
  addVariable("_id") ;
  if (_id->cutOnIsolation()) addVariable("_iso") ;
  if (_id->cutOnEMFraction()) addVariable("_EMfrac") ;
  if (_id->cutOnpT() || _pt>0 || _ptmax>0 ) {
    addVariable("fX") ;
    addVariable("fY") ;
  }
  bool etavar = false ;
  if ( _etaCC>0 || _etaECmin>0 || _etaECmax>0 ||
       _id->cutOnEtaCC() || _id->cutOnEtaECmin() || _id->cutOnEtaECmax() ) {
    addVariable("_CalDetectorEta") ;
    etavar = true ;
  }

  bool hmx7var = false ;
  if (_id->cutOnHMx7CC()) {
    if (!etavar) {
      addVariable("_CalDetectorEta") ;
      etavar = true ;
    }
    addVariable("_HMx7") ;
    hmx7var = true ;
  }

  if (_id->cutOnHMx8EC()) {
    if (!etavar) {
      addVariable("_CalDetectorEta") ;
      etavar = true ;
    }
    addVariable("_HMx8") ;
  }

  if (_id->cutOnHMx7()) {
    if (!hmx7var) {
      addVariable("_HMx7") ;
      hmx7var = true ;
    }
  }

  if (_id->cutOnTrkMatchChi2EOP())  {
    addVariable("_TrMatchChi2ProbBest") ;
    addVariable("_SpatialTrMatchChi2ProbBest") ;
  } if (_id->cutOnTrkMatchChi2()) addVariable("_SpatialTrMatchChi2ProbBest") ;


  if (_id->cutOnLHood()) {
    addVariable("_Lhood8") ;
    addVariable("_chptr") ;
  }
  if (_id->cutOnTrackPt()) {
    addVariable("_chptr") ;
    if (_id->cutOnTrkMatchChi2EOP() ) addVariable("_chptrBest") ;
    _finder = new TRefFinder() ;
    _finder->SetVariablesTrack(Variables("fUniqueID", "fX", "fY")) ;
  }
  _id->saveCuts() ;

  if(!readCELLs())
        printf("\n ### NO TheCellLocTable.dat file! ###\n");
  else  printf("\n ### File TheCellLocTable.dat is loaded. ###\n");

  lHood = new EMLikelihood("./caf_util/configs/templates_ZeeDatEMincDat.root");

  }


bool ElectronSelector::selectObject(const TMBEMCluster &electron)
{

  Collection<TMBPrimaryVertex> vertices = _event->getPrimaryVertices();
  if(vertices.size() > 0 )
    zvtx = (vertices.begin())->vz();
  else
    zvtx = -1000;

  float trkpt = 1000.0 ;
  if (_finder) {
    const TMBTrack* trk = 0 ;
    //    cout << (electron.GetChargedTrackRef().GetUniqueID() & 0xffffff )
    //   << " " << (electron.getPtrChpSpatialRef().GetUniqueID() & 0xffffff )
    //	 << " " << (electron.getPtrChpRef().GetUniqueID() & 0xffffff)  << endl ;
    if (_id->cutOnTrkMatchChi2EOP() && electron.track_match_chi2prob()>0) {
      trk = _finder->FindTrack(*_event, electron.getPtrChpRef()) ;
    } else if (electron.track_match_spatialchi2prob()>0) {
      trk = _finder->FindTrack(*_event, electron.getPtrChpSpatialRef()) ;
      //      trk = _finder->FindTrack(*_event, electron.GetChargedTrackRef()) ;
    }

    if (!trk) 
      trkpt = -1.0 ;
    else 
      trkpt =trk->Pt() ;
  }

  double IsoHC4 = trk_iso_ann(electron, zvtx, 0.05, 0.4);
  double Isolation_n = -1000.;
  if (_ver >= 10) {
    Isolation_n = Isolation(electron);
    Isolation_n = Isolation_n / electron.CalE();
  }
  double emhits_e_f = electron.emhits_e_f_discriminant();
  double Sigphi = electron.flrS1(3);

 // ANN 
  bool inCC = fabs(electron.CalDetectorEta()) < 1.1;
  bool inEC = fabs(electron.CalDetectorEta()) > 1.5 && fabs(electron.CalDetectorEta()) < 3.3;  

  double NNout3_cc = -1;
  double NNout4_cc = -1;
  double NNout5_cc = -1;
  double NNout7_cc = -1;
  double NNout3_ec = -1;

  double hmx8 = electron.HMx8();

  if(inCC) { 
 
     double _Eratio_ = 0.004;
     double _Ethreshold_ = 0.25;

     double threshold = _Eratio_*( electron.Pt()) + _Ethreshold_;
     if( threshold < 0.45 ) threshold = 0.45;
     double EM1fr = 0.;
     int EM1ncells = 0;   
     EMClusterNcell( electron, threshold, EM1ncells, EM1fr);  
     double EM1ncells_02 = cone_cell(electron, threshold, 0.2);
     double EM1ncells_04 = cone_cell(electron, threshold, 0.4);
     double EM1ncells_0204 = EM1ncells_04 - EM1ncells_02;
     int ntrks = cone5_ntrks(electron, zvtx, 0.05);  

      double e_RMS_CPS  = electron.E_RMS_CPS();
      double e2_RMS_CPS = electron.E2_RMS_CPS();
      int cps_total = 0;
      int cps_tight = 0;
      CPS_cluster(electron, cps_total, cps_tight);


//     double VARIABLE[11] ={EM1fr,EM1ncells,EM1ncells_0204,0.0,
//                            0.0,0.0,0.0,IsoHC4,
//                            ntrks, 0.0, 0.0};
//     NNout3_cc = eCCANN_3VAL(VARIABLE);
//     NNout4_cc = eCCANN_4VAL(VARIABLE);
//     NNout5_cc = eCCANN_5VAL(VARIABLE);

      double VARIABLE_ann7[11] = {EM1fr,EM1ncells,EM1ncells_0204,0.0,cps_total,
                                  e_RMS_CPS,e2_RMS_CPS,IsoHC4,ntrks, 0.0, 0.0};

      NNout7_cc = eCCANN_7VAL(VARIABLE_ann7);   

  }

  else if(inEC) {  

     double EM1fr_ec = 0.;
     int EM1ncells_ec = 0;
     double threshold_ec = 0.200;   
     EMClusterNcell_EC(electron, threshold_ec, EM1ncells_ec, EM1fr_ec);    
     int EM1ncells_02ec = cone_cell_EC(electron, threshold_ec, 0.2);    
     int EM1ncells_04ec = cone_cell_EC(electron, threshold_ec, 0.4);
     int EM1ncells_0204ec = EM1ncells_04ec - EM1ncells_02ec;

     double VARIABLE_EC[11] = {0.0,EM1ncells_ec,0.0,0.0,
                               0.0,hmx8,0.0,IsoHC4,
                               0.0,0.0,0.0};

     NNout3_ec = eECANN_3VAL(VARIABLE_EC);    

  }
     
  //LHood 
  double Lhood8_p20 = -1; 
  Bool_t cL = lHood->calcLhood( &electron, _event); 
  if (cL) Lhood8_p20 = lHood->getLhood8();  


  bool pass = _id->select(1, 
			  electron.id(), 
			  electron.iso(), 
			  electron.emfrac(), 
			  electron.Pt(), 
			  electron.CalDetectorEta(), 
			  electron.HMx7(), 
			  electron.HMx8(), 
			  electron.track_match_spatialchi2prob(), 
			  electron.Lhood8(),
			  trkpt,
			  electron.track_match_chi2prob(),
			  // new vars
                          Isolation_n, IsoHC4, NNout7_cc, NNout3_ec, Lhood8_p20, 
                          electron.emhits_e_f_discriminant(), Sigphi
			  ) ;

  const vector<string>* cuts = _id->passCuts() ;
  if (cuts && _nelectrons > 0 ) 
    for (vector<string>::const_iterator it = cuts->begin();
       it != cuts->end(); it++) _stat.EventSelected(*it) ;

  if (!pass) return false; 

  if( _doNewFiducial ) {
    const TMBVertex * vtx = electron.GetVertex();
    double vertexz = 0.;
    if( vtx )
      vertexz = vtx->z();
 
    emIDFiducial::FidCutResults fidresults = _fid->Select( electron.CalDetectorEta(), electron.Phi(),vertexz, electron.Pt(),"");
    if( _fidCutsToApply.size() == 0)
      pass = pass && fidresults;
    else
      pass = pass && fidresults.accept(_fidCutsToApply);

    if (!pass) return false;
    if ( _nelectrons > 0 ) _stat.EventSelected("Fiducial cuts [" + _fiducial + "]");
  }
  else { 
    /// Additional Pt cut
    if ( _pt > 0 && electron.Pt() < _pt ) return false;
    if ( _pt > 0 && _nelectrons > 0 ) _stat.EventSelected(_ptCutString);
    
    if ( _ptmax > 0 && electron.Pt() > _ptmax ) return false;
    if ( _ptmax > 0 && _nelectrons > 0 ) _stat.EventSelected(_ptmaxCutString);
    
    /// Additional Eta cut
    if ( _etaCC > 0 || ( _etaECmin>0 && _etaECmax>0 ) ) {
      bool inCC_ = _etaCC > 0 ? fabs(electron.CalDetectorEta()) < _etaCC : false;
      bool inEC_ = ( _etaECmin>0 && _etaECmax>0 ) ? 
	fabs(electron.CalDetectorEta()) > _etaECmin && fabs(electron.CalDetectorEta()) < _etaECmax : false;
      if ( !inCC_ && !inEC_  ) return false;
      if ( _nelectrons > 0 ) _stat.EventSelected(_etaCutString);
    }
  }
  _nselected++ ;

  if (debug()>1) 
    cout << "PROCESSOR \"" << name() 
	 << "\" SELECTED OBJECT: pT = " << electron.Pt() 
	 << " eta = " << electron.Eta() 
	 << " phi = " << electron.Phi() 
	 << " charge = " << electron.charge() 
	 << endl ;
  return true ;
}

void ElectronSelector::finish() {
  if (debug() >=1 ) _id->Stat() ;
  _id->Stat();

}

  bool ElectronSelector::processEvent(cafe::Event &event)
  {

    //  cout << " New event ! " << endl;

  Collection<TMBEMCluster> allEMs = event.getEMscone();
  Collection<TMBEMCluster>::iterator emcl;

  const TMBCellContainer*  clist  = event.getCaloCells();

//   for ( emcl=allEMs.begin(); emcl!=allEMs.end(); ++emcl) {

//      bool probe_quality =
//      fabs(emcl->CalDetectorEta())<2.5                      &&
//      emcl->iso()<0.2                                       &&
//      emcl->emfrac()>0.80                                   &&
//      abs(emcl->id()) < 12                                  &&
//      emcl->CalE()*sin(kinem::theta(emcl->CalEta())) > 15.0 &&
//      emcl->is_in_eta_fiducial();
//      if(!probe_quality) continue;
//      //    cout << " processEvent:emcl.ncells() = " << emcl->ncells() << endl;
//      for(UInt_t c = 0; c < emcl->ncells(); c++) {
//        //          cout << "processEvent:cell 1: " << c << endl;
//           const TMBCaloCell *cell = emcl->GetCaloCell(c);
// 	  //          cout << "processEvent:cell 2: " << c << endl;
//  //          if(cell==0) cout << "cell is not found!!!" << endl;
// //           if(cell==0) continue;
//        //    bool normal = cell->isNormalCell();
//           if(!cell->isNormalCell()) continue;
//           //cout << " normal = "<< normal << endl;
//           //if( !normal )       continue;
//      } 
//   }


    //get pointer to statistics collector
    event.get("StatPointer", _stat)  ;
    _event = &event ;
    _nselected = 0 ;
    SelectUserObjects<TMBEMCluster>::processEvent(event);

    if (debug() >=2 ) 
      out() << " ELECTRON SELECTED: " << _nselected 
	    << endl ;
    
    for (int i = 1; i <= _nelectrons; i++) {
      if (_nselected < i) return false ; 
      ostringstream st ;
      st << " N electrons >= "  << i;
      _stat.EventSelected(st.str()) ;
    }

    if (_nelectronsmax >= 0) {
      if (_nselected > _nelectronsmax) return false ; 
      ostringstream st ;
      st << "Number of electrons <= " <<  _nelectronsmax ;
      _stat.EventSelected(st.str()) ;      
    }   

    return true  ;	
  }

  void ElectronSelector::before(cafe::Collection<TMBEMCluster>& from) {
      if(_sort) {
	  from.sort(from.begin(),from.end(),::moreThan) ;
      } 
     
  }

  void ElectronSelector::after(cafe::Collection<TMBEMCluster>& accepted,
			       cafe::Collection<TMBEMCluster>& rejected) {
      if ( _treehandler && _createLeadingBranch ) {
	  Collection<TMBEMCluster> leadingEM;
	  if (accepted.size()>0) leadingEM.push_back(accepted[0]);
	  _treehandler->FillBranch(_leadingBranchName,leadingEM);
      }
  }

  void ElectronSelector::inputFileOpened(TFile *file) {
      SelectUserObjects<TMBEMCluster>::inputFileOpened(file);
      if (_treehandler) {
	  _treehandler->InitTree(_toTreeName,file);
	  if (_createLeadingBranch) _treehandler->AddBranch(_leadingBranchName,"TMBEMCluster");
      }
  }

  void ElectronSelector::inputFileClosing(TFile *file) {
      if (_treehandler) _treehandler->ClearIt();
      SelectUserObjects<TMBEMCluster>::inputFileClosing(file);
  }


 double ElectronSelector::Isolation(const TMBEMCluster &ele)
 {
    double lumi_factor_CC = 0.00685;
    double lumi_factor_EC = 0.02250;
    double IsoE   = ele.EisoCore()*ele.iso();
    const TMBGlobal* glob = _event->getGlobal();
    double instLum = glob->instlum()*36.;
    if(fabs(ele.CalDetectorEta()) < 1.1)
      IsoE   = IsoE   - instLum * lumi_factor_CC;
    else
      IsoE   = IsoE   - instLum * lumi_factor_EC;

    if(IsoE < 0.) IsoE = 0.;

    return IsoE;

 }

// CALO cells:
//**********************************************************
// getCELLcoord()
//**********************************************************
      
 bool ElectronSelector::getCELLcoord(int iEta, int iPhi, int iLayer, double& x, double& y, double& z)
 {
  if (iEta < minEta     || iEta > maxEta ||
      iLayer < minLayer || iLayer > maxLayer ||
      iPhi < minPhi     || iPhi > maxPhi)
    return false;

  long mangled = iLayer*10000+iPhi*100+iEta;
  std::map<long, xyz>::const_iterator it = CELLtable.find(mangled);
  if (it == CELLtable.end()) return false; // can't find it.
  const xyz& ref = it->second;
  x = ref.x; y=ref.y; z=ref.z;
  return true;
}

//**********************************************************
// readCELLs()
//**********************************************************
bool ElectronSelector::readCELLs()
{
  ifstream infile("./caf_util/configs/TheCellLocTable.dat");
  if ( !infile) return false;

  while (infile) {
    long key;
    double x,y,z;
    infile >> key >> x >> y >> z;
    xyz e = {x, y, z};
    CELLtable [key] = e;
  }

  infile.close();
  return true;
}

//=======================
// ANN methods
//=======================
//
//-----------------------------
// Parametrized ANN functions:
//----------------------------

//namespace caf_util {

// Neural Network output variable
// use EM1frac+NCells+TrkIso+ntrks+conecells+nCPS+cps_rms2
double ElectronSelector::eCCANN_7VAL(double pattern[11])
{

  // NN node numbers: Input, hidden, and output
  int MaxIndim = 11;
  int Indim = 7;
  int Hidden = 4;

  int ivalid = 0;
  double I_SUM = 0.0;
  double H_O_SUM=0.0;

   // Input scales
 float scale[] = {1.000000,100.000000,10.000000,0.000000,100.000000
                        ,0.100000,0.100000,1000.000000,10.000000,0.000000
                        ,0.000000};

// veto input patterns (0 == DO NOT USE)
 int veto[] = {1,1,1,0,1,0,1,1,1,0,0};

// NN node weights and threshholds
 float I_H_Weight[] = {3.997103,4.439505,-2.596792,1.793211,-5.961729
                ,16.274263,16.936739,-15.617092,2.901612,-0.062824
                ,2.832576,-57.351852,4.814333,4.869636,3.653307
                ,1.672297,-12.886433,8.774302,8.917918,8.684820
                ,-114.349709,117.369705,124.744110,-165.534042,8.223936
                ,5.301783,4.564135,141.896973};
 float H_Thresh[] = {-3.241092,-3.305562,-1.086682,-1.311124};
 float H_O_Weight[] = {3.222120,-2.726241,-5.399126,3.660837};
 float O_Thresh = -0.347845;
 //============================================

 for (int h = 0; h<Hidden; h++)
   {
     I_SUM = 0.;
     ivalid = 0;
     for (int i = 0; i<MaxIndim; i++)
       {
         if(veto[i])
           {
             I_SUM += I_H_Weight[Hidden*ivalid+h] * pattern[i]/scale[i];
             ivalid++;
           }
       }
     I_SUM += H_Thresh[h];
     H_O_SUM += H_O_Weight[h]*( 1.0 / ( 1.0 + exp(-2.0*I_SUM)));
   }
 H_O_SUM += O_Thresh;
 return ( 1.0 / ( 1.0 + exp(-2.0*H_O_SUM) ) );

};


// Neural Network output variable
// use EM1frac+NCells+TrkIso
double ElectronSelector::eCCANN_3VAL(double pattern[11])
{
  // NN node numbers: Input, hidden, and output
  int MaxIndim = 11;
  int Indim = 3;
  int Hidden = 4;
  
  int ivalid = 0;
  double I_SUM = 0.0;
  double H_O_SUM = 0.0;

  // Input scales
 double scale[] = {1.000000,100.000000,100.000000,1000.000000,1000.000000
                        ,1000.000000,10.000000,10000.000000,10.000000,100.000000
                        ,0.100000};

// veto input patterns (0 == DO NOT USE)
 int veto[] = {1,1,0,0,0,0,0,1,0,0
                        ,0};

// NN node weights and threshholds
 double I_H_Weight[] = {-0.245486,0.314878,-0.099278,-0.209881,3.577051
                ,-3.727264,4.423425,4.074361,284.924347,-284.052216
                ,284.775299,284.440765};
 double H_Thresh[] = {-0.297218,0.178977,-0.451276,-0.406109};
 double H_O_Weight[] = {-7.482326,9.414640,-8.031065,-7.231708};
 double O_Thresh = 3.480176;


// Loop over output, hidden, and input nodes
 // to produce the neural net output 
 //============================================
 
 for (int h = 0; h<Hidden; h++)
   {
     I_SUM = 0.;
     ivalid = 0;
     for (int i = 0; i<MaxIndim; i++)
       {
	 if(veto[i])
	   {
	     I_SUM += I_H_Weight[Hidden*ivalid+h] * pattern[i]/scale[i];
	     ivalid++;
	   }
       }
     I_SUM += H_Thresh[h];
     H_O_SUM += H_O_Weight[h]*( 1.0 / ( 1.0 + exp(-2.0*I_SUM)));
   }
 H_O_SUM += O_Thresh;
 return ( 1.0 / ( 1.0 + exp(-2.0*H_O_SUM) ) );
};


// Neural Network output variable
// use EM1frac+NCells+TrkIso+ntrks
double ElectronSelector::eCCANN_4VAL(double pattern[11])
{
  // NN node numbers: Input, hidden, and output
  int MaxIndim = 11;
  int Indim = 4;
  int Hidden = 4;

  int ivalid = 0;
  double I_SUM = 0.0;
  double H_O_SUM = 0.0;

  // Input scales
 double scale[] = {1.000000,100.000000,100.000000,1000.000000,1000.000000
                        ,1000.000000,10.000000,10000.000000,10.000000,100.000000
                        ,0.100000};

// veto input patterns (0 == DO NOT USE)
 int veto[] = {1,1,0,0,0,0,0,1,1,0
                        ,0};

// NN node weights and threshholds
 double I_H_Weight[] = {1.307330,0.576619,-0.193091,0.267955,-6.417233
                ,-5.363276,6.153428,-5.523644,-391.844940,-388.295166
                ,383.850983,-363.240845,2.313047,1.866142,1.464222
                ,41.607906};
 double H_Thresh[] = {-1.986848,-0.602294,-0.547508,0.503302};
 double H_O_Weight[] = {4.850314,4.805796,-12.817219,6.291244};
 double O_Thresh = -1.915041;


 //============================================

 for (int h = 0; h<Hidden; h++)
   {
     I_SUM = 0.;
     ivalid = 0;
     for (int i = 0; i<MaxIndim; i++)
       {
         if(veto[i])
           {
             I_SUM += I_H_Weight[Hidden*ivalid+h] * pattern[i]/scale[i];
             ivalid++;
           }
       }
     I_SUM += H_Thresh[h];
     H_O_SUM += H_O_Weight[h]*( 1.0 / ( 1.0 + exp(-2.0*I_SUM)));
   }
 H_O_SUM += O_Thresh;
 return ( 1.0 / ( 1.0 + exp(-2.0*H_O_SUM) ) );
};

// Neural Network output variable
// use EM1frac+NCells+TrkIso+ntrks+conecells
double ElectronSelector::eCCANN_5VAL(double pattern[11])
{
  // NN node numbers: Input, hidden, and output
  int MaxIndim = 11;
  int Indim = 5;
  int Hidden = 4;

  int ivalid = 0;
  double I_SUM = 0.0;
  double H_O_SUM = 0.0;

  // Input scales
 double scale[] = {1.000000,100.000000,100.000000,1000.000000,1000.000000
                        ,1000.000000,10.000000,10000.000000,10.000000,100.000000
                        ,0.100000};

// veto input patterns (0 == DO NOT USE)
 int veto[] = {1,1,1,0,0,0,0,1,1,0
                        ,0};

// NN node weights and threshholds
 double I_H_Weight[] = {-0.579286,0.044048,0.598504,-0.001557,7.065582
                ,-5.502957,-5.922191,-5.555092,9.488820,-7.120522
                ,-36.014442,-6.800014,423.613800,-425.063293,-397.159332
                ,-426.368958,2.037131,2.429470,38.408684,2.513462};
 double H_Thresh[] = {-0.625721,-0.886000,0.328361,-1.086039};
 double H_O_Weight[] = {-11.077871,4.973232,5.643337,4.736973};
 double O_Thresh = -2.018292;

 //============================================

 for (int h = 0; h<Hidden; h++)
   {
     I_SUM = 0.;
     ivalid = 0;
     for (int i = 0; i<MaxIndim; i++)
       {
         if(veto[i])
           {
             I_SUM += I_H_Weight[Hidden*ivalid+h] * pattern[i]/scale[i];
             ivalid++;
           }
       }
     I_SUM += H_Thresh[h];
     H_O_SUM += H_O_Weight[h]*( 1.0 / ( 1.0 + exp(-2.0*I_SUM)));
   }
 H_O_SUM += O_Thresh;
 return ( 1.0 / ( 1.0 + exp(-2.0*H_O_SUM) ) );
};

// ncells02 + trkIso005_04 + HMx8
double ElectronSelector::eECANN_3VAL(double pattern[11])
{
  // NN node numbers: Input, hidden, and output 
  int MaxIndim = 11;
  int Indim = 3;
  int Hidden = 4;

  int ivalid = 0;
  double I_SUM = 0.0;
  double H_O_SUM = 0.0;

  // Input scales
 double scale[] = {1.000000,100.000000,100.000000,100.000000,10000.000000
                        ,10000.000000,10.000000,1000.000000,100.000000,0.000000
                        ,0.000000};

// veto input patterns (0 == DO NOT USE)
 int veto[] = {0,1,0,0,0,1,0,1,0,0
                        ,0};

// NN node weights and threshholds
 double I_H_Weight[] = {5.382861,-3.808220,3.319416,3.094138,75.685478
                ,-77.542099,77.989510,77.157494,57.442909,-54.978638
                ,55.918785,55.673016};
 double H_Thresh[] = {-0.953798,0.306949,-0.277013,-0.261774};
 double H_O_Weight[] = {-4.459886,6.262339,-4.066019,-3.770394};
 double O_Thresh = 2.158080;

  // Loop over output, hidden, and input nodes 
  // to produce the neural net output  
  //============================================ 

  for (int h = 0; h<Hidden; h++)
    {
      I_SUM = 0.;
      ivalid = 0;
      for (int i = 0; i<MaxIndim; i++)
        {
          if(veto[i])
            {
              I_SUM += I_H_Weight[Hidden*ivalid+h] * pattern[i]/scale[i];
              ivalid++;
            }
        }
      I_SUM += H_Thresh[h];
      H_O_SUM += H_O_Weight[h]*( 1.0 / ( 1.0 + exp(-2.0*I_SUM)));
    }
  H_O_SUM += O_Thresh;
  return ( 1.0 / ( 1.0 + exp(-2.0*H_O_SUM) ) );
};

//------------------------
// Auxiliary ANN methods:
//------------------------

//calculate the number of CPS 3D clusters in dR=0.1 around an emcluster
int ElectronSelector::CPS_cluster(const  TMBEMCluster &gam, int &N1, int &N2)
{
 // CPS stuff
  Collection<TMBCps>  cps_3d   = _event->getCollection<TMBCps>("Cps");
  Collection<TMBCps>::iterator cps_iter;
  Collection<TMBCps> cps_my(0);
  for( cps_iter=cps_3d.begin(); cps_iter != cps_3d.end(); ++cps_iter ) {

      double TempR = (*cps_iter).r();
      double TempPhi = (*cps_iter).phi();
      double cps_xyz[3];
      cps_xyz[0]= TempR * cos(TempPhi);
      cps_xyz[1]= TempR * sin(TempPhi);
      cps_xyz[2]= (*cps_iter).z();
      double TempEta =kinem::eta(cps_xyz[0],cps_xyz[1],cps_xyz[2]); // 
      //double TempPhi = kinem::phi(cps_xyz[0],cps_xyz[1]);
      double dPhi = kinem::delta_phi(TempPhi,gam.CalDetectorPhi());
      double dEta = TempEta - gam.CalDetectorEta();
      double Dist = sqrt(dPhi*dPhi + dEta*dEta);
      // reject clusters with R > 0.1
      if ( Dist > 0.10 ) continue;
     
      cps_my.push_back(*cps_iter);
  } // loop over CPS clusters
  //cout<<"CPS size: "<<cps_my.size()<<endl;
 for ( int ie=0; ie != cps_my.size(); ++ie ) 
    {
     //all 3D clusters

     ++N1;
      if ( cps_my[ie].isTight() ) 
      {
        ++N2;
      }
  }

return 0;
};


//calculate the number of tracks in cone0.05 around an emcluster
int ElectronSelector::cone5_ntrks(const TMBEMCluster &em, double zvtx, double conesize) 
{ 
  int n_trks = 0;
  double dz_cut = 1.0;//cm 
  double emEta = em.CalEta();
  double emPhi = em.CalPhi(); 
 
  Collection<TMBTrack> tracks = _event->getTracks(); 
  Collection<TMBTrack>::iterator trk; 
 
  for( trk=tracks.begin(); trk != tracks.end(); ++trk )
    { 
      if ( fabs(trk->z() - zvtx) > dz_cut )             continue; 
      double dr = kinem::delta_R(trk->Eta(), trk->Phi(), emEta, emPhi);
      if( dr > conesize )        continue;
      if( trk->Pt() < 0.5 ) continue;
      n_trks++; 
    }
  return n_trks;
}

 // calculate the EM1 energy fraction and No. of EM1 cells in CC region
void ElectronSelector::EMClusterNcell(const TMBEMCluster &gam, double Thres, int& GamNcell,double& GamEM1fr)
{
   _event->getCaloCells();
  // Get Cells for EM cluster
  int ncell = 0;
  double Etot = 0;
  double E_EM1 = 0;
  for(UInt_t c = 0; c < gam.ncells(); c++) {
    const TMBCaloCell *cell = gam.GetCaloCell(c);
    int ilay = cell->layer();
    double Ecell = cell->E();
    Etot += Ecell;
    if ( ilay != 1 )    continue;
    bool normal = cell->isNormalCell();
    if( !normal )       continue;
    if (Ecell < Thres ) continue;
    ncell++;
    E_EM1 += Ecell;
  }
  GamNcell = ncell;
  //GamEM1fr = 0;
  if ( Etot>0 ) GamEM1fr = E_EM1 / Etot;
  //printf("\n gam.ncells() =%d",gam.ncells());
  //printf("\n EMClusterNcell: %d %f %f %f",GamNcell,GamEM1fr,E_EM1,Etot);
  // if (GamNcell < 1 && GamEM1fr < 0.00001 )  
//     printf("\n EMClusterNcell: %d %f %f %f   %d  %f %f\n",GamNcell,GamEM1fr,E_EM1,Etot,gam.ncells(),  
// 	   gam.CalE(), gam.CalDetectorEta());

  return;

}


// calculate the EM1 energy fraction and No. of EM1 cells in EC region
void ElectronSelector::EMClusterNcell_EC(const TMBEMCluster &gam,double Thres,int& GamNcell,double& GamEM1fr)
{
  _event->getCaloCells();
  // Get Cells for EM cluster
  int ncell = 0;
  double Etot = 0.;
  double E_EM1 = 0.;
  double detEta = gam.CalDetectorEta();
  for(UInt_t c = 0; c < gam.ncells(); c++) 
    {
      const TMBCaloCell *cell = gam.GetCaloCell(c);
      int ilay = cell->layer();
      double Ecell = cell->E();
      Etot += Ecell;
      bool normal = cell->isNormalCell();
      if( !normal )       continue;
      double cellET = Ecell/cosh(detEta);
      if (cellET < Thres ) continue;
      if( ilay == 1 )
      { 
        ncell++;
        E_EM1 += Ecell;
      }
    }
  GamNcell = ncell;
  GamEM1fr = 0.;
  if(Etot>0.) GamEM1fr = E_EM1/Etot;  
  
  return;
}

// calculate # of EM1 cells in cone 0.2 or 0.4 in CC region
int ElectronSelector::cone_cell(const TMBEMCluster &em, double Thres, double conesize)
{
  double _PI_ = TMath::Pi();
  double _TWOPI_ = 2.*TMath::Pi();
  double eEta = em.CalDetectorEta();
  double ePhi = em.CalDetectorPhi();
  const TMBCellContainer *cells = _event->get<TMBCellContainer>("CaloCells");
  
  int n_cell = 0;
  if(CaloGeometryManager::get_instance())
    {
      const map<int, CELLXYZ_STRUCT>& calogeometry_map = 
      CaloGeometryManager::get_instance()->getCaloGeometry2005Map();
      for(UInt_t c = 0; c < cells->NumCells(); c++)
	{
          const TMBCaloCell *cell = cells->GetCell(c);
	  //const TMBCaloCell *cell = em.GetCaloCell(c);
	  double Ecell = cell->E();
	  if( Ecell < Thres ) continue;
	  int ilay = cell->layer();
	  if( ilay != 1 )     continue;
	  bool normal = cell->isNormalCell();
	  if( !normal )       continue;
	  int iEta = cell->ieta();
          int iPhi = cell->iphi();
          int cell_index = wz_utils::EncodeCellIndex(iEta, iPhi, ilay);
          // calorimeter cell position
          double cell_pos[3] = {0., 0., 0.};
          map<int, CELLXYZ_STRUCT>::const_iterator iter = calogeometry_map.find(cell_index);
          if( iter != calogeometry_map.end() )
            {
              cell_pos[0] = (*iter).second.x;
              cell_pos[1] = (*iter).second.y;
              cell_pos[2] = (*iter).second.z;
	      double eta_cell = kinem::eta(cell_pos[0],cell_pos[1],cell_pos[2]);
              double phi_cell = kinem::phi(cell_pos[0], cell_pos[1]);
              double deta = fabs(eta_cell - eEta);
              double dphi = fabs(phi_cell - ePhi);
              if( dphi > _TWOPI_ ) dphi -= _TWOPI_;
              if( dphi > _PI_ ) dphi = _TWOPI_ - dphi;
              double dR = sqrt(deta*deta + dphi*dphi);
	      if( dR < conesize )
		++n_cell;
	    }

	}//cell loop
    }//get_instance
  return n_cell;
};

// calculate # of EM1 cells in cone 0.2 or 0.4 in EC region
int ElectronSelector::cone_cell_EC(const TMBEMCluster &em, double Thres, double conesize)
{

   double _PI_ = TMath::Pi();
  double _TWOPI_ = 2.*TMath::Pi();
  double eEta = em.CalDetectorEta();
  double ePhi = em.CalDetectorPhi();  
  const TMBCellContainer *cells = _event->get<TMBCellContainer>("CaloCells");  
  int n_cell = 0;
  //double detEta = em->CalDetectorEta();
  if(CaloGeometryManager::get_instance())
    {      
      const map<int, CELLXYZ_STRUCT>& calogeometry_map =
       CaloGeometryManager::get_instance()->getCaloGeometry2005Map();
      for(UInt_t c = 0; c < cells->NumCells(); c++)
        {
          const TMBCaloCell *cell = cells->GetCell(c);
          //const TMBCaloCell *cell = em.GetCaloCell(c);
          double Ecell = cell->E();
	  double cellET = Ecell/cosh(eEta);
          if( cellET < Thres ) continue;
          int ilay = cell->layer();
          if( ilay != 1 )     continue;
          bool normal = cell->isNormalCell();
          if( !normal )       continue;
          int iEta = cell->ieta();
          int iPhi = cell->iphi();
	  int cell_index = wz_utils::EncodeCellIndex(iEta, iPhi, ilay);
          // calorimeter cell position
          double cell_pos[3] = {0., 0., 0.};
          map<int, CELLXYZ_STRUCT>::const_iterator iter = calogeometry_map.find(cell_index);
          if( iter != calogeometry_map.end() )
            {
              cell_pos[0] = (*iter).second.x;
              cell_pos[1] = (*iter).second.y;
              cell_pos[2] = (*iter).second.z;
              double eta_cell = kinem::eta(cell_pos[0],cell_pos[1],cell_pos[2]);
              double phi_cell = kinem::phi(cell_pos[0], cell_pos[1]);
              double deta = fabs(eta_cell - eEta);
              double dphi = fabs(phi_cell - ePhi);
              if( dphi > _TWOPI_ ) dphi -= _TWOPI_;
              if( dphi > _PI_ ) dphi = _TWOPI_ - dphi;
              double dR = sqrt(deta*deta + dphi*dphi);
              if( dR < conesize )
                ++n_cell;
            }

        }//cell loop
    }//get_instance
  return n_cell;
};

                           
double ElectronSelector::trk_iso_ann(const TMBEMCluster &em, double zvtx, double rmin, double rmax)
{
  double eta = em.CalEta();
  double phi = em.CalPhi();
  double dz_cut = 1.0;
  double iso = 0.0;
  //cout<<"COMPARISON: "<<em.Pt();
  const TMBTrack *electron = em.getPtrChp();
  Collection<TMBTrack> tracks = _event->getTracks();
  Collection<TMBTrack>::iterator trk;
  vector<TMBTrack*> final_tracks;
  final_tracks.clear();
  for( trk=tracks.begin(); trk != tracks.end(); ++trk )
   {
    if(electron) if(size_t(&(*electron)) == size_t(&(*trk))) continue;
    if ( fabs(trk->z() - zvtx) > dz_cut ) continue;
    double dr = kinem::delta_R(trk->Eta(), trk->Phi(), eta, phi);
    double trkPT = trk->Pt();
    if( trkPT<0.5) continue;

    if((dr>rmin) && (dr<rmax))
       {
        iso += trkPT;
       //final_tracks.push_back((TMBTrack*)&(*trk));
       }//dr
   }//trk
 
   return iso;

}


 //} // namespace



} // namespace

ClassImp(caf_util::ElectronSelector) ;