MuonSeedPtExtractor.cc

Go to the documentation of this file.

#include "RecoMuon/MuonSeedGenerator/src/MuonSeedPtExtractor.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "DataFormats/MuonDetId/interface/CSCDetId.h"
#include "DataFormats/MuonDetId/interface/ME0DetId.h"

#include "TMath.h"
#include <sstream>

MuonSeedPtExtractor::MuonSeedPtExtractor(const edm::ParameterSet& par)
: theBeamSpot(0.,0.,0.),
  scaleDT_( par.getParameter<bool>("scaleDT") )
{
  init(par);
}


void MuonSeedPtExtractor::init(const edm::ParameterSet& par)
{
  // load pT seed parameters
  // DT combinations
  fillParametersForCombo("DT_12", par);
  fillParametersForCombo("DT_13", par);
  fillParametersForCombo("DT_14", par);
  fillParametersForCombo("DT_23", par);
  fillParametersForCombo("DT_24", par);
  fillParametersForCombo("DT_34", par);
  // CSC combinations
  fillParametersForCombo("CSC_01", par);
  fillParametersForCombo("CSC_12", par);
  fillParametersForCombo("CSC_02", par);
  fillParametersForCombo("CSC_13", par);
  fillParametersForCombo("CSC_03", par);
  fillParametersForCombo("CSC_14", par);
  fillParametersForCombo("CSC_23", par);
  fillParametersForCombo("CSC_24", par);
  fillParametersForCombo("CSC_34", par);

  // Overlap combinations
  fillParametersForCombo("OL_1213", par);
  fillParametersForCombo("OL_1222", par);
  fillParametersForCombo("OL_1232", par);
  fillParametersForCombo("OL_2213", par);
  fillParametersForCombo("OL_2222", par);

  // Single segments (CSC)
  fillParametersForCombo("SME_11", par);
  fillParametersForCombo("SME_12", par);
  fillParametersForCombo("SME_13", par);
  fillParametersForCombo("SME_21", par);
  fillParametersForCombo("SME_22", par);
  fillParametersForCombo("SME_31", par);
  fillParametersForCombo("SME_32", par);
  fillParametersForCombo("SME_41", par);
  fillParametersForCombo("SME_42", par);

  // Single segments (DT)
  fillParametersForCombo("SMB_10", par);
  fillParametersForCombo("SMB_11", par);
  fillParametersForCombo("SMB_12", par);
  fillParametersForCombo("SMB_20", par);
  fillParametersForCombo("SMB_21", par);
  fillParametersForCombo("SMB_22", par);
  fillParametersForCombo("SMB_30", par);
  fillParametersForCombo("SMB_31", par);
  fillParametersForCombo("SMB_32", par);

  fillScalesForCombo("CSC_01_1_scale", par);
  fillScalesForCombo("CSC_12_1_scale", par);
  fillScalesForCombo("CSC_12_2_scale", par);
  fillScalesForCombo("CSC_12_3_scale", par);
  fillScalesForCombo("CSC_13_2_scale", par);
  fillScalesForCombo("CSC_13_3_scale", par);
  fillScalesForCombo("CSC_14_3_scale", par);
  fillScalesForCombo("CSC_23_1_scale", par);
  fillScalesForCombo("CSC_23_2_scale", par);
  fillScalesForCombo("CSC_24_1_scale", par);
  fillScalesForCombo("CSC_34_1_scale", par);
  fillScalesForCombo("DT_12_1_scale", par);
  fillScalesForCombo("DT_12_2_scale", par);
  fillScalesForCombo("DT_13_1_scale", par);
  fillScalesForCombo("DT_13_2_scale", par);
  fillScalesForCombo("DT_14_1_scale", par);
  fillScalesForCombo("DT_14_2_scale", par);
  fillScalesForCombo("DT_23_1_scale", par);
  fillScalesForCombo("DT_23_2_scale", par);
  fillScalesForCombo("DT_24_1_scale", par);
  fillScalesForCombo("DT_24_2_scale", par);
  fillScalesForCombo("DT_34_1_scale", par);
  fillScalesForCombo("DT_34_2_scale", par);
  fillScalesForCombo("OL_1213_0_scale", par);
  fillScalesForCombo("OL_1222_0_scale", par);
  fillScalesForCombo("OL_1232_0_scale", par);
  fillScalesForCombo("OL_2213_0_scale", par);
  fillScalesForCombo("OL_2222_0_scale", par);
  fillScalesForCombo("SMB_10_0_scale", par);
  fillScalesForCombo("SMB_11_0_scale", par);
  fillScalesForCombo("SMB_12_0_scale", par);
  fillScalesForCombo("SMB_20_0_scale", par);
  fillScalesForCombo("SMB_21_0_scale", par);
  fillScalesForCombo("SMB_22_0_scale", par);
  fillScalesForCombo("SMB_30_0_scale", par);
  fillScalesForCombo("SMB_31_0_scale", par);
  fillScalesForCombo("SMB_32_0_scale", par);
  fillScalesForCombo("SME_11_0_scale", par);
  fillScalesForCombo("SME_12_0_scale", par);
  fillScalesForCombo("SME_13_0_scale", par);
  fillScalesForCombo("SME_21_0_scale", par);
  fillScalesForCombo("SME_22_0_scale", par);

}


