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)
 
virtual ~CSCSegAlgoRU ()
 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 36 of file CSCSegAlgoRU.h.

Member Typedef Documentation

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

Definition at line 63 of file CSCSegAlgoRU.h.

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

Definition at line 55 of file CSCSegAlgoRU.h.

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

Definition at line 56 of file CSCSegAlgoRU.h.

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

Definition at line 54 of file CSCSegAlgoRU.h.

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

Typedefs.

Definition at line 52 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, edm::ParameterSet::getParameter(), LogDebug, minLayersApart, myName, and wideSeg.

21  : CSCSegmentAlgorithm(ps), myName("CSCSegAlgoRU"){
22  doCollisions = ps.getParameter<bool>("doCollisions");
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:147
float dPhiIntMax
Definition: CSCSegAlgoRU.h:158
float chi2Norm_2D_
Definition: CSCSegAlgoRU.h:161
virtual CSCSegAlgoRU::~CSCSegAlgoRU ( )
inlinevirtual

Destructor.

Definition at line 68 of file CSCSegAlgoRU.h.

References buildSegments(), and TrackInfoProducer_cfi::rechits.

68 {};

Member Function Documentation

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

Utility functions.

Definition at line 521 of file CSCSegAlgoRU.cc.

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

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

521  {
522  // Return true if hit was added successfully
523  // (and then parameters are updated).
524  // Return false if there is already a hit on the same layer, or insert failed.
525  ChamberHitContainer::const_iterator it;
526  for(it = aState.proto_segment.begin(); it != aState.proto_segment.end(); it++)
527  if (((*it)->cscDetId().layer() == layer) && (aHit != (*it)))
528  return false;
529  aState.proto_segment.push_back(aHit);
530  // make a fit
531  updateParameters(aState);
532  return true;
533 }
void updateParameters(AlgoState &aState) const
bool CSCSegAlgoRU::areHitsCloseInGlobalPhi ( const AlgoState aState,
const CSCRecHit2D h1,
const CSCRecHit2D h2 
) const
private

Definition at line 400 of file CSCSegAlgoRU.cc.

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

Referenced by buildSegments().

400  {
401  float strip_width[10] = {0.003878509, 0.002958185, 0.002327105, 0.00152552, 0.00465421, 0.002327105, 0.00465421, 0.002327105, 0.00465421, 0.002327105};//in rad
402  const CSCLayer* l1 = aState.aChamber->layer(h1->cscDetId().layer());
403  GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
404  const CSCLayer* l2 = aState.aChamber->layer(h2->cscDetId().layer());
405  GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
406  float err_stpos_h1 = h1->errorWithinStrip();
407  float err_stpos_h2 = h2->errorWithinStrip();
408  CSCDetId id = h1->cscDetId();
409  int iStn = id.iChamberType()-1;
410  float dphi_incr = 0;
411  if(err_stpos_h1>0.25*strip_width[iStn] || err_stpos_h2>0.25*strip_width[iStn])dphi_incr = 0.5*strip_width[iStn];
412  float dphi12 = deltaPhi(gp1.barePhi(),gp2.barePhi());
413  return (fabs(dphi12) < (aState.dPhiMax*aState.strip_iadd+dphi_incr))? true:false; // +v
414 }
CSCDetId cscDetId() const
Definition: CSCRecHit2D.h:52
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:54
int layer() const
Definition: CSCDetId.h:61
T barePhi() const
Definition: PV3DBase.h:68
unsigned short iChamberType() const
Definition: CSCDetId.h:107
float errorWithinStrip() const
The uncertainty of the estimated position within the strip.
Definition: CSCRecHit2D.h:79
LocalPoint localPosition() const
Definition: CSCRecHit2D.h:50
bool CSCSegAlgoRU::areHitsCloseInR ( const AlgoState aState,
const CSCRecHit2D h1,
const CSCRecHit2D h2 
) const
private

Utility functions.

Definition at line 360 of file CSCSegAlgoRU.cc.

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

Referenced by buildSegments().

360  {
361  float maxWG_width[10] = {0, 0, 4.1, 5.69, 2.49, 5.06, 2.49, 5.06, 1.87, 5.06};
362  CSCDetId id = h1->cscDetId();
363  int iStn = id.iChamberType()-1;
364  //find maxWG_width for ME11 (tilt = 29deg)
365  int wg_num = h2->hitWire();
366  if(iStn == 0 || iStn == 1){
367  if (wg_num == 1){
368  maxWG_width[0] = 9.25;
369  maxWG_width[1] = 9.25;
370  }
371  if (wg_num > 1 && wg_num < 48){
372  maxWG_width[0] = 3.14;
373  maxWG_width[1] = 3.14;
374  }
375  if (wg_num == 48){
376  maxWG_width[0] = 10.75;
377  maxWG_width[1] = 10.75;
378  }
379  }
380  const CSCLayer* l1 = aState.aChamber->layer(h1->cscDetId().layer());
381  GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
382  const CSCLayer* l2 = aState.aChamber->layer(h2->cscDetId().layer());
383  GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
384  //find z to understand the direction
385  float h1z = gp1.z();
386  float h2z = gp2.z();
387  //switch off the IP check for non collisions case
388  if (!aState.doCollisions){
389  h1z = 1;
390  h2z = 1;
391  }
392  if (gp2.perp() > ((gp1.perp() - aState.dRMax*maxWG_width[iStn])*h2z)/h1z && gp2.perp() < ((gp1.perp() + aState.dRMax*maxWG_width[iStn])*h2z)/h1z){
393  return true;
394  }
395  else{
396  return false;
397  }
398 }
CSCDetId cscDetId() const
Definition: CSCRecHit2D.h:52
T perp() const
Definition: PV3DBase.h:72
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:54
int layer() const
Definition: CSCDetId.h:61
T z() const
Definition: PV3DBase.h:64
unsigned short iChamberType() const
Definition: CSCDetId.h:107
short int hitWire() const
L1A.
Definition: CSCRecHit2D.h:62
LocalPoint localPosition() const
Definition: CSCRecHit2D.h:50
void CSCSegAlgoRU::baseline ( AlgoState aState,
int  n_seg_min 
) const
private

(nhits-2)

(nhits-3)

Definition at line 561 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, init_size, gen::k, kLayer(), CSCDetId::layer(), min(), nhits, CSCRecHit2D::nStrips(), CSCSegAlgoRU::AlgoState::proto_segment, CSCDetId::ring(), slope, CSCDetId::station(), and z.

Referenced by buildSegments().

561  {
562  int nhits = aState.proto_segment.size();
563  ChamberHitContainer::const_iterator iRH_worst;
564  //initialise vectors for strip position and error within strip
565  SVector6 sp;
566  SVector6 se;
567  unsigned int init_size = aState.proto_segment.size();
569  buffer.clear();
570  buffer.reserve(init_size);
571  while (buffer.size()< init_size){
572  ChamberHitContainer::iterator min;
573  int min_layer = 10;
574  for(ChamberHitContainer::iterator k = aState.proto_segment.begin(); k != aState.proto_segment.end(); k++){
575  const CSCRecHit2D* iRHk = *k;
576  CSCDetId idRHk = iRHk->cscDetId();
577  int kLayer = idRHk.layer();
578  if(kLayer < min_layer){
579  min_layer = kLayer;
580  min = k;
581  }
582  }
583  buffer.push_back(*min);
584  aState.proto_segment.erase(min);
585  }//while
586 
587  aState.proto_segment.clear();
588  for (ChamberHitContainer::const_iterator cand = buffer.begin(); cand != buffer.end(); cand++) {
589  aState.proto_segment.push_back(*cand);
590  }
591 
592  for(ChamberHitContainer::const_iterator iRH = aState.proto_segment.begin(); iRH != aState.proto_segment.end(); iRH++){
593  const CSCRecHit2D* iRHp = *iRH;
594  CSCDetId idRH = iRHp->cscDetId();
595  int kRing = idRH.ring();
596  int kStation = idRH.station();
597  int kLayer = idRH.layer();
598  // Find the strip containing this hit
599  int centerid = iRHp->nStrips()/2;
600  int centerStrip = iRHp->channels(centerid);
601  float stpos = (*iRHp).positionWithinStrip();
602  se(kLayer-1) = (*iRHp).errorWithinStrip();
603  // Take into account half-strip staggering of layers (ME1/1 has no staggering)
604  if (kStation == 1 && (kRing == 1 || kRing == 4)) sp(kLayer-1) = stpos + centerStrip;
605  else{
606  if (kLayer == 1 || kLayer == 3 || kLayer == 5) sp(kLayer-1) = stpos + centerStrip;
607  if (kLayer == 2 || kLayer == 4 || kLayer == 6) sp(kLayer-1) = stpos - 0.5 + centerStrip;
608  }
609  }
610  float chi2_str;
611  fitX(aState, sp, se, -1, -1, chi2_str);
612 
613  //-----------------------------------------------------
614  // Optimal point rejection method
615  //-----------------------------------------------------
616  float minSum = 1000;
617  int i1b = 0;
618  int i2b = 0;
619  int iworst = -1;
620  int bad_layer = -1;
621  ChamberHitContainer::const_iterator rh_to_be_deleted_1;
622  ChamberHitContainer::const_iterator rh_to_be_deleted_2;
623  if ( (chi2_str) > aState.chi2_str_*aState.chi2D_iadd){
624  for (ChamberHitContainer::const_iterator i1 = aState.proto_segment.begin(); i1 != aState.proto_segment.end();++i1) {
625  ++i1b;
626  const CSCRecHit2D* i1_1 = *i1;
627  CSCDetId idRH1 = i1_1->cscDetId();
628  int z1 = idRH1.layer();
629  i2b = i1b;
630  for (ChamberHitContainer::const_iterator i2 = i1+1; i2 != aState.proto_segment.end(); ++i2) {
631  ++i2b;
632  const CSCRecHit2D* i2_1 = *i2;
633  CSCDetId idRH2 = i2_1->cscDetId();
634  int z2 = idRH2.layer();
635  int irej = 0;
636  for ( ChamberHitContainer::const_iterator ir = aState.proto_segment.begin(); ir != aState.proto_segment.end(); ++ir) {
637  ++irej;
638  if (ir == i1 || ir == i2) continue;
639  float dsum = 0;
640  int hit_nr = 0;
641  const CSCRecHit2D* ir_1 = *ir;
642  CSCDetId idRH = ir_1->cscDetId();
643  int worst_layer = idRH.layer();
644  for (ChamberHitContainer::const_iterator i = aState.proto_segment.begin(); i != aState.proto_segment.end(); ++i) {
645  ++hit_nr;
646  const CSCRecHit2D* i_1 = *i;
647  if (i == i1 || i == i2 || i == ir) continue;
648  float slope = (sp(z2-1)-sp(z1-1))/(z2-z1);
649  float intersept = sp(z1-1) - slope*z1;
650  CSCDetId idRH = i_1->cscDetId();
651  int z = idRH.layer();
652  float di = fabs(sp(z-1) - intersept - slope*z);
653  dsum = dsum + di;
654  }//i
655  if (dsum < minSum){
656  minSum = dsum;
657  bad_layer = worst_layer;
658  iworst = irej;
659  rh_to_be_deleted_1 = ir;
660  }
661  }//ir
662  }//i2
663  }//i1
664  fitX(aState, sp, se, bad_layer, -1, chi2_str);
665  }//if chi2prob<1.0e-4
666 
667  //find worst from n-1 hits
668  int iworst2 = -1;
669  int bad_layer2 = -1;
670  if (iworst > -1 && (nhits-1) > n_seg_min && (chi2_str) > aState.chi2_str_*aState.chi2D_iadd){
671  iworst = -1;
672  float minSum = 1000;
673  int i1b = 0;
674  int i2b = 0;
675  for (ChamberHitContainer::const_iterator i1 = aState.proto_segment.begin(); i1 != aState.proto_segment.end();++i1) {
676  ++i1b;
677  const CSCRecHit2D* i1_1 = *i1;
678  CSCDetId idRH1 = i1_1->cscDetId();
679  int z1 = idRH1.layer();
680  i2b = i1b;
681  for ( ChamberHitContainer::const_iterator i2 = i1+1; i2 != aState.proto_segment.end(); ++i2) {
682  ++i2b;
683  const CSCRecHit2D* i2_1 = *i2;
684  CSCDetId idRH2 = i2_1->cscDetId();
685  int z2 = idRH2.layer();
686  int irej = 0;
687  for ( ChamberHitContainer::const_iterator ir = aState.proto_segment.begin(); ir != aState.proto_segment.end(); ++ir) {
688  ++irej;
689  int irej2 = 0;
690  if (ir == i1 || ir == i2 ) continue;
691  const CSCRecHit2D* ir_1 = *ir;
692  CSCDetId idRH = ir_1->cscDetId();
693  int worst_layer = idRH.layer();
694  for ( ChamberHitContainer::const_iterator ir2 = aState.proto_segment.begin(); ir2 != aState.proto_segment.end(); ++ir2) {
695  ++irej2;
696  if (ir2 == i1 || ir2 == i2 || ir2 ==ir ) continue;
697  float dsum = 0;
698  int hit_nr = 0;
699  const CSCRecHit2D* ir2_1 = *ir2;
700  CSCDetId idRH = ir2_1->cscDetId();
701  int worst_layer2 = idRH.layer();
702  for ( ChamberHitContainer::const_iterator i = aState.proto_segment.begin(); i != aState.proto_segment.end(); ++i) {
703  ++hit_nr;
704  const CSCRecHit2D* i_1 = *i;
705  if (i == i1 || i == i2 || i == ir|| i == ir2 ) continue;
706  float slope = (sp(z2-1)-sp(z1-1))/(z2-z1);
707  float intersept = sp(z1-1) - slope*z1;
708  CSCDetId idRH = i_1->cscDetId();
709  int z = idRH.layer();
710  float di = fabs(sp(z-1) - intersept - slope*z);
711  dsum = dsum + di;
712  }//i
713  if (dsum < minSum){
714  minSum = dsum;
715  iworst2 = irej2;
716  iworst = irej;
717  bad_layer = worst_layer;
718  bad_layer2 = worst_layer2;
719  rh_to_be_deleted_1 = ir;
720  rh_to_be_deleted_2 = ir2;
721  }
722  }//ir2
723  }//ir
724  }//i2
725  }//i1
726  fitX(aState, sp, se, bad_layer ,bad_layer2, chi2_str);
727  }//if prob(n-1)<e-4
728 
729  //----------------------------------
730  //erase bad_hits
731  //----------------------------------
732  if( iworst2-1 >= 0 && iworst2 <= int(aState.proto_segment.size()) ) {
733  aState.proto_segment.erase( rh_to_be_deleted_2);
734  }
735  if( iworst-1 >= 0 && iworst <= int(aState.proto_segment.size()) ){
736  aState.proto_segment.erase(rh_to_be_deleted_1);
737  }
738 }
CSCDetId cscDetId() const
Definition: CSCRecHit2D.h:52
ROOT::Math::SVector< double, 6 > SVector6
Typedefs.
Definition: CSCSegAlgoRU.h:52
static const double slope[3]
int layer() const
Definition: CSCDetId.h:61
int channels(unsigned int i) const
Extracting strip channel numbers comprising the rechit - low.
Definition: CSCRecHit2D.h:55
unsigned int nStrips() const
Definition: CSCRecHit2D.h:56
T min(T a, T b)
Definition: MathUtil.h:58
static const std::string kLayer("layer")
std::vector< const CSCRecHit2D * > ChamberHitContainer
Definition: CSCSegAlgoRU.h:55
int k[5][pyjets_maxn]
int ring() const
Definition: CSCDetId.h:75
float fitX(const AlgoState &aState, SVector6 points, SVector6 errors, int ir, int ir2, float &chi2_str) const
const int init_size
int station() const
Definition: CSCDetId.h:86
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, flagHitsAsUsed(), CSCRecHit2D::hitWire(), mps_fire::i, cuy::ib, isSegmentGood(), gen::k, CSCChamber::layer(), CSCDetId::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().

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

Definition at line 792 of file CSCSegAlgoRU.cc.

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

Referenced by tryAddingHitsToSegment().

792  {
793  // Copy the input CSCSegFit
794  std::unique_ptr<CSCSegFit> oldfit;
795  oldfit.reset(new CSCSegFit( aState.aChamber, aState.proto_segment ));
796  oldfit->fit();
797  ChamberHitContainer oldproto = aState.proto_segment;
798 
799  // May create a new fit
800  bool ok = replaceHit(aState, h, layer);
801  if ( (aState.sfit->chi2() >= oldfit->chi2() ) || !ok ) {
802  // keep original fit
803  aState.proto_segment = oldproto;
804  aState.sfit = std::move(oldfit); // reset to the original input fit
805  }
806 }
bool replaceHit(AlgoState &aState, const CSCRecHit2D *h, int layer) const
std::vector< const CSCRecHit2D * > ChamberHitContainer
Definition: CSCSegAlgoRU.h:55
def move(src, dest)
Definition: eostools.py:510
float CSCSegAlgoRU::fit_r_phi ( const AlgoState aState,
const SVector6 points,
int  layer 
) const
private

Definition at line 542 of file CSCSegAlgoRU.cc.

References delta, mps_fire::i, hiPixelPairStep_cff::points, and slope.

Referenced by isHitNearSegment().

542  {
543  //find R or Phi on the given layer using the given points for the interpolation
544  float Sx = 0;
545  float Sy = 0;
546  float Sxx = 0;
547  float Sxy = 0;
548  for (int i=1;i<7;i++){
549  if (points(i-1)== 0.) continue;
550  Sy = Sy + (points(i-1));
551  Sx = Sx + i;
552  Sxx = Sxx + (i*i);
553  Sxy = Sxy + ((i)*points(i-1));
554  }
555  float delta = 2*Sxx - Sx*Sx;
556  float intercept = (Sxx*Sy - Sx*Sxy)/delta;
557  float slope = (2*Sxy - Sx*Sy)/delta;
558  return (intercept + slope*layer);
559 }
dbl * delta
Definition: mlp_gen.cc:36
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 740 of file CSCSegAlgoRU.cc.

References delta, benchmark_cfg::errors, mps_fire::i, hiPixelPairStep_cff::points, S(), and slope.

Referenced by baseline().

740  {
741  float S = 0;
742  float Sx = 0;
743  float Sy = 0;
744  float Sxx = 0;
745  float Sxy = 0;
746  float sigma2 = 0;
747  for (int i=1;i<7;i++){
748  if (i == ir || i == ir2 || points(i-1) == 0.) continue;
749  sigma2 = errors(i-1)*errors(i-1);
750  float i1 = i - 3.5;
751  S = S + (1/sigma2);
752  Sy = Sy + (points(i-1)/sigma2);
753  Sx = Sx + ((i1)/sigma2);
754  Sxx = Sxx + (i1*i1)/sigma2;
755  Sxy = Sxy + (((i1)*points(i-1))/sigma2);
756  }
757  float delta = S*Sxx - Sx*Sx;
758  float intercept = (Sxx*Sy - Sx*Sxy)/delta;
759  float slope = (S*Sxy - Sx*Sy)/delta;
760  float chi_str = 0;
761  chi2_str = 0;
762  // calculate chi2_str
763  for (int i=1;i<7;i++){
764  if (i == ir || i == ir2 || points(i-1) == 0.) continue;
765  chi_str = (points(i-1) - intercept - slope*(i-3.5))/(errors(i-1));
766  chi2_str = chi2_str + chi_str*chi_str;
767  }
768  return (intercept + slope*0);
769 }
dbl * delta
Definition: mlp_gen.cc:36
static const double slope[3]
double S(const TLorentzVector &, const TLorentzVector &)
Definition: Particle.cc:99
void CSCSegAlgoRU::flagHitsAsUsed ( const AlgoState aState,
const ChamberHitContainer rechitsInChamber,
BoolContainer used 
) const
private

Flag hits on segment as used

Definition at line 508 of file CSCSegAlgoRU.cc.

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

Referenced by buildSegments().

509  {
510  // Flag hits on segment as used
511  ChamberHitContainerCIt ib = rechitsInChamber.begin();
513  for(hi = aState.proto_segment.begin(); hi != aState.proto_segment.end(); ++hi) {
514  for(iu = ib; iu != rechitsInChamber.end(); ++iu) {
515  if(*hi == *iu)
516  used[iu-ib] = true;
517  }
518  }
519 }
std::vector< const CSCRecHit2D * >::const_iterator ChamberHitContainerCIt
Definition: CSCSegAlgoRU.h:56
ib
Definition: cuy.py:660
bool CSCSegAlgoRU::hasHitOnLayer ( const AlgoState aState,
int  layer 
) const
private

Definition at line 771 of file CSCSegAlgoRU.cc.

References CSCSegAlgoRU::AlgoState::proto_segment.

Referenced by tryAddingHitsToSegment().

771  {
772  // Is there is already a hit on this layer?
774  for(it = aState.proto_segment.begin(); it != aState.proto_segment.end(); it++)
775  if ((*it)->cscDetId().layer() == layer)
776  return true;
777  return false;
778 }
std::vector< const CSCRecHit2D * >::const_iterator ChamberHitContainerCIt
Definition: CSCSegAlgoRU.h:56
void CSCSegAlgoRU::increaseProtoSegment ( AlgoState aState,
const CSCRecHit2D h,
int  layer,
int  chi2_factor 
) const
private

Definition at line 808 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().

808  {
809  // Creates a new fit
810  std::unique_ptr<CSCSegFit> oldfit;
811  ChamberHitContainer oldproto = aState.proto_segment;
812  oldfit.reset(new CSCSegFit( aState.aChamber, aState.proto_segment ));
813  oldfit->fit();
814 
815  bool ok = addHit(aState, h, layer);
816  //@@ TEST ON ndof<=0 IS JUST TO ACCEPT nhits=2 CASE??
817  if ( !ok || ( (aState.sfit->ndof() > 0) && (aState.sfit->chi2()/aState.sfit->ndof() >= aState.chi2Max)) ) {
818  // reset to original fit
819  aState.proto_segment = oldproto;
820  aState.sfit = std::move(oldfit);
821  }
822 }
std::vector< const CSCRecHit2D * > ChamberHitContainer
Definition: CSCSegAlgoRU.h:55
bool addHit(AlgoState &aState, const CSCRecHit2D *hit, int layer) const
Utility functions.
def move(src, dest)
Definition: eostools.py:510
bool CSCSegAlgoRU::isHitNearSegment ( const AlgoState aState,
const CSCRecHit2D h 
) const
private

Definition at line 416 of file CSCSegAlgoRU.cc.

References CSCSegAlgoRU::AlgoState::aChamber, CSCRecHit2D::cscDetId(), CSCSegAlgoRU::AlgoState::dPhiIntMax, runTauDisplay::dr, CSCSegAlgoRU::AlgoState::dRIntMax, CSCRecHit2D::errorWithinStrip(), fit_r_phi(), CSCRecHit2D::hitWire(), AnalysisDataFormats_SUSYBSMObjects::hp, CSCDetId::iChamberType(), checklumidiff::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().

416  {
417  // Is hit near segment?
418  // Requires deltaphi and deltaR within ranges specified in parameter set.
419  // Note that to make intuitive cuts on delta(phi) one must work in
420  // phi range (-pi, +pi] not [0, 2pi)
421  float strip_width[10] = {0.003878509, 0.002958185, 0.002327105, 0.00152552, 0.00465421, 0.002327105, 0.00465421, 0.002327105, 0.00465421, 0.002327105};//in rad
422  const CSCLayer* l1 = aState.aChamber->layer((*(aState.proto_segment.begin()))->cscDetId().layer());
423  GlobalPoint gp1 = l1->toGlobal((*(aState.proto_segment.begin()))->localPosition());
424  const CSCLayer* l2 = aState.aChamber->layer((*(aState.proto_segment.begin()+1))->cscDetId().layer());
425  GlobalPoint gp2 = l2->toGlobal((*(aState.proto_segment.begin()+1))->localPosition());
426  float err_stpos_h1 = (*(aState.proto_segment.begin()))->errorWithinStrip();
427  float err_stpos_h2 = (*(aState.proto_segment.begin()+1))->errorWithinStrip();
428  const CSCLayer* l = aState.aChamber->layer(h->cscDetId().layer());
430  float err_stpos_h = h->errorWithinStrip();
431  float hphi = hp.phi(); // in (-pi, +pi]
432  if (hphi < 0.)
433  hphi += 2.*M_PI; // into range [0, 2pi)
434  float sphi = phiAtZ(aState, hp.z()); // in [0, 2*pi)
435  float phidif = sphi-hphi;
436  if (phidif < 0.)
437  phidif += 2.*M_PI; // into range [0, 2pi)
438  if (phidif > M_PI)
439  phidif -= 2.*M_PI; // into range (-pi, pi]
440  SVector6 r_glob;
441  CSCDetId id = h->cscDetId();
442  int iStn = id.iChamberType()-1;
443  float dphi_incr = 0;
444  float pos_str = 1;
445  //increase dPhi cut if the hit is on the edge of the strip
446  float stpos = (*h).positionWithinStrip();
447  bool centr_str = false;
448  if(iStn != 0 && iStn != 1){
449  if (stpos > -0.25 && stpos < 0.25) centr_str = true;
450  }
451  if(err_stpos_h1<0.25*strip_width[iStn] || err_stpos_h2<0.25*strip_width[iStn] || err_stpos_h < 0.25*strip_width[iStn]){
452  dphi_incr = 0.5*strip_width[iStn];
453  }else{
454  if(centr_str) pos_str = 1.3;
455  }
456  r_glob((*(aState.proto_segment.begin()))->cscDetId().layer()-1) = gp1.perp();
457  r_glob((*(aState.proto_segment.begin()+1))->cscDetId().layer()-1) = gp2.perp();
458  float R = hp.perp();
459  int layer = h->cscDetId().layer();
460  float r_interpolated = fit_r_phi(aState, r_glob,layer);
461  float dr = fabs(r_interpolated - R);
462  float maxWG_width[10] = {0, 0, 4.1, 5.69, 2.49, 5.06, 2.49, 5.06, 1.87, 5.06};
463  //find maxWG_width for ME11 (tilt = 29deg)
464  int wg_num = h->hitWire();
465  if(iStn == 0 || iStn == 1){
466  if (wg_num == 1){
467  maxWG_width[0] = 9.25;
468  maxWG_width[1] = 9.25;
469  }
470  if (wg_num > 1 && wg_num < 48){
471  maxWG_width[0] = 3.14;
472  maxWG_width[1] = 3.14;
473  }
474  if (wg_num == 48){
475  maxWG_width[0] = 10.75;
476  maxWG_width[1] = 10.75;
477  }
478  }
479  return (fabs(phidif) < aState.dPhiIntMax*aState.strip_iadd*pos_str+dphi_incr && fabs(dr) < aState.dRIntMax*maxWG_width[iStn])? true:false;
480 }
CSCDetId cscDetId() const
Definition: CSCRecHit2D.h:52
T perp() const
Definition: PV3DBase.h:72
ROOT::Math::SVector< double, 6 > SVector6
Typedefs.
Definition: CSCSegAlgoRU.h:52
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:54
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
int layer() const
Definition: CSCDetId.h:61
susybsm::HSCParticleRefProd hp
Definition: classes.h:27
T z() const
Definition: PV3DBase.h:64
unsigned short iChamberType() const
Definition: CSCDetId.h:107
#define M_PI
float fit_r_phi(const AlgoState &aState, const SVector6 &points, int layer) const
short int hitWire() const
L1A.
Definition: CSCRecHit2D.h:62
float errorWithinStrip() const
The uncertainty of the estimated position within the strip.
Definition: CSCRecHit2D.h:79
LocalPoint localPosition() const
Definition: CSCRecHit2D.h:50
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 495 of file CSCSegAlgoRU.cc.

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

Referenced by buildSegments().

495  {
496  // If the chamber has 20 hits or fewer, require at least 3 hits on segment
497  // If the chamber has >20 hits require at least 4 hits
498  //@@ THESE VALUES SHOULD BECOME PARAMETERS?
499  bool ok = false;
500  unsigned int iadd = ( rechitsInChamber.size() > 20)? 1 : 0;
501  if (aState.windowScale > 1.)
502  iadd = 1;
503  if (aState.proto_segment.size() >= 3+iadd)
504  ok = true;
505  return ok;
506 }
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 482 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().

482  {
483  if ( !aState.sfit ) return 0.;
484  // Returns a phi in [ 0, 2*pi )
485  const CSCLayer* l1 = aState.aChamber->layer((*(aState.proto_segment.begin()))->cscDetId().layer());
486  GlobalPoint gp = l1->toGlobal(aState.sfit->intercept());
487  GlobalVector gv = l1->toGlobal(aState.sfit->localdir());
488  float x = gp.x() + (gv.x()/gv.z())*(z - gp.z());
489  float y = gp.y() + (gv.y()/gv.z())*(z - gp.z());
490  float phi = atan2(y, x);
491  if (phi < 0.f ) phi += 2. * M_PI;
492  return phi ;
493 }
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:54
T y() const
Definition: PV3DBase.h:63
T z() const
Definition: PV3DBase.h:64
double f[11][100]
#define M_PI
T x() const
Definition: PV3DBase.h:62
bool CSCSegAlgoRU::replaceHit ( AlgoState aState,
const CSCRecHit2D h,
int  layer 
) const
private

Definition at line 780 of file CSCSegAlgoRU.cc.

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

Referenced by compareProtoSegment().

780  {
781  // replace a hit from a layer
782  ChamberHitContainer::const_iterator it;
783  for (it = aState.proto_segment.begin(); it != aState.proto_segment.end();) {
784  if ((*it)->cscDetId().layer() == layer)
785  it = aState.proto_segment.erase(it);
786  else
787  ++it;
788  }
789  return addHit(aState, h, layer);
790 }
bool addHit(AlgoState &aState, const CSCRecHit2D *hit, int layer) const
Utility functions.
std::vector<CSCSegment> CSCSegAlgoRU::run ( const CSCChamber aChamber,
const ChamberHitContainer rechits 
)
inline

Here we must implement the algorithm

Definition at line 81 of file CSCSegAlgoRU.h.

References buildSegments().

81 { return buildSegments(aChamber, rechits); }
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 327 of file CSCSegAlgoRU.cc.

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

Referenced by buildSegments().

329  {
330  // Iterate over the layers with hits in the chamber
331  // Skip the layers containing the segment endpoints
332  // Test each hit on the other layers to see if it is near the segment
333  // If it is, see whether there is already a hit on the segment from the same layer
334  // - if so, and there are more than 2 hits on the segment, copy the segment,
335  // replace the old hit with the new hit. If the new segment chi2 is better
336  // then replace the original segment with the new one (by swap)
337  // - if not, copy the segment, add the hit. If the new chi2/dof is still satisfactory
338  // then replace the original segment with the new one (by swap)
340  ChamberHitContainerCIt ie = rechits.end();
341  for (ChamberHitContainerCIt i = ib; i != ie; ++i) {
342  int layer = layerIndex[i-ib];
343  if (hasHitOnLayer(aState, layer) && aState.proto_segment.size() <= 2)continue;
344  if (layerIndex[i-ib] == layerIndex[i1-ib] || layerIndex[i-ib] == layerIndex[i2-ib] || used[i-ib])continue;
345 
346  const CSCRecHit2D* h = *i;
347  if (isHitNearSegment(aState, h)) {
348  // Don't consider alternate hits on layers holding the two starting points
349  if (hasHitOnLayer(aState, layer)) {
350  if (aState.proto_segment.size() <= 2)continue;
351  compareProtoSegment(aState, h, layer);
352  }
353  else{
354  increaseProtoSegment(aState, h, layer, aState.chi2D_iadd);
355  }
356  } // h & seg close
357  } // i
358 }
std::vector< const CSCRecHit2D * >::const_iterator ChamberHitContainerCIt
Definition: CSCSegAlgoRU.h:56
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:660
void CSCSegAlgoRU::updateParameters ( AlgoState aState) const
private

Definition at line 535 of file CSCSegAlgoRU.cc.

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

Referenced by addHit(), and buildSegments().

535  {
536  // Delete input CSCSegFit, create a new one and make the fit
537  // delete sfit;
538  aState.sfit.reset(new CSCSegFit( aState.aChamber, aState.proto_segment ));
539  aState.sfit->fit();
540 }

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().

int CSCSegAlgoRU::minLayersApart
private

Definition at line 163 of file CSCSegAlgoRU.h.

Referenced by CSCSegAlgoRU().

const std::string CSCSegAlgoRU::myName
private

Definition at line 147 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().