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
313 for (
unsigned l=0;
l<seg.size();
l++) {
334 unsigned size = seg.size();
335 if (size < 2)
return;
343 std::vector<double> ptEstimate;
344 std::vector<double> sptEstimate;
353 int layer0 = layers[0];
354 segPos[0] = seg[0]->globalPosition();
355 float eta = fabs( segPos[0].
eta() );
371 unsigned idx1 =
size;
375 int layer1 = layers[idx1];
376 if (layer0 == layer1)
continue;
377 segPos[1] = seg[idx1]->globalPosition();
379 double dphi = segPos[0].
phi() - segPos[1].
phi();
381 double temp_dphi = dphi;
384 if (temp_dphi < 0.) {
385 temp_dphi = -1.0*temp_dphi;
390 if (temp_dphi < 0.0001 ) {
394 ptEstimate.push_back( pt*sign );
395 sptEstimate.push_back( spt );
398 if ( layer0 == 0 && temp_dphi >= 0.0001 ) {
407 else if ( layer1 == 2 ) {
413 else if ( layer1 == 3 ) {
424 ptEstimate.push_back( pt*sign );
425 sptEstimate.push_back( spt );
439 else if ( layer1 == 3 ) {
450 ptEstimate.push_back( pt*sign );
451 sptEstimate.push_back( spt );
455 if ( layer0 == 2 && temp_dphi > 0.0001 ) {
475 ptEstimate.push_back( pt*sign );
476 sptEstimate.push_back( spt );
480 if ( layer0 == 3 && temp_dphi > 0.0001 ) {
485 ptEstimate.push_back( pt*sign );
486 sptEstimate.push_back( spt );
493 if ( ptEstimate.size() > 0 )
weightedPt( ptEstimate, sptEstimate, thePt, theSpt);
506 unsigned size = seg.size();
507 if (size < 2)
return;
509 std::vector<double> ptEstimate;
510 std::vector<double> sptEstimate;
519 int layer0 = layers[0];
520 segPos[0] = seg[0]->globalPosition();
521 float eta = fabs(segPos[0].
eta());
530 for (
unsigned idx1 = 1; idx1 <
size ; ++idx1 ) {
532 int layer1 = layers[idx1];
533 segPos[1] = seg[idx1]->globalPosition();
538 double dphi = segPos[0].
phi() - segPos[1].
phi();
539 double temp_dphi = dphi;
544 if (temp_dphi < 0.) {
545 temp_dphi = -temp_dphi;
549 if (temp_dphi < 0.0001 ) {
553 ptEstimate.push_back( pt*sign );
554 sptEstimate.push_back( spt );
559 if (layer0 == -1 && temp_dphi > 0.0001 ) {
570 else if (layer1 == -3) {
587 ptEstimate.push_back( pt*sign );
588 sptEstimate.push_back( spt );
592 if (layer0 == -2 && temp_dphi > 0.0001 ) {
609 ptEstimate.push_back( pt*sign );
610 sptEstimate.push_back( spt );
614 if (layer0 == -3 && temp_dphi > 0.0001 ) {
621 ptEstimate.push_back( pt*sign );
622 sptEstimate.push_back( spt );
629 if (ptEstimate.size() > 0 )
weightedPt( ptEstimate, sptEstimate, thePt, theSpt);
640 int size = layers.size();
646 std::vector<int> layersCSC;
648 std::vector<int> layersDT;
651 for (
unsigned j = 0;
j < layers.size(); ++
j ) {
652 if ( layers[
j] > -1 ) {
653 segCSC.push_back(seg[
j]);
654 layersCSC.push_back(layers[j]);
657 segDT.push_back(seg[
j]);
658 layersDT.push_back(layers[j]);
662 std::vector<double> ptEstimate;
663 std::vector<double> sptEstimate;
666 int layer0 = layers[0];
667 segPos[0] = seg[0]->globalPosition();
668 float eta = fabs(segPos[0].
eta());
671 if ( segDT.size() > 0 && segCSC.size() > 0 ) {
672 int layer1 = layers[size-1];
673 segPos[1] = seg[size-1]->globalPosition();
675 double dphi = segPos[0].
phi() - segPos[1].
phi();
676 double temp_dphi = dphi;
681 if (temp_dphi < 0.) {
682 temp_dphi = -temp_dphi;
686 if (temp_dphi < 0.0001 ) {
690 ptEstimate.push_back( thePt*sign );
691 sptEstimate.push_back( theSpt );
695 if ( layer0 == -1 && temp_dphi > 0.0001 ) {
703 else if ( layer1 == 2) {
714 ptEstimate.push_back(thePt*sign);
715 sptEstimate.push_back(theSpt);
718 if ( layer0 == -2 && temp_dphi > 0.0001 ) {
724 ptEstimate.push_back(thePt*sign);
725 sptEstimate.push_back(theSpt);
736 if ( segDT.size() > 1 ) {
738 ptEstimate.push_back(thePt);
739 sptEstimate.push_back(theSpt);
759 if (ptEstimate.size() > 0 )
weightedPt( ptEstimate, sptEstimate, thePt, theSpt);
773 double eta = segPos.
eta();
778 double cosDpsi = (gv.
x()*segPos.
x() + gv.
y()*segPos.
y());
779 cosDpsi /=
sqrt(segPos.
x()*segPos.
x() + segPos.
y()*segPos.
y());
780 cosDpsi /=
sqrt(gv.
x()*gv.
x() + gv.
y()*gv.
y());
782 double axb = ( segPos.
x()*gv.
y() ) - ( segPos.
y()*gv.
x() ) ;
783 double sign = (axb < 0.) ? 1.0 : -1.0;
785 double dpsi = fabs(acos(cosDpsi)) ;
786 if ( dpsi > 1.570796 ) {
787 dpsi = 3.141592 - dpsi;
790 if (fabs(dpsi) < 0.00005) {
795 if ( layers[0] == -1 ) {
797 if ( fabs(eta) < 0.3 ) {
803 if ( fabs(eta) >= 0.3 && fabs(eta) < 0.82 ) {
809 if ( fabs(eta) >= 0.82 && fabs(eta) < 1.2 ) {
815 if ( layers[0] == 1 ) {
817 if ( fabs(eta) > 0.92 && fabs(eta) < 1.16 ) {
823 if ( fabs(eta) >= 1.16 && fabs(eta) <= 1.6 ) {
829 if ( layers[0] == 0 ) {
831 if ( fabs(eta) > 1.6 ) {
838 if ( layers[0] == -2 ) {
840 if ( fabs(eta) < 0.25 ) {
846 if ( fabs(eta) >= 0.25 && fabs(eta) < 0.72 ) {
852 if ( fabs(eta) >= 0.72 && fabs(eta) < 1.04 ) {
858 if ( layers[0] == 2 ) {
860 if ( fabs(eta) > 0.95 && fabs(eta) <= 1.6 ) {
866 if ( fabs(eta) > 1.6 && fabs(eta) < 2.45 ) {
874 if ( layers[0] == -3 ) {
876 if ( fabs(eta) <= 0.22 ) {
882 if ( fabs(eta) > 0.22 && fabs(eta) <= 0.6 ) {
888 if ( fabs(eta) > 0.6 && fabs(eta) < 0.95 ) {
894 thePt = fabs(thePt)*
sign;
895 theSpt = fabs(theSpt);
903 if ( NShowers > 2 && thePt < 300. ) {
907 if ( NShowers == 2 && NShowerSegments > 11 && thePt < 150. ) {
911 if ( NShowers == 2 && NShowerSegments <= 11 && thePt < 50. ) {
915 if ( NShowers == 1 && NShowerSegments <= 5 && thePt < 10. ) {
931 int size = ptEstimate.size();
935 thePt = ptEstimate[0];
936 theSpt = sptEstimate[0];
944 for (
unsigned j = 0;
j < ptEstimate.size();
j++ ) {
946 if ( ptEstimate[
j] < 0. ) {
948 charge -= 1. * (ptEstimate[
j]*ptEstimate[
j]) / (sptEstimate[
j]*sptEstimate[
j] );
950 charge += 1. * (ptEstimate[
j]*ptEstimate[
j]) / (sptEstimate[
j]*sptEstimate[
j] );
962 double weightPtSum = 0.;
963 double sigmaSqr_sum = 0.;
967 for (
unsigned j = 0;
j < ptEstimate.size(); ++
j ) {
971 sigmaSqr_sum += 1.0 / (sptEstimate[
j]*sptEstimate[
j]);
972 weightPtSum += fabs(ptEstimate[
j])/(sptEstimate[
j]*sptEstimate[
j]);
984 thePt = (charge*weightPtSum) / sigmaSqr_sum;
985 theSpt =
sqrt( 1.0 / sigmaSqr_sum ) ;
994 double h = fabs(eta);
995 double estPt = ( vPara[0] + vPara[1]*h + vPara[2]*h*
h ) / dPhi;
996 double estSPt = ( vPara[3] + vPara[4]*h + vPara[5]*h*
h ) * estPt;
997 std::vector<double> paraPt ;
998 paraPt.push_back( estPt );
999 paraPt.push_back( estSPt ) ;
1009 double oPhi = 1./dphi ;
1010 dphi = dphi /( 1. + t1/( oPhi + 10. ) ) ;
std::vector< double > OL2213
T getParameter(std::string const &) const
std::vector< double > DT13_2
std::vector< double > SMB32
std::vector< double > CSC23_1
std::vector< double > SMB_30S
std::vector< double > SMB_21S
std::vector< double > DT14_2
std::vector< LayerSetAndLayers > layers(const SeedingLayerSetsHits &sets)
std::vector< double > DT23_2
std::vector< double > DT34_1
std::vector< double > DT14_1
std::vector< double > SME_22S
std::vector< double > SMB22
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
std::vector< double > OL1232
std::vector< double > DT34_2
std::vector< double > CSC24
void estimatePtDT(const SegmentContainer &seg, const std::vector< int > &layers, double &pt, double &spt)
Estimate transverse momentum of track from DT measurements.
void weightedPt(const std::vector< double > &ptEstimate, const std::vector< double > &sptEstimate, double &ptAvg, double &sptAvg)
Compute weighted mean pt from different pt estimators.
std::vector< double > CSC23_2
std::vector< double > SMB20
Geom::Phi< T > phi() const
Global3DPoint GlobalPoint
std::vector< double > CSC14
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
TrajectorySeed createSeed(int type, const SegmentContainer &seg, const std::vector< int > &layers, int NShower, int NShowerSeg)
Create a seed from set of segments.
MuonTransientTrackingRecHit::MuonRecHitContainer SegmentContainer
std::vector< double > SME21
void estimatePtOverlap(const SegmentContainer &seg, const std::vector< int > &layers, double &pt, double &spt)
Estimate transverse momentum of track from CSC + DT measurements.
std::vector< double > CSC13_3
MuonSeedCreator(const edm::ParameterSet &pset)
Constructor.
const MagneticField * BField
void estimatePtSingle(const SegmentContainer &seg, const std::vector< int > &layers, double &pt, double &spt)
Estimate transverse momentum of track from single segment.
std::vector< double > CSC34
void estimatePtCSC(const SegmentContainer &seg, const std::vector< int > &layers, double &pt, double &spt)
Estimate transverse momentum of track from CSC measurements.
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
std::vector< double > SMB_20S
std::vector< double > SMB_11S
std::vector< double > OL1213
std::vector< double > getPt(const std::vector< double > &vParameters, double eta, double dPhi)
Compute pt from parameters.
CLHEP::HepVector AlgebraicVector
std::vector< double > SMB12
std::vector< double > SME_21S
std::vector< double > SME_12S
std::vector< double > OL2222
TEveGeoShape * clone(const TEveElement *element, TEveElement *parent)
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
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