MuonSeedPtExtractor::~MuonSeedPtExtractor(){
}


void MuonSeedPtExtractor::fillParametersForCombo(const std::string & name, const edm::ParameterSet & pset)
{
  theParametersForCombo[name] = pset.getParameter<std::vector<double> >(name);
}


void MuonSeedPtExtractor::fillScalesForCombo(const std::string & name, const edm::ParameterSet & pset)
{
  theScalesForCombo[name] = pset.getParameter<std::vector<double> >(name);
}


std::vector<double> MuonSeedPtExtractor::pT_extract(MuonTransientTrackingRecHit::ConstMuonRecHitPointer firstHit,
                                                    MuonTransientTrackingRecHit::ConstMuonRecHitPointer secondHit) const
{
  GlobalPoint innerPoint = firstHit->globalPosition() - theBeamSpot;
  GlobalPoint outerPoint = secondHit->globalPosition() - theBeamSpot;
  MuonTransientTrackingRecHit::ConstMuonRecHitPointer innerHit = firstHit;
  MuonTransientTrackingRecHit::ConstMuonRecHitPointer outerHit = secondHit;

  // ways in which the hits could be in the wrong order
  if( (outerHit->isDT() && innerHit->isCSC())
   || (outerHit->isDT() && innerHit->isDT() && (innerPoint.perp() > outerPoint.perp())) 
   || (outerHit->isCSC() && innerHit->isCSC() && (fabs(innerPoint.z()) > fabs(outerPoint.z()))) )
  {
    innerHit = secondHit;
    outerHit = firstHit;
    innerPoint = innerHit->globalPosition() - theBeamSpot;
    outerPoint = outerHit->globalPosition() - theBeamSpot;
  } 
  
  double phiInner = innerPoint.phi();
  double phiOuter = outerPoint.phi();

  double etaInner = innerPoint.eta();
  double etaOuter = outerPoint.eta();
  //std::cout<<" inner pos = "<< innerPoint << " phi eta " << phiInner << " " << etaInner << std::endl;
  //std::cout<<" outer pos = "<< outerPoint << " phi eta " << phiOuter << " " << etaOuter << std::endl;
  //double thetaInner = firstHit->globalPosition().theta();
  // if some of the segments is missing r-phi measurement then we should
  // use only the 4D phi estimate (also use 4D eta estimate only)
  // the direction is not so important (it will be corrected) 
  /*
    bool firstOK = (4==allValidSets[iSet][firstMeasurement]->hit->dimension());
    bool lastOK = (4==allValidSets[iSet][lastMeasurement]->hit->dimension());
    if(!(firstOK * lastOK)){
    if(!firstOK){
    }
    if(!firstOK){
    }
    }
  */
  double dPhi =  phiInner - phiOuter;
  if(dPhi < -TMath::Pi()){
    dPhi += 2*TMath::Pi();
  }
  else if(dPhi >  TMath::Pi()){
    dPhi -= 2*TMath::Pi();
  }
  int sign = 1;
  if ( dPhi< 0.) {
    dPhi = -dPhi;
    sign = -1;
  }

  if (dPhi < 1.0e-6){
    dPhi = 1.0e-6;
  }
  double eta = fabs(etaOuter);// what if it is 2D segment? use inner?


  std::vector <int> stationCoded(2);
  stationCoded[0] = stationCoded[1] = 999;

  DetId  detId_inner = innerHit->hit()->geographicalId();
  DetId  detId_outer = outerHit->hit()->geographicalId();

  stationCoded[0] = stationCode(innerHit);
  stationCoded[1] = stationCode(outerHit);

  std::ostringstream os0 ;
  std::ostringstream os1 ;
  os0 << abs(stationCoded[0]);
  os1 << abs(stationCoded[1]);

  //std::cout<<" st1 = "<<stationCoded[0]<<" st2 = "<<stationCoded[1]<<std::endl;
  //std::cout<<" detId_inner = "<<detId_inner.rawId()<<"  detId_outer = "<< detId_outer.rawId()<<std::endl;
  std::string  combination = "0";
  std::string init_combination = combination;
  bool singleSegment = false;
  //if(detId_first == detId_second){
  if( stationCoded[0] == stationCoded[1]){
    // single segment - DT or CSC
    singleSegment = true;
    //eta = innerPoint.eta();
    GlobalVector gv = innerHit->globalDirection();
    double gvPerpPos = gv.x()*gv.x() + gv.y()*gv.y();
    if (gvPerpPos < 1e-32) gvPerpPos = 1e-32;
    gvPerpPos=sqrt(gvPerpPos);
    // Psi is angle between the segment origin and segment direction
    // Use dot product between two vectors to get Psi in global x-y plane
    double cosDpsi  = (gv.x()*innerPoint.x() + gv.y()*innerPoint.y());
    if (cosDpsi!=0){
      cosDpsi /= sqrt(innerPoint.x()*innerPoint.x() + innerPoint.y()*innerPoint.y());
      cosDpsi /= gvPerpPos;
      cosDpsi = cosDpsi > 1 ? 1 : cosDpsi;
      cosDpsi = cosDpsi < -1 ? -1 : cosDpsi;
    }

    double axb = ( innerPoint.x()*gv.y() ) - ( innerPoint.y()*gv.x() ) ;
    sign = (axb < 0.) ? 1 : -1;

    double dpsi = fabs(acos(cosDpsi)) ;
    if ( dpsi > TMath::Pi()/2.) {
      dpsi = TMath::Pi() - dpsi;
    }
    if (fabs(dpsi) < 0.00005) {
      dpsi = 0.00005;
    }
    dPhi = dpsi;

    if(innerHit->isDT())
    {
      DTChamberId dtCh(detId_inner);
      std::ostringstream os;
      os <<  dtCh.station() << abs(dtCh.wheel());
      combination = "SMB_"+os.str();
    }
    else if(innerHit->isCSC())
    {
      CSCDetId cscId(detId_inner);
      std::ostringstream os;
      int ring = cscId.ring();
      if(ring == 4) ring = 1;
      os << cscId.station() << ring;
      combination = "SME_"+os.str();
    }
    else if(innerHit->isME0())
    {
      ME0DetId me0Id(detId_inner);
      std::ostringstream os;
      int ring = 1;//me0Id.ring(); me0 only in ring 1
      os << me0Id.station() << ring;
      combination = "SME_"+os.str();
    }
    else
    {
      throw cms::Exception("MuonSeedPtExtractor") << "Bad hit DetId";
    }
  }
  else{
    if(stationCoded[0]<0){
      if(stationCoded[1]<0){
        // DT-DT
        combination = "DT_" + os0.str() + os1.str();
      }
      else{
        // DT-CSC
        eta = fabs(etaInner);
        if(-1==stationCoded[0]){
          switch (stationCoded[1]){
          case 1:
            combination = "OL_1213";
            break;
          case 2:
            combination = "OL_1222";
            break;
          case 3:
            combination = "OL_1232";
            break;
          default:
            // can not be 
            break;
          }
        }
        else if (-2==stationCoded[0]){
          if(1==stationCoded[1]){
            combination = "OL_2213";
          }
          else{
            // can not be (not coded?)
        combination = "OL_2222";// in case
          }
        }
        else{
          // can not be
        }
      }
    }
    else{
      if(stationCoded[1]<0){
        // CSC-DT
        // can not be
      }
      else{
        // CSC-CSC
        combination = "CSC_" + os0.str() + os1.str();
        if("CSC_04" == combination){
          combination = "CSC_14";
        }
      }
    }
  }

  std::vector<double> pTestimate(2);// 
  //std::cout<<" combination = "<<combination<<std::endl;
  if(init_combination!=combination){
        //std::cout<<" combination = "<<combination<<" eta = "<<eta<<" dPhi = "<<dPhi<<std::endl;
    ParametersMap::const_iterator parametersItr = theParametersForCombo.find(combination);
    if(parametersItr == theParametersForCombo.end()) {
      //      edm::LogWarning("RecoMuon|MuonSeedGenerator|MuonSeedPtExtractor") << "Cannot find parameters for combo " << combination;
      edm::LogWarning("BadSegmentCombination") << "Cannot find parameters for combo " << combination;
      pTestimate[0] = pTestimate[1] = 100;
      //       throw cms::Exception("MuonSeedPtEstimator") << "Cannot find parameters for combo " << combination;
    } else {

      if(scaleDT_ && outerHit->isDT() )
    {
      pTestimate = getPt(parametersItr->second, eta, dPhi, combination, detId_outer);
    }
      else 
    {
      pTestimate = getPt(parametersItr->second, eta, dPhi);
    }
      
      if(singleSegment){
    pTestimate[0] = fabs(pTestimate[0]);
    pTestimate[1] = fabs(pTestimate[1]);
      }
      pTestimate[0] *= double(sign);
    }
  }
  else{
    // often a MB3 - ME1/3 seed
    pTestimate[0] = pTestimate[1] = 100;
    // hmm
  }
  // std::cout<<" pTestimate[0] = "<<pTestimate[0]<<" pTestimate[1] = "<<pTestimate[1]<<std::endl;
  /*
                              MuonRecHitContainer_clusters[cluster][iHit+1]->isDT());
          if(specialCase){
            DTChamberId dtCh(detId);
            DTChamberId dtCh_2(detId_2);
            specialCase =  (dtCh.station() == dtCh_2.station());
          }
  */
  //return vPara;
  return pTestimate;
}


