/*
>>>>>>>>  Met Resolutio Model
>>>>>>>> Author: Sasha Pronko (pronko@fnal.gov)
 Met Model was designed to separate events with fake and true MET. 
 It avaluates a probability that a measured value of MET arises purely 
 from energy fluctuations in Calorimeter.
 Met Model can also provide a 2-dimensional PDF of fake MET for a particular event.   
 */
#ifndef MetModel_hh
#define MetModel_hh

#include <TH1.h>
#include <TF1.h>
#include <TH2.h>
#include <TProfile.h>
#include <TMath.h>
#include <TLorentzVector.h>
#include <TVector3.h>
#include <vector>

#include <Stntuple/loop/TStnModule.hh>
#include <Stntuple/alg/TStntuple.hh>
#include <Stntuple/photon/TPhotonUtil.hh>
#include "JetUser/JetEnergyCorrections.hh"
//------declaration of blocks to be read
#include <Stntuple/obj/TStnEvent.hh>
#include <Stntuple/obj/TCalDataBlock.hh>
#include <Stntuple/obj/TStnTrackBlock.hh>
#include <Stntuple/obj/TStnPhotonBlock.hh>
#include <Stntuple/obj/TStnElectronBlock.hh>
#include <Stntuple/obj/TStnHeaderBlock.hh>
#include <Stntuple/obj/TStnVertexBlock.hh>
#include <Stntuple/obj/TStnJetBlock.hh>
#include <Stntuple/obj/TStnMetBlock.hh>
#include <CalorGeometry/CalConstants.hh>

class TStnMetModel: public TObject {
protected:

  enum { MetsigArray = 15 };

  //-------- container for met-significance stuff
  struct MetsigStuff {
    double metsig_Def; // default metsig
    double metsig_syst[MetsigArray]; // for metsig systematics

    std::vector<double> metsig_DefPseudo; // default metsig for pseudo-experiments
    std::vector<double> metsig_systPseudo[MetsigArray]; //  metsig systematics for pseudo-experiments
    std::vector<TVector2> MetPseudo; // MET from pseudo-experiments
    std::vector<TVector2> MetPseudoSyst[MetsigArray]; // systematics for MET from pseudo-experiments
  };

  //------------------- parameters for Jet Correction initialization
  struct CorInit {
    int level;
    int nvx;
    int cone;
    int version;
    int sys;
    int Nrun;
    int imode;
  };

  struct JetStuff {
    
    TVector2 myMETcorr_th15; // my MEt(4) after corrections for jets (muons), using jets with Et>15 GeV
    double mySumEtCorr_th15; // SumEt after corrections for photons, electrons, JetClu-0.4 jets (muons), using jets with Et>15 GeV  
    double mySumEtJet_th15;  // SumEt of JetClu-0.4 jets (lev6), using jets with Et>15 GeV   
    //_____________________________________________________________
    //----- my new output params (old one's to be removed later)
    
    int myNjet; // number of clusters ("raw jets")
    int myNjet_th15; // number of level-6 jets with Et>15, after matching with EM object
    
