/*
  Met Resolution Model 
*/
#include <iomanip>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <Stntuple/loop/TStnAna.hh>
#include <Stntuple/jet/TStnMetModel.hh>
#include <TRandom.h>

//------- need this for jet-EM matching
bool SortTowersByEnergy(TCalTower* c1, TCalTower* c2){
  if(c1->Energy() > c2->Energy()) return true;
  else return false;
}

//------- Function for MET resolution due to unclustered energy
double UnclMetResolution(double* x, double* par){
  double value=0.0;
  double mean=par[0];
  double sigma1=par[1];
  double sigma2=par[3]*sigma1;
  double arg=x[0];
  value=(TMath::Gaus(arg,mean,sigma1)+par[2]*TMath::Gaus(arg,mean,sigma2))/(1.0+par[2]);
  return value;
}

//-------------- raw MetSig fit function: p0*exp(-p1*x+p2*x*x)+p3*exp(-p4*x)
double rawMetSigFunc(double *x, double *par)
{
  double fitval=0.0;
  if(x[0]>0.0)
    {
      fitval=par[0]*exp(-1.0*par[1]*x[0]+par[2]*x[0]*x[0])+par[3]*exp(-1.0*par[4]*x[0]);
    }
  else fitval=-1.0E6;
  return fitval;
}

//-------- Jet Energy Resolution: centered at zero 
// //-------------- fit function: [Const*Gaus(-x/(1+x))+Landau(-x/(1+x))]/(1+Const)
double MyJER(double* x, double* par) {
  double value=0.0;
  double arg=-x[0]/(x[0]+1.0);
  double arg_L=-x[0]/(x[0]+1.0);
  double mean=par[0];
  double sigmaG=par[1];
  double mpv=par[2];
  double sigmaL=par[3];
  double normL=par[4];
  if(normL<0.0) normL=0.0;
  double f1=normL/(1.0+normL)*TMath::Gaus(arg,mean,sigmaG);
  double f2=TMath::Landau(arg_L,mpv,sigmaL)/(1.0+normL);
  value=f1+f2;
  if(value<0.0) return 0.0;
  return value;
}


ClassImp(TStnMetModel)
//_____________________________________________________________________________
TStnMetModel::TStnMetModel(){
  DebugMode=0; 
  //------ CEM Photon ID parameters
  fPhoMinEt=10.0;
  fCemPhoMaxHADEM[0]=0.125;        
  fCemPhoMaxHADEM[1]=0.055;        
  fCemPhoMaxHADEM[2]=0.00045;                                      
  fCemPhoMaxCesChi2=20.0;
  fCemPhoMaxN3D=1;  
  fCemPhoMaxPtmax[0]=1.0;     
  fCemPhoMaxPtmax[1]=0.005;     
  fCemPhoMaxCalIso4_corr[0]=0.1;    
  fCemPhoMaxCalIso4_corr[1]=2.0;    
  fCemPhoMaxCalIso4_corr[2]=0.02;    
  fCemPhoMaxTrkIso4[0]=2.0;
  fCemPhoMaxTrkIso4[1]=0.005;
  fCemPhoMax2ndCesEt[0]=0.14; 
  fCemPhoMax2ndCesEt[1]=2.4;      
  fCemPhoMax2ndCesEt[2]=0.01;         
  fCemPhoMinCesX=0.0;           
  fCemPhoMaxCesX=21.0;
  fCemPhoMinCesZ=9.0; 
  fCemPhoMaxCesZ=230.0;
  //------- CEM Electron ID parameters
  fEleMinEt=10.0;    
  fEleTrkMinPt=5.0; 
  fCemEleMaxHADEM[0]=0.055;                      
  fCemEleMaxHADEM[1]=0.00045;                      
  fCemEleMaxIsoFr=0.1;  
  fCemEleMaxTrkZ0=60.0;
  fCemEleFidCut=1; // cut on CEM fiduciality
  fCemEleMinCotAxSl=1; // cut on number of COT Axial SL with 5 or more hits 
  fCemEleMinCotStSl=1; // cut on number of COT Stereo SL with 5 or more hits 

  //-------- cuts for PEM electrons
  fPemEleMinEtaPes=1.2; 
  fPemEleMaxEtaPes=3.0;
  fPemEleMaxHADEM=0.05;       
  fPemEle3x3FitCut=0; // cut on PEM 3x3 tower fit status code
  fPemEleMax3x3Chi2=10.0;
  fPemEleMinPes5x9U=0.65; 
  fPemEleMinPes5x9V=0.65; 
  fPemEleMaxIsoFr=0.1;
  fPemEleMaxPesdR=3.0;

  //------- Jet energy corrections;
  fJTC_systcode=0; // default
  //------ Other jet stuff
  MinJetEta=0.0;
  MaxJetEta=3.0;
  //------ Metsig related stuff
  fSampleID=10; // by default, not metsig calibartion
  fUnclParamSwitch=1; // use Zee parameterization by default; 
  Npseudo=0; // no pseudo-experiments by default

  ClearExternalInput(); // have to call it here for initial cleanup
//   std::cout<<"Hi, entering TStnMetModel"<<std::endl;
}

//_____________________________________________________________________________
TStnMetModel::~TStnMetModel() {
}


//_____________________________________________________________________________
// Main Method to run MetModel
int TStnMetModel::RunMetModel(TStnEvent* event) {

  int code=UnpackBlocks(event);
  //--------- clearing input stuff
  ClearCommonStuff(allstuff); // clears common stuff
  ClearJetStuff(jet04stuff); // clears JetClu-0.4 stuff 
  ClearMetsigStuff(metsigstuff);
  //---- accessing vertex info
  Zvx_best=0.0;
  myNvx_class12=GetMyNvx12(Zvx_best);
  //_________________________ Initializing parameters
  DoCommonStuff(allstuff); // does rawSumEt, MET, EM objects, muons/tracks
  DoMyJet(fJetBlockClu04,allstuff,jet04stuff); // does my JetClu04 jets
  DoMyMet(allstuff,jet04stuff); // corrects Met & SumEt for jets, EM objects, muons/tracks
  DoMyMetSignificanceData(allstuff,jet04stuff,metsigstuff); // calculates MET-sig for data and MC
  DoMyMetSignificancePseudo(allstuff,jet04stuff,metsigstuff); // calculates MET-sig for data and MC

  return 0;
}


//---------------------------------------------------------------------------------------------------------
//
//   >>>>>>>>>>>>>>>>>>>>>>>>   ALL Additional Methods are described HERE <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
//
//_________________________________________________________________________________________________________

//_____________________________________________________________________________
int TStnMetModel::UnpackBlocks(TStnEvent* event) {

  fHeaderBlock=(TStnHeaderBlock*)event->UnpackDataBlock("HeaderBlock");
  if(!fHeaderBlock) return 1;

  fCalData= (TCalDataBlock*)event->UnpackDataBlock("CalDataBlock");
  if(!fCalData) return 2;

  fZVertexBlock=(TStnVertexBlock*)event->UnpackDataBlock("ZVertexBlock");
  if(!fZVertexBlock) return 3;

  fJetBlockClu04=(TStnJetBlock*)event->UnpackDataBlock("JetBlock");
  if(!fJetBlockClu04) return 4;

  fPhotonBlock=(TStnPhotonBlock*)event->UnpackDataBlock("PhotonBlock");
  if(!fPhotonBlock) return 5;

  fElectronBlock=(TStnElectronBlock*)event->UnpackDataBlock("ElectronBlock");
  if(!fElectronBlock) return 6;

  fMetBlock=(TStnMetBlock*)event->UnpackDataBlock("MetBlock");
  if(!fMetBlock) return 7;

  fTrackBlock=(TStnTrackBlock*)event->UnpackDataBlock("TrackBlock");
  if(!fTrackBlock) return 8;

  return 0;
}

//________________________ accessing vertex information
int TStnMetModel::GetMyNvx12(double &zvx)
{
  if(!fZVertexBlock) return 0;
  //---- accessing vertex info
  int myNvx_all=fZVertexBlock->NVertices();
  int myNvx12=0;
  double sumPt_max=0.0;
  for (int j=0; j<myNvx_all; j++) 
    {
      TStnVertex *vertex= fZVertexBlock->Vertex(j);
      if(vertex->VClass()>=12) myNvx12++;
      if(vertex->SumPt()>sumPt_max)
	{
	  sumPt_max=vertex->SumPt();
	  zvx=vertex->Z();
	}
    }
  return myNvx12;
}

//________ prints out event info
void TStnMetModel::PrintEvent() {
  if(DebugMode==1)
    {
      std::cout<<"From TStnMetModel::PrintEvent(): --------- Debugging Mode -----------------"<<std::endl;
      std::cout<<"________________________________________________________________________"<<std::endl;
      std::cout<<"From TStnMetModel::PrintEvent():    Run="<<fHeaderBlock->RunNumber()<<"  Event="<<fHeaderBlock->EventNumber()<<std::endl;
      std::cout<<"From TStnMetModel::PrintEvent(): raw SumEt="<<allstuff.mySumEt_raw<<std::endl;
      std::cout<<"From TStnMetModel::PrintEvent(): raw MET0="<<allstuff.myMET0_raw.Mod()<<"  MET0_phi="<<TVector2::Phi_0_2pi(allstuff.myMET0_raw.Phi())<<std::endl;
      std::cout<<"From TStnMetModel::PrintEvent(): raw MET4="<<allstuff.myMET_raw.Mod()<<"  MET4_phi="<<TVector2::Phi_0_2pi(allstuff.myMET_raw.Phi())<<std::endl;
      std::cout<<"From TStnMetModel::PrintEvent(): "<<std::endl;
      std::cout<<"From TStnMetModel::PrintEvent(): unclustered SumEt="<<jet04stuff.mySumEtCorr_th15<<std::endl;
      std::cout<<"From TStnMetModel::PrintEvent(): jet sumEt="<<jet04stuff.mySumEtJet_th15<<std::endl;
      std::cout<<"From TStnMetModel::PrintEvent(): corrected MET="<<jet04stuff.myMETcorr_th15.Mod()
	       <<"  MET_phi="<<TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())<<std::endl;

      std::cout<<"From TStnMetModel::PrintEvent(): --- printing photons: Npho="<<allstuff.myCorrPhoton.size()<<std::endl;
      //_________________ printing photons
      for(int i=0; i<allstuff.myCorrPhoton.size(); i++)
	{
	  std::cout<<"From TStnMetModel::PrintEvent(): --- cor Photon-"<<i<<": rawEt, Et, Phi, Eta, EtaDet "
		   <<allstuff.myRawPhoton[i].Pt()<<" "<<allstuff.myCorrPhoton[i].Pt()<<" "<<TVector2::Phi_0_2pi(allstuff.myCorrPhoton[i].Phi())
		   <<" "<<allstuff.myCorrPhoton[i].Eta()<<" "<<allstuff.myPhoEtaDet[i]<<std::endl;
	  std::cout<<"_______________ List of towers for Photon-"<<i<<std::endl;
	  PrintPhoTowers(allstuff.myPhoInd[i],fPhotonBlock,fCalData);
	}
      //_________________ printing electrons
      std::cout<<"From TStnMetModel::PrintEvent(): --- printing electrons: Nele="<<allstuff.myCorrElectron.size()<<std::endl;
      for(int i=0; i<allstuff.myCorrElectron.size(); i++)
	{
	  std::cout<<"From TStnMetModel::PrintEvent(): --- cor Electron-"<<i<<": rawEt, Et, Phi, Eta, EtaDet "
		   <<allstuff.myRawElectron[i].Pt()<<" "<<allstuff.myCorrElectron[i].Pt()<<" "<<TVector2::Phi_0_2pi(allstuff.myCorrElectron[i].Phi())
		   <<" "<<allstuff.myCorrElectron[i].Eta()<<" "<<allstuff.myEleEtaDet[i]<<std::endl;
	  std::cout<<"_______________ List of towers for Electron-"<<i<<std::endl;
	  PrintEleTowers(allstuff.myEleInd[i],fElectronBlock,fCalData);
	}
      //_________________ printing jets
      std::cout<<"From TStnMetModel::PrintEvent(): --- printing jets: Njet="<<jet04stuff.myNjet<<std::endl;
      std::cout<<"___ printing corrected jets after removing EM object ____"<<std::endl;
      for(int i=0; i<jet04stuff.myNjet; i++)
	{
	  std::cout<<"...... jet-"<<i+1
		   <<" Et="<<jet04stuff.Jet_lev6_noEMobj[i].Pt()
		   <<" Phi="<<TVector2::Phi_0_2pi(jet04stuff.Jet_lev6_noEMobj[i].Phi())
		   <<" eta="<<jet04stuff.Jet_lev6_noEMobj[i].Eta()
		   <<" cor_eta_det="<<jet04stuff.EtaDetCorr[i]
		   <<" raw_em_fr="<<jet04stuff.EmFrRaw[i]
		   <<" match_ID="<<jet04stuff.Nobj_match[i]<<std::endl;
	  std::cout<<"_______________ List of towers for Jet-"<<i+1<<std::endl;
	  PrintJetTowers(jet04stuff.JetBlockInd[i],fJetBlockClu04,fCalData);
	}
      std::cout<<"___ printing raw jets after removing EM object ____"<<std::endl;
      for(int i=0; i<jet04stuff.myNjet; i++)
	{
	  std::cout<<"...... jet-"<<i+1
		   <<" Et="<<jet04stuff.Jet_raw_noEMobj[i].Pt()
		   <<" Phi="<<TVector2::Phi_0_2pi(jet04stuff.Jet_raw_noEMobj[i].Phi())
		   <<" eta="<<jet04stuff.Jet_raw_noEMobj[i].Eta()
		   <<" cor_eta_det="<<jet04stuff.EtaDetCorr[i]
		   <<" raw_em_fr="<<jet04stuff.EmFrRaw[i]<<std::endl;
	}
      
      std::cout<<"___ printing corrected jets before removing EM object ____"<<std::endl;
      for(int i=0; i<jet04stuff.myNjet; i++)
	{
	  std::cout<<"...... jet-"<<i+1
		   <<" Et="<<jet04stuff.Jet_lev6[i].Pt()
		   <<" Phi="<<TVector2::Phi_0_2pi(jet04stuff.Jet_lev6[i].Phi())
		   <<" eta="<<jet04stuff.Jet_lev6[i].Eta()
		   <<" raw_eta_det="<<jet04stuff.EtaDet[i]
		   <<" raw_em_fr="<<jet04stuff.EmFrRaw[i]<<std::endl;
	}
      std::cout<<"___ printing raw jets before removing EM object ____"<<std::endl;
      for(int i=0; i<jet04stuff.myNjet; i++)
	{
	  std::cout<<"...... jet-"<<i+1
		   <<" Et="<<jet04stuff.Jet_raw[i].Pt()
		   <<" Phi="<<TVector2::Phi_0_2pi(jet04stuff.Jet_raw[i].Phi())
		   <<" eta="<<jet04stuff.Jet_raw[i].Eta()
		   <<" raw_eta_det="<<jet04stuff.EtaDet[i]
		   <<" raw_em_fr="<<jet04stuff.EmFrRaw[i]<<std::endl;
	}
      
      std::cout<<"........................................................................"<<std::endl;
      std::cout<<"From TStnMetModel::PrintEvent(): --- printing MET-significance ------------"<<std::endl;
      std::cout<<"________ MET-significance: default="<<metsigstuff.metsig_Def<<std::endl;
      for(int i=0; i<MetsigArray; i++)
	{
	  std::cout<<"________ MET-significance: syst="<<metsigstuff.metsig_syst[i]<<std::endl;
	}

      std::cout<<"........................................................................"<<std::endl;
      std::cout<<"From TStnMetModel::PrintEvent(): --- Printing  Pseudo-experiments ---------"<<std::endl;
      for(int i=0; i<metsigstuff.metsig_DefPseudo.size(); i++)
	{
	  std::cout<<"________ MET-sig Pseudo-exp: default="<<metsigstuff.metsig_DefPseudo[i]<<std::endl;
	  std::cout<<"________ MET Pseudo-exp: default="<<metsigstuff.MetPseudo[i].Mod()<<std::endl;
	  for(int j=0; j<MetsigArray; j++)
	    {
	      std::cout<<"________ Pseudo-exp syst-"<<j<<": metsig="<<metsigstuff.metsig_systPseudo[j][i]
		       <<"  MET="<<metsigstuff.MetPseudoSyst[j][i].Mod()
		       <<"  phi="<<metsigstuff.MetPseudoSyst[j][i].Phi()<<std::endl;
	    }

	}

      std::cout<<"__________ End of Event _________________________________________________"<<std::endl;
    }
  return;
}


//_________________________ clears parameters which do not depend on Jet Cone size
void TStnMetModel::ClearCommonStuff(CommonStuff &stuff) { // clears CommonStuff structures 
  stuff.myRawPhoton.clear();    
  stuff.myCorrPhoton.clear();
  stuff.myPhoInd.clear();
  stuff.myPhoEtaDet.clear();
  stuff.myRawElectron.clear();
  stuff.myEleInd.clear();
  stuff.myCorrElectron.clear(); 
  stuff.myEleEtaDet.clear();
  stuff.myElePhiDet.clear();
  stuff.myElePprEt.clear();
  stuff.myRawMuon.clear();      
  stuff.myCorrMuon.clear();
  stuff.myCalMuon.clear(); 
  stuff.myMuonZvx.clear(); 
  stuff.myMuonAlgo.clear();
  stuff.myMuonTrkInd.clear(); 
  stuff.myTrkStereo.clear(); 
  stuff.myTrkAxial.clear(); 
  stuff.myMET_raw.Set(0.0,0.0); 
  stuff.myMET0_raw.Set(0.0,0.0); 
  stuff.mySumEt_raw=0.0;

  return;
}

//__________________________ clears the JetStuff structures
void TStnMetModel::ClearJetStuff(JetStuff &stuff) { 
  
  stuff.myMETcorr_th15.Set(-1.0E6,-1.0E6);
  stuff.mySumEtCorr_th15=-1.0E6;
  stuff.mySumEtJet_th15=-1.0E6;
  stuff.myNjet=0; 
  stuff.myNjet_th15=0;
  stuff.Jet_raw.clear();      
  stuff.Jet_lev1.clear();     
  stuff.Jet_lev4.clear();     
  stuff.Jet_lev5.clear();     
  stuff.Jet_lev6.clear();     
  stuff.Jet_lev7.clear();     
  stuff.Jet_raw_noEMobj.clear();      
  stuff.Jet_lev1_noEMobj.clear();     
  stuff.Jet_lev4_noEMobj.clear();     
  stuff.Jet_lev5_noEMobj.clear();     
  stuff.Jet_lev6_noEMobj.clear();     
  stuff.Jet_lev7_noEMobj.clear();
  stuff.newJetLev6.clear();     
  stuff.JetNtrk.clear();  
  stuff.JetNtwr.clear();  
  stuff.EtaDet.clear();   
  stuff.EtaDetCorr.clear();   
  stuff.EmFrRaw.clear();  
  stuff.Nobj_match.clear(); 
  stuff.Npho_match.clear(); 
  stuff.Nele_match.clear(); 
  stuff.Nmu_match.clear(); 
  stuff.JetBlockInd.clear();
  return; 		 
}

//__________________________ clears external input
void TStnMetModel::ClearExternalInput() {
  input_Nele=0;
  input_EleInd.clear();
  input_CorrEle.clear(); 
  input_Nmuon=0;
  input_CorrMuon.clear();
  input_RawMuon.clear();
  input_trkInd.clear();
  input_TrkAlgoID.clear();
  input_TrkZvx.clear();
  input_CalCluster.clear(); 

  return;
}

//____________________ clears MetsigStuff
void TStnMetModel::ClearMetsigStuff(MetsigStuff &sigstuff) {
  sigstuff.metsig_Def=0.0;
  sigstuff.metsig_DefPseudo.clear();
  sigstuff.MetPseudo.clear();
  for(int i=0; i<MetsigArray; i++)
    {
      sigstuff.metsig_syst[i]=0.0;
      sigstuff.metsig_systPseudo[i].clear();
      sigstuff.MetPseudoSyst[i].clear();
    }
  return;
}


//___________________________________________ ****** external setters for intput params
void TStnMetModel::SetNele(int nele)
{
  input_Nele=nele;
  return;
}
void TStnMetModel::SetEleInd(int i, int ind)
{
  if(i<input_EleInd.size()) input_EleInd[i]=ind;
  else input_EleInd.push_back(ind);
  return;
}
void TStnMetModel::SetCorrEle(int i, TLorentzVector vec)
{
  if(i<input_CorrEle.size()) input_CorrEle[i]=vec;
  else input_CorrEle.push_back(vec);
  return;
}
void TStnMetModel::SetNmuon(int nmu)
{
  input_Nmuon=nmu;
  return;
}
void TStnMetModel::SetCorrMuon(int i, TLorentzVector vec)
{
  if(i<input_CorrMuon.size()) input_CorrMuon[i]=vec;
  else input_CorrMuon.push_back(vec);
  return;
}
void TStnMetModel::SetRawMuon(int i, TLorentzVector vec)
{
  if(i<input_RawMuon.size()) input_RawMuon[i]=vec;
  else input_RawMuon.push_back(vec);
  return;
}
void TStnMetModel::SetTrkInd(int i, int ind)
{
  if(i<input_trkInd.size()) input_trkInd[i]=ind;
  else input_trkInd.push_back(ind);
  return;
}
void TStnMetModel::SetTrkAlgoID(int i, int id_ind)
{
  if(i<input_TrkAlgoID.size()) input_TrkAlgoID[i]=id_ind;
  else input_TrkAlgoID.push_back(id_ind);
  return;
}
void TStnMetModel::SetTrkZvx(int i, double zvx) 
{
  if(i<input_TrkZvx.size()) input_TrkZvx[i]=zvx;
  else input_TrkZvx.push_back(zvx);
  return;
}
void TStnMetModel::SetCalCluster(int i, TLorentzVector vec)
{
  if(i<input_CalCluster.size()) input_CalCluster[i]=vec;
  else input_CalCluster.push_back(vec);
  return;
} 


//______________________________________ reads raw Met, pho, ele, muon info and fills CommonStuff
void TStnMetModel::DoCommonStuff(CommonStuff &miscstuff) {

  //_______________________________________ accessing MetBlock info
  if(!fMetBlock) return;
  miscstuff.mySumEt_raw=fMetBlock->Sumet(2);
  miscstuff.myMET_raw.SetMagPhi(fMetBlock->Met(4),fMetBlock->MetPhi(4));
  miscstuff.myMET0_raw.SetMagPhi(fMetBlock->Met(0),fMetBlock->MetPhi(0));

  //>>>> Relative order of the following methods is important
  //     1) do photons
  //     2) do electrons
  //     3) do muons/tracks
  //     4) do jets
  //     5) do MET
  //     6) do MetModel part
  //______________________________________ accessing Photon information
  DoPhotons(miscstuff); 
  //______________________________________ accessing Electron information
  DoElectrons(miscstuff);
  AddExternalElectrons(miscstuff);
  //______________________________________ accessing Muon/Track information
  DoMuons(miscstuff);
  AddExternalMuons(miscstuff);


  return;
}

//______________________________________ accessing Photon information
void TStnMetModel::DoPhotons(CommonStuff &miscstuff)
{
  if(!fPhotonBlock) return;
  TPhotonUtil::CorrectPhotonIsolation(fHeaderBlock->GetEvent());
  int phoCor_stat=TStntuple::CorrectPhotonEnergy(fHeaderBlock->GetEvent(),Zvx_best);
  int Npho=fPhotonBlock->NPhotons();
  std::vector<int> ind_passed;
  ind_passed.clear();
  for(int i=0; i<Npho; i++)
    {
      TStnPhoton* _Pho=fPhotonBlock->Photon(i);
      if(CEMPhotonIDcut(_Pho)==1) ind_passed.push_back(i); // found CEM photon
    }
  // re-ordering photons, if needed
  int reorder_word=0;
  int reorder_word_last=0;
  int Npho_passed=ind_passed.size();
  for(int j=0; j<Npho_passed-1; j++)
    {
      reorder_word_last=reorder_word;
      for(int i=1; i<Npho_passed; i++)
	{
	  TStnPhoton* _myPhoThis = fPhotonBlock->Photon(ind_passed[i]);
	  TStnPhoton* _myPhoLast = fPhotonBlock->Photon(ind_passed[i-1]);   
	  if(_myPhoThis->Etc()>_myPhoLast->Etc()) 
	    {
	      int serv_ind=ind_passed[i-1];
	      ind_passed[i-1]=ind_passed[i];
	      ind_passed[i]=serv_ind;
	      reorder_word++;
	    }
	}
      if(reorder_word_last==reorder_word) break;
    }
  for(int i=0; i<Npho_passed; i++)
    {
      TStnPhoton* _Pho=fPhotonBlock->Photon(ind_passed[i]);
      TLorentzVector _pho_raw(0.0,0.0,0.0,0.0);
      TLorentzVector _pho_cor(0.0,0.0,0.0,0.0);
      _pho_cor.SetPx(_Pho->Momentum()->Px());
      _pho_cor.SetPy(_Pho->Momentum()->Py());
      _pho_cor.SetPz(_Pho->Momentum()->Pz());
      _pho_cor.SetE(_Pho->Momentum()->E());
      _pho_raw.SetPtEtaPhiM(_Pho->Et(),_Pho->Eta(),_Pho->Phi(),0.0);
      miscstuff.myPhoInd.push_back(i);
      miscstuff.myRawPhoton.push_back(_pho_raw);
      miscstuff.myCorrPhoton.push_back(_pho_cor);
      miscstuff.myPhoEtaDet.push_back(_Pho->DetEta());
    }

  return;
}

