31 #include "gsl/gsl_statistics.h"
65 OL1213 = pset.
getParameter<std::vector<double> >(
"OL_1213");
168 if (type == 3 )
estimatePtDT(seg, layers, ptmean, sptmean);
180 if (ptmean < 0.) charge = -1.0;
195 double chi2_dof = 9999.0;
196 unsigned int ini_seg = 0;
198 if ( type == 5 ) ini_seg = 1;
199 for (
size_t i = ini_seg ;
i < seg.size();
i++) {
200 double dof =
static_cast<double>(seg[
i]->degreesOfFreedom());
201 if ( chi2_dof < ( seg[
i]->
chi2()/dof ) )
continue;
202 chi2_dof = seg[
i]->chi2() / dof ;
203 best_seg =
static_cast<int>(
i);
207 if ( type==1 || type==5 || type== 4) {
213 segPos = seg[
last]->localPosition();
219 polar *= fabs(ptmean)/polar.
perp();
222 int chargeI =
static_cast<int>(
charge);
225 p_err = (sptmean*sptmean)/(polar.mag()*polar.mag()*ptmean*ptmean) ;
226 mat = seg[
last]->parametersError().similarityT( seg[last]->projectionMatrix() );
229 mat[0][0] = mat[0][0]/fabs(
tan(mom.
theta()) );
230 mat[1][1] = mat[1][1]/fabs(
tan(mom.
theta()) );
231 mat[3][3] = 2.25*mat[3][3];
232 mat[4][4] = 2.25*mat[4][4];
235 mat[0][0] = mat[0][0]/fabs(
tan(mom.
theta()) );
236 mat[1][1] = mat[1][1]/fabs(
tan(mom.
theta()) );
237 mat[2][2] = 2.25*mat[2][2];
238 mat[3][3] = 2.25*mat[3][3];
239 mat[4][4] = 2.25*mat[4][4];
241 double dh = fabs( seg[last]->globalPosition().
eta() ) - 1.6 ;
242 if ( fabs(dh) < 0.1 && type == 1 ) {
243 mat[1][1] = 4.*mat[1][1];
244 mat[2][2] = 4.*mat[2][2];
245 mat[3][3] = 9.*mat[3][3];
246 mat[4][4] = 9.*mat[4][4];
259 segPos = seg[
last]->localPosition();
270 polar *= fabs(ptmean)/polar.
perp();
273 int chargeI =
static_cast<int>(
charge);
276 p_err = (sptmean*sptmean)/(polar.mag()*polar.mag()*ptmean*ptmean) ;
277 mat = seg[
last]->parametersError().similarityT( seg[last]->projectionMatrix() );
284 float Geta = gp.
eta();
286 std::cout <<
"Type " << type <<
" Nsegments " << layers.size() <<
" ";
287 std::cout <<
"pt " << ptmean <<
" errpt " << sptmean
288 <<
" eta " << Geta <<
" charge " << charge
311 std::auto_ptr<PTrajectoryStateOnDet> seedTSOS(tsTransform.
persistentState( tsos,
id.rawId()));
314 for (
unsigned l=0;
l<seg.size();
l++) {
335 unsigned size = seg.size();
336 if (size < 2)
return;
344 std::vector<double> ptEstimate;
345 std::vector<double> sptEstimate;
354 int layer0 = layers[0];
355 segPos[0] = seg[0]->globalPosition();
356 float eta = fabs( segPos[0].
eta() );
372 unsigned idx1 =
size;
376 int layer1 = layers[idx1];
377 if (layer0 == layer1)
continue;
378 segPos[1] = seg[idx1]->globalPosition();
380 double dphi = segPos[0].
phi() - segPos[1].
phi();
382 double temp_dphi = dphi;
385 if (temp_dphi < 0.) {
386 temp_dphi = -1.0*temp_dphi;
391 if (temp_dphi < 0.0001 ) {
395 ptEstimate.push_back( pt*sign );
396 sptEstimate.push_back( spt );
399 if ( layer0 == 0 && temp_dphi >= 0.0001 ) {
408 else if ( layer1 == 2 ) {
414 else if ( layer1 == 3 ) {
425 ptEstimate.push_back( pt*sign );
426 sptEstimate.push_back( spt );
440 else if ( layer1 == 3 ) {
451 ptEstimate.push_back( pt*sign );
452 sptEstimate.push_back( spt );
456 if ( layer0 == 2 && temp_dphi > 0.0001 ) {
476 ptEstimate.push_back( pt*sign );
477 sptEstimate.push_back( spt );
481 if ( layer0 == 3 && temp_dphi > 0.0001 ) {
486 ptEstimate.push_back( pt*sign );
487 sptEstimate.push_back( spt );
494 if ( ptEstimate.size() > 0 )
weightedPt( ptEstimate, sptEstimate, thePt, theSpt);
507 unsigned size = seg.size();
508 if (size < 2)
return;
510 std::vector<double> ptEstimate;
511 std::vector<double> sptEstimate;
520 int layer0 = layers[0];
521 segPos[0] = seg[0]->globalPosition();
522 float eta = fabs(segPos[0].
eta());
531 for (
unsigned idx1 = 1; idx1 <
size ; ++idx1 ) {
533 int layer1 = layers[idx1];
534 segPos[1] = seg[idx1]->globalPosition();
539 double dphi = segPos[0].
phi() - segPos[1].
phi();
540 double temp_dphi = dphi;
545 if (temp_dphi < 0.) {
546 temp_dphi = -temp_dphi;
550 if (temp_dphi < 0.0001 ) {
554 ptEstimate.push_back( pt*sign );
555 sptEstimate.push_back( spt );
560 if (layer0 == -1 && temp_dphi > 0.0001 ) {
571 else if (layer1 == -3) {
588 ptEstimate.push_back( pt*sign );
589 sptEstimate.push_back( spt );
593 if (layer0 == -2 && temp_dphi > 0.0001 ) {
610 ptEstimate.push_back( pt*sign );
611 sptEstimate.push_back( spt );
615 if (layer0 == -3 && temp_dphi > 0.0001 ) {
622 ptEstimate.push_back( pt*sign );
623 sptEstimate.push_back( spt );
630 if (ptEstimate.size() > 0 )
weightedPt( ptEstimate, sptEstimate, thePt, theSpt);
641 int size = layers.size();
647 std::vector<int> layersCSC;
649 std::vector<int> layersDT;
652 for (
unsigned j = 0;
j < layers.size(); ++
j ) {
653 if ( layers[
j] > -1 ) {
654 segCSC.push_back(seg[
j]);
655 layersCSC.push_back(layers[j]);
658 segDT.push_back(seg[
j]);
659 layersDT.push_back(layers[j]);
663 std::vector<double> ptEstimate;
664 std::vector<double> sptEstimate;
667 int layer0 = layers[0];
668 segPos[0] = seg[0]->globalPosition();
669 float eta = fabs(segPos[0].
eta());
672 if ( segDT.size() > 0 && segCSC.size() > 0 ) {
673 int layer1 = layers[size-1];
674 segPos[1] = seg[size-1]->globalPosition();
676 double dphi = segPos[0].
phi() - segPos[1].
phi();
677 double temp_dphi = dphi;
682 if (temp_dphi < 0.) {
683 temp_dphi = -temp_dphi;
687 if (temp_dphi < 0.0001 ) {
691 ptEstimate.push_back( thePt*sign );
692 sptEstimate.push_back( theSpt );
696 if ( layer0 == -1 && temp_dphi > 0.0001 ) {
704 else if ( layer1 == 2) {
715 ptEstimate.push_back(thePt*sign);
716 sptEstimate.push_back(theSpt);
719 if ( layer0 == -2 && temp_dphi > 0.0001 ) {
725 ptEstimate.push_back(thePt*sign);
726 sptEstimate.push_back(theSpt);
737 if ( segDT.size() > 1 ) {
739 ptEstimate.push_back(thePt);
740 sptEstimate.push_back(theSpt);
760 if (ptEstimate.size() > 0 )
weightedPt( ptEstimate, sptEstimate, thePt, theSpt);
774 double eta = segPos.
eta();
779 double cosDpsi = (gv.
x()*segPos.
x() + gv.
y()*segPos.
y());
780 cosDpsi /=
sqrt(segPos.
x()*segPos.
x() + segPos.
y()*segPos.
y());
781 cosDpsi /=
sqrt(gv.
x()*gv.
x() + gv.
y()*gv.
y());
783 double axb = ( segPos.
x()*gv.
y() ) - ( segPos.
y()*gv.
x() ) ;
784 double sign = (axb < 0.) ? 1.0 : -1.0;
786 double dpsi = fabs(acos(cosDpsi)) ;
787 if ( dpsi > 1.570796 ) {
788 dpsi = 3.141592 - dpsi;
791 if (fabs(dpsi) < 0.00005) {
796 if ( layers[0] == -1 ) {
798 if ( fabs(eta) < 0.3 ) {
804 if ( fabs(eta) >= 0.3 && fabs(eta) < 0.82 ) {
810 if ( fabs(eta) >= 0.82 && fabs(eta) < 1.2 ) {
816 if ( layers[0] == 1 ) {
818 if ( fabs(eta) > 0.92 && fabs(eta) < 1.16 ) {
824 if ( fabs(eta) >= 1.16 && fabs(eta) <= 1.6 ) {
830 if ( layers[0] == 0 ) {
832 if ( fabs(eta) > 1.6 ) {
839 if ( layers[0] == -2 ) {
841 if ( fabs(eta) < 0.25 ) {
847 if ( fabs(eta) >= 0.25 && fabs(eta) < 0.72 ) {
853 if ( fabs(eta) >= 0.72 && fabs(eta) < 1.04 ) {
859 if ( layers[0] == 2 ) {
861 if ( fabs(eta) > 0.95 && fabs(eta) <= 1.6 ) {
867 if ( fabs(eta) > 1.6 && fabs(eta) < 2.45 ) {
875 if ( layers[0] == -3 ) {
877 if ( fabs(eta) <= 0.22 ) {
883 if ( fabs(eta) > 0.22 && fabs(eta) <= 0.6 ) {
889 if ( fabs(eta) > 0.6 && fabs(eta) < 0.95 ) {
895 thePt = fabs(thePt)*sign;
896 theSpt = fabs(theSpt);
904 if ( NShowers > 2 && thePt < 300. ) {
908 if ( NShowers == 2 && NShowerSegments > 11 && thePt < 150. ) {
912 if ( NShowers == 2 && NShowerSegments <= 11 && thePt < 50. ) {
916 if ( NShowers == 1 && NShowerSegments <= 5 && thePt < 10. ) {
932 int size = ptEstimate.size();
936 thePt = ptEstimate[0];
937 theSpt = sptEstimate[0];
945 for (
unsigned j = 0;
j < ptEstimate.size();
j++ ) {
947 if ( ptEstimate[
j] < 0. ) {
949 charge -= 1. * (ptEstimate[
j]*ptEstimate[
j]) / (sptEstimate[
j]*sptEstimate[
j] );
951 charge += 1. * (ptEstimate[
j]*ptEstimate[
j]) / (sptEstimate[
j]*sptEstimate[
j] );
963 double weightPtSum = 0.;
964 double sigmaSqr_sum = 0.;
968 for (
unsigned j = 0;
j < ptEstimate.size(); ++
j ) {
972 sigmaSqr_sum += 1.0 / (sptEstimate[
j]*sptEstimate[
j]);
973 weightPtSum += fabs(ptEstimate[
j])/(sptEstimate[
j]*sptEstimate[
j]);
985 thePt = (charge*weightPtSum) / sigmaSqr_sum;
986 theSpt =
sqrt( 1.0 / sigmaSqr_sum ) ;
995 double h = fabs(eta);
996 double estPt = ( vPara[0] + vPara[1]*h + vPara[2]*h*
h ) / dPhi;
997 double estSPt = ( vPara[3] + vPara[4]*h + vPara[5]*h*
h ) * estPt;
998 std::vector<double> paraPt ;
999 paraPt.push_back( estPt );
1000 paraPt.push_back( estSPt ) ;
1010 double oPhi = 1./dphi ;
1011 dphi = dphi /( 1. + t1/( oPhi + 10. ) ) ;
std::vector< double > OL2213
T getParameter(std::string const &) const
std::vector< double > DT13_2
TrajectorySeed createSeed(int type, SegmentContainer seg, std::vector< int > layers, int NShower, int NShowerSeg)
Create a seed from set of segments.
std::vector< double > SMB32
void estimatePtCSC(SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
Estimate transverse momentum of track from CSC measurements.
std::vector< double > CSC23_1
std::vector< double > SMB_30S
std::vector< double > SMB_21S
std::vector< double > DT14_2
std::vector< double > DT23_2
std::vector< double > DT34_1
std::vector< double > DT14_1
std::vector< double > SME_22S
std::vector< double > SMB22
std::vector< double > OL1232
std::vector< double > DT34_2
std::vector< double > CSC24
std::vector< double > CSC23_2
std::vector< double > SMB20
Geom::Phi< T > phi() const
std::vector< double > getPt(std::vector< double > vParameters, double eta, double dPhi)
Compute pt from parameters.
Global3DPoint GlobalPoint
std::vector< double > CSC14
void estimatePtSingle(SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
Estimate transverse momentum of track from single segment.
void estimatePtOverlap(SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
Estimate transverse momentum of track from CSC + DT measurements.
std::vector< double > DT34
std::vector< double > CSC12_2
std::vector< double > CSC13_2
std::vector< double > CSC12_1
std::vector< double > SME12
std::vector< double > SME22
void weightedPt(std::vector< double > ptEstimate, std::vector< double > sptEstimate, double &ptAvg, double &sptAvg)
Compute weighted mean pt from different pt estimators.
MuonTransientTrackingRecHit::MuonRecHitContainer SegmentContainer
std::vector< double > SME21
std::vector< double > CSC13_3
MuonSeedCreator(const edm::ParameterSet &pset)
Constructor.
const MagneticField * BField
std::vector< double > CSC34
std::vector< double > SME13
Geom::Theta< T > theta() const
std::vector< double > DT13
std::vector< double > CSC12_3
std::vector< double > SMB10
std::vector< double > CSC12
std::vector< double > DT12_2
std::vector< double > OL_1222
std::vector< double > OL1222
double dPhi(double phi1, double phi2)
void estimatePtShowering(int &NShowers, int &NShowerSeg, double &pt, double &spt)
Estimate transverse momentum of track from showering segment.
std::vector< double > CSC23
std::vector< double > CSC13
std::vector< double > DT23_1
Tan< T >::type tan(const T &t)
std::vector< double > SMB30
std::vector< double > SMB_10S
std::vector< double > DT24_2
std::vector< double > SMB_31S
std::vector< double > SMB_32S
std::vector< double > DT12_1
std::vector< double > CSC14_3
std::vector< double > DT12
void estimatePtDT(SegmentContainer seg, std::vector< int > layers, double &pt, double &spt)
Estimate transverse momentum of track from DT measurements.
std::vector< double > SMB_20S
std::vector< double > SMB_11S
std::vector< double > OL1213
CLHEP::HepVector AlgebraicVector
std::vector< double > SMB12
std::vector< double > SME_21S
std::vector< double > SME_12S
std::vector< double > OL2222
std::vector< double > SME_11S
std::vector< double > DT14
std::vector< double > DT13_1
std::vector< double > SME31
std::vector< double > DT24
std::vector< double > OL_2213
std::vector< double > DT23
std::vector< double > SME32
std::vector< double > CSC24_1
std::vector< double > SMB21
CLHEP::HepSymMatrix AlgebraicSymMatrix
std::vector< double > CSC01_1
double scaledPhi(double dphi, double t1)
Scale the dPhi from segment position.
~MuonSeedCreator()
Destructor.
std::vector< double > SME_13S
std::vector< double > OL_2222
std::vector< double > SMB31
std::vector< double > CSC03
std::vector< double > CSC02
std::vector< double > DT24_1
std::vector< double > CSC34_1
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
std::vector< double > SME11
std::vector< double > SMB11
std::vector< double > OL_1232
std::vector< double > SMB_12S
std::vector< double > SME41
tuple size
Write out results.
std::vector< double > OL_1213
std::vector< double > SMB_22S