mc08.106445.Pythia_RPV_bil1.recon.AOD.e362_s462_r563_tid028084

My first approach

Set up athena in the release corresponding to the AOD reconstruction. After that, I used a code for top physics and did a trigger study. AnalysisSkeleton.h

Header file.

#include "GaudiKernel/ToolHandle.h"
#include "GaudiKernel/Algorithm.h"
#include "GaudiKernel/ObjectVector.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include "StoreGate/StoreGateSvc.h"
#include "GaudiKernel/ITHistSvc.h"

#include "AnalysisTools/AnalysisTools.h"

#include "UserAnalysisUtils/UserAnalysisSelectionTool.h"
#include "UserAnalysisUtils/UserAnalysisPreparationTool.h"
#include "UserAnalysisUtils/UserAnalysisOverlapCheckingTool.h"
#include "UserAnalysisUtils/UserAnalysisOverlapRemovalTool.h"

#include "TrigDecision/TrigDecisionTool.h"

//#include "TLorentzVector.h"
//#include "CLHEP/Vector/LorentzVector.h"

#include "CBNT_Utils/CBNT_AthenaAwareBase.h"

#include <string>

#include "TH1.h"

class JetCollection;

class AnalysisSkeleton : public CBNT_AthenaAwareBase  {

 public:

   AnalysisSkeleton(const std::string& name, ISvcLocator* pSvcLocator);
   ~AnalysisSkeleton();

   virtual StatusCode CBNT_initializeBeforeEventLoop();
   virtual StatusCode CBNT_initialize();
   virtual StatusCode CBNT_finalize();
   virtual StatusCode CBNT_execute();
   virtual StatusCode CBNT_clear();

 private:

   /** methods called by CBNT_execute() */
   StatusCode electronSkeleton();
   StatusCode triggerSkeleton();

   /** an example of pre-selection, overlap-checking and overlap removal */
   StatusCode analysisPreparation();

   /** look at b-jet tagging information */
   StatusCode bjetInfo();

   /** get quark flavour of jets */
   int getQuarkJetFlavour(JetCollection::const_iterator jetItr);

   /** get missing ET information */
   StatusCode getMissingET();

   /** make plots for SUSY studies */
   StatusCode SusyStudies();

   /** get pT of top quarks */
   StatusCode getTopQpT(int &, double&, double&);

 private:

   /** get a handle to the tool helper */
   ToolHandle<AnalysisTools> m_analysisTools;

   /** get a handle on the user tool for pre-selection and overlap removal */
   ToolHandle<UserAnalysisSelectionTool>       m_analysisSelectionTool;
   ToolHandle<UserAnalysisPreparationTool>     m_analysisPreparationTool;
   ToolHandle<UserAnalysisOverlapCheckingTool> m_analysisOverlapCheckingTool;
   ToolHandle<UserAnalysisOverlapRemovalTool>  m_analysisOverlapRemovalTool;

   /** tool to access the trigger decision */
   ToolHandle<TrigDec::TrigDecisionTool> m_trigDec;

   /** a handle on the Hist/TTree registration service */
   ITHistSvc * m_thistSvc;

   /** a handle on Store Gate for access to the Event Store */
   StoreGateSvc* m_storeGate;

   /** the key of the Electron Container to retrieve from the AOD */
   std::string m_electronContainerName; 

   /** name of the AOD truth particle container to retrieve from StoreGate */
   std::string m_truthParticleContainerName;

   /** key to get missing ET information */
   std::string m_missingETObjectName;

   /// The missing ET object
   const MissingET * m_pMissing;
   double m_pxMiss;
   double m_pyMiss;
   double m_ptMiss;   

   /** additional user cuts after pre-selections */ 
 
   double m_deltaRMatchCut;
   double m_maxDeltaR;

   /** electron specific cuts */
   double m_etElecCut;
   double m_elecCone;
   double m_etaElecCut;

   /** bjet specific cuts */
   double m_bjetWt_ip3dsv1Cut;
   double m_bjet_etaCut;
   double m_bjet_etCut;

   /** missing ET cuts */
   double m_missingETCut;

   /** Atlfast data? */
   bool m_isAtlFastData;

   /** min Jet ET cut for SUSY studies */
   double m_SusyJetMinEt;
   
  /** Histograms */
  TH1F* m_h_elecpt;
  TH1F* m_h_eleceta;
  TH1F* m_h_elec_deltaRMatch;

  TH1F* m_h_jet_eta_beforeOR;
  TH1F* m_h_jet_et_beforeOR;
  TH1F* m_h_jet_ip3dsv1Wt_beforeOR;
  TH1F* m_h_jet_label_beforeOR;
  TH1F* m_h_jet_ip3dsv1Wt_bjet_beforeOR;
  TH1F* m_h_jet_ip3dsv1Wt_ujet_beforeOR;

  TH1F* m_h_jet_eta_afterOR;
  TH1F* m_h_jet_et_afterOR;
  TH1F* m_h_jet_ip3dsv1Wt_afterOR;
  TH1F* m_h_jet_label_afterOR;
  TH1F* m_h_jet_ip3dsv1Wt_bjet_afterOR;
  TH1F* m_h_jet_ip3dsv1Wt_ujet_afterOR;