//_____________________________________ CEM photon ID
int TStnMetModel::CEMPhotonIDcut(TStnPhoton* Pho)
{
  int passcode=1;
  //______________________________________ only CEM photons with Et>7 GeV
  if((Pho->Detector())!=0 || (Pho->Etc())<fPhoMinEt) return 0;
  //_______________________________________ Fiducial cuts
  if(fabs(Pho->XCes())<fCemPhoMinCesX || fabs(Pho->XCes())>fCemPhoMaxCesX) return 0;
  if(fabs(Pho->ZCes())<fCemPhoMinCesZ || fabs(Pho->ZCes())>fCemPhoMaxCesZ) return 0;
  //______________________________________ HADEM cut using 3 towers
  float cutMax_HADEM=fCemPhoMaxHADEM[1]+fCemPhoMaxHADEM[2]*(Pho->Momentum()->E());
  if(cutMax_HADEM<fCemPhoMaxHADEM[0]) cutMax_HADEM=fCemPhoMaxHADEM[0];
  if((Pho->HadEm())>cutMax_HADEM) return 0;
  //______________________________________ CES Chi^2 cut 
  if((Pho->Chi2())>fCemPhoMaxCesChi2) return 0;
  //______________________________________ N3D cut
  if((Pho->N3d())>fCemPhoMaxN3D) return 0;
  //______________________________________ cut on max Pt of track in cluster if N3D=1
  if((Pho->N3d())>0)
    {
      float trkPtcut_max=fCemPhoMaxPtmax[0]+fCemPhoMaxPtmax[1]*(Pho->Etc());
      if((Pho->Pt())>trkPtcut_max) return 0;
    }
  //______________________________________ CalIso4 cut
  float cutMax_CalIso=0.0;
  if((Pho->Etc())<20.0) cutMax_CalIso=fCemPhoMaxCalIso4_corr[0]*(Pho->Etc());
  else cutMax_CalIso=fCemPhoMaxCalIso4_corr[1]+fCemPhoMaxCalIso4_corr[2]*(Pho->Etc()-20.0);
  if((Pho->EIso4(2))>cutMax_CalIso) return 0; // test
  //______________________________________ TrkIso4 cut
  if((Pho->SumPt4())>(fCemPhoMaxTrkIso4[0]+fCemPhoMaxTrkIso4[1]*(Pho->Etc()))) return 0;
  //______________________________________ Energy of 2nd CES cluster (Wire & Strip)
  float _Et2ndCES=0.0;
  if((Pho->CesStripE2())>(Pho->CesWireE2())) _Et2ndCES=(Pho->CesStripE2())*(Pho->SinTheta());
  else _Et2ndCES=(Pho->CesWireE2())*(Pho->SinTheta());
  float cutMax_2ndCes=0.0;  
  if((Pho->Etc())<18.0) cutMax_2ndCes=fCemPhoMax2ndCesEt[0]*(Pho->Etc());
  else cutMax_2ndCes=fCemPhoMax2ndCesEt[1]+fCemPhoMax2ndCesEt[2]*(Pho->Etc());
  if(_Et2ndCES>cutMax_2ndCes) return 0;

  return passcode;
}

//_______________________________ fills out electrons _____________________________
void TStnMetModel::DoElectrons(CommonStuff &miscstuff) {
  int eleCor_stat=0;
  eleCor_stat=TStntuple::CorrectElectronEnergy(fHeaderBlock->GetEvent(),Zvx_best);
  if(!fTrackBlock) return;
  if(!fElectronBlock) return;
  int Nele=fElectronBlock->NElectrons();
  std::vector<int> ind_passed;
  ind_passed.clear();
  for(int i=0; i<Nele; i++)
    {
      TStnElectron* _Ele=fElectronBlock->Electron(i);
      int _trk_ind=_Ele->TrackNumber();
      TStnTrack* _myTrk = NULL; // this needs a double-check; potential source of problems
      if(_trk_ind>-1) _myTrk = fTrackBlock->Track(_trk_ind); // getting "regular" track
      if(CEMElectronLooseIDcut(_Ele,_myTrk)==1) 
	{
	  // to prevent ele==pho double count
	  if(IsPhoEle(_Ele)==0) ind_passed.push_back(i); // found CEM electron
	}
      if(PEMElectronLooseIDcut(_Ele)==1) 
	{
	  // to prevent ele==pho double count
	  if(IsPhoEle(_Ele)==0) ind_passed.push_back(i); // found PEM electron
	}      
    }
  //_______________________________________________________________________________
  // re-ordering electrons, if needed
  int reorder_word=0;
  int reorder_word_last=0;
  int Nele_passed=ind_passed.size();
  for(int j=0; j<Nele_passed-1; j++)
    {
      reorder_word_last=reorder_word;
      for(int i=1; i<Nele_passed; i++)
	{
	  TStnElectron* _myEleThis = fElectronBlock->Electron(ind_passed[i]);
	  TStnElectron* _myEleLast = fElectronBlock->Electron(ind_passed[i-1]);   
	  if(myCorrEtEle(_myEleThis)>myCorrEtEle(_myEleLast)) 
	    {
	      int serv_ind=ind_passed[i-1];   // sorting index
	      ind_passed[i-1]=ind_passed[i];
	      ind_passed[i]=serv_ind;
	      reorder_word++;
	    }
	}
      if(reorder_word_last==reorder_word) break;
    }
  //___________________________________________________________________________________________
  for(int i=0; i<Nele_passed; i++)
    {
      TStnElectron* _myEle = fElectronBlock->Electron(ind_passed[i]);
      const TLorentzVector* ele_momentumvec = _myEle->Momentum(); 

      miscstuff.myEleInd.push_back(ind_passed[i]);
      miscstuff.myEleEtaDet.push_back(_myEle->DetEta());
      miscstuff.myElePhiDet.push_back(_myEle->EmClusPhi());

      double e_mass=0.000511;
      double e_eta=ele_momentumvec->Eta();
      double e_phi=ele_momentumvec->Phi();
      double e_pt=myCorrEtEle(_myEle);
      double clus_eta=_myEle->EmClusEvEta();
      double clus_phi=_myEle->EmClusPhi();
      double clus_pt=_myEle->TotalEt();
      TLorentzVector vec;
      //____________________ corrected electron
      vec.SetPtEtaPhiM(e_pt,e_eta,e_phi,e_mass);
      miscstuff.myCorrElectron.push_back(vec); 
      //____________________ raw electron (to be the same as raw photon)
      vec.SetPtEtaPhiM(clus_pt,clus_eta,clus_phi,0.0);
      miscstuff.myRawElectron.push_back(vec);      

      if(_myEle->DetectorCode()==0) miscstuff.myElePprEt.push_back(0.0);
      if(_myEle->DetectorCode()==1) 
	{
	  double factor=_myEle->Etcor();
	  double extra_PemEleEt=(factor-1.0)*_myEle->TotalEt();
	  double pprEt=extra_PemEleEt+factor*(_myEle->PprEnergy())*sin(miscstuff.myRawElectron[i].Theta());
	  if(pprEt<0.0) pprEt=0.0; 
	  miscstuff.myElePprEt.push_back(pprEt);
	}
    }
  
  return;
}

//________________ adds external electrons to the list
void TStnMetModel::AddExternalElectrons(CommonStuff &miscstuff) {
  if(input_EleInd.size()!=input_CorrEle.size()) 
    {
      std::cout<<"WARNING from TStnMetModel::AddExternalElectrons --->>> You are not passing correct electron information"<<std::endl;
      return;
    }
  else
    {
      if(input_EleInd.size()>0)
	{
	  int match_stat=0; // 0=all matched; >0 mismatch
	  std::vector<int> input_match;
	  input_match.clear();
	  
	  for(int i=0; i<input_EleInd.size(); i++)
	    {
	      input_match.push_back(0); // assign zero before matching
	      for(int j=0; j<miscstuff.myCorrElectron.size(); j++)
		{
		  double dR=miscstuff.myCorrElectron[j].DeltaR(input_CorrEle[i]);
		  if(dR<0.2 && input_EleInd[i]==miscstuff.myEleInd[j])
		    {
		      match_stat++;
		      input_match[i]=1; // re-assign 1 if matched
		    }
		}
	    }
	  if(match_stat!=input_EleInd.size())
	    {
	      for(int i=0; i<input_match.size(); i++)
		{
		  int Nele=fElectronBlock->NElectrons();
		  if(input_match[i]==0 && input_EleInd[i]<Nele)
		    {
		      TStnElectron* _myEle = fElectronBlock->Electron(input_EleInd[i]);
		      if(IsPhoEle(_myEle)==0) // to prevent ele==pho double-count
			{
			  miscstuff.myCorrElectron.push_back(input_CorrEle[i]);
			  miscstuff.myEleInd.push_back(input_EleInd[i]);
			  miscstuff.myEleEtaDet.push_back(_myEle->DetEta());
			  miscstuff.myElePhiDet.push_back(_myEle->EmClusPhi());
			  
			  double clus_eta=_myEle->EmClusEvEta();
			  double clus_phi=_myEle->EmClusPhi();
			  double clus_pt=_myEle->TotalEt();
			  TLorentzVector vec;
			  vec.SetPtEtaPhiM(clus_pt,clus_eta,clus_phi,0.0);
			  miscstuff.myRawElectron.push_back(vec);      
			  if(_myEle->DetectorCode()==0) miscstuff.myElePprEt.push_back(0.0);
			  if(_myEle->DetectorCode()==1) 
			    {
			      double factor=_myEle->Etcor();
			      double extra_PemEleEt=(factor-1.0)*_myEle->TotalEt();
			      double pprEt=extra_PemEleEt+factor*(_myEle->PprEnergy())*sin(vec.Theta());
			      if(pprEt<0.0) pprEt=0.0; 
			      miscstuff.myElePprEt.push_back(pprEt);
			    }
			}
		    }
		}
	    }
	}
    }
  return; 
}

//______ prevents double count of ele==pho: returns 1 if ele==pho, 0--otherwise
int TStnMetModel::IsPhoEle(TStnElectron* Ele) {
  int passcode=0;
  if(allstuff.myRawPhoton.size()>0)
    {
      TStnPhoton* myPho=TPhotonUtil::MatchPhoton(fHeaderBlock->GetEvent(),Ele);
      TLorentzVector pho_vec(0.0,0.0,0.0,0.0);
      pho_vec.SetPtEtaPhiM(myPho->Et(),myPho->Eta(),myPho->Phi(),0.0);
      for(int j=0; j<allstuff.myRawPhoton.size(); j++)
	{
	  double dr=pho_vec.DeltaR(allstuff.myRawPhoton[j]);
	  if(dr<0.2) return 1;
	}
    }
  return passcode;
}


//_______________________________ loose CEM electron ID ____________________________ 
int TStnMetModel::CEMElectronLooseIDcut(TStnElectron* Ele, TStnTrack* Trk) {
  int passcode=1;
  if(Trk==NULL) return 0;
  if((Ele->DetectorCode())!=0) return 0;
  if((Ele->FidEleSmx())!=fCemEleFidCut) return 0;
  if((Ele->TrackNumber())==-1) return 0; // these are junk CEM electrons
  if((myCorrEtEle(Ele))<fEleMinEt) return 0;
  if(fabs(Ele->Z0())>fCemEleMaxTrkZ0) return 0; // note, not beam-constrained track !!!
  if((Trk->Algorithm())==2)
    {
      if((Ele->TrackBcPt())<fEleTrkMinPt) return 0; // note, beam-constrained track !!! 
    }
  else
    {
      if((Ele->TrackPt())<fEleTrkMinPt) return 0; // note, not beam-constrained track !!! 
    }
  if((Ele->Nasl())<fCemEleMinCotAxSl) return 0;
  if((Ele->Nssl())<fCemEleMinCotStSl) return 0;
  if((Ele->HadEm())>(fCemEleMaxHADEM[0]+fCemEleMaxHADEM[1]*(myCorrEnergyEle(Ele)))) return 0;
  if((myCorrIsoFr(Ele))>fCemEleMaxIsoFr) return 0;
  
  return passcode;
}

//_______________________________ loose PEM electron ID ____________________________ 
int TStnMetModel::PEMElectronLooseIDcut(TStnElectron* Ele) {
  int passcode=1;
  double Et_em=myCorrEtEle(Ele);
  if((Ele->DetectorCode())!=1) return 0;
  if(Et_em<fEleMinEt) return 0;   
  if(fabs(Ele->PesEta())<fPemEleMinEtaPes || fabs(Ele->PesEta())>fPemEleMaxEtaPes) return 0;
  if((Ele->HadEm())>fPemEleMaxHADEM) return 0;
  if((Ele->Pem3x3FitTower())==fPemEle3x3FitCut) return 0;
  if((Ele->Chi2Three())>fPemEleMax3x3Chi2) return 0;
  if((Ele->Pes5x9(0))<fPemEleMinPes5x9U) return 0;
  if((Ele->Pes5x9(1))<fPemEleMinPes5x9V) return 0;
  if((myCorrIsoFr(Ele))>fPemEleMaxIsoFr) return 0;
  TStnPhoton* myPho=TPhotonUtil::MatchPhoton(fHeaderBlock->GetEvent(),Ele);
  if((myPho->PesDeltaR())>fPemEleMaxPesdR) return 0;

  return passcode;
}

//____________________________________________________________________________
// returns Iso corrected for leakage & nvx12
double TStnMetModel::myCorrIso(TStnElectron* ele) {
  double Iso=ele->Iso();
  Iso=Iso - (myNvx_class12>1? 0.3563*(myNvx_class12-1) : 0.0) - ele->LeakageEnergy();
  return Iso;
}
//____________________________________________________________________________
// returns Iso/Et(ele) corrected for leakage & nvx12
double TStnMetModel::myCorrIsoFr(TStnElectron* ele) {
  return myCorrIso(ele)/ele->Et();
}
//____________returns corrected E of electrons (includes all corrections)_________ 
double TStnMetModel::myCorrEnergyEle(TStnElectron* ele) {
  return (ele->Et()>0.0? ele->E()*myCorrEtEle(ele)/(ele->Et()) : ele->E());
}

//____________returns corrected Et of electrons (includes all corrections)_________ 
double TStnMetModel::myCorrEtEle(TStnElectron* ele) {
  double et=ele->Et();
  double factor=1.0;
  factor=ele->Etcor(); // provide that fEtcor has already been corrected by Ray's CorrectElectronEnergy()
  et *= factor;
  return et;
}

//______________ fills out muons
void TStnMetModel::DoMuons(CommonStuff &miscstuff) {
  // empty for a moment
  if(!fTrackBlock) return;
  return; 
}

//______________ adds external muons to the list
void TStnMetModel::AddExternalMuons(CommonStuff &miscstuff) {
  if(input_Nmuon!=input_CorrMuon.size() || input_trkInd.size()!=input_CorrMuon.size()) 
    {
      std::cout<<"WARNING from TMetModel::AddExternalMuons --->>> You are not passing correct muon information"<<std::endl;
      return;
    }
  else
    {
      if(input_trkInd.size()>0 || miscstuff.myCorrMuon.size()>0)
	{
	  int match_stat=0; // 0=all matched; >0 mismatch
	  std::vector<int> input_match;
	  input_match.clear();
	  
	  for(int i=0; i<input_trkInd.size(); i++)
	    {
	      input_match.push_back(0); // assign zero before matching
	      for(int j=0; j<miscstuff.myCorrMuon.size(); j++)
		{
		  double dR=miscstuff.myCorrMuon[j].DeltaR(input_CorrMuon[i]);
		  if(dR<0.2)
		    {
		      match_stat++;
		      input_match[i]=1; // re-assign 1 if matched
		    }
		}
	    }
	  if(match_stat!=input_trkInd.size())
	    {
	      for(int i=0; i<input_match.size(); i++)
		{
		  if(input_match[i]==0)
		    {
// 		      miscstuff.myRawMuon.push_back(input_RawMuon[i]);
		      miscstuff.myCorrMuon.push_back(input_CorrMuon[i]);
		      miscstuff.myCalMuon.push_back(input_CalCluster[i]);
// 		      miscstuff.myMuonZvx.push_back(input_TrkZvx[i]);
// 		      miscstuff.myMuonAlgo.push_back(input_TrkAlgoID[i]);
		      miscstuff.myMuonTrkInd.push_back(input_trkInd[i]);
		      TStnTrack* _myTrk = fTrackBlock->Track(input_trkInd[i]);
		      int COTSeg_St=_myTrk->NCotStSeg(7); 
		      int COTSeg_Ax=_myTrk->NCotAxSeg(7);
		      miscstuff.myTrkStereo.push_back(COTSeg_St); 
		      miscstuff.myTrkAxial.push_back(COTSeg_Ax);
		    }
		}
	    }
	}
    }
  
  return;
}


//----------------------------------------------------------------------------------------
//
// >>>>>>>>>>>>>>>>>>>>>  HERE come JETS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
//
//________________________________________________________________________________________
//__________________________________ does my jets; main routine  
void TStnMetModel::DoMyJet(TStnJetBlock* fJetBlock, CommonStuff miscstuff, JetStuff &jetstuff) {
  myJetRun=fHeaderBlock->RunNumber(); // obtaining run number, myJetRun is globaly defined 
  fJTC_coneSize=0; // cone 0.4
  fJTC_imode = 1-fHeaderBlock->McFlag(); // 1==data; 0==MC
  fJTC_version=5;
  DoMyJetNoMatch(fJetBlock,jetstuff);
  DoMyJetWithMatch(fJetBlock,miscstuff,jetstuff);
  ReorderMyJets(jetstuff);
  MyNewJetWithMet(jetstuff); 
  return;
}

//__________________________________ does my jets after removing EM objetcs, to be called after DoMyJetNoMatch
void TStnMetModel::DoMyJetWithMatch(TStnJetBlock* fJetBlock, CommonStuff miscstuff, JetStuff &jetstuff) {

  //_________________________ Initializing Jet corrections
  CorInit lev1;
  CorInit lev4;
  CorInit lev5;
  CorInit lev6;
  CorInit lev7;
  lev1.level=1;
  lev1.nvx=myNvx_class12;
  lev1.cone=fJTC_coneSize;
  lev1.version=fJTC_version;
  lev1.sys=fJTC_systcode; 
  lev1.imode=fJTC_imode;  
  lev1.Nrun=myJetRun;
  lev4.level=4;
  lev4.nvx=myNvx_class12;
  lev4.cone=fJTC_coneSize;
  lev4.version=fJTC_version;
  lev4.sys=fJTC_systcode;  
  lev4.imode=fJTC_imode;  
  lev4.Nrun=myJetRun;
  lev5.level=5;
  lev5.nvx=myNvx_class12;
  lev5.cone=fJTC_coneSize;
  lev5.version=fJTC_version;
  lev5.sys=fJTC_systcode;
  lev5.imode=fJTC_imode;    
  lev5.Nrun=myJetRun; 
  lev6.level=6;
  lev6.nvx=myNvx_class12;
  lev6.cone=fJTC_coneSize;
  lev6.version=fJTC_version;
  lev6.sys=fJTC_systcode;
  lev6.imode=fJTC_imode;  
  lev6.Nrun=myJetRun;
  lev7.level=7;
  lev7.nvx=myNvx_class12;
  lev7.cone=fJTC_coneSize;
  lev7.version=fJTC_version;
  lev7.sys=fJTC_systcode;
  lev7.imode=fJTC_imode;  
  lev7.Nrun=myJetRun;

  int _Npho_match=0;
  int _Nele_match=0;
  double epsilon=1.0E-10;
  for(int i=0; i<jetstuff.myNjet; i++)
    {
      float _etadet_tmp=jetstuff.EtaDet[i]; // new line, 10/20/05 
      TLorentzVector *_jet=&jetstuff.Jet_raw[i];
      //_____ removing photons
      for(int j=0; j<miscstuff.myRawPhoton.size() && _Npho_match<miscstuff.myRawPhoton.size(); j++)
	{
	  TLorentzVector *_pho=&miscstuff.myRawPhoton[j];
	  int match_stat=MyMatchPhoJet(jetstuff.JetBlockInd[i],miscstuff.myPhoInd[j],-1,fJetBlock,_pho,_jet);
	  if(match_stat==1)
	    {
	      _Npho_match++;
	      jetstuff.Nobj_match[i]=jetstuff.Nobj_match[i]+1;
	      jetstuff.Npho_match[i]=jetstuff.Npho_match[i]+1;
	      // corrected EtaDet is to be calculated before removing EM object 
	      _etadet_tmp=CorrectEtaDet(_pho,&jetstuff.Jet_raw_noEMobj[i], 
					jetstuff.EtaDetCorr[i],miscstuff.myPhoEtaDet[j]); // new line, 10/20/05
	      jetstuff.EtaDetCorr[i]=_etadet_tmp;  // new line, 10/20/05
	      if(jetstuff.Jet_raw_noEMobj[i].E()>miscstuff.myRawPhoton[j].E()) 
		jetstuff.Jet_raw_noEMobj[i]=jetstuff.Jet_raw_noEMobj[i]-miscstuff.myRawPhoton[j]; // removing raw photon from raw jet
	      else jetstuff.Jet_raw_noEMobj[i].SetPxPyPzE(epsilon,epsilon,epsilon,sqrt(3.0)*epsilon);
	      if(jetstuff.Jet_raw_noEMobj[i].Pt()<epsilon) jetstuff.Jet_raw_noEMobj[i].SetPxPyPzE(epsilon,epsilon,epsilon,sqrt(3.0)*epsilon); // new, 05/28/08
	      if(jetstuff.Jet_raw_noEMobj[i].E()<jetstuff.Jet_raw_noEMobj[i].P()
		 || jetstuff.Jet_raw_noEMobj[i].E()<=0.0
		 || (jetstuff.Jet_raw[i].DeltaR(jetstuff.Jet_raw_noEMobj[i]))>2.0*(fJetBlock->ConeSize())) 
		{
		  jetstuff.Jet_raw_noEMobj[i].SetPxPyPzE(epsilon,epsilon,epsilon,sqrt(3.0)*epsilon); // in case removing EM object leads to M()<0.0
		}
	    }
	}
      //_____ removing electrons
      for(int j=0; j<miscstuff.myRawElectron.size() && _Nele_match<miscstuff.myRawElectron.size(); j++)
	{
	  TLorentzVector *_ele=&miscstuff.myRawElectron[j];
	  int match_stat=MyMatchPhoJet(jetstuff.JetBlockInd[i],-1,miscstuff.myEleInd[j],fJetBlock,_ele,_jet);
	  if(match_stat==1)
	    {
	      _Nele_match++;
	      jetstuff.Nobj_match[i]=jetstuff.Nobj_match[i]+1;
	      jetstuff.Nele_match[i]=jetstuff.Nele_match[i]+1;
	      // corrected EtaDet is to be calculated before removing EM object 
	      _etadet_tmp=CorrectEtaDet(_ele,&jetstuff.Jet_raw_noEMobj[i], 
					jetstuff.EtaDetCorr[i],miscstuff.myEleEtaDet[j]); // new line, 10/20/05
	      jetstuff.EtaDetCorr[i]=_etadet_tmp;  // new line, 10/20/05
	      if(jetstuff.Jet_raw_noEMobj[i].E()>miscstuff.myRawElectron[j].E()) 
		jetstuff.Jet_raw_noEMobj[i]=jetstuff.Jet_raw_noEMobj[i]-miscstuff.myRawElectron[j]; // removing raw electron from raw jet
	      else jetstuff.Jet_raw_noEMobj[i].SetPxPyPzE(epsilon,epsilon,epsilon,sqrt(3.0)*epsilon);
	      if(jetstuff.Jet_raw_noEMobj[i].Pt()<epsilon) jetstuff.Jet_raw_noEMobj[i].SetPxPyPzE(epsilon,epsilon,epsilon,sqrt(3.0)*epsilon); // new, 05/28/08
	      if(jetstuff.Jet_raw_noEMobj[i].E()<jetstuff.Jet_raw_noEMobj[i].P() 
		 || jetstuff.Jet_raw_noEMobj[i].E()<=0.0
		 || (jetstuff.Jet_raw[i].DeltaR(jetstuff.Jet_raw_noEMobj[i]))>2.0*(fJetBlock->ConeSize())) 
		{
		  jetstuff.Jet_raw_noEMobj[i].SetPxPyPzE(epsilon,epsilon,epsilon,sqrt(3.0)*epsilon); // in case removing EM object leads to M()<0.0
		}
	    }
	}
      if((jetstuff.Nele_match[i]+jetstuff.Npho_match[i])>0) // correcting jet energy after removing EM objetcs
	{
	  TLorentzVector _rawJet=jetstuff.Jet_raw_noEMobj[i];
	  TLorentzVector _rawJet_tmp;
	  _rawJet_tmp=_rawJet;
	  float _myEmFr_tmp=jetstuff.EmFrRaw[i];
	  double corr7=GetCorrection(lev7,_rawJet_tmp,_myEmFr_tmp,_etadet_tmp);
	  _myEmFr_tmp=jetstuff.EmFrRaw[i]; // do this again because EmFr can be potentially changed in GetCorrection
	  _rawJet_tmp=_rawJet; // do this again because _rawJet_tmp can be potentially changed in GetCorrection
	  double corr6=GetCorrection(lev6,_rawJet_tmp,_myEmFr_tmp,_etadet_tmp);
	  _myEmFr_tmp=jetstuff.EmFrRaw[i]; // do this again because EmFr can be potentially changed in GetCorrection
	  _rawJet_tmp=_rawJet; // do this again because _rawJet_tmp can be potentially changed in GetCorrection
	  double corr5=GetCorrection(lev5,_rawJet_tmp,_myEmFr_tmp,_etadet_tmp);
	  _myEmFr_tmp=jetstuff.EmFrRaw[i]; // do this again because EmFr can be potentially changed in GetCorrection
	  _rawJet_tmp=_rawJet; // do this again because _rawJet_tmp can be potentially changed in GetCorrection
	  double corr4=GetCorrection(lev4,_rawJet_tmp,_myEmFr_tmp,_etadet_tmp);
	  _myEmFr_tmp=jetstuff.EmFrRaw[i]; // do this again because EmFr can be potentially changed in GetCorrection
	  _rawJet_tmp=_rawJet; // do this again because _rawJet_tmp can be potentially changed in GetCorrection
	  double corr1=GetCorrection(lev1,_rawJet_tmp,_myEmFr_tmp,_etadet_tmp);	  
	  if(corr1<0.0) corr1=epsilon;
	  if(corr4<0.0) corr4=epsilon;
	  if(corr5<0.0) corr5=epsilon;
	  if(corr6<0.0) corr6=epsilon;
	  if(corr7<0.0) corr7=epsilon;
	  jetstuff.Jet_lev1_noEMobj[i]=_rawJet*corr1; 
	  jetstuff.Jet_lev4_noEMobj[i]=_rawJet*corr4; 
	  jetstuff.Jet_lev5_noEMobj[i]=_rawJet*corr5; 
	  jetstuff.Jet_lev6_noEMobj[i]=_rawJet*corr6; 
	  jetstuff.Jet_lev7_noEMobj[i]=_rawJet*corr7; 	  
	}
      if(fabs(jetstuff.Jet_lev6_noEMobj[i].Eta())>=MinJetEta && fabs(jetstuff.Jet_lev6_noEMobj[i].Eta())<=MaxJetEta)
	{
	  if(jetstuff.Jet_lev6_noEMobj[i].Pt()>15.0) jetstuff.myNjet_th15++;
	}
    }
  return;
}

