#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;
   }
   
   
   // ============================================
   // ============================================
   /** @brief Contains resolution term values in a given eta range.
    *
    *  For a given muon 4-vector, constant and pT terms are used to
    *  compute a pT resolution.
    */
 
   
   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);
 
   }
   /** Calculates constant and pT terms to be used in resolution calculation */ 
   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;
       
     }    
   }
   /** Takes muon 4-vector and calculates pT resolution */
   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;
   }
  /** Dumps EtaBin properties to standard output */
   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;
   }
   
   // ============================================
   /** @brief Contains all EtaBins for a given range in phi.
    *
    *  Functions called on a PhiBin are delegated to the relevant
    *  EtaBin.
    */
   
   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();
   }
    /** Sets phi limits on this bin  */  
   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);
     
   }
   /** Delegates combine method to relevant EtaBin */
   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);
     }
     
   }
   /** Delegates calculateResolution method to relevant EtaBin */
   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);    
     
   }
   /** Dumps PhiBin properties to standard output and calls dump for
     *  all constituent EtaBins */
   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));
   }
 
   // ============================================
   // ============================================
   /** @brief Contains PhiBins for a given spectrometer sector
    *
    *  Functions called on a Sector are delegated to the relevant
    *  PhiBin.
    */
 
   /** Constructor takes xml DOMNodes for detector boundaries and sector index */
   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();
   }
   /** Delegates combine method to relevant PhiBin */
   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);
     }
   }
   
   /** Delegates calculateResolution method to relevant PhiBin.
      *  Phi values is translated to local system */
   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);    
     
   }
   /** Dumps Sector properties to standard output and calls dump
     *  for all constituent PhiBins*/
   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));
   }
 
   // ============================================
   /** @brief Spectrometer class containing exact resolutions at given pT values
    *
    *  These are combined (via the combine method) to produce MuonSpectrometers
    *  which contain the constant and pT terms necessary to compute a resolution
    *  at any pT.
    */
 
 
   /** Constructor takes xml DOMNodes for detector boundaries and pt resolutions and
    *  indicators for whether this is the lowest pT or highest pT ProtoSpectrometer. */
   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;
     
   }
 
   // ============================================
 
 /** @brief Contains collection of spectrometer Sectors.
    *
    *  Functions called on a MuonSpectrometer are delegated to the relevant
    *  Sector.
    */
 
   /** Constructor takes upper and lower pT values for range covered and
      *  constituent Sectors */
   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();
   }
   /** Delegates calculateResolution method to relevant Sector */
   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);    
 
   }
   /** Dumps MuonSpectrometer properties to standard output and calls
      *  dump for all constituent Sectors */
   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));
   }
 
   // ============================================
 
   /** The constructor takes the two main nodes in the input xml file
      *  and uses them to instantiate some ProtoSpectrometers. The
      *  ProtoSpectrometers contain exact resolution values that are
      *  later used to extract pT and constant terms per eta and phi
      *  range (represented by EtaBin and PhiBin respectively).
      *
      *  The two main nodes in the input xml file are:
      *
      *   <interpretation>   Contains information to define the spatial
      *                      boundaries on components of the muon spectrometer
      *
      *   <resolutions>      Contains the resolutions relevant to the
      *                      components defined in <interpretation>
      */
 
 
   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
 

Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems
This page was automatically generated by the LXR engine.