  TH1F* m_pxMis;
  TH1F* m_pyMis;
  TH1F* m_ptMis;

  /** Athena-Aware Ntuple (AAN) variables - branches of the AAN TTree */

  /** Simple variables by Ketevi */
  int m_aan_size;
  std::vector<double> * m_aan_eta;
  std::vector<double> * m_aan_pt;
  std::vector<double> * m_aan_elecetres;

  /** Variables by VJ */
  double m_aan_ptMiss;
  int    m_aan_njets;
  int    m_aan_njets_etaLT25;
  int    m_aan_njets_SusyETCut;
  double m_aan_effmass;
  double m_aan_ht;
  double m_aan_maxJetET;
  int    m_aan_nbjets;

  std::vector<double>* m_aan_JetEta;
  std::vector<double>* m_aan_JetEt;
  std::vector<double>* m_aan_JetBTagWt;

  /** Look at final electrons/muons */
  int m_aan_NFinalEl;
  int m_aan_NFinalMu;

  std::vector<double>* m_aan_FinalElEta;
  std::vector<double>* m_aan_FinalElPt;
  std::vector<double>* m_aan_FinalElEtCone20;
  std::vector<double>* m_aan_FinalElPtrat;

  std::vector<double>* m_aan_FinalMuEta;
  std::vector<double>* m_aan_FinalMuPt;
  std::vector<double>* m_aan_FinalMuEtCone20;
  std::vector<int>*    m_aan_FinalMuBestMat;
  std::vector<double>* m_aan_FinalMuMatChi2;

  double m_aan_FinalLepEtSum;
  double m_aan_FinalElEtSum;
  double m_aan_FinalMuEtSum;

  /** number top quarks */

  int m_aan_NumTopQ;
  double m_aan_pTtop1;
  double m_aan_pTtop2;

  bool m_doTrigger;
};

#endif // ANALYSIS_SKELETON_H

Here the AnalysisSkeleton.cxx code.

Implementation file

#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, m_truthParticleContainerName);
  if( sc.isFailure()  ||  !mcpartTES ) {
     mLog << MSG::WARNING
          << "No AOD MC truth particle container found in TDS"
          << endreq; 
     return StatusCode::SUCCESS;
  }
  mLog <<MSG::DEBUG << "MC Truth Container Successfully Retrieved" << endreq;
  
  /** get the container of the original AOD electron - without any selection */
  /** get the AOD electron container for TES */
  const ElectronContainer* elecTES = 0;
  sc=m_storeGate->retrieve( elecTES, m_electronContainerName);
  if( sc.isFailure()  ||  !elecTES ) {
     mLog << MSG::WARNING
          << "No AOD electron container found in TDS"
          << endreq; 
     return StatusCode::FAILURE;
  }  
  mLog << MSG::DEBUG << "ElectronContainer successfully retrieved - size is " << elecTES->size() << " electrons " << endreq;


  /** iterators over the electron container - the pre-selected ones or the original ones */ 
  ElectronContainer::const_iterator elecItr  = elecTES->begin();
  ElectronContainer::const_iterator elecItrE = elecTES->end();

  for (; elecItr != elecItrE; ++elecItr) {

    /** apply further selections if necessary */
    /** check for the author of the electron */
    bool author = (*elecItr)->author(egammaParameters::AuthorElectron);
    if ( !author || (*elecItr)->pt() < m_etElecCut ) continue;

    m_aan_size++;


    /** fill histograms */
    m_h_elecpt->Fill( (*elecItr)->pt(), 1.);
    m_h_eleceta->Fill( (*elecItr)->eta(), 1.);

    /** fill Athena-Aware NTuple */
    m_aan_eta->push_back((*elecItr)->eta());
    m_aan_pt->push_back((*elecItr)->pt());

    /** find a match to this electron in the MC truth container
        the index and deltaR are returned */
    int index = -1;
    double deltaRMatch;
    if( (*elecItr)->trackParticle() && (*elecItr)->pt()> m_etElecCut ) {
       const TruthParticleContainer * truthContainer = mcpartTES;
       bool findAMatch = m_analysisTools->matchR((*elecItr), truthContainer, 
             index, deltaRMatch, (*elecItr)->pdgId());
       if (findAMatch) {
          deltaRMatch = (deltaRMatch > m_maxDeltaR) ? m_maxDeltaR : deltaRMatch;

          m_h_elec_deltaRMatch->Fill(deltaRMatch);
          mLog << MSG::DEBUG << "Electron: MC/Reco DeltaR " << deltaRMatch << endreq;
          /** check for good match */
          if ( deltaRMatch < m_deltaRMatchCut) {
             const TruthParticle*  electronMCMatch = (*mcpartTES)[index]; 
             double res = (*elecItr)->pt() / electronMCMatch->pt();
             m_aan_elecetres->push_back(res);
          }
       }
    }    
  }

  mLog << MSG::DEBUG << "electronSkeleton() succeeded" << endreq;
        
  return StatusCode::SUCCESS;
}

