#include <MuonSeedCreator.h>

Public Types
typedef MuonTransientTrackingRecHit::MuonRecHitContainer	SegmentContainer
Public Member Functions
TrajectorySeed	createSeed (int type, SegmentContainer seg, std::vector< int > layers, int NShower, int NShowerSeg)
	Create a seed from set of segments.
	MuonSeedCreator (const edm::ParameterSet &pset)
	Constructor.
void	setBField (const MagneticField *theField)
	Cache Magnetic Field for current event.
	~MuonSeedCreator ()
	Destructor.
Private Member Functions
void	estimatePtCSC (SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
	Estimate transverse momentum of track from CSC measurements.
void	estimatePtDT (SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
	Estimate transverse momentum of track from DT measurements.
void	estimatePtOverlap (SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
	Estimate transverse momentum of track from CSC + DT measurements.
void	estimatePtShowering (int &NShowers, int &NShowerSeg, double &pt, double &spt)
	Estimate transverse momentum of track from showering segment.
void	estimatePtSingle (SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
	Estimate transverse momentum of track from single segment.
std::vector< double >	getPt (std::vector< double > vParameters, double eta, double dPhi)
	Compute pt from parameters.
double	scaledPhi (double dphi, double t1)
	Scale the dPhi from segment position.
void	weightedPt (std::vector< double > ptEstimate, std::vector< double > sptEstimate, double &ptAvg, double &sptAvg)
	Compute weighted mean pt from different pt estimators.
Private Attributes
const MagneticField *	BField
std::vector< double >	CSC01
std::vector< double >	CSC01_1
std::vector< double >	CSC02
std::vector< double >	CSC03
std::vector< double >	CSC12
std::vector< double >	CSC12_1
std::vector< double >	CSC12_2
std::vector< double >	CSC12_3
std::vector< double >	CSC13
std::vector< double >	CSC13_2
std::vector< double >	CSC13_3
std::vector< double >	CSC14
std::vector< double >	CSC14_3
std::vector< double >	CSC23
std::vector< double >	CSC23_1
std::vector< double >	CSC23_2
std::vector< double >	CSC24
std::vector< double >	CSC24_1
std::vector< double >	CSC34
std::vector< double >	CSC34_1
bool	debug
float	defaultMomentum
std::vector< double >	DT12
std::vector< double >	DT12_1
std::vector< double >	DT12_2
std::vector< double >	DT13
std::vector< double >	DT13_1
std::vector< double >	DT13_2
std::vector< double >	DT14
std::vector< double >	DT14_1
std::vector< double >	DT14_2
std::vector< double >	DT23
std::vector< double >	DT23_1
std::vector< double >	DT23_2
std::vector< double >	DT24
std::vector< double >	DT24_1
std::vector< double >	DT24_2
std::vector< double >	DT34
std::vector< double >	DT34_1
std::vector< double >	DT34_2
std::vector< double >	OL1213
std::vector< double >	OL1222
std::vector< double >	OL1232
std::vector< double >	OL2213
std::vector< double >	OL2222
std::vector< double >	OL_1213
std::vector< double >	OL_1222
std::vector< double >	OL_1232
std::vector< double >	OL_2213
std::vector< double >	OL_2222
std::vector< double >	SMB10
std::vector< double >	SMB11
std::vector< double >	SMB12
std::vector< double >	SMB20
std::vector< double >	SMB21
std::vector< double >	SMB22
std::vector< double >	SMB30
std::vector< double >	SMB31
std::vector< double >	SMB32
std::vector< double >	SMB_10S
std::vector< double >	SMB_11S
std::vector< double >	SMB_12S
std::vector< double >	SMB_20S
std::vector< double >	SMB_21S
std::vector< double >	SMB_22S
std::vector< double >	SMB_30S
std::vector< double >	SMB_31S
std::vector< double >	SMB_32S
std::vector< double >	SME11
std::vector< double >	SME12
std::vector< double >	SME13
std::vector< double >	SME21
std::vector< double >	SME22
std::vector< double >	SME31
std::vector< double >	SME32
std::vector< double >	SME41
std::vector< double >	SME_11S
std::vector< double >	SME_12S
std::vector< double >	SME_13S
std::vector< double >	SME_21S
std::vector< double >	SME_22S
double	sysError
float	theMaxMomentum
float	theMinMomentum

Detailed Description

Creates seed from vector of segment

Determine pt of seed using various combination of segments from different layers (stations) Parameterization used to determine pt between layers i and j:

pt = [ c_0 + c_1 * (Delta phi_ij) + c_2 * (Delta phi_ij)^2 ] / eta

Author:: Dominique Fortin - UCR

Definition at line 30 of file MuonSeedCreator.h.

Member Typedef Documentation

typedef MuonTransientTrackingRecHit::MuonRecHitContainer MuonSeedCreator::SegmentContainer

Definition at line 34 of file MuonSeedCreator.h.

Constructor & Destructor Documentation

MuonSeedCreator::MuonSeedCreator ( const edm::ParameterSet & pset ) [explicit]

Constructor.

See header file for a description of this class.

Author:: : Shih-Chuan Kao, Dominique Fortin - UCR

Definition at line 37 of file MuonSeedCreator.cc.

References CSC01_1, CSC02, CSC03, CSC12, CSC12_1, CSC12_2, CSC12_3, CSC13, CSC13_2, CSC13_3, CSC14, CSC14_3, CSC23, CSC23_1, CSC23_2, CSC24, CSC24_1, CSC34, CSC34_1, debug, defaultMomentum, DT12, DT12_1, DT12_2, DT13, DT13_1, DT13_2, DT14, DT14_1, DT14_2, DT23, DT23_1, DT23_2, DT24, DT24_1, DT24_2, DT34, DT34_1, DT34_2, edm::ParameterSet::getParameter(), OL1213, OL1222, OL1232, OL2222, OL_1213, OL_1222, OL_1232, OL_2213, OL_2222, SMB10, SMB11, SMB12, SMB20, SMB21, SMB22, SMB30, SMB31, SMB32, SMB_10S, SMB_11S, SMB_12S, SMB_20S, SMB_21S, SMB_22S, SMB_30S, SMB_31S, SMB_32S, SME11, SME12, SME13, SME21, SME22, SME31, SME32, SME41, SME_11S, SME_12S, SME_13S, SME_21S, SME_22S, sysError, theMaxMomentum, and theMinMomentum.

                                                           {
  
  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");

}

MuonSeedCreator::~MuonSeedCreator ( )

Destructor.

Definition at line 141 of file MuonSeedCreator.cc.

                                 {
  
}

Member Function Documentation

TrajectorySeed MuonSeedCreator::createSeed	(	int	type,
		SegmentContainer	seg,
		std::vector< int >	layers,
		int	NShower,
		int	NShowerSeg
	)

Create a seed from set of segments.

get the Global position

get the Global direction

scale the magnitude of total momentum

Trasfer into local direction

get the Global position

get the Global direction

count the energy loss - from parameterization

scale the magnitude of total momentum

Trasfer into local direction

Definition at line 153 of file MuonSeedCreator.cc.

References alongMomentum, BField, DeDxDiscriminatorTools::charge(), clone(), gather_cfg::cout, debug, error, estimatePtCSC(), estimatePtDT(), estimatePtOverlap(), estimatePtShowering(), estimatePtSingle(), eta(), PV3DBase< T, PVType, FrameType >::eta(), i, prof2calltree::l, prof2calltree::last, PV3DBase< T, PVType, FrameType >::perp(), trajectoryStateTransform::persistentState(), edm::OwnVector< T, P >::push_back(), lumiQTWidget::t, funct::tan(), theMaxMomentum, theMinMomentum, and PV3DBase< T, PVType, FrameType >::theta().

Referenced by MuonSeedBuilder::build().

                                                                                                                                    {

  // 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;
}

void MuonSeedCreator::estimatePtCSC	(	SegmentContainer	seg,
		std::vector< int >	layers,
		double &	pt,
		double &	spt
	)		`[private]`

Estimate transverse momentum of track from CSC measurements.

Definition at line 332 of file MuonSeedCreator.cc.

References CSC02, CSC03, CSC12, CSC13, CSC13_2, CSC14, CSC14_3, CSC23, CSC23_1, CSC23_2, CSC24, CSC24_1, CSC34, CSC34_1, defaultMomentum, eta(), getPt(), PV3DBase< T, PVType, FrameType >::phi(), scaledPhi(), findQualityFiles::size, theMaxMomentum, and weightedPt().

Referenced by createSeed().

                                                                                                               {

  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
        bool ME4av =false;
        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];
          ME4av = true;
        }
        // ME3 is outer-most
        else {
          // if ME2-4 is availabe , discard ME2-3 
          if ( !ME4av ) {
            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.size() > 0 ) weightedPt( ptEstimate, sptEstimate, thePt, theSpt);

}

void MuonSeedCreator::estimatePtDT	(	SegmentContainer	seg,
		std::vector< int >	layers,
		double &	pt,
		double &	spt
	)		`[private]`

Estimate transverse momentum of track from DT measurements.

Definition at line 504 of file MuonSeedCreator.cc.

References defaultMomentum, DT12, DT12_1, DT12_2, DT13, DT13_1, DT13_2, DT14, DT14_1, DT14_2, DT23, DT23_1, DT23_2, DT24, DT24_1, DT24_2, DT34, DT34_1, DT34_2, eta(), getPt(), PV3DBase< T, PVType, FrameType >::phi(), scaledPhi(), findQualityFiles::size, theMaxMomentum, and weightedPt().

Referenced by createSeed(), and estimatePtOverlap().

                                                                                                             {

  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
      bool MB23av = false;
      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];
          MB23av = true;
        }
        // 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];
          MB23av = true;
        }
        // MB4 is outer-most
        else {
          if ( !MB23av ) {
             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.size() > 0 ) weightedPt( ptEstimate, sptEstimate, thePt, theSpt);

}

void MuonSeedCreator::estimatePtOverlap	(	SegmentContainer	seg,
		std::vector< int >	layers,
		double &	pt,
		double &	spt
	)		`[private]`

Estimate transverse momentum of track from CSC + DT measurements.

Definition at line 638 of file MuonSeedCreator.cc.

References defaultMomentum, estimatePtDT(), eta(), getPt(), j, OL1213, OL1222, OL1232, OL2213, OL2222, OL_1213, OL_1222, OL_1232, OL_2213, OL_2222, PV3DBase< T, PVType, FrameType >::phi(), scaledPhi(), findQualityFiles::size, theMaxMomentum, and weightedPt().

Referenced by createSeed().

                                                                                                                  {

  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.size() > 0 && segCSC.size() > 0 ) {
    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.size() > 0 ) weightedPt( ptEstimate, sptEstimate, thePt, theSpt);

}

