MuonSeedCreator.cc

Go to the documentation of this file.

#include "RecoMuon/MuonSeedGenerator/src/MuonSeedCreator.h"

#include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"

#include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"

#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"

#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/DetLayers/interface/DetLayer.h"

#include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
#include "DataFormats/Common/interface/OwnVector.h"
#include "DataFormats/MuonDetId/interface/DTChamberId.h"
#include "DataFormats/MuonDetId/interface/CSCDetId.h"
#include "DataFormats/MuonDetId/interface/RPCDetId.h"

#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "gsl/gsl_statistics.h"

/*
 * Constructor
 */
MuonSeedCreator::MuonSeedCreator(const edm::ParameterSet& pset) {
  theMinMomentum = pset.getParameter<double>("minimumSeedPt");
  theMaxMomentum = pset.getParameter<double>("maximumSeedPt");
  defaultMomentum = pset.getParameter<double>("defaultSeedPt");
  debug = pset.getParameter<bool>("DebugMuonSeed");
  sysError = pset.getParameter<double>("SeedPtSystematics");
  // load seed PT parameters
  DT12 = pset.getParameter<std::vector<double> >("DT_12");
  DT13 = pset.getParameter<std::vector<double> >("DT_13");
  DT14 = pset.getParameter<std::vector<double> >("DT_14");
  DT23 = pset.getParameter<std::vector<double> >("DT_23");
  DT24 = pset.getParameter<std::vector<double> >("DT_24");
  DT34 = pset.getParameter<std::vector<double> >("DT_34");

  CSC01 = pset.getParameter<std::vector<double> >("CSC_01");
  CSC12 = pset.getParameter<std::vector<double> >("CSC_12");
  CSC02 = pset.getParameter<std::vector<double> >("CSC_02");
  CSC13 = pset.getParameter<std::vector<double> >("CSC_13");
  CSC03 = pset.getParameter<std::vector<double> >("CSC_03");
  CSC14 = pset.getParameter<std::vector<double> >("CSC_14");
  CSC23 = pset.getParameter<std::vector<double> >("CSC_23");
  CSC24 = pset.getParameter<std::vector<double> >("CSC_24");
  CSC34 = pset.getParameter<std::vector<double> >("CSC_34");

  OL1213 = pset.getParameter<std::vector<double> >("OL_1213");
  OL1222 = pset.getParameter<std::vector<double> >("OL_1222");
  OL1232 = pset.getParameter<std::vector<double> >("OL_1232");
  OL1213 = pset.getParameter<std::vector<double> >("OL_1213");
  OL2222 = pset.getParameter<std::vector<double> >("OL_1222");

  SME11 = pset.getParameter<std::vector<double> >("SME_11");
  SME12 = pset.getParameter<std::vector<double> >("SME_12");
  SME13 = pset.getParameter<std::vector<double> >("SME_13");
  SME21 = pset.getParameter<std::vector<double> >("SME_21");
  SME22 = pset.getParameter<std::vector<double> >("SME_22");
  SME31 = pset.getParameter<std::vector<double> >("SME_31");
  SME32 = pset.getParameter<std::vector<double> >("SME_32");
  SME41 = pset.getParameter<std::vector<double> >("SME_41");

  SMB10 = pset.getParameter<std::vector<double> >("SMB_10");
  SMB11 = pset.getParameter<std::vector<double> >("SMB_11");
  SMB12 = pset.getParameter<std::vector<double> >("SMB_12");
  SMB20 = pset.getParameter<std::vector<double> >("SMB_20");
  SMB21 = pset.getParameter<std::vector<double> >("SMB_21");
  SMB22 = pset.getParameter<std::vector<double> >("SMB_22");
  SMB30 = pset.getParameter<std::vector<double> >("SMB_30");
  SMB31 = pset.getParameter<std::vector<double> >("SMB_31");
  SMB32 = pset.getParameter<std::vector<double> >("SMB_32");

  // Load dphi scale parameters
  CSC01_1 = pset.getParameter<std::vector<double> >("CSC_01_1_scale");
  CSC12_1 = pset.getParameter<std::vector<double> >("CSC_12_1_scale");
  CSC12_2 = pset.getParameter<std::vector<double> >("CSC_12_2_scale");
  CSC12_3 = pset.getParameter<std::vector<double> >("CSC_12_3_scale");
  CSC13_2 = pset.getParameter<std::vector<double> >("CSC_13_2_scale");
  CSC13_3 = pset.getParameter<std::vector<double> >("CSC_13_3_scale");
  CSC14_3 = pset.getParameter<std::vector<double> >("CSC_14_3_scale");
  CSC23_1 = pset.getParameter<std::vector<double> >("CSC_23_1_scale");
  CSC23_2 = pset.getParameter<std::vector<double> >("CSC_23_2_scale");
  CSC24_1 = pset.getParameter<std::vector<double> >("CSC_24_1_scale");
  CSC34_1 = pset.getParameter<std::vector<double> >("CSC_34_1_scale");

  DT12_1 = pset.getParameter<std::vector<double> >("DT_12_1_scale");
  DT12_2 = pset.getParameter<std::vector<double> >("DT_12_2_scale");
  DT13_1 = pset.getParameter<std::vector<double> >("DT_13_1_scale");
  DT13_2 = pset.getParameter<std::vector<double> >("DT_13_2_scale");
  DT14_1 = pset.getParameter<std::vector<double> >("DT_14_1_scale");
  DT14_2 = pset.getParameter<std::vector<double> >("DT_14_2_scale");
  DT23_1 = pset.getParameter<std::vector<double> >("DT_23_1_scale");
  DT23_2 = pset.getParameter<std::vector<double> >("DT_23_2_scale");
  DT24_1 = pset.getParameter<std::vector<double> >("DT_24_1_scale");
  DT24_2 = pset.getParameter<std::vector<double> >("DT_24_2_scale");
  DT34_1 = pset.getParameter<std::vector<double> >("DT_34_1_scale");
  DT34_2 = pset.getParameter<std::vector<double> >("DT_34_2_scale");

  OL_1213 = pset.getParameter<std::vector<double> >("OL_1213_0_scale");
  OL_1222 = pset.getParameter<std::vector<double> >("OL_1222_0_scale");
  OL_1232 = pset.getParameter<std::vector<double> >("OL_1232_0_scale");
  OL_2213 = pset.getParameter<std::vector<double> >("OL_2213_0_scale");
  OL_2222 = pset.getParameter<std::vector<double> >("OL_2222_0_scale");

  SMB_10S = pset.getParameter<std::vector<double> >("SMB_10_0_scale");
  SMB_11S = pset.getParameter<std::vector<double> >("SMB_11_0_scale");
  SMB_12S = pset.getParameter<std::vector<double> >("SMB_12_0_scale");
  SMB_20S = pset.getParameter<std::vector<double> >("SMB_20_0_scale");
  SMB_21S = pset.getParameter<std::vector<double> >("SMB_21_0_scale");
  SMB_22S = pset.getParameter<std::vector<double> >("SMB_22_0_scale");
  SMB_30S = pset.getParameter<std::vector<double> >("SMB_30_0_scale");
  SMB_31S = pset.getParameter<std::vector<double> >("SMB_31_0_scale");
  SMB_32S = pset.getParameter<std::vector<double> >("SMB_32_0_scale");

  SME_11S = pset.getParameter<std::vector<double> >("SME_11_0_scale");
  SME_12S = pset.getParameter<std::vector<double> >("SME_12_0_scale");
  SME_13S = pset.getParameter<std::vector<double> >("SME_13_0_scale");
  SME_21S = pset.getParameter<std::vector<double> >("SME_21_0_scale");
  SME_22S = pset.getParameter<std::vector<double> >("SME_22_0_scale");
}

