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
 
bool enlarge
 
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

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

Definition at line 59 of file CSCSegAlgoRU.h.

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

Definition at line 51 of file CSCSegAlgoRU.h.

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

Definition at line 52 of file CSCSegAlgoRU.h.

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

Definition at line 50 of file CSCSegAlgoRU.h.

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

Typedefs.

Definition at line 48 of file CSCSegAlgoRU.h.

Constructor & Destructor Documentation

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

Constructor.

Definition at line 20 of file CSCSegAlgoRU.cc.

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

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

Destructor.

Definition at line 64 of file CSCSegAlgoRU.h.

References buildSegments(), and TrackInfoProducer_cfi::rechits.

64 {};

Member Function Documentation

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

Utility functions.

Definition at line 606 of file CSCSegAlgoRU.cc.

References CSCSegAlgoRU::AlgoState::proto_segment, and updateParameters().

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

606  {
607  // Return true if hit was added successfully
608  // (and then parameters are updated).
609  // Return false if there is already a hit on the same layer, or insert failed.
610  ChamberHitContainer::const_iterator it;
611  for (it = aState.proto_segment.begin(); it != aState.proto_segment.end(); it++)
612  if (((*it)->cscDetId().layer() == layer) && (aHit != (*it)))
613  return false;
614  aState.proto_segment.push_back(aHit);
615  // make a fit
616  updateParameters(aState);
617  return true;
618 }
void updateParameters(AlgoState &aState) const
bool CSCSegAlgoRU::areHitsCloseInGlobalPhi ( const AlgoState aState,
const CSCRecHit2D h1,
const CSCRecHit2D h2 
) const
private

Definition at line 442 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(), CSCChamber::layer(), CSCDetId::layer(), CSCRecHit2D::localPosition(), CSCSegAlgoRU::AlgoState::strip_iadd, and GeomDet::toGlobal().

Referenced by buildSegments().

444  {
445  float strip_width[10] = {0.003878509,
446  0.002958185,
447  0.002327105,
448  0.00152552,
449  0.00465421,
450  0.002327105,
451  0.00465421,
452  0.002327105,
453  0.00465421,
454  0.002327105}; //in rad
455  const CSCLayer* l1 = aState.aChamber->layer(h1->cscDetId().layer());
456  GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
457  const CSCLayer* l2 = aState.aChamber->layer(h2->cscDetId().layer());
458  GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
459  float err_stpos_h1 = h1->errorWithinStrip();
460  float err_stpos_h2 = h2->errorWithinStrip();
461  CSCDetId id = h1->cscDetId();
462  int iStn = id.iChamberType() - 1;
463  float dphi_incr = 0;
464  if (err_stpos_h1 > 0.25 * strip_width[iStn] || err_stpos_h2 > 0.25 * strip_width[iStn])
465  dphi_incr = 0.5 * strip_width[iStn];
466  float dphi12 = deltaPhi(gp1.barePhi(), gp2.barePhi());
467  return (fabs(dphi12) < (aState.dPhiMax * aState.strip_iadd + dphi_incr)) ? true : false; // +v
468 }
CSCDetId cscDetId() const
Definition: CSCRecHit2D.h:58
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
int layer() const
Definition: CSCDetId.h:56
T barePhi() const
Definition: PV3DBase.h:65
unsigned short iChamberType() const
Definition: CSCDetId.h:96
float errorWithinStrip() const
The uncertainty of the estimated position within the strip.
Definition: CSCRecHit2D.h:85
bool CSCSegAlgoRU::areHitsCloseInR ( const AlgoState aState,
const CSCRecHit2D h1,
const CSCRecHit2D h2 
) const
private

Utility functions.

Definition at line 395 of file CSCSegAlgoRU.cc.

References CSCSegAlgoRU::AlgoState::aChamber, CSCRecHit2D::cscDetId(), CSCSegAlgoRU::AlgoState::doCollisions, CSCSegAlgoRU::AlgoState::dRMax, CSCSegAlgoRU::AlgoState::enlarge, 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().

