CMS 3D CMS Logo

List of all members | Classes | Public Types | Public Member Functions | Private Member Functions | Private Attributes
CSCSegAlgoRU Class Reference

#include <CSCSegAlgoRU.h>

Inheritance diagram for CSCSegAlgoRU:
CSCSegmentAlgorithm

Classes

struct  AlgoState
 

Public Types

typedef std::vector< bool > BoolContainer
 
typedef std::vector< const CSCRecHit2D * > ChamberHitContainer
 
typedef std::vector< const CSCRecHit2D * >::const_iterator ChamberHitContainerCIt
 
typedef std::vector< int > LayerIndex
 
typedef ROOT::Math::SVector< double, 6 > SVector6
 Typedefs. More...
 

Public Member Functions

std::vector< CSCSegmentbuildSegments (const CSCChamber *aChamber, const ChamberHitContainer &rechits) const
 
 CSCSegAlgoRU (const edm::ParameterSet &ps)
 Constructor. More...
 
std::vector< CSCSegmentrun (const CSCChamber *aChamber, const ChamberHitContainer &rechits) override
 
 ~CSCSegAlgoRU () override
 Destructor. More...
 
- Public Member Functions inherited from CSCSegmentAlgorithm
 CSCSegmentAlgorithm (const edm::ParameterSet &)
 Constructor. More...
 
virtual std::vector< CSCSegmentrun (const CSCChamber *chamber, const std::vector< const CSCRecHit2D *> &rechits)=0
 
virtual ~CSCSegmentAlgorithm ()
 Destructor. More...
 

Private Member Functions

bool addHit (AlgoState &aState, const CSCRecHit2D *hit, int layer) const
 Utility functions. More...
 
bool areHitsCloseInGlobalPhi (const AlgoState &aState, const CSCRecHit2D *h1, const CSCRecHit2D *h2) const
 
bool areHitsCloseInR (const AlgoState &aState, const CSCRecHit2D *h1, const CSCRecHit2D *h2) const
 Utility functions. More...
 
void baseline (AlgoState &aState, int n_seg_min) const
 
void compareProtoSegment (AlgoState &aState, const CSCRecHit2D *h, int layer) const
 
float fit_r_phi (const AlgoState &aState, const SVector6 &points, int layer) const
 
float fitX (const AlgoState &aState, SVector6 points, SVector6 errors, int ir, int ir2, float &chi2_str) const
 
void flagHitsAsUsed (const AlgoState &aState, const ChamberHitContainer &rechitsInChamber, BoolContainer &used) const
 
bool hasHitOnLayer (const AlgoState &aState, int layer) const
 
void increaseProtoSegment (AlgoState &aState, const CSCRecHit2D *h, int layer, int chi2_factor) const
 
bool isHitNearSegment (const AlgoState &aState, const CSCRecHit2D *h) const
 
bool isSegmentGood (const AlgoState &aState, const ChamberHitContainer &rechitsInChamber) const
 
float phiAtZ (const AlgoState &aState, float z) const
 
bool replaceHit (AlgoState &aState, const CSCRecHit2D *h, int layer) const
 
void tryAddingHitsToSegment (AlgoState &aState, const ChamberHitContainer &rechitsInChamber, const BoolContainer &used, const LayerIndex &layerIndex, const ChamberHitContainerCIt i1, const ChamberHitContainerCIt i2) const
 
void updateParameters (AlgoState &aState) const
 

Private Attributes

float chi2_str_
 
float chi2Max
 
float chi2Norm_2D_
 
bool debugInfo
 
bool doCollisions
 
float dPhiIntMax
 
float dPhiMax
 
float dRIntMax
 
float dRMax
 
int minLayersApart
 
const std::string myName
 
double theChi2
 
LocalVector theDirection
 
LocalPoint theOrigin
 
float uz
 
float vz
 
float wideSeg
 

Detailed Description

This is the original algorithm for building endcap muon track segments out of the rechit's in a CSCChamber 'RU' = 'RUssia' = Road Usage

A CSCSegment is a RecSegment4D, and is built from CSCRecHit2D objects, each of which is a RecHit2DLocalPos.

This class is used by the CSCSegmentAlgorithm.
Alternative algorithms can be used for the segment building by writing classes like this, and then selecting which one is actually used via the CSCSegmentBuilder.

developed and implemented by Vladimir Palichik Vladi.nosp@m.mir..nosp@m.Paltc.nosp@m.hik@.nosp@m.cern..nosp@m.ch and Nikolay Voytishin nikol.nosp@m.ay.v.nosp@m.oytis.nosp@m.hin@.nosp@m.cern..nosp@m.ch

Definition at line 35 of file CSCSegAlgoRU.h.

Member Typedef Documentation

◆ BoolContainer

typedef std::vector<bool> CSCSegAlgoRU::BoolContainer

Definition at line 59 of file CSCSegAlgoRU.h.

◆ ChamberHitContainer

typedef std::vector<const CSCRecHit2D*> CSCSegAlgoRU::ChamberHitContainer

Definition at line 51 of file CSCSegAlgoRU.h.

◆ ChamberHitContainerCIt

typedef std::vector<const CSCRecHit2D*>::const_iterator CSCSegAlgoRU::ChamberHitContainerCIt

Definition at line 52 of file CSCSegAlgoRU.h.

◆ LayerIndex

typedef std::vector<int> CSCSegAlgoRU::LayerIndex

Definition at line 50 of file CSCSegAlgoRU.h.

◆ SVector6

typedef ROOT::Math::SVector<double, 6> CSCSegAlgoRU::SVector6

Typedefs.

Definition at line 48 of file CSCSegAlgoRU.h.

Constructor & Destructor Documentation

◆ CSCSegAlgoRU()

CSCSegAlgoRU::CSCSegAlgoRU ( const edm::ParameterSet ps)
explicit

Constructor.

Definition at line 22 of file CSCSegAlgoRU.cc.

References chi2_str_, chi2Max, chi2Norm_2D_, doCollisions, dPhiIntMax, dPhiMax, dRIntMax, dRMax, edm::ParameterSet::getParameter(), LogDebug, minLayersApart, myName, and wideSeg.

22  : CSCSegmentAlgorithm(ps), myName("CSCSegAlgoRU") {
23  doCollisions = ps.getParameter<bool>("doCollisions");
24  chi2_str_ = ps.getParameter<double>("chi2_str");
25  chi2Norm_2D_ = ps.getParameter<double>("chi2Norm_2D_");
26  dRMax = ps.getParameter<double>("dRMax");
27  dPhiMax = ps.getParameter<double>("dPhiMax");
28  dRIntMax = ps.getParameter<double>("dRIntMax");
29  dPhiIntMax = ps.getParameter<double>("dPhiIntMax");
30  chi2Max = ps.getParameter<double>("chi2Max");
31  wideSeg = ps.getParameter<double>("wideSeg");
32  minLayersApart = ps.getParameter<int>("minLayersApart");
33 
34  LogDebug("CSC") << myName << " has algorithm cuts set to: \n"
35  << "--------------------------------------------------------------------\n"
36  << "dRMax = " << dRMax << '\n'
37  << "dPhiMax = " << dPhiMax << '\n'
38  << "dRIntMax = " << dRIntMax << '\n'
39  << "dPhiIntMax = " << dPhiIntMax << '\n'
40  << "chi2Max = " << chi2Max << '\n'
41  << "wideSeg = " << wideSeg << '\n'
42  << "minLayersApart = " << minLayersApart << std::endl;
43 
44  //reset the thresholds for non-collision data
45  if (!doCollisions) {
46  dRMax = 2.0;
47  dPhiMax = 2 * dPhiMax;
48  dRIntMax = 2 * dRIntMax;
49  dPhiIntMax = 2 * dPhiIntMax;
51  chi2_str_ = 100;
52  chi2Max = 2 * chi2Max;
53  }
54 }
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
CSCSegmentAlgorithm(const edm::ParameterSet &)
Constructor.
const std::string myName
Definition: CSCSegAlgoRU.h:147
float dPhiIntMax
Definition: CSCSegAlgoRU.h:157
float chi2Norm_2D_
Definition: CSCSegAlgoRU.h:160
#define LogDebug(id)