/*
 * Destructor
 */
MuonSeedCreator::~MuonSeedCreator() {}

/*
 * createSeed
 *
 * Note type = 1 --> CSC
 *           = 2 --> Overlap
 *           = 3 --> DT
 */
TrajectorySeed MuonSeedCreator::createSeed(
    int type, const SegmentContainer& seg, const std::vector<int>& layers, int NShowers, int NShowerSegments) {
  // The index of the station closest to the IP
  int last = 0;

  double ptmean = theMinMomentum;
  double sptmean = theMinMomentum;

  AlgebraicVector t(4);
  AlgebraicSymMatrix mat(5, 0);
  LocalPoint segPos;

  // Compute the pt according to station types used;
  if (type == 1)
    estimatePtCSC(seg, layers, ptmean, sptmean);
  if (type == 2)
    estimatePtOverlap(seg, layers, ptmean, sptmean);
  if (type == 3)
    estimatePtDT(seg, layers, ptmean, sptmean);
  if (type == 4)
    estimatePtSingle(seg, layers, ptmean, sptmean);
  // type 5 are the seeding for ME1/4
  if (type == 5)
    estimatePtCSC(seg, layers, ptmean, sptmean);

  // for certain clear showering case, set-up the minimum value
  if (NShowers > 0)
    estimatePtShowering(NShowers, NShowerSegments, ptmean, sptmean);
  //if ( NShowers > 0 ) std::cout<<" Showering happened "<<NShowers<<" times w/ "<< NShowerSegments<<std::endl; ;

  // Minimal pt
  double charge = 1.0;
  if (ptmean < 0.)
    charge = -1.0;
  if ((charge * ptmean) < theMinMomentum) {
    ptmean = theMinMomentum * charge;
    sptmean = theMinMomentum;
  } else if ((charge * ptmean) > theMaxMomentum) {
    ptmean = theMaxMomentum * charge;
    sptmean = theMaxMomentum * 0.25;
  }

  LocalTrajectoryParameters param;

  double p_err = 0.0;
  // determine the seed layer
  int best_seg = 0;
  double chi2_dof = 9999.0;
  unsigned int ini_seg = 0;
  // avoid generating seed from  1st layer(ME1/1)
  if (type == 5)
    ini_seg = 1;
  for (size_t i = ini_seg; i < seg.size(); i++) {
    double dof = static_cast<double>(seg[i]->degreesOfFreedom());
    if (chi2_dof < (seg[i]->chi2() / dof))
      continue;
    chi2_dof = seg[i]->chi2() / dof;
    best_seg = static_cast<int>(i);
  }

  if (type == 1 || type == 5 || type == 4) {
    // Fill the LocalTrajectoryParameters
    last = best_seg;
    // last = 0;
    GlobalVector mom = seg[last]->globalPosition() - GlobalPoint();
    segPos = seg[last]->localPosition();

    GlobalVector polar(GlobalVector::Spherical(mom.theta(), seg[last]->globalDirection().phi(), 1.));

    polar *= fabs(ptmean) / polar.perp();
    LocalVector segDirFromPos = seg[last]->det()->toLocal(polar);
    int chargeI = static_cast<int>(charge);
    LocalTrajectoryParameters param1(segPos, segDirFromPos, chargeI);
    param = param1;
    p_err = (sptmean * sptmean) / (polar.mag() * polar.mag() * ptmean * ptmean);
    mat = seg[last]->parametersError().similarityT(seg[last]->projectionMatrix());
    mat[0][0] = p_err;
    if (type == 5) {
      mat[0][0] = mat[0][0] / fabs(tan(mom.theta()));
      mat[1][1] = mat[1][1] / fabs(tan(mom.theta()));
      mat[3][3] = 2.25 * mat[3][3];
      mat[4][4] = 2.25 * mat[4][4];
    }
    if (type == 4) {
      mat[0][0] = mat[0][0] / fabs(tan(mom.theta()));
      mat[1][1] = mat[1][1] / fabs(tan(mom.theta()));
      mat[2][2] = 2.25 * mat[2][2];
      mat[3][3] = 2.25 * mat[3][3];
      mat[4][4] = 2.25 * mat[4][4];
    }
    double dh = fabs(seg[last]->globalPosition().eta()) - 1.6;
    if (fabs(dh) < 0.1 && type == 1) {
      mat[1][1] = 4. * mat[1][1];
      mat[2][2] = 4. * mat[2][2];
      mat[3][3] = 9. * mat[3][3];
      mat[4][4] = 9. * mat[4][4];
    }

    //if ( !highPt && type != 1 ) mat[1][1]= 2.25*mat[1][1];
    //if (  highPt && type != 1 ) mat[3][3]= 2.25*mat[1][1];
    //mat[2][2]= 3.*mat[2][2];
    //mat[3][3]= 2.*mat[3][3];
    //mat[4][4]= 2.*mat[4][4];
  } else {
    // Fill the LocalTrajectoryParameters
    last = 0;
    segPos = seg[last]->localPosition();
    GlobalVector mom = seg[last]->globalPosition() - GlobalPoint();
    GlobalVector polar(GlobalVector::Spherical(mom.theta(), seg[last]->globalDirection().phi(), 1.));
    //GlobalVector polar(GlobalVector::Spherical(seg[last]->globalDirection().theta(),seg[last]->globalDirection().phi(),1.));

    //double ptRatio = 1.0 - (2.808/(fabs(ptmean) -1)) + (4.546/( (fabs(ptmean)-1)*(fabs(ptmean)-1)) );
    //ptmean = ptmean*ptRatio ;

    polar *= fabs(ptmean) / polar.perp();
    LocalVector segDirFromPos = seg[last]->det()->toLocal(polar);
    int chargeI = static_cast<int>(charge);
    LocalTrajectoryParameters param1(segPos, segDirFromPos, chargeI);
    param = param1;
    p_err = (sptmean * sptmean) / (polar.mag() * polar.mag() * ptmean * ptmean);
    mat = seg[last]->parametersError().similarityT(seg[last]->projectionMatrix());
    //mat[0][0]= 1.44 * p_err;
    mat[0][0] = p_err;
  }

  if (debug) {
    GlobalPoint gp = seg[last]->globalPosition();
    float Geta = gp.eta();

    std::cout << "Type " << type << " Nsegments " << layers.size() << " ";
    std::cout << "pt " << ptmean << " errpt    " << sptmean << " eta " << Geta << " charge " << charge << std::endl;
  }

  // this perform H.T() * parErr * H, which is the projection of the
  // the measurement error (rechit rf) to the state error (TSOS rf)
  // Legend:
  // H => is the 4x5 projection matrix
  // parError the 4x4 parameter error matrix of the RecHit

  LocalTrajectoryError error(asSMatrix<5>(mat));

  // Create the TrajectoryStateOnSurface
  TrajectoryStateOnSurface tsos(param, error, seg[last]->det()->surface(), &*BField);

  // Take the DetLayer on which relies the segment
  DetId id = seg[last]->geographicalId();

  // Transform it in a TrajectoryStateOnSurface

  PTrajectoryStateOnDet seedTSOS = trajectoryStateTransform::persistentState(tsos, id.rawId());

  edm::OwnVector<TrackingRecHit> container;
  for (unsigned l = 0; l < seg.size(); l++) {
    container.push_back(seg[l]->hit()->clone());
    //container.push_back(seg[l]->hit());
  }

  TrajectorySeed theSeed(seedTSOS, container, alongMomentum);
  return theSeed;
}