//__________________________________ does my jets, but doesn't remove EM objetcs
void TStnMetModel::DoMyJetNoMatch(TStnJetBlock* fJetBlock, JetStuff &jetstuff) {

  //_________________________ Initializing Jet corrections
  CorInit lev1;
  CorInit lev4;
  CorInit lev5;
  CorInit lev6;
  CorInit lev7;
  lev1.level=1;
  lev1.nvx=myNvx_class12;
  lev1.cone=fJTC_coneSize;
  lev1.version=fJTC_version;
  lev1.sys=fJTC_systcode; 
  lev1.imode=fJTC_imode;  
  lev1.Nrun=myJetRun;
  lev4.level=4;
  lev4.nvx=myNvx_class12;
  lev4.cone=fJTC_coneSize;
  lev4.version=fJTC_version;
  lev4.sys=fJTC_systcode;  
  lev4.imode=fJTC_imode;  
  lev4.Nrun=myJetRun;
  lev5.level=5;
  lev5.nvx=myNvx_class12;
  lev5.cone=fJTC_coneSize;
  lev5.version=fJTC_version;
  lev5.sys=fJTC_systcode;
  lev5.imode=fJTC_imode;    
  lev5.Nrun=myJetRun; 
  lev6.level=6;
  lev6.nvx=myNvx_class12;
  lev6.cone=fJTC_coneSize;
  lev6.version=fJTC_version;
  lev6.sys=fJTC_systcode;
  lev6.imode=fJTC_imode;  
  lev6.Nrun=myJetRun;
  lev7.level=7;
  lev7.nvx=myNvx_class12;
  lev7.cone=fJTC_coneSize;
  lev7.version=fJTC_version;
  lev7.sys=fJTC_systcode;
  lev7.imode=fJTC_imode;  
  lev7.Nrun=myJetRun;

  jetstuff.myNjet=fJetBlock->NJets();
  double epsilon=1.0E-10;
  for(int i=0; i<jetstuff.myNjet; i++)
    {
      TStnJet* jet = fJetBlock->Jet(i);
      jetstuff.JetBlockInd.push_back(i);
      jetstuff.JetNtrk.push_back(jet->NTracks());
      jetstuff.JetNtwr.push_back(jet->NTowers());
      jetstuff.EtaDet.push_back(jet->DetEta());
      jetstuff.EtaDetCorr.push_back(jet->DetEta()); // for a moment, make it the same as "EtaDet"
      jetstuff.EmFrRaw.push_back(jet->Emfr());

      jetstuff.Nobj_match.push_back(0); // these are packed with zero's for a moment 
      jetstuff.Npho_match.push_back(0);  
      jetstuff.Nele_match.push_back(0);  
      jetstuff.Nmu_match.push_back(0);   

      //__________________ getting raw jets
      TLorentzVector _rawJet;
      TLorentzVector _rawJet_tmp;
      _rawJet.SetPx(jet->Momentum()->Px());
      _rawJet.SetPy(jet->Momentum()->Py());
      _rawJet.SetPz(jet->Momentum()->Pz());
      _rawJet.SetE(jet->Momentum()->E());
      jetstuff.Jet_raw.push_back(_rawJet);
      jetstuff.Jet_raw_noEMobj.push_back(_rawJet); // for a moment, make it the same as "Jet_raw"
      
      float _myEmFr_tmp=jet->Emfr();
      float _etadet_tmp=jet->DetEta();
      _rawJet_tmp=_rawJet;
      double corr7=GetCorrection(lev7,_rawJet_tmp,_myEmFr_tmp,_etadet_tmp);
      _myEmFr_tmp=jet->Emfr(); // do this again because EmFr can be potentially changed in GetCorrection
      _rawJet_tmp=_rawJet; // do this again because _rawJet_tmp can be potentially changed in GetCorrection
      double corr6=GetCorrection(lev6,_rawJet_tmp,_myEmFr_tmp,_etadet_tmp);
      _myEmFr_tmp=jet->Emfr(); // do this again because EmFr can be potentially changed in GetCorrection
      _rawJet_tmp=_rawJet; // do this again because _rawJet_tmp can be potentially changed in GetCorrection
      double corr5=GetCorrection(lev5,_rawJet_tmp,_myEmFr_tmp,_etadet_tmp);
      _myEmFr_tmp=jet->Emfr(); // do this again because EmFr can be potentially changed in GetCorrection
      _rawJet_tmp=_rawJet; // do this again because _rawJet_tmp can be potentially changed in GetCorrection
      double corr4=GetCorrection(lev4,_rawJet_tmp,_myEmFr_tmp,_etadet_tmp);
      _myEmFr_tmp=jet->Emfr(); // do this again because EmFr can be potentially changed in GetCorrection
      _rawJet_tmp=_rawJet; // do this again because _rawJet_tmp can be potentially changed in GetCorrection
      double corr1=GetCorrection(lev1,_rawJet_tmp,_myEmFr_tmp,_etadet_tmp);

      if(corr1<0.0) corr1=epsilon;
      if(corr4<0.0) corr4=epsilon;
      if(corr5<0.0) corr5=epsilon;
      if(corr6<0.0) corr6=epsilon;
      if(corr7<0.0) corr7=epsilon;
      jetstuff.Jet_lev1.push_back(_rawJet*corr1);
      jetstuff.Jet_lev1_noEMobj.push_back(_rawJet*corr1); // for a moment, make it the same as "Jet_lev1"
      jetstuff.Jet_lev4.push_back(_rawJet*corr4);
      jetstuff.Jet_lev4_noEMobj.push_back(_rawJet*corr4); // for a moment, make it the same as "Jet_lev4"
      jetstuff.Jet_lev5.push_back(_rawJet*corr5);
      jetstuff.Jet_lev5_noEMobj.push_back(_rawJet*corr5); // for a moment, make it the same as "Jet_lev5"
      jetstuff.Jet_lev6.push_back(_rawJet*corr6);
      jetstuff.Jet_lev6_noEMobj.push_back(_rawJet*corr6); // for a moment, make it the same as "Jet_lev6"
      jetstuff.Jet_lev7.push_back(_rawJet*corr7);
      jetstuff.Jet_lev7_noEMobj.push_back(_rawJet*corr7); // for a moment, make it the same as "Jet_lev7"
            
    }
  return;
}

//_________ This function adds Met to closest jet, the output is "new" JetLev6
void TStnMetModel::MyNewJetWithMet(JetStuff &jetstuff) {
  double phi_met=jetstuff.myMETcorr_th15.Phi();
  int closejet_indJ=-1; // for jets not matched to EM objects
  int NrawJet=0;
  double dphi_min=10.0;
  double maxNewPt=0.0;
  for(int i=0; i<jetstuff.Jet_lev6_noEMobj.size(); i++) // packing "jetstuff.newJetLev6" with "old" jets 
    {
      jetstuff.newJetLev6.push_back(jetstuff.Jet_lev6_noEMobj[i]);
      NrawJet++;
      if(jetstuff.myMETcorr_th15.Mod()>0.0 && jetstuff.Jet_raw_noEMobj[i].Pt()>3.0) 
	{
	  double phi_jet=jetstuff.Jet_lev6_noEMobj[i].Phi();
	  double dphi=fabs(TVector2::Phi_mpi_pi(phi_jet-phi_met));
	  double metjetRatio=jetstuff.myMETcorr_th15.Mod()*cos(dphi)/jetstuff.Jet_lev6_noEMobj[i].Pt();
	  double low_lim1=0.0; // based on MET*cos(dPhi)/Et(lev6) studies for QCD H.F & ZH
	  double up_lim=20.0; //  based on MET*cos(dPhi)/Et(lev6) studies for QCD H.F & ZH
	  if(dphi<0.4 && metjetRatio>low_lim1 && metjetRatio<up_lim && (metjetRatio+1.0)*jetstuff.Jet_lev6_noEMobj[i].Pt()>maxNewPt)
	    {
	      maxNewPt=(metjetRatio+1.0)*jetstuff.Jet_lev6_noEMobj[i].Pt();
	      closejet_indJ=i;
	    }
	}
    }
  if(closejet_indJ<0) return;
  double phi_jet=jetstuff.Jet_lev6_noEMobj[closejet_indJ].Phi();
  double eta_jet=jetstuff.Jet_lev6_noEMobj[closejet_indJ].Eta();
  double pt_jet=jetstuff.Jet_lev6_noEMobj[closejet_indJ].Pt();
  double raw_pt_jet=jetstuff.Jet_raw_noEMobj[closejet_indJ].Pt();
  double e_jet=jetstuff.Jet_lev6_noEMobj[closejet_indJ].E();
  double dphi=fabs(TVector2::Phi_mpi_pi(phi_jet-phi_met));
  double pt_scale=(jetstuff.myMETcorr_th15.Mod()*cos(dphi)+pt_jet)/pt_jet;
  if(pt_scale<=0.0) pt_scale=1.0; 
  jetstuff.newJetLev6[closejet_indJ].SetPtEtaPhiE(pt_jet*pt_scale,eta_jet,phi_jet,e_jet*pt_scale);
  return;
}


//--the following three functions could have been replaced by one which works with generic data type (to be done)
void TStnMetModel::myExchange_tlv(TLorentzVector& val1, TLorentzVector& val2) { //exchanges val1 and val2
  TLorentzVector dummy=val1;
  val1=val2;
  val2=dummy;
  return;
}
void TStnMetModel::myExchange_dbl(double& val1, double& val2) { //exchanges val1 and val2
  double dummy=val1;
  val1=val2;
  val2=dummy;
  return;
}
void TStnMetModel::myExchange_int(int& val1, int& val2) { //exchanges val1 and val2
  int dummy=val1;
  val1=val2;
  val2=dummy;
  return;
}
//_________________________________________ reorders jets after removing EM objects
void TStnMetModel::ReorderMyJets(JetStuff &jetstuff) {
  int reorder_word=0;
  int reorder_word_last=0;
  for(int j=0; j<jetstuff.myNjet-1; j++)
    {
      reorder_word_last=reorder_word;
      for(int i=1; i<jetstuff.myNjet; i++)
	{
	  if(jetstuff.Jet_lev6_noEMobj[i].Pt()>jetstuff.Jet_lev6_noEMobj[i-1].Pt()) 
	    {
	      myExchange_tlv(jetstuff.Jet_raw[i-1],jetstuff.Jet_raw[i]);
	      myExchange_tlv(jetstuff.Jet_lev1[i-1],jetstuff.Jet_lev1[i]);
	      myExchange_tlv(jetstuff.Jet_lev4[i-1],jetstuff.Jet_lev4[i]);
	      myExchange_tlv(jetstuff.Jet_lev5[i-1],jetstuff.Jet_lev5[i]);
	      myExchange_tlv(jetstuff.Jet_lev6[i-1],jetstuff.Jet_lev6[i]);
	      myExchange_tlv(jetstuff.Jet_lev7[i-1],jetstuff.Jet_lev7[i]);
	      myExchange_tlv(jetstuff.Jet_raw_noEMobj[i-1],jetstuff.Jet_raw_noEMobj[i]);
	      myExchange_tlv(jetstuff.Jet_lev1_noEMobj[i-1],jetstuff.Jet_lev1_noEMobj[i]);
	      myExchange_tlv(jetstuff.Jet_lev4_noEMobj[i-1],jetstuff.Jet_lev4_noEMobj[i]);
	      myExchange_tlv(jetstuff.Jet_lev5_noEMobj[i-1],jetstuff.Jet_lev5_noEMobj[i]);
	      myExchange_tlv(jetstuff.Jet_lev6_noEMobj[i-1],jetstuff.Jet_lev6_noEMobj[i]);
	      myExchange_tlv(jetstuff.Jet_lev7_noEMobj[i-1],jetstuff.Jet_lev7_noEMobj[i]);
	      
	      myExchange_dbl(jetstuff.EtaDet[i-1],jetstuff.EtaDet[i]);
	      myExchange_dbl(jetstuff.EtaDetCorr[i-1],jetstuff.EtaDetCorr[i]);
	      myExchange_dbl(jetstuff.EmFrRaw[i-1],jetstuff.EmFrRaw[i]);

	      myExchange_int(jetstuff.JetNtrk[i-1],jetstuff.JetNtrk[i]);
	      myExchange_int(jetstuff.JetNtwr[i-1],jetstuff.JetNtwr[i]);
	      myExchange_int(jetstuff.Nobj_match[i-1],jetstuff.Nobj_match[i]); 
	      myExchange_int(jetstuff.Npho_match[i-1],jetstuff.Npho_match[i]); 
	      myExchange_int(jetstuff.Nele_match[i-1],jetstuff.Nele_match[i]); 
	      myExchange_int(jetstuff.Nmu_match[i-1],jetstuff.Nmu_match[i]);
	      myExchange_int(jetstuff.JetBlockInd[i-1],jetstuff.JetBlockInd[i]);

	      reorder_word++;
	    }
	}
      if(reorder_word_last==reorder_word) break;
    }  
  return;
}

//____________________________________ re-calculates jet EtaDet after removing EM object
//  this function uses Snowmass convention to calculate jet centroid (true for JetClu only)
double TStnMetModel::CorrectEtaDet(TLorentzVector *vec_pho, TLorentzVector *vec_old, 
				double jetetadet_old, double pho_etadet) {
  double E_old=vec_old->E();
  double E_pho=vec_pho->E();
  double theta_det_old=2.0*atan(exp(-jetetadet_old));
  double theta_det_pho=2.0*atan(exp(-pho_etadet));
  double etadet_new=jetetadet_old;
  double det=E_old*fabs(sin(theta_det_old))-E_pho*fabs(sin(theta_det_pho));
  double num=jetetadet_old*E_old*fabs(sin(theta_det_old))-pho_etadet*E_pho*fabs(sin(theta_det_pho));
  if(det>0.1) etadet_new=num/det;
  return etadet_new;
}

double TStnMetModel::GetCorrection(CorInit settings, TLorentzVector vec, 
					 float& emf, float etad) {
  int nrun = settings.Nrun; 
  int nVertex = settings.nvx; 
  int coneSize=settings.cone;  
  int version=settings.version;  
  int syscode=settings.sys; 
  int level=settings.level;
  int imode=settings.imode;

  JetEnergyCorrections myJetEnergyCorrections=
    JetEnergyCorrections("JetCorrections","JetCorrections",
			 level,nVertex,coneSize,version,syscode,nrun,imode); 

//   HepLorentzVector P4Jet; 
//   P4Jet.setPx(vec.Px()); 
//   P4Jet.setPy(vec.Py()); 
//   P4Jet.setPz(vec.Pz()); 
//   P4Jet.setE(vec.E());
  if(abs(syscode)>0) 
    {
      int n_sigma= (syscode>0) ? 1 : -1 ;
      myJetEnergyCorrections.setTotalSysUncertainties(n_sigma);
    }
  double scaleFactor = myJetEnergyCorrections.doEnergyCorrections(vec.Pt(),emf,etad); 
  return scaleFactor; 
}

double TStnMetModel::GetCorrection(TLorentzVector vec, float& emf, float etad) {

  int nrun = myJetRun; 
  int nVertex = myNvx_class12; 
  int coneSize=fJTC_coneSize;  
  int version=fJTC_version;  
  int syscode=fJTC_systcode; 
  int level=fJTC_level;
  int imode=fJTC_imode;

  JetEnergyCorrections myJetEnergyCorrections=
    JetEnergyCorrections("JetCorrections","JetCorrections",
			 level,nVertex,coneSize,version,syscode,nrun,imode); 

//   HepLorentzVector P4Jet; 
//   P4Jet.setPx(vec.Px()); 
//   P4Jet.setPy(vec.Py()); 
//   P4Jet.setPz(vec.Pz()); 
//   P4Jet.setE(vec.E()); 
  if(abs(syscode)>0) 
    {
      int n_sigma= (syscode>0) ? 1 : -1 ;
      myJetEnergyCorrections.setTotalSysUncertainties(n_sigma);
    }
  double scaleFactor = myJetEnergyCorrections.doEnergyCorrections(vec.Pt(),emf,etad); 
  return scaleFactor; 
}


//________ creates a list of towers for all jets
void TStnMetModel::MatchCalorTowers(int jet_ind,
				 TStnJetBlock* fJetBlock,
				 TCalDataBlock *fCalDataBlock,
				 CalDataArray* cdholder) {
  cdholder->clear();
  TStnLinkBlock* links = fJetBlock->TowerLinkList();
  int nptow = links->NLinks(jet_ind);

  TCalTower* ctower = NULL;
  for(int j=0; j<nptow; j++) 
    {
      int iptow = links->Index(jet_ind,j);
      int iphi = iptow & 0x3F;
      int ieta = (iptow>>8) & 0x3F;
      
      ctower = fCalDataBlock->Tower(ieta,iphi);
      cdholder->push_back(ctower);
    }

  std::sort(cdholder->begin(), cdholder->end(),SortTowersByEnergy);
  
  return;
}

//________ creates a list of towers for all EM objects
void TStnMetModel::MatchCalorTowers(int pho_ind,
				 TStnPhotonBlock* fPhotonBlock,
				 TCalDataBlock *fCalDataBlock,
				 CalDataArray* cdholder) {
  cdholder->clear();
  TStnLinkBlock* links = fPhotonBlock->TowerLinkList();
  int nptow = links->NLinks(pho_ind);
  TCalTower* ctower = NULL;
  for(int j=0; j<nptow; j++) 
    {
      int iptow = links->Index(pho_ind,j);
      int iphi = iptow & 0x3F;
      int ieta = (iptow>>8) & 0x3F;
      
      ctower = fCalDataBlock->Tower(ieta,iphi);
      cdholder->push_back(ctower);
    }
  
  std::sort(cdholder->begin(),cdholder->end(),SortTowersByEnergy);
  
  return;
}

//________ creates a list of towers for electrons
void TStnMetModel::MatchCalorTowers(int ele_ind,
				 TStnElectronBlock* fElectronBlock,
				 TCalDataBlock *fCalDataBlock,
				 CalDataArray* cdholder) {
  cdholder->clear();
  TStnLinkBlock* links = fElectronBlock->TowerLinkList();
  int nptow = links->NLinks(ele_ind);
  TCalTower* ctower = NULL;
  for(int j=0; j<nptow; j++) 
    {
      int iptow = links->Index(ele_ind,j);
      int iphi = iptow & 0x3F;
      int ieta = (iptow>>8) & 0x3F;
      
      ctower = fCalDataBlock->Tower(ieta,iphi);
      cdholder->push_back(ctower);
    }
  
  std::sort(cdholder->begin(),cdholder->end(),SortTowersByEnergy);
  
  return;
}

//__________________________________________________ returns 1 if jet is matched  
int TStnMetModel::MyMatchPhoJet(int jet_ind, int pho_ind, int ele_ind,TStnJetBlock* fJetBlock,
			     TLorentzVector *mypho,TLorentzVector *myjet)
{
  int match_code=0;
  
  if(jet_ind<0) return match_code;
  if(pho_ind<0 && ele_ind<0) return match_code;
  
  if(pho_ind>=0)
    {
      CalDataArray jet_towers;
      CalDataArray pho_towers;
      MatchCalorTowers(jet_ind,fJetBlock,fCalData,&jet_towers);
      MatchCalorTowers(pho_ind,fPhotonBlock,fCalData,&pho_towers);
      CalDataArrayI tt_em;
      CalDataArrayI tt_jet;
      for(tt_em = pho_towers.begin(); tt_em != pho_towers.end(); tt_em++)
	{
	  TCalTower* calt_em = *tt_em;
	  int iEta_em = calt_em->IEta();
	  int iPhi_em = calt_em->IPhi();
	  for(tt_jet = jet_towers.begin(); tt_jet != jet_towers.end(); tt_jet++)
	    {
	      TCalTower* calt_jet = *tt_jet;
	      int iEta_jet = calt_jet->IEta();
	      int iPhi_jet = calt_jet->IPhi();
	      if(iEta_em==iEta_jet && iPhi_em==iPhi_jet) match_code=1;
	      if(match_code==1) break;
	    } 
	  if(match_code==1) break; 
	}  
    }
  
  if(ele_ind>=0)
    {
      CalDataArray jet_towers;
      CalDataArray ele_towers;
      MatchCalorTowers(jet_ind,fJetBlock,fCalData,&jet_towers);
      MatchCalorTowers(ele_ind,fElectronBlock,fCalData,&ele_towers);
      CalDataArrayI tt_em;
      CalDataArrayI tt_jet;
      for(tt_em = ele_towers.begin(); tt_em != ele_towers.end(); tt_em++)
	{
	  TCalTower* calt_em = *tt_em;
	  int iEta_em = calt_em->IEta();
	  int iPhi_em = calt_em->IPhi();
	  for(tt_jet = jet_towers.begin(); tt_jet != jet_towers.end(); tt_jet++)
	    {
	      TCalTower* calt_jet = *tt_jet;
	      int iEta_jet = calt_jet->IEta();
	      int iPhi_jet = calt_jet->IPhi();
	      if(iEta_em==iEta_jet && iPhi_em==iPhi_jet) match_code=1;
	      if(match_code==1) break;
	    }
	  if(match_code==1) break;
	}  
    }
  return match_code;
}

//__________________ prints photon towers
void TStnMetModel::PrintPhoTowers(int pho_ind,TStnPhotonBlock* fPhotonBlock,TCalDataBlock *fCalDataBlock) {
  if(pho_ind>=0)
    {
      CalDataArray pho_towers;
      MatchCalorTowers(pho_ind,fPhotonBlock,fCalData,&pho_towers);
      CalDataArrayI tt_em;
      for(tt_em = pho_towers.begin(); tt_em != pho_towers.end(); tt_em++)
	{
	  TCalTower* calt_em = *tt_em;
	  int iEta_em = calt_em->IEta();
	  int iPhi_em = calt_em->IPhi();
	  std::cout<<"---- Photon tower: iEta="<<iEta_em<<" IPhi="<<iPhi_em<<" Et="<<calt_em->Et()<<std::endl;
	}  
    }
  return;  
}

//__________________ prints electron towers
void TStnMetModel::PrintEleTowers(int ele_ind,TStnElectronBlock* fElectronBlock,TCalDataBlock *fCalDataBlock) {
  if(ele_ind>=0)
    {
      CalDataArray ele_towers;
      MatchCalorTowers(ele_ind,fElectronBlock,fCalData,&ele_towers);
      CalDataArrayI tt_em;
      for(tt_em = ele_towers.begin(); tt_em != ele_towers.end(); tt_em++)
	{
	  TCalTower* calt_em = *tt_em;
	  int iEta_em = calt_em->IEta();
	  int iPhi_em = calt_em->IPhi();
	  std::cout<<"---- Electron tower: iEta="<<iEta_em<<" IPhi="<<iPhi_em<<" Et="<<calt_em->Et()<<std::endl;
	}  
    }
  return;
}

