MuonSpectrometer.cxx

Go to the documentation of this file.

#include <iostream>
#include <cmath>
#include <cstdlib>

#include "AtlfastAlgs/MuonSpectrometer.h"

namespace Atlfast{

  // ============================================

  double dEnergyLoss(double eta, double pT){
    
    // Contribution from energy loss fluctuations
    // in the colorimeters
    // Unitless fraction
    
    eta = fabs(eta);
    double theta = 2.*atan(exp(-eta));
    double sintheta = sin(theta);
    double energy = fabs(pT)/sintheta;
    double factor = (eta < 1.1) ? 1./sqrt(sintheta) : sqrt(14./11.);
    double dEloss = (eta < 1.1) ? 
      factor*(240. + 0.105*energy/100.)/energy : 
      factor*(480. + 0.105*energy/100.)/energy;

    return dEloss;

  }
  
  double dEnergyLoss(const HepLorentzVector& avec){
    return dEnergyLoss(avec.pseudoRapidity(),avec.perp());
  }
  
  // ============================================
  // Predicates
  // ============================================
  
  IsCorrectEtaBin::IsCorrectEtaBin(double eta):
    m_eta(eta){}
  bool IsCorrectEtaBin::operator()(EtaBin& eb){
    // returns true if eta value is within this EtaBin
    return ( m_eta >= eb.getEtaMin() && m_eta < eb.getEtaMax() ) ? true : false;
  }
  
  // ============================================
  
  IsCorrectPhiBin::IsCorrectPhiBin(double phi):
    m_phi(phi){}
  bool IsCorrectPhiBin::operator()(PhiBin& pb){
    // returns true if phi value is within this PhiBin
    return ( m_phi >= pb.getPhiMin() && m_phi < pb.getPhiMax() ) ? true : false;
  }
  
  // ============================================
  
  IsCorrectSector::IsCorrectSector(std::string sectortype):
    m_sectortype(sectortype){}
  bool IsCorrectSector::operator()(Sector& sec){
    // returns true if this Sector is of the correct type
    return ( m_sectortype == sec.getSectorType() ) ? true : false;
  }
  
  // ============================================
  
  IsCorrectMuonSpectrometer::IsCorrectMuonSpectrometer(double pT):
    m_pT(pT){}
  bool IsCorrectMuonSpectrometer::operator()(MuonSpectrometer& ms){
    // returns true if pt value is within this MuonSpectrometer
    return ( m_pT >= ms.getPtMin() && m_pT < ms.getPtMax() ) ? true : false;
  }
  
  
  // ============================================
  
  EtaBin::EtaBin(double deltaponp, int etabin):
    m_etabin(etabin), m_etamin(0.), m_etamax(0.), m_deltaponp(0.01*deltaponp){}
  
  void EtaBin::setEtaStuff(double etaFirst, double deltaEta){
    
    m_etavalue = deltaEta*double(m_etabin);
    m_etamin = m_etabin ? deltaEta*double(m_etabin-0.5) + etaFirst: etaFirst;
    m_etamax = deltaEta*double(m_etabin+0.5);

  }
  
  void EtaBin::combine(EtaBin& otherEtaBin, double thispt, double nextpt){

    
    m_constterm = 1.;
    m_ptterm = 0;

    if (m_deltaponp < 1.){
      
      double resthis = m_deltaponp;
      double resnext = otherEtaBin.m_deltaponp;
      double delossthis = dEnergyLoss(m_etavalue,thispt);
      double delossnext = dEnergyLoss(m_etavalue,nextpt);
      
      m_ptterm = ( resnext*resnext - resthis*resthis
                   - (delossnext*delossnext - delossthis*delossthis ) ) /
        (nextpt*nextpt - thispt*thispt);
      
      m_constterm = resthis*resthis - m_ptterm*thispt*thispt - delossthis*delossthis;
      
    }    
  }
  
