TrackRecordAnalysis.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
#include "TrackRecordAnalysis.h"
 
// Section of includes for TrackRecord tests
#include "TrackRecord/TrackRecord.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "StoreGate/ReadHandle.h"
 
 
StatusCode TrackRecordAnalysis::initialize() {
  ATH_MSG_DEBUG( "Initializing TrackRecordAnalysis" );
 
  std::string detName ="CaloEntry";
  std::string ntupName = "TrackRecordCaloEntry";
  ATH_CHECK(m_readKey.initialize());
  if (m_readKey.key() == "CaloEntryLayer") {
    detName="CaloEntry";
    ntupName="TrackRecordCaloEntry";
  }
  else if (m_readKey.key( )== "MuonEntryLayer") {
    detName="MuonEntry";
    ntupName="TrackRecordMuonEntry";
  }
  else if (m_readKey.key() == "MuonExitLayer") {
    detName="MuonExit";
    ntupName="TrackRecordMuonExit";
  }
  else {
    ATH_MSG_ERROR("TrackRecordAnalysis for "<< name() << "not supported !!! \n");
    return StatusCode::FAILURE;
  }
  
 
  float x_down = -5000;
  float x_up =  5000;
  float radius =  4500;
  float eta_down =  -5.8;
  float eta_up = 5.8;
  float z_down = -7000;
  float z_up = 7000;
  if (detName == "CaloEntry") {
    x_down = -1180;
    x_up =  1180;
    radius = 1200;
    eta_down = - 5.6;
    eta_up = 5.6;
    z_down = -3700;
    z_up = 3700;
  }
 
  m_h_hits_x = new TH1D((detName+"_x").c_str(),"hits_x", 100,x_down, x_up);
  m_h_hits_x->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_hits_x->GetName(), m_h_hits_x));
 
  m_h_hits_y = new TH1D((detName+"_y").c_str(), "hits_y", 100,x_down,x_up);
  m_h_hits_y->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_hits_y->GetName(), m_h_hits_y));
 
  m_h_hits_z = new TH1D((detName+"_z").c_str(), "hits_z", 100,z_down,z_up);
  m_h_hits_z->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_hits_z->GetName(), m_h_hits_z));
 
  m_h_hits_r = new TH1D((detName+"_r").c_str(), "hits_r", 100,0,radius);
  m_h_hits_r->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_hits_r->GetName(), m_h_hits_r));
 
  m_h_xy = new TH2D((detName+"_xy").c_str(), "xy", 100,x_down,x_up,100, x_down, x_up);
  m_h_xy->StatOverflows();
  ATH_CHECK(histSvc()->regHist( m_path+m_h_xy->GetName(), m_h_xy));
 
  m_h_zr = new TH2D((detName+"_zr").c_str(), "zr", 100,z_down,z_up,100, 0., radius);
  m_h_zr->StatOverflows();
  ATH_CHECK(histSvc()->regHist( m_path+m_h_zr->GetName(), m_h_zr));
 
  m_h_hits_eta = new TH1D((detName+"_eta").c_str(), "hits_eta", 100,eta_down,eta_up);
  m_h_hits_eta->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_hits_eta->GetName(), m_h_hits_eta));
 
  m_h_hits_phi = new TH1D((detName+"_phi").c_str(), "hits_phi", 100,-3.1416,3.1416);
  m_h_hits_phi->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_hits_phi->GetName(), m_h_hits_phi));
 
  m_h_hits_px = new TH1D((detName+"_px").c_str(),"Px", 100, -1500, 1500);
  m_h_hits_px->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_hits_px->GetName(), m_h_hits_px));
 
  m_h_hits_py = new TH1D((detName+"_py").c_str(), "Py", 100,-1500,1500);
  m_h_hits_py->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_hits_py->GetName(), m_h_hits_py));
 
  m_h_hits_pz = new TH1D((detName+"_pz").c_str(), "Pz", 100,-1500,1500);
  m_h_hits_pz->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_hits_pz->GetName(), m_h_hits_pz));
 
  m_h_hits_pt = new TH1D((detName+"_pt").c_str(), "hits_pt", 100,0,2500);
  m_h_hits_pt->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_hits_pt->GetName(), m_h_hits_pt));
 
  m_h_time = new TH1D((detName+"_time").c_str(),"time", 100,0, 250);
  m_h_time->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_time->GetName(), m_h_time));
 
  m_h_edep = new TH1D((detName+"_edep").c_str(), "energy", 100,0,7500);
  m_h_edep->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_edep->GetName(), m_h_edep));
 
  m_h_pdg = new TH1D((detName+"_pdg").c_str(), "pdg", 100,-1000,1000);
  m_h_pdg->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_pdg->GetName(), m_h_pdg));
 
  m_tree = new TTree(ntupName.c_str(), ntupName.c_str());
  std::string fullNtupleName =  "/" + m_ntupleFileName + "/" + detName;
  ATH_CHECK(histSvc()->regTree(fullNtupleName,m_tree));
    
  m_tree->Branch("x", &m_x);
  m_tree->Branch("y", &m_y);
  m_tree->Branch("z", &m_z);
  m_tree->Branch("r", &m_r);
  m_tree->Branch("eta", &m_eta);
  m_tree->Branch("phi", &m_phi);
  m_tree->Branch("px", &m_px);
  m_tree->Branch("py", &m_py);
  m_tree->Branch("pz", &m_pz);
  m_tree->Branch("pt", &m_pt);
  m_tree->Branch("time", &m_time);
  m_tree->Branch("energy", &m_edep);
  m_tree->Branch("pdg", &m_pdg);
  
  return StatusCode::SUCCESS;
}        
 
 
StatusCode TrackRecordAnalysis::execute() {
  ATH_MSG_DEBUG( "In TrackRecordAnalysis::execute()" );
 
  m_x->clear();
  m_y->clear();
  m_z->clear();
  m_r->clear();
  m_eta->clear();
  m_phi->clear();
  m_px->clear();
  m_py->clear();
  m_pz->clear();
  m_pt->clear();
  m_time->clear();
  m_edep->clear();
  m_pdg->clear();
  
  const EventContext& ctx{Gaudi::Hive::currentContext()};
  const TrackRecordCollection* TRcoll{nullptr};
  ATH_CHECK(SG::get(TRcoll,m_readKey ,ctx));
  for (TrackRecordCollection::const_iterator track = TRcoll->begin(); track != TRcoll->end(); ++track) {
      //TrackRecordCollection::const_iterator track;
      //for(auto track : *TRcoll){
      //std::cout<<"Entra en el loop"<<std::endl;
      
      CLHEP::Hep3Vector p =(*track).GetPosition();
      m_h_hits_x->Fill(p.x());
      m_h_hits_y->Fill(p.y());
      m_h_hits_z->Fill(p.z());
      m_h_hits_r->Fill(p.perp());
      m_h_xy->Fill(p.x(), p.y());
      m_h_zr->Fill(p.z(),p.perp());
      m_h_hits_eta->Fill(p.eta());
      m_h_hits_phi->Fill(p.phi());
      
      CLHEP::Hep3Vector mom = (*track).GetMomentum();
      m_h_hits_px->Fill(mom.x());
      m_h_hits_py->Fill(mom.y());
      m_h_hits_pz->Fill(mom.z());
      m_h_hits_pt->Fill(mom.perp());
      m_h_edep->Fill((*track).GetEnergy());
      m_h_time->Fill((*track).GetTime());
      m_h_pdg->Fill((*track).GetPDGCode());
  
      m_x->push_back(p.x());
      m_y->push_back(p.y());
      m_z->push_back(p.z());
      m_r->push_back(p.perp());
      m_eta->push_back(p.eta());
      m_phi->push_back(p.phi());
      m_px->push_back(mom.x());
      m_py->push_back(mom.y());
      m_pz->push_back(mom.z());
      m_pt->push_back(mom.perp());
      m_edep->push_back((*track).GetEnergy());
      m_time->push_back((*track).GetTime());
      m_pdg->push_back((*track).GetPDGCode());
  } // End while hits 
 
  m_tree->Fill();
 
  return StatusCode::SUCCESS;
}