/*
 * estimatePtCSC
 *
 * Look at delta phi to determine pt as:
 * pt = (c_1 * eta + c_2) / dphi
 *
 * Note that only segment pairs with one segment in the ME1 layers give sensitive results
 *
 */
void MuonSeedCreator::estimatePtCSC(const SegmentContainer& seg,
                                    const std::vector<int>& layers,
                                    double& thePt,
                                    double& theSpt) {
  unsigned size = seg.size();
  if (size < 2)
    return;

  // reverse the segment and layer container first for pure CSC case
  //if ( layers[0] > layers[ layers.size()-1 ] ) {
  //   reverse( layers.begin(), layers.end() );
  //   reverse( seg.begin(), seg.end() );
  //}

  std::vector<double> ptEstimate;
  std::vector<double> sptEstimate;

  thePt = defaultMomentum;
  theSpt = defaultMomentum;

  double pt = 0.;
  double spt = 0.;
  GlobalPoint segPos[2];

  int layer0 = layers[0];
  segPos[0] = seg[0]->globalPosition();
  float eta = fabs(segPos[0].eta());
  //float corr = fabs( tan(segPos[0].theta()) );
  // use pt from vertex information
  /*
  if ( layer0 == 0 ) {
     SegmentContainer seg0;
     seg0.push_back(seg[0]);
     std::vector<int> lyr0(1,0);
     estimatePtSingle( seg0, lyr0, thePt, theSpt);
     ptEstimate.push_back( thePt );
     sptEstimate.push_back( theSpt );
  }
  */

  //std::cout<<" estimate CSC "<<std::endl;

  unsigned idx1 = size;
  if (size > 1) {
    while (idx1 > 1) {
      idx1--;
      int layer1 = layers[idx1];
      if (layer0 == layer1)
        continue;
      segPos[1] = seg[idx1]->globalPosition();

      double dphi = segPos[0].phi() - segPos[1].phi();
      //double temp_dphi = dphi/corr;
      double temp_dphi = dphi;

      double sign = 1.0;
      if (temp_dphi < 0.) {
        temp_dphi = -1.0 * temp_dphi;
        sign = -1.0;
      }

      // Ensure that delta phi is not too small to prevent pt from blowing up
      if (temp_dphi < 0.0001) {
        temp_dphi = 0.0001;
        pt = theMaxMomentum;
        spt = theMaxMomentum * 0.25;
        ptEstimate.push_back(pt * sign);
        sptEstimate.push_back(spt);
      }
      // ME1 is inner-most
      if (layer0 == 0 && temp_dphi >= 0.0001) {
        // ME1/2 is outer-most
        if (layer1 == 1) {
          //temp_dphi = scaledPhi(temp_dphi, CSC01_1[3] );
          pt = getPt(CSC01, eta, temp_dphi)[0];
          spt = getPt(CSC01, eta, temp_dphi)[1];
        }
        // ME2 is outer-most
        else if (layer1 == 2) {
          //temp_dphi = scaledPhi(temp_dphi, CSC12_3[3] );
          pt = getPt(CSC02, eta, temp_dphi)[0];
          spt = getPt(CSC02, eta, temp_dphi)[1];
        }
        // ME3 is outer-most
        else if (layer1 == 3) {
          //temp_dphi = scaledPhi(temp_dphi, CSC13_3[3] );
          pt = getPt(CSC03, eta, temp_dphi)[0];
          spt = getPt(CSC03, eta, temp_dphi)[1];
        }
        // ME4 is outer-most
        else {
          //temp_dphi = scaledPhi(temp_dphi, CSC14_3[3]);
          pt = getPt(CSC14, eta, temp_dphi)[0];
          spt = getPt(CSC14, eta, temp_dphi)[1];
        }
        ptEstimate.push_back(pt * sign);
        sptEstimate.push_back(spt);
      }

      // ME1/2,ME1/3 is inner-most
      if (layer0 == 1) {
        // ME2 is outer-most
        if (layer1 == 2) {
          //if ( eta <= 1.2 )  {   temp_dphi = scaledPhi(temp_dphi, CSC12_1[3]); }
          //if ( eta >  1.2 )  {   temp_dphi = scaledPhi(temp_dphi, CSC12_2[3]); }
          pt = getPt(CSC12, eta, temp_dphi)[0];
          spt = getPt(CSC12, eta, temp_dphi)[1];
        }
        // ME3 is outer-most
        else if (layer1 == 3) {
          temp_dphi = scaledPhi(temp_dphi, CSC13_2[3]);
          pt = getPt(CSC13, eta, temp_dphi)[0];
          spt = getPt(CSC13, eta, temp_dphi)[1];
        }
        // ME4 is outer-most
        else {
          temp_dphi = scaledPhi(temp_dphi, CSC14_3[3]);
          pt = getPt(CSC14, eta, temp_dphi)[0];
          spt = getPt(CSC14, eta, temp_dphi)[1];
        }
        ptEstimate.push_back(pt * sign);
        sptEstimate.push_back(spt);
      }

      // ME2 is inner-most
      if (layer0 == 2 && temp_dphi > 0.0001) {
        // ME4 is outer-most
        if (layer1 == 4) {
          temp_dphi = scaledPhi(temp_dphi, CSC24_1[3]);
          pt = getPt(CSC24, eta, temp_dphi)[0];
          spt = getPt(CSC24, eta, temp_dphi)[1];
        }
        // ME3 is outer-most
        else {
          // if ME2-4 is availabe , discard ME2-3
          if (eta <= 1.7) {
            temp_dphi = scaledPhi(temp_dphi, CSC23_1[3]);
          }
          if (eta > 1.7) {
            temp_dphi = scaledPhi(temp_dphi, CSC23_2[3]);
          }
          pt = getPt(CSC23, eta, temp_dphi)[0];
          spt = getPt(CSC23, eta, temp_dphi)[1];
        }
        ptEstimate.push_back(pt * sign);
        sptEstimate.push_back(spt);
      }

      // ME3 is inner-most
      if (layer0 == 3 && temp_dphi > 0.0001) {
        temp_dphi = scaledPhi(temp_dphi, CSC34_1[3]);
        pt = getPt(CSC34, eta, temp_dphi)[0];
        spt = getPt(CSC34, eta, temp_dphi)[1];
        ptEstimate.push_back(pt * sign);
        sptEstimate.push_back(spt);
      }
    }
  }

  // Compute weighted average if have more than one estimator
  if (!ptEstimate.empty())
    weightedPt(ptEstimate, sptEstimate, thePt, theSpt);
}