  double EtaBin::calculateResolution(const HepLorentzVector& avec){

    double resolution = 1.;
    if (m_constterm < 1.){ // Otherwise return resolution of 1.0
      double dEloss = dEnergyLoss(avec);
      double pT = avec.perp();
      double resolution_squared = m_constterm + m_ptterm*pT*pT + dEloss*dEloss;
      if (resolution_squared > 0. && resolution_squared < 1.)
        resolution = sqrt(resolution_squared);
    }
    return resolution;
  }

  void EtaBin::dump(std::string indent){
    indent += " ";
    std::cout<<std::endl;
    std::cout<<indent<<"\nDump for EtaBin:\n";
    std::cout<<indent<<"m_etabin "<<m_etabin<<std::endl;
    std::cout<<indent<<"m_etamin "<<m_etamin<<std::endl;
    std::cout<<indent<<"m_etamax "<<m_etamax<<std::endl;
    std::cout<<indent<<"m_deltaponp "<<m_deltaponp<<std::endl;
    std::cout<<indent<<"m_constterm "<<m_constterm<<std::endl;
    std::cout<<indent<<"m_ptterm "<<m_ptterm<<std::endl;
  }
  
  // ============================================
  
  PhiBin::PhiBin(DOMNode* boundaries, DOMNode* phi):
    m_phibin(0), m_phimin(0.), m_phimax(0.){
    
    // phibin goes from 0 to 90
    getFirstValue(phi, "phibinnumber", m_phibin);
    
    std::vector<double> deltaponps;
    getAllValues(phi, "deltaponp",deltaponps);
    SpectrometerComponentBinInstantiator<EtaBin> instantiator(boundaries);
    m_etas = std::for_each(deltaponps.begin(),
                           deltaponps.end(),             
                           instantiator).objects();
    
    double etafirst;
    double etalast;
    
    DOMNode* etaInfoNode =  getFirstChildTagByName(boundaries, "etainfo");
    getFirstValue(etaInfoNode, "etafirst", etafirst);
    getFirstValue(etaInfoNode, "etalast", etalast);
    
    int etaPoints = m_etas.size();
    double deltaEta = (etalast-etafirst)/(etaPoints-1);
    
    std::vector<EtaBin>::iterator itr = m_etas.begin();
    for(;itr != m_etas.end(); ++itr){
      (*itr).setEtaStuff(etafirst,deltaEta);
    }
    
  }
  

  PhiBin::~PhiBin(){
    m_etas.clear();
  }
  
  void PhiBin::setPhiStuff(double sectorPhimin, double sectorPhimax, double deltaPhi, int phipoints){

    // these are phi values within a sector    
    m_phimin = m_phibin ? deltaPhi*double(m_phibin-0.5) + sectorPhimin : sectorPhimin;
    if (m_phimin > M_PI) m_phimin -= (2.*M_PI);
    
    m_phimax = (m_phibin != (phipoints-1)) ? deltaPhi*double(m_phibin+0.5) : sectorPhimax;
    if (m_phimax > M_PI) m_phimax -= (2.*M_PI);
    
  }
  
  void PhiBin::combine(PhiBin &otherPhiBin, double thispt, double nextpt){
    
    std::vector<EtaBin>::iterator thisEtaItr = m_etas.begin();
    std::vector<EtaBin>::iterator otherEtaItr = otherPhiBin.m_etas.begin();
    for(;thisEtaItr != m_etas.end(); ++thisEtaItr, ++otherEtaItr){

      (*thisEtaItr).combine(*otherEtaItr, thispt, nextpt);
    }
    
  }

  double PhiBin::calculateResolution(const HepLorentzVector& avec){
    
    IsCorrectEtaBin iceb(fabs(avec.pseudoRapidity()));
    std::vector<EtaBin>::iterator ebitr = find_if(m_etas.begin(),
                                                  m_etas.end(),
                                                  iceb);
    if ( ebitr == m_etas.end() ){
      // Muon outside eta range, returning resolution of 1.0
      return 1.;
    }
    
    return (*ebitr).calculateResolution(avec);    
    
  }