395  {
396  float maxWG_width[10] = {0, 0, 4.1, 5.69, 2.49, 5.06, 2.49, 5.06, 1.87, 5.06};
397  CSCDetId id = h1->cscDetId();
398  int iStn = id.iChamberType() - 1;
399  //find maxWG_width for ME11 (tilt = 29deg)
400  int wg_num = h2->hitWire();
401  if (iStn == 0 || iStn == 1) {
402  if (wg_num == 1) {
403  maxWG_width[0] = 9.25;
404  maxWG_width[1] = 9.25;
405  }
406  if (wg_num > 1 && wg_num < 48) {
407  maxWG_width[0] = 3.14;
408  maxWG_width[1] = 3.14;
409  }
410  if (wg_num == 48) {
411  maxWG_width[0] = 10.75;
412  maxWG_width[1] = 10.75;
413  }
414  }
415  const CSCLayer* l1 = aState.aChamber->layer(h1->cscDetId().layer());
416  GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
417  const CSCLayer* l2 = aState.aChamber->layer(h2->cscDetId().layer());
418  GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
419  //find z to understand the direction
420  float h1z = gp1.z();
421  float h2z = gp2.z();
422  //switch off the IP check for non collisions case
423  if (!aState.doCollisions) {
424  h1z = 1;
425  h2z = 1;
426  }
427 
428  if (aState.enlarge) {
429  return (gp2.perp() > ((gp1.perp() - aState.dRMax * aState.strip_iadd * maxWG_width[iStn]) * h2z) / h1z &&
430  gp2.perp() < ((gp1.perp() + aState.dRMax * aState.strip_iadd * maxWG_width[iStn]) * h2z) / h1z)
431  ? true
432  : false;
433 
434  } else {
435  return (gp2.perp() > ((gp1.perp() - aState.dRMax * maxWG_width[iStn]) * h2z) / h1z &&
436  gp2.perp() < ((gp1.perp() + aState.dRMax * maxWG_width[iStn]) * h2z) / h1z)
437  ? true
438  : false;
439  }
440 }
CSCDetId cscDetId() const
Definition: CSCRecHit2D.h:58
T perp() const
Definition: PV3DBase.h:69
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
int layer() const
Definition: CSCDetId.h:56
T z() const
Definition: PV3DBase.h:61
unsigned short iChamberType() const
Definition: CSCDetId.h:96
short int hitWire() const
L1A.
Definition: CSCRecHit2D.h:68
void CSCSegAlgoRU::baseline ( AlgoState aState,
int  n_seg_min 
) const
private

(nhits-2)

(nhits-3)

Definition at line 647 of file CSCSegAlgoRU.cc.

References edmScanValgrind::buffer, CSCRecHit2D::channels(), CSCSegmentAlgorithmRU_cfi::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(), min(), nhits, CSCRecHit2D::nStrips(), CSCSegAlgoRU::AlgoState::proto_segment, CSCDetId::ring(), slope, CSCDetId::station(), z, and testProducerWithPsetDescEmpty_cfi::z2.

Referenced by buildSegments().