//___________________ prints jet towers
void TStnMetModel::PrintJetTowers(int jet_ind,TStnJetBlock* fJetBlock,TCalDataBlock *fCalDataBlock) {
  if(jet_ind>=0)
    {
      CalDataArray jet_towers;
      MatchCalorTowers(jet_ind,fJetBlock,fCalData,&jet_towers);
      CalDataArrayI tt_jet;
      for(tt_jet = jet_towers.begin(); tt_jet != jet_towers.end(); tt_jet++)
	{
	  TCalTower* calt_jet = *tt_jet;
	  int iEta_jet = calt_jet->IEta();
	  int iPhi_jet = calt_jet->IPhi();
	  std::cout<<"---- Jet tower: iEta="<<iEta_jet<<" IPhi="<<iPhi_jet<<" Et="<<calt_jet->Et()<<std::endl;
	}
    }
  return;
}

//___________ corrects Met & SumEt for jets, EM objetcs, and muons/tracks
void TStnMetModel::DoMyMet(CommonStuff miscstuff, JetStuff &jetstuff) {
  //-----correcting SumET
  jetstuff.mySumEtCorr_th15=miscstuff.mySumEt_raw; // make it the same for a moment
  jetstuff.mySumEtJet_th15=0.0; 
  //
  TLorentzVector MetPho(0.0,0.0,0.0,0.0);
  TLorentzVector MetEle(0.0,0.0,0.0,0.0);
  //_________ correcting all SumEt & MET for photons 
  for(int i=0; i<miscstuff.myRawPhoton.size(); i++)
    {
      jetstuff.mySumEtCorr_th15=jetstuff.mySumEtCorr_th15-miscstuff.myRawPhoton[i].Pt(); 
      MetPho=MetPho+miscstuff.myCorrPhoton[i]-miscstuff.myRawPhoton[i];
    }
  //_________ correcting all SumEt & MET for electrons 
  for(int i=0; i<miscstuff.myRawElectron.size(); i++)
    {
      jetstuff.mySumEtCorr_th15=jetstuff.mySumEtCorr_th15-miscstuff.myRawElectron[i].Pt(); 
      MetEle=MetEle+miscstuff.myCorrElectron[i]-miscstuff.myRawElectron[i];
    }
  if(MetPho.E()<MetPho.P() || MetPho.E()<0.0) MetPho.SetPxPyPzE(0.0,0.0,0.0,0.0);
  if(MetEle.E()<MetEle.P() || MetEle.E()<0.0) MetEle.SetPxPyPzE(0.0,0.0,0.0,0.0);
  //_________ correcting all SumEt and calculating Met correction for jets
  TLorentzVector MetJet15(0.0,0.0,0.0,0.0);
  for(int i=0; i<jetstuff.Jet_raw_noEMobj.size(); i++)
    {
      //--------- may need an eta cut here
      if(fabs(jetstuff.Jet_lev6_noEMobj[i].Eta())>=MinJetEta && fabs(jetstuff.Jet_lev6_noEMobj[i].Eta())<=MaxJetEta)
	{
	  if(jetstuff.Jet_lev6_noEMobj[i].Pt()>15.0) 
	    {
	      // corr=lev5-raw+M.I.=lev5-raw+lev1-lev4, 
	      jetstuff.mySumEtCorr_th15=jetstuff.mySumEtCorr_th15-jetstuff.Jet_raw_noEMobj[i].Pt()
		+jetstuff.Jet_lev5_noEMobj[i].Pt()
		+jetstuff.Jet_lev1_noEMobj[i].Pt()
		-jetstuff.Jet_lev4_noEMobj[i].Pt()
		-jetstuff.Jet_lev6_noEMobj[i].Pt();
	      jetstuff.mySumEtJet_th15=jetstuff.mySumEtJet_th15+jetstuff.Jet_lev5_noEMobj[i].Pt()
		+jetstuff.Jet_lev1_noEMobj[i].Pt()
		-jetstuff.Jet_lev4_noEMobj[i].Pt();
	      MetJet15=MetJet15+jetstuff.Jet_lev5_noEMobj[i]
		+jetstuff.Jet_lev1_noEMobj[i]
		-jetstuff.Jet_lev4_noEMobj[i]
		-jetstuff.Jet_raw_noEMobj[i];
	    } 
	}
    }

  //_________ correcting all MET for muons/tracks 
  TLorentzVector MetMuon(0.0,0.0,0.0,0.0);
  for(int i=0; i<miscstuff.myCorrMuon.size(); i++)
    {
      MetMuon=MetMuon+miscstuff.myCorrMuon[i]-miscstuff.myCalMuon[i];
    }
  if(MetMuon.E()<MetMuon.P() || MetMuon.E()<0.0) MetEle.SetPxPyPzE(0.0,0.0,0.0,0.0);

  //---- making sure all MetJet makes sense:
  if(MetJet15.E()<MetJet15.P()) MetJet15.SetPxPyPzE(0.0,0.0,0.0,0.0); 
  //---- this part was modified on 09/20/06
  MetJet15=MetJet15+MetEle+MetPho+MetMuon; // modified on 09/20/06 
  jetstuff.myMETcorr_th15.Set(MetJet15.Px(),MetJet15.Py());
  //---- finally, calculating corrected MET
  jetstuff.myMETcorr_th15=miscstuff.myMET_raw-jetstuff.myMETcorr_th15;
  
  return;
}

//_____________ main routine that does metsig for data & MC
void TStnMetModel::DoMyMetSignificanceData(CommonStuff miscstuff, JetStuff jetstuff, MetsigStuff &sigstuff) {

  if(fJTC_imode==1) // data
    {
      sigstuff.metsig_Def=MyTotalMetSignificance(jetstuff,jetstuff.myMETcorr_th15,0,0);
      for(int i=0; i<MetsigArray; i++)
	{
	  sigstuff.metsig_syst[i]=sigstuff.metsig_Def; // no systematics for data
	}
    }
  if(fJTC_imode==0) // MC
    {
      sigstuff.metsig_Def=MyTotalMetSignificance(jetstuff,jetstuff.myMETcorr_th15,0,0);
      double def=sigstuff.metsig_Def;
      //-----------------------------------------------------------------------------------------------------
      //-- syst: different U.E. parameterization
      int jer_code=0;
      int ue_code=1;
      sigstuff.metsig_syst[0]=MyTotalMetSignificance(jetstuff,jetstuff.myMETcorr_th15,jer_code,ue_code);
      //-----------------------------------------------------------------------------------------------------
      //-- syst: Uncl. En. systematics: mean-G 
      jer_code=0;
      ue_code=-2;
      double syst_dn=MyTotalMetSignificance(jetstuff,jetstuff.myMETcorr_th15,jer_code,ue_code);
      //-- syst: Uncl. En. systematics: mean+G 
      jer_code=0;
      ue_code=2;
      double syst_up=MyTotalMetSignificance(jetstuff,jetstuff.myMETcorr_th15,jer_code,ue_code);
      sigstuff.metsig_syst[1]= fabs(syst_up-def) > fabs(syst_dn-def) ? syst_up : syst_dn;
      //-----------------------------------------------------------------------------------------------------
      //-- syst: Uncl. En. systematics: sigma-G
      jer_code=0;
      ue_code=-3;
      syst_dn=MyTotalMetSignificance(jetstuff,jetstuff.myMETcorr_th15,jer_code,ue_code);
      //-- syst: Uncl. En. systematics: sigma+G
      jer_code=0;
      ue_code=3;
      syst_up=MyTotalMetSignificance(jetstuff,jetstuff.myMETcorr_th15,jer_code,ue_code);
      sigstuff.metsig_syst[2]= fabs(syst_up-def) > fabs(syst_dn-def) ? syst_up : syst_dn;
      //-----------------------------------------------------------------------------------------------------
      //-- syst: Uncl. En. systematics: scale-G
      jer_code=0;
      ue_code=-4;
      syst_dn=MyTotalMetSignificance(jetstuff,jetstuff.myMETcorr_th15,jer_code,ue_code);      
      //-- syst: Uncl. En. systematics: scale+G
      jer_code=0;
      ue_code=4;
      syst_up=MyTotalMetSignificance(jetstuff,jetstuff.myMETcorr_th15,jer_code,ue_code);
      sigstuff.metsig_syst[3]= fabs(syst_up-def) > fabs(syst_dn-def) ? syst_up : syst_dn;

      //-----------------------------------------------------------------------------------------------------
      //-- syst: Uncl. En. systematics: norm-G
      jer_code=0;
      ue_code=-5;
      syst_dn=MyTotalMetSignificance(jetstuff,jetstuff.myMETcorr_th15,jer_code,ue_code);      
      //-- syst: Uncl. En. systematics: norm+G
      jer_code=0;
      ue_code=5;
      syst_up=MyTotalMetSignificance(jetstuff,jetstuff.myMETcorr_th15,jer_code,ue_code);
      sigstuff.metsig_syst[4]= fabs(syst_up-def) > fabs(syst_dn-def) ? syst_up : syst_dn;
      //-----------------------------------------------------------------------------------------------------
      //-- syst: JER systematics due to Data-MC difference: sigma-G
      jer_code=-6;
      ue_code=0;
      syst_dn=MyTotalMetSignificance(jetstuff,jetstuff.myMETcorr_th15,jer_code,ue_code);      
      //-- syst: JER systematics due to Data-MC difference: sigma+G
      jer_code=6;
      ue_code=0;
      syst_up=MyTotalMetSignificance(jetstuff,jetstuff.myMETcorr_th15,jer_code,ue_code);
      sigstuff.metsig_syst[5]= fabs(syst_up-def) > fabs(syst_dn-def) ? syst_up : syst_dn;

      //------------------ remaining syst
      for(int i=6; i<MetsigArray; i++)
	{
	  sigstuff.metsig_syst[i]=sigstuff.metsig_Def; // no systematics for data
	}
    }
  return;
}

//______ Main routine that does metsig stuff based on pseudo-experiments
void TStnMetModel::DoMyMetSignificancePseudo(CommonStuff miscstuff, JetStuff jetstuff, MetsigStuff &sigstuff) {
  for(int i=0; i<Npseudo; i++)
    {
      //______________________ Met Model default prediction
      int rnd_seed=gRandom->GetSeed();
      TVector2 myGenMet(0.0,0.0);
      int jer_code=0; // JER systematics code
      int ue_code=0; // unclustered energy parameterization systematics 
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet);
      sigstuff.MetPseudo.push_back(myGenMet);
      double def=MyTotalMetSignificance(jetstuff,myGenMet,jer_code,ue_code);
      sigstuff.metsig_DefPseudo.push_back(def);
      //______________________ Met Model systematics: z->ee vs. dipho sideband parameterization
      myGenMet.Set(0.0,0.0);
      jer_code=0;
      ue_code=1;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet);
      sigstuff.MetPseudoSyst[0].push_back(myGenMet);
      sigstuff.metsig_systPseudo[0].push_back(MyTotalMetSignificance(jetstuff,myGenMet,jer_code,ue_code));
      //______________________ Met Model Uncl. En. systematics: mean
      TVector2 myGenMet_up(0.0,0.0);
      TVector2 myGenMet_dn(0.0,0.0);
      jer_code=0;
      ue_code=-2;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_dn);
      double syst_dn=MyTotalMetSignificance(jetstuff,myGenMet_dn,jer_code,ue_code);
      jer_code=0;
      ue_code=2;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_up);
      double syst_up=MyTotalMetSignificance(jetstuff,myGenMet_up,jer_code,ue_code);
      if(fabs(syst_up-def)>fabs(syst_dn-def))
	{
	  sigstuff.MetPseudoSyst[1].push_back(myGenMet_up);
	  sigstuff.metsig_systPseudo[1].push_back(syst_up);
	}
      else
	{
	  sigstuff.MetPseudoSyst[1].push_back(myGenMet_dn);
	  sigstuff.metsig_systPseudo[1].push_back(syst_dn);
	}
      //______________________ Met Model Uncl. En. systematics: sigma
      myGenMet_up.Set(0.0,0.0);
      myGenMet_dn.Set(0.0,0.0);
      jer_code=0;
      ue_code=-3;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_dn);
      syst_dn=MyTotalMetSignificance(jetstuff,myGenMet_dn,jer_code,ue_code);
      jer_code=0;
      ue_code=3;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_up);
      syst_up=MyTotalMetSignificance(jetstuff,myGenMet_up,jer_code,ue_code);
      if(fabs(syst_up-def)>fabs(syst_dn-def))
	{
	  sigstuff.MetPseudoSyst[2].push_back(myGenMet_up);
	  sigstuff.metsig_systPseudo[2].push_back(syst_up);
	}
      else
	{
	  sigstuff.MetPseudoSyst[2].push_back(myGenMet_dn);
	  sigstuff.metsig_systPseudo[2].push_back(syst_dn);
	}
      //______________________ Met Model Uncl. En. systematics: scale
      myGenMet_up.Set(0.0,0.0);
      myGenMet_dn.Set(0.0,0.0);
      jer_code=0;
      ue_code=-4;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_dn);
      syst_dn=MyTotalMetSignificance(jetstuff,myGenMet_dn,jer_code,ue_code);
      jer_code=0;
      ue_code=4;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_up);
      syst_up=MyTotalMetSignificance(jetstuff,myGenMet_up,jer_code,ue_code);
      if(fabs(syst_up-def)>fabs(syst_dn-def))
	{
	  sigstuff.MetPseudoSyst[3].push_back(myGenMet_up);
	  sigstuff.metsig_systPseudo[3].push_back(syst_up);
	}
      else
	{
	  sigstuff.MetPseudoSyst[3].push_back(myGenMet_dn);
	  sigstuff.metsig_systPseudo[3].push_back(syst_dn);
	}
      //______________________ Met Model Uncl. En. systematics: norm
      myGenMet_up.Set(0.0,0.0);
      myGenMet_dn.Set(0.0,0.0);
      jer_code=0;
      ue_code=-5;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_dn);
      syst_dn=MyTotalMetSignificance(jetstuff,myGenMet_dn,jer_code,ue_code);
      jer_code=0;
      ue_code=5;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_up);
      syst_up=MyTotalMetSignificance(jetstuff,myGenMet_up,jer_code,ue_code);
      if(fabs(syst_up-def)>fabs(syst_dn-def))
	{
	  sigstuff.MetPseudoSyst[4].push_back(myGenMet_up);
	  sigstuff.metsig_systPseudo[4].push_back(syst_up);
	}
      else
	{
	  sigstuff.MetPseudoSyst[4].push_back(myGenMet_dn);
	  sigstuff.metsig_systPseudo[4].push_back(syst_dn);
	}
      //-----------------------------------------------------------------------------------------------------
      //-- syst: JER systematics due to Data-MC difference
      myGenMet_up.Set(0.0,0.0);
      myGenMet_dn.Set(0.0,0.0);
      jer_code=-6;
      ue_code=0;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_dn);
      syst_dn=MyTotalMetSignificance(jetstuff,myGenMet_dn,jer_code,ue_code);
      jer_code=6;
      ue_code=0;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_up);
      syst_up=MyTotalMetSignificance(jetstuff,myGenMet_up,jer_code,ue_code);
      if(fabs(syst_up-def)>fabs(syst_dn-def))
	{
	  sigstuff.MetPseudoSyst[5].push_back(myGenMet_up);
	  sigstuff.metsig_systPseudo[5].push_back(syst_up);
	}
      else
	{
	  sigstuff.MetPseudoSyst[5].push_back(myGenMet_dn);
	  sigstuff.metsig_systPseudo[5].push_back(syst_dn);
	}
      //______________________ Met Model JER systematics: meanG
      myGenMet_up.Set(0.0,0.0);
      myGenMet_dn.Set(0.0,0.0);
      jer_code=-1;
      ue_code=0;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_dn);
      syst_dn=MyTotalMetSignificance(jetstuff,myGenMet_dn,jer_code,ue_code);
      jer_code=1;
      ue_code=0;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_up);
      syst_up=MyTotalMetSignificance(jetstuff,myGenMet_up,jer_code,ue_code);
      if(fabs(syst_up-def)>fabs(syst_dn-def))
	{
	  sigstuff.MetPseudoSyst[6].push_back(myGenMet_up);
	  sigstuff.metsig_systPseudo[6].push_back(syst_up);
	}
      else
	{
	  sigstuff.MetPseudoSyst[6].push_back(myGenMet_dn);
	  sigstuff.metsig_systPseudo[6].push_back(syst_dn);
	}
      //______________________ Met Model JER systematics: sigmaG
      myGenMet_up.Set(0.0,0.0);
      myGenMet_dn.Set(0.0,0.0);
      jer_code=-2;
      ue_code=0;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_dn);
      syst_dn=MyTotalMetSignificance(jetstuff,myGenMet_dn,jer_code,ue_code);
      jer_code=2;
      ue_code=0;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_up);
      syst_up=MyTotalMetSignificance(jetstuff,myGenMet_up,jer_code,ue_code);
      if(fabs(syst_up-def)>fabs(syst_dn-def))
	{
	  sigstuff.MetPseudoSyst[7].push_back(myGenMet_up);
	  sigstuff.metsig_systPseudo[7].push_back(syst_up);
	}
      else
	{
	  sigstuff.MetPseudoSyst[7].push_back(myGenMet_dn);
	  sigstuff.metsig_systPseudo[7].push_back(syst_dn);
	}
      //______________________ Met Model JER systematics: mpvL
      myGenMet_up.Set(0.0,0.0);
      myGenMet_dn.Set(0.0,0.0);
      jer_code=-3;
      ue_code=0;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_dn);
      syst_dn=MyTotalMetSignificance(jetstuff,myGenMet_dn,jer_code,ue_code);
      jer_code=3;
      ue_code=0;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_up);
      syst_up=MyTotalMetSignificance(jetstuff,myGenMet_up,jer_code,ue_code);
      if(fabs(syst_up-def)>fabs(syst_dn-def))
	{
	  sigstuff.MetPseudoSyst[8].push_back(myGenMet_up);
	  sigstuff.metsig_systPseudo[8].push_back(syst_up);
	}
      else
	{
	  sigstuff.MetPseudoSyst[8].push_back(myGenMet_dn);
	  sigstuff.metsig_systPseudo[8].push_back(syst_dn);
	}
      //______________________ Met Model JER systematics: sigmaL
      myGenMet_up.Set(0.0,0.0);
      myGenMet_dn.Set(0.0,0.0);
      jer_code=-4;
      ue_code=0;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_dn);
      syst_dn=MyTotalMetSignificance(jetstuff,myGenMet_dn,jer_code,ue_code);
      jer_code=4;
      ue_code=0;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_up);
      syst_up=MyTotalMetSignificance(jetstuff,myGenMet_up,jer_code,ue_code);
      if(fabs(syst_up-def)>fabs(syst_dn-def))
	{
	  sigstuff.MetPseudoSyst[9].push_back(myGenMet_up);
	  sigstuff.metsig_systPseudo[9].push_back(syst_up);
	}
      else
	{
	  sigstuff.MetPseudoSyst[9].push_back(myGenMet_dn);
	  sigstuff.metsig_systPseudo[9].push_back(syst_dn);
	}
      //______________________ Met Model JER systematics: norm
      myGenMet_up.Set(0.0,0.0);
      myGenMet_dn.Set(0.0,0.0);
      jer_code=-5;
      ue_code=0;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_dn);
      syst_dn=MyTotalMetSignificance(jetstuff,myGenMet_dn,jer_code,ue_code);
      jer_code=5;
      ue_code=0;
      GenerateMyTotalMet(jetstuff,miscstuff,jer_code,ue_code,rnd_seed,myGenMet_up);
      syst_up=MyTotalMetSignificance(jetstuff,myGenMet_up,jer_code,ue_code);
      if(fabs(syst_up-def)>fabs(syst_dn-def))
	{
	  sigstuff.MetPseudoSyst[10].push_back(myGenMet_up);
	  sigstuff.metsig_systPseudo[10].push_back(syst_up);
	}
      else
	{
	  sigstuff.MetPseudoSyst[10].push_back(myGenMet_dn);
	  sigstuff.metsig_systPseudo[10].push_back(syst_dn);
	}

      //__________ remaining syst. to be filled with default values
      for(int j=11; j<MetsigArray; j++)
	{
	  sigstuff.MetPseudoSyst[j].push_back(sigstuff.MetPseudo[i]);
	  sigstuff.metsig_systPseudo[j].push_back(sigstuff.metsig_DefPseudo[i]);
	}
    }
  return;
}



//___ This function returns my Met significance for all events
double TStnMetModel::MyTotalMetSignificance(JetStuff jetstuff, TVector2 myMetVec, int systcode, int systcodeUncl) {

  double max_MetProbInteg=1.0;
  InitMetProbStuff(max_MetProbInteg,jetstuff,myMetVec,systcode,systcodeUncl); //inits params for MetProb

  double metSig_est=0.0;
  if(max_MetProbInteg>0.0 && max_MetProbInteg<1.0) metSig_est=-1.0*log10(1.0-max_MetProbInteg);
  if(max_MetProbInteg<=0.0) metSig_est=0.0;
  if(max_MetProbInteg>=1.0) metSig_est=19.0;
  metSig_est=MetSigCorrection(metSig_est);
  return metSig_est;
}

//------- New version as of 06/23/08. Takes into account MET direction (to some extent) by projecting jets on MET
//___ This function inits params for MetProb
void TStnMetModel::InitMetProbStuff(double &metprobInt,JetStuff jetstuff, TVector2 myMetVec, int systcode, int systcodeUncl) {

  double jet_EtCut=15.0;
  double _met=myMetVec.Mod();
  double _metX=fabs(myMetVec.X());
  double _metY=fabs(myMetVec.Y());

  //-------------------------- Getting unclustered energy parameters
  double _meanX=0.0;
  double _sigmaX=0.0;
  double _meanY=0.0;
  double _sigmaY=0.0;
  double _normX=0.0;
  double _scaleX=0.0;
  double _normY=0.0;
  double _scaleY=0.0;
  GetMyUnclMetResolution(jetstuff.mySumEtCorr_th15,fJTC_imode,fUnclParamSwitch,systcodeUncl,_sigmaX,_sigmaY,_meanX,_meanY,_normX,_normY,_scaleX,_scaleY);

  //------------------------------------------------------------
  // calculating precision for multi-D integration
  //______ uncl. energy part
  TF1* ueFun=new TF1("ue",UnclMetResolution,-1000.0,1000.0,4);
  ueFun->SetParameter(0,_meanX);
  ueFun->SetParameter(1,_sigmaX);
  ueFun->SetParameter(2,_normX);
  ueFun->SetParameter(3,_scaleX);
  double ue_lim=5.0*_scaleX*_sigmaX>_met ? 5.0*_scaleX*_sigmaX : 2.0*_met;
  double scaleF=ueFun->Integral(-1.0*ue_lim,1.0*ue_lim);
  double test_prob;
  if(scaleF>0.0) test_prob=1.0-(ueFun->Integral(-1.0*_met,_met))/scaleF;
  else test_prob=1.0;
  if(test_prob<=0.0) test_prob=0.0;
  if(test_prob>1.0) test_prob=1.0; 
  metprobInt=metprobInt*(1.0-test_prob); // calculating upper limit on MetProb
  ueFun->SetParameter(0,_meanY);
  ueFun->SetParameter(1,_sigmaY);
  ueFun->SetParameter(2,_normY);
  ueFun->SetParameter(3,_scaleY);
  ue_lim=5.0*_scaleY*_sigmaY>_met ? 5.0*_scaleY*_sigmaY : 2.0*_met;
  scaleF=ueFun->Integral(-1.0*ue_lim,1.0*ue_lim);
  if(scaleF>0.0) test_prob=1.0-(ueFun->Integral(-1.0*_met,_met))/scaleF;
  else test_prob=1.0;
  if(test_prob<=0.0) test_prob=0.0;
  if(test_prob>1.0) test_prob=1.0; 
  metprobInt=metprobInt*(1.0-test_prob); // calculating upper limit on MetProb
  delete ueFun;

  //------------------------------------------------------------
  //______ jet part
  for(int i=0; i<jetstuff.Jet_lev6_noEMobj.size(); i++)
    {
      if(jetstuff.Jet_lev6_noEMobj[i].Pt()>3.0 && jetstuff.Jet_raw_noEMobj[i].Pt()>3.0 && fabs(jetstuff.EtaDet[i])<=MaxJetEta)
	{
	  double dPhi_jetmet=fabs(TVector2::Phi_mpi_pi(jetstuff.Jet_lev6_noEMobj[i].Phi()-myMetVec.Phi()));
	  double trueEjet=jetstuff.Jet_lev6_noEMobj[i].E();
	  double scaleE=jetstuff.Jet_lev6_noEMobj[i].E()/jetstuff.Jet_lev6_noEMobj[i].Pt();
	  double cos_dphi=cos(dPhi_jetmet);
	  if(fabs(cos_dphi)>0.0) // only consider these cases (no contribution is dPhi_jetmet=pi/2)
	    {
	      trueEjet=trueEjet+scaleE*myMetVec.Mod()/cos_dphi;
	      if(trueEjet>3.0 && trueEjet<1000.0)
		{
		  double parJER[5];
		  parJER[0]=MyJer_meanG(trueEjet,jetstuff.EtaDetCorr[i],0,systcode); // for now systcode is used as jer_stat_code
		  parJER[1]=MyJer_sigmaG(trueEjet,jetstuff.EtaDetCorr[i],0,systcode);
		  parJER[2]=MyJer_mpvL(trueEjet,jetstuff.EtaDetCorr[i],0,systcode);
		  parJER[3]=MyJer_sigmaL(trueEjet,jetstuff.EtaDetCorr[i],0,systcode);
		  parJER[4]=MyJer_normG(trueEjet,jetstuff.EtaDetCorr[i],0,systcode);

		  double jer_limit=MyIntegralUpLimit(trueEjet,jetstuff.EtaDetCorr[i],0,systcode);
		  TF1* JerFun=new TF1("jer",MyJER,-1.0,jer_limit,5);
		  for(int par_ind=0; par_ind<5; par_ind++)
		    {
		      JerFun->SetParameter(par_ind,parJER[par_ind]);
		    }
		  double off_set=JerFun->GetMaximumX(-1.0,1.0); // returns mpv of JER function;		  
		  if(trueEjet>jetstuff.Jet_lev6_noEMobj[i].E()) // case-1: MET points along the jet
		    {
		      double x1=jetstuff.Jet_lev6_noEMobj[i].E()/trueEjet-1.0+off_set;
		      if(x1<-1.0) x1=-1.0;
		      if(x1>jer_limit) x1=jer_limit;
		      test_prob=JerFun->Integral(-1.0,x1);
		      if(test_prob<=0.0) test_prob=0.0;
		      double tot_jetInt=test_prob+(JerFun->Integral(x1,jer_limit));
		      if(tot_jetInt>0.0)
			{
			  test_prob=test_prob/tot_jetInt;
			  double prob_factor=1.0-test_prob;
			  if(prob_factor>1.0) prob_factor=1.0;
			  if(prob_factor<=0.0) prob_factor=1.0E-50;
			  metprobInt=metprobInt*prob_factor; // calculating upper limit on MetProb
			}
		    }
		  else // case-2: MET points opposite to the jet
		    {
		      double x1=jetstuff.Jet_lev6_noEMobj[i].E()/trueEjet-1.0+off_set;
		      if(x1<-1.0) x1=-1.0;
		      if(x1>jer_limit) x1=jer_limit;
		      test_prob=JerFun->Integral(x1,jer_limit);
		      if(test_prob<=0.0) test_prob=0.0;
		      double tot_jetInt=test_prob+(JerFun->Integral(-1.0,x1));
		      if(tot_jetInt>0.0) 
			{
			  test_prob=test_prob/tot_jetInt;
			  double prob_factor=1.0-test_prob;
			  if(prob_factor>1.0) prob_factor=1.0;
			  if(prob_factor<=0.0) prob_factor=1.0E-50;
			  metprobInt=metprobInt*prob_factor; // calculating upper limit on MetProb
			}
		    }
		  delete JerFun;
		}
	    }
	}
    }
  return;
}