int MuonSeedPtExtractor::stationCode(MuonTransientTrackingRecHit::ConstMuonRecHitPointer hit) const
{
  DetId detId(hit->hit()->geographicalId());
  int result = -999;
  if(hit->isDT() ){
    DTChamberId dtCh(detId);
    //std::cout<<"first (DT) St/W/S = "<<dtCh.station()<<"/"<<dtCh.wheel()<<"/"<<dtCh.sector()<<"/"<<std::endl;
    result = -1 * dtCh.station();
  }
  else if( hit->isCSC() ){
    CSCDetId cscID(detId);
    //std::cout<<"first (CSC) E/S/R/C = "<<cscID.endcap()<<"/"<<cscID.station()<<"/"<<cscID.ring()<<"/"<<cscID.chamber()<<std::endl;
    result = cscID.station();
    if(result == 1 && (1 == cscID.ring() ||  4 == cscID.ring()) )
       result = 0;
  }
  else if( hit->isME0() ){
    result = 0;
  }
  else if(hit->isRPC()){
  }
  return result;
}


std::vector<double> MuonSeedPtExtractor::getPt(const std::vector<double> & vPara, double eta, double dPhi ) const {
   //std::cout<<" eta = "<<eta<<" dPhi = "<<dPhi<<" vPara[0] = "<<vPara[0]<<" vPara[1] = "<<vPara[1]<<" vPara[2] = "<<vPara[2]<<std::endl;
  double h  = fabs(eta);
  double estPt  = ( vPara[0] + vPara[1]*h + vPara[2]*h*h ) / dPhi;
  double estSPt = ( vPara[3] + vPara[4]*h + vPara[5]*h*h ) / dPhi;
  // std::cout<<"estPt = "<<estPt<<std::endl;
  std::vector<double> paraPt ;
  paraPt.push_back( estPt );
  paraPt.push_back( estSPt ) ;
  return paraPt ;
}

// combined pt parameterization and pt scale 
std::vector<double> MuonSeedPtExtractor::getPt(const std::vector<double> & vPara, double eta, double dPhi, 
     const std::string & combination, const DTChamberId & outerDetId ) const {
  double h  = fabs(eta);
  double estPt  = ( vPara[0] + vPara[1]*h + vPara[2]*h*h ) / dPhi;
  // changed by S.C.
  double estSPt = ( vPara[3] + vPara[4]*h + vPara[5]*h*h ) / dPhi;

  // scale the pt and spt , changed by S.C.
  int wheel = 0;
  if(combination[0] == 'D') {
    wheel = abs(outerDetId.wheel());
  }

  std::ostringstream os;
  os << combination << "_" << wheel << "_scale";

  ScalesMap::const_iterator scalesItr = theScalesForCombo.find(os.str());
  if ( scalesItr != theScalesForCombo.end()) {
    double t1 = scalesItr->second[0];
    double scaleFactor = 1./( 1. + t1/( estPt + 10. ) ) ;

    estSPt = estSPt * scaleFactor  ;
    estPt  = estPt  * scaleFactor  ;
  }

  // std::cout<<"estPt = "<<estPt<<std::endl;
  std::vector<double> paraPt ;
  paraPt.push_back( estPt );
  paraPt.push_back( estSPt ) ;
  return paraPt ;
}