    std::vector<TLorentzVector> Jet_raw;      // uncorrected jets 
    std::vector<TLorentzVector> Jet_lev1;     // jets corrected to level-1
    std::vector<TLorentzVector> Jet_lev4;     // jets corrected to level-4
    std::vector<TLorentzVector> Jet_lev5;     // jets corrected to level-5
    std::vector<TLorentzVector> Jet_lev6;     // jets corrected to level-6
    std::vector<TLorentzVector> Jet_lev7;     // jets corrected to level-6
    std::vector<TLorentzVector> Jet_raw_noEMobj;      // uncorrected jets, matched EM object removed 
    std::vector<TLorentzVector> Jet_lev1_noEMobj;     // jets corrected to level-1, matched EM object removed
    std::vector<TLorentzVector> Jet_lev4_noEMobj;     // jets corrected to level-4, matched EM object removed
    std::vector<TLorentzVector> Jet_lev5_noEMobj;     // jets corrected to level-5, matched EM object removed
    std::vector<TLorentzVector> Jet_lev6_noEMobj;     // jets corrected to level-6, matched EM object removed
    std::vector<TLorentzVector> Jet_lev7_noEMobj;     // jets corrected to level-7, matched EM object removed
    std::vector<TLorentzVector> newJetLev6;           // pseudo-jets for MetModel: noEM_lev6+MET*cos(dPhi)
    std::vector<int> JetNtrk;  // number of tracks in jet (Ntrk as defined in Stntuple)
    std::vector<int> JetNtwr;  // number of towers in jet (Ntwr as defined in Stntuple)
    std::vector<double> EtaDet;   // jet detector eta (before removing matched EM objects)
    std::vector<double> EtaDetCorr;   // jet detector eta (after removing matched EM objects)
    std::vector<double> EmFrRaw;  // EM fraction (before removing matched EM objects)
    std::vector<int> Nobj_match; // total number of matched objects (pho,ele,mu,tau,btag)
    std::vector<int> Npho_match; // number of matched photons
    std::vector<int> Nele_match; // number of matched electrons
    std::vector<int> Nmu_match; // number of matched muons
    std::vector<int> JetBlockInd; // original jet index in JetBlock, need this after jet reordering    
  };


  struct CommonStuff { // params which do not depend on Jet Cone size
    std::vector<TLorentzVector> myRawPhoton;    // raw photons
    std::vector<TLorentzVector> myCorrPhoton;   // corrected photons (leakage, em scale,etc.)
    std::vector<int> myPhoInd; // original photon index
    std::vector<double> myPhoEtaDet;    // detector eta for photons
    std::vector<TLorentzVector> myRawElectron;  // raw electrons
    std::vector<TLorentzVector> myCorrElectron; // corrected electrons (leakage, em scale,etc.)
    std::vector<int> myEleInd; // original electron index
    std::vector<double> myEleEtaDet;    // detector eta for electrons
    std::vector<double> myElePhiDet;    // detector eta for electrons
    std::vector<double> myElePprEt;    // Ppr Et for electrons
    std::vector<TLorentzVector> myRawMuon;      // raw muons
    std::vector<TLorentzVector> myCorrMuon;     // corrected muons (corrections???)
    std::vector<TLorentzVector> myCalMuon;     // CAL energy cluster for muons
    std::vector<double> myMuonZvx;     // zvx of muons and tracks
    std::vector<int> myMuonAlgo;    // Track Algorithm ID for muons and tracks
    std::vector<int> myMuonTrkInd; // muon ind in TStnTrack block
    std::vector<int> myTrkStereo;  // number of COT stereo SL
    std::vector<int> myTrkAxial;  // number of COT axial SL
    TVector2 myMET_raw; // MEt(4) before any corrections
    TVector2 myMET0_raw; // MEt(0) before any corrections
    double mySumEt_raw; // SumEt before any corrections    
  };

  //_______________ data blocks
  TStnHeaderBlock*   fHeaderBlock;
  TCalDataBlock*     fCalData;					
  TStnVertexBlock*   fZVertexBlock;
  TStnJetBlock*      fJetBlockClu04;
  TStnPhotonBlock*   fPhotonBlock;
  TStnElectronBlock* fElectronBlock;
  TStnMetBlock*      fMetBlock;
  TStnTrackBlock*    fTrackBlock;
  typedef std::vector<TCalTower*> CalDataArray; // need for EM-jet matching
  typedef std::vector<TCalTower*>::iterator CalDataArrayI; // need for EM-jet matching


  //----- declaring container for pho,ele,muon,jets, etc.
  CommonStuff allstuff;
  JetStuff jet04stuff;
  MetsigStuff metsigstuff;