//______________ This functio returns corrected (true) MetSig
double TStnMetModel::MetSigCorrection(double rawSig) {
  double corr=1.0;
  if(fSampleID>=0 && fSampleID<=4)
    {
      float cut1=0.0;
      float cut2=rawSig<20.0 ? 20.0 : rawSig+1.0;
      TF1 *fitFun=new TF1("fit",rawMetSigFunc,cut1,cut2,5);
      int njet=jet04stuff.myNjet_th15;
      fitFun->SetParameter(0,MetSigCorrPAR(0,njet,fJTC_imode,fSampleID));
      fitFun->SetParameter(1,MetSigCorrPAR(1,njet,fJTC_imode,fSampleID));
      fitFun->SetParameter(2,MetSigCorrPAR(2,njet,fJTC_imode,fSampleID));
      fitFun->SetParameter(3,MetSigCorrPAR(3,njet,fJTC_imode,fSampleID));
      fitFun->SetParameter(4,MetSigCorrPAR(4,njet,fJTC_imode,fSampleID));
      
      double bin_int=0.0;
      double bin_int0=0.0;
      double last_bin_slope=1.0;
      double tot_int=fitFun->Integral(0.0,cut2);
      double par5=MetSigCorrPAR(5,njet,fJTC_imode,fSampleID);
      
      bin_int=fitFun->Integral(0.0,rawSig);
      if(tot_int>0.0 && par5>0.0) bin_int=bin_int/tot_int;
      else 
	{
	  delete fitFun;
	  return rawSig;
	}
      if(bin_int>0.0 && bin_int<1.0 && rawSig<par5 && rawSig>0.0) corr=-log10(1-bin_int)/rawSig;
      else 
	{
	  bin_int0=fitFun->Integral(0.0,par5);
	  bin_int0=bin_int0/tot_int;
	  last_bin_slope=-log10(1-bin_int0)/par5;      
	  corr=last_bin_slope;
	}
      if(corr<0.0) corr=1.0;
      delete fitFun;
    }
  return rawSig*corr;
}

//___ This function returns input parameters for MetSig correction
double TStnMetModel::MetSigCorrPAR(int ind, int njet15, int isMC, int sample_id) {
  double param=0.0;
  double par_arr[3][6];
  int jetind=0;
  if(njet15==0) jetind=0;
  if(njet15==1) jetind=1;
  if(njet15>=2) jetind=2;
  if(ind<0 || ind>5) return 0.0;
  if(isMC==0) // for MC
    {
      if(sample_id==0) // Pythia ggX signal
	{
	  // Pythia calibration params are the same as in the e-log entry on 06/23/08
	  par_arr[0][0]=0.008;  par_arr[1][0]=0.1586;  par_arr[2][0]=0.1438;
	  par_arr[0][1]=-3.088; par_arr[1][1]=1.515;   par_arr[2][1]=1.4; 
	  par_arr[0][2]=-2.625;	par_arr[1][2]=-0.104;  par_arr[2][2]=-0.081;
	  par_arr[0][3]=0.1807;	par_arr[1][3]=0.00076; par_arr[2][3]=0.000179;
	  par_arr[0][4]=2.247;  par_arr[1][4]=0.507;   par_arr[2][4]=0.143;  
	  par_arr[0][5]=6.17;   par_arr[1][5]=12.0;   par_arr[2][5]=9.0;  
	}
      if(sample_id==2) // Pythia cem-cem Z->ee
	{
	  // fixed some entires on 06/20/08 (made the same as in 01/20/08 -- "Re-calibration of MetSig")
	  par_arr[0][0]=0.00986; par_arr[1][0]=0.6469;  par_arr[2][0]=0.5971;
	  par_arr[0][1]=-2.19; 	 par_arr[1][1]=1.507;   par_arr[2][1]=1.382; 
	  par_arr[0][2]=-2.10;	 par_arr[1][2]=-0.1073; par_arr[2][2]=-0.0857;
	  par_arr[0][3]=0.1774;	 par_arr[1][3]=0.00564; par_arr[2][3]=0.00121;
	  par_arr[0][4]=2.168;   par_arr[1][4]=0.6112;  par_arr[2][4]=0.359;  
	  par_arr[0][5]=5.97;    par_arr[1][5]=12.0;   par_arr[2][5]=9.0;  
	}
      if(sample_id==3) // Pythia pho-jet(Pt_hat>22) (if Z->ee parameterization is used) ATTN::need to be replaced  
	{
	  par_arr[0][0]=0.00986; par_arr[1][0]=0.1624;    par_arr[2][0]=0.1418;
	  par_arr[0][1]=-2.19; 	 par_arr[1][1]=1.641;     par_arr[2][1]=1.405; 
	  par_arr[0][2]=-2.10;	 par_arr[1][2]=-0.0276;   par_arr[2][2]=-0.06082;
	  par_arr[0][3]=0.1774;	 par_arr[1][3]=0.000395;  par_arr[2][3]=0.000426;
	  par_arr[0][4]=2.168;   par_arr[1][4]=0.4546;    par_arr[2][4]=0.40;  
	  par_arr[0][5]=5.97;    par_arr[1][5]=19.01;     par_arr[2][5]=19.01;  
	}
    }
  if(isMC==1) // for DATA
    {
      if(sample_id==0) // ggX signal
	{
	  // fixed some entires on 06/20/08 (made the same as in 01/20/08 -- "Re-calibration of MetSig")
	  par_arr[0][0]=0.00693; par_arr[1][0]=0.6237; par_arr[2][0]=0.587;
	  par_arr[0][1]=-3.37; 	 par_arr[1][1]=1.387;  par_arr[2][1]=1.316; 
	  par_arr[0][2]=-2.79;	 par_arr[1][2]=-0.180; par_arr[2][2]=-0.127;
	  par_arr[0][3]=0.1838;	 par_arr[1][3]=0.0083; par_arr[2][3]=0.0022;
	  par_arr[0][4]=2.26;    par_arr[1][4]=0.624;  par_arr[2][4]=0.36;  
	  par_arr[0][5]=5.59;    par_arr[1][5]=12.0;  par_arr[2][5]=9.0;  
	}
      if(sample_id==1) // ggX sideband
	{
	  par_arr[0][0]=0.00849; par_arr[1][0]=0.5939; par_arr[2][0]=0.570;
	  par_arr[0][1]=-3.04; 	 par_arr[1][1]=1.295;  par_arr[2][1]=1.273; 
	  par_arr[0][2]=-2.61;	 par_arr[1][2]=-0.227; par_arr[2][2]=-0.141;
	  par_arr[0][3]=0.1813;	 par_arr[1][3]=0.0161; par_arr[2][3]=0.0051;
	  par_arr[0][4]=2.277;   par_arr[1][4]=0.724;  par_arr[2][4]=0.50;  
	  par_arr[0][5]=5.39;    par_arr[1][5]=12.0;   par_arr[2][5]=9.0;  
	}
      if(sample_id==2) // cem-cem Z->ee
	{
	  // fixed some entires on 06/20/08 (made the same as in 01/20/08 -- "Re-calibration of MetSig")
	  par_arr[0][0]=0.0065;  par_arr[1][0]=0.6395; par_arr[2][0]=0.611;
	  par_arr[0][1]=-4.03; 	 par_arr[1][1]=1.504;  par_arr[2][1]=1.466; 
	  par_arr[0][2]=-3.70;	 par_arr[1][2]=-0.140; par_arr[2][2]=-0.041;
	  par_arr[0][3]=0.1794;	 par_arr[1][3]=0.0142; par_arr[2][3]=0.00056;
	  par_arr[0][4]=2.151;   par_arr[1][4]=0.675;  par_arr[2][4]=0.17;  
	  par_arr[0][5]=6.11;    par_arr[1][5]=12.0;  par_arr[2][5]=9.0;  
	}
      if(sample_id==3) // Use Pythia pho-jet(Pt_hat>22) for a moment (if Z->ee parameterization is used) ATTN::need to be replaced  
	{
	  par_arr[0][0]=0.00986; par_arr[1][0]=0.1624;    par_arr[2][0]=0.1418;
	  par_arr[0][1]=-2.19; 	 par_arr[1][1]=1.641;     par_arr[2][1]=1.405; 
	  par_arr[0][2]=-2.10;	 par_arr[1][2]=-0.0276;   par_arr[2][2]=-0.06082;
	  par_arr[0][3]=0.1774;	 par_arr[1][3]=0.000395;  par_arr[2][3]=0.000426;
	  par_arr[0][4]=2.168;   par_arr[1][4]=0.4546;    par_arr[2][4]=0.40;  
	  par_arr[0][5]=5.97;    par_arr[1][5]=19.01;     par_arr[2][5]=19.01;  
	}
    }
  param=par_arr[jetind][ind];
  return param;
}

//___ This function returns upper limit of allowed jet energy fluctuation.
double TStnMetModel::MyIntegralUpLimit(double jet_E, double eta_det, int stat_code, int syst_code) {
  double max_value=10.0; // maximum allowed: E(det)/E(true)-1<=10.0
  double value=10.0; 
  double parJER[5];
  parJER[0]=MyJer_meanG(jet_E,eta_det,0,syst_code); // for now systcode is used as jer_stat_code
  parJER[1]=MyJer_sigmaG(jet_E,eta_det,0,syst_code);
  parJER[2]=MyJer_mpvL(jet_E,eta_det,0,syst_code);
  parJER[3]=MyJer_sigmaL(jet_E,eta_det,0,syst_code);
  parJER[4]=MyJer_normG(jet_E,eta_det,0,syst_code);

  double lim_L=max_value;
  double lim_G=max_value;
  //------ calculating limit for Landau part
  //--- This limit is roughly equivalent to 1.0E-6 significance level of JER Landau part if it is found in range [-1.0;10.0] 
  if(fabs(1+parJER[2]-3.7*parJER[3])>0.0) lim_L=(3.7*parJER[3]-parJER[2])/(1+parJER[2]-3.7*parJER[3]);
  if(lim_L>max_value || lim_L<parJER[3]) lim_L=max_value;
  //------ calculating limit for Gaussian part taking into account relative normalization of Gaussian
  if(parJER[4]>0.0 && parJER[4]<=1.0)
    {
      double sig_tmp=-log10(parJER[4]/(1+parJER[4]));
      double sig=4.2-sig_tmp; // 4.2 roughly corresponds to 4*Sigma significance level (~1.0E-5)
      if(sig<0.0) sig=0.0;
      double n_sigma=0.6807+1.035*sig-0.05574*sig*sig; // crude parametrization of N_sigma vs. significance
      if(fabs(1+parJER[0]-n_sigma*parJER[1])>0.0) lim_G=(n_sigma*parJER[1]-parJER[0])/(1+parJER[0]-n_sigma*parJER[1]);
      else // just to make a "smooth" transition between Landau and Gauss
	{
	  if(parJER[0]<1.0 && parJER[0]>0.0) lim_G=1.0/parJER[0]-1.0;
	  else lim_G=-1.0;
	}
      if(lim_G>max_value || lim_G<parJER[1]) lim_G=max_value;     
    }
  //------ deciding what limit should be used: Landau or Gauss
  if(parJER[4]>0.0) // gauss is non-zero
    {
      if(lim_L<lim_G) value=lim_G; // Landau lim smaller 
      else value=lim_L;// Landau lim larger 
    }
  else value=lim_L; // gauss is zero
  if(value>max_value) value=max_value;
  return value;
}


