MuonSegmentPlots.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
 
#include "MuonHistUtils/MuonSegmentPlots.h"
#include "xAODTracking/TrackingPrimitives.h"
#include "MuonStationIndex/MuonStationIndex.h"
#include "GeoPrimitives/GeoPrimitives.h"
#include <cmath>
 
namespace Muon{
 
using ChIdx = MuonStationIndex::ChIndex;
constexpr int chIdxMax = static_cast<int>(ChIdx::ChIndexMax);
 
MuonSegmentPlots::MuonSegmentPlots(PlotBase* pParent, const std::string& sDir): PlotBase(pParent, sDir) {  
 
  //booking histograms
  segmentfitChi2 = Book1D("segmentfitChi2", "Segment fit #chi^{2};#chi^{2};Entries", 100,0,200);
  segmentfitNdof = Book1D("segmentfitNdof", "Segment fit N_{dof};N_{dof};Entries", 100,0,20);
  segmentfitChi2oNdof = Book1D("segmentfitChi2oNdof", "Segment fit #chi^{2}/N_{dof};Segment fit #chi^{2}/N_{dof};Entries", 100,0,20);
  
  t0 = Book1D("t0","Segment Refit t_{0};t_{0};Entries",100,-25,25);
  t0_top = Book1D("t0_top","Segment Refit t_{0} for y>0;t_{0};Entries",100,-25,25);
  t0_bottom = Book1D("t0_bottom","Segment Refit t_{0} for y<0;t_{0};Entries",100,-25,25);
  t0err = Book1D("t0err","Segment Refit t_{0} error;t_{0} error;Entries",100,0,10);
  t0err_top = Book1D("t0err_top","Segment Refit t_{0} error for y>0;t_{0} error;Entries",100,0,10);
  t0err_bottom = Book1D("t0err_bottom","Segment Refit t_{0} error for y<0;t_{0} error;Entries",100,0,10);
 
  nPrecisionHits = Book1D("nPrecisionHits","Segment precision hits;hits;Entries",20,0,20);
  nPhiLayers = Book1D("nPhiLayers","Segment phi layers;#phi layers;Entries",10,0,10);
  nTrigEtaLayers = Book1D("nTrigEtaLayers","Segment eta trigger layers;#eta trigger layers;Entries",10,0,10);
  nPrecisionHits_nTriggerHits = Book2D("nPrecisionHits_nTriggerHits", "Number of MDT hits vs Tigger station hits; MDT hits; Trigger hits", 20, -0.5, 19.5, 20, -0.5, 19.5);
 
  etaIndex = Book1D("etaIndex","Segment #eta Index ;#eta index;Entries",21,-10.5,10.5);
  sector = Book1D("sector","Segment phi sector;#phi sector;Entries",16,0.5,16.5);
  
 
  xypos_barrel = Book2D("xypos_barrel","Segment position x-y, barrel;x_{pos};y_{pos}",150,-14000,14000,150,-14000,14000);
  xypos_endcap = Book2D("xypos_endcap","Segment position x-y, endcap;x_{pos};y_{pos}",150,-14000,14000,150,-14000,14000);
  rzpos_sectorLarge = Book2D("rzpos_sectorLarge","Segment position r-z, large sectors normalized by solid angle;z_{pos};r_{pos}",100,-23000,23000,75,0,14000);
 
  rzpos_sectorSmall = Book2D("rzpos_sectorSmall","Segment position r-z, small sectors normalized by solid angle;z_{pos};r_{pos}",100,-23000,23000,75,0,14000);
  
 
  etadir = Book1D("etadir","Segment pointing direction eta;#eta_{dir};Entries",100,-5,5);
  etadir_barrel = Book1D("etadir_barrel","Segment pointing direction eta, barrel;#eta_{dir};Entries",100,-5,5);
  etadir_endcap = Book1D("etadir_endcap","Segment pointing direction eta, endcap;#eta_{dir};Entries",100,-5,5);
  phidir = Book1D("phidir","Segment pointing direction phi;#phi_{dir};Entries",64,-3.2,3.2);
  etaphidir = Book2D("etaphidir","Segment pointing direction phi vs eta;#eta_{dir};#phi_{dir}",64,-3.2,3.2,64,-3.2,3.2);
 
  chamberIndex = Book1D("chamberIndex","Chamber index; Chamber Index",chIdxMax,0,chIdxMax);
  chamberIndex_perSector = Book2D("chamberIndex_perSector","Number of Segments per Chamber, normalized by solid angle; Sector; Chamber Index ", 33, -16.5, 16.5, chIdxMax,0,chIdxMax);
  eff_chamberIndex_perSector_numerator = Book2D("eff_chamberIndex_perSector_numerator","Number of expected hits for Segments per Chamber; Sector; Chamber Index ", 33, -16.5, 16.5, chIdxMax,0,chIdxMax);
  eff_chamberIndex_perSector_denominator = Book2D("eff_chamberIndex_perSector_denominator","Number of recorded precision hits for Segments per Chamber; Sector; Chamber Index ", 33, -16.5, 16.5, chIdxMax,0,chIdxMax);
  eff_chamberIndex_perSector = Book2D("eff_chamberIndex_perSector","precision layer hit efficiency per chamber; Sector; Chamber Index ", 33, -16.5, 16.5, chIdxMax,0,chIdxMax);
  //chamberIndex_dtheta = Book2D("chamberIndex_dtheta","Segment #Delta#theta between position and momentum; #Delta#theta; Chamber Index ", 180, -90.0, 90.0, chIdxMax,0,chIdxMax);
  for (int i=1; i<=chamberIndex->GetXaxis()->GetNbins(); i++){
    const char *temp_chambername = MuonStationIndex::chName(static_cast<MuonStationIndex::ChIndex>(chamberIndex->GetBinLowEdge(i))).c_str();
    chamberIndex->GetXaxis()->SetBinLabel(i, temp_chambername);
    chamberIndex_perSector->GetYaxis()->SetBinLabel(i, temp_chambername);
    eff_chamberIndex_perSector_numerator->GetYaxis()->SetBinLabel(i, temp_chambername);
    eff_chamberIndex_perSector_denominator->GetYaxis()->SetBinLabel(i, temp_chambername);
    eff_chamberIndex_perSector->GetYaxis()->SetBinLabel(i, temp_chambername);
 
    sector_etaIndex.push_back(Book2D(Form("%s_etastation", temp_chambername), Form("Number of Segment in %s; #phi Sector; #eta Index", temp_chambername), 18, -0.5, 17.5, 19, -9.5, 9.5));
    sector_etaIndex_nPrechit.push_back(Book2D(Form("%s_etastation_nPrechit", temp_chambername), Form("Number of Segment Prec hit in %s; #phi Sector; #eta Index", temp_chambername), 18, -0.5, 17.5, 19, -9.5, 9.5));
    sector_etaIndex_nTrighit.push_back(Book2D(Form("%s_etastation_nTrighit", temp_chambername), Form("Number of Segment Phi + Trigeta hit in %s; #phi Sector; #eta Index", temp_chambername), 18, -0.5, 17.5, 19, -9.5, 9.5));
    eff_sector_etaIndex_nPrechit.push_back(Book2D(Form("eff_%s_etastation_nPrechit", temp_chambername), Form("Segment Prec hit eff in %s; #phi Sector; #eta Index", temp_chambername), 18, -0.5, 17.5, 19, -9.5, 9.5));
    eff_sector_etaIndex_nTrighit.push_back(Book2D(Form("eff_%s_etastation_nTrighit", temp_chambername), Form("Segment Phi + Trigeta hit eff in %s; #phi Sector; #eta Index", temp_chambername), 18, -0.5, 17.5, 19, -9.5, 9.5));
  }
}
 
MuonSegmentPlots::~MuonSegmentPlots()  = default;
 