/*
 * estimatePtDT
 *
 * Look at delta phi between segments to determine pt as:
 * pt = (c_1 * eta + c_2) / dphi
 */
void MuonSeedCreator::estimatePtDT(const SegmentContainer& seg,
                                   const std::vector<int>& layers,
                                   double& thePt,
                                   double& theSpt) {
  unsigned size = seg.size();
  if (size < 2)
    return;

  std::vector<double> ptEstimate;
  std::vector<double> sptEstimate;

  thePt = defaultMomentum;
  theSpt = defaultMomentum;

  double pt = 0.;
  double spt = 0.;
  GlobalPoint segPos[2];

  int layer0 = layers[0];
  segPos[0] = seg[0]->globalPosition();
  float eta = fabs(segPos[0].eta());

  //std::cout<<" estimate DT "<<std::endl;
  // inner-most layer
  //for ( unsigned idx0 = 0; idx0 < size-1; ++idx0 ) {
  //  layer0 = layers[idx0];
  //  segPos[0]  = seg[idx0]->globalPosition();
  // outer-most layer
  // for ( unsigned idx1 = idx0+1; idx1 < size; ++idx1 ) {
  for (unsigned idx1 = 1; idx1 < size; ++idx1) {
    int layer1 = layers[idx1];
    segPos[1] = seg[idx1]->globalPosition();

    //eta = fabs(segPos[1].eta());
    //if (layer1 == -4) eta = fabs(segPos[0].eta());

    double dphi = segPos[0].phi() - segPos[1].phi();
    double temp_dphi = dphi;

    // Ensure that delta phi is not too small to prevent pt from blowing up

    double sign = 1.0;
    if (temp_dphi < 0.) {
      temp_dphi = -temp_dphi;
      sign = -1.0;
    }

    if (temp_dphi < 0.0001) {
      temp_dphi = 0.0001;
      pt = theMaxMomentum;
      spt = theMaxMomentum * 0.25;
      ptEstimate.push_back(pt * sign);
      sptEstimate.push_back(spt);
    }

    // MB1 is inner-most
    if (layer0 == -1 && temp_dphi > 0.0001) {
      // MB2 is outer-most
      if (layer1 == -2) {
        if (eta <= 0.7) {
          temp_dphi = scaledPhi(temp_dphi, DT12_1[3]);
        }
        if (eta > 0.7) {
          temp_dphi = scaledPhi(temp_dphi, DT12_2[3]);
        }
        pt = getPt(DT12, eta, temp_dphi)[0];
        spt = getPt(DT12, eta, temp_dphi)[1];
      }
      // MB3 is outer-most
      else if (layer1 == -3) {
        if (eta <= 0.6) {
          temp_dphi = scaledPhi(temp_dphi, DT13_1[3]);
        }
        if (eta > 0.6) {
          temp_dphi = scaledPhi(temp_dphi, DT13_2[3]);
        }
        pt = getPt(DT13, eta, temp_dphi)[0];
        spt = getPt(DT13, eta, temp_dphi)[1];
      }
      // MB4 is outer-most
      else {
        if (eta <= 0.52) {
          temp_dphi = scaledPhi(temp_dphi, DT14_1[3]);
        }
        if (eta > 0.52) {
          temp_dphi = scaledPhi(temp_dphi, DT14_2[3]);
        }
        pt = getPt(DT14, eta, temp_dphi)[0];
        spt = getPt(DT14, eta, temp_dphi)[1];
      }
      ptEstimate.push_back(pt * sign);
      sptEstimate.push_back(spt);
    }

    // MB2 is inner-most
    if (layer0 == -2 && temp_dphi > 0.0001) {
      // MB3 is outer-most
      if (layer1 == -3) {
        if (eta <= 0.6) {
          temp_dphi = scaledPhi(temp_dphi, DT23_1[3]);
        }
        if (eta > 0.6) {
          temp_dphi = scaledPhi(temp_dphi, DT23_2[3]);
        }
        pt = getPt(DT23, eta, temp_dphi)[0];
        spt = getPt(DT23, eta, temp_dphi)[1];
      }
      // MB4 is outer-most
      else {
        if (eta <= 0.52) {
          temp_dphi = scaledPhi(temp_dphi, DT24_1[3]);
        }
        if (eta > 0.52) {
          temp_dphi = scaledPhi(temp_dphi, DT24_2[3]);
        }
        pt = getPt(DT24, eta, temp_dphi)[0];
        spt = getPt(DT24, eta, temp_dphi)[1];
      }
      ptEstimate.push_back(pt * sign);
      sptEstimate.push_back(spt);
    }

    // MB3 is inner-most    -> only marginally useful to pick up the charge
    if (layer0 == -3 && temp_dphi > 0.0001) {
      // MB4 is outer-most

      if (eta <= 0.51) {
        temp_dphi = scaledPhi(temp_dphi, DT34_1[3]);
      }
      if (eta > 0.51) {
        temp_dphi = scaledPhi(temp_dphi, DT34_2[3]);
      }
      pt = getPt(DT34, eta, temp_dphi)[0];
      spt = getPt(DT34, eta, temp_dphi)[1];
      ptEstimate.push_back(pt * sign);
      sptEstimate.push_back(spt);
    }
  }
  //}

  // Compute weighted average if have more than one estimator
  if (!ptEstimate.empty())
    weightedPt(ptEstimate, sptEstimate, thePt, theSpt);
}