647  {
648  int nhits = aState.proto_segment.size();
649  //initialise vectors for strip position and error within strip
650  SVector6 sp;
651  SVector6 se;
652  unsigned int init_size = aState.proto_segment.size();
654  buffer.clear();
655  buffer.reserve(init_size);
656  while (buffer.size() < init_size) {
657  ChamberHitContainer::iterator min;
658  int min_layer = 10;
659  for (ChamberHitContainer::iterator k = aState.proto_segment.begin(); k != aState.proto_segment.end(); k++) {
660  const CSCRecHit2D* iRHk = *k;
661  CSCDetId idRHk = iRHk->cscDetId();
662  int kLayer = idRHk.layer();
663  if (kLayer < min_layer) {
664  min_layer = kLayer;
665  min = k;
666  }
667  }
668  buffer.push_back(*min);
669  aState.proto_segment.erase(min);
670  } //while
671 
672  aState.proto_segment.clear();
673  for (ChamberHitContainer::const_iterator cand = buffer.begin(); cand != buffer.end(); cand++) {
674  aState.proto_segment.push_back(*cand);
675  }
676 
677  for (ChamberHitContainer::const_iterator iRH = aState.proto_segment.begin(); iRH != aState.proto_segment.end();
678  iRH++) {
679  const CSCRecHit2D* iRHp = *iRH;
680  CSCDetId idRH = iRHp->cscDetId();
681  int kRing = idRH.ring();
682  int kStation = idRH.station();
683  int kLayer = idRH.layer();
684  // Find the strip containing this hit
685  int centerid = iRHp->nStrips() / 2;
686  int centerStrip = iRHp->channels(centerid);
687  float stpos = (*iRHp).positionWithinStrip();
688  se(kLayer - 1) = (*iRHp).errorWithinStrip();
689  // Take into account half-strip staggering of layers (ME1/1 has no staggering)
690  if (kStation == 1 && (kRing == 1 || kRing == 4))
691  sp(kLayer - 1) = stpos + centerStrip;
692  else {
693  if (kLayer == 1 || kLayer == 3 || kLayer == 5)
694  sp(kLayer - 1) = stpos + centerStrip;
695  if (kLayer == 2 || kLayer == 4 || kLayer == 6)
696  sp(kLayer - 1) = stpos - 0.5 + centerStrip;
697  }
698  }
699  float chi2_str;
700  fitX(aState, sp, se, -1, -1, chi2_str);
701 
702  //-----------------------------------------------------
703  // Optimal point rejection method
704  //-----------------------------------------------------
705  float minSum = 1000;
706  int i1b = 0;
707  int i2b = 0;
708  int iworst = -1;
709  int bad_layer = -1;
710  ChamberHitContainer::const_iterator rh_to_be_deleted_1;
711  ChamberHitContainer::const_iterator rh_to_be_deleted_2;
712  if ((chi2_str) > aState.chi2_str_ * aState.chi2D_iadd) {
713  for (ChamberHitContainer::const_iterator i1 = aState.proto_segment.begin(); i1 != aState.proto_segment.end();
714  ++i1) {
715  ++i1b;
716  const CSCRecHit2D* i1_1 = *i1;
717  CSCDetId idRH1 = i1_1->cscDetId();
718  int z1 = idRH1.layer();
719  i2b = i1b;
720  for (ChamberHitContainer::const_iterator i2 = i1 + 1; i2 != aState.proto_segment.end(); ++i2) {
721  ++i2b;
722  const CSCRecHit2D* i2_1 = *i2;
723  CSCDetId idRH2 = i2_1->cscDetId();
724  int z2 = idRH2.layer();
725  int irej = 0;
726  for (ChamberHitContainer::const_iterator ir = aState.proto_segment.begin(); ir != aState.proto_segment.end();
727  ++ir) {
728  ++irej;
729  if (ir == i1 || ir == i2)
730  continue;
731  float dsum = 0;
732  int hit_nr = 0;
733  const CSCRecHit2D* ir_1 = *ir;
734  CSCDetId idRH = ir_1->cscDetId();
735  int worst_layer = idRH.layer();
736  for (ChamberHitContainer::const_iterator i = aState.proto_segment.begin(); i != aState.proto_segment.end();
737  ++i) {
738  ++hit_nr;
739  const CSCRecHit2D* i_1 = *i;
740  if (i == i1 || i == i2 || i == ir)
741  continue;
742  float slope = (sp(z2 - 1) - sp(z1 - 1)) / (z2 - z1);
743  float intersept = sp(z1 - 1) - slope * z1;
744  CSCDetId idRH = i_1->cscDetId();
745  int z = idRH.layer();
746  float di = fabs(sp(z - 1) - intersept - slope * z);
747  dsum = dsum + di;
748  } //i
749  if (dsum < minSum) {
750  minSum = dsum;
751  bad_layer = worst_layer;
752  iworst = irej;
753  rh_to_be_deleted_1 = ir;
754  }
755  } //ir
756  } //i2
757  } //i1
758  fitX(aState, sp, se, bad_layer, -1, chi2_str);
759  } //if chi2prob<1.0e-4
760 
761  //find worst from n-1 hits
762  int iworst2 = -1;
763  int bad_layer2 = -1;
764  if (iworst > -1 && (nhits - 1) > n_seg_min && (chi2_str) > aState.chi2_str_ * aState.chi2D_iadd) {
765  iworst = -1;
766  float minSum = 1000;
767  int i1b = 0;
768  int i2b = 0;
769  for (ChamberHitContainer::const_iterator i1 = aState.proto_segment.begin(); i1 != aState.proto_segment.end();
770  ++i1) {
771  ++i1b;
772  const CSCRecHit2D* i1_1 = *i1;
773  CSCDetId idRH1 = i1_1->cscDetId();
774  int z1 = idRH1.layer();
775  i2b = i1b;
776  for (ChamberHitContainer::const_iterator i2 = i1 + 1; i2 != aState.proto_segment.end(); ++i2) {
777  ++i2b;
778  const CSCRecHit2D* i2_1 = *i2;
779  CSCDetId idRH2 = i2_1->cscDetId();
780  int z2 = idRH2.layer();
781  int irej = 0;
782  for (ChamberHitContainer::const_iterator ir = aState.proto_segment.begin(); ir != aState.proto_segment.end();
783  ++ir) {
784  ++irej;
785  int irej2 = 0;
786  if (ir == i1 || ir == i2)
787  continue;
788  const CSCRecHit2D* ir_1 = *ir;
789  CSCDetId idRH = ir_1->cscDetId();
790  int worst_layer = idRH.layer();
791  for (ChamberHitContainer::const_iterator ir2 = aState.proto_segment.begin();
792  ir2 != aState.proto_segment.end();
793  ++ir2) {
794  ++irej2;
795  if (ir2 == i1 || ir2 == i2 || ir2 == ir)
796  continue;
797  float dsum = 0;
798  int hit_nr = 0;
799  const CSCRecHit2D* ir2_1 = *ir2;
800  CSCDetId idRH = ir2_1->cscDetId();
801  int worst_layer2 = idRH.layer();
802  for (ChamberHitContainer::const_iterator i = aState.proto_segment.begin(); i != aState.proto_segment.end();
803  ++i) {
804  ++hit_nr;
805  const CSCRecHit2D* i_1 = *i;
806  if (i == i1 || i == i2 || i == ir || i == ir2)
807  continue;
808  float slope = (sp(z2 - 1) - sp(z1 - 1)) / (z2 - z1);
809  float intersept = sp(z1 - 1) - slope * z1;
810  CSCDetId idRH = i_1->cscDetId();
811  int z = idRH.layer();
812  float di = fabs(sp(z - 1) - intersept - slope * z);
813  dsum = dsum + di;
814  } //i
815  if (dsum < minSum) {
816  minSum = dsum;
817  iworst2 = irej2;
818  iworst = irej;
819  bad_layer = worst_layer;
820  bad_layer2 = worst_layer2;
821  rh_to_be_deleted_1 = ir;
822  rh_to_be_deleted_2 = ir2;
823  }
824  } //ir2
825  } //ir
826  } //i2
827  } //i1
828  fitX(aState, sp, se, bad_layer, bad_layer2, chi2_str);
829  } //if prob(n-1)<e-4
830 
831  //----------------------------------
832  //erase bad_hits
833  //----------------------------------
834  if (iworst2 - 1 >= 0 && iworst2 <= int(aState.proto_segment.size())) {
835  aState.proto_segment.erase(rh_to_be_deleted_2);
836  }
837  if (iworst - 1 >= 0 && iworst <= int(aState.proto_segment.size())) {
838  aState.proto_segment.erase(rh_to_be_deleted_1);
839  }
840 }
CSCDetId cscDetId() const
Definition: CSCRecHit2D.h:58
static const double slope[3]
int layer() const
Definition: CSCDetId.h:56
int channels(unsigned int i) const
Extracting strip channel numbers comprising the rechit - low.
Definition: CSCRecHit2D.h:61
unsigned int nStrips() const
Definition: CSCRecHit2D.h:62
T min(T a, T b)
Definition: MathUtil.h:58
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 ring() const
Definition: CSCDetId.h:68
float fitX(const AlgoState &aState, SVector6 points, SVector6 errors, int ir, int ir2, float &chi2_str) const
int station() const
Definition: CSCDetId.h:79
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 55 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, enlarge, flagHitsAsUsed(), CSCRecHit2D::hitWire(), mps_fire::i, testProducerWithPsetDescEmpty_cfi::i1, testProducerWithPsetDescEmpty_cfi::i2, cuy::ib, isSegmentGood(), dqmiolumiharvest::j, dqmdumpme::k, CSCDetId::layer(), CSCChamber::layer(), CSCRecHit2D::nStrips(), GeomDet::position(), TrackInfoProducer_cfi::rechits, groupFilesInBlocks::reverse, groupFilesInBlocks::temp, tryAddingHitsToSegment(), updateParameters(), wideSeg, and PV3DBase< T, PVType, FrameType >::z().