  void PhiBin::dump(std::string indent){
    indent += "  ";
    std::cout<<indent<<"Dump for PhiBin:\n";
    std::cout<<indent<<"m_phibin "<<m_phibin<<std::endl;
    std::cout<<indent<<"m_phimin "<<m_phimin<<std::endl;
    std::cout<<indent<<"m_phimax "<<m_phimax<<std::endl;
    std::cout<<indent<<"Dump of my etas:"<<std::endl;
    std::for_each(m_etas.begin(), m_etas.end(), AddToDump(indent));
  }

  // ============================================

  Sector::Sector(DOMNode* boundaries, DOMNode* sector):m_phimin(0.){

    m_sectorType = "Unknown";
    
    getFirstValue(sector, "sectortype", m_sectorType);
    
    std::vector<DOMNode*> phis = getAllChildTagsByName(sector, "phibinres");
    SpectrometerComponentInstantiator<PhiBin> instantiator(boundaries);
    m_phis = std::for_each(phis.begin(),
                           phis.end(),
                           instantiator).objects();
    
    DOMNode* phiInfoNode =  getFirstChildTagByName(boundaries, "phiinfo");
    getFirstValue(phiInfoNode, "phifirst", m_phimin);
    getFirstValue(phiInfoNode, "philast", m_phimax);
    
    // Calculate phi info once here and give it to the PhiBins
    
    int phipoints = m_phis.size();
    double deltaphi = (m_phimax-m_phimin)/(phipoints-1);


    std::vector<PhiBin>::iterator itr = m_phis.begin();
    for(;itr != m_phis.end(); ++itr){
      (*itr).setPhiStuff(m_phimin, m_phimax, deltaphi, phipoints);
    }
    
  }
  
  Sector::~Sector(){
    m_phis.clear();
  }

  void Sector::combine(Sector& otherSec, double thispt, double nextpt){
    
    std::vector<PhiBin>::iterator thisPhiItr = m_phis.begin();
    std::vector<PhiBin>::iterator otherPhiItr = otherSec.m_phis.begin();
    for(;thisPhiItr != m_phis.end(); ++thisPhiItr, ++otherPhiItr){
      (*thisPhiItr).combine(*otherPhiItr, thispt, nextpt);
    }
  }
  
  double Sector::calculateResolution(const HepLorentzVector& avec){

    // check for appropriate PhiBin
    // First calculate the local phi (zeroed at the sector phimin)
    
    double localphi = avec.phi();
    while (localphi < 0.) localphi += 2.*M_PI;
    double global_phi = localphi;
    double sectorsize = m_phimax - m_phimin;
    while (localphi > sectorsize ) localphi -= sectorsize;
    if ( global_phi >= 3.926990815 && global_phi < 4.712388978 )
      localphi = sectorsize - localphi;
    
    IsCorrectPhiBin icpb(localphi);
    std::vector<PhiBin>::iterator pbitr = find_if(m_phis.begin(),
                                                  m_phis.end(),
                                                  icpb);
    if ( pbitr == m_phis.end() ){
      // Could not find correct PhiBin, returning resolution of 1.0
      return 1.;
    }
    
    return (*pbitr).calculateResolution(avec);    
    
  }

  void Sector::dump(std::string indent){
    indent += "  ";
    std::cout<<indent<<"Dump for Sector:\n";
    std::cout<<indent<<"m_sectorType "<<m_sectorType<<std::endl;
    std::cout<<indent<<"m_phimin     "<<m_phimin<<std::endl;
    std::for_each(m_phis.begin(), m_phis.end(), AddToDump(indent));
  }

  // ============================================
  