/*
 * estimatePtOverlap
 *
 */
void MuonSeedCreator::estimatePtOverlap(const SegmentContainer& seg,
                                        const std::vector<int>& layers,
                                        double& thePt,
                                        double& theSpt) {
  int size = layers.size();

  thePt = defaultMomentum;
  theSpt = defaultMomentum;

  SegmentContainer segCSC;
  std::vector<int> layersCSC;
  SegmentContainer segDT;
  std::vector<int> layersDT;

  // DT layers are numbered as -4 to -1, whereas CSC layers go from 0 to 4:
  for (unsigned j = 0; j < layers.size(); ++j) {
    if (layers[j] > -1) {
      segCSC.push_back(seg[j]);
      layersCSC.push_back(layers[j]);
    } else {
      segDT.push_back(seg[j]);
      layersDT.push_back(layers[j]);
    }
  }

  std::vector<double> ptEstimate;
  std::vector<double> sptEstimate;

  GlobalPoint segPos[2];
  int layer0 = layers[0];
  segPos[0] = seg[0]->globalPosition();
  float eta = fabs(segPos[0].eta());
  //std::cout<<" estimate OL "<<std::endl;

  if (!segDT.empty() && !segCSC.empty()) {
    int layer1 = layers[size - 1];
    segPos[1] = seg[size - 1]->globalPosition();

    double dphi = segPos[0].phi() - segPos[1].phi();
    double temp_dphi = dphi;

    // Ensure that delta phi is not too small to prevent pt from blowing up

    double sign = 1.0;
    if (temp_dphi < 0.) {
      temp_dphi = -temp_dphi;
      sign = -1.0;
    }

    if (temp_dphi < 0.0001) {
      temp_dphi = 0.0001;
      thePt = theMaxMomentum;
      theSpt = theMaxMomentum * 0.25;
      ptEstimate.push_back(thePt * sign);
      sptEstimate.push_back(theSpt);
    }

    // MB1 is inner-most
    if (layer0 == -1 && temp_dphi > 0.0001) {
      // ME1/3 is outer-most
      if (layer1 == 1) {
        temp_dphi = scaledPhi(temp_dphi, OL_1213[3]);
        thePt = getPt(OL1213, eta, temp_dphi)[0];
        theSpt = getPt(OL1213, eta, temp_dphi)[1];
      }
      // ME2 is outer-most
      else if (layer1 == 2) {
        temp_dphi = scaledPhi(temp_dphi, OL_1222[3]);
        thePt = getPt(OL1222, eta, temp_dphi)[0];
        theSpt = getPt(OL1222, eta, temp_dphi)[1];
      }
      // ME3 is outer-most
      else {
        temp_dphi = scaledPhi(temp_dphi, OL_1232[3]);
        thePt = getPt(OL1232, eta, temp_dphi)[0];
        theSpt = getPt(OL1232, eta, temp_dphi)[1];
      }
      ptEstimate.push_back(thePt * sign);
      sptEstimate.push_back(theSpt);
    }
    // MB2 is inner-most
    if (layer0 == -2 && temp_dphi > 0.0001) {
      // ME1/3 is outer-most
      if (layer1 == 1) {
        temp_dphi = scaledPhi(temp_dphi, OL_2213[3]);
        thePt = getPt(OL2213, eta, temp_dphi)[0];
        theSpt = getPt(OL2213, eta, temp_dphi)[1];
        ptEstimate.push_back(thePt * sign);
        sptEstimate.push_back(theSpt);
      }
      // ME2 is outer-most
      if (layer1 == 2) {
        temp_dphi = scaledPhi(temp_dphi, OL_2222[3]);
        thePt = getPt(OL2222, eta, temp_dphi)[0];
        theSpt = getPt(OL2222, eta, temp_dphi)[1];
      }
    }
  }

  if (segDT.size() > 1) {
    estimatePtDT(segDT, layersDT, thePt, theSpt);
    ptEstimate.push_back(thePt);
    sptEstimate.push_back(theSpt);
  }

  /*
  // not useful ....and pt estimation is bad
  if ( segCSC.size() > 1 ) {
    // don't estimate pt without ME1 information
    bool CSCLayer1=false;
    for (unsigned i=0; i< layersCSC.size(); i++) {
        if ( layersCSC[i]==0 || layersCSC[i]==1 ) CSCLayer1 = true;
    }
    if (CSCLayer1) {
       estimatePtCSC(segCSC, layersCSC, thePt, theSpt);
       ptEstimate.push_back(thePt);
       sptEstimate.push_back(theSpt);
    }
  }
  */

  // Compute weighted average if have more than one estimator
  if (!ptEstimate.empty())
    weightedPt(ptEstimate, sptEstimate, thePt, theSpt);
}
/*
 *
 *   estimate Pt for single segment events
 *
 */
void MuonSeedCreator::estimatePtSingle(const SegmentContainer& seg,
                                       const std::vector<int>& layers,
                                       double& thePt,
                                       double& theSpt) {
  thePt = defaultMomentum;
  theSpt = defaultMomentum;

  GlobalPoint segPos = seg[0]->globalPosition();
  double eta = segPos.eta();
  GlobalVector gv = seg[0]->globalDirection();

  // 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() * segPos.x() + gv.y() * segPos.y());
  cosDpsi /= sqrt(segPos.x() * segPos.x() + segPos.y() * segPos.y());
  cosDpsi /= sqrt(gv.x() * gv.x() + gv.y() * gv.y());

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

  double dpsi = fabs(acos(cosDpsi));
  if (dpsi > 1.570796) {
    dpsi = 3.141592 - dpsi;
    sign = -1. * sign;
  }
  if (fabs(dpsi) < 0.00005) {
    dpsi = 0.00005;
  }

  // the 1st layer
  if (layers[0] == -1) {
    // MB10
    if (fabs(eta) < 0.3) {
      dpsi = scaledPhi(dpsi, SMB_10S[3]);
      thePt = getPt(SMB10, eta, dpsi)[0];
      theSpt = getPt(SMB10, eta, dpsi)[1];
    }
    // MB11
    if (fabs(eta) >= 0.3 && fabs(eta) < 0.82) {
      dpsi = scaledPhi(dpsi, SMB_11S[3]);
      thePt = getPt(SMB11, eta, dpsi)[0];
      theSpt = getPt(SMB11, eta, dpsi)[1];
    }
    // MB12
    if (fabs(eta) >= 0.82 && fabs(eta) < 1.2) {
      dpsi = scaledPhi(dpsi, SMB_12S[3]);
      thePt = getPt(SMB12, eta, dpsi)[0];
      theSpt = getPt(SMB12, eta, dpsi)[1];
    }
  }
  if (layers[0] == 1) {
    // ME13
    if (fabs(eta) > 0.92 && fabs(eta) < 1.16) {
      dpsi = scaledPhi(dpsi, SME_13S[3]);
      thePt = getPt(SME13, eta, dpsi)[0];
      theSpt = getPt(SME13, eta, dpsi)[1];
    }
    // ME12
    if (fabs(eta) >= 1.16 && fabs(eta) <= 1.6) {
      dpsi = scaledPhi(dpsi, SME_12S[3]);
      thePt = getPt(SME12, eta, dpsi)[0];
      theSpt = getPt(SME12, eta, dpsi)[1];
    }
  }
  if (layers[0] == 0) {
    // ME11
    if (fabs(eta) > 1.6) {
      dpsi = scaledPhi(dpsi, SMB_11S[3]);
      thePt = getPt(SME11, eta, dpsi)[0];
      theSpt = getPt(SME11, eta, dpsi)[1];
    }
  }
  // the 2nd layer
  if (layers[0] == -2) {
    // MB20
    if (fabs(eta) < 0.25) {
      dpsi = scaledPhi(dpsi, SMB_20S[3]);
      thePt = getPt(SMB20, eta, dpsi)[0];
      theSpt = getPt(SMB20, eta, dpsi)[1];
    }
    // MB21
    if (fabs(eta) >= 0.25 && fabs(eta) < 0.72) {
      dpsi = scaledPhi(dpsi, SMB_21S[3]);
      thePt = getPt(SMB21, eta, dpsi)[0];
      theSpt = getPt(SMB21, eta, dpsi)[1];
    }
    // MB22
    if (fabs(eta) >= 0.72 && fabs(eta) < 1.04) {
      dpsi = scaledPhi(dpsi, SMB_22S[3]);
      thePt = getPt(SMB22, eta, dpsi)[0];
      theSpt = getPt(SMB22, eta, dpsi)[1];
    }
  }
  if (layers[0] == 2) {
    // ME22
    if (fabs(eta) > 0.95 && fabs(eta) <= 1.6) {
      dpsi = scaledPhi(dpsi, SME_22S[3]);
      thePt = getPt(SME22, eta, dpsi)[0];
      theSpt = getPt(SME22, eta, dpsi)[1];
    }
    // ME21
    if (fabs(eta) > 1.6 && fabs(eta) < 2.45) {
      dpsi = scaledPhi(dpsi, SME_21S[3]);
      thePt = getPt(SME21, eta, dpsi)[0];
      theSpt = getPt(SME21, eta, dpsi)[1];
    }
  }

  // the 3rd layer
  if (layers[0] == -3) {
    // MB30
    if (fabs(eta) <= 0.22) {
      dpsi = scaledPhi(dpsi, SMB_30S[3]);
      thePt = getPt(SMB30, eta, dpsi)[0];
      theSpt = getPt(SMB30, eta, dpsi)[1];
    }
    // MB31
    if (fabs(eta) > 0.22 && fabs(eta) <= 0.6) {
      dpsi = scaledPhi(dpsi, SMB_31S[3]);
      thePt = getPt(SMB31, eta, dpsi)[0];
      theSpt = getPt(SMB31, eta, dpsi)[1];
    }
    // MB32
    if (fabs(eta) > 0.6 && fabs(eta) < 0.95) {
      dpsi = scaledPhi(dpsi, SMB_32S[3]);
      thePt = getPt(SMB32, eta, dpsi)[0];
      theSpt = getPt(SMB32, eta, dpsi)[1];
    }
  }
  thePt = fabs(thePt) * sign;
  theSpt = fabs(theSpt);

  return;
}

// setup the minimum value for obvious showering cases
void MuonSeedCreator::estimatePtShowering(int& NShowers, int& NShowerSegments, double& thePt, double& theSpt) {
  if (NShowers > 2 && thePt < 300.) {
    thePt = 800.;
    theSpt = 200.;
  }
  if (NShowers == 2 && NShowerSegments > 11 && thePt < 150.) {
    thePt = 280.;
    theSpt = 70.;
  }
  if (NShowers == 2 && NShowerSegments <= 11 && thePt < 50.) {
    thePt = 80.;
    theSpt = 40.;
  }
  if (NShowers == 1 && NShowerSegments <= 5 && thePt < 10.) {
    thePt = 16.;
    theSpt = 8.;
  }
}

/*
 * weightedPt
 *
 * Look at delta phi between segments to determine pt as:
 * pt = (c_1 * eta + c_2) / dphi
 */
void MuonSeedCreator::weightedPt(const std::vector<double>& ptEstimate,
                                 const std::vector<double>& sptEstimate,
                                 double& thePt,
                                 double& theSpt) {
  int size = ptEstimate.size();

  // If only one element, by-pass computation below
  if (size < 2) {
    thePt = ptEstimate[0];
    theSpt = sptEstimate[0];
    return;
  }

  double charge = 0.;
  // If have more than one pt estimator, first look if any estimator is biased
  // by comparing the charge of each estimator

  for (unsigned j = 0; j < ptEstimate.size(); j++) {
    //std::cout<<" weighting pt: "<< ptEstimate[j] <<std::endl;
    if (ptEstimate[j] < 0.) {
      // To prevent from blowing up, add 0.1
      // weight by relative error on pt
      charge -= 1. * (ptEstimate[j] * ptEstimate[j]) / (sptEstimate[j] * sptEstimate[j]);
    } else {
      // weight by relative error on pt
      charge += 1. * (ptEstimate[j] * ptEstimate[j]) / (sptEstimate[j] * sptEstimate[j]);
    }
  }

  // No need to normalize as we want to know only sign ( + or - )
  if (charge < 0.) {
    charge = -1.;
  } else {
    charge = 1.;
  }

  //int n = 0;
  double weightPtSum = 0.;
  double sigmaSqr_sum = 0.;

  // Now, we want to compute average Pt using estimators with "correct" charge
  // This is to remove biases
  for (unsigned j = 0; j < ptEstimate.size(); ++j) {
    //if ( (minpt_ratio < 0.5) && (fabs(ptEstimate[j]) < 5.0) ) continue;
    //if ( ptEstimate[j] * charge > 0. ) {
    //n++;
    sigmaSqr_sum += 1.0 / (sptEstimate[j] * sptEstimate[j]);
    weightPtSum += fabs(ptEstimate[j]) / (sptEstimate[j] * sptEstimate[j]);
    //}
  }
  /*  
  if (n < 1) {
    thePt  = defaultMomentum*charge;
    theSpt = defaultMomentum; 
    return;
  } 
  */
  // Compute weighted mean and error

  thePt = (charge * weightPtSum) / sigmaSqr_sum;
  theSpt = sqrt(1.0 / sigmaSqr_sum);

  //std::cout<<" final weighting : "<< thePt <<" ~ "<< fabs( theSpt/thePt ) <<std::endl;

  return;
}

std::vector<double> MuonSeedCreator::getPt(const std::vector<double>& vPara, double eta, double dPhi) {
  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) * estPt;
  std::vector<double> paraPt;
  paraPt.push_back(estPt);
  paraPt.push_back(estSPt);

  //std::cout<<"      pt:"<<estPt<<" +/-"<< estSPt<<"   h:"<<eta<<"  df:"<<dPhi<<std::endl;
  return paraPt;
}

double MuonSeedCreator::scaledPhi(double dphi, double t1) {
  if (dphi != 0.) {
    double oPhi = 1. / dphi;
    dphi = dphi / (1. + t1 / (oPhi + 10.));
    return dphi;

  } else {
    return dphi;
  }
}