//////////////////////////////////////////////////////////////////////////////////
/// Analysis Preparation method - called by execute() on event by event
/// A lot of the AOD container are read in
/// pre-selection is done using the UserAnalysisSelectionTool
/// An example of overlap checking is demonstrated
/// An example of overlap removal is demonstration
StatusCode AnalysisSkeleton::analysisPreparation() {

  MsgStream mLog( messageService(), name() );

  mLog << MSG::DEBUG << "in analysisPreparation()" << endreq;

  /** loop over Electrons from the AOD and see which pass the recommended electron selection 
      These selections are defined in m_analysisSelectionTool - to be changed if necessary */
  const ElectronContainer* elecTES = 0;
  StatusCode sc=m_storeGate->retrieve( elecTES, m_electronContainerName);
  if( sc.isFailure()  ||  !elecTES ) {
     mLog << MSG::WARNING
          << "No AOD electron container found in TDS"
          << endreq;
  }
  ElectronContainer::const_iterator elecItr  = elecTES->begin();
  ElectronContainer::const_iterator elecItrE = elecTES->end();
  for (; elecItr != elecItrE; ++elecItr) {
     bool passedSelection = m_analysisSelectionTool->isSelected( *elecItr );
     if ( passedSelection ) mLog << MSG::DEBUG<<"Found a potential good Electron " << endreq;
  }

  /** do analysis preparation using the AnalysisPreparationTool
      selections based or recommended selections from performance groups  
      The tool outputs various containers of pre-selected objects */
  sc = m_analysisPreparationTool->execute();
  if ( sc.isFailure() ) {
    mLog << MSG::WARNING << "AnalysisPreparation Failed - selection " << endreq;
    return StatusCode::SUCCESS;
  }


  /** get the pre-selected Electrons - given by the AnalysisPreparationTool */
  const ElectronContainer* preselectedElecTES = m_analysisPreparationTool->selectedElectrons();
  if ( !preselectedElecTES ) {
    mLog << MSG::WARNING << "Selected Electrons Not Found " << endreq;
  }
  mLog << MSG::DEBUG << "Pre-selected Electrons successfully retrieved - size is " << preselectedElecTES->size() << " electrons " << endreq;

  /** get the pre-selected Muons - given by the AnalysisPreparationTool */
  const MuonContainer* preselectedMuonTES = m_analysisPreparationTool->selectedMuons();
  if ( !preselectedMuonTES ) {
    mLog << MSG::WARNING << "Selected Muons Not Found " << endreq;
  }
  mLog << MSG::DEBUG << "Pre-selected Muons successfully retrieved - size is " << preselectedMuonTES->size() << " muons " << endreq;

  /** Check if the leading Electron and the Leadign Muon overlap or not */
  double deltaR = -1.0;
  if ( preselectedElecTES->size() > 0 && preselectedMuonTES->size() > 0) {
     const Electron * leadingElectron = preselectedElecTES->at(0);
     const Muon     * leadingMuon     = preselectedMuonTES->at(0);
     bool areOverlapping = m_analysisOverlapCheckingTool->overlap( leadingElectron, leadingMuon, deltaR ); 
     if ( areOverlapping ) mLog << MSG::INFO << "Leading Electron and Leading Muon overlap - deltaR = " << endreq;
  }

  /** now remove the overlaps using this tool 
      The input to the tool is the collection of pre-selected obeject obtained from the AnalysisPreparationTool 
      The output is various collections of final state non-overlapping particles 
      You can change the order of the overlap removal by change the input to the tool in job options */
  sc = m_analysisOverlapRemovalTool->execute();
  if ( sc.isFailure() ) {
    mLog << MSG::WARNING << "AnalysisPreparation Failed - overlap removal " << endreq;
    return StatusCode::SUCCESS;
  }

  /** get the final state Electrons - given by the AnalysisOverlapRemovalTool */
  const ElectronContainer* finalStateElecTES = m_analysisOverlapRemovalTool->finalStateElectrons();
  if ( !finalStateElecTES ) {

    mLog << MSG::WARNING << "Final State Electrons Not Found " << endreq;
  }
  mLog << MSG::DEBUG << "Final State Electrons successfully retrieved - size is " << finalStateElecTES->size() << " electrons " << endreq;




  mLog << MSG::DEBUG << "AnalysisPreparation() succeeded" << endreq;
 
  return StatusCode::SUCCESS;

}
//////////////////////////////////////////////////////////////////////////////////
/// Method to look at bjetInfo - called by execute() on event by event
/// Use jets after overlap removal and look at the b-tagging weights within
///
//////////////////////////////////////////////////////////////////////////////////
StatusCode AnalysisSkeleton::bjetInfo() {

  MsgStream mLog( messageService(), name() );

  mLog << MSG::DEBUG << "in bjetInfo()" << endreq;

  /** loop over jet container after overlap removal 
      As a check first get the container after selection cuts, but BEFORE Overlap Removal */

  const JetCollection* selectedJetTES = m_analysisPreparationTool->selectedJets();
  if ( !selectedJetTES ) {
    mLog << MSG::WARNING << "Selected Particle Jets Not Found " << endreq;
  }
  else{
    mLog << MSG::DEBUG << "Selected Jets successfully retrieved - size is " << selectedJetTES->size() << " jets " << endreq;}
  // 
  const JetCollection* finalStateJetTES = m_analysisOverlapRemovalTool->finalStateJets();
  if ( !finalStateJetTES ) {
    mLog << MSG::WARNING << "Final State Particle Jets Not Found " << endreq;
    return StatusCode::SUCCESS;
  }
  mLog << MSG::DEBUG << "Final State Jets successfully retrieved - size is " << finalStateJetTES->size() << " jets " << endreq;

  /** now look at some variables before and after overlap removal */

  int iflav;
  JetCollection::const_iterator jetItr_sel  = selectedJetTES->begin();
  JetCollection::const_iterator jetItrE_sel = selectedJetTES->end();
  for (; jetItr_sel != jetItrE_sel; ++jetItr_sel) {

    HepLorentzVector p4((*jetItr_sel)->px(),
         (*jetItr_sel)->py(),
         (*jetItr_sel)->pz(),
         (*jetItr_sel)->e());

    m_h_jet_eta_beforeOR->Fill(p4.pseudoRapidity());
    m_h_jet_et_beforeOR->Fill(p4.et()/1000.);

    /** get b-tagging info */

    double w_cmb = (*jetItr_sel)->getFlavourTagWeight(); // weight for IP3DSV1
    m_h_jet_ip3dsv1Wt_beforeOR->Fill(w_cmb);

    /** get quark flavour that originates this jet */
    iflav = getQuarkJetFlavour(jetItr_sel);
    m_h_jet_label_beforeOR->Fill((float) iflav);

    if(p4.et() > m_bjet_etCut && fabs(p4.pseudoRapidity()) < m_bjet_etaCut) {
      if(iflav==5) m_h_jet_ip3dsv1Wt_bjet_beforeOR->Fill(w_cmb);
      if(iflav==0) m_h_jet_ip3dsv1Wt_ujet_beforeOR->Fill(w_cmb);
    }

  }
  /** after overlapRemoval */
  JetCollection::const_iterator jetItr_fin  = finalStateJetTES->begin();
  JetCollection::const_iterator jetItrE_fin = finalStateJetTES->end();
  for (; jetItr_fin != jetItrE_fin; ++jetItr_fin) {

    HepLorentzVector p4((*jetItr_fin)->px(),
         (*jetItr_fin)->py(),
         (*jetItr_fin)->pz(),
         (*jetItr_fin)->e());

    m_h_jet_eta_afterOR->Fill(p4.pseudoRapidity());
    m_h_jet_et_afterOR->Fill(p4.et()/1000.);

    /** get b-tagging info */



    if(p4.et() > m_bjet_etCut && fabs(p4.pseudoRapidity()) < m_bjet_etaCut) {

      double w_cmb = (*jetItr_fin)->getFlavourTagWeight(); // weight for IP3DSV1
      m_h_jet_ip3dsv1Wt_afterOR->Fill(w_cmb);

      /** get quark flavour that originates this jet */
      iflav = getQuarkJetFlavour(jetItr_fin);
      m_h_jet_label_afterOR->Fill((float) iflav);
      //

      if(w_cmb > m_bjetWt_ip3dsv1Cut) m_aan_nbjets++; // count # of bjets in event

      if(iflav==5) m_h_jet_ip3dsv1Wt_bjet_afterOR->Fill(w_cmb);
      if(iflav==0) m_h_jet_ip3dsv1Wt_ujet_afterOR->Fill(w_cmb);
    }

  }

  return StatusCode::SUCCESS;

}
////////