  ProtoSpectrometer::ProtoSpectrometer(DOMNode* boundaries, DOMNode* ptbinres, bool lowestpT, bool highestpT):
    m_lowestpT(lowestpT), m_highestpT(highestpT){
    
    getFirstValue(ptbinres, "ptbinnumber", m_ptbin);
    
    DOMNode* node = getFirstChildTagByName(boundaries, "ptpoints");
    
    std::vector<DOMNode*> ptBoundaryNodes = 
      getAllChildTagsByName(node, "ptpoint");
    
    std::vector<double> ptBoundaries = 
      for_each(
               ptBoundaryNodes.begin(),
               ptBoundaryNodes.end(),
               GetDoublesFromNodes()
               ).doubles();
    
    if(m_ptbin < static_cast<int>(ptBoundaries.size())){
      m_ptvalue = ptBoundaries[m_ptbin];
    }
    std::vector<DOMNode*> sectors = getAllChildTagsByName(ptbinres, "sector");
    SpectrometerComponentInstantiator<Sector> instantiator(boundaries);
    m_sectors = std::for_each(sectors.begin(),
                              sectors.end(),
                              instantiator).objects();
  }
  
  ProtoSpectrometer::~ProtoSpectrometer(){
    m_sectors.clear();
  }

  void ProtoSpectrometer::dump(std::string indent){
    indent += "  ";
    std::cout<<indent<<"Dump for Spectrometer:\n";
    std::cout<<indent<<"m_bin "<<m_ptbin<<std::endl;
    std::cout<<indent<<"m_ptvalue "<<m_ptvalue<<std::endl;
    std::cout<<indent<<"Dump of my sectors:"<<std::endl;
    std::for_each(m_sectors.begin(), m_sectors.end(), AddToDump(indent));
  }

  MuonSpectrometer ProtoSpectrometer::combine(ProtoSpectrometer& otherPS){

    std::vector<Sector>::iterator thisSecItr = m_sectors.begin();
    std::vector<Sector>::iterator otherSecItr = otherPS.m_sectors.begin();
    
    for(; thisSecItr!= m_sectors.end(); ++thisSecItr, ++otherSecItr){
      (*thisSecItr).combine(*otherSecItr, m_ptvalue, otherPS.m_ptvalue);
    }
    
    double ptlow = m_lowestpT ? 0.0 : m_ptvalue;
    double pthigh = otherPS.m_highestpT ? 7000000.0 : otherPS.m_ptvalue;
    
    MuonSpectrometer thisMS(ptlow, pthigh, m_sectors);
    
    return thisMS;
    
  }

  // ============================================

  MuonSpectrometer::MuonSpectrometer(double ptmin, double ptmax, std::vector<Sector> sectors):
    m_ptmin(ptmin), m_ptmax(ptmax), m_sectors(sectors){
  }

  void MuonSpectrometer::setSectorMap(std::map<int,std::string> sectortypes){
    m_sectortypes = sectortypes;
    m_sectorsize = 2.*M_PI/m_sectortypes.size();
  }
  
  double MuonSpectrometer::calculateResolution(const HepLorentzVector& avec){

    // Choose sector number based on phi
    // -pi < HepLorentzVector.phi() < +pi
    double phi = avec.phi();
    while (phi < 0.) phi += 2.*M_PI;
    // Is this compiler-dependent rounding-up? Must check
    int sectorindex = int(phi/m_sectorsize) + 1;

    // Translate this to a sector type
    std::string sectortype = m_sectortypes[sectorindex];
    
    // Now find the corresponding sector
    IsCorrectSector ics(sectortype);
    std::vector<Sector>::iterator secitr = find_if(m_sectors.begin(),
                                                    m_sectors.end(),
                                                    ics);
    if ( secitr == m_sectors.end() ){
      // Could not find correct Sector, returning resolution of 1.0
      return 1.;
    }
    
    return (*secitr).calculateResolution(avec);    

  }

  void MuonSpectrometer::dump(std::string indent){
    indent += "  ";
    std::cout<<indent<<"Dump for MuonSpectrometer:\n";
    std::cout<<indent<<"m_ptmin "<<m_ptmin<<std::endl;
    std::cout<<indent<<"m_ptmax "<<m_ptmax<<std::endl;
    std::cout<<indent<<"Dump of my sectors:"<<std::endl;
    std::for_each(m_sectors.begin(), m_sectors.end(), AddToDump(indent));
  }