//___ generates Sigma and Mean for "unclustered" Met
//___ takes into account correlations between params
//___ systcode=0 -default param
//___ systcode=1 -choice of parameterization: Z->ee vs. dipho-sideband
//___ systcode=2(+/-) -Met Model Uncl. En. systematics: mean+/-G
//___ systcode=3(+/-) -Met Model Uncl. En. systematics: sigma+/-G
//___ systcode=4(+/-) -Met Model Uncl. En. systematics: scale+/-G
//___ systcode=5(+/-) -Met Model Uncl. En. systematics: norm+/-G   
void TStnMetModel::GetMyUnclMetResolution(double sumEt, int isData, int ParamSwitch, int systcode, 
			    double &sigmaMx, double &sigmaMy, double &meanMx, double &meanMy,
			    double &normX, double &normY, double &scaleX, double &scaleY) {
  
  double sigmaX_err=0.0;
  double meanX_err=0.0;
  double sigmaY_err=0.0;
  double meanY_err=0.0;
  double normX_err=0.0;
  double normY_err=0.0;
  double scaleX_err=0.0;
  double scaleY_err=0.0;
  double dummypar=0.0;

  if(sumEt<0.0) sumEt=0.0;

//--------------------------------------------------------------------------------
//_______________ ***** met model parametrization.
//                for MetX & MetY: sigma=p0+p1*sqrt(SumEt)
//                for MetX & MetY: mean=p0+p1*SumEt
//                first index is for di-pho sideband ([0]) & central Z->ee ([1]) parametrizations   
  double metmodel_sigmaX[2][2];    // sigmaX=p[0]+p[1]*sqrt(corrSumEt)
  double metmodel_sigmaX_er[2][2];
  double metmodel_sigmaY[2][2];    // sigmaY=p[0]+p[1]*sqrt(corrSumEt) 
  double metmodel_sigmaY_er[2][2];
  double metmodel_meanX[2][2];     // meanX=p[0]+p[1]*corrSumEt
  double metmodel_meanX_er[2][2];
  double metmodel_meanY[2][2];     // meanY=p[0]+p[1]*corrSumEt
  double metmodel_meanY_er[2][2];

  double metmodel_normX[2]; // normX=p[0]
  double metmodel_normX_er[2]; 
  double metmodel_normY[2]; // normY=p[0]
  double metmodel_normY_er[2]; 
  double metmodel_sigmaScaleX[2][2];    // sigmaScaleX=p[0]+p[1]*sqrt(corrSumEt)
  double metmodel_sigmaScaleX_er[2][2];
  double metmodel_sigmaScaleY[2][2];    // sigmaScaleY=p[0]+p[1]*sqrt(corrSumEt) 
  double metmodel_sigmaScaleY_er[2][2];

  //-------------- For unclustered MetModel with Njet(cut)=0
  double metmodel_sigmaXcorr[2]; // correlation coefficient for sigmaX
  double metmodel_sigmaYcorr[2]; // correlation coefficient for sigmaY
  double metmodel_meanXcorr[2]; // correlation coefficient for meanX
  double metmodel_meanYcorr[2]; // correlation coefficient for meanY
  double metmodel_sigmaScaleXcorr[2]; // correlation coefficient for sigmaScaleX
  double metmodel_sigmaScaleYcorr[2]; // correlation coefficient for sigmaScaleY

  if(isData==1)   //---------- Data 
    {
      //____________________________ Data di-pho sideband parameterization for events with Njet15 = 0
      //------------- parametrization after 09/12/07
      metmodel_sigmaX[0][0]=0.82;
      metmodel_sigmaX[0][1]=0.372;
      metmodel_sigmaX_er[0][0]=0.25;
      metmodel_sigmaX_er[0][1]=0.031;
      metmodel_sigmaY[0][0]=0.60;
      metmodel_sigmaY[0][1]=0.387;
      metmodel_sigmaY_er[0][0]=0.25;
      metmodel_sigmaY_er[0][1]=0.022;
      metmodel_meanX[0][0]=-0.022;
      metmodel_meanX[0][1]=0.00647;
      metmodel_meanX_er[0][0]=0.057;
      metmodel_meanX_er[0][1]=0.00076;
      metmodel_meanY[0][0]=-0.016;
      metmodel_meanY[0][1]=-0.00337;
      metmodel_meanY_er[0][0]=0.038;
      metmodel_meanY_er[0][1]=0.00033;
      
      metmodel_sigmaXcorr[0]=-0.911;
      metmodel_sigmaYcorr[0]=-0.906;
      metmodel_meanXcorr[0]=-0.814;
      metmodel_meanYcorr[0]=-0.764;
      
      metmodel_normX[0]=0.080;
      metmodel_normX_er[0]=0.022; 
      metmodel_normY[0]=0.147;
      metmodel_normY_er[0]=0.036; 
      metmodel_sigmaScaleX[0][0]=2.16;  
      metmodel_sigmaScaleX[0][1]=-0.064;
      metmodel_sigmaScaleX_er[0][0]=0.17;
      metmodel_sigmaScaleX_er[0][1]=0.020;
      metmodel_sigmaScaleY[0][0]=1.99; 
      metmodel_sigmaScaleY[0][1]=-0.046;
      metmodel_sigmaScaleY_er[0][0]=0.17;
      metmodel_sigmaScaleY_er[0][1]=0.020;
      metmodel_sigmaScaleXcorr[0]=-0.973;
      metmodel_sigmaScaleYcorr[0]=-0.976;
      
      //------------- Data Z->ee (CC) parameterization (as of 09/12/07)
      metmodel_sigmaX[1][0]=1.03;
      metmodel_sigmaX[1][1]=0.371;
      metmodel_sigmaX_er[1][0]=0.36;
      metmodel_sigmaX_er[1][1]=0.042;
      metmodel_sigmaY[1][0]=1.04;
      metmodel_sigmaY[1][1]=0.389;
      metmodel_sigmaY_er[1][0]=0.32;
      metmodel_sigmaY_er[1][1]=0.034;
      metmodel_meanX[1][0]=-0.048;
      metmodel_meanX[1][1]=0.00674;
      metmodel_meanX_er[1][0]=0.057;
      metmodel_meanX_er[1][1]=0.00073;
      metmodel_meanY[1][0]=-0.117;
      metmodel_meanY[1][1]=-0.00323;
      metmodel_meanY_er[1][0]=0.053;
      metmodel_meanY_er[1][1]=0.00073;
      
      metmodel_sigmaXcorr[1]=-0.921;
      metmodel_sigmaYcorr[1]=-0.903;
      metmodel_meanXcorr[1]=-0.774;
      metmodel_meanYcorr[1]=-0.827;
      
      metmodel_normX[1]=0.281;
      metmodel_normX_er[1]=0.076; 
      metmodel_normY[1]=0.235;
      metmodel_normY_er[1]=0.078; 
      metmodel_sigmaScaleX[1][0]=1.94;  
      metmodel_sigmaScaleX[1][1]=-0.051;
      metmodel_sigmaScaleX_er[1][0]=0.11;
      metmodel_sigmaScaleX_er[1][1]=0.013;
      metmodel_sigmaScaleY[1][0]=2.19; 
      metmodel_sigmaScaleY[1][1]=-0.079;
      metmodel_sigmaScaleY_er[1][0]=0.16;
      metmodel_sigmaScaleY_er[1][1]=0.019;
      metmodel_sigmaScaleXcorr[1]=-0.944;
      metmodel_sigmaScaleYcorr[1]=-0.965;
    }
  else //--------- MC 
    {
      //____________________________ Pythia di-pho signal parametrization for events with Njet15 = 0
      //------------- parametrization after 01/31/07
      metmodel_sigmaX[0][0]=0.15;
      metmodel_sigmaX[0][1]=0.502;
      metmodel_sigmaX_er[0][0]=0.46;
      metmodel_sigmaX_er[0][1]=0.073;
      metmodel_sigmaY[0][0]=0.11;
      metmodel_sigmaY[0][1]=0.526;
      metmodel_sigmaY_er[0][0]=0.37;
      metmodel_sigmaY_er[0][1]=0.045;
      metmodel_meanX[0][0]=0.044;
      metmodel_meanX[0][1]=0.0068;
      metmodel_meanX_er[0][0]=0.048;
      metmodel_meanX_er[0][1]=0.0013;
      metmodel_meanY[0][0]=-0.195;
      metmodel_meanY[0][1]=-0.00927;
      metmodel_meanY_er[0][0]=0.050;
      metmodel_meanY_er[0][1]=0.00096;
      
      metmodel_sigmaXcorr[0]=-0.943;
      metmodel_sigmaYcorr[0]=-0.945;
      metmodel_meanXcorr[0]=-0.812;
      metmodel_meanYcorr[0]=-0.800;
      
      metmodel_normX[0]=0.101;
      metmodel_normX_er[0]=0.024; 
      metmodel_normY[0]=0.096;
      metmodel_normY_er[0]=0.029; 
      metmodel_sigmaScaleX[0][0]=1.98;  
      metmodel_sigmaScaleX[0][1]=-0.052;
      metmodel_sigmaScaleX_er[0][0]=0.16;
      metmodel_sigmaScaleX_er[0][1]=0.026;
      metmodel_sigmaScaleY[0][0]=2.18; 
      metmodel_sigmaScaleY[0][1]=-0.090;
      metmodel_sigmaScaleY_er[0][0]=0.33;
      metmodel_sigmaScaleY_er[0][1]=0.053;
      metmodel_sigmaScaleXcorr[0]=-0.969;
      metmodel_sigmaScaleYcorr[0]=-0.988;
      
      //------------- Pythia Z->ee (CC) parametrization 
      metmodel_sigmaX[1][0]=0.98;
      metmodel_sigmaX[1][1]=0.457;
      metmodel_sigmaX_er[1][0]=0.25;
      metmodel_sigmaX_er[1][1]=0.027;
      metmodel_sigmaY[1][0]=1.66;
      metmodel_sigmaY[1][1]=0.367;
      metmodel_sigmaY_er[1][0]=0.17;
      metmodel_sigmaY_er[1][1]=0.027;
      metmodel_meanX[1][0]=0.155;
      metmodel_meanX[1][1]=0.00137;
      metmodel_meanX_er[1][0]=0.020;
      metmodel_meanX_er[1][1]=0.00035;
      metmodel_meanY[1][0]=-0.179;
      metmodel_meanY[1][1]=-0.00333;
      metmodel_meanY_er[1][0]=0.026;
      metmodel_meanY_er[1][1]=0.00044;
      
      metmodel_sigmaXcorr[1]=-0.893;
      metmodel_sigmaYcorr[1]=-0.857;
      metmodel_meanXcorr[1]=-0.857;
      metmodel_meanYcorr[1]=-0.807;
      
      metmodel_normX[1]=0.110;
      metmodel_normX_er[1]=0.025; 
      metmodel_normY[1]=0.134;
      metmodel_normY_er[1]=0.029; 
      metmodel_sigmaScaleX[1][0]=1.755;  
      metmodel_sigmaScaleX[1][1]=-0.046;
      metmodel_sigmaScaleX_er[1][0]=0.091;
      metmodel_sigmaScaleX_er[1][1]=0.011;
      metmodel_sigmaScaleY[1][0]=1.697; 
      metmodel_sigmaScaleY[1][1]=-0.035;
      metmodel_sigmaScaleY_er[1][0]=0.083;
      metmodel_sigmaScaleY_er[1][1]=0.011;
      metmodel_sigmaScaleXcorr[1]=-0.966;
      metmodel_sigmaScaleYcorr[1]=-0.961;
    }

//___ systcode=4(+/-) -Met Model Uncl. En. systematics: scale+/-G
//___ systcode=5(+/-) -Met Model Uncl. En. systematics: norm+/-G   
//_____________________ systcode=0 -default param
  if(systcode==0) 
    {
      //___________ X-axis
      //-sigma
      sigmaMx=metmodel_sigmaX[ParamSwitch][0]+metmodel_sigmaX[ParamSwitch][1]*sqrt(sumEt);
      if(sigmaMx<0.0) sigmaMx=0.0;      
      //-mean
      meanMx=metmodel_meanX[ParamSwitch][0]+metmodel_meanX[ParamSwitch][1]*sumEt;
      //-norm
      normX=metmodel_normX[ParamSwitch];
      if(normX<0.0) normX=0.0;
      //-scale
      scaleX=metmodel_sigmaScaleX[ParamSwitch][0]+metmodel_sigmaScaleX[ParamSwitch][1]*sqrt(sumEt);
      if(scaleX<1.0) scaleX=1.0;
      //___________ Y-axis
      //-sigma
      sigmaMy=metmodel_sigmaY[ParamSwitch][0]+metmodel_sigmaY[ParamSwitch][1]*sqrt(sumEt);
      if(sigmaMy<0.0) sigmaMy=0.0;
      //-mean
      meanMy=metmodel_meanY[ParamSwitch][0]+metmodel_meanY[ParamSwitch][1]*sumEt;
      //-norm
      normY=metmodel_normY[ParamSwitch];
      if(normY<0.0) normY=0.0;
      //-scale
      scaleY=metmodel_sigmaScaleY[ParamSwitch][0]+metmodel_sigmaScaleY[ParamSwitch][1]*sqrt(sumEt);
      if(scaleY<1.0) scaleY=1.0;
    }
  //________________________ systcode=1 -choice of parameterization: Z->ee vs. dipho-sideband
  if(abs(systcode)==1)
    {
      //___________ X-axis
      //-sigma
      sigmaMx=metmodel_sigmaX[1-ParamSwitch][0]+metmodel_sigmaX[1-ParamSwitch][1]*sqrt(sumEt);
      if(sigmaMx<0.0) sigmaMx=0.0;      
      //-mean
      meanMx=metmodel_meanX[1-ParamSwitch][0]+metmodel_meanX[1-ParamSwitch][1]*sumEt;
      //-norm
      normX=metmodel_normX[1-ParamSwitch];
      if(normX<0.0) normX=0.0;
      //-scale
      scaleX=metmodel_sigmaScaleX[1-ParamSwitch][0]+metmodel_sigmaScaleX[1-ParamSwitch][1]*sqrt(sumEt);
      if(scaleX<1.0) scaleX=1.0;
      //___________ Y-axis
      //-sigma
      sigmaMy=metmodel_sigmaY[1-ParamSwitch][0]+metmodel_sigmaY[1-ParamSwitch][1]*sqrt(sumEt);
      if(sigmaMy<0.0) sigmaMy=0.0;
      //-mean
      meanMy=metmodel_meanY[1-ParamSwitch][0]+metmodel_meanY[1-ParamSwitch][1]*sumEt;
      //-norm
      normY=metmodel_normY[1-ParamSwitch];
      if(normY<0.0) normY=0.0;
      //-scale
      scaleY=metmodel_sigmaScaleY[1-ParamSwitch][0]+metmodel_sigmaScaleY[1-ParamSwitch][1]*sqrt(sumEt);
      if(scaleY<1.0) scaleY=1.0;
    }
  
  //___________________ calculating systematics
  if(abs(systcode)>1)
    {
      //___________ X-axis
      //-sigma
      dummypar=metmodel_sigmaX_er[ParamSwitch][0]*metmodel_sigmaX_er[ParamSwitch][0]
	+ metmodel_sigmaX_er[ParamSwitch][1]*metmodel_sigmaX_er[ParamSwitch][1]*sumEt
	+ 2.0*metmodel_sigmaXcorr[ParamSwitch]*metmodel_sigmaX_er[ParamSwitch][0]*metmodel_sigmaX_er[ParamSwitch][1]*sqrt(sumEt);
      if(dummypar<0.0) dummypar=0.0;
      sigmaX_err=sqrt(dummypar);
      //-mean
      dummypar=metmodel_meanX_er[ParamSwitch][0]*metmodel_meanX_er[ParamSwitch][0]
	+ metmodel_meanX_er[ParamSwitch][1]*metmodel_meanX_er[ParamSwitch][1]*sumEt*sumEt
	+ 2.0*metmodel_meanXcorr[ParamSwitch]*metmodel_meanX_er[ParamSwitch][0]*metmodel_meanX_er[ParamSwitch][1]*sumEt;
      if(dummypar<0.0) dummypar=0.0;
      meanX_err=sqrt(dummypar);
      //-norm
      dummypar=metmodel_normX_er[ParamSwitch]*metmodel_normX_er[ParamSwitch];
      normX_err=sqrt(dummypar);
      //-scale
      dummypar=metmodel_sigmaScaleX_er[ParamSwitch][0]*metmodel_sigmaScaleX_er[ParamSwitch][0]
	+ metmodel_sigmaScaleX_er[ParamSwitch][1]*metmodel_sigmaScaleX_er[ParamSwitch][1]*sumEt
	+ 2.0*metmodel_sigmaScaleXcorr[ParamSwitch]*metmodel_sigmaScaleX_er[ParamSwitch][0]*metmodel_sigmaScaleX_er[ParamSwitch][1]*sqrt(sumEt);
      if(dummypar<0.0) dummypar=0.0;
      scaleX_err=sqrt(dummypar);
      //___________ Y-axis
      //-sigma
      dummypar=metmodel_sigmaY_er[ParamSwitch][0]*metmodel_sigmaY_er[ParamSwitch][0]
	+ metmodel_sigmaY_er[ParamSwitch][1]*metmodel_sigmaY_er[ParamSwitch][1]*sumEt
	+ 2.0*metmodel_sigmaYcorr[ParamSwitch]*metmodel_sigmaY_er[ParamSwitch][0]*metmodel_sigmaY_er[ParamSwitch][1]*sqrt(sumEt);
      if(dummypar<0.0) dummypar=0.0;
      sigmaY_err=sqrt(dummypar);
      //-mean
      dummypar=metmodel_meanY_er[ParamSwitch][0]*metmodel_meanY_er[ParamSwitch][0]
	+ metmodel_meanY_er[ParamSwitch][1]*metmodel_meanY_er[ParamSwitch][1]*sumEt*sumEt
	+ 2.0*metmodel_meanYcorr[ParamSwitch]*metmodel_meanY_er[ParamSwitch][0]*metmodel_meanY_er[ParamSwitch][1]*sumEt;
      if(dummypar<0.0) dummypar=0.0;
      meanY_err=sqrt(dummypar);  
      //-norm
      dummypar=metmodel_normY_er[ParamSwitch]*metmodel_normY_er[ParamSwitch];
      normY_err=sqrt(dummypar);
      //-scale
      dummypar=metmodel_sigmaScaleY_er[ParamSwitch][0]*metmodel_sigmaScaleY_er[ParamSwitch][0]
	+ metmodel_sigmaScaleY_er[ParamSwitch][1]*metmodel_sigmaScaleY_er[ParamSwitch][1]*sumEt
	+ 2.0*metmodel_sigmaScaleYcorr[ParamSwitch]*metmodel_sigmaScaleY_er[ParamSwitch][0]*metmodel_sigmaScaleY_er[ParamSwitch][1]*sqrt(sumEt);
      if(dummypar<0.0) dummypar=0.0;
      scaleY_err=sqrt(dummypar);
    }
  //____________________ systcode=2(+/-) -Met Model Uncl. En. systematics: mean+/-G
  if(abs(systcode)==2)
    {
      int int_systcode=(systcode>0) ? 1 : -1 ;
      //___________ X-axis
      //-sigma
      sigmaMx=metmodel_sigmaX[ParamSwitch][0]+metmodel_sigmaX[ParamSwitch][1]*sqrt(sumEt);
      if(sigmaMx<0.0) sigmaMx=0.0;      
      //-mean
      meanMx=metmodel_meanX[ParamSwitch][0]+metmodel_meanX[ParamSwitch][1]*sumEt+int_systcode*meanX_err;
      //-norm
      normX=metmodel_normX[ParamSwitch];
      if(normX<0.0) normX=0.0;
      //-scale
      scaleX=metmodel_sigmaScaleX[ParamSwitch][0]+metmodel_sigmaScaleX[ParamSwitch][1]*sqrt(sumEt);
      if(scaleX<1.0) scaleX=1.0;
      //___________ Y-axis
      //-sigma
      sigmaMy=metmodel_sigmaY[ParamSwitch][0]+metmodel_sigmaY[ParamSwitch][1]*sqrt(sumEt);
      if(sigmaMy<0.0) sigmaMy=0.0;
      //-mean
      meanMy=metmodel_meanY[ParamSwitch][0]+metmodel_meanY[ParamSwitch][1]*sumEt+int_systcode*meanY_err;
      //-norm
      normY=metmodel_normY[ParamSwitch];
      if(normY<0.0) normY=0.0;
      //-scale
      scaleY=metmodel_sigmaScaleY[ParamSwitch][0]+metmodel_sigmaScaleY[ParamSwitch][1]*sqrt(sumEt);
      if(scaleY<1.0) scaleY=1.0;
    }
  //____________________ systcode=3(+/-) -Met Model Uncl. En. systematics: sigma+/-G
  if(abs(systcode)==3)
    {
      int int_systcode=(systcode>0) ? 1 : -1 ;
      //___________ X-axis
      //-sigma
      sigmaMx=metmodel_sigmaX[ParamSwitch][0]+metmodel_sigmaX[ParamSwitch][1]*sqrt(sumEt)+int_systcode*sigmaX_err;
      if(sigmaMx<0.0) sigmaMx=0.0;      
      //-mean
      meanMx=metmodel_meanX[ParamSwitch][0]+metmodel_meanX[ParamSwitch][1]*sumEt;
      //-norm
      normX=metmodel_normX[ParamSwitch];
      if(normX<0.0) normX=0.0;
      //-scale
      scaleX=metmodel_sigmaScaleX[ParamSwitch][0]+metmodel_sigmaScaleX[ParamSwitch][1]*sqrt(sumEt);
      if(scaleX<1.0) scaleX=1.0;
      //___________ Y-axis
      //-sigma
      sigmaMy=metmodel_sigmaY[ParamSwitch][0]+metmodel_sigmaY[ParamSwitch][1]*sqrt(sumEt)+int_systcode*sigmaY_err;
      if(sigmaMy<0.0) sigmaMy=0.0;
      //-mean
      meanMy=metmodel_meanY[ParamSwitch][0]+metmodel_meanY[ParamSwitch][1]*sumEt;
      //-norm
      normY=metmodel_normY[ParamSwitch];
      if(normY<0.0) normY=0.0;
      //-scale
      scaleY=metmodel_sigmaScaleY[ParamSwitch][0]+metmodel_sigmaScaleY[ParamSwitch][1]*sqrt(sumEt);
      if(scaleY<1.0) scaleY=1.0;
    }
  //____________________ systcode=4(+/-) -Met Model Uncl. En. systematics: scale+/-G
  if(abs(systcode)==4)
    {
      int int_systcode=(systcode>0) ? 1 : -1 ;
      //___________ X-axis
      //-sigma
      sigmaMx=metmodel_sigmaX[ParamSwitch][0]+metmodel_sigmaX[ParamSwitch][1]*sqrt(sumEt);
      if(sigmaMx<0.0) sigmaMx=0.0;      
      //-mean
      meanMx=metmodel_meanX[ParamSwitch][0]+metmodel_meanX[ParamSwitch][1]*sumEt;
      //-norm
      normX=metmodel_normX[ParamSwitch];
      if(normX<0.0) normX=0.0;
      //-scale
      scaleX=metmodel_sigmaScaleX[ParamSwitch][0]+metmodel_sigmaScaleX[ParamSwitch][1]*sqrt(sumEt)+int_systcode*scaleX_err;
      if(scaleX<1.0) scaleX=1.0;
      //___________ Y-axis
      //-sigma
      sigmaMy=metmodel_sigmaY[ParamSwitch][0]+metmodel_sigmaY[ParamSwitch][1]*sqrt(sumEt);
      if(sigmaMy<0.0) sigmaMy=0.0;
      //-mean
      meanMy=metmodel_meanY[ParamSwitch][0]+metmodel_meanY[ParamSwitch][1]*sumEt;
      //-norm
      normY=metmodel_normY[ParamSwitch];
      if(normY<0.0) normY=0.0;
      //-scale
      scaleY=metmodel_sigmaScaleY[ParamSwitch][0]+metmodel_sigmaScaleY[ParamSwitch][1]*sqrt(sumEt)+int_systcode*scaleY_err;
      if(scaleY<1.0) scaleY=1.0;
    }
  //____________________ systcode=5(+/-) -Met Model Uncl. En. systematics: scale+/-G
  if(abs(systcode)==5)
    {
      int int_systcode=(systcode>0) ? 1 : -1 ;
      //___________ X-axis
      //-sigma
      sigmaMx=metmodel_sigmaX[ParamSwitch][0]+metmodel_sigmaX[ParamSwitch][1]*sqrt(sumEt);
      if(sigmaMx<0.0) sigmaMx=0.0;      
      //-mean
      meanMx=metmodel_meanX[ParamSwitch][0]+metmodel_meanX[ParamSwitch][1]*sumEt;
      //-norm
      normX=metmodel_normX[ParamSwitch]+int_systcode*normX_err;
      if(normX<0.0) normX=0.0;
      //-scale
      scaleX=metmodel_sigmaScaleX[ParamSwitch][0]+metmodel_sigmaScaleX[ParamSwitch][1]*sqrt(sumEt);
      if(scaleX<1.0) scaleX=1.0;
      //___________ Y-axis
      //-sigma
      sigmaMy=metmodel_sigmaY[ParamSwitch][0]+metmodel_sigmaY[ParamSwitch][1]*sqrt(sumEt);
      if(sigmaMy<0.0) sigmaMy=0.0;
      //-mean
      meanMy=metmodel_meanY[ParamSwitch][0]+metmodel_meanY[ParamSwitch][1]*sumEt;
      //-norm
      normY=metmodel_normY[ParamSwitch]+int_systcode*normY_err;
      if(normY<0.0) normY=0.0;
      //-scale
      scaleY=metmodel_sigmaScaleY[ParamSwitch][0]+metmodel_sigmaScaleY[ParamSwitch][1]*sqrt(sumEt);
      if(scaleY<1.0) scaleY=1.0;
    }
  return;
}



//________________ returns correlation coefficients for JER params
double TStnMetModel::JERCorrCoeff(int i, int j)
{
  double corr_param[15][11]; // [eta][param]
  corr_param[0][0]=-0.953359; // Gauss mean: c01
  corr_param[0][1]=-0.883424; // Gauss mean: c02
  corr_param[0][2]=0.783976; // Gauss mean: c12
  corr_param[0][3]=-0.769484; // Gauss sigma: c01
  corr_param[0][4]=-0.916005; // Landau mpv: c01
  corr_param[0][5]=-0.920646; // Landau mpv: c02
  corr_param[0][6]=0.774913; // Landau mpv: c12
  corr_param[0][7]=-0.798127; // Landau sigma: c01
  corr_param[0][8]=-0.979048; // norm: c01
  corr_param[0][9]=0.914083; // norm: c02
  corr_param[0][10]=-0.972488; // norm: c12
  corr_param[1][0]=-0.930448; // Gauss mean: c01
  corr_param[1][1]=-0.893613; // Gauss mean: c02
  corr_param[1][2]=0.743464; // Gauss mean: c12
  corr_param[1][3]=-0.855253; // Gauss sigma: c01
  corr_param[1][4]=-0.93506; // Landau mpv: c01
  corr_param[1][5]=-0.790079; // Landau mpv: c02
  corr_param[1][6]=0.633928; // Landau mpv: c12
  corr_param[1][7]=-0.791488; // Landau sigma: c01
  corr_param[1][8]=-0.979059; // norm: c01
  corr_param[1][9]=0.935759; // norm: c02
  corr_param[1][10]=-0.978165; // norm: c12
  corr_param[2][0]=-0.933467; // Gauss mean: c01
  corr_param[2][1]=-0.903664; // Gauss mean: c02
  corr_param[2][2]=0.769636; // Gauss mean: c12
  corr_param[2][3]=-0.834006; // Gauss sigma: c01
  corr_param[2][4]=-0.941456; // Landau mpv: c01
  corr_param[2][5]=-0.863316; // Landau mpv: c02
  corr_param[2][6]=0.730447; // Landau mpv: c12
  corr_param[2][7]=-0.795485; // Landau sigma: c01
  corr_param[2][8]=-0.973561; // norm: c01
  corr_param[2][9]=0.922559; // norm: c02
  corr_param[2][10]=-0.979269; // norm: c12
  corr_param[3][0]=-0.92146; // Gauss mean: c01
  corr_param[3][1]=-0.875047; // Gauss mean: c02
  corr_param[3][2]=0.720247; // Gauss mean: c12
  corr_param[3][3]=-0.820492; // Gauss sigma: c01
  corr_param[3][4]=-0.932934; // Landau mpv: c01
  corr_param[3][5]=-0.806167; // Landau mpv: c02
  corr_param[3][6]=0.666462; // Landau mpv: c12
  corr_param[3][7]=-0.752656; // Landau sigma: c01
  corr_param[3][8]=-0.97604; // norm: c01
  corr_param[3][9]=0.916931; // norm: c02
  corr_param[3][10]=-0.970166; // norm: c12
  corr_param[4][0]=-0.921373; // Gauss mean: c01
  corr_param[4][1]=-0.865997; // Gauss mean: c02
  corr_param[4][2]=0.713135; // Gauss mean: c12
  corr_param[4][3]=-0.812212; // Gauss sigma: c01
  corr_param[4][4]=-0.933552; // Landau mpv: c01
  corr_param[4][5]=-0.826206; // Landau mpv: c02
  corr_param[4][6]=0.696959; // Landau mpv: c12
  corr_param[4][7]=-0.731748; // Landau sigma: c01
  corr_param[4][8]=-0.972365; // norm: c01
  corr_param[4][9]=0.900603; // norm: c02
  corr_param[4][10]=-0.96338; // norm: c12
  corr_param[5][0]=-0.945056; // Gauss mean: c01
  corr_param[5][1]=-0.918674; // Gauss mean: c02
  corr_param[5][2]=0.808282; // Gauss mean: c12
  corr_param[5][3]=-0.811252; // Gauss sigma: c01
  corr_param[5][4]=-0.941345; // Landau mpv: c01
  corr_param[5][5]=-0.908828; // Landau mpv: c02
  corr_param[5][6]=0.805962; // Landau mpv: c12
  corr_param[5][7]=-0.737379; // Landau sigma: c01
  corr_param[5][8]=-0.974342; // norm: c01
  corr_param[5][9]=0.878683; // norm: c02
  corr_param[5][10]=-0.952568; // norm: c12
  corr_param[6][0]=-0.863692; // Gauss mean: c01
  corr_param[6][1]=-0.898828; // Gauss mean: c02
  corr_param[6][2]=0.683443; // Gauss mean: c12
  corr_param[6][3]=-0.864674; // Gauss sigma: c01
  corr_param[6][4]=-0.910948; // Landau mpv: c01
  corr_param[6][5]=-0.878636; // Landau mpv: c02
  corr_param[6][6]=0.709797; // Landau mpv: c12
  corr_param[6][7]=-0.790837; // Landau sigma: c01
  corr_param[6][8]=-0.976028; // norm: c01
  corr_param[6][9]=0.920123; // norm: c02
  corr_param[6][10]=-0.976345; // norm: c12
  corr_param[7][0]=-0.901329; // Gauss mean: c01
  corr_param[7][1]=-0.868814; // Gauss mean: c02
  corr_param[7][2]=0.664341; // Gauss mean: c12
  corr_param[7][3]=-0.879818; // Gauss sigma: c01
  corr_param[7][4]=-0.910909; // Landau mpv: c01
  corr_param[7][5]=-0.848758; // Landau mpv: c02
  corr_param[7][6]=0.654596; // Landau mpv: c12
  corr_param[7][7]=-0.854805; // Landau sigma: c01
  corr_param[7][8]=-0.98475; // norm: c01
  corr_param[7][9]=0.944593; // norm: c02
  corr_param[7][10]=-0.983111; // norm: c12
  corr_param[8][0]=-0.93546; // Gauss mean: c01
  corr_param[8][1]=-0.919387; // Gauss mean: c02
  corr_param[8][2]=0.775089; // Gauss mean: c12
  corr_param[8][3]=-0.885525; // Gauss sigma: c01
  corr_param[8][4]=-0.902932; // Landau mpv: c01
  corr_param[8][5]=-0.858334; // Landau mpv: c02
  corr_param[8][6]=0.664962; // Landau mpv: c12
  corr_param[8][7]=-0.83488; // Landau sigma: c01
  corr_param[8][8]=-0.989491; // norm: c01
  corr_param[8][9]=0.946333; // norm: c02
  corr_param[8][10]=-0.980041; // norm: c12
  corr_param[9][0]=-0.958622; // Gauss mean: c01
  corr_param[9][1]=-0.92783; // Gauss mean: c02
  corr_param[9][2]=0.813899; // Gauss mean: c12
  corr_param[9][3]=-0.918357; // Gauss sigma: c01
  corr_param[9][4]=-0.946212; // Landau mpv: c01
  corr_param[9][5]=-0.874383; // Landau mpv: c02
  corr_param[9][6]=0.732471; // Landau mpv: c12
  corr_param[9][7]=-0.890866; // Landau sigma: c01
  corr_param[9][8]=-0.993186; // norm: c01
  corr_param[9][9]=0.965285; // norm: c02
  corr_param[9][10]=-0.987298; // norm: c12
  corr_param[10][0]=-0.959116; // Gauss mean: c01
  corr_param[10][1]=-0.934235; // Gauss mean: c02
  corr_param[10][2]=0.822224; // Gauss mean: c12
  corr_param[10][3]=-0.929578; // Gauss sigma: c01
  corr_param[10][4]=-0.945125; // Landau mpv: c01
  corr_param[10][5]=-0.889286; // Landau mpv: c02
  corr_param[10][6]=0.741957; // Landau mpv: c12
  corr_param[10][7]=-0.900695; // Landau sigma: c01
  corr_param[10][8]=-0.992615; // norm: c01
  corr_param[10][9]=0.967849; // norm: c02
  corr_param[10][10]=-0.989424; // norm: c12
  corr_param[11][0]=-0.961204; // Gauss mean: c01
  corr_param[11][1]=-0.933701; // Gauss mean: c02
  corr_param[11][2]=0.822536; // Gauss mean: c12
  corr_param[11][3]=-0.938163; // Gauss sigma: c01
  corr_param[11][4]=-0.95; // Landau mpv: c01
  corr_param[11][5]=-0.875959; // Landau mpv: c02
  corr_param[11][6]=0.737455; // Landau mpv: c12
  corr_param[11][7]=-0.895806; // Landau sigma: c01
  corr_param[11][8]=-0.993612; // norm: c01
  corr_param[11][9]=0.968624; // norm: c02
  corr_param[11][10]=-0.988782; // norm: c12
  corr_param[12][0]=-0.961494; // Gauss mean: c01
  corr_param[12][1]=-0.888003; // Gauss mean: c02
  corr_param[12][2]=0.756614; // Gauss mean: c12
  corr_param[12][3]=-0.940026; // Gauss sigma: c01
  corr_param[12][4]=-0.971454; // Landau mpv: c01
  corr_param[12][5]=-0.907369; // Landau mpv: c02
  corr_param[12][6]=0.819359; // Landau mpv: c12
  corr_param[12][7]=-0.890781; // Landau sigma: c01
  corr_param[12][8]=-0.995606; // norm: c01
  corr_param[12][9]=0.976755; // norm: c02
  corr_param[12][10]=-0.991439; // norm: c12
  corr_param[13][0]=-0.979353; // Gauss mean: c01
  corr_param[13][1]=-0.901815; // Gauss mean: c02
  corr_param[13][2]=0.819153; // Gauss mean: c12
  corr_param[13][3]=-0.906452; // Gauss sigma: c01
  corr_param[13][4]=-0.977947; // Landau mpv: c01
  corr_param[13][5]=-0.902621; // Landau mpv: c02
  corr_param[13][6]=0.821221; // Landau mpv: c12
  corr_param[13][7]=-0.907155; // Landau sigma: c01
  corr_param[13][8]=-0.996171; // norm: c01
  corr_param[13][9]=0.98616; // norm: c02
  corr_param[13][10]=-0.996056; // norm: c12
  corr_param[14][0]=0; // Gauss mean: c01
  corr_param[14][1]=0; // Gauss mean: c02
  corr_param[14][2]=0; // Gauss mean: c12
  corr_param[14][3]=0; // Gauss sigma: c01
  corr_param[14][4]=-0.993129; // Landau mpv: c01
  corr_param[14][5]=-0.979168; // Landau mpv: c02
  corr_param[14][6]=0.957156; // Landau mpv: c12
  corr_param[14][7]=-0.895607; // Landau sigma: c01
  corr_param[14][8]=0; // norm: c01
  corr_param[14][9]=0; // norm: c02
  corr_param[14][10]=0; // norm: c12

  double value=0.0;
  if(i<15 && j<11) value=corr_param[i][j];

  return value;
}


