MuonSeedCreator Class Reference

#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. More...
	MuonSeedCreator (const edm::ParameterSet &pset)
	Constructor. More...
void	setBField (const MagneticField *theField)
	Cache Magnetic Field for current event. More...
	~MuonSeedCreator ()
	Destructor. More...

Private Member Functions
void	estimatePtCSC (SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
	Estimate transverse momentum of track from CSC measurements. More...
void	estimatePtDT (SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
	Estimate transverse momentum of track from DT measurements. More...
void	estimatePtOverlap (SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
	Estimate transverse momentum of track from CSC + DT measurements. More...
void	estimatePtShowering (int &NShowers, int &NShowerSeg, double &pt, double &spt)
	Estimate transverse momentum of track from showering segment. More...
void	estimatePtSingle (SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
	Estimate transverse momentum of track from single segment. More...
std::vector< double >	getPt (std::vector< double > vParameters, double eta, double dPhi)
	Compute pt from parameters. More...
double	scaledPhi (double dphi, double t1)
	Scale the dPhi from segment position. More...
void	weightedPt (std::vector< double > ptEstimate, std::vector< double > sptEstimate, double &ptAvg, double &sptAvg)
	Compute weighted mean pt from different pt estimators. More...

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(), PV3DBase< T, PVType, FrameType >::eta(), eta(), i, prof2calltree::l, prof2calltree::last, PV3DBase< T, PVType, FrameType >::perp(), TrajectoryStateTransform::persistentState(), edm::OwnVector< T, P >::push_back(), matplotRender::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(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
  TrajectoryStateTransform tsTransform;
  
  //PTrajectoryStateOnDet *seedTSOS = tsTransform.persistentState( tsos, id.rawId());
  std::auto_ptr<PTrajectoryStateOnDet> seedTSOS(tsTransform.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 333 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 505 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 639 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 902 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 768 of file MuonSeedCreator.cc.

References defaultMomentum, PV3DBase< T, PVType, FrameType >::eta(), 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 993 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 1006 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 929 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;
}

Member Data Documentation

const MagneticField* MuonSeedCreator::BField

private

Definition at line 90 of file MuonSeedCreator.h.

Referenced by createSeed(), and setBField().

std::vector<double> MuonSeedCreator::CSC01

private

Definition at line 100 of file MuonSeedCreator.h.

std::vector<double> MuonSeedCreator::CSC01_1

private

Definition at line 137 of file MuonSeedCreator.h.

Referenced by MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC02

private

Definition at line 102 of file MuonSeedCreator.h.

Referenced by estimatePtCSC(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC03

private

Definition at line 104 of file MuonSeedCreator.h.

Referenced by estimatePtCSC(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC12

private

Definition at line 101 of file MuonSeedCreator.h.

Referenced by estimatePtCSC(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC12_1

private

Definition at line 138 of file MuonSeedCreator.h.

Referenced by MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC12_2

private

Definition at line 139 of file MuonSeedCreator.h.

Referenced by MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC12_3

private

Definition at line 140 of file MuonSeedCreator.h.

Referenced by MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC13

private

Definition at line 103 of file MuonSeedCreator.h.

Referenced by estimatePtCSC(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC13_2

private

Definition at line 141 of file MuonSeedCreator.h.

Referenced by estimatePtCSC(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC13_3

private

Definition at line 142 of file MuonSeedCreator.h.

Referenced by MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC14

private

Definition at line 105 of file MuonSeedCreator.h.

Referenced by estimatePtCSC(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC14_3

private

Definition at line 143 of file MuonSeedCreator.h.

Referenced by estimatePtCSC(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC23

private

Definition at line 106 of file MuonSeedCreator.h.

Referenced by estimatePtCSC(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC23_1

private

Definition at line 144 of file MuonSeedCreator.h.

Referenced by estimatePtCSC(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC23_2

private

Definition at line 145 of file MuonSeedCreator.h.

Referenced by estimatePtCSC(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC24

private

Definition at line 107 of file MuonSeedCreator.h.

Referenced by estimatePtCSC(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC24_1

private

Definition at line 146 of file MuonSeedCreator.h.

Referenced by estimatePtCSC(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC34

private

Definition at line 108 of file MuonSeedCreator.h.

Referenced by estimatePtCSC(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::CSC34_1

private

Definition at line 147 of file MuonSeedCreator.h.

Referenced by estimatePtCSC(), and MuonSeedCreator().

bool MuonSeedCreator::debug

private

Definition at line 87 of file MuonSeedCreator.h.

Referenced by createSeed(), and MuonSeedCreator().

float MuonSeedCreator::defaultMomentum

private

Definition at line 81 of file MuonSeedCreator.h.

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

std::vector<double> MuonSeedCreator::DT12

private

Definition at line 93 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT12_1

private

Definition at line 149 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT12_2

private

Definition at line 150 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT13

private

Definition at line 94 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT13_1

private

Definition at line 151 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT13_2

private

Definition at line 152 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT14

private

Definition at line 95 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT14_1

private

Definition at line 153 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT14_2

private

Definition at line 154 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT23

private

Definition at line 96 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT23_1

private

Definition at line 155 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT23_2

private

Definition at line 156 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT24

private

Definition at line 97 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT24_1

private

Definition at line 157 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT24_2

private

Definition at line 158 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT34

private

Definition at line 98 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT34_1

private

Definition at line 159 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::DT34_2

private

Definition at line 160 of file MuonSeedCreator.h.

Referenced by estimatePtDT(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::OL1213

private

Definition at line 110 of file MuonSeedCreator.h.

Referenced by estimatePtOverlap(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::OL1222

private

Definition at line 111 of file MuonSeedCreator.h.

Referenced by estimatePtOverlap(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::OL1232

private

Definition at line 112 of file MuonSeedCreator.h.

Referenced by estimatePtOverlap(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::OL2213

private

Definition at line 113 of file MuonSeedCreator.h.

Referenced by estimatePtOverlap().

std::vector<double> MuonSeedCreator::OL2222

private

Definition at line 114 of file MuonSeedCreator.h.

Referenced by estimatePtOverlap(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::OL_1213

private

Definition at line 162 of file MuonSeedCreator.h.

Referenced by estimatePtOverlap(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::OL_1222

private

Definition at line 163 of file MuonSeedCreator.h.

Referenced by estimatePtOverlap(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::OL_1232

private

Definition at line 164 of file MuonSeedCreator.h.

Referenced by estimatePtOverlap(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::OL_2213

private

Definition at line 165 of file MuonSeedCreator.h.

Referenced by estimatePtOverlap(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::OL_2222

private

Definition at line 166 of file MuonSeedCreator.h.

Referenced by estimatePtOverlap(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB10

private

Definition at line 125 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB11

private

Definition at line 126 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB12

private

Definition at line 127 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB20

private

Definition at line 128 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB21

private

Definition at line 129 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB22

private

Definition at line 130 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB30

private

Definition at line 131 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB31

private

Definition at line 132 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB32

private

Definition at line 133 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB_10S

private

Definition at line 168 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB_11S

private

Definition at line 169 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB_12S

private

Definition at line 170 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB_20S

private

Definition at line 171 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB_21S

private

Definition at line 172 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB_22S

private

Definition at line 173 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB_30S

private

Definition at line 174 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB_31S

private

Definition at line 175 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SMB_32S

private

Definition at line 176 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SME11

private

Definition at line 116 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SME12

private

Definition at line 117 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SME13

private

Definition at line 118 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SME21

private

Definition at line 119 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SME22

private

Definition at line 120 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SME31

private

Definition at line 121 of file MuonSeedCreator.h.

Referenced by MuonSeedCreator().

std::vector<double> MuonSeedCreator::SME32

private

Definition at line 122 of file MuonSeedCreator.h.

Referenced by MuonSeedCreator().

std::vector<double> MuonSeedCreator::SME41

private

Definition at line 123 of file MuonSeedCreator.h.

Referenced by MuonSeedCreator().

std::vector<double> MuonSeedCreator::SME_11S

private

Definition at line 178 of file MuonSeedCreator.h.

Referenced by MuonSeedCreator().

std::vector<double> MuonSeedCreator::SME_12S

private

Definition at line 179 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SME_13S

private

Definition at line 180 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SME_21S

private

Definition at line 181 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

std::vector<double> MuonSeedCreator::SME_22S

private

Definition at line 182 of file MuonSeedCreator.h.

Referenced by estimatePtSingle(), and MuonSeedCreator().

double MuonSeedCreator::sysError

private

Definition at line 84 of file MuonSeedCreator.h.

Referenced by MuonSeedCreator().

float MuonSeedCreator::theMaxMomentum

private

Definition at line 80 of file MuonSeedCreator.h.

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

float MuonSeedCreator::theMinMomentum

private

Definition at line 79 of file MuonSeedCreator.h.

Referenced by createSeed(), and MuonSeedCreator().