TruthHitAnalysis.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
*/
 
#include "TruthHitAnalysis.h"
 
// Section of includes for Truth tests
#include "AtlasHepMC/GenEvent.h"
#include "TruthUtils/MagicNumbers.h"
 
 
 
StatusCode TruthHitAnalysis::initialize() {
  ATH_MSG_DEBUG( "Initializing TruthHitAnalysis" );
 
  // Grab the Ntuple and histogramming service for the tree
  ATH_CHECK(m_readKey.initialize());

  m_h_n_vert = new TH1D("h_n_vert","n_vert", 100,200, 1500);
  m_h_n_vert->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_n_vert->GetName(), m_h_n_vert));
 
  m_h_n_part = new TH1D("h_n_part","n_part", 100,1000, 10000);
  m_h_n_part->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_n_part->GetName(), m_h_n_part));
 
  m_h_n_vert_prim = new TH1D("h_n_vert_prim","n_vert prim", 100,0, 1000);
  m_h_n_vert_prim->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_n_vert_prim->GetName(), m_h_n_vert_prim));
 
  m_h_n_part_prim = new TH1D("h_n_part_prim","n_part prim", 100,200, 1500);
  m_h_n_part_prim->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_n_part_prim->GetName(), m_h_n_part_prim));
 
  m_h_n_vert_sec = new TH1D("h_n_vert_sec","n_vert sec", 100,0, 1000);
  m_h_n_vert_sec->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_n_vert_sec->GetName(), m_h_n_vert_sec));
 
  m_h_n_part_sec = new TH1D("h_n_part_sec","n_part sec", 100,0, 5000);
  m_h_n_part_sec->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_n_part_sec->GetName(), m_h_n_part_sec));
 
  m_h_vtx_x = new TH1D("h_vtx_x","vtx_x", 100,-1300, 1300);
  m_h_vtx_x->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_x->GetName(), m_h_vtx_x));
 
  m_h_vtx_y = new TH1D("h_vtx_y","vtx_y", 100,-1200, 1200);
  m_h_vtx_y->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_y->GetName(), m_h_vtx_y));
 
  m_h_vtx_z = new TH1D("h_vtx_z","vtx_z", 100,-5000, 5000);
  m_h_vtx_z->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_z->GetName(), m_h_vtx_z));
 
  m_h_vtx_r = new TH1D("h_vtx_r","vtx_r", 100,0, 1160);
  m_h_vtx_r->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_r->GetName(), m_h_vtx_r));
 
  m_h_vtx_prim_xy = new TH2D("h_vtx_prim_xy","vtx_prim_xy", 100,-100, 100, 100,-100, 100);
  m_h_vtx_prim_xy->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_prim_xy->GetName(), m_h_vtx_prim_xy));
 
  m_h_vtx_prim_zr = new TH2D("h_vtx_prim_zr","vtx_prim_zr", 100,-1500, 1500, 100,0, 150);
  m_h_vtx_prim_zr->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_prim_zr->GetName(), m_h_vtx_prim_zr));
 
  m_h_vtx_sec_xy = new TH2D("h_vtx_sec_xy","vtx_sec_xy", 100,-1200, 1200, 100,-1200, 1200);
  m_h_vtx_sec_xy->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_sec_xy->GetName(), m_h_vtx_sec_xy));
 
  m_h_vtx_sec_zr = new TH2D("h_vtx_sec_zr","vtx_sec_zr", 100,-6000, 6000, 100,0, 1160);
  m_h_vtx_sec_zr->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_sec_zr->GetName(), m_h_vtx_sec_zr));
 
  m_h_n_generations = new TH1D("h_n_generations","h_generations", 100,0, 25);
  m_h_n_generations->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_n_generations->GetName(), m_h_n_generations));
 
  m_h_truth_px = new TH1D("h_truth_px","truth_px", 100,0, 4000);
  m_h_truth_px->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_px->GetName(), m_h_truth_px));
 
  m_h_truth_py = new TH1D("h_truth_py","truth_py", 100,0, 4000);
  m_h_truth_py->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_py->GetName(), m_h_truth_py));
 
  m_h_truth_pz = new TH1D("h_truth_pz","truth_pz", 100,0, 4000);
  m_h_truth_pz->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_pz->GetName(), m_h_truth_pz));
 
  m_h_truth_pt = new TH1D("h_truth_pt","truth_pt", 100,0, 4000);
  m_h_truth_pt->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_pt->GetName(), m_h_truth_pt));
 
  m_h_truth_eta = new TH1D("h_truth_eta","truth_eta", 50,-10, 10);
  m_h_truth_eta->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_eta->GetName(), m_h_truth_eta));
 
  m_h_truth_phi = new TH1D("h_truth_phi","truth_phi", 25,-3.1416, 3.1416);
  m_h_truth_phi->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_phi->GetName(), m_h_truth_phi));
 
  m_h_barcode = new TH1D("h_truth_barcode","truth_barcode", 100,0, 300000);
  m_h_barcode->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_barcode->GetName(), m_h_barcode));
 
  m_h_part_status = new TH1D("h_part_status","part status", 100,0,50);
  m_h_part_status->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_part_status->GetName(), m_h_part_status));
 
  m_h_part_pdgid = new TH1D("h_part_pdgid","part pdgid", 100,-5000, 5000);
  m_h_part_pdgid->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_part_pdgid->GetName(), m_h_part_pdgid));
 
  m_h_part_pdgid_sec = new TH1D("h_part_pdgid_sec","part pdgid sec", 100,-5000, 5000);
  m_h_part_pdgid_sec->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_part_pdgid_sec->GetName(), m_h_part_pdgid_sec));
 
  m_h_part_eta = new TH1D("h_part_eta","part eta", 100,-10, 10);
  m_h_part_eta->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_part_eta->GetName(), m_h_part_eta));
 
  m_h_part_phi = new TH1D("h_part_phi","part phi", 100,-3.2, 3.2);
  m_h_part_phi->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_part_phi->GetName(), m_h_part_phi));
 
  m_h_part_p = new TH1D("h_part_p","part p", 100,0, 5000);
  m_h_part_p->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_part_p->GetName(), m_h_part_p));

  m_tree = new TTree("Truth", "Truth");
  std::string fullNtupleName =  "/" + m_ntupleFileName + "/";
  ATH_CHECK(histSvc()->regTree(fullNtupleName,m_tree));
  
  m_tree->Branch("vtx_x", &m_vtx_x);
  m_tree->Branch("vtx_y", &m_vtx_y);
  m_tree->Branch("vtx_z", &m_vtx_z);
  m_tree->Branch("vtx_r", &m_vtx_r);
  m_tree->Branch("vtx_barcode", &m_vtx_barcode);
  m_tree->Branch("truth_px", &m_truth_px);
  m_tree->Branch("truth_py", &m_truth_py);
  m_tree->Branch("truth_pz", &m_truth_pz);
  m_tree->Branch("truth_pt", &m_truth_pt);
  m_tree->Branch("truth_eta", &m_truth_eta);
  m_tree->Branch("truth_phi", &m_truth_phi);
  m_tree->Branch("barcode", &m_barcode);
  m_tree->Branch("status", &m_status);
  m_tree->Branch("pdg_id", &m_pdgid); 
 
  return StatusCode::SUCCESS;
}
 
 
StatusCode TruthHitAnalysis::execute() {
  ATH_MSG_DEBUG( "In TruthHitAnalysis::execute()" );
 
  m_vtx_x->clear();
  m_vtx_y->clear();
  m_vtx_z->clear();
  m_vtx_r->clear();
  m_vtx_barcode->clear();
  m_truth_px->clear();
  m_truth_py->clear();
  m_truth_pz->clear();
  m_truth_pt->clear();
  m_truth_eta->clear();
  m_truth_phi->clear();
  m_barcode->clear();
  m_status->clear();
  m_pdgid->clear();
  
  const EventContext& ctx{Gaudi::Hive::currentContext()};
  const McEventCollection* mcCollection{nullptr};
  ATH_CHECK(SG::get(mcCollection, m_readKey, ctx));
  
 
  McEventCollection::const_iterator currentGenEventIter = mcCollection->begin(); 
    if (currentGenEventIter != mcCollection->end()) {
      int nvtx = 0;
      int nvtx_sec=0;
#ifdef HEPMC3
    const auto &barcodes = (*currentGenEventIter)->attribute<HepMC::GenEventBarcodes> ("barcodes");
    std::map<int,int> id_to_barcode_map;
    if (barcodes) id_to_barcode_map = barcodes->id_to_barcode_map();
    for (const auto& vtx: (*currentGenEventIter)->vertices()) {
    int bcode = id_to_barcode_map[vtx->id()];
#else
    for (HepMC::GenEvent::vertex_const_iterator vtxit=(*currentGenEventIter)->vertices_begin(); vtxit!=(*currentGenEventIter)->vertices_end(); ++vtxit) {
    auto vtx=*vtxit;
        int bcode = HepMC::barcode(*vtx);
#endif
    double x = vtx->position().x();
    double y = vtx->position().y();
    double z = vtx->position().z();
    double r = std::sqrt(x*x+y*y);
    m_h_vtx_x->Fill(x);
    m_h_vtx_y->Fill(y);
    m_h_vtx_r->Fill(r);
    m_h_vtx_z->Fill(z);
 
    m_vtx_x->push_back(x);
    m_vtx_y->push_back(y);
    m_vtx_r->push_back(r);
    m_vtx_z->push_back(z);
    m_vtx_barcode->push_back(bcode);
 
    if (!HepMC::is_simulation_vertex(vtx)) {
      m_h_vtx_prim_xy->Fill(x,y);
      m_h_vtx_prim_zr->Fill(z,r);
      ++nvtx;
    }
    else {
      m_h_vtx_sec_xy->Fill(x,y);
      m_h_vtx_sec_zr->Fill(z,r);
      ++nvtx_sec; 
    }
      } //End iteration over vertices
 
      m_h_n_vert->Fill(nvtx+nvtx_sec);
      m_h_n_vert_prim->Fill(nvtx);
      m_h_n_vert_sec->Fill(nvtx_sec);
 
      int npart_prim=0; 
      int npart_sec=0;
 
      for (auto currentGenParticle: *(*currentGenEventIter)) {
    const HepMC::FourVector mom = currentGenParticle->momentum();
#ifdef HEPMC3
    int currentGenParticlebarcode = id_to_barcode_map[currentGenParticle->id()];
#else
        int currentGenParticlebarcode=HepMC::barcode(currentGenParticle);
#endif
    m_h_truth_px->Fill(mom.x());
    m_h_truth_py->Fill(mom.y());
    m_h_truth_pz->Fill(mom.z());
    m_h_truth_pt->Fill(mom.perp());
    m_h_truth_eta->Fill(mom.eta());
    m_h_truth_phi->Fill(mom.phi());
    m_h_barcode->Fill(currentGenParticlebarcode);
    m_h_part_status->Fill(currentGenParticle->status());
    m_truth_px->push_back(mom.x());
    m_truth_py->push_back(mom.y());
    m_truth_pz->push_back(mom.z());
    m_truth_pt->push_back(mom.perp());
    m_truth_eta->push_back(mom.eta());
    m_truth_phi->push_back(mom.phi());
    m_barcode->push_back(currentGenParticlebarcode);       
    m_status->push_back(currentGenParticle->status());
    
    int pdg = currentGenParticle->pdg_id();
    m_pdgid->push_back(pdg);
    
    if (!HepMC::is_simulation_particle(currentGenParticle)) {
      m_h_part_pdgid->Fill(pdg);
      m_h_part_p->Fill(std::sqrt(mom.x()*mom.x()+mom.y()*mom.y()+mom.z()*mom.z()));
      m_h_part_eta->Fill(mom.eta());
      m_h_part_phi->Fill(mom.phi());
      ++npart_prim; 
      m_h_n_generations->Fill(1);
/* The conditions don't make any sense.
      if (currentGenParticlebarcode < HepMC::PHOTOSMIN) {
        m_h_n_generations->Fill(0);
      }
      else {
        m_h_n_generations->Fill(1);
      }
*/
    } //End simulation particle
    else {
      m_h_part_pdgid_sec->Fill(pdg);
      ++npart_sec;
      const int gen = HepMC::generations(currentGenParticle) + 2; 
      m_h_n_generations->Fill(gen);  
    }
      } // End iteration over particles
 
      m_h_n_part_prim->Fill(npart_prim);
      m_h_n_part_sec->Fill(npart_sec);
      m_h_n_part->Fill(npart_prim+npart_sec); 
    } // End mcCollection 
 
  m_tree->Fill();
 
  return StatusCode::SUCCESS;
}