  //_____________________________________________ External input params;
  int input_Nele;
  std::vector<int> input_EleInd;
  std::vector<TLorentzVector> input_CorrEle; 
  int input_Nmuon;
  std::vector<TLorentzVector> input_CorrMuon;
  std::vector<TLorentzVector> input_RawMuon;
  std::vector<int> input_TrkAlgoID;
  std::vector<int> input_trkInd;
  std::vector<double> input_TrkZvx;
  std::vector<TLorentzVector> input_CalCluster; 
  

  //_____________________________________________________________________ cuts & service params
  double Zvx_best;  // global variable
  int myNvx_class12; // global variable
  int myJetRun; // global variable
  int DebugMode;  // debugging mode
  int fSampleID;  // sample ID for Metsig calibation
  int fUnclParamSwitch; // switch for Unclustered parameterization: 0==di-pho; 1==Zee
  int Npseudo; // number of pseudo-experiments for MetModel
  //__________________________________________________ parameters for Jet Corrections
  int fJTC_coneSize; 
  int fJTC_version;  // for 5.3.1 and later it's recommended to use version 5  
  int fJTC_level;    // Level of corrections: 
                       //       Level=1 Relative Energy Corrections 
                       //       Level=2 Relative Energy Corrections and Time Dependence 
                       //       Level=3 Raw Energy Scale corrections, Relative Energy 
                       //               Corrections and Time Dependence 
                       //       Level=4 Multiple Interaction corrections, and all of above 
                       //       Level=5 Absolute Energy Scale and all of above 
                       //       Level=6 Underlying Event corrections and all of above. 
                       //       Level=7 Out-of-cone and all of above.
  int fJTC_systcode; // code for systematics in corrections:
                       //     0--default corrections                       
                       //     1--Relative Correction
                       //     2--Central Cal Stability
                       //     3--Scale Correction
                       //     4--Multiple Interaction Correction
                       //     5--Absolute Correction
                       //     6--Underlying Event Correction
                       //     7--Out-Of-Cone Correction
                       //     8--Splash Out Correction
                       //     100--Total Uncertainty
                       //--- if systcode>0 -- plus sigma deviation(?)
                       //--- if systcode<0 -- minos sigma deviation(?)
                       //-------- All these codes have to be double-checked!!!!
  int fJTC_imode;      // correction mode in V5* data and MC: 0=MC, 1=data


  //---- Photon Cuts
  double fPhoMinEt; // min Et of photons
  double fCemPhoMaxHADEM[3];        // HADEM cut using 3towers : [0]-value of fixed cut 
                              // **** [1]-const term in sliding cut; [2]-linear term in sliding cut
  double fCemPhoMaxCesChi2;         // CES Chi^2 (wire+strip)
  int fCemPhoMaxN3D;                // cut on number of 3D tracks in the cluster
  double fCemPhoMaxPtmax[2];     // cut on track Pt pointing to the cluster: [0]-const term; [1]-linear term
  double fCemPhoMaxCalIso4_corr[3];   // cut on fully corr. CalIso4: [0]-fixed value
                                // **** [1]-const term in sliding cut, [2]-linear term 
  double fCemPhoMaxTrkIso4[2];        // cut on TrkIso4: [0]-const term, [1]-linear term 
  double fCemPhoMax2ndCesEt[3];        // **** cut on energy of 2nd CES cluster: [0]-fixed value; [1]-const term; [2]-linear term
  double fCemPhoMinCesX;           // cut on CES X
  double fCemPhoMaxCesX; 
  double fCemPhoMinCesZ;           // cut on CES Z
  double fCemPhoMaxCesZ;           

  //---- cuts for CEM electrons
  double fEleMinEt; // min Et of  electrons
  double fEleTrkMinPt; // min Pt of Tracks corresponding to CEM electrons
  double fCemEleMaxHADEM[2];                      
  double fCemEleMaxIsoFr;  
  double fCemEleMaxTrkZ0; // cut on Trk Z0, beam-constrained track
  int fCemEleFidCut; // cut on CEM fiduciality
  int fCemEleMinCotAxSl; // cut on number of COT Axial SL with 5 or more hits 
  int fCemEleMinCotStSl; // cut on number of COT Stereo SL with 5 or more hits 