void MuonSeedCreator::estimatePtShowering	(	int &	NShowers,
		int &	NShowerSeg,
		double &	pt,
		double &	spt
	)		`[private]`

Estimate transverse momentum of track from showering segment.

Definition at line 901 of file MuonSeedCreator.cc.

Referenced by createSeed().

                                                                                                             {

  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. ; 
  }

}

void MuonSeedCreator::estimatePtSingle	(	SegmentContainer	seg,
		std::vector< int >	layers,
		double &	pt,
		double &	spt
	)		`[private]`

Estimate transverse momentum of track from single segment.

Definition at line 767 of file MuonSeedCreator.cc.

References defaultMomentum, eta(), PV3DBase< T, PVType, FrameType >::eta(), getPt(), scaledPhi(), SMB10, SMB11, SMB12, SMB20, SMB21, SMB22, SMB30, SMB31, SMB32, SMB_10S, SMB_11S, SMB_12S, SMB_20S, SMB_21S, SMB_22S, SMB_30S, SMB_31S, SMB_32S, SME11, SME12, SME13, SME21, SME22, SME_12S, SME_13S, SME_21S, SME_22S, mathSSE::sqrt(), PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::y().

Referenced by createSeed().

                                                                                                                 {

  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;
}

std::vector< double > MuonSeedCreator::getPt	(	std::vector< double >	vParameters,
		double	eta,
		double	dPhi
	)		`[private]`