Referenced by run(), and ~CSCSegAlgoRU().

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

Definition at line 897 of file CSCSegAlgoRU.cc.

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

Referenced by tryAddingHitsToSegment().

897  {
898  // Copy the input CSCSegFit
899  std::unique_ptr<CSCSegFit> oldfit;
900  oldfit.reset(new CSCSegFit(aState.aChamber, aState.proto_segment));
901  oldfit->fit();
902  ChamberHitContainer oldproto = aState.proto_segment;
903 
904  // May create a new fit
905  bool ok = replaceHit(aState, h, layer);
906  if ((aState.sfit->chi2() >= oldfit->chi2()) || !ok) {
907  // keep original fit
908  aState.proto_segment = oldproto;
909  aState.sfit = std::move(oldfit); // reset to the original input fit
910  }
911 }
bool replaceHit(AlgoState &aState, const CSCRecHit2D *h, int layer) const
std::vector< const CSCRecHit2D * > ChamberHitContainer
Definition: CSCSegAlgoRU.h:51
def move(src, dest)
Definition: eostools.py:511
float CSCSegAlgoRU::fit_r_phi ( const AlgoState aState,
const SVector6 points,
int  layer 
) const
private

Definition at line 627 of file CSCSegAlgoRU.cc.

References dumpMFGeometry_cfg::delta, mps_fire::i, HLT_2018_cff::points, and slope.

Referenced by isHitNearSegment().

627  {
628  //find R or Phi on the given layer using the given points for the interpolation
629  float Sx = 0;
630  float Sy = 0;
631  float Sxx = 0;
632  float Sxy = 0;
633  for (int i = 1; i < 7; i++) {
634  if (points(i - 1) == 0.)
635  continue;
636  Sy = Sy + (points(i - 1));
637  Sx = Sx + i;
638  Sxx = Sxx + (i * i);
639  Sxy = Sxy + ((i)*points(i - 1));
640  }
641  float delta = 2 * Sxx - Sx * Sx;
642  float intercept = (Sxx * Sy - Sx * Sxy) / delta;
643  float slope = (2 * Sxy - Sx * Sy) / delta;
644  return (intercept + slope * layer);
645 }
static const double slope[3]
float CSCSegAlgoRU::fitX ( const AlgoState aState,
SVector6  points,
SVector6  errors,
int  ir,
int  ir2,
float &  chi2_str 
) const
private

Definition at line 842 of file CSCSegAlgoRU.cc.

References dumpMFGeometry_cfg::delta, MessageLogger_cfi::errors, mps_fire::i, testProducerWithPsetDescEmpty_cfi::i1, HLT_2018_cff::points, and slope.

Referenced by baseline().

843  {
844  float S = 0;
845  float Sx = 0;
846  float Sy = 0;
847  float Sxx = 0;
848  float Sxy = 0;
849  float sigma2 = 0;
850  for (int i = 1; i < 7; i++) {
851  if (i == ir || i == ir2 || points(i - 1) == 0.)
852  continue;
853  sigma2 = errors(i - 1) * errors(i - 1);
854  float i1 = i - 3.5;
855  S = S + (1 / sigma2);
856  Sy = Sy + (points(i - 1) / sigma2);
857  Sx = Sx + ((i1) / sigma2);
858  Sxx = Sxx + (i1 * i1) / sigma2;
859  Sxy = Sxy + (((i1)*points(i - 1)) / sigma2);
860  }
861  float delta = S * Sxx - Sx * Sx;
862  float intercept = (Sxx * Sy - Sx * Sxy) / delta;
863  float slope = (S * Sxy - Sx * Sy) / delta;
864  float chi_str = 0;
865  chi2_str = 0;
866  // calculate chi2_str
867  for (int i = 1; i < 7; i++) {
868  if (i == ir || i == ir2 || points(i - 1) == 0.)
869  continue;
870  chi_str = (points(i - 1) - intercept - slope * (i - 3.5)) / (errors(i - 1));
871  chi2_str = chi2_str + chi_str * chi_str;
872  }
873  return (intercept + slope * 0);
874 }
static const double slope[3]
void CSCSegAlgoRU::flagHitsAsUsed ( const AlgoState aState,
const ChamberHitContainer rechitsInChamber,
BoolContainer used 
) const
private

Flag hits on segment as used

Definition at line 592 of file CSCSegAlgoRU.cc.

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

Referenced by buildSegments().

594  {
595  // Flag hits on segment as used
596  ChamberHitContainerCIt ib = rechitsInChamber.begin();
598  for (hi = aState.proto_segment.begin(); hi != aState.proto_segment.end(); ++hi) {
599  for (iu = ib; iu != rechitsInChamber.end(); ++iu) {
600  if (*hi == *iu)
601  used[iu - ib] = true;
602  }
603  }
604 }
std::vector< const CSCRecHit2D * >::const_iterator ChamberHitContainerCIt
Definition: CSCSegAlgoRU.h:52
ib
Definition: cuy.py:662
bool CSCSegAlgoRU::hasHitOnLayer ( const AlgoState aState,
int  layer 
) const
private

Definition at line 876 of file CSCSegAlgoRU.cc.

References CSCSegAlgoRU::AlgoState::proto_segment.

Referenced by tryAddingHitsToSegment().

876  {
877  // Is there is already a hit on this layer?
879  for (it = aState.proto_segment.begin(); it != aState.proto_segment.end(); it++)
880  if ((*it)->cscDetId().layer() == layer)
881  return true;
882  return false;
883 }
std::vector< const CSCRecHit2D * >::const_iterator ChamberHitContainerCIt
Definition: CSCSegAlgoRU.h:52
void CSCSegAlgoRU::increaseProtoSegment ( AlgoState aState,
const CSCRecHit2D h,
int  layer,
int  chi2_factor 
) const
private

Definition at line 913 of file CSCSegAlgoRU.cc.

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

Referenced by tryAddingHitsToSegment().

913  {
914  // Creates a new fit
915  std::unique_ptr<CSCSegFit> oldfit;
916  ChamberHitContainer oldproto = aState.proto_segment;
917  oldfit.reset(new CSCSegFit(aState.aChamber, aState.proto_segment));
918  oldfit->fit();
919 
920  bool ok = addHit(aState, h, layer);
921  //@@ TEST ON ndof<=0 IS JUST TO ACCEPT nhits=2 CASE??
922  if (!ok || ((aState.sfit->ndof() > 0) && (aState.sfit->chi2() / aState.sfit->ndof() >= aState.chi2Max))) {
923  // reset to original fit
924  aState.proto_segment = oldproto;
925  aState.sfit = std::move(oldfit);
926  }
927 }
std::vector< const CSCRecHit2D * > ChamberHitContainer
Definition: CSCSegAlgoRU.h:51
bool addHit(AlgoState &aState, const CSCRecHit2D *hit, int layer) const
Utility functions.
def move(src, dest)
Definition: eostools.py:511
bool CSCSegAlgoRU::isHitNearSegment ( const AlgoState aState,
const CSCRecHit2D h 
) const
private

Definition at line 470 of file CSCSegAlgoRU.cc.

References CSCSegAlgoRU::AlgoState::aChamber, CSCRecHit2D::cscDetId(), CSCSegAlgoRU::AlgoState::dPhiIntMax, flavorHistoryFilter_cfi::dr, CSCSegAlgoRU::AlgoState::dRIntMax, CSCSegAlgoRU::AlgoState::enlarge, CSCRecHit2D::errorWithinStrip(), fit_r_phi(), CSCRecHit2D::hitWire(), AnalysisDataFormats_SUSYBSMObjects::hp, CSCDetId::iChamberType(), cmsLHEtoEOSManager::l, CSCChamber::layer(), CSCDetId::layer(), CSCRecHit2D::localPosition(), M_PI, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), phiAtZ(), CSCSegAlgoRU::AlgoState::proto_segment, dttmaxenums::R, CSCSegAlgoRU::AlgoState::strip_iadd, GeomDet::toGlobal(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by tryAddingHitsToSegment().

470  {
471  // Is hit near segment?
472  // Requires deltaphi and deltaR within ranges specified in parameter set.
473  // Note that to make intuitive cuts on delta(phi) one must work in
474  // phi range (-pi, +pi] not [0, 2pi)
475  float strip_width[10] = {0.003878509,
476  0.002958185,
477  0.002327105,
478  0.00152552,
479  0.00465421,
480  0.002327105,
481  0.00465421,
482  0.002327105,
483  0.00465421,
484  0.002327105}; //in rad
485  const CSCLayer* l1 = aState.aChamber->layer((*(aState.proto_segment.begin()))->cscDetId().layer());
486  GlobalPoint gp1 = l1->toGlobal((*(aState.proto_segment.begin()))->localPosition());
487  const CSCLayer* l2 = aState.aChamber->layer((*(aState.proto_segment.begin() + 1))->cscDetId().layer());
488  GlobalPoint gp2 = l2->toGlobal((*(aState.proto_segment.begin() + 1))->localPosition());
489  float err_stpos_h1 = (*(aState.proto_segment.begin()))->errorWithinStrip();
490  float err_stpos_h2 = (*(aState.proto_segment.begin() + 1))->errorWithinStrip();
491  const CSCLayer* l = aState.aChamber->layer(h->cscDetId().layer());
493  float err_stpos_h = h->errorWithinStrip();
494  float hphi = hp.phi(); // in (-pi, +pi]
495  if (hphi < 0.)
496  hphi += 2. * M_PI; // into range [0, 2pi)
497  float sphi = phiAtZ(aState, hp.z()); // in [0, 2*pi)
498  float phidif = sphi - hphi;
499  if (phidif < 0.)
500  phidif += 2. * M_PI; // into range [0, 2pi)
501  if (phidif > M_PI)
502  phidif -= 2. * M_PI; // into range (-pi, pi]
503  SVector6 r_glob;
504  CSCDetId id = h->cscDetId();
505  int iStn = id.iChamberType() - 1;
506  float dphi_incr = 0;
507  float pos_str = 1;
508  //increase dPhi cut if the hit is on the edge of the strip
509  float stpos = (*h).positionWithinStrip();
510  bool centr_str = false;
511  if (iStn != 0 && iStn != 1) {
512  if (stpos > -0.25 && stpos < 0.25)
513  centr_str = true;
514  }
515  if (err_stpos_h1 < 0.25 * strip_width[iStn] || err_stpos_h2 < 0.25 * strip_width[iStn] ||
516  err_stpos_h < 0.25 * strip_width[iStn]) {
517  dphi_incr = 0.5 * strip_width[iStn];
518  } else {
519  if (centr_str)
520  pos_str = 1.3;
521  }
522  r_glob((*(aState.proto_segment.begin()))->cscDetId().layer() - 1) = gp1.perp();
523  r_glob((*(aState.proto_segment.begin() + 1))->cscDetId().layer() - 1) = gp2.perp();
524  float R = hp.perp();
525  int layer = h->cscDetId().layer();
526  float r_interpolated = fit_r_phi(aState, r_glob, layer);
527  float dr = fabs(r_interpolated - R);
528  float maxWG_width[10] = {0, 0, 4.1, 5.69, 2.49, 5.06, 2.49, 5.06, 1.87, 5.06};
529  //find maxWG_width for ME11 (tilt = 29deg)
530  int wg_num = h->hitWire();
531  if (iStn == 0 || iStn == 1) {
532  if (wg_num == 1) {
533  maxWG_width[0] = 9.25;
534  maxWG_width[1] = 9.25;
535  }
536  if (wg_num > 1 && wg_num < 48) {
537  maxWG_width[0] = 3.14;
538  maxWG_width[1] = 3.14;
539  }
540  if (wg_num == 48) {
541  maxWG_width[0] = 10.75;
542  maxWG_width[1] = 10.75;
543  }
544  }
545 
546  if (aState.enlarge) {
547  return (fabs(phidif) < aState.dPhiIntMax * aState.strip_iadd * pos_str + dphi_incr &&
548  fabs(dr) < aState.dRIntMax * aState.strip_iadd * maxWG_width[iStn])
549  ? true
550  : false;
551 
552  } else {
553  return (fabs(phidif) < aState.dPhiIntMax * aState.strip_iadd * pos_str + dphi_incr &&
554  fabs(dr) < aState.dRIntMax * maxWG_width[iStn])
555  ? true
556  : false;
557  }
558 }
CSCDetId cscDetId() const
Definition: CSCRecHit2D.h:58
T perp() const
Definition: PV3DBase.h:69
float phiAtZ(const AlgoState &aState, float z) const
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
LocalPoint localPosition() const override
Definition: CSCRecHit2D.h:56
int layer() const
Definition: CSCDetId.h:56
susybsm::HSCParticleRefProd hp
Definition: classes.h:27
T z() const
Definition: PV3DBase.h:61
unsigned short iChamberType() const
Definition: CSCDetId.h:96
#define M_PI
ROOT::Math::SVector< double, 6 > SVector6
Typedefs.
Definition: CSCSegAlgoRU.h:48
float fit_r_phi(const AlgoState &aState, const SVector6 &points, int layer) const
short int hitWire() const
L1A.
Definition: CSCRecHit2D.h:68
float errorWithinStrip() const
The uncertainty of the estimated position within the strip.
Definition: CSCRecHit2D.h:85
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 575 of file CSCSegAlgoRU.cc.

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

Referenced by buildSegments().

575  {
576  // If the chamber has 20 hits or fewer, require at least 3 hits on segment
577  // If the chamber has >20 hits require at least 4 hits
578  // If it's the second cycle of the builder and there are <= 12 hits in chamber, require at least 3 hits on segment
579  //@@ THESE VALUES SHOULD BECOME PARAMETERS?
580  bool ok = false;
581  unsigned int iadd = (rechitsInChamber.size() > 20) ? 1 : 0;
582  if (aState.windowScale > 1.) {
583  iadd = 1;
584  if (rechitsInChamber.size() <= 12 && aState.enlarge)
585  iadd = 0;
586  }
587  if (aState.proto_segment.size() >= 3 + iadd)
588  ok = true;
589  return ok;
590 }
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 560 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(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by isHitNearSegment().

560  {
561  if (!aState.sfit)
562  return 0.;
563  // Returns a phi in [ 0, 2*pi )
564  const CSCLayer* l1 = aState.aChamber->layer((*(aState.proto_segment.begin()))->cscDetId().layer());
565  GlobalPoint gp = l1->toGlobal(aState.sfit->intercept());
566  GlobalVector gv = l1->toGlobal(aState.sfit->localdir());
567  float x = gp.x() + (gv.x() / gv.z()) * (z - gp.z());
568  float y = gp.y() + (gv.y() / gv.z()) * (z - gp.z());
569  float phi = atan2(y, x);
570  if (phi < 0.f)
571  phi += 2. * M_PI;
572  return phi;
573 }
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
T y() const
Definition: PV3DBase.h:60
T z() const
Definition: PV3DBase.h:61
double f[11][100]
#define M_PI
T x() const
Definition: PV3DBase.h:59
bool CSCSegAlgoRU::replaceHit ( AlgoState aState,
const CSCRecHit2D h,
int  layer 
) const
private

Definition at line 885 of file CSCSegAlgoRU.cc.

References addHit(), and CSCSegAlgoRU::AlgoState::proto_segment.

Referenced by compareProtoSegment().

885  {
886  // replace a hit from a layer
887  ChamberHitContainer::const_iterator it;
888  for (it = aState.proto_segment.begin(); it != aState.proto_segment.end();) {
889  if ((*it)->cscDetId().layer() == layer)
890  it = aState.proto_segment.erase(it);
891  else
892  ++it;
893  }
894  return addHit(aState, h, layer);
895 }
bool addHit(AlgoState &aState, const CSCRecHit2D *hit, int layer) const
Utility functions.
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:55
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 357 of file CSCSegAlgoRU.cc.

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

Referenced by buildSegments().

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

Definition at line 620 of file CSCSegAlgoRU.cc.

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

Referenced by addHit(), and buildSegments().

620  {
621  // Delete input CSCSegFit, create a new one and make the fit
622  // delete sfit;
623  aState.sfit.reset(new CSCSegFit(aState.aChamber, aState.proto_segment));
624  aState.sfit->fit();
625 }

Member Data Documentation

float CSCSegAlgoRU::chi2_str_
private

Definition at line 160 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

float CSCSegAlgoRU::chi2Max
private

Definition at line 159 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

float CSCSegAlgoRU::chi2Norm_2D_
private

Definition at line 161 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

bool CSCSegAlgoRU::debugInfo
private

Definition at line 164 of file CSCSegAlgoRU.h.

bool CSCSegAlgoRU::doCollisions
private

Definition at line 154 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

float CSCSegAlgoRU::dPhiIntMax
private

Definition at line 158 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

float CSCSegAlgoRU::dPhiMax
private

Definition at line 156 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

float CSCSegAlgoRU::dRIntMax
private

Definition at line 157 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

float CSCSegAlgoRU::dRMax
private

Definition at line 155 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

bool CSCSegAlgoRU::enlarge
private

Definition at line 165 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().

int CSCSegAlgoRU::minLayersApart
private

Definition at line 163 of file CSCSegAlgoRU.h.

Referenced by CSCSegAlgoRU().

const std::string CSCSegAlgoRU::myName
private

Definition at line 148 of file CSCSegAlgoRU.h.

Referenced by CSCSegAlgoRU().

double CSCSegAlgoRU::theChi2
private

Definition at line 150 of file CSCSegAlgoRU.h.

LocalVector CSCSegAlgoRU::theDirection
private

Definition at line 152 of file CSCSegAlgoRU.h.

LocalPoint CSCSegAlgoRU::theOrigin
private

Definition at line 151 of file CSCSegAlgoRU.h.

float CSCSegAlgoRU::uz
private

Definition at line 153 of file CSCSegAlgoRU.h.

float CSCSegAlgoRU::vz
private

Definition at line 153 of file CSCSegAlgoRU.h.

float CSCSegAlgoRU::wideSeg
private

Definition at line 162 of file CSCSegAlgoRU.h.

Referenced by buildSegments(), and CSCSegAlgoRU().