  // ============================================

  MuonResolutionCalculator::MuonResolutionCalculator(const char* filename){
    
    DOMDocument* dom = parseFileForDocument(filename);
    
    if(dom == 0){
      std::cout<<"Could not parse " <<filename<<std::endl;
    }
    
    DOMElement* docElement = dom->getDocumentElement();
    
    if(docElement == 0){
      std::cout<<"Could not get DOMElement from DOMDocument"<<std::endl;
    }
    
    // get the node containing information for interpreting the resolutions
    XMLCh* s = XMLString::transcode("interpretation");
    DOMNodeList* binTags = 
      docElement->getElementsByTagName(s);
    XMLString::release(&s);

    DOMNode* boundaries = binTags->item(0);
    
    // get the node needed to build a spectrometer 
    s = XMLString::transcode("ptbinres");
    DOMNodeList* ptbinresTags = 
      docElement->getElementsByTagName(s);
    XMLString::release(&s);

    // build the spectrometers
    
    for(XMLSize_t ptr=0; ptr<ptbinresTags->getLength();++ptr){
      
      ProtoSpectrometer thisPS(boundaries,
                               ptbinresTags->item(ptr),
                               ptr==0,
                               ptr==(ptbinresTags->getLength()-1));
      
      m_protoSpectrometers.push_back(thisPS);
      
    }
    
    std::map<int,std::string> sectorMap = makeSectorMap(boundaries);
    
    std::vector<ProtoSpectrometer>::iterator psItr = m_protoSpectrometers.begin();
    for (; psItr!=m_protoSpectrometers.end()-1; ++psItr){
      
      MuonSpectrometer newMS = (*psItr).combine(*(psItr+1));
      newMS.setSectorMap(sectorMap);
      m_muonSpectrometers.push_back(newMS);
      
    }
    
    // Free up memory from the ProtoSpectrometers as they're no longer needed
    m_protoSpectrometers.clear();

    // Delete the DOMDocument
    delete dom;

    // Clean up xerces-style. This should deallocate all memory used by Xerces
    // but doesn't seem to work. Nice one.
    XMLPlatformUtils::Terminate();

  }
  
  std::map<int,std::string> MuonResolutionCalculator::makeSectorMap(DOMNode* boundaries){
    


    DOMNode* node = getFirstChildTagByName(boundaries, "sectortypes");

    std::vector<DOMNode*> sectorTypeNodes =
      getAllChildTagsByName(node, "sectortype");

    AddToMap atm;
    std::map<int,std::string> sectorMap = std::for_each(sectorTypeNodes.begin(),
                                                        sectorTypeNodes.end(),
                                                        atm).getMap();
    return sectorMap;
    
  }

  double MuonResolutionCalculator::calculateResolution(const HepLorentzVector& avec){

    // check for appropriate MuonSpectrometer
    IsCorrectMuonSpectrometer icms(avec.perp());
    std::vector<MuonSpectrometer>::iterator msitr = find_if(m_muonSpectrometers.begin(),
                                                             m_muonSpectrometers.end(),
                                                             icms);
    if ( msitr == m_muonSpectrometers.end() ){
      // Could not find correct MuonSpectrometer, returning resolution of 1.0
      return 1.;
    }
    
    return (*msitr).calculateResolution(avec);    
  }

  void MuonResolutionCalculator::dump(std::string indent){
    indent += " ";
    std::cout<<indent<<"\nDump for MuonResolutionCalculator:\n";
    std::for_each(m_protoSpectrometers.begin(), 
                  m_protoSpectrometers.end(), 
                  AddToDump(indent));
    std::for_each(m_muonSpectrometers.begin(), 
                  m_muonSpectrometers.end(), 
                  AddToDump(indent));
    
  }      
  
} // namespace Atlfast

---