sTGCHitAnalysis.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
*/
 
#include "sTGCHitAnalysis.h"
 
 
/*Includes specific to sTGCs*/
#include "GeoAdaptors/GeoMuonHits.h"
#include "MuonReadoutGeometry/MuonDetectorManager.h"
#include "MuonReadoutGeometry/sTgcReadoutElement.h"
#include "MuonSimEvent/sTgcHitIdHelper.h"
#include "MuonSimEvent/sTGCSimHit.h"
#include "CLHEP/Vector/LorentzVector.h"
 
#include "TH1.h"
#include "TTree.h"
#include "TString.h"
 
#include <algorithm>
#include <math.h>
#include <functional>
#include <iostream>
#include <stdio.h>
 
 
StatusCode sTGCHitAnalysis::initialize() {
  ATH_MSG_DEBUG( "Initializing sTGCHitAnalysis" );
 
  // Grab the Ntuple and histogramming service for the tree
  ATH_CHECK(m_readKey.initialize());
 
 
  m_h_rz_A = new TH2D("h_stgc_rz_A","rz_A", 2000,6500., 8500., 1000, -5000., 5000.);
  m_h_rz_A->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_rz_A->GetName(), m_h_rz_A));
 
  m_h_rz_C = new TH2D("h_stgc_rz_C","rz_C", 2000,-8500., -6500., 5000, 0., 5000.);
  m_h_rz_C->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_rz_C->GetName(), m_h_rz_C));
 
  m_h_L_rz_A = new TH2D("h_stgc_L_rz_A","L_rz_A", 2000, 6500., 8500., 5000, 0., 5000.);
  m_h_L_rz_A->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_L_rz_A->GetName(), m_h_L_rz_A));
 
  m_h_L_rz_C = new TH2D("h_stgc_L_rz_C","L_rz_C", 2000,-8500., -6500., 5000, 0., 5000.);
  m_h_L_rz_C->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_L_rz_C->GetName(), m_h_L_rz_C));
 
  m_h_S_rz_A = new TH2D("h_stgc_S_rz_A","S_rz_A", 2000,6500., 8500., 5000, 0., 5000.);
  m_h_S_rz_A->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_S_rz_A->GetName(), m_h_S_rz_A));
 
  m_h_S_rz_C = new TH2D("h_stgc_S_rz_C","S_rz_C", 2000, -8500., -6500., 5000, 0., 5000.);
  m_h_S_rz_C->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_S_rz_C->GetName(), m_h_S_rz_C));
 
  
  m_h_xy_A = new TH2D("h_stgc_xy_A","xy_A", 1200,-6000., 6000., 1200, -6000., 6000.);
  m_h_xy_A->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_xy_A->GetName(), m_h_xy_A));
 
  m_h_xy_C = new TH2D("h_stgc_xy_C","xy_C", 1200,-6000., 6000., 1200, -6000., 6000.);
  m_h_xy_C->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_xy_C->GetName(), m_h_xy_C));
 
  m_h_L_xy_A = new TH2D("h_stgc_L_xy_A","L_xy_A", 1200,-6000., 6000., 1200, -6000., 6000.);
  m_h_L_xy_A->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_L_xy_A->GetName(), m_h_L_xy_A));
 
  m_h_L_xy_C = new TH2D("h_stgc_L_xy_C","L_xy_C", 1200,-6000., 6000., 1200, -6000., 6000.);
  m_h_L_xy_C->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_L_xy_C->GetName(), m_h_L_xy_C));
 
  m_h_S_xy_A = new TH2D("h_stgc_S_xy_A","S_xy_A", 1200,-6000., 6000., 1200, -6000., 6000.);
  m_h_S_xy_A->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_S_xy_A->GetName(), m_h_S_xy_A));
 
  m_h_S_xy_C = new TH2D("h_stgc_S_xy_C","S_xy_C", 1200,-6000., 6000., 1200, -6000., 6000.);
  m_h_S_xy_C->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_S_xy_C->GetName(), m_h_S_xy_C));
 
  m_h_r_A = new TH1D("h_stgc_r_A","r_A", 1000,0., 14000.);
  m_h_r_A->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_r_A->GetName(), m_h_r_A));
  
  m_h_r_C = new TH1D("h_stgc_r_C","r_C", 1000,0., 14000.);
  m_h_r_C->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_r_C->GetName(), m_h_r_C));
 
  
  return StatusCode::SUCCESS;
}
 
StatusCode sTGCHitAnalysis::execute() {
  ATH_MSG_DEBUG( "In sTGCHitAnalysis::execute()" );
 
  const EventContext& ctx{Gaudi::Hive::currentContext()};
  const sTGCSimHitCollection* p_collection{nullptr};
  ATH_CHECK(SG::get(p_collection, m_readKey, ctx));
  for (sTGCSimHitCollection::const_iterator i_hit = p_collection->begin(); i_hit != p_collection->end(); ++i_hit){
 
      Amg::Vector3D p = (*i_hit).globalPosition();
      //Get station names and make plots for each wedge
      const sTgcHitIdHelper* hitHelper = sTgcHitIdHelper::GetHelper();
      int simId = (*i_hit).sTGCId();
      std::string sim_stationName = hitHelper->GetStationName(simId);
 
      static const std::string QS1C("QS1C");
      static const std::string QS2C("QS2C");
      static const std::string QS3C("QS3C");
      static const std::string QL1P("QL1P");
      static const std::string QL2P("QL2P");
      static const std::string QL3P("QL3P");
 
      int sim_side = hitHelper->GetSide(simId);
 
      if ( sim_side==1){ //Both sectors A side
        double r_A = std::sqrt(p.x()*p.x()+p.y()*p.y()); //Evaluate r 
        m_h_rz_A->Fill(p.z(), r_A);
          m_h_xy_A->Fill(p.x(), p.y());
          m_h_r_A->Fill(r_A);
      }
 
      if ( sim_side==-1){ //Both sectors C side
        double r_C = std::sqrt(p.x()*p.x()+p.y()*p.y()); //Evaluate r
        m_h_rz_C->Fill(p.z(), r_C);
          m_h_xy_C->Fill(p.x(), p.y());
          m_h_r_C->Fill(r_C);
      }
 
      if ((sim_stationName==QS1C || sim_stationName==QS2C || sim_stationName==QS3C) && sim_side==1){ //Small sectors A side 
        double r_S_A = std::sqrt(p.x()*p.x()+p.y()*p.y()); //Evaluate r
        m_h_S_rz_A->Fill(p.z(), r_S_A);
          m_h_S_xy_A->Fill(p.x(), p.y());
      }
 
      if ((sim_stationName==QS1C || sim_stationName==QS2C || sim_stationName==QS3C) && sim_side==-1){ //Small sectors C side 
        double r_S_C = std::sqrt(p.x()*p.x()+p.y()*p.y()); //Evaluate r
        m_h_S_rz_C->Fill(p.z(), r_S_C);
        m_h_S_xy_C->Fill(p.x(), p.y());
      }
 
      if ((sim_stationName==QL1P || sim_stationName==QL2P || sim_stationName==QL3P) && sim_side==1){ //Large sectors A side
        double r_L_A = std::sqrt(p.x()*p.x()+p.y()*p.y()); //Evaluate r
        m_h_L_rz_A->Fill(p.z(), r_L_A);
        m_h_L_xy_A->Fill(p.x(), p.y());
      }
 
      if ((sim_stationName==QL1P || sim_stationName==QL2P || sim_stationName==QL3P) && sim_side==-1){ //Large sectors C side
        double r_L_C = std::sqrt(p.x()*p.x()+p.y()*p.y()); //Evaluate r
        m_h_L_rz_C->Fill(p.z(), r_L_C);
        m_h_L_xy_C->Fill(p.x(), p.y());
      }
 
  }
  return StatusCode::SUCCESS;
}