36 <<
"--------------------------------------------------------------------\n"
37 <<
"dRMax = " << dRMax <<
'\n'
38 <<
"dPhiMax = " << dPhiMax <<
'\n'
39 <<
"dRIntMax = " << dRIntMax <<
'\n'
40 <<
"dPhiIntMax = " << dPhiIntMax <<
'\n'
41 <<
"chi2Max = " << chi2Max <<
'\n'
42 <<
"wideSeg = " << wideSeg <<
'\n'
43 <<
"minLayersApart = " << minLayersApart << std::endl;
65 int recHits_per_layer[6] = {0,0,0,0,0,0};
68 if (rechits.size()>150){
69 return std::vector<CSCSegment>();
73 for(
unsigned int i = 0;
i < rechits.size();
i++) {
74 recHits_per_layer[rechits[
i]->cscDetId().layer()-1]++;
75 layerIndex[
i] = rechits[
i]->cscDetId().layer();
82 reverse(layerIndex.begin(), layerIndex.end());
83 reverse(rechits.begin(), rechits.end());
86 if (rechits.size() < 2) {
87 return std::vector<CSCSegment>();
112 std::vector<CSCSegment> segments;
120 bool search_disp =
false;
121 int npass = (
wideSeg > 1.)? 3 : 2;
123 for (
int ipass = 0; ipass < npass; ++ipass) {
132 if(used[i1-ib])used_rh++;
136 if(
doCollisions && search_disp &&
int(rechits.size()-used_rh)>2){
148 for(
unsigned int n_seg_min = 6u; n_seg_min > 2u + iadd; --n_seg_min){
150 std::array<BoolContainer, 120> common_used_it = {};
151 for (
unsigned int i = 0;
i < common_used_it.size();
i++) {
152 common_used_it[
i] = common_used;
155 float min_chi[120] = {9999};
157 bool first_proto_segment =
true;
162 if(used[i1-ib] || recHits_per_layer[
int(layerIndex[i1-ib])-1]>25 || (n_seg_min == 3 && used3p[i1-
ib]))
continue;
164 int layer1 = layerIndex[i1-
ib];
169 if(used[i2-ib] || recHits_per_layer[
int(layerIndex[i2-ib])-1]>25 || (n_seg_min == 3 && used3p[i2-
ib]))
continue;
171 int layer2 = layerIndex[i2-
ib];
172 if((
abs(layer2 - layer1) + 1) <
int(n_seg_min))
break;
177 if (!
addHit(h1, layer1))
continue;
178 if (!
addHit(h2, layer2))
continue;
196 if(first_proto_segment){
199 min_chi[0] = sfit_->chi2()/sfit_->ndof();
201 first_proto_segment =
false;
206 min_chi[common_it] = sfit_->chi2()/sfit_->ndof();
209 int iter = common_it;
210 for(iu = ib; iu != ie; ++iu) {
214 for(
int k = 0;
k < iter+1;
k++){
215 if(common_used_it[
k][iu-ib] ==
true){
219 for(ik = ib; ik != ie; ++ik) {
220 if(common_used_it[
k][ik-ib] ==
true){
221 common_used_it[merge_nr][ik-
ib] =
true;
222 common_used_it[
k][ik-
ib] =
false;
226 if(min_chi[
k] < min_chi[merge_nr]){
227 min_chi[merge_nr] = min_chi[
k];
228 best_proto_segment[merge_nr] = best_proto_segment[
k];
229 best_proto_segment[
k].clear();
252 for(
int j = 0;
j < common_it+1;
j++){
255 best_proto_segment[
j].clear();
264 sfit_->localdir(), sfit_->covarianceMatrix(), sfit_->chi2());
266 segments.push_back(
temp);
292 std::vector<CSCSegment>::iterator it =segments.begin();
293 bool good_segs =
false;
294 while(it != segments.end()) {
295 if ((*it).nRecHits() > 3){
313 std::vector<CSCSegment>::iterator it =segments.begin();
314 while(it != segments.end()) {
315 if((*it).nRecHits() == 3){
316 bool found_common =
false;
317 const std::vector<CSCRecHit2D>& theseRH = (*it).specificRecHits();
319 if(used[i1-ib] && used3p[i1-ib]){
323 int sh1layer =
id.
layer();
324 int RH_centerid = sh1->
nStrips()/2;
325 int RH_centerStrip = sh1->
channels(RH_centerid);
327 std::vector<CSCRecHit2D>::const_iterator sh;
328 for(sh = theseRH.begin(); sh != theseRH.end(); ++sh){
332 int shlayer = idRH.
layer();
333 int SegRH_centerid = sh->nStrips()/2;
334 int SegRH_centerStrip = sh->channels(SegRH_centerid);
335 int SegRH_wg = sh->hitWire();
337 if(sh1layer == shlayer && SegRH_centerStrip == RH_centerStrip && SegRH_wg == RH_wg){
339 segments.erase(it,(it+1));
345 if(found_common)
break;
347 if(!found_common)++it;
375 if(layerIndex[i1-ib]<layerIndex[i2-ib]){
376 if (layerIndex[
i-ib] <= layerIndex[i1-ib] || layerIndex[
i-ib] >= layerIndex[i2-ib] ||
i == i1 ||
i == i2 || used[
i-ib]){
377 if (
i == i1 ||
i == i2 || used[
i-ib])
382 if (layerIndex[
i-ib] >= layerIndex[i1-ib] || layerIndex[
i-ib] <= layerIndex[i2-ib] ||
i == i1 ||
i == i2 || used[
i-ib]){
383 if (
i == i1 ||
i == i2 || used[
i-ib])
387 int layer = layerIndex[
i-
ib];
403 float maxWG_width[10] = {0, 0, 4.1, 5.69, 2.49, 5.06, 2.49, 5.06, 1.87, 5.06};
409 if(iStn == 0 || iStn == 1){
411 maxWG_width[0] = 9.25;
412 maxWG_width[1] = 9.25;
414 if (wg_num > 1 && wg_num < 48){
415 maxWG_width[0] = 3.14;
416 maxWG_width[1] = 3.14;
419 maxWG_width[0] = 10.75;
420 maxWG_width[1] = 10.75;
436 if (gp2.
perp() > ((gp1.
perp() -
dRMax*maxWG_width[iStn])*h2z)/h1z && gp2.
perp() < ((gp1.
perp() +
dRMax*maxWG_width[iStn])*h2z)/h1z){
446 float strip_width[10] = {0.003878509, 0.002958185, 0.002327105, 0.00152552, 0.00465421, 0.002327105, 0.00465421, 0.002327105, 0.00465421, 0.002327105};
460 if(err_stpos_h1>0.25*strip_width[iStn] || err_stpos_h2>0.25*strip_width[iStn])dphi_incr = 0.5*strip_width[iStn];
473 float strip_width[10] = {0.003878509, 0.002958185, 0.002327105, 0.00152552, 0.00465421, 0.002327105, 0.00465421, 0.002327105, 0.00465421, 0.002327105};
479 float err_stpos_h1 = (*(
proto_segment.begin()))->errorWithinStrip();
480 float err_stpos_h2 = (*(
proto_segment.begin()+1))->errorWithinStrip();
486 float hphi = hp.
phi();
490 float phidif = sphi-hphi;
503 float stpos = (*h).positionWithinStrip();
504 bool centr_str =
false;
505 if(iStn != 0 && iStn != 1){
506 if (stpos > -0.25 && stpos < 0.25) centr_str =
true;
508 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]){
509 dphi_incr = 0.5*strip_width[iStn];
511 if(centr_str) pos_str = 1.3;
518 float r_interpolated =
fit_r_phi(r_glob,layer);
519 float dr = fabs(r_interpolated - R);
521 float maxWG_width[10] = {0, 0, 4.1, 5.69, 2.49, 5.06, 2.49, 5.06, 1.87, 5.06};
525 if(iStn == 0 || iStn == 1){
527 maxWG_width[0] = 9.25;
528 maxWG_width[1] = 9.25;
530 if (wg_num > 1 && wg_num < 48){
531 maxWG_width[0] = 3.14;
532 maxWG_width[1] = 3.14;
535 maxWG_width[0] = 10.75;
536 maxWG_width[1] = 10.75;
544 if ( !
sfit_ )
return 0.;
551 float x = gp.
x() + (gv.
x()/gv.
z())*(z - gp.
z());
552 float y = gp.
y() + (gv.
y()/gv.
z())*(z - gp.
z());
553 float phi = atan2(y, x);
554 if (phi < 0.
f ) phi += 2. *
M_PI;
566 unsigned int iadd = ( rechitsInChamber.size() > 20)? 1 : 0;
584 for(iu = ib; iu != rechitsInChamber.end(); ++iu) {
597 ChamberHitContainer::const_iterator it;
600 if (((*it)->cscDetId().layer() == layer) && (aHit != (*it)))
622 for (
int i=1;
i<7;
i++){
623 if (points(
i-1)== 0.)
continue;
624 Sy = Sy + (points(
i-1));
627 Sxy = Sxy + ((
i)*points(i-1));
629 float delta = 2*Sxx - Sx*Sx;
630 float intercept = (Sxx*Sy - Sx*Sxy)/delta;
631 float slope = (2*Sxy - Sx*Sy)/delta;
632 return (intercept + slope*layer);
638 ChamberHitContainer::const_iterator iRH_worst;
648 buffer.reserve(init_size);
651 ChamberHitContainer::iterator
min;
657 if(kLayer < min_layer){
662 buffer.push_back(*min);
667 for (ChamberHitContainer::const_iterator cand = buffer.begin(); cand != buffer.end(); cand++) {
674 int kRing = idRH.
ring();
678 int centerid = iRHp->
nStrips()/2;
679 int centerStrip = iRHp->
channels(centerid);
680 float stpos = (*iRHp).positionWithinStrip();
681 se(kLayer-1) = (*iRHp).errorWithinStrip();
683 if (kStation == 1 && (kRing == 1 || kRing == 4)) sp(kLayer-1) = stpos + centerStrip;
685 if (kLayer == 1 || kLayer == 3 || kLayer == 5) sp(kLayer-1) = stpos + centerStrip;
686 if (kLayer == 2 || kLayer == 4 || kLayer == 6) sp(kLayer-1) = stpos - 0.5 + centerStrip;
691 fitX(sp, se, -1, -1, chi2_str);
702 ChamberHitContainer::const_iterator rh_to_be_deleted_1;
703 ChamberHitContainer::const_iterator rh_to_be_deleted_2;
709 int z1 = idRH1.
layer();
711 for (ChamberHitContainer::const_iterator i2 = i1+1; i2 !=
proto_segment.end(); ++i2) {
715 int z2 = idRH2.
layer();
721 if (ir == i1 || ir == i2)
continue;
726 int worst_layer = idRH.
layer();
730 if (
i == i1 ||
i == i2 ||
i == ir)
continue;
731 float slope = (sp(z2-1)-sp(z1-1))/(z2-z1);
732 float intersept = sp(z1-1) - slope*z1;
735 float di = fabs(sp(z-1) - intersept - slope*z);
740 bad_layer = worst_layer;
742 rh_to_be_deleted_1 = ir;
747 fitX(sp, se, bad_layer, -1, chi2_str);
753 if (iworst > -1 && (nhits-1) > n_seg_min && (chi2_str/(nhits-3)) >
chi2_str_*chi2D_iadd){
762 int z1 = idRH1.
layer();
764 for ( ChamberHitContainer::const_iterator i2 = i1+1; i2 !=
proto_segment.end(); ++i2) {
769 int z2 = idRH2.
layer();
776 if (ir == i1 || ir == i2 )
continue;
779 int worst_layer = idRH.
layer();
783 if (ir2 == i1 || ir2 == i2 || ir2 ==ir )
continue;
788 int worst_layer2 = idRH.
layer();
792 if (
i == i1 ||
i == i2 ||
i == ir||
i == ir2 )
continue;
793 float slope = (sp(z2-1)-sp(z1-1))/(z2-z1);
794 float intersept = sp(z1-1) - slope*z1;
797 float di = fabs(sp(z-1) - intersept - slope*z);
804 bad_layer = worst_layer;
805 bad_layer2 = worst_layer2;
806 rh_to_be_deleted_1 = ir;
807 rh_to_be_deleted_2 = ir2;
814 fitX(sp, se, bad_layer ,bad_layer2, chi2_str);
821 if( iworst2-1 >= 0 && iworst2 <=
int(
proto_segment.size()) ) {
839 for (
int i=1;
i<7;
i++){
840 if (
i == ir ||
i == ir2 || points(
i-1) == 0.)
continue;
844 Sy = Sy + (points(
i-1)/sigma2);
845 Sx = Sx + ((i1)/sigma2);
846 Sxx = Sxx + (i1*i1)/sigma2;
847 Sxy = Sxy + (((i1)*points(
i-1))/sigma2);
850 float delta = S*Sxx - Sx*Sx;
851 float intercept = (Sxx*Sy - Sx*Sxy)/delta;
852 float slope = (S*Sxy - Sx*Sy)/delta;
858 for (
int i=1;
i<7;
i++){
859 if (
i == ir ||
i == ir2 || points(
i-1) == 0.)
continue;
860 chi_str = (points(
i-1) - intercept - slope*(
i-3.5))/(
errors(
i-1));
861 chi2_str = chi2_str + chi_str*chi_str;
864 return (intercept + slope*0);
873 if ((*it)->cscDetId().layer() == layer)
882 ChamberHitContainer::const_iterator it;
884 if ((*it)->cscDetId().layer() == layer)
901 if ( (
sfit_->chi2() >= oldfit->chi2() ) || !ok ) {
T getParameter(std::string const &) const
bool areHitsCloseInR(const CSCRecHit2D *h1, const CSCRecHit2D *h2) const
Utility functions.
CSCDetId cscDetId() const
ROOT::Math::SVector< double, 6 > SVector6
Typedefs.
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
void flagHitsAsUsed(const ChamberHitContainer &rechitsInChamber, BoolContainer &used) const
std::vector< const CSCRecHit2D * >::const_iterator ChamberHitContainerCIt
static const double slope[3]
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Geom::Phi< T > phi() const
bool isSegmentGood(const ChamberHitContainer &rechitsInChamber) const
void updateParameters(void)
std::vector< CSCSegment > run(const CSCChamber *aChamber, const ChamberHitContainer &rechits)
bool isHitNearSegment(const CSCRecHit2D *h) const
int channels(unsigned int i) const
Extracting strip channel numbers comprising the rechit - low.
const Surface::PositionType & position() const
The position (origin of the R.F.)
bool replaceHit(const CSCRecHit2D *h, int layer)
ChamberHitContainer proto_segment
float fit_r_phi(SVector6 points, int layer) const
std::vector< CSCSegment > buildSegments(const ChamberHitContainer &rechits)
const CSCChamber * theChamber
unsigned int nStrips() const
Abs< T >::type abs(const T &t)
void tryAddingHitsToSegment(const ChamberHitContainer &rechitsInChamber, const BoolContainer &used, const LayerIndex &layerIndex, const ChamberHitContainerCIt i1, const ChamberHitContainerCIt i2)
const CSCLayer * layer(CSCDetId id) const
Return the layer corresponding to the given id.
void compareProtoSegment(const CSCRecHit2D *h, int layer)
static const std::string kLayer("layer")
CSCSegAlgoRU(const edm::ParameterSet &ps)
Constructor.
std::vector< const CSCRecHit2D * > ChamberHitContainer
std::vector< bool > BoolContainer
std::vector< int > LayerIndex
float phiAtZ(float z) const
unsigned short iChamberType()
bool areHitsCloseInGlobalPhi(const CSCRecHit2D *h1, const CSCRecHit2D *h2) const
float fitX(SVector6 points, SVector6 errors, int ir, int ir2, float &chi2_str)
bool hasHitOnLayer(int layer) const
double S(const TLorentzVector &, const TLorentzVector &)
short int hitWire() const
L1A.
void baseline(int n_seg_min)
float errorWithinStrip() const
The uncertainty of the estimated position within the strip.
LocalPoint localPosition() const
bool addHit(const CSCRecHit2D *hit, int layer)
Utility functions.
std::unique_ptr< CSCSegFit > sfit_
void increaseProtoSegment(const CSCRecHit2D *h, int layer, int chi2_factor)