|
|
|
#include <MuonDTSeedFromRecHits.h>
Public Member Functions | |
| ConstMuonRecHitPointer | bestBarrelHit (const MuonRecHitContainer &barrelHits) const |
| MuonDTSeedFromRecHits () | |
| virtual TrajectorySeed | seed () const |
Private Member Functions | |
| float | bestEta () const |
| void | computeBestPt (double *pt, double *spt, float &ptmean, float &sptmean) const |
| void | computeMean (const double *pt, const double *weights, int sz, bool tossOutlyers, float &ptmean, float &sptmean) const |
| void | computePtWithoutVtx (double *pt, double *spt) const |
| void | computePtWithVtx (double *pt, double *spt) const |
Generate a seed starting from a list of RecHits make use of TrajectorySeed from CommonDet
Definition at line 20 of file MuonDTSeedFromRecHits.h.
| MuonDTSeedFromRecHits::MuonDTSeedFromRecHits | ( | ) |
Definition at line 29 of file MuonDTSeedFromRecHits.cc.
: MuonSeedFromRecHits() { }
| MuonDTSeedFromRecHits::ConstMuonRecHitPointer MuonDTSeedFromRecHits::bestBarrelHit | ( | const MuonRecHitContainer & | barrelHits | ) | const |
Definition at line 82 of file MuonDTSeedFromRecHits.cc.
Referenced by MuonOverlapSeedFromRecHits::bestHit(), and seed().
{
int alt_npt = 0;
int best_npt = 0;
int cur_npt = 0;
MuonRecHitPointer best = 0;
MuonRecHitPointer alter=0;
for (MuonRecHitContainer::const_iterator iter=barrelHits.begin();
iter!=barrelHits.end(); iter++ ) {
bool hasZed = ((*iter)->projectionMatrix()[1][1]==1);
cur_npt = 0;
vector<TrackingRecHit*> slrhs=(*iter)->recHits();
for (vector<TrackingRecHit*>::const_iterator slrh=slrhs.begin(); slrh!=slrhs.end(); ++slrh ) {
cur_npt+=(*slrh)->recHits().size();
}
float radius1 = (*iter)->det()->position().perp();
if (hasZed) {
if ( cur_npt > best_npt ) {
best = (*iter);
best_npt = cur_npt;
}
else if ( best && cur_npt == best_npt ) {
float radius2 = best->det()->position().perp();
if (radius1 < radius2) {
best = (*iter);
best_npt = cur_npt;
}
}
}
if ( cur_npt > alt_npt ) {
alter = (*iter);
alt_npt = cur_npt;
}
else if ( alter && cur_npt == alt_npt ) {
float radius2 = alter->det()->position().perp();
if (radius1 < radius2) {
alter = (*iter);
alt_npt = cur_npt;
}
}
}
if ( !best ) best = alter;
return best;
}
| float MuonDTSeedFromRecHits::bestEta | ( | ) | const [private] |
Definition at line 135 of file MuonDTSeedFromRecHits.cc.
References query::result, and MuonSeedFromRecHits::theRhits.
Referenced by computePtWithoutVtx(), and computePtWithVtx().
{
int Maxseg = 0;
float Msdeta = 100.;
float result = (*theRhits.begin())->globalPosition().eta();
for (MuonRecHitContainer::const_iterator iter=theRhits.begin(); iter!=theRhits.end(); iter++ ) {
float eta1= (*iter)->globalPosition().eta();
int Nseg = 0;
float sdeta = .0;
for (MuonRecHitContainer::const_iterator iter2=theRhits.begin(); iter2!=theRhits.end(); iter2++ ) {
if ( iter2 == iter ) continue;
float eta2 = (*iter2)->globalPosition().eta();
if ( fabs (eta1-eta2) > .2 ) continue;
Nseg++;
sdeta += fabs (eta1-eta2);
if ( Nseg > Maxseg || (Nseg == Maxseg && sdeta < Msdeta) ) {
Maxseg = Nseg;
Msdeta = sdeta;
result = eta1;
}
}
} // +v.
return result;
}
| void MuonDTSeedFromRecHits::computeBestPt | ( | double * | pt, |
| double * | spt, | ||
| float & | ptmean, | ||
| float & | sptmean | ||
| ) | const [private] |
just one pt estimate: use it.
more than a pt estimate, do all the job.
Definition at line 260 of file MuonDTSeedFromRecHits.cc.
References computeMean(), i, LogTrace, metname, and mathSSE::sqrt().
Referenced by seed().
{
const std::string metname = "Muon|RecoMuon|MuonDTSeedFromRecHits";
int nTotal=8;
int igood=-1;
for (int i=0;i<=7;i++) {
if(pt[i]==0) {
spt[i]=0.;
nTotal--;
} else {
igood=i;
}
}
if (nTotal==1) {
ptmean=pt[igood];
sptmean=1/sqrt(spt[igood]);
} else if (nTotal==0) {
// No estimate (e.g. only one rechit!)
ptmean=50;
sptmean=30;
} else {
// calculate pt with vertex
float ptvtx=0.0;
float sptvtx=0.0;
computeMean(pt, spt, 2, false, ptvtx, sptvtx);
LogTrace(metname) << " GSL: Pt w vtx :" << ptvtx << "+/-" <<
sptvtx << endl;
// FIXME: temp hack
if(ptvtx != 0.) {
ptmean = ptvtx;
sptmean = sptvtx;
return;
//
}
// calculate pt w/o vertex
float ptMB=0.0;
float sptMB=0.0;
computeMean(pt+2, spt+2, 6, false, ptMB, sptMB);
LogTrace(metname) << " GSL: Pt w/o vtx :" << ptMB << "+/-" <<
sptMB << endl;
// decide wheter the muon comes or not from vertex
float ptpool=0.0;
if((ptvtx+ptMB)!=0.0) ptpool=(ptvtx-ptMB)/(ptvtx+ptMB);
bool fromvtx=true;
if(fabs(ptpool)>0.2) fromvtx=false;
LogTrace(metname) << "From vtx? " <<fromvtx << " ptpool "<< ptpool << endl;
// now choose the "right" pt => weighted mean
computeMean(pt, spt, 8, true, ptmean, sptmean);
LogTrace(metname) << " GSL Pt :" << ptmean << "+/-" << sptmean << endl;
}
}
| void MuonDTSeedFromRecHits::computeMean | ( | const double * | pt, |
| const double * | weights, | ||
| int | sz, | ||
| bool | tossOutlyers, | ||
| float & | ptmean, | ||
| float & | sptmean | ||
| ) | const [private] |
Definition at line 324 of file MuonDTSeedFromRecHits.cc.
References i, n, and mathSSE::sqrt().
Referenced by computeBestPt().
{
int n=0;
double ptTmp[8];
double wtTmp[8];
assert(sz<=8);
for (int i=0; i<sz; ++i) {
ptTmp[i] = 0.;
wtTmp[i] = 0;
if (pt[i]!=0) {
ptTmp[n]=pt[i];
wtTmp[n]=weights[i];
++n;
}
}
if(n != 0)
{
if (n==1) {
ptmean=ptTmp[0];
sptmean=1/sqrt(wtTmp[0]);
} else {
ptmean = gsl_stats_wmean(wtTmp, 1, ptTmp, 1, n);
sptmean = sqrt( gsl_stats_wvariance_m (wtTmp, 1, ptTmp, 1, n, ptmean) );
}
if(tossOutlyers)
{
// Recompute mean with a cut at 3 sigma
for ( int nm =0; nm<8; nm++ ){
if ( ptTmp[nm]!=0 && fabs(ptTmp[nm]-ptmean)>3*(sptmean) ) {
wtTmp[nm]=0.;
}
}
ptmean = gsl_stats_wmean(wtTmp, 1, ptTmp, 1, n);
sptmean = sqrt( gsl_stats_wvariance_m (wtTmp, 1, ptTmp, 1, n, ptmean) );
}
}
}
| void MuonDTSeedFromRecHits::computePtWithoutVtx | ( | double * | pt, |
| double * | spt | ||
| ) | const [private] |
Definition at line 202 of file MuonDTSeedFromRecHits.cc.
References bestEta(), DTChamberId, MuonSeedPtExtractor::pT_extract(), relativeConstraints::station, MuonSeedFromRecHits::thePtExtractor, MuonSeedFromRecHits::theRhits, and tmp.
Referenced by seed().
{
float eta0 = bestEta();
for (MuonRecHitContainer::const_iterator iter=theRhits.begin();
iter!=theRhits.end(); iter++ ) {
//+vvp !:
float eta1= (*iter)->globalPosition().eta();
if ( fabs (eta1-eta0) > .2 ) continue; // !!! +vvp
for (MuonRecHitContainer::const_iterator iter2=theRhits.begin();
iter2!=iter; iter2++ ) {
//+vvp !:
float eta2= (*iter2)->globalPosition().eta();
if ( fabs (eta2-eta0) > .2 ) continue; // !!! +vvp
unsigned int stat1 = DTChamberId((*iter)->geographicalId()).station();
unsigned int stat2 = DTChamberId((*iter2)->geographicalId()).station();
if ( stat1>stat2) {
int tmp = stat1;
stat1 = stat2;
stat2 = tmp;
}
unsigned int st = stat1*10+stat2;
double thispt = thePtExtractor->pT_extract(*iter, *iter2)[0];
switch (st) {
case 12 : {//MB1-MB2
pt[2] = thispt;
break;
}
case 13 : {//MB1-MB3
pt[3] = thispt;
break;
}
case 14 : {//MB1-MB4
pt[5] = thispt;
break;
}
case 23 : {//MB2-MB3
pt[4] = thispt;
break;
}
case 24 : {//MB2-MB4
pt[6] = thispt;
break;
}
case 34 : {//MB3-MB4
pt[7] = thispt;
break;
}
default: {
break;
}
}
}
}
}
| void MuonDTSeedFromRecHits::computePtWithVtx | ( | double * | pt, |
| double * | spt | ||
| ) | const [private] |
Definition at line 171 of file MuonDTSeedFromRecHits.cc.
References bestEta(), DTChamberId, MuonSeedPtExtractor::pT_extract(), relativeConstraints::station, MuonSeedFromRecHits::thePtExtractor, and MuonSeedFromRecHits::theRhits.
Referenced by seed().
{
float eta0 = bestEta();
for (MuonRecHitContainer::const_iterator iter=theRhits.begin(); iter!=theRhits.end(); iter++ ) {
//+vvp !:
float eta1= (*iter)->globalPosition().eta();
if ( fabs (eta1-eta0) > .2 ) continue; // !!! +vvp
// assign Pt from MB1 & vtx
unsigned int stat = DTChamberId((**iter).geographicalId()).station();
if(stat > 2) continue;
double thispt = thePtExtractor->pT_extract(*iter, *iter)[0];
float ptmax = 2000.;
if(thispt > ptmax) thispt = ptmax;
if(thispt < -ptmax) thispt = -ptmax;
if( stat==1 ) {
pt[0] = thispt;
}
// assign Pt from MB2 & vtx
if( stat==2 ) {
pt[1] = thispt;
}
}
}
| TrajectorySeed MuonDTSeedFromRecHits::seed | ( | ) | const [virtual] |
compute pts with vertex constraint
now w/o vertex constrain
now combine all pt estimate
Definition at line 35 of file MuonDTSeedFromRecHits.cc.
References bestBarrelHit(), computeBestPt(), computePtWithoutVtx(), computePtWithVtx(), MuonSeedFromRecHits::createSeed(), i, prof2calltree::last, LogTrace, metname, mathSSE::sqrt(), and MuonSeedFromRecHits::theRhits.
Referenced by MuonSeedFinder::seeds().
{
double pt[16] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. };
// these weights are supposed to be 1/sigma^2(pt), but they're so small.
// Instead of the 0.2-0.5 GeV here, we'll add something extra in quadrature later
double spt[8] = { 1/0.048 , 1/0.075 , 1/0.226 , 1/0.132 , 1/0.106 , 1/0.175 , 1/0.125 , 1/0.185 };
const std::string metname = "Muon|RecoMuon|MuonDTSeedFromRecHits";
computePtWithVtx(pt, spt);
computePtWithoutVtx(pt, spt);
// some dump...
LogTrace(metname) << " Pt MB1 @vtx: " << pt[0] << " w: " << spt[0] << "\n"
<< " Pt MB2 @vtx: " << pt[1] << " w: " << spt[1]<< endl ;
LogTrace(metname) << " Pt MB2-MB1 " << pt[2] << " w: " << spt[2]<< "\n"
<< " Pt MB3-MB1 " << pt[3] << " w: " << spt[3]<< "\n"
<< " Pt MB3-MB2 " << pt[4] << " w: " << spt[4]<< "\n"
<< " Pt MB4-MB1 " << pt[5] << " w: " << spt[5]<< "\n"
<< " Pt MB4-MB2 " << pt[6] << " w: " << spt[6]<< "\n"
<< " Pt MB4-MB3 " << pt[7] << " w: " << spt[7]<< endl ;
float ptmean=0.;
float sptmean=0.;
//@@ Use Shih-Chuan's
computeBestPt(pt, spt, ptmean, sptmean);
// add an extra term to the error in quadrature, 30% of pT per point
int npoints = 0;
for(int i = 0; i < 8; ++i)
if(pt[i] != 0) ++npoints;
if(npoints != 0) {
sptmean = sqrt(sptmean*sptmean + 0.09*ptmean*ptmean/npoints);
}
LogTrace(metname) << " Seed Pt: " << ptmean << " +/- " << sptmean << endl;
// take the best candidate
ConstMuonRecHitPointer last = bestBarrelHit(theRhits);
return createSeed(ptmean, sptmean,last);
}
1.7.3