diff --git a/Generators/HepMCAnalysis_i/CMakeLists.txt b/Generators/HepMCAnalysis_i/CMakeLists.txt
index 65867868c356622011ff25a2950ec5f3e0684630..abbed775bb511f15534b554be57e941e7bcb3cd5 100644
--- a/Generators/HepMCAnalysis_i/CMakeLists.txt
+++ b/Generators/HepMCAnalysis_i/CMakeLists.txt
@@ -23,6 +23,7 @@ find_package( HEPUtils )
find_package( HepMCAnalysis )
find_package( ROOT COMPONENTS Core MathCore Hist RIO )
find_package( FastJet )
+find_package( HepMC )
# Component(s) in the package:
atlas_add_component( HepMCAnalysis_i
@@ -30,6 +31,7 @@ atlas_add_component( HepMCAnalysis_i
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${HEPMCANALYSIS_INCLUDE_DIRS}
${HEPUTILS_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS}
${FASTJET_INCLUDE_DIRS}
+ ${HEPMC_INCLUDE_DIRS}
LINK_LIBRARIES ${ROOT_LIBRARIES} ${HEPMCANALYSIS_LIBRARIES}
${HEPUTILS_LIBRARIES} ${CLHEP_LIBRARIES} AtlasHepMCLib
${FASTJET_LIBRARIES} AthenaBaseComps GaudiKernel StoreGateLib EventInfo
diff --git a/Generators/HepMCAnalysis_i/HepMCAnalysis_i/HepMCAnalysis.h b/Generators/HepMCAnalysis_i/HepMCAnalysis_i/HepMCAnalysis.h
index 114515b787949cb8713f62e18312383ce1f8349f..7b23626d5b23d03c14a9433f967f6e47e6330e3f 100644
--- a/Generators/HepMCAnalysis_i/HepMCAnalysis_i/HepMCAnalysis.h
+++ b/Generators/HepMCAnalysis_i/HepMCAnalysis_i/HepMCAnalysis.h
@@ -19,7 +19,6 @@
#include "GaudiKernel/ServiceHandle.h"
#include "TH1.h"
-#include "AtlasHepMC/GenEvent.h"
// forward declarations
class ITHistSvc;
@@ -66,8 +65,6 @@ class HepMCAnalysis
std::string m_key;
std::string m_infokey;
- //convert unit MeV to GeV for energy and momenta
- void MeVToGeV (HepMC::GenEvent* evt);
bool m_JetFinder;
bool m_JetAnalysis;
diff --git a/Generators/HepMCAnalysis_i/HepMCAnalysis_i/LeptonJetAnalysis.h b/Generators/HepMCAnalysis_i/HepMCAnalysis_i/LeptonJetAnalysis.h
index 6c1f44dffacc25e1beeb46f8afbe9b1373bbc6d0..50bd4735a5a06edc345d87c8950d4650b437b08b 100644
--- a/Generators/HepMCAnalysis_i/HepMCAnalysis_i/LeptonJetAnalysis.h
+++ b/Generators/HepMCAnalysis_i/HepMCAnalysis_i/LeptonJetAnalysis.h
@@ -7,8 +7,6 @@
#include "baseAnalysis.h"
-// forward declarations
-#include "AtlasHepMC/GenEvent_fwd.h"
class TH1D;
class LeptonJetAnalysis: public baseAnalysis
diff --git a/Generators/HepMCAnalysis_i/HepMCAnalysis_i/ParticleContentAnalysis.h b/Generators/HepMCAnalysis_i/HepMCAnalysis_i/ParticleContentAnalysis.h
index 45cfd6ce142dd8fa01990f53c9eb7e09c44090d9..8cb78563a7972969de7f150aabb54020041b75e5 100644
--- a/Generators/HepMCAnalysis_i/HepMCAnalysis_i/ParticleContentAnalysis.h
+++ b/Generators/HepMCAnalysis_i/HepMCAnalysis_i/ParticleContentAnalysis.h
@@ -7,11 +7,6 @@
#include "baseAnalysis.h"
-// forward declarations
-// forward declarations
-#include "AtlasHepMC/GenEvent_fwd.h"
-#include "AtlasHepMC/GenParticle_fwd.h"
-
class TH1D;
/*
diff --git a/Generators/HepMCAnalysis_i/HepMCAnalysis_i/PdfAnalysis.h b/Generators/HepMCAnalysis_i/HepMCAnalysis_i/PdfAnalysis.h
index 2cecf588b5a6f99210325a08397830465bc1186a..95a20911cfb4d59de98aaf26c8fc699cfb815e82 100644
--- a/Generators/HepMCAnalysis_i/HepMCAnalysis_i/PdfAnalysis.h
+++ b/Generators/HepMCAnalysis_i/HepMCAnalysis_i/PdfAnalysis.h
@@ -7,10 +7,6 @@
#include "baseAnalysis.h"
-// forward declarations
-#include "AtlasHepMC/GenEvent_fwd.h"
-#include "AtlasHepMC/PdfInfo.h"
-
class TH1D;
class PdfAnalysis: public baseAnalysis
diff --git a/Generators/HepMCAnalysis_i/HepMCAnalysis_i/UserAnalysis.h b/Generators/HepMCAnalysis_i/HepMCAnalysis_i/UserAnalysis.h
index a81836c552559fac4febb8826c0019a890ddbcc2..112850fccf5d528c12602e3a26f8ce53046a9c39 100644
--- a/Generators/HepMCAnalysis_i/HepMCAnalysis_i/UserAnalysis.h
+++ b/Generators/HepMCAnalysis_i/HepMCAnalysis_i/UserAnalysis.h
@@ -7,8 +7,6 @@
#include "baseAnalysis.h"
-// forward declarations
-#include "AtlasHepMC/GenEvent_fwd.h"
class TH1D;
class UserAnalysis: public baseAnalysis
diff --git a/Generators/HepMCAnalysis_i/src/GetEvents.cxx b/Generators/HepMCAnalysis_i/src/GetEvents.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..693bc84995c1ea53b51b49f48b88a6580fbffd03
--- /dev/null
+++ b/Generators/HepMCAnalysis_i/src/GetEvents.cxx
@@ -0,0 +1,55 @@
+/**
+ * Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ * Author: Andrii Verbytsky andrii.verbytskyi@mpp.mpg.de
+ *
+ * This file implements helper functions that prepare the ATLAS MC events
+ * to be used with the HepMCAnalysis package.
+ */
+#include <vector>
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "StoreGate/StoreGateSvc.h"
+#include "GaudiKernel/ITHistSvc.h"
+#include "GaudiKernel/Property.h"
+
+#include "EventInfo/EventInfo.h"
+#include "EventInfo/EventID.h"
+#include "GeneratorObjects/McEventCollection.h"
+
+//class HepMCAnalysisGenEvent;
+using HepMCAnalysisGenEvent=HepMC::GenEvent;
+HepMCAnalysisGenEvent* PrepareHepMCAnalysisGenEvent(const HepMC::GenEvent* cevent);
+
+StatusCode GetRunEventNumber(AthAlgorithm* a,int& runNumber,int& evtNumber, const std::string im)
+{
+ // load event info
+ const EventInfo* mcInfoptr;
+ if ( a->evtStore()->retrieve( mcInfoptr, im ).isFailure() ) {
+ //ATH_MSG_ERROR( "Could not retrieve EventInfo" );
+ return StatusCode::FAILURE;
+
+ } else{
+ runNumber = mcInfoptr->event_ID()->run_number();
+ evtNumber = mcInfoptr->event_ID()->event_number();
+ //ATH_MSG_DEBUG( "run: " << runNumber << " event: " << evtNumber );
+ //ATH_MSG_DEBUG("Event info loaded");
+ }
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode GetEvents(AthAlgorithm* a,std::vector<HepMCAnalysisGenEvent*>& evts, const std::string km)
+{
+ // load HepMC info
+ const McEventCollection* mcCollptr;
+ if ( a->evtStore()->retrieve( mcCollptr, km ).isFailure() ) {
+ //ATH_MSG_WARNING( "Could not retrieve McEventCollection" );
+ } else {
+ //ATH_MSG_DEBUG("HepMC info loaded");
+
+ // loop over all events in McEventCollection
+ for ( McEventCollection::const_iterator itr = mcCollptr->begin(); itr != mcCollptr->end(); ++itr ) {
+ evts.push_back(PrepareHepMCAnalysisGenEvent(*itr));
+}
+}
+ return StatusCode::SUCCESS;
+}
\ No newline at end of file
diff --git a/Generators/HepMCAnalysis_i/src/HepMCAnalysis.cxx b/Generators/HepMCAnalysis_i/src/HepMCAnalysis.cxx
index 30d84aad66f36f7ff56ef5e871f77714b0791c38..1b4b25141e91424ae390b05d021f69b53a0184b1 100644
--- a/Generators/HepMCAnalysis_i/src/HepMCAnalysis.cxx
+++ b/Generators/HepMCAnalysis_i/src/HepMCAnalysis.cxx
@@ -13,12 +13,7 @@
#include "GaudiKernel/ITHistSvc.h"
#include "GaudiKernel/Property.h"
-#include "EventInfo/EventInfo.h"
-#include "EventInfo/EventID.h"
-#include "GeneratorObjects/McEventCollection.h"
-
-#include "AtlasHepMC/GenEvent.h"
-
+// ROOT
#include "TFile.h"
#include "TH1.h"
#include "TString.h"
@@ -45,13 +40,10 @@
#include "HepMCAnalysis_i/ParticleContentAnalysis.h"
#include "HepMCAnalysis_i/PdfAnalysis.h"
-using std::cout;
-using std::endl;
-using std::string;
-using std::vector;
+
// ----------------------------------------------------------------------
-HepMCAnalysis::HepMCAnalysis(const string &name, ISvcLocator *pSvcLocator):
+HepMCAnalysis::HepMCAnalysis(const std::string &name, ISvcLocator *pSvcLocator):
AthAlgorithm(name, pSvcLocator), m_histSvc("ITHistSvc", name)
{
declareProperty("McEventKey", m_key = "GEN_EVENT");
@@ -95,119 +87,106 @@ HepMCAnalysis::~HepMCAnalysis()
}
-// ----------------------------------------------------------------------
-void HepMCAnalysis::MeVToGeV(HepMC::GenEvent *event)
-{
- for (HepMC::GenEvent::particle_iterator p = event->particles_begin(); p != event->particles_end(); ++p) {
- const HepMC::FourVector fv((*p)->momentum().px() / 1000.0,
- (*p)->momentum().py() / 1000.0,
- (*p)->momentum().pz() / 1000.0,
- (*p)->momentum().e() / 1000.0);
-
- (*p)->set_momentum(fv);
- }
-}
-
// ----------------------------------------------------------------------
StatusCode HepMCAnalysis::initialize()
{
ATH_MSG_DEBUG("Initializing " << name() << "...");
if (m_JetFinder) {
- cout << "Adding module JetFinder" << endl;
+ std::cout << "Adding module JetFinder" << std::endl;
m_jetfinder.push_back(new JetFinder);
}
// create analyses and initialize them
if (m_JetAnalysis) {
- cout << "Adding module JetAnalysis" << endl;
+ std::cout << "Adding module JetAnalysis" << std::endl;
m_analysis.push_back(new JetAnalysis);
}
if (m_WplusJetAnalysis) {
- cout << "Adding module WplusJetAnalysis" << endl;
+ std::cout << "Adding module WplusJetAnalysis" << std::endl;
m_analysis.push_back(new WplusJetAnalysis);
}
if (m_ZAnalysis) {
- cout << "Adding module ZAnalysis" << endl;
+ std::cout << "Adding module ZAnalysis" << std::endl;
m_analysis.push_back(new ZAnalysis);
}
if (m_ZtautauAnalysis) {
- cout << "Adding module ZtautauAnalysis" << endl;
+ std::cout << "Adding module ZtautauAnalysis" << std::endl;
m_analysis.push_back(new ZtautauAnalysis);
}
if (m_WtaunuAnalysis) {
- cout << "Adding module WtaunuAnalysis" << endl;
+ std::cout << "Adding module WtaunuAnalysis" << std::endl;
m_analysis.push_back(new WtaunuAnalysis);
}
if (m_ttbarAnalysis) {
- cout << "Adding module ttbarAnalysis" << endl;
+ std::cout << "Adding module ttbarAnalysis" << std::endl;
m_analysis.push_back(new ttbarAnalysis);
}
if (m_bbbarAnalysis) {
- cout << "Adding module bbbarAnalysis" << endl;
+ std::cout << "Adding module bbbarAnalysis" << std::endl;
m_analysis.push_back(new bbbarAnalysis);
}
if (m_UEAnalysis) {
- cout << "Adding module UEAnalysis" << endl;
+ std::cout << "Adding module UEAnalysis" << std::endl;
m_analysis.push_back(new UEAnalysis);
}
if (m_EtmissAnalysis) {
- cout << "Adding module EtMissAnalysis" << endl;
+ std::cout << "Adding module EtMissAnalysis" << std::endl;
m_analysis.push_back(new EtMissAnalysis);
}
if (m_ElasScatAnalysis) {
- cout << "Adding module ElasScatAnalysis" << endl;
+ std::cout << "Adding module ElasScatAnalysis" << std::endl;
m_analysis.push_back(new ElasScatAnalysis);
}
if (m_UserAnalysis) {
- cout << "Adding module UserAnalysis" << endl;
+ std::cout << "Adding module UserAnalysis" << std::endl;
m_analysis.push_back(new UserAnalysis);
}
ATH_MSG_DEBUG("Analyses added");
if (m_LeptonJetAnalysis) {
- cout << "Adding module LeptonJetAnalysis" << endl;
+ std::cout << "Adding module LeptonJetAnalysis" << std::endl;
m_analysis.push_back(new LeptonJetAnalysis);
}
ATH_MSG_DEBUG("Analyses added");
if (m_ParticleContentAnalysis) {
- cout << "Adding module ParticleContentAnalysis" << endl;
+ std::cout << "Adding module ParticleContentAnalysis" << std::endl;
m_analysis.push_back(new ParticleContentAnalysis);
}
ATH_MSG_DEBUG("Analyses added");
if (m_PdfAnalysis) {
- cout << "Adding module PdfAnalysis" << endl;
+ std::cout << "Adding module PdfAnalysis" << std::endl;
m_analysis.push_back(new PdfAnalysis);
}
ATH_MSG_DEBUG("Analyses added");
- for (vector<baseAnalysis*>::const_iterator i(m_jetfinder.begin()); i != m_jetfinder.end(); ++i) {
+ for (std::vector<baseAnalysis*>::const_iterator i(m_jetfinder.begin()); i != m_jetfinder.end(); ++i) {
(*i)->Init(m_max_eta,m_min_pt);
(*i)->InitJetFinder(m_jet_coneRadius,m_jet_overlapThreshold,m_jet_ptmin, m_lepton_ptmin, m_DeltaR_lepton_track);
}
- for (vector<baseAnalysis*>::const_iterator i(m_analysis.begin()); i != m_analysis.end(); ++i) {
+ for (std::vector<baseAnalysis*>::const_iterator i(m_analysis.begin()); i != m_analysis.end(); ++i) {
(*i)->Init(m_max_eta,m_min_pt);
(*i)->InitJetFinder(m_jet_coneRadius,m_jet_overlapThreshold,m_jet_ptmin, m_lepton_ptmin, m_DeltaR_lepton_track);
}
//TH1::SetDirectory(0) prevent histograms to be added to the list of histiograms of currently opened file (solve conflict with TestHepMC)
- for ( vector< baseAnalysis* >::const_iterator ana = m_analysis.begin(); ana != m_analysis.end(); ++ana )
+ for ( std::vector< baseAnalysis* >::const_iterator ana = m_analysis.begin(); ana != m_analysis.end(); ++ana )
{
- vector<TH1D*> hists = (*ana)->Histograms();
- for (vector<TH1D* >::iterator h = hists.begin(); h != hists.end(); h++)
+ std::vector<TH1D*> hists = (*ana)->Histograms();
+ for (std::vector<TH1D* >::iterator h = hists.begin(); h != hists.end(); h++)
(*h)->SetDirectory(0);
}
@@ -225,6 +204,8 @@ StatusCode HepMCAnalysis::initialize()
return StatusCode::SUCCESS;
}
+StatusCode GetRunEventNumber(AthAlgorithm* a,int& runNumber,int& evtNumber, const std::string im);
+StatusCode GetEvents( AthAlgorithm* a, std::vector<HepMC::GenEvent*>& evts, const std::string km);
// ----------------------------------------------------------------------
StatusCode HepMCAnalysis::execute()
{
@@ -232,71 +213,32 @@ StatusCode HepMCAnalysis::execute()
int runNumber = 0, evtNumber = 0;
- // load event info
- const EventInfo* mcInfoptr;
- if ( evtStore()->retrieve( mcInfoptr, m_infokey ).isFailure() ) {
- ATH_MSG_ERROR( "Could not retrieve EventInfo" );
- return StatusCode::FAILURE;
-
- } else{
- runNumber = mcInfoptr->event_ID()->run_number();
- evtNumber = mcInfoptr->event_ID()->event_number();
- ATH_MSG_DEBUG( "run: " << runNumber << " event: " << evtNumber );
- }
- ATH_MSG_DEBUG("Event info loaded");
-
- // load HepMC info
- const McEventCollection* mcCollptr;
- if ( evtStore()->retrieve( mcCollptr, m_key ).isFailure() ) {
- ATH_MSG_WARNING( "Could not retrieve McEventCollection" );
- return StatusCode::SUCCESS;
-
- } else {
- ATH_MSG_DEBUG("HepMC info loaded");
-
- // loop over all events in McEventCollection
- for ( McEventCollection::const_iterator itr = mcCollptr->begin(); itr != mcCollptr->end(); ++itr ) {
- // FIXME: This gets an object from SG, which we are not allowed to modify.
- // But we do modify it in MeVToGeV.
- // This needs to change.
- // Further, we need a non-const pointer to call the Process() and
- // ClearEvent() methods of baseAnalysis, although from a spot check,
- // it appears that these do not actually change the object.
- // If that is the case, then these interfaces should really be changed
- // to take a const pointer.
- HepMC::GenEvent *evt = const_cast<HepMC::GenEvent*>(*itr);
-
- // convert units
- MeVToGeV( evt );
-
+ StatusCode runeventnumber_statuscode =GetRunEventNumber(this,runNumber,evtNumber,m_infokey);
+ if (runeventnumber_statuscode==StatusCode::FAILURE) return StatusCode::FAILURE;
+ std::vector<HepMC::GenEvent*> evts;
+ StatusCode events_statuscode=GetEvents(this,evts,m_key);
+ if (events_statuscode==StatusCode::FAILURE) return StatusCode::FAILURE;
+ for (auto evt: evts) {
// call JetFinder
- vector< fastjet::PseudoJet > inclusive_jets;
- for ( vector< baseAnalysis* >::const_iterator i( m_jetfinder.begin() ); i != m_jetfinder.end(); ++i ) {
+ std::vector< fastjet::PseudoJet > inclusive_jets;
+ for ( std::vector< baseAnalysis* >::const_iterator i( m_jetfinder.begin() ); i != m_jetfinder.end(); ++i ) {
inclusive_jets = (*i)->GetJet(evt);
}
ATH_MSG_DEBUG( "after jet finder" );
-
-// int ret_all = true;
- for ( vector< baseAnalysis* >::const_iterator analysis = m_analysis.begin(); analysis != m_analysis.end(); ++analysis ) {
+ for ( std::vector< baseAnalysis* >::const_iterator analysis = m_analysis.begin(); analysis != m_analysis.end(); ++analysis ) {
(*analysis)->SetJet( &inclusive_jets );
-// int ret =
(*analysis)->Process( evt ); // call analysis
-
-// if ( !ret ) {
-// ret_all = false;
-// }
}
ATH_MSG_DEBUG("After process function");
// delete jet object, missing Et etc.
- for ( vector<baseAnalysis*>::const_iterator i( m_jetfinder.begin() ); i != m_jetfinder.end(); ++i ) {
+ for ( std::vector<baseAnalysis*>::const_iterator i( m_jetfinder.begin() ); i != m_jetfinder.end(); ++i ) {
(*i)->ClearEvent(evt);
}
ATH_MSG_DEBUG("After jetfinder mem cleanup");
}
- }
-
+ for (size_t i=0;i<evts.size();i++) delete evts[i];
ATH_MSG_DEBUG("Successful execution");
return StatusCode::SUCCESS;
}
@@ -314,18 +256,18 @@ StatusCode HepMCAnalysis::finalize()
}
// average histograms
- for ( vector< baseAnalysis* >::const_iterator analysis = m_analysis.begin(); analysis != m_analysis.end(); ++analysis ) {
+ for ( std::vector< baseAnalysis* >::const_iterator analysis = m_analysis.begin(); analysis != m_analysis.end(); ++analysis ) {
(*analysis)->averagedHistograms();
}
// register histos with histogram service
TH1D *hist = 0;
- for ( vector< baseAnalysis* >::const_iterator ana = m_analysis.begin(); ana != m_analysis.end(); ++ana ) {
+ for ( std::vector< baseAnalysis* >::const_iterator ana = m_analysis.begin(); ana != m_analysis.end(); ++ana ) {
while ( ( hist =(*ana)->popHisto() ) ) {
TString histname("/hepmcanalysis/");
- string name = hist->GetName();
- string dir = name.substr(0, name.find("_"));
+ std::string name = hist->GetName();
+ std::string dir = name.substr(0, name.find("_"));
histname += dir;
histname += "/";
histname += hist->GetName();
@@ -342,12 +284,12 @@ StatusCode HepMCAnalysis::finalize()
ATH_MSG_INFO( "Histograms registered successfully" );
// clean up analysis modules
- for ( vector< baseAnalysis* >::const_iterator j = m_jetfinder.begin(); j != m_jetfinder.end(); ++j ) {
+ for ( std::vector< baseAnalysis* >::const_iterator j = m_jetfinder.begin(); j != m_jetfinder.end(); ++j ) {
delete (*j);
}
m_jetfinder.clear();
- for ( vector< baseAnalysis* >::const_iterator ana = m_analysis.begin(); ana != m_analysis.end(); ++ana ) {
+ for ( std::vector< baseAnalysis* >::const_iterator ana = m_analysis.begin(); ana != m_analysis.end(); ++ana ) {
delete (*ana);
}
m_analysis.clear();
diff --git a/Generators/HepMCAnalysis_i/src/LeptonJetAnalysis.cxx b/Generators/HepMCAnalysis_i/src/LeptonJetAnalysis.cxx
index afe57c9276bfacecbc241f773ed1b35052c54c56..dca2138cce8ab55d433a017dad53b005bffebb5e 100644
--- a/Generators/HepMCAnalysis_i/src/LeptonJetAnalysis.cxx
+++ b/Generators/HepMCAnalysis_i/src/LeptonJetAnalysis.cxx
@@ -6,16 +6,6 @@
#include <iostream>
-#include "AtlasHepMC/GenEvent.h"
-#include "AtlasHepMC/IO_GenEvent.h"
-#include "AtlasHepMC/GenParticle.h"
-#include "AtlasHepMC/GenVertex.h"
-#include "AtlasHepMC/IO_AsciiParticles.h"
-#include "AtlasHepMC/SimpleVector.h"
-#include "AtlasHepMC/WeightContainer.h"
-#include "CLHEP/Vector/LorentzVector.h"
-
-
#include "TH1.h"
#include "TMath.h"
#include "TLorentzVector.h"
@@ -25,12 +15,22 @@
#include "fastjet/JetDefinition.hh"
#include "fastjet/SISConePlugin.hh"
-#include "TruthUtils/HepMCHelpers.h" //MCUtils.h
+/** Here we state explicitely what is used */
+#include "MCUtils/PIDUtils.h"
+#include "MCUtils/HepMCEventUtils.h"
+#include "MCUtils/HepMCParticleUtils.h"
+#include "MCUtils/HepMCVectors.h"
+#include "MCUtils/HepMCParticleClassifiers.h"
#include "HEPUtils/FastJet.h"
#include "../HepMCAnalysis_i/LeptonJetAnalysis.h"
-using namespace std;
+namespace MC {
+using namespace MCUtils::PID;
+using MCUtils::ptGtr;
+using MCUtils::get_weight;
+using MCUtils::isLastReplica;
+}
// ----------------------------------------------------------------------
LeptonJetAnalysis::LeptonJetAnalysis()
@@ -298,26 +298,26 @@ int LeptonJetAnalysis::Init(double maxEta, double minPt)
for(UInt_t i = 0; i < nleptons; ++i){
m_histName.Form("lepton%d_pT", i+1);
- m_leptonPt[i] = initHist(string( m_histName ), string( m_histName ), string( "p_{T} [GeV]") , 100, 0., 300.);
+ m_leptonPt[i] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "p_{T} [GeV]") , 100, 0., 300.);
m_leptonPt[i] -> Sumw2();
m_histName.Form("lepton%d_eta", i+1);
- m_leptonEta[i] = initHist(string( m_histName ), string( m_histName ), string( "#eta") , 50, -4, 4);
+ m_leptonEta[i] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "#eta") , 50, -4, 4);
m_leptonEta[i] -> Sumw2();
}
for(UInt_t i = 0; i < nleptons-1; ++i){
for (UInt_t j=i+1; j < nleptons;++j){
m_histName.Form("lepton%d%d_deltaR", i+1,j+1);
- m_leptondR[i][j] = initHist(string( m_histName ), string( m_histName ), string( "#Delta R" ), 64, 0, 3.2);
+ m_leptondR[i][j] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "#Delta R" ), 64, 0, 3.2);
m_leptondR[i][j] -> Sumw2();
m_histName.Form("lepton%d%d_deltaPhi", i+1,j+1);
- m_leptondPhi[i][j] = initHist(string( m_histName ), string( m_histName ), string( "#Delta #phi" ), 40, -3.14, 3.14);
+ m_leptondPhi[i][j] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "#Delta #phi" ), 40, -3.14, 3.14);
m_leptondPhi[i][j] -> Sumw2();
m_histName.Form("lepton%d%d_Mass", i+1,j+1);
- m_leptonMass[i][j] = initHist(string( m_histName ), string( m_histName ), string( "mass [GeV]" ), 40, 40.0, 120.0);
+ m_leptonMass[i][j] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "mass [GeV]" ), 40, 40.0, 120.0);
m_leptonMass[i][j] -> Sumw2();
m_histName.Form("lepton%d%d_lowMass", i+1,j+1);
- m_leptonLowMass[i][j] = initHist(string( m_histName ), string( m_histName ), string( "mass [GeV]" ), 20, 0.0, 40.0);
+ m_leptonLowMass[i][j] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "mass [GeV]" ), 20, 0.0, 40.0);
m_leptonLowMass[i][j] -> Sumw2();
}
}
@@ -351,16 +351,16 @@ int LeptonJetAnalysis::Init(double maxEta, double minPt)
for(UInt_t i = 0; i < njets; ++i){
m_histName.Form("jet%d_pT", i+1);
- m_jetPt[i] = initHist(string( m_histName ), string( m_histName ), string( "p_{T} [GeV]") , 100, 0., 300.);
+ m_jetPt[i] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "p_{T} [GeV]") , 100, 0., 300.);
m_jetPt[i] -> Sumw2();
m_histName.Form("jet%d_pT_highrange", i+1);
- m_jetPtHighRange[i] = initHist(string( m_histName ), string( m_histName ), string( "p_{T} [GeV]") , 300, 0., 1000.);
+ m_jetPtHighRange[i] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "p_{T} [GeV]") , 300, 0., 1000.);
m_jetPtHighRange[i] -> Sumw2();
m_histName.Form("jet%d_eta", i+1);
- m_jetEta[i] = initHist(string( m_histName ), string( m_histName ), string( "#eta") , 50, -4, 4);
+ m_jetEta[i] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "#eta") , 50, -4, 4);
m_jetEta[i] -> Sumw2();
m_histName.Form("jet%d_MassPt", i+1);
- m_jetMassPt[i] = initHist(string( m_histName ), string( m_histName ), string( "mass^{2}/p_{T}^{2}") , 100, 0.0, 0.1);
+ m_jetMassPt[i] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "mass^{2}/p_{T}^{2}") , 100, 0.0, 0.1);
m_jetMassPt[i] -> Sumw2();
}
@@ -368,13 +368,13 @@ int LeptonJetAnalysis::Init(double maxEta, double minPt)
for(UInt_t i = 0; i < njets-1; ++i){
for (UInt_t j=i+1; j < njets;++j){
m_histName.Form("jet%d%d_deltaR", i+1,j+1);
- m_jetdR[i][j] = initHist(string( m_histName ), string( m_histName ), string( "#Delta R" ), 64, 0, 3.2);
+ m_jetdR[i][j] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "#Delta R" ), 64, 0, 3.2);
m_jetdR[i][j] -> Sumw2();
m_histName.Form("jet%d%d_deltaPhi", i+1,j+1);
- m_jetdPhi[i][j] = initHist(string( m_histName ), string( m_histName ), string( "#Delta #phi" ), 40, -3.14, 3.14);
+ m_jetdPhi[i][j] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "#Delta #phi" ), 40, -3.14, 3.14);
m_jetdPhi[i][j] -> Sumw2();
m_histName.Form("jet%d%d_Mass", i+1,j+1);
- m_jetMass[i][j] = initHist(string( m_histName ), string( m_histName ), string( "mass [GeV]" ), 50, 0.0, 50.0);
+ m_jetMass[i][j] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "mass [GeV]" ), 50, 0.0, 50.0);
m_jetMass[i][j] -> Sumw2();
}
}
@@ -388,29 +388,29 @@ int LeptonJetAnalysis::Init(double maxEta, double minPt)
for(UInt_t i = 0; i < njets; ++i){
m_histName.Form("jet%d_pT_nocuts", i+1);
- m_jetPt_nocuts[i] = initHist(string( m_histName ), string( m_histName ), string( "p_{T} [GeV]") , 100, 0., 300.);
+ m_jetPt_nocuts[i] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "p_{T} [GeV]") , 100, 0., 300.);
m_jetPt_nocuts[i] -> Sumw2();
m_histName.Form("jet%d_pT_highrange_nocuts", i+1);
- m_jetPtHighRange_nocuts[i] = initHist(string( m_histName ), string( m_histName ), string( "p_{T} [GeV]") , 300, 0., 1000.);
+ m_jetPtHighRange_nocuts[i] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "p_{T} [GeV]") , 300, 0., 1000.);
m_jetPtHighRange_nocuts[i] -> Sumw2();
m_histName.Form("jet%d_eta_nocuts", i+1);
- m_jetEta_nocuts[i] = initHist(string( m_histName ), string( m_histName ), string( "#eta") , 50, -4, 4);
+ m_jetEta_nocuts[i] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "#eta") , 50, -4, 4);
m_jetEta_nocuts[i] -> Sumw2();
m_histName.Form("jet%d_MassPt_nocuts", i+1);
- m_jetMassPt_nocuts[i] = initHist(string( m_histName ), string( m_histName ), string( "mass^{2}/p_{T}^{2}") , 100, 0.0, 0.1);
+ m_jetMassPt_nocuts[i] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "mass^{2}/p_{T}^{2}") , 100, 0.0, 0.1);
m_jetMassPt_nocuts[i] -> Sumw2();
}
for(UInt_t i = 0; i < njets-1; ++i){
for (UInt_t j=i+1; j < njets;++j){
m_histName.Form("jet%d%d_deltaR_nocuts", i+1,j+1);
- m_jetdR_nocuts[i][j] = initHist(string( m_histName ), string( m_histName ), string( "#Delta R" ), 64, 0, 3.2);
+ m_jetdR_nocuts[i][j] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "#Delta R" ), 64, 0, 3.2);
m_jetdR_nocuts[i][j] -> Sumw2();
m_histName.Form("jet%d%d_deltaPhi_nocuts", i+1,j+1);
- m_jetdPhi_nocuts[i][j] = initHist(string( m_histName ), string( m_histName ), string( "#Delta #phi" ), 40, -3.14, 3.14);
+ m_jetdPhi_nocuts[i][j] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "#Delta #phi" ), 40, -3.14, 3.14);
m_jetdPhi_nocuts[i][j] -> Sumw2();
m_histName.Form("jet%d%d_Mass_nocuts", i+1,j+1);
- m_jetMass_nocuts[i][j] = initHist(string( m_histName ), string( m_histName ), string( "mass [GeV]" ), 50, 0.0, 50.0);
+ m_jetMass_nocuts[i][j] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "mass [GeV]" ), 50, 0.0, 50.0);
m_jetMass_nocuts[i][j] -> Sumw2();
}
}
@@ -424,29 +424,29 @@ int LeptonJetAnalysis::Init(double maxEta, double minPt)
for(UInt_t i = 0; i < njets; ++i){
m_histName.Form("jet%d_pT_forward", i+1);
- m_jetPt_forward[i] = initHist(string( m_histName ), string( m_histName ), string( "p_{T} [GeV]") , 100, 0., 300.);
+ m_jetPt_forward[i] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "p_{T} [GeV]") , 100, 0., 300.);
m_jetPt_forward[i] -> Sumw2();
m_histName.Form("jet%d_pT_highrange_forward", i+1);
- m_jetPtHighRange_forward[i] = initHist(string( m_histName ), string( m_histName ), string( "p_{T} [GeV]") , 300, 0., 1000.);
+ m_jetPtHighRange_forward[i] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "p_{T} [GeV]") , 300, 0., 1000.);
m_jetPtHighRange_forward[i] -> Sumw2();
m_histName.Form("jet%d_eta_forward", i+1);
- m_jetEta_forward[i] = initHist(string( m_histName ), string( m_histName ), string( "#eta") , 50, -4, 4);
+ m_jetEta_forward[i] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "#eta") , 50, -4, 4);
m_jetEta_forward[i] -> Sumw2();
m_histName.Form("jet%d_MassPt_forward", i+1);
- m_jetMassPt_forward[i] = initHist(string( m_histName ), string( m_histName ), string( "mass^{2}/p_{T}^{2}") , 100, 0.0, 0.1);
+ m_jetMassPt_forward[i] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "mass^{2}/p_{T}^{2}") , 100, 0.0, 0.1);
m_jetMassPt_forward[i] -> Sumw2();
}
for(UInt_t i = 0; i < njets-1; ++i){
for (UInt_t j=i+1; j < njets;++j){
m_histName.Form("jet%d%d_deltaR_forward", i+1,j+1);
- m_jetdR_forward[i][j] = initHist(string( m_histName ), string( m_histName ), string( "#Delta R" ), 64, 0, 3.2);
+ m_jetdR_forward[i][j] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "#Delta R" ), 64, 0, 3.2);
m_jetdR_forward[i][j] -> Sumw2();
m_histName.Form("jet%d%d_deltaPhi_forward", i+1,j+1);
- m_jetdPhi_forward[i][j] = initHist(string( m_histName ), string( m_histName ), string( "#Delta #phi" ), 40, -3.14, 3.14);
+ m_jetdPhi_forward[i][j] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "#Delta #phi" ), 40, -3.14, 3.14);
m_jetdPhi_forward[i][j] -> Sumw2();
m_histName.Form("jet%d%d_Mass_forward", i+1,j+1);
- m_jetMass_forward[i][j] = initHist(string( m_histName ), string( m_histName ), string( "mass [GeV]" ), 50, 0.0, 50.0);
+ m_jetMass_forward[i][j] = initHist(std::string( m_histName ), std::string( m_histName ), std::string( "mass [GeV]" ), 50, 0.0, 50.0);
m_jetMass_forward[i][j] -> Sumw2();
}
}
@@ -520,7 +520,7 @@ int LeptonJetAnalysis::Process(HepMC::GenEvent *event)
{
double weight;
- HepMC::GenParticle *lepton = 0;
+ HepMC::GenParticle* lepton = nullptr;
std::vector<fastjet::PseudoJet> selected_jets;
std::vector<fastjet::PseudoJet> forward_jets;
@@ -538,7 +538,7 @@ int LeptonJetAnalysis::Process(HepMC::GenEvent *event)
Int_t NumLeptons = 0, NumElectrons = 0, NumMuons = 0, NumTaus = 0, NumZ = 0, NumW = 0, NumGamma = 0;
if (event -> event_number() % 1000 == 0) {
- cout << "Processing event " << event -> event_number() << endl;
+ std::cout << "Processing event " << event -> event_number() << std::endl;
}
weight = MC::get_weight(event, 0 );
m_event_weight -> Fill(weight);
@@ -547,7 +547,7 @@ int LeptonJetAnalysis::Process(HepMC::GenEvent *event)
for (HepMC::GenEvent::particle_const_iterator p = event -> particles_begin(); p != event -> particles_end(); ++p){
- if ( abs((*p) -> pdg_id()) == 23 && (MC::isLastReplica(*p))){
+ if ( std::abs((*p) -> pdg_id()) == 23 && (MC::isLastReplica(*p))){
NumZ++;
m_ZNum->Fill(weight);
m_Z1pt->Fill((*p) -> momentum().perp(), weight);
@@ -555,7 +555,7 @@ int LeptonJetAnalysis::Process(HepMC::GenEvent *event)
}
- if ( abs((*p) -> pdg_id()) == 24 && (MC::isLastReplica(*p))){
+ if ( std::abs((*p) -> pdg_id()) == 24 && (MC::isLastReplica(*p))){
NumW++;
m_WNum->Fill(weight);
m_W1pt->Fill((*p) -> momentum().perp(), weight);
@@ -563,7 +563,7 @@ int LeptonJetAnalysis::Process(HepMC::GenEvent *event)
}
- if ( abs((*p) -> pdg_id()) == 22 && (MC::isLastReplica(*p))){
+ if ( std::abs((*p) -> pdg_id()) == 22 && (MC::isLastReplica(*p))){
NumGamma++;
m_gammaNum->Fill(weight);
m_gamma1pt->Fill((*p) -> momentum().perp(), weight);
@@ -573,32 +573,32 @@ int LeptonJetAnalysis::Process(HepMC::GenEvent *event)
- if ((MC::PID::isLepton((*p) -> pdg_id())) && (MC::isLastReplica(*p))) // last replica of charged lepton selected
+ if ((MC::isLepton((*p) -> pdg_id())) && (MC::isLastReplica(*p))) // last replica of charged lepton selected
{
lepton = (*p);
if (!MC::ptGtr(lepton, ptmin_l)) continue;
- if (!MC::isStable(lepton)) continue;
- if (!MC::isCharged(lepton)) continue;
+ if (!(lepton->status()==1)) continue;
+ if (!MC::isCharged(lepton->pdg_id())) continue;
int charge = ((*p) -> pdg_id() > 0) ? 1 : (((*p) -> pdg_id() < 0) ? -1 : 0);
m_leptonInclPt -> Fill(lepton -> momentum().perp(), weight);
//fill plots for electrons
- if (abs(lepton -> pdg_id()) == 11) {
+ if (std::abs(lepton -> pdg_id()) == 11) {
m_electronPt -> Fill(lepton -> momentum().perp(), weight);
m_electronCharge -> Fill(charge, weight);
NumElectrons++;
}
//fill plots for muons
- if (abs(lepton -> pdg_id()) == 13) {
+ if (std::abs(lepton -> pdg_id()) == 13) {
m_muonPt -> Fill(lepton -> momentum().perp(), weight);
m_muonCharge -> Fill(charge, weight);
NumMuons++;
}
//fill plots for taus
- //pdgid 17??? seems wrong if ((abs(lepton -> pdg_id()) == 15) || (abs(lepton -> pdg_id()) == 17)) {
- if (abs(lepton -> pdg_id()) == 15) {
+ //pdgid 17??? seems wrong if ((std::abs(lepton -> pdg_id()) == 15) || (std::abs(lepton -> pdg_id()) == 17)) {
+ if (std::abs(lepton -> pdg_id()) == 15) {
m_tauPt -> Fill(lepton -> momentum().perp(), weight);
m_tauCharge -> Fill(charge, weight);
NumTaus++;
@@ -610,7 +610,7 @@ int LeptonJetAnalysis::Process(HepMC::GenEvent *event)
const HEPUtils::P4 leptonMomentum(p4.px(), p4.py(), p4.pz(), p4.e());
unsorted_leptons.push_back(HEPUtils::mk_pseudojet(leptonMomentum));
- if( abs(lepton -> momentum().eta()) < etamax_l){
+ if( std::abs(lepton -> momentum().eta()) < etamax_l){
// unsorted_tight_leptons.push_back(MC::mk_pseudojet(lepton));
unsorted_tight_leptons.push_back(HEPUtils::mk_pseudojet(leptonMomentum));
}
@@ -649,13 +649,13 @@ int LeptonJetAnalysis::Process(HepMC::GenEvent *event)
//loop over jets and fill vectors
for (size_t i = 0; i < m_inclusive_jets.size(); ++i) {
- if ((abs(m_inclusive_jets[i].eta()) < etamax_jet) && (m_inclusive_jets[i].perp() > ptmin_jet)){
+ if ((std::abs(m_inclusive_jets[i].eta()) < etamax_jet) && (m_inclusive_jets[i].perp() > ptmin_jet)){
selected_jets.push_back(m_inclusive_jets[i]);
}
- if (abs(m_inclusive_jets[i].eta()) < etamax_jet){
+ if (std::abs(m_inclusive_jets[i].eta()) < etamax_jet){
no_cuts_jets.push_back(m_inclusive_jets[i]);
}
- if ( (abs(m_inclusive_jets[i].eta()) >= etamin_jetf) && (m_inclusive_jets[i].perp() > ptmin_jet)){
+ if ( (std::abs(m_inclusive_jets[i].eta()) >= etamin_jetf) && (m_inclusive_jets[i].perp() > ptmin_jet)){
forward_jets.push_back(m_inclusive_jets[i]);
}
}
@@ -671,7 +671,7 @@ int LeptonJetAnalysis::Process(HepMC::GenEvent *event)
}
for (UInt_t j = i + 1; j < selected_jets.size() ;++j){
if ((i < njets-1)&&(j < njets)){
- m_jetdR[i][j]->Fill(sqrt(selected_jets[i].squared_distance(selected_jets[j])), weight);
+ m_jetdR[i][j]->Fill(std::sqrt(selected_jets[i].squared_distance(selected_jets[j])), weight);
m_jetMass[i][j] -> Fill(selected_jets[i].m() + selected_jets[j].m(), weight);
m_jetdPhi[i][j]->Fill(selected_jets[i].delta_phi_to(selected_jets[j]), weight);
}
@@ -689,7 +689,7 @@ int LeptonJetAnalysis::Process(HepMC::GenEvent *event)
}
for (UInt_t j = i + 1; j < no_cuts_jets.size() ;++j){
if ((i < njets-1)&&(j < njets)){
- m_jetdR_nocuts[i][j]->Fill(sqrt(no_cuts_jets[i].squared_distance(no_cuts_jets[j])), weight);
+ m_jetdR_nocuts[i][j]->Fill(std::sqrt(no_cuts_jets[i].squared_distance(no_cuts_jets[j])), weight);
m_jetMass_nocuts[i][j] -> Fill(no_cuts_jets[i].m() + no_cuts_jets[j].m(), weight);
m_jetdPhi_nocuts[i][j]->Fill(no_cuts_jets[i].delta_phi_to(no_cuts_jets[j]), weight);
}
@@ -707,7 +707,7 @@ int LeptonJetAnalysis::Process(HepMC::GenEvent *event)
}
for (UInt_t j = i + 1; j < forward_jets.size() ;++j){
if ((i < njets-1)&&(j < njets)){
- m_jetdR_forward[i][j]->Fill(sqrt(forward_jets[i].squared_distance(forward_jets[j])), weight);
+ m_jetdR_forward[i][j]->Fill(std::sqrt(forward_jets[i].squared_distance(forward_jets[j])), weight);
m_jetMass_forward[i][j] -> Fill(forward_jets[i].m() + forward_jets[j].m(), weight);
m_jetdPhi_forward[i][j]->Fill(forward_jets[i].delta_phi_to(forward_jets[j]), weight);
}
@@ -729,7 +729,7 @@ int LeptonJetAnalysis::Process(HepMC::GenEvent *event)
m_jet1Pt_tight -> Fill(selected_jets[0].perp(), weight);
m_jet1PtHighRange_tight -> Fill(selected_jets[0].perp(), weight);
m_jet1Eta_tight -> Fill(selected_jets[0].eta(), weight);
- m_dR_jet1_to_lepton_tight -> Fill(sqrt(selected_jets[0].squared_distance(tight_leptons[0])), weight);
+ m_dR_jet1_to_lepton_tight -> Fill(std::sqrt(selected_jets[0].squared_distance(tight_leptons[0])), weight);
m_dPhi_jet1_to_lepton_tight -> Fill(selected_jets[0].delta_phi_to(tight_leptons[0]), weight);
dR1 = selected_jets[0].squared_distance (tight_leptons[0]);
@@ -740,15 +740,15 @@ int LeptonJetAnalysis::Process(HepMC::GenEvent *event)
m_jet2Pt_tight -> Fill(selected_jets[1].perp(), weight);
m_jet2PtHighRange_tight -> Fill(selected_jets[1].perp(), weight);
m_jet2Eta_tight -> Fill(selected_jets[1].eta(), weight);
- m_dR_jet2_to_lepton_tight -> Fill(sqrt(selected_jets[1].squared_distance(tight_leptons[0])), weight);
+ m_dR_jet2_to_lepton_tight -> Fill(std::sqrt(selected_jets[1].squared_distance(tight_leptons[0])), weight);
m_dPhi_jet2_to_lepton_tight -> Fill(selected_jets[1].delta_phi_to(tight_leptons[0]), weight);
dR2 = selected_jets[1].squared_distance(tight_leptons[0]);
dPhi2 = selected_jets[1].delta_phi_to(tight_leptons[0]);
}
- dRmin = ( abs(dR2) > abs(dR1)) ? dR1 :dR2;
- dPhimin = (abs(dPhi2) > abs(dPhi1)) ? dPhi1 :dPhi2;
+ dRmin = ( std::abs(dR2) > std::abs(dR1)) ? dR1 :dR2;
+ dPhimin = (std::abs(dPhi2) > std::abs(dPhi1)) ? dPhi1 :dPhi2;
m_dR_lepton_to_closest_jet_tight -> Fill(dRmin, weight);
m_dPhi_lepton_to_closest_jet_tight -> Fill(dPhimin, weight);
}
diff --git a/Generators/HepMCAnalysis_i/src/ParticleContentAnalysis.cxx b/Generators/HepMCAnalysis_i/src/ParticleContentAnalysis.cxx
index fe44ac3d0c57029b69a3cff7b3b94d99c381a832..c5a24da478e3381250b19b31f82bee6fea50d889 100644
--- a/Generators/HepMCAnalysis_i/src/ParticleContentAnalysis.cxx
+++ b/Generators/HepMCAnalysis_i/src/ParticleContentAnalysis.cxx
@@ -5,12 +5,6 @@
#include <iostream>
#include <sstream>
#include <stdio.h>
-#include "AtlasHepMC/GenEvent.h"
-#include "AtlasHepMC/IO_GenEvent.h"
-#include "AtlasHepMC/GenParticle.h"
-#include "AtlasHepMC/GenVertex.h"
-#include "AtlasHepMC/IO_AsciiParticles.h"
-#include "AtlasHepMC/SimpleVector.h"
#include "CLHEP/Vector/LorentzVector.h"
// ROOT headers
@@ -18,13 +12,23 @@
#include "TH2.h"
#include "TFile.h"
#include "TMath.h"
-#include "TLorentzVector.h"
-#include "TruthUtils/HepMCHelpers.h"
+/** Here we state explicitely what is used */
+#include "MCUtils/PIDUtils.h"
+#include "MCUtils/HepMCEventUtils.h"
+#include "MCUtils/HepMCParticleUtils.h"
+#include "MCUtils/HepMCVectors.h"
+#include "MCUtils/HepMCParticleClassifiers.h"
#include "../HepMCAnalysis_i/ParticleContentAnalysis.h"
-using namespace std;
+namespace MC {
+using namespace MCUtils::PID;
+using MCUtils::get_weight;
+using MCUtils::flightLength;
+using MCUtils::isLastReplica;
+}
+
// empty default constructor
ParticleContentAnalysis::ParticleContentAnalysis()
@@ -208,83 +212,83 @@ int ParticleContentAnalysis::Init(double tr_max_eta, double tr_min_pt)
m_outputRootDir = "PartCont";
//declaration of histograms
- m_evtnr = initHist( string( "EventNumber" ), string( "Event number" ), string("Eventnumber"), 100, 0., 1000. );
- m_evtweight = initHist( string( "EventWeight"), string( "event weight"), string("weight"), 1000, -5000000.0, 5000000.0);
- m_evtweight_zoom = initHist(string("EventWeight_Zoom"), string("event weight"), string("weight"), 41, -20.5, 20.5);
+ m_evtnr = initHist( std::string( "EventNumber" ), std::string( "Event number" ), std::string("Eventnumber"), 100, 0., 1000. );
+ m_evtweight = initHist( std::string( "EventWeight"), std::string( "event weight"), std::string("weight"), 1000, -5000000.0, 5000000.0);
+ m_evtweight_zoom = initHist(std::string("EventWeight_Zoom"), std::string("event weight"), std::string("weight"), 41, -20.5, 20.5);
//all particles in the event
- m_initial_particle_pid = initHist( string( "Initial_Particle_PID" ) , string( "PID of initial particles") , string(" pid" ), 30001., -15000.5, 15000.5 );
- m_intermediate_particle_pid = initHist( string( "Intermediate_Particle_PID" ) , string( "PID of intermediate particles" ), string(" pid" ), 30001., -15000.5, 15000.5 );
- m_final_particle_pid = initHist( string( "Final_Particle_PID") , string( "PID of final particles" ), string(" pid" ), 30001., -15000.5, 15000.5 );
+ m_initial_particle_pid = initHist( std::string( "Initial_Particle_PID" ) , std::string( "PID of initial particles") , std::string(" pid" ), 30001., -15000.5, 15000.5 );
+ m_intermediate_particle_pid = initHist( std::string( "Intermediate_Particle_PID" ) , std::string( "PID of intermediate particles" ), std::string(" pid" ), 30001., -15000.5, 15000.5 );
+ m_final_particle_pid = initHist( std::string( "Final_Particle_PID") , std::string( "PID of final particles" ), std::string(" pid" ), 30001., -15000.5, 15000.5 );
- m_initial_particle_pid_zoom = initHist( string( "Initial_Particle_PID_zoom" ) , string( "PID of initial particles") , string(" pid" ), 61., -30.5, 30.5 );
- m_intermediate_particle_pid_zoom = initHist( string( "Intermediate_Particle_PID_zoom" ) , string( "PID of intermediate particles" ), string(" pid" ), 61., -30.5, 30.5 );
- m_final_particle_pid_zoom = initHist( string( "Final_Particle_PID_zoom") , string( "PID of final particles" ), string(" pid" ), 61., -30.5, 30.5 );
+ m_initial_particle_pid_zoom = initHist( std::string( "Initial_Particle_PID_zoom" ) , std::string( "PID of initial particles") , std::string(" pid" ), 61., -30.5, 30.5 );
+ m_intermediate_particle_pid_zoom = initHist( std::string( "Intermediate_Particle_PID_zoom" ) , std::string( "PID of intermediate particles" ), std::string(" pid" ), 61., -30.5, 30.5 );
+ m_final_particle_pid_zoom = initHist( std::string( "Final_Particle_PID_zoom") , std::string( "PID of final particles" ), std::string(" pid" ), 61., -30.5, 30.5 );
- m_particle_status = initHist( string( "Particle_Status" ), string( "status of particles" ), string("status" ), 151., -0.5, 150.5 );
+ m_particle_status = initHist( std::string( "Particle_Status" ), std::string( "status of particles" ), std::string("status" ), 151., -0.5, 150.5 );
//bottom hadron plots
- m_BottomHadron_pt = initHist( string( "BottomHadron_Pt" ), string( "p_{T} of bottom hadron" ), string( "p_{T} [GeV]" ), 200., 0., 400. );
- m_BottomHadron_eta = initHist( string( "BottomHadron_Eta" ), string( "#eta of bottom hadron" ), string(" #eta" ), 60, -6., 6. );
- m_BottomHadron_phi = initHist( string( "BottomHadron_Phi" ), string( "#Phi of bottom hadron" ), string( "#Phi" ), 48, -3.15, 3.15 );
- m_BottomHadron_flightlength = initHist( string( "BottomHadron_FlightLength" ), string( "flight length of bottom hadron" ), string( "length" ), 100., -0.05, 150.05 );
- m_BottomHadron_status = initHist( string( "BottomHadron_Status" ), string( "status of bottom hadron" ), string("status" ), 151., -0.5, 150.5 );
- m_BottomHadron_decay_multiplicity = initHist( string( "BottomHadron_decay_multiplicity" ), string( "multiplicity of bottom hadron decay" ), string("N" ), 11., -0.5, 10.5 );
- m_BottomHadron_decay_charged_multiplicity = initHist( string( "BottomHadron_decay_charged_multiplicity" ), string( "charged multiplicity of bottom hadron decay" ), string("N" ), 11., -0.5, 10.5 );
-
- m_BottomHadron_tight_pt = initHist( string( "BottomHadron_tight_Pt" ), string( "p_{T} of bottom hadron" ), string( "p_{T} [GeV]" ), 200., 0., 400. );
- m_BottomHadron_tight_eta = initHist( string( "BottomHadron_tight_Eta" ), string( "#eta of bottom hadron" ), string(" #eta" ), 60, -6., 6. );
- m_BottomHadron_tight_phi = initHist( string( "BottomHadron_tight_Phi" ), string( "#Phi of bottom hadron" ), string( "#Phi" ), 48, -3.15, 3.15 );
- m_BottomHadron_tight_flightlength = initHist( string( "BottomHadron_tight_FlightLength" ), string( "flight length of bottom hadron" ), string( "length" ), 100., -0.05, 150.05 );
- m_BottomHadron_tight_status = initHist( string( "BottomHadron_tight_Status" ), string( "status of bottom hadron" ), string("status" ), 151., -0.5, 150.5 );
- m_BottomHadron_tight_decay_multiplicity = initHist( string( "BottomHadron_tight_decay_multiplicity" ), string( "multiplicity of bottom hadron decay" ), string("N" ), 11., -0.5, 10.5 );
- m_BottomHadron_tight_decay_charged_multiplicity = initHist( string( "BottomHadron_tight_decay_charged_multiplicity" ), string( "charged multiplicity of bottom hadron decay" ), string("N" ), 11., -0.5, 10.5 );
-
- m_BottomMeson_number = initHist( string( "NumBottomMeson"), string( "number of bottom meson"), string( "number"), 81, -0.5, 80.5);
- m_BottomBaryon_number = initHist( string( "NumBottomBaryon"), string( "number of bottom baryon"), string( "number"), 81, -0.5, 80.5);
+ m_BottomHadron_pt = initHist( std::string( "BottomHadron_Pt" ), std::string( "p_{T} of bottom hadron" ), std::string( "p_{T} [GeV]" ), 200., 0., 400. );
+ m_BottomHadron_eta = initHist( std::string( "BottomHadron_Eta" ), std::string( "#eta of bottom hadron" ), std::string(" #eta" ), 60, -6., 6. );
+ m_BottomHadron_phi = initHist( std::string( "BottomHadron_Phi" ), std::string( "#Phi of bottom hadron" ), std::string( "#Phi" ), 48, -3.15, 3.15 );
+ m_BottomHadron_flightlength = initHist( std::string( "BottomHadron_FlightLength" ), std::string( "flight length of bottom hadron" ), std::string( "length" ), 100., -0.05, 150.05 );
+ m_BottomHadron_status = initHist( std::string( "BottomHadron_Status" ), std::string( "status of bottom hadron" ), std::string("status" ), 151., -0.5, 150.5 );
+ m_BottomHadron_decay_multiplicity = initHist( std::string( "BottomHadron_decay_multiplicity" ), std::string( "multiplicity of bottom hadron decay" ), std::string("N" ), 11., -0.5, 10.5 );
+ m_BottomHadron_decay_charged_multiplicity = initHist( std::string( "BottomHadron_decay_charged_multiplicity" ), std::string( "charged multiplicity of bottom hadron decay" ), std::string("N" ), 11., -0.5, 10.5 );
+
+ m_BottomHadron_tight_pt = initHist( std::string( "BottomHadron_tight_Pt" ), std::string( "p_{T} of bottom hadron" ), std::string( "p_{T} [GeV]" ), 200., 0., 400. );
+ m_BottomHadron_tight_eta = initHist( std::string( "BottomHadron_tight_Eta" ), std::string( "#eta of bottom hadron" ), std::string(" #eta" ), 60, -6., 6. );
+ m_BottomHadron_tight_phi = initHist( std::string( "BottomHadron_tight_Phi" ), std::string( "#Phi of bottom hadron" ), std::string( "#Phi" ), 48, -3.15, 3.15 );
+ m_BottomHadron_tight_flightlength = initHist( std::string( "BottomHadron_tight_FlightLength" ), std::string( "flight length of bottom hadron" ), std::string( "length" ), 100., -0.05, 150.05 );
+ m_BottomHadron_tight_status = initHist( std::string( "BottomHadron_tight_Status" ), std::string( "status of bottom hadron" ), std::string("status" ), 151., -0.5, 150.5 );
+ m_BottomHadron_tight_decay_multiplicity = initHist( std::string( "BottomHadron_tight_decay_multiplicity" ), std::string( "multiplicity of bottom hadron decay" ), std::string("N" ), 11., -0.5, 10.5 );
+ m_BottomHadron_tight_decay_charged_multiplicity = initHist( std::string( "BottomHadron_tight_decay_charged_multiplicity" ), std::string( "charged multiplicity of bottom hadron decay" ), std::string("N" ), 11., -0.5, 10.5 );
+
+ m_BottomMeson_number = initHist( std::string( "NumBottomMeson"), std::string( "number of bottom meson"), std::string( "number"), 81, -0.5, 80.5);
+ m_BottomBaryon_number = initHist( std::string( "NumBottomBaryon"), std::string( "number of bottom baryon"), std::string( "number"), 81, -0.5, 80.5);
//charmed hadron plots
- m_CharmHadron_pt = initHist( string( "CharmHadron_Pt" ), string( "p_{T} of charmed hadron" ), string( "p_{T} [GeV]" ), 200., 0., 400.);
- m_CharmHadron_eta = initHist( string( "CharmHadron_Eta" ), string( "#eta of charmed hadron" ), string( "#eta" ), 60, -6., 6.);
- m_CharmHadron_phi = initHist( string( "CharmHadron_Phi" ), string( "#Phi of charmed hadron" ), string( "#Phi" ), 48, -3.15, 3.15);
- m_CharmHadron_flightlength = initHist( string( "CharmHadron_FlightLength" ), string( "flight length of charmed hadron" ), string( "length" ), 100., -0.05, 150.05 );
- m_CharmHadron_status = initHist( string( "CharmHadron_Status" ), string( "status of charmed hadron" ), string("status" ), 151., -0.5, 150.5 );
- m_CharmHadron_decay_multiplicity = initHist( string( "CharmHadron_decay_multiplicity" ), string( "multiplicity of charmed hadron decay" ), string("N" ), 11., -0.5, 10.5 );
- m_CharmHadron_decay_charged_multiplicity = initHist( string( "CharmHadron_decay_charged_multiplicity" ), string( "charged multiplicity of charmed hadron decay" ), string("N" ), 11., -0.5, 10.5 );
- m_CharmMeson_number = initHist( string( "NumCharmMeson"), string( "number of charm meson"), string( "number"), 81, -0.5, 80.5);
- m_CharmBaryon_number = initHist( string( "NumCharmBaryon"), string( "number of charm baryon"), string( "number"), 81, -0.5, 80.5);
+ m_CharmHadron_pt = initHist( std::string( "CharmHadron_Pt" ), std::string( "p_{T} of charmed hadron" ), std::string( "p_{T} [GeV]" ), 200., 0., 400.);
+ m_CharmHadron_eta = initHist( std::string( "CharmHadron_Eta" ), std::string( "#eta of charmed hadron" ), std::string( "#eta" ), 60, -6., 6.);
+ m_CharmHadron_phi = initHist( std::string( "CharmHadron_Phi" ), std::string( "#Phi of charmed hadron" ), std::string( "#Phi" ), 48, -3.15, 3.15);
+ m_CharmHadron_flightlength = initHist( std::string( "CharmHadron_FlightLength" ), std::string( "flight length of charmed hadron" ), std::string( "length" ), 100., -0.05, 150.05 );
+ m_CharmHadron_status = initHist( std::string( "CharmHadron_Status" ), std::string( "status of charmed hadron" ), std::string("status" ), 151., -0.5, 150.5 );
+ m_CharmHadron_decay_multiplicity = initHist( std::string( "CharmHadron_decay_multiplicity" ), std::string( "multiplicity of charmed hadron decay" ), std::string("N" ), 11., -0.5, 10.5 );
+ m_CharmHadron_decay_charged_multiplicity = initHist( std::string( "CharmHadron_decay_charged_multiplicity" ), std::string( "charged multiplicity of charmed hadron decay" ), std::string("N" ), 11., -0.5, 10.5 );
+ m_CharmMeson_number = initHist( std::string( "NumCharmMeson"), std::string( "number of charm meson"), std::string( "number"), 81, -0.5, 80.5);
+ m_CharmBaryon_number = initHist( std::string( "NumCharmBaryon"), std::string( "number of charm baryon"), std::string( "number"), 81, -0.5, 80.5);
//strange hadron plots
- m_StrangeHadron_pt = initHist( string( "StrangeHadron_Pt" ), string( "p_{T} of strange hadron" ), string( "p_{T} [GeV]" ), 200., 0., 400. );
- m_StrangeHadron_unstable_pt = initHist( string( "StrangeHadron_Unstable_Pt" ), string("p_{T} of the last strange hadron (unstable)" ), string( "p_{T} [GeV]" ), 200.,0.,400. );
- m_StrangeHadron_eta = initHist( string( "StrangeHadron_Eta" ), string( "#eta of strange hadron" ), string( "#eta" ), 60, -6., 6.);
- m_StrangeHadron_phi = initHist( string( "StrangeHadron_Phi" ), string( "#Phi of strange hadron" ), string( "#Phi" ), 48, -3.15, 3.15);
- m_StrangeHadron_flightlength = initHist( string( "StrangeHadron_FlightLength" ), string( "flight length of strange hadron" ), string( "length" ), 100., -0.05, 150.05 );
- m_StrangeHadron_status = initHist( string( "StrangeHadron_Status" ), string( "status of strange hadron" ), string("status" ), 151., -0.5, 150.5 );
- m_StrangeHadron_decay_multiplicity = initHist( string( "StrangeHadron_decay_multiplicity" ), string( "multiplicity of strange hadron decay" ), string("N" ), 11., -0.5, 10.5 );
- m_StrangeHadron_decay_charged_multiplicity = initHist( string( "StrangeHadron_decay_charged_multiplicity" ), string( "charged multiplicity of strange hadron decay" ), string("N" ), 11., -0.5, 10.5 );
-
- m_StrangeMeson_number = initHist( string( "NumStrangeMeson"), string( "number of strange meson"), string( "number"), 81, -0.5, 80.5);
- m_StrangeBaryon_number = initHist( string( "NumStrangeBaryon"), string( "number of strange baryon"), string( "number"), 81, -0.5, 80.5);
+ m_StrangeHadron_pt = initHist( std::string( "StrangeHadron_Pt" ), std::string( "p_{T} of strange hadron" ), std::string( "p_{T} [GeV]" ), 200., 0., 400. );
+ m_StrangeHadron_unstable_pt = initHist( std::string( "StrangeHadron_Unstable_Pt" ), std::string("p_{T} of the last strange hadron (unstable)" ), std::string( "p_{T} [GeV]" ), 200.,0.,400. );
+ m_StrangeHadron_eta = initHist( std::string( "StrangeHadron_Eta" ), std::string( "#eta of strange hadron" ), std::string( "#eta" ), 60, -6., 6.);
+ m_StrangeHadron_phi = initHist( std::string( "StrangeHadron_Phi" ), std::string( "#Phi of strange hadron" ), std::string( "#Phi" ), 48, -3.15, 3.15);
+ m_StrangeHadron_flightlength = initHist( std::string( "StrangeHadron_FlightLength" ), std::string( "flight length of strange hadron" ), std::string( "length" ), 100., -0.05, 150.05 );
+ m_StrangeHadron_status = initHist( std::string( "StrangeHadron_Status" ), std::string( "status of strange hadron" ), std::string("status" ), 151., -0.5, 150.5 );
+ m_StrangeHadron_decay_multiplicity = initHist( std::string( "StrangeHadron_decay_multiplicity" ), std::string( "multiplicity of strange hadron decay" ), std::string("N" ), 11., -0.5, 10.5 );
+ m_StrangeHadron_decay_charged_multiplicity = initHist( std::string( "StrangeHadron_decay_charged_multiplicity" ), std::string( "charged multiplicity of strange hadron decay" ), std::string("N" ), 11., -0.5, 10.5 );
+
+ m_StrangeMeson_number = initHist( std::string( "NumStrangeMeson"), std::string( "number of strange meson"), std::string( "number"), 81, -0.5, 80.5);
+ m_StrangeBaryon_number = initHist( std::string( "NumStrangeBaryon"), std::string( "number of strange baryon"), std::string( "number"), 81, -0.5, 80.5);
//tau plots
- m_Tau_decay_multiplicity = initHist( string( "Tau_decay_multiplicity") , string("multiplicity of tau decay") , string("N" ), 11., -0.5, 10.5 );
- m_Tau_decay_charged_multiplicity = initHist( string( "Tau_charged_decay_multiplicity") , string("charged multiplicity of tau decay") , string("N" ), 11., -0.5, 10.5 );
- m_Tau_cos_theta = initHist( string( "Tau_cos_theta") , string("cos(#theta)") , string("cos(#theta)") , 80., -1.0, 1.0 );
- m_Tau_charged_energy_fraction = initHist( string( "Tau_charged_energy_fraction") , string("charged energy fraction") , string("E_{charged}/E") , 40., -1.0, 1.0 );
- m_Tau_flightlength = initHist( string( "Tau_FlightLength") , string( "flight length of tau") , string( "length" ), 100., -0.05, 150.05 );
-
- m_uquark_mass = initHist( string( "uquark_mass" ), string( "mass of u-quark" ), string( "mass [GeV]" ), 101., -0.25, 50.25 );
- m_dquark_mass = initHist( string( "dquark_mass" ), string( "mass of d-quark" ), string( "mass [GeV]" ), 101., -0.25, 50.25 );
- m_cquark_mass = initHist( string( "cquark_mass" ), string( "mass of c-quark" ), string( "mass [GeV]" ), 101., -0.25, 50.25 );
- m_squark_mass = initHist( string( "squark_mass" ), string( "mass of s-quark" ), string( "mass [GeV]" ), 101., -0.25, 50.25 );
- m_tquark_mass = initHist( string( "tquark_mass" ), string( "mass of t-quark" ), string( "mass [GeV]" ), 201., 99.5, 200.5 );
- m_bquark_mass = initHist( string( "bquark_mass" ), string( "mass of b-quark" ), string( "mass [GeV]" ), 101., -0.25, 50.25 );
- m_gluon_mass = initHist( string( "gluon_mass" ), string( "mass of gluon" ), string( "mass [GeV]" ), 101., -0.25, 50.25 );
- m_photon_mass = initHist( string( "photon_mass" ), string( "mass of photon" ), string( "mass [GeV]" ), 101., -0.25, 50.25 );
- m_Zboson_mass = initHist( string( "Zboson_mass" ), string( "mass of Z boson" ), string( "mass [GeV]" ), 100., 50., 150 );
- m_Wboson_mass = initHist( string( "Wboson_mass" ), string( "mass of W boson" ), string( "mass [GeV]" ), 100., 50., 150 );
+ m_Tau_decay_multiplicity = initHist( std::string( "Tau_decay_multiplicity") , std::string("multiplicity of tau decay") , std::string("N" ), 11., -0.5, 10.5 );
+ m_Tau_decay_charged_multiplicity = initHist( std::string( "Tau_charged_decay_multiplicity") , std::string("charged multiplicity of tau decay") , std::string("N" ), 11., -0.5, 10.5 );
+ m_Tau_cos_theta = initHist( std::string( "Tau_cos_theta") , std::string("cos(#theta)") , std::string("cos(#theta)") , 80., -1.0, 1.0 );
+ m_Tau_charged_energy_fraction = initHist( std::string( "Tau_charged_energy_fraction") , std::string("charged energy fraction") , std::string("E_{charged}/E") , 40., -1.0, 1.0 );
+ m_Tau_flightlength = initHist( std::string( "Tau_FlightLength") , std::string( "flight length of tau") , std::string( "length" ), 100., -0.05, 150.05 );
+
+ m_uquark_mass = initHist( std::string( "uquark_mass" ), std::string( "mass of u-quark" ), std::string( "mass [GeV]" ), 101., -0.25, 50.25 );
+ m_dquark_mass = initHist( std::string( "dquark_mass" ), std::string( "mass of d-quark" ), std::string( "mass [GeV]" ), 101., -0.25, 50.25 );
+ m_cquark_mass = initHist( std::string( "cquark_mass" ), std::string( "mass of c-quark" ), std::string( "mass [GeV]" ), 101., -0.25, 50.25 );
+ m_squark_mass = initHist( std::string( "squark_mass" ), std::string( "mass of s-quark" ), std::string( "mass [GeV]" ), 101., -0.25, 50.25 );
+ m_tquark_mass = initHist( std::string( "tquark_mass" ), std::string( "mass of t-quark" ), std::string( "mass [GeV]" ), 201., 99.5, 200.5 );
+ m_bquark_mass = initHist( std::string( "bquark_mass" ), std::string( "mass of b-quark" ), std::string( "mass [GeV]" ), 101., -0.25, 50.25 );
+ m_gluon_mass = initHist( std::string( "gluon_mass" ), std::string( "mass of gluon" ), std::string( "mass [GeV]" ), 101., -0.25, 50.25 );
+ m_photon_mass = initHist( std::string( "photon_mass" ), std::string( "mass of photon" ), std::string( "mass [GeV]" ), 101., -0.25, 50.25 );
+ m_Zboson_mass = initHist( std::string( "Zboson_mass" ), std::string( "mass of Z boson" ), std::string( "mass [GeV]" ), 100., 50., 150 );
+ m_Wboson_mass = initHist( std::string( "Wboson_mass" ), std::string( "mass of W boson" ), std::string( "mass [GeV]" ), 100., 50., 150 );
m_evtnr -> Sumw2();
@@ -409,7 +413,7 @@ int ParticleContentAnalysis::Process(HepMC::GenEvent* hepmcevt)
m_particle_status -> Fill((*p) -> status(), weight);
- if ((MC::isBottomHadron(*p)) && MC::isLastReplica(*p))
+ if ((MC::isBottomHadron((*p)->pdg_id())) && MC::isLastReplica(*p))
{
m_BottomHadron_pt -> Fill((*p) -> momentum().perp(), weight);
m_BottomHadron_eta -> Fill((*p) -> momentum().eta(), weight);
@@ -418,7 +422,7 @@ int ParticleContentAnalysis::Process(HepMC::GenEvent* hepmcevt)
m_BottomHadron_flightlength -> Fill(bflightlength, weight);
m_BottomHadron_status -> Fill((*p) -> status(), weight);
- if((*p) -> momentum().perp() > tight_bottom_pt && fabs((*p) -> momentum().eta())< tight_bottom_eta) isTightBottomHadron = true;
+ if((*p) -> momentum().perp() > tight_bottom_pt && std::abs((*p) -> momentum().eta())< tight_bottom_eta) isTightBottomHadron = true;
if(isTightBottomHadron){
m_BottomHadron_tight_pt -> Fill((*p) -> momentum().perp(), weight);
m_BottomHadron_tight_eta -> Fill((*p) -> momentum().eta(), weight);
@@ -437,7 +441,7 @@ int ParticleContentAnalysis::Process(HepMC::GenEvent* hepmcevt)
current_daughter != decay_vertex -> particles_end(HepMC::children); ++current_daughter)
{
NBD++;
- if (MC::isCharged(*current_daughter))NBCD++;
+ if (MC::isCharged((*current_daughter)->pdg_id()))NBCD++;
}
m_BottomHadron_decay_multiplicity -> Fill(NBD, weight);
m_BottomHadron_decay_charged_multiplicity -> Fill(NBCD, weight);
@@ -446,11 +450,11 @@ int ParticleContentAnalysis::Process(HepMC::GenEvent* hepmcevt)
}
}
- if ((MC::isBottomMeson(*p)) && MC::isLastReplica(*p)) NumBottomMeson++;
- if ((MC::isBottomBaryon(*p)) && MC::isLastReplica(*p)) NumBottomBaryon++;
+ if ((MC::isBottomMeson((*p)->pdg_id())) && MC::isLastReplica(*p)) NumBottomMeson++;
+ if ((MC::isBottomBaryon((*p)->pdg_id())) && MC::isLastReplica(*p)) NumBottomBaryon++;
- if ((MC::isCharmHadron(*p)) && MC::isLastReplica(*p))
+ if ((MC::isCharmHadron((*p)->pdg_id())) && MC::isLastReplica(*p))
{
m_CharmHadron_pt -> Fill((*p) -> momentum().perp(), weight);
m_CharmHadron_eta -> Fill((*p) -> momentum().eta(), weight);
@@ -467,20 +471,20 @@ int ParticleContentAnalysis::Process(HepMC::GenEvent* hepmcevt)
for (HepMC::GenVertex::particle_iterator current_daughter = decay_vertex -> particles_begin(HepMC::children); current_daughter != decay_vertex -> particles_end(HepMC::children); ++current_daughter)
{
NDD++;
- if (MC::isCharged(*current_daughter))NDCD++;
+ if (MC::isCharged((*current_daughter)->pdg_id()))NDCD++;
}
m_CharmHadron_decay_multiplicity -> Fill(NDD, weight);
m_CharmHadron_decay_charged_multiplicity -> Fill(NDCD, weight);
}
}
- if ((MC::isCharmMeson(*p)) && MC::isLastReplica(*p)) NumCharmMeson++;
- if ((MC::isCharmBaryon(*p)) && MC::isLastReplica(*p)) NumCharmBaryon++;
+ if ((MC::isCharmMeson((*p)->pdg_id())) && MC::isLastReplica(*p)) NumCharmMeson++;
+ if ((MC::isCharmBaryon((*p)->pdg_id())) && MC::isLastReplica(*p)) NumCharmBaryon++;
- if ((MC::isStrangeHadron(*p)) && MC::isLastReplica(*p))
+ if ((MC::isStrangeHadron((*p)->pdg_id())) && MC::isLastReplica(*p))
{
- if (!(MC::isStable(*p))){
+ if (!((*p)->status()==1)){
m_StrangeHadron_unstable_pt -> Fill((*p) -> momentum().perp(), weight);
}
m_StrangeHadron_pt -> Fill((*p) -> momentum().perp(), weight);
@@ -499,7 +503,7 @@ int ParticleContentAnalysis::Process(HepMC::GenEvent* hepmcevt)
for (HepMC::GenVertex::particle_iterator current_daughter = decay_vertex -> particles_begin(HepMC::children); current_daughter != decay_vertex -> particles_end(HepMC::children); ++current_daughter)
{
NSD++;
- if (MC::isCharged(*current_daughter))NSCD++;
+ if (MC::isCharged((*current_daughter)->pdg_id()))NSCD++;
}
m_StrangeHadron_decay_multiplicity -> Fill(NSD, weight);
m_StrangeHadron_decay_charged_multiplicity -> Fill(NSCD, weight);
@@ -507,11 +511,11 @@ int ParticleContentAnalysis::Process(HepMC::GenEvent* hepmcevt)
}
- if ((MC::isStrangeMeson(*p)) && MC::isLastReplica(*p)) NumStrangeMeson++;
- if ((MC::isStrangeBaryon(*p)) && MC::isLastReplica(*p)) NumStrangeBaryon++;
+ if ((MC::isStrangeMeson((*p)->pdg_id())) && MC::isLastReplica(*p)) NumStrangeMeson++;
+ if ((MC::isStrangeBaryon((*p)->pdg_id())) && MC::isLastReplica(*p)) NumStrangeBaryon++;
- if ( abs((*p) -> pdg_id()) == 15 && MC::isLastReplica(*p) )
+ if ( std::abs((*p) -> pdg_id()) == 15 && MC::isLastReplica(*p) )
{
Int_t NTauDecay = 0;
Int_t NTauChargedDecay = 0;
@@ -525,12 +529,12 @@ int ParticleContentAnalysis::Process(HepMC::GenEvent* hepmcevt)
if (decay_vertex){
for (HepMC::GenVertex::particle_iterator current_daughter = decay_vertex -> particles_begin(HepMC::children); current_daughter != decay_vertex -> particles_end(HepMC::children); ++current_daughter)
{
- if (abs((*current_daughter) -> pdg_id()) == 16) tau_neutrino = (*current_daughter);
- if (abs((*current_daughter) -> pdg_id()) == 211) charged_pion = (*current_daughter);
- if (abs((*current_daughter) -> pdg_id()) == 111) uncharged_pion = (*current_daughter);
+ if (std::abs((*current_daughter) -> pdg_id()) == 16) tau_neutrino = (*current_daughter);
+ if (std::abs((*current_daughter) -> pdg_id()) == 211) charged_pion = (*current_daughter);
+ if (std::abs((*current_daughter) -> pdg_id()) == 111) uncharged_pion = (*current_daughter);
NTauDecay++;
- if (MC::isCharged(*current_daughter)) NTauChargedDecay++;
+ if (MC::isCharged((*current_daughter)->pdg_id())) NTauChargedDecay++;
}
m_Tau_decay_multiplicity -> Fill(NTauDecay, weight);
@@ -548,16 +552,16 @@ int ParticleContentAnalysis::Process(HepMC::GenEvent* hepmcevt)
}
}
- if ( abs((*p) -> pdg_id()) == 2) m_uquark_mass -> Fill((*p) -> momentum().m(), weight);
- if ( abs((*p) -> pdg_id()) == 1) m_dquark_mass -> Fill((*p) -> momentum().m(), weight);
- if ( abs((*p) -> pdg_id()) == 4) m_cquark_mass -> Fill((*p) -> momentum().m(), weight);
- if ( abs((*p) -> pdg_id()) == 3) m_squark_mass -> Fill((*p) -> momentum().m(), weight);
- if ( abs((*p) -> pdg_id()) == 6) m_tquark_mass -> Fill((*p) -> momentum().m(), weight);
- if ( abs((*p) -> pdg_id()) == 5) m_bquark_mass -> Fill((*p) -> momentum().m(), weight);
- if ( (abs((*p) -> pdg_id()) == 21) || (abs((*p) -> pdg_id()) == 9) ) m_gluon_mass -> Fill((*p) -> momentum().m(), weight);
- if ( abs((*p) -> pdg_id()) == 22) m_photon_mass -> Fill((*p) -> momentum().m(), weight);
- if ( abs((*p) -> pdg_id()) == 23) m_Zboson_mass -> Fill((*p) -> momentum().m(), weight);
- if ( abs((*p) -> pdg_id()) == 24) m_Wboson_mass -> Fill((*p) -> momentum().m(), weight);
+ if ( std::abs((*p) -> pdg_id()) == 2) m_uquark_mass -> Fill((*p) -> momentum().m(), weight);
+ if ( std::abs((*p) -> pdg_id()) == 1) m_dquark_mass -> Fill((*p) -> momentum().m(), weight);
+ if ( std::abs((*p) -> pdg_id()) == 4) m_cquark_mass -> Fill((*p) -> momentum().m(), weight);
+ if ( std::abs((*p) -> pdg_id()) == 3) m_squark_mass -> Fill((*p) -> momentum().m(), weight);
+ if ( std::abs((*p) -> pdg_id()) == 6) m_tquark_mass -> Fill((*p) -> momentum().m(), weight);
+ if ( std::abs((*p) -> pdg_id()) == 5) m_bquark_mass -> Fill((*p) -> momentum().m(), weight);
+ if ( (std::abs((*p) -> pdg_id()) == 21) || (std::abs((*p) -> pdg_id()) == 9) ) m_gluon_mass -> Fill((*p) -> momentum().m(), weight);
+ if ( std::abs((*p) -> pdg_id()) == 22) m_photon_mass -> Fill((*p) -> momentum().m(), weight);
+ if ( std::abs((*p) -> pdg_id()) == 23) m_Zboson_mass -> Fill((*p) -> momentum().m(), weight);
+ if ( std::abs((*p) -> pdg_id()) == 24) m_Wboson_mass -> Fill((*p) -> momentum().m(), weight);
cmult++;
@@ -602,7 +606,7 @@ inline double ParticleContentAnalysis::getCosTheta(const HepMC::GenParticle* cha
//calculate theta defined as above
double cos_theta;
- cos_theta = cos (charged_pion_vector.angle(tau_vector));
+ cos_theta = std::cos (charged_pion_vector.angle(tau_vector));
return cos_theta;
}
diff --git a/Generators/HepMCAnalysis_i/src/PdfAnalysis.cxx b/Generators/HepMCAnalysis_i/src/PdfAnalysis.cxx
index 1c9814872e9d8788d44db53b63d64dbe7f766a23..9825a45a7a897b7641187ae47418c6e6b738d407 100644
--- a/Generators/HepMCAnalysis_i/src/PdfAnalysis.cxx
+++ b/Generators/HepMCAnalysis_i/src/PdfAnalysis.cxx
@@ -5,32 +5,13 @@
// ----------------------------------------------------------------------
#include <iostream>
-
-#include "AtlasHepMC/GenEvent.h"
-#include "AtlasHepMC/IO_GenEvent.h"
-#include "AtlasHepMC/GenParticle.h"
-#include "AtlasHepMC/GenVertex.h"
-#include "AtlasHepMC/IO_AsciiParticles.h"
-#include "AtlasHepMC/SimpleVector.h"
-#include "AtlasHepMC/WeightContainer.h"
-#include "CLHEP/Vector/LorentzVector.h"
-
-
#include "TH1.h"
#include "TMath.h"
-#include "TLorentzVector.h"
-
-#include "fastjet/PseudoJet.hh"
-#include "fastjet/ClusterSequence.hh"
-#include "fastjet/JetDefinition.hh"
-#include "fastjet/SISConePlugin.hh"
-
-
-#include "TruthUtils/HepMCHelpers.h"
-
+#include "MCUtils/HepMCEventUtils.h"
#include "../HepMCAnalysis_i/PdfAnalysis.h"
-
-using namespace std;
+namespace MC {
+using MCUtils::get_weight;
+}
// ----------------------------------------------------------------------
PdfAnalysis::PdfAnalysis()
diff --git a/Generators/HepMCAnalysis_i/src/PrepareHepMCAnalysisGenEvent.cxx b/Generators/HepMCAnalysis_i/src/PrepareHepMCAnalysisGenEvent.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..63f688bb1284d8e25be5f122e3bd32bc831e04ff
--- /dev/null
+++ b/Generators/HepMCAnalysis_i/src/PrepareHepMCAnalysisGenEvent.cxx
@@ -0,0 +1,22 @@
+/**
+ * Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ * Author: Andrii Verbytsky andrii.verbytskyi@mpp.mpg.de
+ *
+ * This file implements helper functions that prepare the ATLAS MC events
+ * to be used with the HepMCAnalysis package.
+ */
+#include <vector>
+#include "HepMC/GenEvent.h"
+//This tricky convention is needed for usage with HepMC3.
+using HepMCAnalysisGenEvent=HepMC::GenEvent;
+HepMC::GenEvent* PrepareHepMCAnalysisGenEvent(const HepMC::GenEvent* cevent)
+{
+//Note: the deep copy solves the issue described in the previous version.
+HepMC::GenEvent* event = new HepMC::GenEvent(*cevent);
+ for (auto p = event->particles_begin(); p != event->particles_end(); ++p) {
+ auto fv=(*p)->momentum();
+ fv.set(fv.x()*0.001,fv.y()*0.001,fv.z()*0.001,fv.e()*0.001);
+ (*p)->set_momentum(fv);
+ }
+ return event;
+}
diff --git a/Generators/HepMCAnalysis_i/src/UserAnalysis.cxx b/Generators/HepMCAnalysis_i/src/UserAnalysis.cxx
index 64c4867b2c31a756b665d44d5dd661889455d7c7..51f2714a2d9d66945751e1f4899a991c8025ddcf 100644
--- a/Generators/HepMCAnalysis_i/src/UserAnalysis.cxx
+++ b/Generators/HepMCAnalysis_i/src/UserAnalysis.cxx
@@ -6,17 +6,8 @@
#include <iostream>
-#include "AtlasHepMC/GenEvent.h"
-#include "AtlasHepMC/IO_GenEvent.h"
-#include "AtlasHepMC/GenParticle.h"
-#include "AtlasHepMC/GenVertex.h"
-#include "AtlasHepMC/IO_AsciiParticles.h"
-#include "AtlasHepMC/SimpleVector.h"
-#include "CLHEP/Vector/LorentzVector.h"
-
#include "TH1.h"
#include "TMath.h"
-#include "TLorentzVector.h"
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
@@ -25,7 +16,6 @@
#include "../HepMCAnalysis_i/UserAnalysis.h"
-using namespace std;
// ----------------------------------------------------------------------
UserAnalysis::UserAnalysis():
@@ -55,7 +45,7 @@ int UserAnalysis::Init(double maxEta, double minPt)
int UserAnalysis::Process(HepMC::GenEvent *event)
{
if (event->event_number() % 1000 == 0) {
- cout << "Processing event " << event->event_number() << endl;
+ std::cout << "Processing event " << event->event_number() << std::endl;
}
m_h_njets->Fill(m_inclusive_jets.size());