◆ ~CSCSegAlgoRU()

CSCSegAlgoRU::~CSCSegAlgoRU ( )
inlineoverride

Destructor.

Definition at line 64 of file CSCSegAlgoRU.h.

64 {};

Member Function Documentation

◆ addHit()

bool CSCSegAlgoRU::addHit ( AlgoState aState,
const CSCRecHit2D hit,
int  layer 
) const
private

Utility functions.

Definition at line 582 of file CSCSegAlgoRU.cc.

References ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, nano_mu_digi_cff::layer, CSCSegAlgoRU::AlgoState::proto_segment, and updateParameters().

Referenced by buildSegments(), increaseProtoSegment(), and replaceHit().

582  {
583  // Return true if hit was added successfully
584  // (and then parameters are updated).
585  // Return false if there is already a hit on the same layer, or insert failed.
586  ChamberHitContainer::const_iterator it;
587  for (it = aState.proto_segment.begin(); it != aState.proto_segment.end(); it++)
588  if (((*it)->cscDetId().layer() == layer) && (aHit != (*it)))
589  return false;
590  aState.proto_segment.push_back(aHit);
591  // make a fit
592  updateParameters(aState);
593  return true;
594 }
void updateParameters(AlgoState &aState) const

◆ areHitsCloseInGlobalPhi()

bool CSCSegAlgoRU::areHitsCloseInGlobalPhi ( const AlgoState aState,
const CSCRecHit2D h1,
const CSCRecHit2D h2 
) const
private

Definition at line 427 of file CSCSegAlgoRU.cc.

References CSCSegAlgoRU::AlgoState::aChamber, PV3DBase< T, PVType, FrameType >::barePhi(), CSCRecHit2D::cscDetId(), SiPixelRawToDigiRegional_cfi::deltaPhi, CSCSegAlgoRU::AlgoState::dPhiMax, CSCRecHit2D::errorWithinStrip(), CSCDetId::iChamberType(), CSCDetId::layer(), CSCChamber::layer(), CSCRecHit2D::localPosition(), CSCSegAlgoRU::AlgoState::strip_iadd, GeomDet::toGlobal(), and funct::true.

Referenced by buildSegments().

429  {
430  float strip_width[10] = {0.003878509,
431  0.002958185,
432  0.002327105,
433  0.00152552,
434  0.00465421,
435  0.002327105,
436  0.00465421,
437  0.002327105,
438  0.00465421,
439  0.002327105}; //in rad
440  const CSCLayer* l1 = aState.aChamber->layer(h1->cscDetId().layer());
441  GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
442  const CSCLayer* l2 = aState.aChamber->layer(h2->cscDetId().layer());
443  GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
444  float err_stpos_h1 = h1->errorWithinStrip();
445  float err_stpos_h2 = h2->errorWithinStrip();
446  CSCDetId id = h1->cscDetId();
447  int iStn = id.iChamberType() - 1;
448  float dphi_incr = 0;
449  if (err_stpos_h1 > 0.25 * strip_width[iStn] || err_stpos_h2 > 0.25 * strip_width[iStn])
450  dphi_incr = 0.5 * strip_width[iStn];
451  float dphi12 = deltaPhi(gp1.barePhi(), gp2.barePhi());
452  return (fabs(dphi12) < (aState.dPhiMax * aState.strip_iadd + dphi_incr)) ? true : false; // +v
453 }
CSCDetId cscDetId() const
Definition: CSCRecHit2D.h:58
int layer() const
Definition: CSCDetId.h:56
float errorWithinStrip() const
The uncertainty of the estimated position within the strip.
Definition: CSCRecHit2D.h:85
unsigned short iChamberType() const
Definition: CSCDetId.h:96
T barePhi() const
Definition: PV3DBase.h:65
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
LocalPoint localPosition() const override
Definition: CSCRecHit2D.h:56

◆ areHitsCloseInR()

bool CSCSegAlgoRU::areHitsCloseInR ( const AlgoState aState,
const CSCRecHit2D h1,
const CSCRecHit2D h2 
) const
private

Utility functions.

Definition at line 389 of file CSCSegAlgoRU.cc.

References CSCSegAlgoRU::AlgoState::aChamber, CSCRecHit2D::cscDetId(), CSCSegAlgoRU::AlgoState::doCollisions, CSCSegAlgoRU::AlgoState::dRMax, CSCRecHit2D::hitWire(), CSCDetId::iChamberType(), CSCDetId::layer(), CSCChamber::layer(), CSCRecHit2D::localPosition(), PV3DBase< T, PVType, FrameType >::perp(), CSCSegAlgoRU::AlgoState::strip_iadd, GeomDet::toGlobal(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by buildSegments().

389  {
390  float maxWG_width[10] = {0, 0, 4.1, 5.69, 2.49, 5.06, 2.49, 5.06, 1.87, 5.06};
391  CSCDetId id = h1->cscDetId();
392  int iStn = id.iChamberType() - 1;
393  //find maxWG_width for ME11 (tilt = 29deg)
394  int wg_num = h2->hitWire();
395  if (iStn == 0 || iStn == 1) {
396  if (wg_num == 1) {
397  maxWG_width[0] = 9.25;
398  maxWG_width[1] = 9.25;
399  }
400  if (wg_num > 1 && wg_num < 48) {
401  maxWG_width[0] = 3.14;
402  maxWG_width[1] = 3.14;
403  }
404  if (wg_num == 48) {
405  maxWG_width[0] = 10.75;
406  maxWG_width[1] = 10.75;
407  }
408  }
409  const CSCLayer* l1 = aState.aChamber->layer(h1->cscDetId().layer());
410  GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
411  const CSCLayer* l2 = aState.aChamber->layer(h2->cscDetId().layer());
412  GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
413  //find z to understand the direction
414  float h1z = gp1.z();
415  float h2z = gp2.z();
416  //switch off the IP check for non collisions case
417  if (!aState.doCollisions) {
418  h1z = 1;
419  h2z = 1;
420  }
421  return (gp2.perp() > ((gp1.perp() - aState.dRMax * aState.strip_iadd * maxWG_width[iStn]) * h2z) / h1z &&
422  gp2.perp() < ((gp1.perp() + aState.dRMax * aState.strip_iadd * maxWG_width[iStn]) * h2z) / h1z)
423  ? true
424  : false;
425 }
T perp() const
Definition: PV3DBase.h:69
T z() const
Definition: PV3DBase.h:61
CSCDetId cscDetId() const
Definition: CSCRecHit2D.h:58
int layer() const
Definition: CSCDetId.h:56
unsigned short iChamberType() const
Definition: CSCDetId.h:96
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
LocalPoint localPosition() const override
Definition: CSCRecHit2D.h:56
short int hitWire() const
L1A.
Definition: CSCRecHit2D.h:68

◆ baseline()

void CSCSegAlgoRU::baseline ( AlgoState aState,
int  n_seg_min 
) const
private

(nhits-2)

(nhits-3)

Definition at line 623 of file CSCSegAlgoRU.cc.

References edmScanValgrind::buffer, CSCRecHit2D::channels(), HLT_2024v10_cff::chi2_str, CSCSegAlgoRU::AlgoState::chi2_str_, CSCSegAlgoRU::AlgoState::chi2D_iadd, CSCRecHit2D::cscDetId(), fitX(), mps_fire::i, testProducerWithPsetDescEmpty_cfi::i1, testProducerWithPsetDescEmpty_cfi::i2, dqmdumpme::k, kLayer(), CSCDetId::layer(), SiStripPI::min, TrackingDataMCValidation_Standalone_cff::nhits, CSCRecHit2D::nStrips(), CSCSegAlgoRU::AlgoState::proto_segment, CSCDetId::ring(), slope, CSCDetId::station(), z, and testProducerWithPsetDescEmpty_cfi::z2.

Referenced by buildSegments().

623  {
624  int nhits = aState.proto_segment.size();
625  //initialise vectors for strip position and error within strip
626  SVector6 sp;
627  SVector6 se;
628  unsigned int init_size = aState.proto_segment.size();
630  buffer.clear();
631  buffer.reserve(init_size);
632  while (buffer.size() < init_size) {
633  ChamberHitContainer::iterator min;
634  int min_layer = 10;
635  for (ChamberHitContainer::iterator k = aState.proto_segment.begin(); k != aState.proto_segment.end(); k++) {
636  const CSCRecHit2D* iRHk = *k;
637  CSCDetId idRHk = iRHk->cscDetId();
638  int kLayer = idRHk.layer();
639  if (kLayer < min_layer) {
640  min_layer = kLayer;
641  min = k;
642  }
643  }
644  buffer.push_back(*min);
645  aState.proto_segment.erase(min);
646  } //while
647 
648  aState.proto_segment.clear();
649  for (ChamberHitContainer::const_iterator cand = buffer.begin(); cand != buffer.end(); cand++) {
650  aState.proto_segment.push_back(*cand);
651  }
652 
653  for (ChamberHitContainer::const_iterator iRH = aState.proto_segment.begin(); iRH != aState.proto_segment.end();
654  iRH++) {
655  const CSCRecHit2D* iRHp = *iRH;
656  CSCDetId idRH = iRHp->cscDetId();
657  int kRing = idRH.ring();
658  int kStation = idRH.station();
659  int kLayer = idRH.layer();
660  // Find the strip containing this hit
661  int centerid = iRHp->nStrips() / 2;
662  int centerStrip = iRHp->channels(centerid);
663  float stpos = (*iRHp).positionWithinStrip();
664  se(kLayer - 1) = (*iRHp).errorWithinStrip();
665  // Take into account half-strip staggering of layers (ME1/1 has no staggering)
666  if (kStation == 1 && (kRing == 1 || kRing == 4))
667  sp(kLayer - 1) = stpos + centerStrip;
668  else {
669  if (kLayer == 1 || kLayer == 3 || kLayer == 5)
670  sp(kLayer - 1) = stpos + centerStrip;
671  if (kLayer == 2 || kLayer == 4 || kLayer == 6)
672  sp(kLayer - 1) = stpos - 0.5 + centerStrip;
673  }
674  }
675  float chi2_str;
676  fitX(aState, sp, se, -1, -1, chi2_str);
677 
678  //-----------------------------------------------------
679  // Optimal point rejection method
680  //-----------------------------------------------------
681  float minSum = 1000;
682  int iworst = -1;
683  int bad_layer = -1;
684  ChamberHitContainer::const_iterator rh_to_be_deleted_1;
685  ChamberHitContainer::const_iterator rh_to_be_deleted_2;
686  if ((chi2_str) > aState.chi2_str_ * aState.chi2D_iadd) {
687  for (ChamberHitContainer::const_iterator i1 = aState.proto_segment.begin(); i1 != aState.proto_segment.end();
688  ++i1) {
689  const CSCRecHit2D* i1_1 = *i1;
690  CSCDetId idRH1 = i1_1->cscDetId();
691  int z1 = idRH1.layer();
692  for (ChamberHitContainer::const_iterator i2 = i1 + 1; i2 != aState.proto_segment.end(); ++i2) {
693  const CSCRecHit2D* i2_1 = *i2;
694  CSCDetId idRH2 = i2_1->cscDetId();
695  int z2 = idRH2.layer();
696  int irej = 0;
697  for (ChamberHitContainer::const_iterator ir = aState.proto_segment.begin(); ir != aState.proto_segment.end();
698  ++ir) {
699  ++irej;
700  if (ir == i1 || ir == i2)
701  continue;
702  float dsum = 0;
703  const CSCRecHit2D* ir_1 = *ir;
704  CSCDetId idRH = ir_1->cscDetId();
705  int worst_layer = idRH.layer();
706  for (ChamberHitContainer::const_iterator i = aState.proto_segment.begin(); i != aState.proto_segment.end();
707  ++i) {
708  const CSCRecHit2D* i_1 = *i;
709  if (i == i1 || i == i2 || i == ir)
710  continue;
711  float slope = (sp(z2 - 1) - sp(z1 - 1)) / (z2 - z1);
712  float intersept = sp(z1 - 1) - slope * z1;
713  CSCDetId idRH = i_1->cscDetId();
714  int z = idRH.layer();
715  float di = fabs(sp(z - 1) - intersept - slope * z);
716  dsum = dsum + di;
717  } //i
718  if (dsum < minSum) {
719  minSum = dsum;
720  bad_layer = worst_layer;
721  iworst = irej;
722  rh_to_be_deleted_1 = ir;
723  }
724  } //ir
725  } //i2
726  } //i1
727  fitX(aState, sp, se, bad_layer, -1, chi2_str);
728  } //if chi2prob<1.0e-4
729 
730  //find worst from n-1 hits
731  int iworst2 = -1;
732  int bad_layer2 = -1;
733  if (iworst > -1 && (nhits - 1) > n_seg_min && (chi2_str) > aState.chi2_str_ * aState.chi2D_iadd) {
734  iworst = -1;
735  float minSum = 1000;
736  for (ChamberHitContainer::const_iterator i1 = aState.proto_segment.begin(); i1 != aState.proto_segment.end();
737  ++i1) {
738  const CSCRecHit2D* i1_1 = *i1;
739  CSCDetId idRH1 = i1_1->cscDetId();
740  int z1 = idRH1.layer();
741  for (ChamberHitContainer::const_iterator i2 = i1 + 1; i2 != aState.proto_segment.end(); ++i2) {
742  const CSCRecHit2D* i2_1 = *i2;
743  CSCDetId idRH2 = i2_1->cscDetId();
744  int z2 = idRH2.layer();
745  int irej = 0;
746  for (ChamberHitContainer::const_iterator ir = aState.proto_segment.begin(); ir != aState.proto_segment.end();
747  ++ir) {
748  ++irej;
749  int irej2 = 0;
750  if (ir == i1 || ir == i2)
751  continue;
752  const CSCRecHit2D* ir_1 = *ir;
753  CSCDetId idRH = ir_1->cscDetId();
754  int worst_layer = idRH.layer();
755  for (ChamberHitContainer::const_iterator ir2 = aState.proto_segment.begin();
756  ir2 != aState.proto_segment.end();
757  ++ir2) {
758  ++irej2;
759  if (ir2 == i1 || ir2 == i2 || ir2 == ir)
760  continue;
761  float dsum = 0;
762  const CSCRecHit2D* ir2_1 = *ir2;
763  CSCDetId idRH = ir2_1->cscDetId();
764  int worst_layer2 = idRH.layer();
765  for (ChamberHitContainer::const_iterator i = aState.proto_segment.begin(); i != aState.proto_segment.end();
766  ++i) {
767  const CSCRecHit2D* i_1 = *i;
768  if (i == i1 || i == i2 || i == ir || i == ir2)
769  continue;
770  float slope = (sp(z2 - 1) - sp(z1 - 1)) / (z2 - z1);
771  float intersept = sp(z1 - 1) - slope * z1;
772  CSCDetId idRH = i_1->cscDetId();
773  int z = idRH.layer();
774  float di = fabs(sp(z - 1) - intersept - slope * z);
775  dsum = dsum + di;
776  } //i
777  if (dsum < minSum) {
778  minSum = dsum;
779  iworst2 = irej2;
780  iworst = irej;
781  bad_layer = worst_layer;
782  bad_layer2 = worst_layer2;
783  rh_to_be_deleted_1 = ir;
784  rh_to_be_deleted_2 = ir2;
785  }
786  } //ir2
787  } //ir
788  } //i2
789  } //i1
790  fitX(aState, sp, se, bad_layer, bad_layer2, chi2_str);
791  } //if prob(n-1)<e-4
792 
793  //----------------------------------
794  //erase bad_hits
795  //----------------------------------
796  if (iworst2 - 1 >= 0 && iworst2 <= int(aState.proto_segment.size())) {
797  aState.proto_segment.erase(rh_to_be_deleted_2);
798  }
799  if (iworst - 1 >= 0 && iworst <= int(aState.proto_segment.size())) {
800  aState.proto_segment.erase(rh_to_be_deleted_1);
801  }
802 }
unsigned int nStrips() const
Definition: CSCRecHit2D.h:62
static const double slope[3]
CSCDetId cscDetId() const
Definition: CSCRecHit2D.h:58
int layer() const
Definition: CSCDetId.h:56
float fitX(const AlgoState &aState, SVector6 points, SVector6 errors, int ir, int ir2, float &chi2_str) const
static const std::string kLayer("layer")
std::vector< const CSCRecHit2D * > ChamberHitContainer
Definition: CSCSegAlgoRU.h:51
ROOT::Math::SVector< double, 6 > SVector6
Typedefs.
Definition: CSCSegAlgoRU.h:48
int station() const
Definition: CSCDetId.h:79
int channels(unsigned int i) const
Extracting strip channel numbers comprising the rechit - low.
Definition: CSCRecHit2D.h:61
int ring() const
Definition: CSCDetId.h:68

◆ buildSegments()

std::vector< CSCSegment > CSCSegAlgoRU::buildSegments ( const CSCChamber aChamber,
const ChamberHitContainer rechits 
) const

Build track segments in this chamber (this is where the actual segment-building algorithm hides.)

Definition at line 56 of file CSCSegAlgoRU.cc.

References funct::abs(), CSCSegAlgoRU::AlgoState::aChamber, addHit(), areHitsCloseInGlobalPhi(), areHitsCloseInR(), baseline(), CSCRecHit2D::channels(), chi2_str_, chi2Max, chi2Norm_2D_, CSCRecHit2D::cscDetId(), doCollisions, dPhiIntMax, dPhiMax, dRIntMax, dRMax, flagHitsAsUsed(), CSCRecHit2D::hitWire(), mps_fire::i, testProducerWithPsetDescEmpty_cfi::i1, testProducerWithPsetDescEmpty_cfi::i2, cuy::ib, isSegmentGood(), ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, dqmiolumiharvest::j, dqmdumpme::k, CSCDetId::layer(), CSCChamber::layer(), CSCRecHit2D::nStrips(), GeomDet::position(), groupFilesInBlocks::reverse, groupFilesInBlocks::temp, tryAddingHitsToSegment(), updateParameters(), wideSeg, and PV3DBase< T, PVType, FrameType >::z().

Referenced by run().

57  {
58  ChamberHitContainer rechits = urechits;
59  LayerIndex layerIndex(rechits.size());
60  int recHits_per_layer[6] = {0, 0, 0, 0, 0, 0};
61  //skip events with high multiplicity of hits
62  if (rechits.size() > 150) {
63  return std::vector<CSCSegment>();
64  }
65  int iadd = 0;
66  for (unsigned int i = 0; i < rechits.size(); i++) {
67  recHits_per_layer[rechits[i]->cscDetId().layer() - 1]++; //count rh per chamber
68  layerIndex[i] = rechits[i]->cscDetId().layer();
69  }
70  double z1 = aChamber->layer(1)->position().z();
71  double z6 = aChamber->layer(6)->position().z();
72  if (std::abs(z1) > std::abs(z6)) {
73  reverse(layerIndex.begin(), layerIndex.end());
74  reverse(rechits.begin(), rechits.end());
75  }
76  if (rechits.size() < 2) {
77  return std::vector<CSCSegment>();
78  }
79  // We have at least 2 hits. We intend to try all possible pairs of hits to start
80  // segment building. 'All possible' means each hit lies on different layers in the chamber.
81  // after all same size segs are build we get rid of the overcrossed segments using the chi2 criteria
82  // the hits from the segs that are left are marked as used and are not added to segs in future iterations
83  // the hits from 3p segs are marked as used separately in order to try to assamble them in longer segments
84  // in case there is a second pass
85  // Choose first hit (as close to IP as possible) h1 and second hit
86  // (as far from IP as possible) h2 To do this we iterate over hits
87  // in the chamber by layer - pick two layers. Then we
88  // iterate over hits within each of these layers and pick h1 and h2
89  // these. If they are 'close enough' we build an empty
90  // segment. Then try adding hits to this segment.
91  // Initialize flags that a given hit has been allocated to a segment
92  BoolContainer used(rechits.size(), false);
93  BoolContainer used3p(rechits.size(), false);
94  // This is going to point to fits to hits, and its content will be used to create a CSCSegment
95  AlgoState aState;
96  aState.aChamber = aChamber;
97  aState.doCollisions = doCollisions;
98  aState.dRMax = dRMax;
99  aState.dPhiMax = dPhiMax;
100  aState.dRIntMax = dRIntMax;
101  aState.dPhiIntMax = dPhiIntMax;
102  aState.chi2Norm_2D_ = chi2Norm_2D_;
103  aState.chi2_str_ = chi2_str_;
104  aState.chi2Max = chi2Max;
105 
106  // Define buffer for segments we build
107  std::vector<CSCSegment> segments;
108  ChamberHitContainerCIt ib = rechits.begin();
109  ChamberHitContainerCIt ie = rechits.end();
110  // Possibly allow 3 passes, second widening scale factor for cuts, third for segments from displaced vertices
111  aState.windowScale = 1.; // scale factor for cuts
112  bool search_disp = false;
113  aState.strip_iadd = 1;
114  aState.chi2D_iadd = 1;
115  int npass = (wideSeg > 1.) ? 3 : 2;
116  for (int ipass = 0; ipass < npass; ++ipass) {
117  if (aState.windowScale > 1.) {
118  iadd = 1;
119  aState.strip_iadd = 4;
120  aState.chi2D_iadd = 4;
121  aState.chi2Max = 2 * chi2Max;
122  if (rechits.size() <= 12)
123  iadd = 0; //allow 3 hit segments for low hit multiplicity chambers
124  }
125 
126  int used_rh = 0;
127  for (ChamberHitContainerCIt i1 = ib; i1 != ie; ++i1) {
128  if (used[i1 - ib])
129  used_rh++;
130  }
131 
132  //change the tresholds if it's time to look for displaced mu segments
133  if (aState.doCollisions && search_disp &&
134  int(rechits.size() - used_rh) >
135  2) { //check if there are enough recHits left to build a segment from displaced vertices
136  aState.doCollisions = false;
137  aState.windowScale = 1.; // scale factor for cuts
138  aState.dRMax = 4.0;
139  aState.dPhiMax = 4 * aState.dPhiMax;
140  aState.dRIntMax = 4 * aState.dRIntMax;
141  aState.dPhiIntMax = 4 * aState.dPhiIntMax;
142  aState.chi2Norm_2D_ = 10 * aState.chi2Norm_2D_;
143  aState.chi2_str_ = 200;
144  aState.chi2Max = 4 * aState.chi2Max;
145  } else {
146  search_disp = false; //make sure the flag is off
147  }
148 
149  for (unsigned int n_seg_min = 6u; n_seg_min > 2u + iadd; --n_seg_min) {
150  BoolContainer common_used(rechits.size(), false);
151  std::array<BoolContainer, 120> common_used_it = {};
152  for (unsigned int i = 0; i < common_used_it.size(); i++) {
153  common_used_it[i] = common_used;
154  }
155  ChamberHitContainer best_proto_segment[120];
156  float min_chi[120] = {9999};
157  int common_it = 0;
158  bool first_proto_segment = true;
159  for (ChamberHitContainerCIt i1 = ib; i1 != ie; ++i1) {
160  bool segok = false;
161  //skip if rh is used and the layer tat has big rh multiplicity(>25RHs)
162  if (used[i1 - ib] || recHits_per_layer[int(layerIndex[i1 - ib]) - 1] > 25 ||
163  (n_seg_min == 3 && used3p[i1 - ib]))
164  continue;
165  int layer1 = layerIndex[i1 - ib];
166  const CSCRecHit2D* h1 = *i1;
167  for (ChamberHitContainerCIt i2 = ie - 1; i2 != i1; --i2) {
168  if (used[i2 - ib] || recHits_per_layer[int(layerIndex[i2 - ib]) - 1] > 25 ||
169  (n_seg_min == 3 && used3p[i2 - ib]))
170  continue;
171  int layer2 = layerIndex[i2 - ib];
172  if ((abs(layer2 - layer1) + 1) < int(n_seg_min))
173  break; //decrease n_seg_min
174  const CSCRecHit2D* h2 = *i2;
175  if (areHitsCloseInR(aState, h1, h2) && areHitsCloseInGlobalPhi(aState, h1, h2)) {
176  aState.proto_segment.clear();
177  if (!addHit(aState, h1, layer1))
178  continue;
179  if (!addHit(aState, h2, layer2))
180  continue;
181  // Can only add hits if already have a segment
182  if (aState.sfit)
183  tryAddingHitsToSegment(aState, rechits, used, layerIndex, i1, i2);
184  segok = isSegmentGood(aState, rechits);
185  if (segok) {
186  if (aState.proto_segment.size() > n_seg_min) {
187  baseline(aState, n_seg_min);
188  updateParameters(aState);
189  }
190  if (aState.sfit->chi2() > aState.chi2Norm_2D_ * aState.chi2D_iadd ||
191  aState.proto_segment.size() < n_seg_min)
192  aState.proto_segment.clear();
193  if (!aState.proto_segment.empty()) {
194  updateParameters(aState);
195  //add same-size overcrossed protosegments to the collection
196  if (first_proto_segment) {
197  flagHitsAsUsed(aState, rechits, common_used_it[0]);
198  min_chi[0] = aState.sfit->chi2();
199  best_proto_segment[0] = aState.proto_segment;
200  first_proto_segment = false;
201  } else { //for the rest of found proto_segments
202  common_it++;
203  flagHitsAsUsed(aState, rechits, common_used_it[common_it]);
204  min_chi[common_it] = aState.sfit->chi2();
205  best_proto_segment[common_it] = aState.proto_segment;
206  ChamberHitContainerCIt hi, iu, ik;
207  int iter = common_it;
208  for (iu = ib; iu != ie; ++iu) {
209  for (hi = aState.proto_segment.begin(); hi != aState.proto_segment.end(); ++hi) {
210  if (*hi == *iu) {
211  int merge_nr = -1;
212  for (int k = 0; k < iter + 1; k++) {
213  if (common_used_it[k][iu - ib] == true) {
214  if (merge_nr != -1) {
215  //merge the k and merge_nr collections of flaged hits into the merge_nr collection and unmark the k collection hits
216  for (ik = ib; ik != ie; ++ik) {
217  if (common_used_it[k][ik - ib] == true) {
218  common_used_it[merge_nr][ik - ib] = true;
219  common_used_it[k][ik - ib] = false;
220  }
221  }
222  //change best_protoseg if min_chi_2 is smaller
223  if (min_chi[k] < min_chi[merge_nr]) {
224  min_chi[merge_nr] = min_chi[k];
225  best_proto_segment[merge_nr] = best_proto_segment[k];
226  best_proto_segment[k].clear();
227  min_chi[k] = 9999;
228  }
229  common_it--;
230  } else {
231  merge_nr = k;
232  }
233  } //if(common_used[k][iu-ib] == true)
234  } //for k
235  } //if
236  } //for proto_seg
237  } //for rec_hits
238  } //else
239  } //proto seg not empty
240  }
241  } // h1 & h2 close
242  if (segok)
243  break;
244  } // i2
245  } // i1
246 
247  //add the reconstructed segments
248  for (int j = 0; j < common_it + 1; j++) {
249  aState.proto_segment = best_proto_segment[j];
250  best_proto_segment[j].clear();
251  //SKIP empty proto-segments
252  if (aState.proto_segment.empty())
253  continue;
254  updateParameters(aState);
255  // Create an actual CSCSegment - retrieve all info from the fit
256  CSCSegment temp(aState.sfit->hits(),
257  aState.sfit->intercept(),
258  aState.sfit->localdir(),
259  aState.sfit->covarianceMatrix(),
260  aState.sfit->chi2());
261  aState.sfit = nullptr;
262  segments.push_back(temp);
263  //if the segment has 3 hits flag them as used in a particular way
264  if (aState.proto_segment.size() == 3) {
265  flagHitsAsUsed(aState, rechits, used3p);
266  } else {
267  flagHitsAsUsed(aState, rechits, used);
268  }
269  aState.proto_segment.clear();
270  }
271  } //for n_seg_min
272 
273  if (search_disp) {
274  //reset params and flags for the next chamber
275  search_disp = false;
276  aState.doCollisions = true;
277  aState.dRMax = 2.0;
278  aState.chi2_str_ = 100;
279  aState.dPhiMax = 0.25 * aState.dPhiMax;
280  aState.dRIntMax = 0.25 * aState.dRIntMax;
281  aState.dPhiIntMax = 0.25 * aState.dPhiIntMax;
282  aState.chi2Norm_2D_ = 0.1 * aState.chi2Norm_2D_;
283  aState.chi2Max = 0.25 * aState.chi2Max;
284  }
285 
286  std::vector<CSCSegment>::iterator it = segments.begin();
287  bool good_segs = false;
288  while (it != segments.end()) {
289  if ((*it).nRecHits() > 3) {
290  good_segs = true;
291  break;
292  }
293  ++it;
294  }
295  if (good_segs &&
296  aState.doCollisions) { // only change window if not enough good segments were found (bool can be changed to int if a >0 number of good segs is required)
297  search_disp = true;
298  continue; //proceed to search the segs from displaced vertices
299  }
300 
301  // Increase cut windows by factor of wideSeg only for collisions
302  if (!aState.doCollisions && !search_disp)
303  break;
304  aState.windowScale = wideSeg;
305  } // ipass
306 
307  //get rid of enchansed 3p segments
308  std::vector<CSCSegment>::iterator it = segments.begin();
309  while (it != segments.end()) {
310  if ((*it).nRecHits() == 3) {
311  bool found_common = false;
312  const std::vector<CSCRecHit2D>& theseRH = (*it).specificRecHits();
313  for (ChamberHitContainerCIt i1 = ib; i1 != ie; ++i1) {
314  if (used[i1 - ib] && used3p[i1 - ib]) {
315  const CSCRecHit2D* sh1 = *i1;
316  CSCDetId id = sh1->cscDetId();
317  int sh1layer = id.layer();
318  int RH_centerid = sh1->nStrips() / 2;
319  int RH_centerStrip = sh1->channels(RH_centerid);
320  int RH_wg = sh1->hitWire();
321  std::vector<CSCRecHit2D>::const_iterator sh;
322  for (sh = theseRH.begin(); sh != theseRH.end(); ++sh) {
323  CSCDetId idRH = sh->cscDetId();
324  //find segment hit coord
325  int shlayer = idRH.layer();
326  int SegRH_centerid = sh->nStrips() / 2;
327  int SegRH_centerStrip = sh->channels(SegRH_centerid);
328  int SegRH_wg = sh->hitWire();
329  if (sh1layer == shlayer && SegRH_centerStrip == RH_centerStrip && SegRH_wg == RH_wg) {
330  //remove the enchansed 3p segment
331  segments.erase(it, (it + 1));
332  found_common = true;
333  break;
334  }
335  } //theserh
336  }
337  if (found_common)
338  break; //current seg has already been erased
339  } //camber hits
340  if (!found_common)
341  ++it;
342  } //its a 3p seg
343  else {
344  ++it; //go to the next seg
345  }
346  } //while
347  // Give the segments to the CSCChamber
348  return segments;
349 } //build segments
unsigned int nStrips() const
Definition: CSCRecHit2D.h:62
const CSCLayer * layer(CSCDetId id) const
Return the layer corresponding to the given id.
Definition: CSCChamber.cc:30
std::vector< const CSCRecHit2D * >::const_iterator ChamberHitContainerCIt
Definition: CSCSegAlgoRU.h:52
void baseline(AlgoState &aState, int n_seg_min) const
T z() const
Definition: PV3DBase.h:61
bool areHitsCloseInR(const AlgoState &aState, const CSCRecHit2D *h1, const CSCRecHit2D *h2) const
Utility functions.
CSCDetId cscDetId() const
Definition: CSCRecHit2D.h:58
int layer() const
Definition: CSCDetId.h:56
float dPhiIntMax
Definition: CSCSegAlgoRU.h:157
bool areHitsCloseInGlobalPhi(const AlgoState &aState, const CSCRecHit2D *h1, const CSCRecHit2D *h2) const
void flagHitsAsUsed(const AlgoState &aState, const ChamberHitContainer &rechitsInChamber, BoolContainer &used) const
Definition: EPCuts.h:4
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< const CSCRecHit2D * > ChamberHitContainer
Definition: CSCSegAlgoRU.h:51
std::vector< bool > BoolContainer
Definition: CSCSegAlgoRU.h:59
std::vector< int > LayerIndex
Definition: CSCSegAlgoRU.h:50
void tryAddingHitsToSegment(AlgoState &aState, const ChamberHitContainer &rechitsInChamber, const BoolContainer &used, const LayerIndex &layerIndex, const ChamberHitContainerCIt i1, const ChamberHitContainerCIt i2) const
float chi2Norm_2D_
Definition: CSCSegAlgoRU.h:160
const Surface::PositionType & position() const
The position (origin of the R.F.)
Definition: GeomDet.h:43
void updateParameters(AlgoState &aState) const
int channels(unsigned int i) const
Extracting strip channel numbers comprising the rechit - low.
Definition: CSCRecHit2D.h:61
short int hitWire() const
L1A.
Definition: CSCRecHit2D.h:68
bool isSegmentGood(const AlgoState &aState, const ChamberHitContainer &rechitsInChamber) const
ib
Definition: cuy.py:661
bool addHit(AlgoState &aState, const CSCRecHit2D *hit, int layer) const
Utility functions.

◆ compareProtoSegment()

void CSCSegAlgoRU::compareProtoSegment ( AlgoState aState,
const CSCRecHit2D h,
int  layer 
) const
private

Definition at line 859 of file CSCSegAlgoRU.cc.

References CSCSegAlgoRU::AlgoState::aChamber, h, nano_mu_digi_cff::layer, eostools::move(), convertSQLiteXML::ok, CSCSegAlgoRU::AlgoState::proto_segment, replaceHit(), and CSCSegAlgoRU::AlgoState::sfit.

Referenced by tryAddingHitsToSegment().

859  {
860  // Copy the input CSCSegFit
861  std::unique_ptr<CSCSegFit> oldfit;
862  oldfit = std::make_unique<CSCSegFit>(aState.aChamber, aState.proto_segment);
863  oldfit->fit();
864  ChamberHitContainer oldproto = aState.proto_segment;
865 
866  // May create a new fit
867  bool ok = replaceHit(aState, h, layer);
868  if ((aState.sfit->chi2() >= oldfit->chi2()) || !ok) {
869  // keep original fit
870  aState.proto_segment = oldproto;
871  aState.sfit = std::move(oldfit); // reset to the original input fit
872  }
873 }
bool replaceHit(AlgoState &aState, const CSCRecHit2D *h, int layer) const
std::vector< const CSCRecHit2D * > ChamberHitContainer
Definition: CSCSegAlgoRU.h:51
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
def move(src, dest)
Definition: eostools.py:511

◆ fit_r_phi()

float CSCSegAlgoRU::fit_r_phi ( const AlgoState aState,
const SVector6 points,
int  layer 
) const
private

Definition at line 603 of file CSCSegAlgoRU.cc.

References dumpMFGeometry_cfg::delta, mps_fire::i, nano_mu_digi_cff::layer, and slope.

Referenced by isHitNearSegment().

603  {
604  //find R or Phi on the given layer using the given points for the interpolation
605  float Sx = 0;
606  float Sy = 0;
607  float Sxx = 0;
608  float Sxy = 0;
609  for (int i = 1; i < 7; i++) {
610  if (points(i - 1) == 0.)
611  continue;
612  Sy = Sy + (points(i - 1));
613  Sx = Sx + i;
614  Sxx = Sxx + (i * i);
615  Sxy = Sxy + ((i)*points(i - 1));
616  }
617  float delta = 2 * Sxx - Sx * Sx;
618  float intercept = (Sxx * Sy - Sx * Sxy) / delta;
619  float slope = (2 * Sxy - Sx * Sy) / delta;
620  return (intercept + slope * layer);
621 }
static const double slope[3]

◆ fitX()

float CSCSegAlgoRU::fitX ( const AlgoState aState,
SVector6  points,
SVector6  errors,
int  ir,
int  ir2,
float &  chi2_str 
) const
private

Definition at line 804 of file CSCSegAlgoRU.cc.

References HLT_2024v10_cff::chi2_str, dumpMFGeometry_cfg::delta, nano_mu_digi_cff::errors, mps_fire::i, testProducerWithPsetDescEmpty_cfi::i1, and slope.

Referenced by baseline().

805  {
806  float S = 0;
807  float Sx = 0;
808  float Sy = 0;
809  float Sxx = 0;
810  float Sxy = 0;
811  float sigma2 = 0;
812  for (int i = 1; i < 7; i++) {
813  if (i == ir || i == ir2 || points(i - 1) == 0.)
814  continue;
815  sigma2 = errors(i - 1) * errors(i - 1);
816  float i1 = i - 3.5;
817  S = S + (1 / sigma2);
818  Sy = Sy + (points(i - 1) / sigma2);
819  Sx = Sx + ((i1) / sigma2);
820  Sxx = Sxx + (i1 * i1) / sigma2;
821  Sxy = Sxy + (((i1)*points(i - 1)) / sigma2);
822  }
823  float delta = S * Sxx - Sx * Sx;
824  float intercept = (Sxx * Sy - Sx * Sxy) / delta;
825  float slope = (S * Sxy - Sx * Sy) / delta;
826  float chi_str = 0;
827  chi2_str = 0;
828  // calculate chi2_str
829  for (int i = 1; i < 7; i++) {
830  if (i == ir || i == ir2 || points(i - 1) == 0.)
831  continue;
832  chi_str = (points(i - 1) - intercept - slope * (i - 3.5)) / (errors(i - 1));
833  chi2_str = chi2_str + chi_str * chi_str;
834  }
835  return (intercept + slope * 0);
836 }
static const double slope[3]

◆ flagHitsAsUsed()

void CSCSegAlgoRU::flagHitsAsUsed ( const AlgoState aState,
const ChamberHitContainer rechitsInChamber,
BoolContainer used 
) const
private

Flag hits on segment as used

Definition at line 568 of file CSCSegAlgoRU.cc.

References cuy::ib, and CSCSegAlgoRU::AlgoState::proto_segment.

Referenced by buildSegments().

570  {
571  // Flag hits on segment as used
572  ChamberHitContainerCIt ib = rechitsInChamber.begin();
574  for (hi = aState.proto_segment.begin(); hi != aState.proto_segment.end(); ++hi) {
575  for (iu = ib; iu != rechitsInChamber.end(); ++iu) {
576  if (*hi == *iu)
577  used[iu - ib] = true;
578  }
579  }
580 }
std::vector< const CSCRecHit2D * >::const_iterator ChamberHitContainerCIt
Definition: CSCSegAlgoRU.h:52
Definition: EPCuts.h:4
ib
Definition: cuy.py:661

◆ hasHitOnLayer()

bool CSCSegAlgoRU::hasHitOnLayer ( const AlgoState aState,
int  layer 
) const
private

Definition at line 838 of file CSCSegAlgoRU.cc.

References ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, nano_mu_digi_cff::layer, and CSCSegAlgoRU::AlgoState::proto_segment.

Referenced by tryAddingHitsToSegment().

838  {
839  // Is there is already a hit on this layer?
841  for (it = aState.proto_segment.begin(); it != aState.proto_segment.end(); it++)
842  if ((*it)->cscDetId().layer() == layer)
843  return true;
844  return false;
845 }
std::vector< const CSCRecHit2D * >::const_iterator ChamberHitContainerCIt
Definition: CSCSegAlgoRU.h:52

◆ increaseProtoSegment()

void CSCSegAlgoRU::increaseProtoSegment ( AlgoState aState,
const CSCRecHit2D h,
int  layer,
int  chi2_factor 
) const
private

Definition at line 875 of file CSCSegAlgoRU.cc.

References CSCSegAlgoRU::AlgoState::aChamber, addHit(), CSCSegAlgoRU::AlgoState::chi2Max, h, nano_mu_digi_cff::layer, eostools::move(), convertSQLiteXML::ok, CSCSegAlgoRU::AlgoState::proto_segment, and CSCSegAlgoRU::AlgoState::sfit.

Referenced by tryAddingHitsToSegment().

875  {
876  // Creates a new fit
877  std::unique_ptr<CSCSegFit> oldfit;
878  ChamberHitContainer oldproto = aState.proto_segment;
879  oldfit = std::make_unique<CSCSegFit>(aState.aChamber, aState.proto_segment);
880  oldfit->fit();
881 
882  bool ok = addHit(aState, h, layer);
883  //@@ TEST ON ndof<=0 IS JUST TO ACCEPT nhits=2 CASE??
884  if (!ok || ((aState.sfit->ndof() > 0) && (aState.sfit->chi2() / aState.sfit->ndof() >= aState.chi2Max))) {
885  // reset to original fit
886  aState.proto_segment = oldproto;
887  aState.sfit = std::move(oldfit);
888  }
889 }
std::vector< const CSCRecHit2D * > ChamberHitContainer
Definition: CSCSegAlgoRU.h:51
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
def move(src, dest)
Definition: eostools.py:511
bool addHit(AlgoState &aState, const CSCRecHit2D *hit, int layer) const
Utility functions.

◆ isHitNearSegment()

bool CSCSegAlgoRU::isHitNearSegment ( const AlgoState aState,
const CSCRecHit2D h 
) const
private

Definition at line 455 of file CSCSegAlgoRU.cc.

References CSCSegAlgoRU::AlgoState::aChamber, CSCSegAlgoRU::AlgoState::dPhiIntMax, l1ctLayer1_cff::dr, CSCSegAlgoRU::AlgoState::dRIntMax, fit_r_phi(), h, trackingPlots::hp, MainPageGenerator::l, nano_mu_digi_cff::layer, CSCChamber::layer(), M_PI, PV3DBase< T, PVType, FrameType >::perp(), phiAtZ(), CSCSegAlgoRU::AlgoState::proto_segment, dttmaxenums::R, CSCSegAlgoRU::AlgoState::strip_iadd, GeomDet::toGlobal(), and funct::true.

Referenced by tryAddingHitsToSegment().

455  {
456  // Is hit near segment?
457  // Requires deltaphi and deltaR within ranges specified in parameter set.
458  // Note that to make intuitive cuts on delta(phi) one must work in
459  // phi range (-pi, +pi] not [0, 2pi)
460  float strip_width[10] = {0.003878509,
461  0.002958185,
462  0.002327105,
463  0.00152552,
464  0.00465421,
465  0.002327105,
466  0.00465421,
467  0.002327105,
468  0.00465421,
469  0.002327105}; //in rad
470  const CSCLayer* l1 = aState.aChamber->layer((*(aState.proto_segment.begin()))->cscDetId().layer());
471  GlobalPoint gp1 = l1->toGlobal((*(aState.proto_segment.begin()))->localPosition());
472  const CSCLayer* l2 = aState.aChamber->layer((*(aState.proto_segment.begin() + 1))->cscDetId().layer());
473  GlobalPoint gp2 = l2->toGlobal((*(aState.proto_segment.begin() + 1))->localPosition());
474  float err_stpos_h1 = (*(aState.proto_segment.begin()))->errorWithinStrip();
475  float err_stpos_h2 = (*(aState.proto_segment.begin() + 1))->errorWithinStrip();
476  const CSCLayer* l = aState.aChamber->layer(h->cscDetId().layer());
477  GlobalPoint hp = l->toGlobal(h->localPosition());
478  float err_stpos_h = h->errorWithinStrip();
479  float hphi = hp.phi(); // in (-pi, +pi]
480  if (hphi < 0.)
481  hphi += 2. * M_PI; // into range [0, 2pi)
482  float sphi = phiAtZ(aState, hp.z()); // in [0, 2*pi)
483  float phidif = sphi - hphi;
484  if (phidif < 0.)
485  phidif += 2. * M_PI; // into range [0, 2pi)
486  if (phidif > M_PI)
487  phidif -= 2. * M_PI; // into range (-pi, pi]
488  SVector6 r_glob;
489  CSCDetId id = h->cscDetId();
490  int iStn = id.iChamberType() - 1;
491  float dphi_incr = 0;
492  float pos_str = 1;
493  //increase dPhi cut if the hit is on the edge of the strip
494  float stpos = (*h).positionWithinStrip();
495  bool centr_str = false;
496  if (iStn != 0 && iStn != 1) {
497  if (stpos > -0.25 && stpos < 0.25)
498  centr_str = true;
499  }
500  if (err_stpos_h1 < 0.25 * strip_width[iStn] || err_stpos_h2 < 0.25 * strip_width[iStn] ||
501  err_stpos_h < 0.25 * strip_width[iStn]) {
502  dphi_incr = 0.5 * strip_width[iStn];
503  } else {
504  if (centr_str)
505  pos_str = 1.3;
506  }
507  r_glob((*(aState.proto_segment.begin()))->cscDetId().layer() - 1) = gp1.perp();
508  r_glob((*(aState.proto_segment.begin() + 1))->cscDetId().layer() - 1) = gp2.perp();
509  float R = hp.perp();
510  int layer = h->cscDetId().layer();
511  float r_interpolated = fit_r_phi(aState, r_glob, layer);
512  float dr = fabs(r_interpolated - R);
513  float maxWG_width[10] = {0, 0, 4.1, 5.69, 2.49, 5.06, 2.49, 5.06, 1.87, 5.06};
514  //find maxWG_width for ME11 (tilt = 29deg)
515  int wg_num = h->hitWire();
516  if (iStn == 0 || iStn == 1) {
517  if (wg_num == 1) {
518  maxWG_width[0] = 9.25;
519  maxWG_width[1] = 9.25;
520  }
521  if (wg_num > 1 && wg_num < 48) {
522  maxWG_width[0] = 3.14;
523  maxWG_width[1] = 3.14;
524  }
525  if (wg_num == 48) {
526  maxWG_width[0] = 10.75;
527  maxWG_width[1] = 10.75;
528  }
529  }
530  return (fabs(phidif) < aState.dPhiIntMax * aState.strip_iadd * pos_str + dphi_incr &&
531  fabs(dr) < aState.dRIntMax * aState.strip_iadd * maxWG_width[iStn])
532  ? true
533  : false;
534 }
T perp() const
Definition: PV3DBase.h:69
float fit_r_phi(const AlgoState &aState, const SVector6 &points, int layer) const
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
#define M_PI
ROOT::Math::SVector< double, 6 > SVector6
Typedefs.
Definition: CSCSegAlgoRU.h:48
float phiAtZ(const AlgoState &aState, float z) const
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4

◆ isSegmentGood()

bool CSCSegAlgoRU::isSegmentGood ( const AlgoState aState,
const ChamberHitContainer rechitsInChamber 
) const
private

Return true if segment is 'good'. In this algorithm, 'good' means has sufficient hits

Definition at line 551 of file CSCSegAlgoRU.cc.

References convertSQLiteXML::ok, CSCSegAlgoRU::AlgoState::proto_segment, and CSCSegAlgoRU::AlgoState::windowScale.

Referenced by buildSegments().

551  {
552  // If the chamber has 20 hits or fewer, require at least 3 hits on segment
553  // If the chamber has >20 hits require at least 4 hits
554  // If it's the second cycle of the builder and there are <= 12 hits in chamber, require at least 3 hits on segment
555  //@@ THESE VALUES SHOULD BECOME PARAMETERS?
556  bool ok = false;
557  unsigned int iadd = (rechitsInChamber.size() > 20) ? 1 : 0;
558  if (aState.windowScale > 1.) {
559  iadd = 1;
560  if (rechitsInChamber.size() <= 12)
561  iadd = 0;
562  }
563  if (aState.proto_segment.size() >= 3 + iadd)
564  ok = true;
565  return ok;
566 }

◆ phiAtZ()

float CSCSegAlgoRU::phiAtZ ( const AlgoState aState,
float  z 
) const
private

Always enforce direction of segment to point from IP outwards (Incorrect for particles not coming from IP, of course.)

Definition at line 536 of file CSCSegAlgoRU.cc.

References CSCSegAlgoRU::AlgoState::aChamber, f, runTauDisplay::gp, CSCChamber::layer(), M_PI, phi, CSCSegAlgoRU::AlgoState::proto_segment, CSCSegAlgoRU::AlgoState::sfit, GeomDet::toGlobal(), x, PV3DBase< T, PVType, FrameType >::x(), y, PV3DBase< T, PVType, FrameType >::y(), z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by isHitNearSegment().

536  {
537  if (!aState.sfit)
538  return 0.;
539  // Returns a phi in [ 0, 2*pi )
540  const CSCLayer* l1 = aState.aChamber->layer((*(aState.proto_segment.begin()))->cscDetId().layer());
541  GlobalPoint gp = l1->toGlobal(aState.sfit->intercept());
542  GlobalVector gv = l1->toGlobal(aState.sfit->localdir());
543  float x = gp.x() + (gv.x() / gv.z()) * (z - gp.z());
544  float y = gp.y() + (gv.y() / gv.z()) * (z - gp.z());
545  float phi = atan2(y, x);
546  if (phi < 0.f)
547  phi += 2. * M_PI;
548  return phi;
549 }
T z() const
Definition: PV3DBase.h:61
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
double f[11][100]
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
#define M_PI

◆ replaceHit()

bool CSCSegAlgoRU::replaceHit ( AlgoState aState,
const CSCRecHit2D h,
int  layer 
) const
private

Definition at line 847 of file CSCSegAlgoRU.cc.

References addHit(), h, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, nano_mu_digi_cff::layer, and CSCSegAlgoRU::AlgoState::proto_segment.

Referenced by compareProtoSegment().

847  {
848  // replace a hit from a layer
849  ChamberHitContainer::const_iterator it;
850  for (it = aState.proto_segment.begin(); it != aState.proto_segment.end();) {
851  if ((*it)->cscDetId().layer() == layer)
852  it = aState.proto_segment.erase(it);
853  else
854  ++it;
855  }
856  return addHit(aState, h, layer);
857 }
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
bool addHit(AlgoState &aState, const CSCRecHit2D *hit, int layer) const
Utility functions.

◆ run()

std::vector<CSCSegment> CSCSegAlgoRU::run ( const CSCChamber aChamber,
const ChamberHitContainer rechits 
)
inlineoverride

Here we must implement the algorithm

Definition at line 77 of file CSCSegAlgoRU.h.

References buildSegments().

77  {
78  return buildSegments(aChamber, rechits);
79  }
std::vector< CSCSegment > buildSegments(const CSCChamber *aChamber, const ChamberHitContainer &rechits) const
Definition: CSCSegAlgoRU.cc:56

◆ tryAddingHitsToSegment()

void CSCSegAlgoRU::tryAddingHitsToSegment ( AlgoState aState,
const ChamberHitContainer rechitsInChamber,
const BoolContainer used,
const LayerIndex layerIndex,
const ChamberHitContainerCIt  i1,
const ChamberHitContainerCIt  i2 
) const
private

Try adding non-used hits to segment

Definition at line 351 of file CSCSegAlgoRU.cc.

References CSCSegAlgoRU::AlgoState::chi2D_iadd, compareProtoSegment(), h, hasHitOnLayer(), mps_fire::i, testProducerWithPsetDescEmpty_cfi::i1, testProducerWithPsetDescEmpty_cfi::i2, cuy::ib, increaseProtoSegment(), isHitNearSegment(), nano_mu_digi_cff::layer, and CSCSegAlgoRU::AlgoState::proto_segment.

Referenced by buildSegments().

356  {
357  // Iterate over the layers with hits in the chamber
358  // Skip the layers containing the segment endpoints
359  // Test each hit on the other layers to see if it is near the segment
360  // If it is, see whether there is already a hit on the segment from the same layer
361  // - if so, and there are more than 2 hits on the segment, copy the segment,
362  // replace the old hit with the new hit. If the new segment chi2 is better
363  // then replace the original segment with the new one (by swap)
364  // - if not, copy the segment, add the hit. If the new chi2/dof is still satisfactory
365  // then replace the original segment with the new one (by swap)
366  ChamberHitContainerCIt ib = rechits.begin();
367  ChamberHitContainerCIt ie = rechits.end();
368  for (ChamberHitContainerCIt i = ib; i != ie; ++i) {
369  int layer = layerIndex[i - ib];
370  if (hasHitOnLayer(aState, layer) && aState.proto_segment.size() <= 2)
371  continue;
372  if (layerIndex[i - ib] == layerIndex[i1 - ib] || layerIndex[i - ib] == layerIndex[i2 - ib] || used[i - ib])
373  continue;
374 
375  const CSCRecHit2D* h = *i;
376  if (isHitNearSegment(aState, h)) {
377  // Don't consider alternate hits on layers holding the two starting points
378  if (hasHitOnLayer(aState, layer)) {
379  if (aState.proto_segment.size() <= 2)
380  continue;
381  compareProtoSegment(aState, h, layer);
382  } else {
383  increaseProtoSegment(aState, h, layer, aState.chi2D_iadd);
384  }
385  } // h & seg close
386  } // i
387 }
std::vector< const CSCRecHit2D * >::const_iterator ChamberHitContainerCIt
Definition: CSCSegAlgoRU.h:52
bool isHitNearSegment(const AlgoState &aState, const CSCRecHit2D *h) const
void compareProtoSegment(AlgoState &aState, const CSCRecHit2D *h, int layer) const
void increaseProtoSegment(AlgoState &aState, const CSCRecHit2D *h, int layer, int chi2_factor) const
bool hasHitOnLayer(const AlgoState &aState, int layer) const
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
ib
Definition: cuy.py:661

◆ updateParameters()

void CSCSegAlgoRU::updateParameters ( AlgoState aState) const
private

Definition at line 596 of file CSCSegAlgoRU.cc.

References CSCSegAlgoRU::AlgoState::aChamber, CSCSegAlgoRU::AlgoState::proto_segment, and CSCSegAlgoRU::AlgoState::sfit.

Referenced by addHit(), and buildSegments().

596  {
597  // Delete input CSCSegFit, create a new one and make the fit
598  // delete sfit;
599  aState.sfit = std::make_unique<CSCSegFit>(aState.aChamber, aState.proto_segment);
600  aState.sfit->fit();
601 }

Member Data Documentation

◆ chi2_str_

float CSCSegAlgoRU::chi2_str_
private

Definition at line 159 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

◆ chi2Max

float CSCSegAlgoRU::chi2Max
private

Definition at line 158 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

◆ chi2Norm_2D_

float CSCSegAlgoRU::chi2Norm_2D_
private

Definition at line 160 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

◆ debugInfo

bool CSCSegAlgoRU::debugInfo
private

Definition at line 163 of file CSCSegAlgoRU.h.

◆ doCollisions

bool CSCSegAlgoRU::doCollisions
private

Definition at line 153 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

◆ dPhiIntMax

float CSCSegAlgoRU::dPhiIntMax
private

Definition at line 157 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

◆ dPhiMax

float CSCSegAlgoRU::dPhiMax
private

Definition at line 155 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

◆ dRIntMax

float CSCSegAlgoRU::dRIntMax
private

Definition at line 156 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

◆ dRMax

float CSCSegAlgoRU::dRMax
private

Definition at line 154 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

◆ minLayersApart

int CSCSegAlgoRU::minLayersApart
private

Definition at line 162 of file CSCSegAlgoRU.h.

Referenced by CSCSegAlgoRU().

◆ myName

const std::string CSCSegAlgoRU::myName
private

Definition at line 147 of file CSCSegAlgoRU.h.

Referenced by CSCSegAlgoRU().

◆ theChi2

double CSCSegAlgoRU::theChi2
private

Definition at line 149 of file CSCSegAlgoRU.h.

◆ theDirection

LocalVector CSCSegAlgoRU::theDirection
private

Definition at line 151 of file CSCSegAlgoRU.h.

◆ theOrigin

LocalPoint CSCSegAlgoRU::theOrigin
private

Definition at line 150 of file CSCSegAlgoRU.h.

◆ uz

float CSCSegAlgoRU::uz
private

Definition at line 152 of file CSCSegAlgoRU.h.

◆ vz

float CSCSegAlgoRU::vz
private

Definition at line 152 of file CSCSegAlgoRU.h.

◆ wideSeg

float CSCSegAlgoRU::wideSeg
private

Definition at line 161 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().