int AnalysisSkeleton::getQuarkJetFlavour(JetCollection::const_iterator jetItr) {

    MsgStream mLog( messageService(), name() );

    /** flavour of quark that originated this jet */
    // --- get the true label of the jet from MC Truth

    std::string label("N/A");
    
    // FIXME - to be done by Vivek
    //const Analysis::TruthInfo* mcinfo = 0; // (*jetItr)->tagInfo<Analysis::TruthInfo>("TruthInfo");
    //if(mcinfo) {
    //  label = mcinfo->jetTruthLabel();
    //} else {
    //  mLog << MSG::VERBOSE << "could not find TruthInfo for matching jet" << endreq;
    //}
    int iflav(0);
    if(label=="B") {
      return iflav = 5;
    }
    if(label=="C") {
      return iflav = 4;
    }
    if(label=="T") {
      return iflav = 15;
    }

    return iflav;
}
////////////////////////////////////////////////////////////////////////////////////////////////
/// missing Et object

StatusCode AnalysisSkeleton::getMissingET() {

  MsgStream mLog(messageService(), name());
  mLog << MSG::DEBUG << "in getMissingEt()" << endreq;

  StatusCode sc = StatusCode::SUCCESS;

  if (!m_isAtlFastData) {
    /// retrieve the missing Et object from TDS
    sc = m_storeGate->retrieve(m_pMissing,m_missingETObjectName);
    if (sc.isFailure()) {
       mLog << MSG::ERROR << "Could not retrieve MissingET Object" << endreq;
       return StatusCode::SUCCESS;
    }
    else{ mLog<< MSG::DEBUG <<" retreived missing ET from AOD"<<endreq;}

    m_pxMiss = m_pMissing->etx();
    m_pyMiss = m_pMissing->ety();
    m_ptMiss = m_pMissing->et();
  } else {
    /// retrieve the missing Et object from TDS
    sc=m_storeGate->retrieve(m_pMissing, "AtlfastMisingEt");
    if( sc.isFailure()  ||  !m_pMissing ) {
      mLog << MSG::WARNING
      << "No Atlfast missing Et object found in TDS"
      << endreq; 
      return StatusCode::SUCCESS;
    }  
    m_pxMiss = m_pMissing->etx();
    m_pyMiss = m_pMissing->ety();
    m_ptMiss = m_pMissing->et();
  }

  /// fill missing energy histograms
  m_pxMis->Fill(m_pxMiss);
  m_pyMis->Fill(m_pyMiss);
  m_ptMis->Fill(m_ptMiss);

  m_aan_ptMiss = m_ptMiss;

  return sc;

}
////////////////////////////////////////////////////////
StatusCode AnalysisSkeleton::SusyStudies() {

  /** Make some introductory plots */

  MsgStream mLog( messageService(), name() );

  mLog << MSG::DEBUG << "in SusyStudies()" << endreq;

  /** loop over truth container and get pT of the earliest top quarks */

  double pTtop1;
  double pTtop2;
  int numTops;

  StatusCode sc = getTopQpT(numTops, pTtop1, pTtop2);  
  if(!sc) mLog << MSG::DEBUG << "Something wrong with finding top quark pT"<<endreq;
  else{

    m_aan_NumTopQ = 2;
    m_aan_pTtop1 = pTtop1;
    m_aan_pTtop2 = pTtop2;
  }

  /** loop over jet container after overlap removal */

  const JetCollection* finalStateJetTES = m_analysisOverlapRemovalTool->finalStateJets();
  if ( !finalStateJetTES ) {
    mLog << MSG::WARNING << "Final State Particle Jets Not Found " << endreq;
    return StatusCode::SUCCESS;
  }
  mLog << MSG::DEBUG << "Final State Jets successfully retrieved - size is " << finalStateJetTES->size() << " jets " << endreq;

  double w_cmb=-100;
  JetCollection::const_iterator jetItr_fin  = finalStateJetTES->begin();
  JetCollection::const_iterator jetItrE_fin = finalStateJetTES->end();
  for (; jetItr_fin != jetItrE_fin; ++jetItr_fin) {

    HepLorentzVector p4((*jetItr_fin)->px(),
         (*jetItr_fin)->py(),
         (*jetItr_fin)->pz(),
         (*jetItr_fin)->e());

    /** variables for the ntuple  */
    m_aan_njets++;
    if(fabs(p4.pseudoRapidity()) < m_bjet_etaCut ) m_aan_njets_etaLT25++;
    if(p4.et()> m_SusyJetMinEt ) m_aan_njets_SusyETCut++;

    m_aan_JetEta->push_back((*jetItr_fin)->eta());
    m_aan_JetEt->push_back(p4.et());

    w_cmb = -100;
    if(p4.et() > m_bjet_etCut && fabs(p4.pseudoRapidity()) < m_bjet_etaCut) {
      w_cmb = (*jetItr_fin)->getFlavourTagWeight();} // weight for IP3DSV1
    m_aan_JetBTagWt->push_back(w_cmb);

    m_aan_ht += p4.et(); // scalar sum of jet ET
    if(p4.et()>m_aan_maxJetET) m_aan_maxJetET = p4.et(); // Jet with max ET

  } // loop over jets

  /** Get final state leptons. We need these in the determination of effective mass
      We should also store the leptons in the ntuple, so that we can re-do the calculation later */
  
  /** First do electrons */

  const ElectronContainer* finalStateElecTES = m_analysisOverlapRemovalTool->finalStateElectrons();
  if ( !finalStateElecTES ) {

    mLog << MSG::WARNING << "SusyStudies: Final State Electrons Not Found " << endreq;
    return StatusCode::SUCCESS;
  }
  mLog << MSG::DEBUG << "SusyStudies: Final State Electrons successfully retrieved - size is " << finalStateElecTES->size() << " electrons " << endreq;


  /** iterators over the electron container - final state */ 
  ElectronContainer::const_iterator felecItr  = finalStateElecTES->begin();
  ElectronContainer::const_iterator felecItrE = finalStateElecTES->end();

  double sum_elET=0;

  for (; felecItr != felecItrE; ++felecItr) {

    /** apply further selections if necessary */
    /** check for the author of the electron */
    bool author = (*felecItr)->author(egammaParameters::AuthorElectron);
    if ( !author || (*felecItr)->pt() < m_etElecCut ) continue;


    const EMShower* eshow = (*felecItr)->detail<EMShower>("egDetailAOD");
    double etisol = -1;
    if( eshow ) etisol = eshow->parameter(egammaParameters::etcone20);

    mLog << MSG::DEBUG << "SusyStudies:isEM/etisol "<< (*felecItr)->isem()<<","<<etisol<<endreq;

    m_aan_NFinalEl++;
    m_aan_FinalElPt->push_back((*felecItr)->pt());
    m_aan_FinalElEta->push_back((*felecItr)->eta());
    m_aan_FinalElEtCone20->push_back(etisol);
    
    /** isEM cut already picks out isolated electrons, so this is just a sanity check */
    if(etisol/1000.<10) sum_elET += (*felecItr)->pt();
  }


  /** Now look at muons */

  const MuonContainer* finalStateMuonTES = m_analysisOverlapRemovalTool->finalStateMuons();
  if ( !finalStateMuonTES ) {

    mLog << MSG::WARNING << "SusyStudies: Final State Muons Not Found " << endreq;
    return StatusCode::SUCCESS;
  }
  mLog << MSG::DEBUG << "SusyStudies: Final State Muons successfully retrieved - size is " << finalStateMuonTES->size() << " muons " << endreq;


  /** iterators over the muon container - final state  */ 
  MuonContainer::const_iterator fmuonItr  = finalStateMuonTES->begin();
  MuonContainer::const_iterator fmuonItrE = finalStateMuonTES->end();

  double sum_muET = 0;

  for (; fmuonItr != fmuonItrE; ++fmuonItr) {

    /** apply further selections if necessary */
    /** check for the author of the muon */

    double etIsol = (*fmuonItr)->parameter( static_cast<MuonParameters::ParamDef>(1) ); // dR of 0.2

    mLog << MSG::DEBUG << "SusyStudies:Muon etisol/best match "<< etIsol<<","<<(*fmuonItr)->bestMatch()<<endreq;
   
    m_aan_NFinalMu++;
    m_aan_FinalMuPt->push_back( (*fmuonItr)->pt());
    m_aan_FinalMuEta->push_back( (*fmuonItr)->eta());
    m_aan_FinalMuEtCone20->push_back( etIsol);
    m_aan_FinalMuBestMat->push_back( (*fmuonItr)->bestMatch());
    m_aan_FinalMuMatChi2->push_back((*fmuonItr)->matchChi2());

    /** require bestMatch, chi2 and isolation cuts */
    if((*fmuonItr)->bestMatch()==1 && (*fmuonItr)->matchChi2() <100. && etIsol/1000. < 10) sum_muET +=(*fmuonItr)->pt();

  }


  /** now calculate effmass */
  m_aan_FinalLepEtSum = sum_elET + sum_muET; // keep leptons separate for now.
  m_aan_FinalElEtSum  = sum_elET;
  m_aan_FinalMuEtSum  = sum_muET;

  m_aan_effmass = m_aan_ptMiss + m_aan_ht; // scalar sum of jet ET + missing ET 

  return StatusCode::SUCCESS;

}
StatusCode AnalysisSkeleton::getTopQpT(int& numTops, double& top1, double& top2) {

  MsgStream mLog( messageService(), name() );

  mLog << MSG::DEBUG << "in getTopQpT()" << endreq;

  //
  top1 = -1; top2 = -1;
  //
  double topPt[2];
  topPt[0]=-1;
  topPt[1]=-1;
  //
  /** get the MC truth particle AOD container from StoreGate */

  const TruthParticleContainer*  mcpartTES = 0;
  StatusCode sc=m_storeGate->retrieve( mcpartTES, m_truthParticleContainerName);
  if( sc.isFailure()  ||  !mcpartTES ) {
     mLog << MSG::WARNING
          << "No AOD MC truth particle container found in TDS"
          << endreq; 
     return StatusCode::SUCCESS;
  }
  mLog <<MSG::DEBUG << "MC Truth Container Successfully Retrieved" << endreq;

  /** loop over particles and get the top quarks produced at the hard scatter */

  TruthParticleContainer::const_iterator mcpItr  = mcpartTES->begin();
  TruthParticleContainer::const_iterator mcpItrE = mcpartTES->end();

  numTops=0;
  for (; mcpItr != mcpItrE; ++mcpItr) {

    const HepMC::GenParticle* part =(*mcpItr)->genParticle();
    int pdgid = part->pdg_id();

    if(numTops==2) break; // quit if we have two already

    if(abs(pdgid)==6) { // it is a top

      HepMC::GenVertex* prod_vtx = part->production_vertex();
      int vtx_barcode = 1;
      if(prod_vtx) vtx_barcode = prod_vtx->barcode();

      if(vtx_barcode == -1) {

        topPt[numTops] = (part->momentum()).perp();
        numTops++;

      }
    

    }
  }

  top1 = topPt[0];
  top2 = topPt[1];
  return StatusCode::SUCCESS;
}
Now the AnalysisSkeleton_topOptions.py