  void MuonSegmentPlots::fill(const xAOD::MuonSegment& muSeg, float weight)
{
  float chi2 = muSeg.chiSquared();
  float ndof = muSeg.numberDoF();
  segmentfitChi2->Fill(chi2,weight);
  segmentfitNdof->Fill(ndof,weight);
  if (ndof>0) segmentfitChi2oNdof->Fill(muSeg.chiSquared()/muSeg.numberDoF(), weight);
 
  float x=muSeg.x();
  float y=muSeg.y();
 
  float segt0 = muSeg.t0();
  float segt0err = muSeg.t0error();
 
  t0->Fill(segt0,weight);
  t0err->Fill(segt0err,weight);
  if (y>0) { 
    t0_top->Fill(segt0,weight);
    t0err_top->Fill(segt0err,weight);
  } else {
    t0_bottom->Fill(segt0,weight);
    t0err_bottom->Fill(segt0err,weight);
  }
 
  sector->Fill(muSeg.sector(),weight);
  etaIndex->Fill(muSeg.etaIndex(),weight);
 
  nPrecisionHits->Fill(muSeg.nPrecisionHits(),weight);
  nPhiLayers->Fill(muSeg.nPhiLayers(),weight);
  nTrigEtaLayers->Fill(muSeg.nTrigEtaLayers(),weight);
  // not sure how to implement weights here JEF 8/4/2021
  nPrecisionHits_nTriggerHits->Fill(muSeg.nPrecisionHits(), muSeg.nPhiLayers() + muSeg.nTrigEtaLayers(),weight); 
 
 
  // not sure how to implement weights here for these chamber/sector Index plots JEF 8/4/2021
  int chIndex = static_cast<int>(muSeg.chamberIndex());
  float chambernorm = 1/Chamberarea[chIndex];//weight of the segment using the chamber eta-phi area
  chamberIndex->Fill(chIndex,weight);
  int sectorIndex = muSeg.sector();
  int etaIndex = muSeg.etaIndex();
  //fill the count of segments; switch the sign here to make the plots
  if (muSeg.z() < 0) { sectorIndex = - sectorIndex;}
  chamberIndex_perSector->Fill(sectorIndex, chIndex, chambernorm*weight);
  eff_chamberIndex_perSector_numerator->Fill(sectorIndex, chIndex, weight*((muSeg.nPrecisionHits() > Chamberexpectedhits[chIndex]) ? Chamberexpectedhits[chIndex]:muSeg.nPrecisionHits()));
  eff_chamberIndex_perSector_denominator->Fill(sectorIndex, chIndex, weight*Chamberexpectedhits[chIndex]);
  //update sector eta index plots; switch the sign back here to make the plots
  sectorIndex = muSeg.sector(); 
  sector_etaIndex[chIndex]->Fill(sectorIndex, etaIndex);//for weighted average
  //double currentfill = sector_etaIndex[chIndex]->GetBinContent(sector_etaIndex[chIndex]->GetXaxis()->FindBin(sectorIndex), sector_etaIndex[chIndex]->GetYaxis()->FindBin(etaIndex));
  sector_etaIndex_nPrechit[chIndex]->Fill(sectorIndex, etaIndex, weight*muSeg.nPrecisionHits());
  sector_etaIndex_nTrighit[chIndex]->Fill(sectorIndex, etaIndex, weight*(muSeg.nPhiLayers() + muSeg.nTrigEtaLayers()));
  eff_sector_etaIndex_nPrechit[chIndex]->Fill(sectorIndex, etaIndex, weight*Chamberexpectedhits[chIndex]);
  eff_sector_etaIndex_nTrighit[chIndex]->Fill(sectorIndex, etaIndex, weight*Chamberexpectedtrighits[chIndex]);
  
  bool isBarrel = (chIndex<static_cast<int>(ChIdx::BEE))? true: false; // BEE -> endcap
  bool isSectorLarge = ( (isBarrel && chIndex%2==1) || (!isBarrel && chIndex%2==0 && chIndex!=static_cast<int>(ChIdx::BEE)) )? true : false; 
 
  
  //position and direction plots
  const Amg::Vector3D globalPos(muSeg.position());
  const Amg::Vector3D globalDir(muSeg.direction());
  
  //protect against cases with no hit information or segments parallell to the beam!
  if ( globalDir.mag() <= FLT_EPSILON  || globalDir.perp() <= FLT_EPSILON) return;
  //set up direction vectors
  float r = globalPos.perp();
  float z = globalPos.z();
  //fill the rz plots
  if (isSectorLarge) {rzpos_sectorLarge->Fill(z,r, chambernorm*weight);}
  else {rzpos_sectorSmall->Fill(z,r, chambernorm*weight);}
 
  float eta = globalDir.eta();
  float phi = globalDir.phi();
  etadir->Fill(eta,weight);
  phidir->Fill(phi,weight);
  etaphidir->Fill(eta,phi);
 
  
  if (isBarrel) {
    xypos_barrel->Fill(x,y, chambernorm*weight);
    etadir_barrel->Fill(eta,weight);
  } else {
    xypos_endcap->Fill(x,y, chambernorm*weight);
    etadir_endcap->Fill(eta,weight);
  }
 
}
 
} // closing namespace Muon