////////////////////////////////////////////////////////////////////////////////// 
 /// Analysis skeleton
 /// also for distributed analysis examples
 /// Author: Ketevi A. Assamagan
 /// BNL, July 22, 2004
 ///
 /// DESCRIPTION:
 ///
 /// This class is an analysis skeleton - The user can implement his analysis here
 /// This class is also used for the demonstration of the distributed analysis
 /// Some electron histograms are used for the distributed case. The user may
 /// remove the histograms and the electron stuff if not needed.
 /// Note t:he single algorithm structure as an analysis code does not scale
 /// For detailed analysis examples, look in CVS: PhysicsAnalysis/AnalysisCommon/AnalysisExamples/
 /// Ketevi A. Assamagan on June 22, 2004
 ///
 ///////////////////////////////////////////////////////////////////////////////////////////////////////
 
 
 //////////////////////////////////////////////////////////////////////////////////
 
 #include "GaudiKernel/MsgStream.h"
 #include "GaudiKernel/AlgFactory.h"
 #include "GaudiKernel/IToolSvc.h"
 
 #include "StoreGate/StoreGateSvc.h"
 
 #include "egammaEvent/ElectronContainer.h"
 #include "egammaEvent/EMShower.h"
 #include "egammaEvent/egammaParamDefs.h"
 
 #include "McParticleEvent/TruthParticleContainer.h"
 
 #include "VxVertex/VxContainer.h"
 #include "Particle/TrackParticleContainer.h"
 #include "CaloEvent/CaloClusterContainer.h"
 
 #include "muonEvent/MuonContainer.h"
 #include "egammaEvent/PhotonContainer.h"
 #include "tauEvent/TauJetContainer.h"
 #include "JetEvent/JetCollection.h"
 #include "MissingETEvent/MissingET.h"
 
 #include "NavFourMom/IParticleContainer.h"
 #include "NavFourMom/INavigable4MomentumCollection.h"
 
 #include "GaudiKernel/ITHistSvc.h"
 #include "TTree.h"
 #include "CLHEP/Vector/LorentzVector.h"
 
 //#include "JetTagInfo/TruthInfo.h"
 
 #include "UserAnalysis/AnalysisSkeleton.h"
 
 #include <algorithm>
 #include <math.h>
 #include <functional>
 #include <iostream>
 
 static const double mZ = 91.19*GeV;
 static const int  MAX_PARTICLES = 20;
 
 using namespace Analysis;
 using namespace Rec;
 
 //////////////////////////////////////////////////////////////////////////////////////
 /// Constructor
 
 AnalysisSkeleton::AnalysisSkeleton(const std::string& name,
   ISvcLocator* pSvcLocator) : CBNT_AthenaAwareBase(name, pSvcLocator),
   m_analysisTools( "AnalysisTools" ),
   m_analysisSelectionTool( "UserAnalysisSelectionTool" ),
   m_analysisPreparationTool( "UserAnalysisPreparationTool" ),
   m_analysisOverlapCheckingTool( "UserAnalysisOverlapCheckingTool" ),
   m_analysisOverlapRemovalTool( "UserAnalysisOverlapRemovalTool" ),
   m_trigDec("TrigDec::TrigDecisionTool"){
 
   /** switches to control the analysis through job options */
 
   declareProperty( "AnalysisTools",               m_analysisTools );
   declareProperty( "AnalysisSelectionTool",       m_analysisSelectionTool);
   declareProperty( "AnalysisPreparationTool",     m_analysisPreparationTool);
   declareProperty( "AnalysisOverlapCheckingTool", m_analysisOverlapCheckingTool);
   declareProperty( "AnalysisOverlapRemovalTool",  m_analysisOverlapRemovalTool);
   declareProperty( "TrigDecisionTool",            m_trigDec, 
                    "The tool to access TrigDecision");
 
   declareProperty("ElectronContainer", m_electronContainerName="ElectronAODCollection"); 
   declareProperty("McParticleContainer", m_truthParticleContainerName = "SpclMC");
   declareProperty("MissingETObject",m_missingETObjectName="MET_RefFinal");
 
   /** the cuts - default values - to be modified in job options */
 
   declareProperty("DeltaRMatchCut", m_deltaRMatchCut = 0.2);
   declareProperty("MaxDeltaR", m_maxDeltaR = 0.9999);
 
   /** additiona cuts for electrons */
   declareProperty("ElectronEtCut", m_etElecCut = 10.0*GeV);
   declareProperty("ElectronEtaCut", m_etaElecCut = 2.5);
   declareProperty("ElectronCone", m_elecCone = 0.9);
 
   /** additional cuts for bjet tagging */
   declareProperty("bjetWt_IP3DSV1Cut", m_bjetWt_ip3dsv1Cut = 6);
   declareProperty("bjet_etaCut", m_bjet_etaCut = 2.5);
   declareProperty("bjet_etCut", m_bjet_etCut = 15.0*GeV);
 
   /** missing ET options */
   declareProperty("MissingETCut",m_missingETCut=20.0*GeV);
 
   /** is this AtlFast */
   declareProperty("IsAtlFastData",m_isAtlFastData="False");
 
   /** count number of jets with ET > min value - for SUSY studies */
   declareProperty("SusyJetMinEt", m_SusyJetMinEt = 50*GeV);
 
 
   declareProperty("DoTrigger", m_doTrigger = true, "enable trigger example");
 }
 
 /////////////////////////////////////////////////////////////////////////////////////
 /// Destructor - check up memory allocation
 /// delete any memory allocation on the heap
 
 AnalysisSkeleton::~AnalysisSkeleton() {}
 
 ////////////////////////////////////////////////////////////////////////////////////
 /// Initialize
 /// initialize StoreGate
 /// get a handle on the analysis tools
 /// book histograms
 
 StatusCode AnalysisSkeleton::CBNT_initializeBeforeEventLoop() {
   MsgStream mLog( messageService(), name() );
 
   mLog << MSG::DEBUG << "Initializing AnalysisSkeleton (before eventloop)" << endreq;
 
   // retrieve trigger decision tool
   // needs to be done before the first run/event since a number of
   // BeginRun/BeginEvents are registered by dependent services
   StatusCode sc = StatusCode::SUCCESS;
 
 /*
   if ( m_doTrigger ) {
      sc = m_trigDec.retrieve();
      if ( sc.isFailure() ){
         mLog << MSG::ERROR << "Can't get handle on TrigDecisionTool" << endreq;
      } else {
        mLog << MSG::DEBUG << "Got handle on TrigDecisionTool" << endreq;
      }
   }
 */
 
   return sc;
 } 
 
 StatusCode AnalysisSkeleton::CBNT_initialize() {
 
   MsgStream mLog( messageService(), name() );
 
   mLog << MSG::DEBUG << "Initializing AnalysisSkeleton" << endreq;
 
   /** get a handle of StoreGate for access to the Event Store */
   StatusCode sc = service("StoreGateSvc", m_storeGate);
   if (sc.isFailure()) {
      mLog << MSG::ERROR
           << "Unable to retrieve pointer to StoreGateSvc"
           << endreq;
      return sc;
   }
   
   /// get a handle on the analysis tools
   sc = m_analysisTools.retrieve();
   if ( sc.isFailure() ) {
       mLog << MSG::ERROR << "Can't get handle on analysis tools" << endreq;
       return sc;
   }
 
   /// get a handle on the preparartion analysis tools
   sc = m_analysisSelectionTool.retrieve();
   if ( sc.isFailure() ) {
       mLog << MSG::ERROR << "Can't get handle on analysis selection tool" << endreq;
       return sc;
   }
 
   sc = m_analysisPreparationTool.retrieve();
   if ( sc.isFailure() ) {
       mLog << MSG::ERROR << "Can't get handle on analysis preparation tool" << endreq;
       return sc;
   }
 
   sc = m_analysisOverlapCheckingTool.retrieve();
   if ( sc.isFailure() ) {
       mLog << MSG::ERROR << "Can't get handle on analysis overlap checking tool" << endreq;
       return sc;
   }
 
   sc = m_analysisOverlapRemovalTool.retrieve();
   if ( sc.isFailure() ) {
       mLog << MSG::ERROR << "Can't get handle on analysis overlap removal tool" << endreq;
       return sc;
   }
 
   if ( m_doTrigger ) {
     sc = m_trigDec.retrieve();
     if ( sc.isFailure() ){
       mLog << MSG::ERROR << "Can't get handle on TrigDecisionTool" << endreq;
       return sc;
     }
   }
 
   /** get a handle on the NTuple and histogramming service */
   sc = service("THistSvc", m_thistSvc);
   if (sc.isFailure()) {
      mLog << MSG::ERROR
           << "Unable to retrieve pointer to THistSvc"
           << endreq;
      return sc;
   }
 
   /** Print out bjet cut values */
   mLog << MSG::DEBUG << "b jet cuts: Et/eta/IP3DSV1 wt. "<<m_bjet_etCut<< ","<<m_bjet_etaCut <<","<<m_bjetWt_ip3dsv1Cut<<endreq;
 
 
   /** now add branches and leaves to the AAN tree */
 
   addBranch("NElectrons",    m_aan_size, "NElectrons/i");
   addBranch("ElectronEta",   m_aan_eta);
   addBranch("ElectronPt",    m_aan_pt);
   addBranch("ElecPtRatio",   m_aan_elecetres);
 
   addBranch("NJets",          m_aan_njets, "NJets/i");
   addBranch("NJets_etaLT25",  m_aan_njets_etaLT25, "NJets_etaLT25/i");
   addBranch("NJets_SusyETCut",m_aan_njets_SusyETCut, "NJets_SusyETCut/i");
   addBranch("JetsEta"        ,m_aan_JetEta);
   addBranch("JetsEt"         ,m_aan_JetEt);
   addBranch("JetsBTagWt"     ,m_aan_JetBTagWt);
   addBranch("MissingET",      m_aan_ptMiss, "MissingET/d");
   addBranch("EffMass",        m_aan_effmass, "EffMass/d");
   addBranch("HT",             m_aan_ht,"HT/d");
   addBranch("NBJets",         m_aan_nbjets, "NBJets/i");
   addBranch("MaxJetET",       m_aan_maxJetET, "MaxJetET/d");
 
   addBranch("NFinalEl",       m_aan_NFinalEl, "NFinalEl/i");
   addBranch("FinalElEta",     m_aan_FinalElEta);
   addBranch("FinalElPt",      m_aan_FinalElPt);
   addBranch("FinalElEtCone20",m_aan_FinalElEtCone20);
 
   addBranch("NFinalMu",       m_aan_NFinalMu, "NFinalMu/i");
   addBranch("FinalMuEta",     m_aan_FinalMuEta);
   addBranch("FinalMuPt",      m_aan_FinalMuPt);
   addBranch("FinalMuEtCone20",m_aan_FinalMuEtCone20);
   addBranch("FinalMuBestMat", m_aan_FinalMuBestMat);
   addBranch("FinalMuMatChi2", m_aan_FinalMuMatChi2);
 
   addBranch("FinalLepEtSum",  m_aan_FinalLepEtSum, "FinalLepEtSum/d");
   addBranch("FinalElEtSum",   m_aan_FinalElEtSum, "FinalElEtSum/d");
   addBranch("FinalMuEtSum",   m_aan_FinalMuEtSum, "FinalMuEtSum/d");
 
   addBranch("NumberTopQ",     m_aan_NumTopQ, "NumberTopQ/i");
   addBranch("pTtop1",         m_aan_pTtop1,  "pTtop1/d");
   addBranch("pTtop2",         m_aan_pTtop2,  "pTtop2/d");
 
 
   /// ROOT histograms ---------------------------------------
 
   /// electrons
   m_h_elecpt     = new TH1F("elec_pt","pt el",50,0,250.*GeV);
   sc = m_thistSvc->regHist("/AANT/Electron/elec_pt",m_h_elecpt);
 
   m_h_eleceta    = new TH1F("elec_eta","eta el",70,-3.5,3.5);
   sc = m_thistSvc->regHist("/AANT/Electron/elec_eta",m_h_eleceta);
 
   m_h_elec_deltaRMatch    = new TH1F("elec_deltaRMatch","elec reco/MC deltaR",50,0.,1.);
   sc = m_thistSvc->regHist("/AANT/Electron/elec_deltaRMatch",m_h_elec_deltaRMatch);
 
   /// jets - before OverlapRemoval
   m_h_jet_eta_beforeOR = new TH1F("jet_eta_beforeOR","jet_eta before OR",50,-5.,5.);
   sc = m_thistSvc->regHist("/AANT/Jet/jet_eta_beforeOR",m_h_jet_eta_beforeOR);
 
   m_h_jet_et_beforeOR = new TH1F("jet_et_beforeOR","jet_et before OR",100,0.,500.);
   sc = m_thistSvc->regHist("/AANT/Jet/jet_et_beforeOR",m_h_jet_et_beforeOR);
 
   m_h_jet_ip3dsv1Wt_beforeOR = new TH1F("jet_ip3dsv1Wt_beforeOR","jet_ip3dsv1Wt before OR",120,-20.,40.);
   sc = m_thistSvc->regHist("/AANT/Jet/jet_ip3dsv1Wt_beforeOR",m_h_jet_ip3dsv1Wt_beforeOR);
 
   m_h_jet_label_beforeOR = new TH1F("jet_label_beforeOR","jet_label before OR",20,0.,20.);
   sc = m_thistSvc->regHist("/AANT/Jet/jet_label_beforeOR",m_h_jet_label_beforeOR);
 
   m_h_jet_ip3dsv1Wt_bjet_beforeOR = new TH1F("jet_ip3dsv1Wt_bjet_beforeOR","b jet_ip3dsv1Wt before OR",120,-20.,40.);
   sc = m_thistSvc->regHist("/AANT/Jet/jet_ip3dsv1Wt_bjet_beforeOR",m_h_jet_ip3dsv1Wt_bjet_beforeOR);
 
   m_h_jet_ip3dsv1Wt_ujet_beforeOR = new TH1F("jet_ip3dsv1Wt_ujet_beforeOR","u jet_ip3dsv1Wt before OR",120,-20.,40.);
   sc = m_thistSvc->regHist("/AANT/Jet/jet_ip3dsv1Wt_ujet_beforeOR",m_h_jet_ip3dsv1Wt_ujet_beforeOR);
 
   /// jets - after OverlapRemoval
   m_h_jet_eta_afterOR = new TH1F("jet_eta_afterOR","jet_eta after OR",50,-5.,5.);
   sc = m_thistSvc->regHist("/AANT/Jet/jet_eta_afterOR",m_h_jet_eta_afterOR);
 
   m_h_jet_et_afterOR = new TH1F("jet_et_afterOR","jet_et after OR",100,0.,500.);
   sc = m_thistSvc->regHist("/AANT/Jet/jet_et_afterOR",m_h_jet_et_afterOR);
 
   m_h_jet_ip3dsv1Wt_afterOR = new TH1F("jet_ip3dsv1Wt_afterOR","jet_ip3dsv1Wt after OR",120,-20.,40.);
   sc = m_thistSvc->regHist("/AANT/Jet/jet_ip3dsv1Wt_afterOR",m_h_jet_ip3dsv1Wt_afterOR);
 
   m_h_jet_label_afterOR = new TH1F("jet_label_afterOR","jet_label after OR",20,0.,20.);
   sc = m_thistSvc->regHist("/AANT/Jet/jet_label_afterOR",m_h_jet_label_afterOR);
 
   m_h_jet_ip3dsv1Wt_bjet_afterOR = new TH1F("jet_ip3dsv1Wt_bjet_afterOR","b jet_ip3dsv1Wt after OR",120,-20.,40.);
   sc = m_thistSvc->regHist("/AANT/Jet/jet_ip3dsv1Wt_bjet_afterOR",m_h_jet_ip3dsv1Wt_bjet_afterOR);
 
   m_h_jet_ip3dsv1Wt_ujet_afterOR = new TH1F("jet_ip3dsv1Wt_ujet_afterOR","u jet_ip3dsv1Wt after OR",120,-20.,40.);
   sc = m_thistSvc->regHist("/AANT/Jet/jet_ip3dsv1Wt_ujet_afterOR",m_h_jet_ip3dsv1Wt_ujet_afterOR);
 
   /// missing ET
 
   m_pxMis   = new TH1F("MissingPx", "MissingPx",200,-500.0*GeV,500.*GeV);
   sc = m_thistSvc->regHist("/AANT/MissingET/MissingPx", m_pxMis);
   m_pyMis   = new TH1F("MissingPy","MissingPy",200,-500.0*GeV,500.*GeV);
   sc = m_thistSvc->regHist("/AANT/MissingET/MissingPy", m_pyMis);
   m_ptMis   = new TH1F("MissingPt","MissingPt",100,0.0,500.*GeV);
   sc = m_thistSvc->regHist("/AANT/MissingET/MissingPt", m_ptMis);
 
 
   if (sc.isFailure()) { 
      mLog << MSG::ERROR << "ROOT Hist registration failed" << endreq; 
      return sc; 
   }
   /// end ROOT Histograms ------------------------------------------
 
   return StatusCode::SUCCESS;
 }                
 
 ///////////////////////////////////////////////////////////////////////////////////
 /// Finalize - delete any memory allocation from the heap
 
 StatusCode AnalysisSkeleton::CBNT_finalize() {
   MsgStream mLog( messageService(), name() );
   
   return StatusCode::SUCCESS;
 
 }
 
 ///////////////////////////////////////////////////////////////////////////////////
 /// Clear - clear CBNT members
 StatusCode AnalysisSkeleton::CBNT_clear() {
   /// For Athena-Aware NTuple
 
   m_aan_size = 0;
   m_aan_eta->clear();
   m_aan_pt->clear();
   m_aan_elecetres->clear();
 
   m_aan_njets=0;
   m_aan_maxJetET = -1000.;
   m_aan_JetEta->clear();
   m_aan_JetEt->clear();
   m_aan_JetBTagWt->clear();  
   //
   m_aan_njets_etaLT25=0;
   m_aan_njets_SusyETCut = 0;
 
   m_aan_ptMiss = -1.;
   m_aan_effmass = 0.;
   m_aan_ht = 0;
   m_aan_nbjets = 0;
 
   //
   m_aan_NFinalEl = 0;
   m_aan_FinalElPt->clear();
   m_aan_FinalElEta->clear();
   m_aan_FinalElEtCone20->clear();
 
 
   m_aan_NFinalMu = 0;
   m_aan_FinalMuPt->clear();
   m_aan_FinalMuBestMat->clear();
   m_aan_FinalMuEta->clear();
   m_aan_FinalMuEtCone20->clear();
   m_aan_FinalMuMatChi2->clear();
 
   //
   m_aan_FinalLepEtSum = 0.;
   m_aan_FinalElEtSum  = 0.;
   m_aan_FinalMuEtSum  = 0.;
 
   // 
   m_aan_NumTopQ=0;
   m_aan_pTtop1=-1;
   m_aan_pTtop2=-1;
 
   return StatusCode::SUCCESS;
 }
 
 //////////////////////////////////////////////////////////////////////////////////
 /// Execute - on event by event
 
 StatusCode AnalysisSkeleton::CBNT_execute() {
   MsgStream mLog( messageService(), name() );
   
   mLog << MSG::DEBUG << "in execute()" << endreq;
 
   StatusCode sc;
 
   /** it shows how to get the Electron and the TruthParticle
       containers shows matching between reconstructed and MC electrons
 
       Can be commented out if you want to do so If you do so, some of
       the electron histograms/ntuple variables will be unfilled */
 
   /*
   sc = electronSkeleton();
   if (sc.isFailure()) {
     mLog << MSG::WARNING << "The method electronSkeleton() failed" << endreq;
     return StatusCode::SUCCESS;
   }
   */
 
   /** an minimal example using the TrigDecisionTool */
   if ( m_doTrigger ) {
     sc = triggerSkeleton();
     if (sc.isFailure()) {
       mLog << MSG::WARNING << "The method triggerSkeleton() failed" << endreq;
       return StatusCode::SUCCESS;
     }
   }
  
   /** an example of analysis preparation, including:
       - pre-selections based on the reocmmendations of performance groups
       - overlap checking
       - overlap removal */
 
   /** Do not comment the next method. This is where we do all the
       selection/overlap removal Those results are then used in the
       methods later on */
 
   sc = analysisPreparation();
   if ( sc.isFailure() ) {
      mLog << MSG::WARNING << "Analysis Preparation Failed " << endreq;
      return StatusCode::SUCCESS;
   }
 
   /** The following methods were added by Vivek Jain They just show
       you how to access different variables and store them in
       histograms and/or ntuples */
 
   /** look at bjet tagging information in the jets after overlap
       removal */
 
   sc = bjetInfo();
   if ( sc.isFailure() ) {
      mLog << MSG::WARNING << "Failure in bjetInfo " << endreq;
      return StatusCode::SUCCESS;
   } 
 
   /** get missing Et information */
 
   sc = getMissingET();
   if( sc.isFailure() ) {
     mLog << MSG::WARNING
          << "Failed to retrieve Et object found in TDS"
          << endreq; 
     return StatusCode::SUCCESS;
   }  
 
   /** do SUSY studies */
 
   sc = SusyStudies();
   if( sc.isFailure() ) {
     mLog << MSG::WARNING
          << "Failed to do SUSY studies"
          << endreq; 
     return StatusCode::SUCCESS;
   }  
 
 
 
   return StatusCode::SUCCESS;
 }
 
 //////////////////////////////////////////////////////////////////////////////////
 /// Trigger method - called by execute() on event by event
 /// to be removed if not needed
 
 StatusCode AnalysisSkeleton::triggerSkeleton() {
   MsgStream mLog( messageService(), name() );
   mLog << MSG::DEBUG << "in triggerSkeleton()" << endreq;
 
   // for example, did event pass Event Filter ?
   // needs to be changed to m_trigDec->isPhysicsPassed
   mLog << MSG::INFO << "Pass state EF = " << m_trigDec->isPassed(TrigDec::EF) << endreq;
 
   // for example, did event pass !L2_e25i chain? 
   // needs to be changed to m_trigDec->isPhysicsPassed
   std::string mychain("L2_e25i");
   const HLT::Chain* chain = m_trigDec->getHLTChain(mychain);
   if (0 == chain){
     mLog << MSG::INFO << "Chain " << mychain << " is not defined";
   } else {
     mLog << MSG::INFO << "Chain " << mychain << ": " << *chain << " passed: " << chain->chainPassed() << endreq;
     
     mLog << MSG::DEBUG << "triggerSkeleton() succeeded" << endreq;
     
     return StatusCode::SUCCESS;
   }
   return StatusCode::SUCCESS;
 }
 
 //////////////////////////////////////////////////////////////////////////////////
 /// Electron method - called by execute() on event by event
 /// to be removed if not needed
 
 StatusCode AnalysisSkeleton::electronSkeleton() {
   MsgStream mLog( messageService(), name() );
   
   mLog << MSG::DEBUG << "in electronSkeleton()" << endreq;
 
   /** get the MC truth particle AOD container from StoreGate */
   const TruthParticleContainer*  mcpartTES = 0;
   StatusCode sc=m_storeGate->retrieve( mcpartTES,