//________________ returns JER param or its uncertainty
double TStnMetModel::JERparam(int code, int i, int j)
{
  double jer_param[15][13]; // [eta][param] 
  double jer_param_er[15][13]; // [eta][param]
  double value=0.0;

  jer_param[0][0]=-0.00505533;  jer_param_er[0][0]=0.00366313;
  jer_param[0][1]=-7.37332E-5;  jer_param_er[0][1]=1.02366E-5;
  jer_param[0][2]=5.06548;  jer_param_er[0][2]=0.198112;
  jer_param[0][3]=1.03604;  jer_param_er[0][3]=0.0162871;
  jer_param[0][4]=0.000363984;  jer_param_er[0][4]=9.99763E-5;
  jer_param[0][5]=0.0404354;  jer_param_er[0][5]=0.00394954;
  jer_param[0][6]=7.70406E-5;  jer_param_er[0][6]=1.80066E-5;
  jer_param[0][7]=-1.05943;  jer_param_er[0][7]=0.148324;
  jer_param[0][8]=0.157408;  jer_param_er[0][8]=0.00948233;
  jer_param[0][9]=0.00745406;  jer_param_er[0][9]=0.000132492;
  jer_param[0][10]=-4.48426;  jer_param_er[0][10]=0.997815;
  jer_param[0][11]=1.14401;  jer_param_er[0][11]=0.269919;
  jer_param[0][12]=0.0459092;  jer_param_er[0][12]=0.0168848;
  jer_param[1][0]=0.0928237;  jer_param_er[1][0]=0.00248223;
  jer_param[1][1]=-0.000201847;  jer_param_er[1][1]=8.82907E-6;
  jer_param[1][2]=0.177168;  jer_param_er[1][2]=0.137301;
  jer_param[1][3]=0.670063;  jer_param_er[1][3]=0.0119732;
  jer_param[1][4]=0.00308811;  jer_param_er[1][4]=9.34083E-5;
  jer_param[1][5]=0.000336138;  jer_param_er[1][5]=0.00152807;
  jer_param[1][6]=5.34059E-6;  jer_param_er[1][6]=5.37503E-6;
  jer_param[1][7]=-1.15979;  jer_param_er[1][7]=0.0761452;
  jer_param[1][8]=0.137285;  jer_param_er[1][8]=0.00491652;
  jer_param[1][9]=0.0020658;  jer_param_er[1][9]=3.46802E-5;
  jer_param[1][10]=-14.7515;  jer_param_er[1][10]=0.612517;
  jer_param[1][11]=4.80581;  jer_param_er[1][11]=0.195418;
  jer_param[1][12]=-0.127518;  jer_param_er[1][12]=0.0109911;
  jer_param[2][0]=0.100634;  jer_param_er[2][0]=0.00241696;
  jer_param[2][1]=-0.000253429;  jer_param_er[2][1]=8.43763E-6;
  jer_param[2][2]=-0.688927;  jer_param_er[2][2]=0.125514;
  jer_param[2][3]=0.655545;  jer_param_er[2][3]=0.0103474;
  jer_param[2][4]=0.00300742;  jer_param_er[2][4]=8.73698E-5;
  jer_param[2][5]=0.0105033;  jer_param_er[2][5]=0.00178783;
  jer_param[2][6]=-3.76786E-5;  jer_param_er[2][6]=5.83178E-6;
  jer_param[2][7]=-1.85956;  jer_param_er[2][7]=0.096736;
  jer_param[2][8]=0.139571;  jer_param_er[2][8]=0.00510727;
  jer_param[2][9]=0.00227237;  jer_param_er[2][9]=3.75247E-5;
  jer_param[2][10]=-19.2885;  jer_param_er[2][10]=1.05958;
  jer_param[2][11]=5.66883;  jer_param_er[2][11]=0.262952;
  jer_param[2][12]=-0.161499;  jer_param_er[2][12]=0.013508;
  jer_param[3][0]=0.0952545;  jer_param_er[3][0]=0.00263113;
  jer_param[3][1]=-0.000189574;  jer_param_er[3][1]=9.63887E-6;
  jer_param[3][2]=0.236287;  jer_param_er[3][2]=0.12846;
  jer_param[3][3]=0.731283;  jer_param_er[3][3]=0.0128143;
  jer_param[3][4]=0.00341218;  jer_param_er[3][4]=0.00011154;
  jer_param[3][5]=0.0258263;  jer_param_er[3][5]=0.0016674;
  jer_param[3][6]=-4.83423E-5;  jer_param_er[3][6]=5.80565E-6;
  jer_param[3][7]=-2.22975;  jer_param_er[3][7]=0.0737835;
  jer_param[3][8]=0.179229;  jer_param_er[3][8]=0.00515981;
  jer_param[3][9]=0.0024204;  jer_param_er[3][9]=4.113E-5;
  jer_param[3][10]=-10.4553;  jer_param_er[3][10]=0.561997;
  jer_param[3][11]=3.36745;  jer_param_er[3][11]=0.171045;
  jer_param[3][12]=-0.0618917;  jer_param_er[3][12]=0.00980421;
  jer_param[4][0]=0.115229;  jer_param_er[4][0]=0.0027871;
  jer_param[4][1]=-0.000157746;  jer_param_er[4][1]=9.95943E-6;
  jer_param[4][2]=-0.0905143;  jer_param_er[4][2]=0.134539;
  jer_param[4][3]=0.795373;  jer_param_er[4][3]=0.016243;
  jer_param[4][4]=0.00635758;  jer_param_er[4][4]=0.000146956;
  jer_param[4][5]=0.0345064;  jer_param_er[4][5]=0.00216488;
  jer_param[4][6]=-5.72569E-5;  jer_param_er[4][6]=7.07075E-6;
  jer_param[4][7]=-2.21206;  jer_param_er[4][7]=0.0958281;
  jer_param[4][8]=0.212192;  jer_param_er[4][8]=0.00725209;
  jer_param[4][9]=0.00349351;  jer_param_er[4][9]=6.06418E-5;
  jer_param[4][10]=-10.6021;  jer_param_er[4][10]=0.660483;
  jer_param[4][11]=3.23948;  jer_param_er[4][11]=0.191871;
  jer_param[4][12]=-0.0236045;  jer_param_er[4][12]=0.0106671;
  jer_param[5][0]=0.0900639;  jer_param_er[5][0]=0.00267652;
  jer_param[5][1]=-8.17704E-5;  jer_param_er[5][1]=7.7302E-6;
  jer_param[5][2]=-0.567021;  jer_param_er[5][2]=0.160166;
  jer_param[5][3]=0.706037;  jer_param_er[5][3]=0.0137402;
  jer_param[5][4]=0.0066424;  jer_param_er[5][4]=0.000110101;
  jer_param[5][5]=0.0161473;  jer_param_er[5][5]=0.00313166;
  jer_param[5][6]=-3.72664E-5;  jer_param_er[5][6]=8.95092E-6;
  jer_param[5][7]=-3.40086;  jer_param_er[5][7]=0.154954;
  jer_param[5][8]=0.163174;  jer_param_er[5][8]=0.00755835;
  jer_param[5][9]=0.00291679;  jer_param_er[5][9]=6.87186E-5;
  jer_param[5][10]=-10.8209;  jer_param_er[5][10]=1.39492;
  jer_param[5][11]=2.05925;  jer_param_er[5][11]=0.352187;
  jer_param[5][12]=0.184539;  jer_param_er[5][12]=0.0192903;
  jer_param[6][0]=0.110669;  jer_param_er[6][0]=0.00336557;
  jer_param[6][1]=-0.000302301;  jer_param_er[6][1]=1.15245E-5;
  jer_param[6][2]=-1.70381;  jer_param_er[6][2]=0.17128;
  jer_param[6][3]=0.84678;  jer_param_er[6][3]=0.0199186;
  jer_param[6][4]=0.00418531;  jer_param_er[6][4]=0.00022283;
  jer_param[6][5]=-0.00347175;  jer_param_er[6][5]=0.0027646;
  jer_param[6][6]=8.50977E-5;  jer_param_er[6][6]=1.10633E-5;
  jer_param[6][7]=-3.2283;  jer_param_er[6][7]=0.126136;
  jer_param[6][8]=0.0758951;  jer_param_er[6][8]=0.00708642;
  jer_param[6][9]=0.00415232;  jer_param_er[6][9]=8.11982E-5;
  jer_param[6][10]=-16.6687;  jer_param_er[6][10]=1.09606;
  jer_param[6][11]=5.18449;  jer_param_er[6][11]=0.288507;
  jer_param[6][12]=-0.190928;  jer_param_er[6][12]=0.016249;
  jer_param[7][0]=0.0658253;  jer_param_er[7][0]=0.00272349;
  jer_param[7][1]=-0.000196113;  jer_param_er[7][1]=1.14286E-5;
  jer_param[7][2]=-1.07207;  jer_param_er[7][2]=0.137955;
  jer_param[7][3]=0.889453;  jer_param_er[7][3]=0.015712;
  jer_param[7][4]=0.00152891;  jer_param_er[7][4]=0.000140114;
  jer_param[7][5]=-0.0287057;  jer_param_er[7][5]=0.00306348;
  jer_param[7][6]=0.000128678;  jer_param_er[7][6]=1.50299E-5;
  jer_param[7][7]=-2.67266;  jer_param_er[7][7]=0.129521;
  jer_param[7][8]=0.170338;  jer_param_er[7][8]=0.00800203;
  jer_param[7][9]=0.00119121;  jer_param_er[7][9]=7.89355E-5;
  jer_param[7][10]=-8.77829;  jer_param_er[7][10]=1.94651;
  jer_param[7][11]=2.5475;  jer_param_er[7][11]=0.581462;
  jer_param[7][12]=0.0673843;  jer_param_er[7][12]=0.0379263;
  jer_param[8][0]=0.058632;  jer_param_er[8][0]=0.00246133;
  jer_param[8][1]=-8.6882E-5;  jer_param_er[8][1]=7.27552E-6;
  jer_param[8][2]=-2.81222;  jer_param_er[8][2]=0.177198;
  jer_param[8][3]=0.87612;  jer_param_er[8][3]=0.0135608;
  jer_param[8][4]=0.000862571;  jer_param_er[8][4]=8.72723E-5;
  jer_param[8][5]=-0.0272999;  jer_param_er[8][5]=0.00336864;
  jer_param[8][6]=5.7274E-5;  jer_param_er[8][6]=1.19158E-5;
  jer_param[8][7]=-3.16373;  jer_param_er[8][7]=0.187887;
  jer_param[8][8]=0.218132;  jer_param_er[8][8]=0.00804389;
  jer_param[8][9]=-0.000115495;  jer_param_er[8][9]=5.26047E-5;
  jer_param[8][10]=6.59552;  jer_param_er[8][10]=4.14997;
  jer_param[8][11]=-3.91415;  jer_param_er[8][11]=1.16124;
  jer_param[8][12]=0.689584;  jer_param_er[8][12]=0.0728198;
  jer_param[9][0]=0.047587;  jer_param_er[9][0]=0.00405247;
  jer_param[9][1]=-3.67703E-5;  jer_param_er[9][1]=1.32232E-5;
  jer_param[9][2]=-2.82013;  jer_param_er[9][2]=0.277916;
  jer_param[9][3]=0.90659;  jer_param_er[9][3]=0.0207463;
  jer_param[9][4]=0.00138775;  jer_param_er[9][4]=0.000128634;
  jer_param[9][5]=-0.0150813;  jer_param_er[9][5]=0.00500665;
  jer_param[9][6]=4.34831E-5;  jer_param_er[9][6]=1.98595E-5;
  jer_param[9][7]=-4.24538;  jer_param_er[9][7]=0.258383;
  jer_param[9][8]=0.235095;  jer_param_er[9][8]=0.0112121;
  jer_param[9][9]=-0.000176939;  jer_param_er[9][9]=7.11021E-5;
  jer_param[9][10]=17.7377;  jer_param_er[9][10]=6.10388;
  jer_param[9][11]=-6.60975;  jer_param_er[9][11]=1.55793;
  jer_param[9][12]=0.757554;  jer_param_er[9][12]=0.0890909;
  jer_param[10][0]=0.0702999;  jer_param_er[10][0]=0.00524409;
  jer_param[10][1]=-0.0001322;  jer_param_er[10][1]=1.44554E-5;
  jer_param[10][2]=-4.64174;  jer_param_er[10][2]=0.433494;
  jer_param[10][3]=1.06671;  jer_param_er[10][3]=0.0316977;
  jer_param[10][4]=0.00105544;  jer_param_er[10][4]=0.000168547;
  jer_param[10][5]=-0.02437;  jer_param_er[10][5]=0.00587413;
  jer_param[10][6]=2.26747E-5;  jer_param_er[10][6]=1.97589E-5;
  jer_param[10][7]=-3.8046;  jer_param_er[10][7]=0.382791;
  jer_param[10][8]=0.268162;  jer_param_er[10][8]=0.0148949;
  jer_param[10][9]=-6.26511E-5;  jer_param_er[10][9]=8.37528E-5;
  jer_param[10][10]=10.2852;  jer_param_er[10][10]=7.77639;
  jer_param[10][11]=-4.34279;  jer_param_er[10][11]=1.80774;
  jer_param[10][12]=0.564928;  jer_param_er[10][12]=0.0947709;
  jer_param[11][0]=0.0276176;  jer_param_er[11][0]=0.00700865;
  jer_param[11][1]=-1.08182E-5;  jer_param_er[11][1]=1.77066E-5;
  jer_param[11][2]=-4.30341;  jer_param_er[11][2]=0.642783;
  jer_param[11][3]=1.15381;  jer_param_er[11][3]=0.0466366;
  jer_param[11][4]=0.000660203;  jer_param_er[11][4]=0.000219101;
  jer_param[11][5]=-0.0365288;  jer_param_er[11][5]=0.007012;
  jer_param[11][6]=5.34815E-5;  jer_param_er[11][6]=2.21572E-5;
  jer_param[11][7]=-4.83064;  jer_param_er[11][7]=0.460712;
  jer_param[11][8]=0.28416;  jer_param_er[11][8]=0.0195007;
  jer_param[11][9]=-4.88533E-7;  jer_param_er[11][9]=0.000101491;
  jer_param[11][10]=14.821;  jer_param_er[11][10]=8.92697;
  jer_param[11][11]=-5.32396;  jer_param_er[11][11]=2.01108;
  jer_param[11][12]=0.568491;  jer_param_er[11][12]=0.10153;
  jer_param[12][0]=-0.00456382;  jer_param_er[12][0]=0.00797029;
  jer_param[12][1]=5.93327E-5;  jer_param_er[12][1]=2.2512E-5;
  jer_param[12][2]=-6.05263;  jer_param_er[12][2]=0.641518;
  jer_param[12][3]=1.01168;  jer_param_er[12][3]=0.0580327;
  jer_param[12][4]=0.00124107;  jer_param_er[12][4]=0.000266261;
  jer_param[12][5]=-0.0672398;  jer_param_er[12][5]=0.0103374;
  jer_param[12][6]=0.000144969;  jer_param_er[12][6]=3.13729E-5;
  jer_param[12][7]=-7.44458;  jer_param_er[12][7]=0.643695;
  jer_param[12][8]=0.236166;  jer_param_er[12][8]=0.0233694;
  jer_param[12][9]=0.000272798;  jer_param_er[12][9]=0.000115376;
  jer_param[12][10]=15.9149;  jer_param_er[12][10]=16.3677;
  jer_param[12][11]=-5.06773;  jer_param_er[12][11]=3.26237;
  jer_param[12][12]=0.515062;  jer_param_er[12][12]=0.14614;
  jer_param[13][0]=-0.0364159;  jer_param_er[13][0]=0.0230243;
  jer_param[13][1]=0.000103308;  jer_param_er[13][1]=6.3942E-5;
  jer_param[13][2]=-2.04414;  jer_param_er[13][2]=1.61387;
  jer_param[13][3]=1.16241;  jer_param_er[13][3]=0.154151;
  jer_param[13][4]=0.00449765;  jer_param_er[13][4]=0.000680983;
  jer_param[13][5]=-0.0713694;  jer_param_er[13][5]=0.0128;
  jer_param[13][6]=6.17634E-5;  jer_param_er[13][6]=3.50159E-5;
  jer_param[13][7]=-6.69936;  jer_param_er[13][7]=0.883717;
  jer_param[13][8]=0.401773;  jer_param_er[13][8]=0.0401819;
  jer_param[13][9]=0.00048766;  jer_param_er[13][9]=0.000171866;
  jer_param[13][10]=-6.84164;  jer_param_er[13][10]=15.4788;
  jer_param[13][11]=1.14003;  jer_param_er[13][11]=2.85757;
  jer_param[13][12]=0.063226;  jer_param_er[13][12]=0.116727;
  jer_param[14][0]=0;  jer_param_er[14][0]=0;
  jer_param[14][1]=0;  jer_param_er[14][1]=0;
  jer_param[14][2]=0;  jer_param_er[14][2]=0;
  jer_param[14][3]=0;  jer_param_er[14][3]=0;
  jer_param[14][4]=0;  jer_param_er[14][4]=0;
  jer_param[14][5]=-0.127208;  jer_param_er[14][5]=0.0191958;
  jer_param[14][6]=5.65717E-5;  jer_param_er[14][6]=3.57586E-5;
  jer_param[14][7]=-4.42422;  jer_param_er[14][7]=1.87829;
  jer_param[14][8]=0.9097;  jer_param_er[14][8]=0.0515384;
  jer_param[14][9]=-0.000149165;  jer_param_er[14][9]=0.000158924;
  jer_param[14][10]=0;  jer_param_er[14][10]=0;
  jer_param[14][11]=0;  jer_param_er[14][11]=0;
  jer_param[14][12]=0;  jer_param_er[14][12]=0;

  if(i<15 && j<13)
    {
      if(code==0) value=jer_param[i][j];
      if(code==1) value=jer_param_er[i][j];
    }
  return value;
}

//________________ returns index for eta_bin in JER
int TStnMetModel::WhatJEReta(double eta_det)
{
  int eta_bin;
  for(int i=0; i<14; i++)
    {
      if(fabs(eta_det)>=i*0.2 && fabs(eta_det)<(i+1)*0.2)
	{
	  eta_bin=i;
	  break;
	}
    }
  if(fabs(eta_det)>=2.8) eta_bin=14;
  return eta_bin;
}

//_________ JER sytematics due to Data-MC difference in JER (see e-log entry as of 01/08/07)
// E-log entry 01/08/07 demonstrates that Data-MC difference is smaller than uncertainty of data fits.
// To be conservative, I'm going to use uncertainty on data as an upper limit on Data-MC differences.
// This scale factor is to be used to scale sigma_Gauss and sigma_Landau
double TStnMetModel::JER_DataMCdiff(double jet_E, double eta_det, int syst_code) {
  double scale=1.0;
  double x=jet_E*fabs(sin(2.0*atan(exp(-eta_det))));
  if(x>0.0)
    {
      double def=sqrt(0.457/x+20.31/(x*x)+0.00834);
      double term1=1.0-sqrt(0.488/x+19.76/(x*x)+0.00814)/sqrt(0.416/x+20.9/(x*x)+0.00858);
      double term2=1.0-sqrt(0.533/x+19.3/(x*x)+0.00831)/def;
      double term3=1.0-sqrt(0.228/x+20.72/(x*x)+0.00880)/def;
      double term4=1.0-(3.942/x+0.08323)/def;
      double syst=sqrt(term1*term1+term2*term2+term3*term3+term4*term4);
      if(syst_code>0.0) scale=scale+syst;
      else scale=scale-syst;
    }
  if(scale<0.0) scale=1.0;
  return scale;
}


//_________________ returns mean of Gauss in JER: A0+A1*Pt+A2/E
double TStnMetModel::MyJer_meanG(double jet_E, double eta_det, int stat_code, int syst_code)
{
  int bin=WhatJEReta(eta_det);
  double p1=JERparam(0,bin,0);
  double p2=JERparam(0,bin,1);
  double p3=JERparam(0,bin,2);
  double val_err=0.0;
  double val=0.0;
  if(jet_E>0.0)
    {
      if(abs(syst_code)==1)
	{
	  double t1=JERparam(1,bin,0);
	  double t2=JERparam(1,bin,1)*jet_E;
	  double t3=JERparam(1,bin,2)/jet_E;
	  double t12=2.0*t1*t2*JERCorrCoeff(bin,0);
	  double t13=2.0*t1*t3*JERCorrCoeff(bin,1);
	  double t23=2.0*t2*t3*JERCorrCoeff(bin,2);
	  double t=t1*t1+t2*t2+t3*t3+t12+t13+t23;
	  if(t>0.0) val_err=(syst_code>0) ? sqrt(t) : -1.0*sqrt(t);
	}
      val=p1+p2*jet_E+p3/jet_E+val_err;
      if(val<-1.0 || val>5.0) val=0.0;
    }
  return val;
}
//_________________ returns sigma of Gauss in JER: sqrt(A3/E+A4)
double TStnMetModel::MyJer_sigmaG(double jet_E, double eta_det, int stat_code, int syst_code)
{
  int bin=WhatJEReta(eta_det);
  double p1=JERparam(0,bin,3);
  double p2=JERparam(0,bin,4);
  double val=0.0;
  double val_err=0.0;
  if(jet_E>0.0) 
    {
      val=p1/jet_E+p2;
      if(val>0.0)
	{
	  if(abs(syst_code)==2)
	    {
	      double t1=JERparam(1,bin,3)/jet_E;
	      double t2=JERparam(1,bin,4);
	      double t12=2.0*t1*t2*JERCorrCoeff(bin,3);
	      double t=0.25*(t1*t1+t2*t2+t12)/val;
	      if(t>0.0) val_err=(syst_code>0) ? sqrt(t) : -1.0*sqrt(t);
	    }
	  if(abs(syst_code)==6) val_err=(JER_DataMCdiff(jet_E,eta_det,syst_code)-1.0)*sqrt(val); // Data-MC difference
	  val=sqrt(val)+val_err;
	}
      if(val<0.0) val=0.0;
    }
  return val;
}
//_________________ returns mpv of Landau in JER: A5+A6*Pt+A7/E
double TStnMetModel::MyJer_mpvL(double jet_E, double eta_det, int stat_code, int syst_code)
{
  int bin=WhatJEReta(eta_det);
  double p1=JERparam(0,bin,5);
  double p2=JERparam(0,bin,6);
  double p3=JERparam(0,bin,7);
  double val_err=0.0;
  double val=0.0;
  if(jet_E>0.0)
    {
      if(abs(syst_code)==3)
	{
	  double t1=JERparam(1,bin,5);
	  double t2=JERparam(1,bin,6)*jet_E;
	  double t3=JERparam(1,bin,7)/jet_E;
	  double t12=2.0*t1*t2*JERCorrCoeff(bin,4);
	  double t13=2.0*t1*t3*JERCorrCoeff(bin,5);
	  double t23=2.0*t2*t3*JERCorrCoeff(bin,6);
	  double t=t1*t1+t2*t2+t3*t3+t12+t13+t23;
	  if(t>0.0) val_err=(syst_code>0) ? sqrt(t) : -1.0*sqrt(t);
	}
      val=p1+p2*jet_E+p3/jet_E+val_err;
      if(val<-1.0 || val>5.0) val=0.0;
    }
  return val;
}
//_________________ returns sigma of Landau in JER: sqrt(A8/E+A9)
double TStnMetModel::MyJer_sigmaL(double jet_E, double eta_det, int stat_code, int syst_code)
{
  int bin=WhatJEReta(eta_det);
  double p1=JERparam(0,bin,8);
  double p2=JERparam(0,bin,9);
  double val=0.0;
  double val_err=0.0;
  if(jet_E>0.0) 
    {
      val=p1/jet_E+p2;
      if(val>0.0)
	{
	  if(abs(syst_code)==4)
	    {
	      double t1=JERparam(1,bin,8)/jet_E;
	      double t2=JERparam(1,bin,9);
	      double t12=2.0*t1*t2*JERCorrCoeff(bin,7);
	      double t=0.25*(t1*t1+t2*t2+t12)/val;
	      if(t>0.0) val_err=(syst_code>0) ? sqrt(t) : -1.0*sqrt(t);
	    }
	  if(abs(syst_code)==6) val_err=(JER_DataMCdiff(jet_E,eta_det,syst_code)-1.0)*sqrt(val); // Data-MC difference
	  val=sqrt(val)+val_err;
	}
      if(val<0.0) val=0.0;
    }
  return val;
}