  //---- cuts for PEM electrons
  double fPemEleMinEtaPes; // min PES Eta of PEM electrons
  double fPemEleMaxEtaPes; 
  double fPemEleMaxHADEM;      
  int fPemEle3x3FitCut; // cut on PEM 3x3 tower fit status code
  double fPemEleMax3x3Chi2;
  double fPemEleMinPes5x9U; // cut on 5x9 shower profile ratios (U)
  double fPemEleMinPes5x9V; // cut on 5x9 shower profile ratios (V)
  double fPemEleMaxIsoFr;  
  double fPemEleMaxPesdR;

  //---- cuts for Jets
  double MinJetEta;
  double MaxJetEta;

  //______________________________________ other membres and functions
  int UnpackBlocks(TStnEvent* event);
  int Event(int ientry);
  int EndJob();

  int GetMyNvx12(double &zvx); // accessing Vertex info
  void ClearCommonStuff(CommonStuff &stuff); // clears CommonStuff structures 
  void ClearExternalInput(); // clears external input
  void ClearJetStuff(JetStuff &stuff);     // clears the JetStuff structures
  void ClearMetsigStuff(MetsigStuff &sigstuff); // clears MetsigStuff
  void DoCommonStuff(CommonStuff &miscstuff); // fills out CommonStuff
  void DoPhotons(CommonStuff &miscstuff); // fills out photons 
  void DoElectrons(CommonStuff &miscstuff); // fills out electrons 
  void AddExternalElectrons(CommonStuff &miscstuff); // adds external electrons to the list
  void DoMuons(CommonStuff &miscstuff); // fills out muons 
  void AddExternalMuons(CommonStuff &miscstuff); // adds external muons to the list

  void DoMyJet(TStnJetBlock* fJetBlock, CommonStuff miscstuff, JetStuff &jetstuff); // does my jets
  void DoMyJetNoMatch(TStnJetBlock* fJetBlock, JetStuff &jetstuff); // does my jets, but doesn't remove EM objetcs
  void DoMyJetWithMatch(TStnJetBlock* fJetBlock, CommonStuff miscstuff, JetStuff &jetstuff); // does my jets after removing EM objects
  void MyNewJetWithMet(JetStuff &jetstuff); // needed for pseudo-exp: adds Met to closest jet, the output is "new" JetLev6. 
  void DoMyMet(CommonStuff miscstuff, JetStuff &jetstuff); // corrects Met & SumEt for jets 
  void DoMyMetSignificanceData(CommonStuff miscstuff, JetStuff jetstuff, MetsigStuff &sigstuff); // does metsig for data & MC 
  void DoMyMetSignificancePseudo(CommonStuff miscstuff, JetStuff jetstuff, MetsigStuff &sigstuff); // does metsig for pseudo-experiments 