Compute pt from parameters.

Definition at line 992 of file MuonSeedCreator.cc.

References h.

Referenced by estimatePtCSC(), estimatePtDT(), estimatePtOverlap(), and estimatePtSingle().

                                                                                          {

       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
	)		`[private]`

Scale the dPhi from segment position.

Definition at line 1005 of file MuonSeedCreator.cc.

Referenced by estimatePtCSC(), estimatePtDT(), estimatePtOverlap(), and estimatePtSingle().

                                                         {

  if (dphi != 0. ) {

    double oPhi = 1./dphi ;
    dphi = dphi /( 1. + t1/( oPhi + 10. ) ) ;
    return dphi ;

  } else {
    return dphi ;
  } 

}

void MuonSeedCreator::setBField ( const MagneticField * theField ) [inline]

Cache Magnetic Field for current event.

Definition at line 45 of file MuonSeedCreator.h.

References BField.

Referenced by MuonSeedBuilder::build().

{ BField = theField; };

void MuonSeedCreator::weightedPt	(	std::vector< double >	ptEstimate,
		std::vector< double >	sptEstimate,
		double &	ptAvg,
		double &	sptAvg
	)		`[private]`

Compute weighted mean pt from different pt estimators.

Definition at line 928 of file MuonSeedCreator.cc.

References DeDxDiscriminatorTools::charge(), j, findQualityFiles::size, and mathSSE::sqrt().

Referenced by estimatePtCSC(), estimatePtDT(), and estimatePtOverlap().

                                                                                                                           {

 
  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 
      charge -= 1. * (ptEstimate[j]*ptEstimate[j]) / (sptEstimate[j]*sptEstimate[j] );  // weight by relative error on pt
    } else {
      charge += 1. * (ptEstimate[j]*ptEstimate[j]) / (sptEstimate[j]*sptEstimate[j] );  // weight by relative error on pt
    }
  }
 
  // 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;
}