//_________________ returns normalization of Gauss: [A10+A11*sqrt(E)]/E+A12
double TStnMetModel::MyJer_normG(double jet_E, double eta_det, int stat_code, int syst_code)
{
  int bin=WhatJEReta(eta_det);
  double p1=JERparam(0,bin,10);
  double p2=JERparam(0,bin,11);
  double p3=JERparam(0,bin,12);
  double val=0.0;
  double val_err=0.0;
  if(jet_E>0.0) 
    {
      if(abs(syst_code)==5)
	{
	  double t1=JERparam(1,bin,10)/jet_E;
	  double t2=JERparam(1,bin,11)/sqrt(jet_E);
	  double t3=JERparam(1,bin,12);
	  double t12=2.0*t1*t2*JERCorrCoeff(bin,8);
	  double t13=2.0*t1*t3*JERCorrCoeff(bin,9);
	  double t23=2.0*t2*t3*JERCorrCoeff(bin,10);
	  double t=t1*t1+t2*t2+t3*t3+t12+t13+t23;
	  if(t>0.0) val_err=(syst_code>0) ? sqrt(t) : -1.0*sqrt(t);
	}
      val=(p1+p2*sqrt(jet_E))/jet_E+p3+val_err;
    }
  if(val<0.0) val=0.0;
  return val;
}


//___ This is new function to generate Met; it takes care of unclustered & jet components of Met.
void TStnMetModel::GenerateMyTotalMet(JetStuff jetstuff, CommonStuff miscstuff, int systcode, int systcodeUncl, int rnd_seed, TVector2 &myGenMet) {
  if(systcode!=0 || systcodeUncl!=0) gRandom->SetSeed(rnd_seed);
  myGenMet.Set(0.0,0.0);
  //___________ generating met due to jets
  double SumEtJet_smear=0.0; // U.E.+M.I. for jets above threshold after smearing
  double SumEtJet=0.0;       // U.E.+M.I. for jets above threshold before smearing 
  double SumEtRawJet=0.0;
  double SumEtRawJet_smear=0.0;
  double Pt_cut=15.0;
  double scale_factor=1.0;
  TLorentzVector jetmetsum(0.0,0.0,0.0,0.0);
  TVector2 myJetMet(0.0,0.0);
  int NjetCut=jetstuff.myNjet_th15;   
  if(NjetCut>0) 
    {
      for(int i=0; i<jetstuff.newJetLev6.size(); i++) // "jetstuff.Jet_lev6_noEMobj" is replaced by "jetstuff.newJetLev6" (11/02/06)
	{
	  //on line below I use exactly the same eta cut as in MET correction
	  if(jetstuff.newJetLev6[i].Pt()>3.0 && fabs(jetstuff.newJetLev6[i].Eta())<=MaxJetEta) // "jetstuff.Jet_lev6_noEMobj" is replaced by "jetstuff.newJetLev6" (11/02/06) 
	    {
	      double parJER[5];
	      parJER[0]=MyJer_meanG(jetstuff.newJetLev6[i].E(),jetstuff.EtaDetCorr[i],0,systcode); // for now systcode is used as jer_stat_code
	      parJER[1]=MyJer_sigmaG(jetstuff.newJetLev6[i].E(),jetstuff.EtaDetCorr[i],0,systcode);
	      parJER[2]=MyJer_mpvL(jetstuff.newJetLev6[i].E(),jetstuff.EtaDetCorr[i],0,systcode);
	      parJER[3]=MyJer_sigmaL(jetstuff.newJetLev6[i].E(),jetstuff.EtaDetCorr[i],0,systcode);
	      parJER[4]=MyJer_normG(jetstuff.newJetLev6[i].E(),jetstuff.EtaDetCorr[i],0,systcode);
	      double jer_limit=MyIntegralUpLimit(jetstuff.newJetLev6[i].E(),jetstuff.EtaDetCorr[i],0,systcode);
	      TF1* JerFun=new TF1("jer",MyJER,-1.0,jer_limit,5);
	      //------ setting Gaussian
	      JerFun->SetParameter(0,parJER[0]); 
	      JerFun->SetParameter(1,parJER[1]); 
	      //------ setting Landau
	      JerFun->SetParameter(2,parJER[2]);
	      JerFun->SetParameter(3,parJER[3]);
	      //------ setting Gaussian normalization
	      JerFun->SetParameter(4,parJER[4]);
	      double off_set=JerFun->GetMaximumX(-1.0,1.0); // returns mpv of JER function;
	      scale_factor=1.0-off_set+JerFun->GetRandom(); // get jet energy scale factor; default=assume jet scale is 1.0
	      if(scale_factor<0.0) scale_factor=0.0; 
	      delete JerFun;
	      if((scale_factor*jetstuff.newJetLev6[i]).Pt()>Pt_cut) // "jetstuff.Jet_lev6_noEMobj" is replaced by "jetstuff.newJetLev6" (11/02/06)
		{
		  jetmetsum=jetmetsum+(1.0-scale_factor)*jetstuff.newJetLev6[i]; // "jetstuff.Jet_lev6_noEMobj" is replaced by "jetstuff.newJetLev6" (11/02/06)
		  // SumEtJet_smear= Sum(U.E.+M.I.)
		  SumEtJet_smear=SumEtJet_smear+jetstuff.Jet_lev5_noEMobj[i].Pt()
		    +jetstuff.Jet_lev1_noEMobj[i].Pt()
		    -jetstuff.Jet_lev4_noEMobj[i].Pt()
		    -jetstuff.Jet_lev6_noEMobj[i].Pt(); // Jet_lev6_noEMobj added on 06/26/06
		  SumEtRawJet_smear=SumEtRawJet_smear+jetstuff.Jet_raw_noEMobj[i].Pt();
		}
	      if(jetstuff.Jet_lev6_noEMobj[i].Pt()>Pt_cut) 
		{
		  SumEtJet=SumEtJet+jetstuff.Jet_lev5_noEMobj[i].Pt()
		    +jetstuff.Jet_lev1_noEMobj[i].Pt()
		    -jetstuff.Jet_lev4_noEMobj[i].Pt()
		    -jetstuff.Jet_lev6_noEMobj[i].Pt(); // Jet_lev6_noEMobj added on 06/26/06
		  SumEtRawJet=SumEtRawJet+jetstuff.Jet_raw_noEMobj[i].Pt();
		}
	    }
	}
      myJetMet.Set(jetmetsum.Px(),jetmetsum.Py());   
    }
  //_____________ generating unclustered component of Met
  double x;
  double y;
  double _meanX=0.0;
  double _sigmaX=0.0;
  double _meanY=0.0;
  double _sigmaY=0.0;
  double _normX=0.0;
  double _scaleX=0.0;
  double _normY=0.0;
  double _scaleY=0.0;
  double sumEt=0.0;
  double sumEt_raw=0.0;  
  //--------------------- this chunk of code is modified on 11/09/06
  sumEt_raw=jetstuff.mySumEtCorr_th15-SumEtJet+SumEtRawJet;
  sumEt=sumEt_raw+SumEtJet_smear-SumEtRawJet_smear;
  //------------------------------------------------ end of 11/09/06 part
  if(sumEt<0.0) sumEt=0.0;
  GetMyUnclMetResolution(sumEt,fJTC_imode,fUnclParamSwitch,systcodeUncl,_sigmaX,_sigmaY,_meanX,_meanY,_normX,_normY,_scaleX,_scaleY);
  //------- new part (05/24/07)
  TF1* UnclMetFun=new TF1("fit",UnclMetResolution,-200.0,200.0,4);
  //------ setting 1st Gaussian
  UnclMetFun->SetParameter(0,_meanX);
  UnclMetFun->SetParameter(1,_sigmaX);  
  //------ setting 2nd Gaussian
  UnclMetFun->SetParameter(2,_normX);  
  UnclMetFun->SetParameter(3,_scaleX);
  x=UnclMetFun->GetRandom();
  UnclMetFun->SetParameter(0,_meanY);
  UnclMetFun->SetParameter(1,_sigmaY);  
  //------ setting 2nd Gaussian
  UnclMetFun->SetParameter(2,_normY);  
  UnclMetFun->SetParameter(3,_scaleY);
  y=UnclMetFun->GetRandom();
  delete UnclMetFun;
  //------- end of new part (05/24/07)
  myGenMet.Set(x,y); 
  myGenMet=myGenMet+myJetMet;
  return;
}



//--------------------------------------------------------------------------------
// >>>>>>>>>>>>>>  Accessors to output params <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
//________________________________________________________________________________

// returns default value of metsig for data event as well as default value for MC event
double TStnMetModel::GetMetsigData() {
  return metsigstuff.metsig_Def;
}

// returns deviated values of metsig for sytematics studies in MC; syst_ind<2*N+1+Nspare
//    2*N-- stands for N=10 sources of syst. (+/-) due to parameterization of energy resolution
//    1 -- differnt schemes of parameterization of unclustered energy
//    Nspare=4 -- reserved for "future" syst. uncertainties we can potentially come up with
double TStnMetModel::GetMetsigSyst(int syst_ind) {
  if(syst_ind>-1 && syst_ind<MetsigArray) return metsigstuff.metsig_syst[syst_ind];
  else metsigstuff.metsig_Def; 
}

// returns corrected MET as it's calculated by Met Model
TVector2 TStnMetModel::GetMetData() {
  return jet04stuff.myMETcorr_th15;
}

// number of pseudo-experiments for MetModel
int TStnMetModel::GetNpseudo() { 
  return Npseudo; 
} 


// returns default MET-significance of generated MET from pseudo-experiments
double TStnMetModel::GetMetsigPseudo(int indPseudo) {
  double sig=0.0;
  if(indPseudo>-1 
     && indPseudo<Npseudo 
     && indPseudo<metsigstuff.metsig_DefPseudo.size()) sig=metsigstuff.metsig_DefPseudo[indPseudo];
  return sig;
}

// returns systematics on MET-significance for generated MET from pseudo-experiments 
double TStnMetModel::GetMetsigPseudoSyst(int indPseudo, int syst_ind) {
  double sig=0.0;
  if(indPseudo>-1 && indPseudo<Npseudo && syst_ind>-1 && syst_ind<MetsigArray
     && indPseudo<metsigstuff.metsig_systPseudo[syst_ind].size()) sig=metsigstuff.metsig_systPseudo[syst_ind][indPseudo];
  return sig;
}

//_____ returns MET predicted by Met Model according to default parameterization
TVector2 TStnMetModel::GetMetPseudo(int indPseudo) {
  TVector2 dummyVec(0.0,0.0);
  if(indPseudo>-1 
     && indPseudo<Npseudo 
     && indPseudo<metsigstuff.MetPseudo.size()) dummyVec=metsigstuff.MetPseudo[indPseudo]; 
  return dummyVec; 
}

//_____ returns deviated MET predicted by Met Model for systematics studies
TVector2 TStnMetModel::GetMetPseudoSyst(int indPseudo, int syst_ind) {
  TVector2 dummyVec(0.0,0.0);
  if(indPseudo>-1 && indPseudo<Npseudo && syst_ind>-1 && syst_ind<MetsigArray
     && indPseudo<metsigstuff.MetPseudoSyst[syst_ind].size()) dummyVec=metsigstuff.MetPseudoSyst[syst_ind][indPseudo];
  return dummyVec; 
}

// _____flag: 0=no bad met; 1=bad met. Output of my MET cleanup cuts tuned to remove events 
// where MET is from badly mis-measured electrons or from Z+pho events with lost photon. You may find it useful.
int TStnMetModel::IsBadMet() {
  int passcode=0;

  double MetEM_dPhi_max=0.4;
  double MetPho_dPhi_max=0.3;
  for(int i=0; i<allstuff.myCorrPhoton.size(); i++) // don't look at 1st photon for twr-9 effect
    {
      double dPhi=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(allstuff.myCorrPhoton[i].Phi())-TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())));
      if(allstuff.myCorrPhoton.size()>1)
	{
	  if(i>0 && (fabs(allstuff.myPhoEtaDet[i])>1.0 || fabs(allstuff.myPhoEtaDet[i])<0.1) && dPhi<MetPho_dPhi_max) passcode++;
	}
      else
	{
	  if(fabs((allstuff.myPhoEtaDet[i])>1.0 || fabs(allstuff.myPhoEtaDet[i])<0.1) && dPhi<MetPho_dPhi_max) passcode++;
	}
    }      
  for(int i=0; i<allstuff.myCorrElectron.size(); i++)
    {
      double dPhi1=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(allstuff.myCorrElectron[i].Phi())-TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())));
      double dPhi2=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(allstuff.myElePhiDet[i])-TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())));
      double dPhi=(dPhi1<dPhi2) ? dPhi1 : dPhi2;
      if(allstuff.myCorrElectron.size()>1)
	{
	  if(i>0 && fabs(fabs(allstuff.myEleEtaDet[i])-1.05)<0.05 && dPhi<MetEM_dPhi_max) passcode++; // new as of 06/25/08
	}
      else
	{
	  if(fabs(fabs(allstuff.myEleEtaDet[i])-1.05)<0.05 && dPhi<MetEM_dPhi_max) passcode++; // new as of 06/25/08	      
	}
    }      
  
  //---------------------- This part was added on 07/03/08, to remove events where a photon got lost in CEM cracks  
  double jet_minPt=5.0;
  double jet_maxEta=1.3;
  double jet_maxEmFr1=0.875;
  double jet_maxEmFr2=0.3;
  double MetJet_dPhi_max=0.3;
  int jet_NtwrMax=10;
  int jet_NtrkMax=5;
  
  double dPhi_min=TMath::Pi();
  double eta_det=10.0;
  double emfr=0.5;
  int ntwr=100;
  int ntrk=100;
  for(int i=0; i<jet04stuff.Jet_lev6_noEMobj.size(); i++)
    {
      if(jet04stuff.Jet_lev6_noEMobj[i].Pt()>jet_minPt && 
	 fabs(jet04stuff.EtaDet[i])<1.3 && 
	 jet04stuff.Nobj_match[i]==0)
	{
	  double dPhi=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(jet04stuff.Jet_lev6_noEMobj[i].Phi())-TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())));
	  if(dPhi<dPhi_min) 
	    {
	      dPhi_min=dPhi;
	      eta_det=jet04stuff.EtaDet[i];
	      emfr=jet04stuff.EmFrRaw[i];
	      ntwr=jet04stuff.JetNtwr[i];
	      ntrk=jet04stuff.JetNtrk[i];
	    }
	}
    }
  if(dPhi_min<MetJet_dPhi_max && emfr>jet_maxEmFr1 && ntwr<jet_NtwrMax && ntrk<jet_NtrkMax && fabs(eta_det)<1.1) passcode++;
  if(dPhi_min<MetJet_dPhi_max && emfr<jet_maxEmFr2 && ntwr<jet_NtwrMax && ntrk<jet_NtrkMax 
     && (fabs(eta_det)<0.1 || (fabs(eta_det)>1.1 && fabs(eta_det)<1.2))) passcode++;
  
  if(passcode>0) passcode=1;   
  return passcode;
}

//_____________ MET*dPhiClosest/Et(closest jet), dPhiClosest=min(dPhi=Phi(met)-Phi(obj))
double TStnMetModel::GetMetJetRatio() {
  double ratio=-1.0;
  double jetEt=0.0;
  double dphi2=CalculateMetDelPhi(2,jetEt);
  if(jetEt>0.0) ratio=cos(dphi2)*jet04stuff.myMETcorr_th15.Mod()/jetEt;
  return ratio;
}

//_____________ Et(closest jet to MET), dPhiClosest=min(dPhi=Phi(met)-Phi(obj))
double TStnMetModel::GetEtClosestJet() {
  double jetEt=0.0;
  double dphi2=CalculateMetDelPhi(2,jetEt);
  return jetEt;
}


//_________ min(dPhi=Phi(met)-Phi(jet))
double TStnMetModel::GetDelPhiClosestJet() {
  return CalculateMetDelPhi(2);
}

//_________ min(dPhi=Phi(met)-Phi(obj))
double TStnMetModel::GetDelPhiClosest() {
  return CalculateMetDelPhi(0);
}

//____ returns default value of "signed" metsig for data event
double TStnMetModel::GetMetsigSignedData() {
  double dphi7=CalculateMetDelPhi(7);
  double _sign=1.0;
  if(dphi7<0.0) _sign=-1.0; 
  return _sign*metsigstuff.metsig_Def;
}

//____ returns deviated values of "signed" metsig for sytematics studies in MC; syst_ind<2*N+1+Nspare
double TStnMetModel::GetMetsigSignedSyst(int syst_ind) {
  double dphi7=CalculateMetDelPhi(7);
  double _sign=1.0;
  if(dphi7<0.0) _sign=-1.0;
  return _sign*GetMetsigSyst(syst_ind);
}


// i=0- closest obj; i=1- closest EM obj; i=2- closest jet; 
// i=3- 1st EM; i=4- 1st jet; i=5- dPhi(jj); i=6- dPhi(jj)+ min(dPhi[met-j])
// i=7 signed dPhi to closest jet
//_____________________________ calculates dPhi for metsig.vs.dPhi studies
double TStnMetModel::CalculateMetDelPhi(int ind) {

  double dphi_min=10.0;
  double dPhi_det[8];
  for(int ij=0; ij<8; ij++)
    {
      dPhi_det[ij]=10.0; // initial dummy value
    }

  for(int ij=0; ij<allstuff.myCorrPhoton.size(); ij++)
    {
      double _dphi_d=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(allstuff.myCorrPhoton[ij].Phi())-TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())));
      if(ij==0) dPhi_det[3]=_dphi_d;
      if(_dphi_d<dPhi_det[1]) dPhi_det[1]=_dphi_d;
    }
  for(int ij=0; ij<allstuff.myCorrElectron.size(); ij++)
    {
      double _dphi_d=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(allstuff.myCorrElectron[ij].Phi())-TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())));
      if(ij==0) 
	{
	  // make sure there are photons before next step	  
	  if(allstuff.myCorrPhoton.size()>0) 
	    {
	      if(allstuff.myCorrElectron[0].Pt()>allstuff.myCorrPhoton[0].Pt()) dPhi_det[3]=_dphi_d;
	    }
	  else dPhi_det[3]=_dphi_d;
	}
      if(_dphi_d<dPhi_det[1]) dPhi_det[1]=_dphi_d;
    }
  double max_jet_Pt=0.0;
  for(int ij=0; ij<jet04stuff.Jet_lev6_noEMobj.size(); ij++)
    {
      if(jet04stuff.Jet_lev6_noEMobj[ij].Pt()>3.0 && (jet04stuff.Npho_match[ij]+jet04stuff.Nele_match[ij])==0)
	{
	  double _dphi_d=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(jet04stuff.Jet_lev6_noEMobj[ij].Phi())-TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())));
	  if(jet04stuff.Jet_lev6_noEMobj[ij].Pt()>max_jet_Pt) 
	    {
	      max_jet_Pt=jet04stuff.Jet_lev6_noEMobj[ij].Pt();
	      dPhi_det[4]=_dphi_d;
	    }
	  if(_dphi_d<dPhi_det[2]) dPhi_det[2]=_dphi_d;
	}
    }
  if(dPhi_det[2]<dPhi_det[1]) dPhi_det[0]=dPhi_det[2];
  else dPhi_det[0]=dPhi_det[1];
  if(jet04stuff.Jet_lev6_noEMobj.size()>=2)
    {
      dPhi_det[5]=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(jet04stuff.Jet_lev6_noEMobj[0].Phi())-TVector2::Phi_0_2pi(jet04stuff.Jet_lev6_noEMobj[1].Phi())));
      dPhi_det[6]=dPhi_det[5]+dPhi_det[2];
      double phi1=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(jet04stuff.Jet_lev6_noEMobj[0].Phi())-TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())));
      double phi2=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(jet04stuff.Jet_lev6_noEMobj[1].Phi())-TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())));
      double phi12=phi1+phi2;
      if(phi12<TMath::Pi()) dPhi_det[7]=-1.0*dPhi_det[2];
      else dPhi_det[7]=dPhi_det[2];
    }
  if(ind>-1 && ind<8) dphi_min=dPhi_det[ind];
  return dphi_min;
}

// i=0- closest obj; i=1- closest EM obj; i=2- closest jet; 
// i=3- 1st EM; i=4- 1st jet; i=5- dPhi(jj); i=6- dPhi(jj)+ min(dPhi[met-j])
// i=7 signed dPhi to closest jet
//_____________________________ calculates dPhi and returns Et(closest object), for metsig.vs.dPhi studies
double TStnMetModel::CalculateMetDelPhi(int ind, double &Et_closest) {

  Et_closest=0.0;
  double dphi_min=10.0;
  double dPhi_det[8];
  double Et_ClosestObj[8];
  for(int ij=0; ij<8; ij++)
    {
      dPhi_det[ij]=10.0; // initial dummy value
      Et_ClosestObj[ij]=0.0;
    }

  for(int ij=0; ij<allstuff.myCorrPhoton.size(); ij++)
    {
      double _dphi_d=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(allstuff.myCorrPhoton[ij].Phi())-TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())));
      if(ij==0) 
	{
	  dPhi_det[3]=_dphi_d;
	  Et_ClosestObj[3]=allstuff.myCorrPhoton[ij].Pt();
	}
      if(_dphi_d<dPhi_det[1]) 
	{
	  dPhi_det[1]=_dphi_d;
	  Et_ClosestObj[1]=allstuff.myCorrPhoton[ij].Pt();
	}
    }
  for(int ij=0; ij<allstuff.myCorrElectron.size(); ij++)
    {
      double _dphi_d=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(allstuff.myCorrElectron[ij].Phi())-TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())));
      if(ij==0) 
	{
	  // make sure there are photons before next step	  
	  if(allstuff.myCorrPhoton.size()>0) 
	    {
	      if(allstuff.myCorrElectron[0].Pt()>allstuff.myCorrPhoton[0].Pt()) 
		{
		  dPhi_det[3]=_dphi_d;
		  Et_ClosestObj[3]=allstuff.myCorrElectron[ij].Pt();
		}
	    }
	  else 
	    {
	      dPhi_det[3]=_dphi_d;
	      Et_ClosestObj[3]=allstuff.myCorrElectron[ij].Pt();
	    }
	}
      if(_dphi_d<dPhi_det[1]) 
	{
	  dPhi_det[1]=_dphi_d;
	  Et_ClosestObj[1]=allstuff.myCorrElectron[ij].Pt();
	}
    }
  double max_jet_Pt=0.0;
  for(int ij=0; ij<jet04stuff.Jet_lev6_noEMobj.size(); ij++)
    {
      if(jet04stuff.Jet_lev6_noEMobj[ij].Pt()>3.0 && (jet04stuff.Npho_match[ij]+jet04stuff.Nele_match[ij])==0)
	{
	  double _dphi_d=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(jet04stuff.Jet_lev6_noEMobj[ij].Phi())-TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())));
	  if(jet04stuff.Jet_lev6_noEMobj[ij].Pt()>max_jet_Pt) 
	    {
	      max_jet_Pt=jet04stuff.Jet_lev6_noEMobj[ij].Pt();
	      dPhi_det[4]=_dphi_d;
	      Et_ClosestObj[4]=jet04stuff.Jet_lev6_noEMobj[ij].Pt();
	    }
	  if(_dphi_d<dPhi_det[2]) 
	    {
	      dPhi_det[2]=_dphi_d;
	      Et_ClosestObj[2]=jet04stuff.Jet_lev6_noEMobj[ij].Pt();
	    }
	}
    }
  if(dPhi_det[2]<dPhi_det[1]) 
    {
      dPhi_det[0]=dPhi_det[2];
      Et_ClosestObj[0]=Et_ClosestObj[2];
    }
  else 
    {
      dPhi_det[0]=dPhi_det[1];
      Et_ClosestObj[0]=Et_ClosestObj[1];
    }
  if(jet04stuff.Jet_lev6_noEMobj.size()>=2)
    {
      dPhi_det[5]=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(jet04stuff.Jet_lev6_noEMobj[0].Phi())-TVector2::Phi_0_2pi(jet04stuff.Jet_lev6_noEMobj[1].Phi())));
      dPhi_det[6]=dPhi_det[5]+dPhi_det[2];
      double phi1=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(jet04stuff.Jet_lev6_noEMobj[0].Phi())-TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())));
      double phi2=fabs(TVector2::Phi_mpi_pi(TVector2::Phi_0_2pi(jet04stuff.Jet_lev6_noEMobj[1].Phi())-TVector2::Phi_0_2pi(jet04stuff.myMETcorr_th15.Phi())));
      double phi12=phi1+phi2;
      if(phi12<TMath::Pi()) dPhi_det[7]=-1.0*dPhi_det[2];
      else dPhi_det[7]=dPhi_det[2];
    }
  if(ind>-1 && ind<8) 
    {
      dphi_min=dPhi_det[ind];
      Et_closest=Et_ClosestObj[ind];
    }
  return dphi_min;
}