  //------------- metsignificance functions
  double MyTotalMetSignificance(JetStuff jetstuff, TVector2 myMetVec, int systcode, int systcodeUncl); // My Met significance for all events
  void InitMetProbStuff(double &metprobInt,JetStuff jetstuff, TVector2 myMetVec, int systcode, int systcodeUncl); // inits params for MetProb;
  void 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); // return parameterization of MET due to Uncl. Energy
  double JERCorrCoeff(int i, int j); // returns correlation coefficients for JER params
  double JERparam(int code, int i, int j); // returns JER param or its uncertainty
  int WhatJEReta(double eta_det); // returns index for eta_bin in JER
  double MyJer_meanG(double jet_E, double eta_det, int stat_code, int syst_code); // returns mean of Gauss in JER: A0+A1*Pt+A2/E
  double MyJer_sigmaG(double jet_E, double eta_det, int stat_code, int syst_code); // returns sigma of Gauss in JER: sqrt(A3/E+A4)
  double MyJer_mpvL(double jet_E, double eta_det, int stat_code, int syst_code); // returns mpv of Landau in JER: A5+A6*Pt+A7/E
  double MyJer_sigmaL(double jet_E, double eta_det, int stat_code, int syst_code); // returns sigma of Landau in JER: sqrt(A8/E+A9)
  double MyJer_normG(double jet_E, double eta_det, int stat_code, int syst_code); // returns normalization of Gauss: [A10+A11*sqrt(E)]/E+A12
  double JER_DataMCdiff(double jet_E, double eta_det, int syst_code); // JER sytematics due to Data-MC difference in JER (see e-log entry as of 01/08/07)
  double MyIntegralUpLimit(double jet_E, double eta_det, int stat_code, int syst_code); // upper integration limit for JER function
  double MetSigCorrection(double rawSig); // returns corrected MetSig
  double MetSigCorrPAR(int ind, int njet15, int isMC, int sample_id); // returns input parameters for MetSig correction 

  //___ This is new function to generate Met; it takes care of unclustered & jet components of Met.
  void GenerateMyTotalMet(JetStuff jetstuff, CommonStuff miscstuff, int systcode, int systcodeUncl, int rnd_seed, TVector2 &myGenMet); 

  void ReorderMyJets(JetStuff &jetstuff); // reorders jets after removing EM objects
  double CorrectEtaDet(TLorentzVector *vec_pho, TLorentzVector *vec_old, 
		       double jetetadet_old, double pho_etadet); // re-calculates jet EtaDet after removing EM object
  void myExchange_tlv(TLorentzVector& val1, TLorentzVector& val2); //exchanges val1 and val2
  void myExchange_dbl(double& val1, double& val2); //exchanges val1 and val2
  void myExchange_int(int& val1, int& val2); //exchanges val1 and val2
  //------------------ need this for jet-EM matching
  void MatchCalorTowers(int jet_ind,TStnJetBlock* fJetBlock,TCalDataBlock *fCalDataBlock,CalDataArray* cdholder); // for jets
  void MatchCalorTowers(int pho_ind,TStnPhotonBlock* fPhotonBlock,TCalDataBlock *fCalDataBlock,CalDataArray* cdholder); // for photons
  void MatchCalorTowers(int ele_ind,TStnElectronBlock* fElectronBlock,TCalDataBlock *fCalDataBlock,CalDataArray* cdholder); // for electrons
  int MyMatchPhoJet(int jet_ind, int pho_ind, int ele_ind, TStnJetBlock* fJetBlock,
		    TLorentzVector *mypho, TLorentzVector *myjet); // matches Jet with EM object
  double GetCorrection(CorInit settings, TLorentzVector vec, float& emf, float etad); /* this function returns correction 
											 based on parameters localy set in the Module */ 
  double GetCorrection(TLorentzVector vec, float& emf, float etad); /* this function returns correction based 
								       on globaly initialized parameters */ 


  int CEMPhotonIDcut(TStnPhoton* Pho); // CEM photon ID
  int CEMElectronLooseIDcut(TStnElectron* Ele, TStnTrack* Trk); // loose CEM electron ID 
  int PEMElectronLooseIDcut(TStnElectron* Ele); // loose PEM electron ID
  int IsPhoEle(TStnElectron* Ele); // prevents double count of ele==pho: returns 1 if ele==pho, 0--otherwise

  double myCorrIso(TStnElectron* ele); // returns Iso corrected for leakage & nvx12
  double myCorrIsoFr(TStnElectron* ele); // returns Iso/Et(ele) corrected for leakage & nvx12
  double myCorrEtEle(TStnElectron* ele); // returns corrected Et of electrons (includes all corrections)
  double myCorrEnergyEle(TStnElectron* ele); // returns corrected E of electrons (includes all corrections)

  void PrintPhoTowers(int pho_ind,TStnPhotonBlock* fPhotonBlock,TCalDataBlock *fCalDataBlock); // prints photon towers
  void PrintEleTowers(int ele_ind,TStnElectronBlock* fElectronBlock,TCalDataBlock *fCalDataBlock); // prints electron towers
  void PrintJetTowers(int jet_ind,TStnJetBlock* fJetBlock,TCalDataBlock *fCalDataBlock); // prints jet towers

  double CalculateMetDelPhi(int ind); // calculates dPhi for metsig.vs.dPhi studies
  double CalculateMetDelPhi(int ind, double &Et_closest); // calculates dPhi and returns Et(closest object), for metsig.vs.dPhi studies

  //-----------------------------------------------------------------------------------------------------------------------------
  //
  //   >>>>>>>>>>>>>   Public methods <<<<<<<<<<<<<<<<<<<<<<< 
  //
  //_____________________________________________________________________________________________________________________________

public:
  TStnMetModel();
  ~TStnMetModel();

  int RunMetModel(TStnEvent* event); // Main Method to run MetModel

  //________ ****** accessors to output			
  double GetMetsigData(); // returns default value of metsig for data event
  double GetMetsigSyst(int syst_ind); // 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
  TVector2 GetMetData(); // returns corrected MET as it's calculated by Met Model
  int GetNpseudo(); // returns number of pseudo-experiments for MetModel
  double GetMetsigPseudo(int indPseudo); // same as GetMetsigData(), but if we choose to use predictions based on pseudo-experiments
  double GetMetsigPseudoSyst(int indPseudo, int syst_ind); // same as GetMetsigSyst, but we choose to use predictions based on pseudo-experiments   
  TVector2 GetMetPseudo(int indPseudo); // returns MET predicted by Met Model according to default parameterization
  TVector2 GetMetPseudoSyst(int indPseudo, int syst_ind); // returns deviated MET predicted by Met Model for systematics studies 
  int IsBadMet(); // flag: 0=no bad met; 1=bad met. Output of my MET cleanup cuts tuned to remove events with lost EM objects
  double GetDelPhiClosest(); // min(dPhi=Phi(met)-Phi(obj))
  double GetDelPhiClosestJet(); // min(dPhi=Phi(met)-Phi(jet))
  double GetMetJetRatio(); // MET*dPhiClosest/Et(closest jet), dPhiClosest=min(dPhi=Phi(met)-Phi(jet))
  double GetEtClosestJet(); // Et(closest jet to MET), dPhiClosest=min(dPhi=Phi(met)-Phi(jet))
  double GetMetsigSignedData(); // returns default value of "signed" metsig for data event
  double GetMetsigSignedSyst(int syst_ind); // returns deviated values of "signed" 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

		
  //________ cuts & input params
  void PrintEvent(); // prints out event info
  void SetDebugMode(int mode) { DebugMode=mode; }
  void SetJTC_systcode(int param) { fJTC_systcode = param; } // for JES systematics studies
  void SetSampleID(int IDind) { fSampleID = IDind; } // sets sample ID to be used in metsig calibration
  void SetUnclParamSwitch(int param) { fUnclParamSwitch = param; } // sets Unclustered Energy parameterization
  void SetNpseudo(int Nexp) { Npseudo=Nexp; } // number of pseudo-experiments for MetModel

  //___________________________________________ ****** external setters for intput params
  void SetNele(int nele);
  void SetEleInd(int i, int ind);
  void SetCorrEle(int i, TLorentzVector vec); 
  void SetNmuon(int nmu);
  void SetTrkInd(int i, int ind);
  void SetCorrMuon(int i, TLorentzVector vec);
  void SetRawMuon(int i, TLorentzVector vec);
  void SetTrkAlgoID(int i, int id_ind);
  void SetTrkZvx(int i, double zvx);
  void SetCalCluster(int i, TLorentzVector vec); 

  ClassDef(TStnMetModel,0)
};

#endif