Python bindings file

from AthenaCommon.AppMgr import ToolSvc

from AthenaCommon.AppMgr import ServiceMgr
# Event selector
import AthenaPoolCnvSvc.ReadAthenaPool
#EventSelector.BackNavigation = True

# Particle Properties
from PartPropSvc.PartPropSvcConf import PartPropSvc

# the POOL converters
include( "ParticleBuilderOptions/ESD_PoolCnv_jobOptions.py" )
include( "ParticleBuilderOptions/AOD_PoolCnv_jobOptions.py")
include( "ParticleBuilderOptions/McAOD_PoolCnv_jobOptions.py")
include( "EventAthenaPool/EventAthenaPool_joboptions.py" )

# The AOD input file
#ServiceMgr.EventSelector.InputCollections = [ "dcache:AOD.pool.root.1" ]
#ServiceMgr.EventSelector.InputCollections = [ "AOD.pool.root" ]
ServiceMgr.EventSelector.InputCollections = [ "/tmp/quinonez/mc08.105200.T1_McAtNlo_Jimmy.recon.AOD.e357_s462_r541/AOD.026357._00491.pool.root.1" ]

# Get the selection, overlap checking and overlap removal tools  
include ( "UserAnalysisUtils/UserAnalysisSelectionTool_jobOptions.py" )
include ( "UserAnalysisUtils/UserAnalysisPreparationTool_jobOptions.py" )
include ( "UserAnalysisUtils/UserAnalysisOverlapCheckingTool_jobOptions.py" )
include ( "UserAnalysisUtils/UserAnalysisOverlapRemovalTool_jobOptions.py" )

# Athena-Aware NTuple making Tools
CBNTAthenaAware = True
include ("CBNT_Athena/CBNT_AthenaAware_jobOptions.py")
include ("CBNT_Athena/CBNT_EventInfo_jobOptions.py")

# list of the algorithms to be executed at run time
from UserAnalysis.UserAnalysisConf import AnalysisSkeleton
topSequence.CBNT_AthenaAware += AnalysisSkeleton() 
AnalysisSkeleton = AnalysisSkeleton()

############# The properties of the AnalysisSkeleton Algorithm
AnalysisSkeleton.AnalysisSelectionTool       = ToolSvc.UserAnalysisSelectionTool
AnalysisSkeleton.AnalysisPreparationTool     = ToolSvc.UserAnalysisPreparationTool
AnalysisSkeleton.AnalysisOverlapCheckingTool = ToolSvc.UserAnalysisOverlapCheckingTool
AnalysisSkeleton.AnalysisOverlapRemovalTool  = ToolSvc.UserAnalysisOverlapRemovalTool

IsAtlfast = False

AnalysisSkeleton.McParticleContainer = "SpclMC"
AnalysisSkeleton.ElectronContainer = "ElectronAODCollection"
AnalysisSkeleton.MissingETObject = "MET_RefFinal"
AnalysisSkeleton.DeltaRMatchCut = 0.2
AnalysisSkeleton.MaxDeltaR = 0.9999
AnalysisSkeleton.ElectronEtCut  = 10.0*GeV
AnalysisSkeleton.ElectronEtaCut = 2.5
AnalysisSkeleton.ElectronCone   = 0.9
AnalysisSkeleton.bjetWt_IP3DSV1Cut = 6
AnalysisSkeleton.bjet_etaCut = 2.5
AnalysisSkeleton.bjet_etCut = 15.0*GeV
AnalysisSkeleton.MissingETCut = 20.0*GeV
AnalysisSkeleton.OutputLevel = INFO
AnalysisSkeleton.IsAtlFastData = IsAtlfast
AnalysisSkeleton.SusyJetMinEt      = 50*GeV

# Change the selections if necesary
# VJ: Feb 28, 2008:
# At the moment, we cannot seem to propagate new settings back to
# the various tools. So, you will have to edit the jO files in
# UserAnalysisUtils/ package
#
ToolSvc.UserAnalysisSelectionTool.IsAtlfastData = IsAtlfast
#AnalysisSkeleton.AnalysisSelectionTool.ElectronIsEMFlag="Loose"
ToolSvc.UserAnalysisSelectionTool.MuonPt=6.0*GeV
ToolSvc.UserAnalysisSelectionTool.JetPt=20.0*GeV

# configure the overlap checking tool
ToolSvc.UserAnalysisOverlapCheckingTool.OverlapDeltaR=0.2
ToolSvc.UserAnalysisOverlapCheckingTool.OverlapDeltaRWithJets=0.3

# Building the  containers of selected obejcts
ToolSvc.UserAnalysisPreparationTool.IsAtlfastData = IsAtlfast

#input cntainer keys to the pre-selection tool
ToolSvc.UserAnalysisPreparationTool.InputContainerKeys = [
    "ElectronAODCollection",
    "StacoMuonCollection",
    "TauRecContainer",
    "Cone4H1TopoJets",
    "PhotonAODCollection",
    "CaloCalTopoCluster",
    "TrackParticleCandidate"
   ]
# Output container keys after the pre-selections
ToolSvc.UserAnalysisPreparationTool.OutputContainerKeys=[ 
    "SelectedElectronCollection",
    "SelectedStacoMuonCollection",
    "SelectedTauRecContainer",
    "SelectedCone4H1TowerJets",
    "SelectedPhotonAODCollection",
    "SelectedCaloCalTopoCluster",
    "SelectedTrackParticleCandidate"
   ]

