MuonDTSeedFromRecHits Class Reference

#include <MuonDTSeedFromRecHits.h>

Inheritance diagram for MuonDTSeedFromRecHits:

Public Member Functions
ConstMuonRecHitPointer	bestBarrelHit (const MuonRecHitContainer &barrelHits) const
	MuonDTSeedFromRecHits ()
virtual TrajectorySeed	seed () const
Public Member Functions inherited from MuonSeedFromRecHits
void	add (MuonTransientTrackingRecHit::MuonRecHitPointer hit)
void	clear ()
TrajectorySeed	createSeed (float ptmean, float sptmean, MuonTransientTrackingRecHit::ConstMuonRecHitPointer last) const
MuonTransientTrackingRecHit::ConstMuonRecHitPointer	firstRecHit () const
	MuonSeedFromRecHits ()
unsigned int	nrhit () const
void	setBField (const MagneticField *field)
void	setPtExtractor (const MuonSeedPtExtractor *extractor)
virtual	~MuonSeedFromRecHits ()

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

Additional Inherited Members
Protected Types inherited from MuonSeedFromRecHits
typedef MuonTransientTrackingRecHit::ConstMuonRecHitPointer	ConstMuonRecHitPointer
typedef MuonTransientTrackingRecHit::MuonRecHitContainer	MuonRecHitContainer
typedef MuonTransientTrackingRecHit::MuonRecHitPointer	MuonRecHitPointer
Protected Attributes inherited from MuonSeedFromRecHits
const MagneticField *	theField
const MuonSeedPtExtractor *	thePtExtractor
MuonTransientTrackingRecHit::MuonRecHitContainer	theRhits

Detailed Description

Author

A. Vitelli - INFN Torino

porting R.Bellan - INFN Torino

Generate a seed starting from a list of RecHits make use of TrajectorySeed from CommonDet

Definition at line 20 of file MuonDTSeedFromRecHits.h.

Constructor & Destructor Documentation

MuonDTSeedFromRecHits::MuonDTSeedFromRecHits

(

)

Definition at line 29 of file MuonDTSeedFromRecHits.cc.

: MuonSeedFromRecHits()
{
}

Member Function Documentation

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);
}