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 19 of file MuonDTSeedFromRecHits.h.

Constructor & Destructor Documentation

MuonDTSeedFromRecHits::MuonDTSeedFromRecHits

(

)

Definition at line 25 of file MuonDTSeedFromRecHits.cc.

: MuonSeedFromRecHits() {}

Member Function Documentation

MuonDTSeedFromRecHits::ConstMuonRecHitPointer MuonDTSeedFromRecHits::bestBarrelHit

(

const MuonRecHitContainer &

barrelHits

)

const

Definition at line 73 of file MuonDTSeedFromRecHits.cc.

Referenced by MuonOverlapSeedFromRecHits::bestHit(), and seed().

                                                 {
  int alt_npt = 0;
  int best_npt = 0;
  int cur_npt = 0;
  MuonRecHitPointer best = nullptr;
  MuonRecHitPointer alter = nullptr;

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

References HLT_2018_cff::eta1, HLT_2018_cff::eta2, mps_fire::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 243 of file MuonDTSeedFromRecHits.cc.

References computeMean(), mps_fire::i, LogTrace, metname, mathSSE::sqrt(), and AlCaHLTBitMon_QueryRunRegistry::string.

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

References mps_fire::i, dqmiodumpmetadata::n, and mathSSE::sqrt().

Referenced by computeBestPt().

                                                                                                             {
  int n = 0;
  double ptTmp[8];
  double wtTmp[8];
  assert(sz <= 8);

  for (int i = 0; i < 8; ++i) {
    ptTmp[i] = 0.;
    wtTmp[i] = 0;
    if (i < sz && 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 187 of file MuonDTSeedFromRecHits.cc.

References bestEta(), HLT_2018_cff::eta1, HLT_2018_cff::eta2, MuonSeedPtExtractor::pT_extract(), relativeConstraints::station, MuonSeedFromRecHits::thePtExtractor, MuonSeedFromRecHits::theRhits, and createJobs::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 155 of file MuonDTSeedFromRecHits.cc.

References bestEta(), HLT_2018_cff::eta1, MuonSeedPtExtractor::pT_extract(), muonTiming_cfi::ptmax, 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 27 of file MuonDTSeedFromRecHits.cc.

References bestBarrelHit(), computeBestPt(), computePtWithoutVtx(), computePtWithVtx(), MuonSeedFromRecHits::createSeed(), mps_fire::i, dqmdumpme::last, LogTrace, metname, npoints, DiDispStaMuonMonitor_cfi::pt, mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, 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);
}