# Use the output containers fromi the selection tool as input to the overalp removal tool
ToolSvc.UserAnalysisOverlapRemovalTool.InputContainerKeys=[  
    "SelectedElectronCollection",
    "SelectedStacoMuonCollection",
    "SelectedTauRecContainer",
    "SelectedCone4H1TowerJets",
    "SelectedPhotonAODCollection",
    "SelectedCaloCalTopoCluster",
    "SelectedTrackParticleCandidate"
   ]

# The output container keys after the overlap-removal
# Note that a container of all leptons is provided on output
# as well as a container of all final state particles
AnalysisSkeleton.AnalysisOverlapRemovalTool.IsAtlfastData          = IsAtlfast
AnalysisSkeleton.AnalysisOverlapRemovalTool.OuputObjectKey         = "FinalStateObjectCollection"
AnalysisSkeleton.AnalysisOverlapRemovalTool.OutputLeptonKey        = "FinalStateLeptonCollection"
AnalysisSkeleton.AnalysisOverlapRemovalTool.OutputPhotonKey        = "FinalStatePhotonCollection"
AnalysisSkeleton.AnalysisOverlapRemovalTool.OutputElectronKey      = "FinalStateElectronCollection"
AnalysisSkeleton.AnalysisOverlapRemovalTool.OutputMuonKey          = "FinalStateMuonCollection"
AnalysisSkeleton.AnalysisOverlapRemovalTool.OutputTauJetKey        = "FinalStateTauJetCollection"
AnalysisSkeleton.AnalysisOverlapRemovalTool.OutputCalloClusterKey  = "FinalStateCaloClusterCollection"
AnalysisSkeleton.AnalysisOverlapRemovalTool.OutputTrackParticleKey = "FinalStateTrackParticleCollection"
AnalysisSkeleton.AnalysisOverlapRemovalTool.OutputJetKey           = "FinalStateJetCollection"
AnalysisSkeleton.AnalysisOverlapRemovalTool.OutputBJetKey          = "FinalStateBJetCollection"
AnalysisSkeleton.AnalysisOverlapRemovalTool.OutputLightJetKey      = "FinalStateLightJetCollection"

##
##########################################
## Set up trigger configuration service and metadata service is relies on, for analysis job without RecExCommon

AnalysisSkeleton.DoTrigger = True

if AnalysisSkeleton.DoTrigger:
   from AthenaCommon.GlobalFlags import GlobalFlags
   GlobalFlags.DetGeo.set_atlas()

   ## set up trigger decision tool
   from TrigDecision.TrigDecisionConfig import TrigDecisionTool
   tdt = TrigDecisionTool()
   ToolSvc += tdt

   # might be needed for TriggerConfigGetter...
   from RecExConfig.RecFlags  import rec
   rec.readAOD=True
   rec.doWriteAOD=False
   rec.doWriteESD=False

   # To read files with trigger config stored as in-file meta-data,
   # i.e. 13.0.40 and above: ds
   # To read AOD produced with 13.0.30 you need to change ds to aod: 
   from TriggerJobOpts.TriggerFlags import TriggerFlags
   TriggerFlags.configurationSourceList = ['ds']

   # set up trigger config service
   from TriggerJobOpts.TriggerConfigGetter import TriggerConfigGetter
   cfg =  TriggerConfigGetter()

## END of trigger setup

##########################################
# setup TTree registration Service
# save ROOT histograms and Tuple
from GaudiSvc.GaudiSvcConf import THistSvc
ServiceMgr += THistSvc()
ServiceMgr.THistSvc.Output = [ "AANT DATAFILE='AnalysisSkeleton.aan.root' OPT='RECREATE'" ]
from AnalysisTools.AnalysisToolsConf import AANTupleStream
topSequence += AANTupleStream()
AANTupleStream = AANTupleStream()
AANTupleStream.ExtraRefNames = [ "StreamESD","Stream1" ]
AANTupleStream.OutputName = 'AnalysisSkeleton.aan.root'
AANTupleStream.WriteInputDataHeader = True
AANTupleStream.OutputLevel = WARNING

# Set output level threshold (2=DEBUG, 3=INFO, 4=WARNING, 5=ERROR, 6=FATAL )
ServiceMgr.MessageSvc.OutputLevel = INFO

# Number of Events to process
theApp.EvtMax = -1
#theApp.EvtMax = 10

###################### For interactive analysis
#include ("PyAnalysisCore/InitPyAnalysisCore.py")

from GaudiSvc.GaudiSvcConf import AuditorSvc
ServiceMgr.AuditorSvc.Auditors  += [ "ChronoAuditor"]

AthenaPoolCnvSvc = Service("AthenaPoolCnvSvc")
AthenaPoolCnvSvc.UseDetailChronoStat = TRUE

QUESTIONS

SHORT-TERM

  • How can I see the pythia output of my AOD?

MEDIUM TERM

  • 10000 events are statistically down?

LONG-TERM

  • Displaced vertices

-- FernandoAndresQuinonezGranados - 30 Apr 2009

Edit | Attach | Watch | Print version | History: r4 < r3 < r2 < r1 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r4 - 2009-05-04 - FernandoAndresQuinonezGranados
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    Main All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright & 2008-2019 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback