CMSSW Reference Manual

MuonDTSeedFromRecHits::MuonDTSeedFromRecHits ( )

Definition at line 35 of file MuonDTSeedFromRecHits.cc.

00036 : MuonSeedFromRecHits()
00037 {
00038 }

MuonDTSeedFromRecHits::ConstMuonRecHitPointer MuonDTSeedFromRecHits::best_cand ( ) const [private]

Definition at line 88 of file MuonDTSeedFromRecHits.cc.

References iter, and MuonSeedFromRecHits::theRhits.

Referenced by seed().

00088                                        {
00089 
00090   int alt_npt = 0;
00091   int best_npt = 0;
00092   int cur_npt = 0;
00093   MuonRecHitPointer best = 0;
00094   MuonRecHitPointer alter=0;
00095 
00096   for (MuonRecHitContainer::const_iterator iter=theRhits.begin();
00097        iter!=theRhits.end(); iter++ ) {
00098 
00099     bool hasZed = ((*iter)->projectionMatrix()[1][1]==1);
00100 
00101     cur_npt = 0;
00102     vector<TrackingRecHit*> slrhs=(*iter)->recHits();
00103     for (vector<TrackingRecHit*>::const_iterator slrh=slrhs.begin(); slrh!=slrhs.end(); ++slrh ) {
00104       cur_npt+=(*slrh)->recHits().size();
00105     }
00106     float radius1 = (*iter)->det()->position().perp();
00107 
00108     if (hasZed) {
00109       if ( cur_npt > best_npt ) {
00110         best = (*iter);
00111         best_npt = cur_npt;
00112       }
00113       else if ( best && cur_npt == best_npt ) {
00114         float radius2 = best->det()->position().perp();
00115        if (radius1 < radius2) {
00116           best = (*iter);
00117           best_npt = cur_npt;
00118        }
00119       }
00120     }
00121 
00122     if ( cur_npt > alt_npt ) {
00123       alter = (*iter);
00124       alt_npt = cur_npt;
00125     }
00126     else if ( alter && cur_npt == alt_npt ) {
00127       float radius2 = alter->det()->position().perp();
00128       if (radius1 < radius2) {
00129         alter = (*iter);
00130         alt_npt = cur_npt;
00131       }
00132     }
00133   }
00134 
00135   if ( !best ) best = alter;
00136 
00137   return best;
00138 }

float MuonDTSeedFromRecHits::bestEta ( ) const [private]

Definition at line 141 of file MuonDTSeedFromRecHits.cc.

References iter, HLT_VtxMuL3::result, and MuonSeedFromRecHits::theRhits.

Referenced by computePtWithoutVtx(), and computePtWithVtx().

00141                                            {
00142 
00143   int Maxseg = 0;
00144   float Msdeta = 100.;
00145   float result = (*theRhits.begin())->globalPosition().eta();
00146 
00147   for (MuonRecHitContainer::const_iterator iter=theRhits.begin(); iter!=theRhits.end(); iter++ ) {
00148 
00149     float eta1= (*iter)->globalPosition().eta(); 
00150 
00151     int Nseg = 0;
00152     float sdeta = .0;
00153 
00154     for (MuonRecHitContainer::const_iterator iter2=theRhits.begin();  iter2!=theRhits.end(); iter2++ ) {
00155 
00156       if ( iter2 == iter )  continue;
00157 
00158       float eta2 = (*iter2)->globalPosition().eta(); 
00159 
00160       if ( fabs (eta1-eta2) > .2 ) continue;
00161 
00162       Nseg++;
00163       sdeta += fabs (eta1-eta2); 
00164 
00165       if ( Nseg > Maxseg ||  Nseg == Maxseg && sdeta < Msdeta ) {
00166         Maxseg = Nseg;
00167         Msdeta = sdeta;
00168         result = eta1;
00169       }
00170 
00171     }
00172   }   //  +v.
00173   return result;
00174 }

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 324 of file MuonDTSeedFromRecHits.cc.

References computeMean(), lat::endl(), i, LogTrace, and funct::sqrt().

Referenced by seed().

00327                                                                  {
00328 
00329   const std::string metname = "Muon|RecoMuon|MuonDTSeedFromRecHits";
00330 
00331   int nTotal=8;
00332   int igood=-1;
00333   for (int i=0;i<=7;i++) {
00334     if(pt[i]==0) {
00335       spt[i]=0.;
00336       nTotal--;
00337     } else {
00338       igood=i;
00339     }
00340   }
00341 
00342   if (nTotal==1) {
00344     ptmean=pt[igood];
00345     sptmean=1/sqrt(spt[igood]);
00346   } else if (nTotal==0) {
00347     // No estimate (e.g. only one rechit!)
00348     ptmean=50;
00349     sptmean=30;
00350   } else {
00352     // calculate pt with vertex
00353     float ptvtx=0.0;
00354     float sptvtx=0.0;
00355     computeMean(pt, spt, 2, false, ptvtx, sptvtx);
00356     LogTrace(metname) << " GSL: Pt w vtx :" << ptvtx << "+/-" <<
00357       sptvtx << endl;
00358       
00359     // FIXME: temp hack
00360     if(ptvtx != 0.) {
00361       ptmean = ptvtx;
00362       sptmean = sptvtx;
00363       return;
00364       //
00365     }
00366     // calculate pt w/o vertex
00367     float ptMB=0.0;
00368     float sptMB=0.0;
00369     computeMean(pt+2, spt+2, 6, false, ptMB, sptMB);
00370     LogTrace(metname) << " GSL: Pt w/o vtx :" << ptMB << "+/-" <<
00371         sptMB << endl;
00372 
00373     // decide wheter the muon comes or not from vertex 
00374     float ptpool=0.0;
00375     if((ptvtx+ptMB)!=0.0) ptpool=(ptvtx-ptMB)/(ptvtx+ptMB);
00376     bool fromvtx=true;
00377     if(fabs(ptpool)>0.2) fromvtx=false; 
00378     LogTrace(metname) << "From vtx? " <<fromvtx << " ptpool "<< ptpool << endl;
00379 
00380     // now choose the "right" pt => weighted mean
00381     computeMean(pt, spt, 8, true, ptmean, sptmean);
00382     LogTrace(metname) << " GSL Pt :" << ptmean << "+/-" << sptmean << endl;
00383 
00384   }
00385 }

void MuonDTSeedFromRecHits::computeMean	(	const double *	pt,
		const double *	weights,
		int	sz,
		bool	tossOutlyers,
		float &	ptmean,
		float &	sptmean
	)			const `[private]`

Definition at line 388 of file MuonDTSeedFromRecHits.cc.

References i, n, and funct::sqrt().

Referenced by computeBestPt().

00390 {
00391   int n=0;
00392   double ptTmp[8];
00393   double wtTmp[8];
00394   assert(sz<=8);
00395 
00396   for (int i=0; i<sz; ++i) {
00397     ptTmp[i] = 0.;
00398     wtTmp[i] = 0;
00399     if (pt[i]!=0) {
00400       ptTmp[n]=pt[i];
00401       wtTmp[n]=weights[i];
00402       ++n;
00403     }
00404   }
00405 
00406   if(n != 0) 
00407   {
00408     if (n==1) {
00409       ptmean=ptTmp[0];
00410       sptmean=1/sqrt(wtTmp[0]);
00411     } else {
00412       ptmean = gsl_stats_wmean(wtTmp, 1, ptTmp, 1, n);
00413       sptmean = sqrt( gsl_stats_wvariance_m (wtTmp, 1, ptTmp, 1, n, ptmean) );
00414     }
00415 
00416     if(tossOutlyers)
00417     {
00418       // Recompute mean with a cut at 3 sigma
00419       for ( int nm =0; nm<8; nm++ ){
00420         if ( ptTmp[nm]!=0 && fabs(ptTmp[nm]-ptmean)>3*(sptmean) ) {
00421           wtTmp[nm]=0.;
00422         }
00423       }
00424       ptmean = gsl_stats_wmean(wtTmp, 1, ptTmp, 1, n);
00425       sptmean = sqrt( gsl_stats_wvariance_m (wtTmp, 1, ptTmp, 1, n, ptmean) );
00426     }
00427   }
00428 
00429 }

void MuonDTSeedFromRecHits::computePtWithoutVtx	(	double *	pt,
		double *	spt
	)			const `[private]`

Definition at line 236 of file MuonDTSeedFromRecHits.cc.

References bestEta(), iter, PV3DBase< T, PVType, FrameType >::phi(), st, MuonSeedFromRecHits::theRhits, tmp, and PV3DBase< T, PVType, FrameType >::z().

Referenced by seed().

00236                                                                              {
00237   float eta0 = bestEta();
00238 
00239   for (MuonRecHitContainer::const_iterator iter=theRhits.begin(); 
00240         iter!=theRhits.end(); iter++ ) {
00241     //+vvp !:
00242     float eta1= (*iter)->globalPosition().eta(); 
00243     if ( fabs (eta1-eta0) > .2 ) continue;  //   !!! +vvp
00244 
00245     float radius1= (*iter)->det()->position().perp(); 
00246 
00247     for (MuonRecHitContainer::const_iterator iter2=theRhits.begin(); 
00248           iter2!=iter; iter2++ ) {
00249 
00250       //+vvp !:
00251       float eta2= (*iter2)->globalPosition().eta(); 
00252       if ( fabs (eta2-eta0) > .2 ) continue;  //   !!! +vvp
00253 
00254       float radius2= (*iter2)->det()->position().perp();
00255 
00256       unsigned int stat1(0), stat2(0);
00257 
00258       if ( radius1>450 && radius1<550 ) stat1=2;
00259       if ( radius1>550 && radius1<650 ) stat1=3;
00260       if ( radius1>650 ) stat1=4;
00261       if ( radius1<450 ) stat1=1;
00262 
00263       if ( radius2>450 && radius2<550 ) stat2=2;
00264       if ( radius2>550 && radius2<650 ) stat2=3;
00265       if ( radius2<450 ) stat2=1;
00266       if ( radius2>650 ) stat2=4;
00267 
00268       GlobalVector globalDir1 = (*iter)->globalDirection();
00269       GlobalVector globalDir2 = (*iter2)->globalDirection();
00270       float dphi = -globalDir1.phi()+globalDir2.phi();
00271       // Maybe these aren't necessary with Geom::Phi
00272       if(dphi>M_PI) dphi -= 2*M_PI;
00273       if(dphi<-M_PI) dphi += 2*M_PI;
00274       // assume we're going inward, so + dphi means + charge
00275       int ch = (dphi > 0) ? 1 : -1;
00276 
00277       if ( stat1>stat2) {
00278         ch = -ch;
00279         int tmp = stat1;
00280         stat1 = stat2;
00281         stat2 = tmp;
00282       }
00283       unsigned int st = stat1*10+stat2;
00284 
00285       if ( dphi ) {
00286         dphi = fabs(dphi);
00287         switch (st) {
00288         case  12 : {//MB1-MB2
00289           pt[2]=(12.802+0.38647/dphi)*ch ; 
00290           GlobalPoint pos_iter = (*iter)->globalPosition();
00291           if (  (*iter)->det()->position().perp() <450 ) {
00292             if ( fabs(pos_iter.z())>500. ) {
00293               pt[2]=(12.802+0.16647/dphi)*ch ; 
00294             }
00295           } else {
00296             if ( fabs(pos_iter.z())>600. ) {
00297               pt[2]=(12.802+0.16647/dphi)*ch ; 
00298             }
00299           }
00300           ;break;
00301         }
00302         case  13 : {//MB1-MB3
00303           pt[3]=(.0307+0.99111/dphi)*ch ; ;break;
00304         }
00305         case  14 : {//MB1-MB4
00306           pt[5]=(2.7947+1.1991/dphi)*ch ; ;break;
00307         }
00308         case  23 : {//MB2-MB3
00309           pt[4]=(2.4583+0.69044/dphi)*ch ;;break; 
00310         }
00311         case  24 : {//MB2-MB4
00312           pt[6]=(2.5267+1.1375/dphi)*ch ; ;break; 
00313         }
00314         case  34 : {//MB3-MB4
00315           pt[7]=(4.06444+0.59189/dphi)*ch ; ;break;
00316         }
00317         default: break;
00318         }  
00319       }
00320     }
00321   }
00322 }

void MuonDTSeedFromRecHits::computePtWithVtx	(	double *	pt,
		double *	spt
	)			const `[private]`

Definition at line 177 of file MuonDTSeedFromRecHits.cc.

References a1, a2, funct::abs(), bestEta(), dir, iter, PV3DBase< T, PVType, FrameType >::phi(), radius(), MuonSeedFromRecHits::theRhits, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by seed().

00177                                                                           {
00178 
00179   float eta0 = bestEta();
00180 
00181   for (MuonRecHitContainer::const_iterator iter=theRhits.begin(); iter!=theRhits.end(); iter++ ) {
00182 
00183  //+vvp !:
00184 
00185     float eta1= (*iter)->globalPosition().eta(); 
00186 
00187     if ( fabs (eta1-eta0) > .2 ) continue;  //   !!! +vvp
00188 
00189     // assign Pt from MB1 & vtx   
00190     float radius = (*iter)->det()->position().perp();
00191     unsigned int stat = 0;
00192     if ( radius>450 && radius<550 ) stat=2;
00193     if ( radius<450 ) stat=1;
00194 
00195     if(stat == 0) continue;
00196 
00197     GlobalPoint pos = (*iter)->globalPosition();
00198     GlobalVector dir = (*iter)->globalDirection();
00199 
00200     float dphi = -pos.phi()+dir.phi();
00201     if(dphi>M_PI) dphi -= 2*M_PI;
00202     if(dphi<-M_PI) dphi += 2*M_PI;
00203     int ch = (dphi<0) ? 1 : -1;
00204 
00205     if( stat==1 ) {
00206       pt[0]=(-1.0+1.46/fabs(dphi)) * ch; 
00207       if ( abs(pos.z()) > 500 ) {
00208         // overlap 
00209         float a1 = dir.y()/dir.x(); float a2 = pos.y()/pos.x();
00210         dphi = fabs((a1-a2)/(1+a1*a2));
00211 
00212         pt[0] = fabs(-3.3104+(1.2373/dphi)) * ch;
00213       }
00214     }
00215     // assign Pt from MB2 & vtx
00216     if( stat==2 ) {
00217       pt[1]=(-1.0+0.9598/fabs(dphi))*ch;
00218       if ( abs(pos.z()) > 600 ) {
00219         // overlap 
00220         float a1 = dir.y()/dir.x(); float a2 = pos.y()/pos.x();
00221         dphi = fabs((a1-a2)/(1+a1*a2));
00222 
00223         pt[1] = fabs(10.236+(0.5766/dphi)) * ch;
00224       }
00225     }
00226     float ptmax = 2000.;
00227     if(pt[0] > ptmax) pt[0] = ptmax;
00228     if(pt[0] < -ptmax) pt[0] = -ptmax;
00229     if(pt[1] > ptmax) pt[1] = ptmax;
00230     if(pt[1] < -ptmax) pt[1] = -ptmax;
00231 
00232   }
00233 }

TrajectorySeed MuonDTSeedFromRecHits::seed ( ) const [virtual]

compute pts with vertex constraint

now w/o vertex constrain

now combine all pt estimate

Definition at line 41 of file MuonDTSeedFromRecHits.cc.

References best_cand(), computeBestPt(), computePtWithoutVtx(), computePtWithVtx(), MuonSeedFromRecHits::createSeed(), lat::endl(), i, prof2calltree::last, LogTrace, and funct::sqrt().

Referenced by MuonSeedFinder::seeds().

00041                                                  {
00042   double pt[8] = { 0.0, 0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 ,0.0 };
00043   // these weights are supposed to be 1/sigma^2(pt), but they're so small.
00044   // Instead of the 0.2-0.5 GeV here, we'll add something extra in quadrature later
00045   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 }; 
00046 
00047   const std::string metname = "Muon|RecoMuon|MuonDTSeedFromRecHits";
00048 
00050   computePtWithVtx(pt, spt);
00051 
00053   computePtWithoutVtx(pt, spt);
00054 
00055   // some dump...
00056   LogTrace(metname) << " Pt MB1 @vtx: " << pt[0] << " w: " << spt[0] << "\n" 
00057                     << " Pt MB2 @vtx: " << pt[1] << " w: " << spt[1]<< endl ;
00058   
00059   LogTrace(metname) << " Pt MB2-MB1 " << pt[2] << " w: " << spt[2]<< "\n" 
00060                     << " Pt MB3-MB1 " << pt[3] << " w: " << spt[3]<< "\n" 
00061                     << " Pt MB3-MB2 " << pt[4] << " w: " << spt[4]<< "\n" 
00062                     << " Pt MB4-MB1 " << pt[5] << " w: " << spt[5]<< "\n" 
00063                     << " Pt MB4-MB2 " << pt[6] << " w: " << spt[6]<< "\n" 
00064                     << " Pt MB4-MB3 " << pt[7] << " w: " << spt[7]<< endl  ;
00065   
00067   float ptmean=0.;
00068   float sptmean=0.;
00069   computeBestPt(pt, spt, ptmean, sptmean);
00070   
00071   // add an extra term to the error in quadrature, 30% of pT per point
00072   int npoints = 0;
00073   for(int i = 0; i < 8; ++i)
00074     if(pt[i] != 0) ++npoints;
00075   if(npoints != 0) {
00076     sptmean = sqrt(sptmean*sptmean + 0.09*ptmean*ptmean/npoints);
00077   }
00078 
00079   LogTrace(metname) << " Seed Pt: " << ptmean << " +/- " << sptmean << endl;
00080   
00081   // take the best candidate
00082   ConstMuonRecHitPointer last = best_cand();
00083   return createSeed(ptmean, sptmean,last);
00084 }


Public Member Functions
	MuonDTSeedFromRecHits ()
virtual TrajectorySeed	seed () const
Private Member Functions
MuonTransientTrackingRecHit::ConstMuonRecHitPointer	best_cand () const
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

MuonDTSeedFromRecHits Class Reference

Public Member Functions

Private Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation