test
CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CSCSegAlgoRU.cc
Go to the documentation of this file.
1 
8 #include "CSCSegAlgoRU.h"
9 #include "CSCSegFit.h"
15 
16 #include <algorithm>
17 #include <cmath>
18 #include <iostream>
19 #include <string>
20 
22  : CSCSegmentAlgorithm(ps), myName("CSCSegAlgoRU"), sfit_(nullptr) {
23 
24  doCollisions = ps.getParameter<bool>("doCollisions");
25  chi2_str_ = ps.getParameter<double>("chi2_str");
26  chi2Norm_2D_ = ps.getParameter<double>("chi2Norm_2D_");
27  dRMax = ps.getParameter<double>("dRMax");
28  dPhiMax = ps.getParameter<double>("dPhiMax");
29  dRIntMax = ps.getParameter<double>("dRIntMax");
30  dPhiIntMax = ps.getParameter<double>("dPhiIntMax");
31  chi2Max = ps.getParameter<double>("chi2Max");
32  wideSeg = ps.getParameter<double>("wideSeg");
33  minLayersApart = ps.getParameter<int>("minLayersApart");
34 
35  LogDebug("CSC") << myName << " has algorithm cuts set to: \n"
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;
44 
45  //reset the thresholds for non-collision data
46  if(!doCollisions){
47  dRMax = 2.0;
48  dPhiMax = 2*dPhiMax;
49  dRIntMax = 2*dRIntMax;
50  dPhiIntMax = 2*dPhiIntMax;
51  chi2Norm_2D_ = 5*chi2Norm_2D_;
52  chi2_str_ = 100;
53  chi2Max = 2*chi2Max;
54  }
55 }
56 
57 std::vector<CSCSegment> CSCSegAlgoRU::run(const CSCChamber* aChamber, const ChamberHitContainer& rechits){
58  theChamber = aChamber;
59  return buildSegments(rechits);
60 }
61 
62 std::vector<CSCSegment> CSCSegAlgoRU::buildSegments(const ChamberHitContainer& urechits) {
63  ChamberHitContainer rechits = urechits;
64  LayerIndex layerIndex(rechits.size());
65  int recHits_per_layer[6] = {0,0,0,0,0,0};
66 
67  //skip events with high multiplicity of hits
68  if (rechits.size()>150){
69  return std::vector<CSCSegment>();
70  }
71 
72  int iadd = 0;
73  for(unsigned int i = 0; i < rechits.size(); i++) {
74  recHits_per_layer[rechits[i]->cscDetId().layer()-1]++;//count rh per chamber
75  layerIndex[i] = rechits[i]->cscDetId().layer();
76  }
77 
78  double z1 = theChamber->layer(1)->position().z();
79  double z6 = theChamber->layer(6)->position().z();
80 
81  if (std::abs(z1) > std::abs(z6)){
82  reverse(layerIndex.begin(), layerIndex.end());
83  reverse(rechits.begin(), rechits.end());
84  }
85 
86  if (rechits.size() < 2) {
87  return std::vector<CSCSegment>();
88  }
89 
90  // We have at least 2 hits. We intend to try all possible pairs of hits to start
91  // segment building. 'All possible' means each hit lies on different layers in the chamber.
92  // after all same size segs are build we get rid of the overcrossed segments using the chi2 criteria
93  // the hits from the segs that are left are marked as used and are not added to segs in future iterations
94  // the hits from 3p segs are marked as used separately in order to try to assamble them in longer segments
95  // in case there is a second pass
96 
97  // Choose first hit (as close to IP as possible) h1 and second hit
98  // (as far from IP as possible) h2 To do this we iterate over hits
99  // in the chamber by layer - pick two layers. Then we
100  // iterate over hits within each of these layers and pick h1 and h2
101  // these. If they are 'close enough' we build an empty
102  // segment. Then try adding hits to this segment.
103 
104  // Initialize flags that a given hit has been allocated to a segment
105  BoolContainer used(rechits.size(), false);
106  BoolContainer used3p(rechits.size(), false);
107 
108  // This is going to point to fits to hits, and its content will be used to create a CSCSegment
109  sfit_ = 0;
110 
111  // Define buffer for segments we build
112  std::vector<CSCSegment> segments;
113 
114  ChamberHitContainerCIt ib = rechits.begin();
115  ChamberHitContainerCIt ie = rechits.end();
116 
117  // Possibly allow 3 passes, second widening scale factor for cuts, third for segments from displaced vertices
118  windowScale = 1.; // scale factor for cuts
119 
120  bool search_disp = false;
121  int npass = (wideSeg > 1.)? 3 : 2;
122 
123  for (int ipass = 0; ipass < npass; ++ipass) {
124  if(windowScale >1.){
125  iadd = 1;
126  strip_iadd = 2;
127  chi2D_iadd = 2;
128  }
129 
130  int used_rh = 0;
131  for (ChamberHitContainerCIt i1 = ib; i1 != ie; ++i1) {
132  if(used[i1-ib])used_rh++;
133  }
134 
135  //change the tresholds if it's time to look for displaced mu segments
136  if(doCollisions && search_disp && int(rechits.size()-used_rh)>2){//check if there are enough recHits left to build a segment from displaced vertices
137  doCollisions = false;
138  windowScale = 1.; // scale factor for cuts
139  dRMax = 2.0;
140  dPhiMax = 2*dPhiMax;
141  dRIntMax = 2*dRIntMax;
144  chi2_str_ = 100;
145  chi2Max = 2*chi2Max;
146  }
147 
148  for(unsigned int n_seg_min = 6u; n_seg_min > 2u + iadd; --n_seg_min){
149  BoolContainer common_used(rechits.size(),false);
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;
153  }
154  ChamberHitContainer best_proto_segment[120];
155  float min_chi[120] = {9999};
156  int common_it = 0;
157  bool first_proto_segment = true;
158  for (ChamberHitContainerCIt i1 = ib; i1 != ie; ++i1) {
159  bool segok = false;
160 
161  //skip if rh is used and the layer tat has big rh multiplicity(>25RHs)
162  if(used[i1-ib] || recHits_per_layer[int(layerIndex[i1-ib])-1]>25 || (n_seg_min == 3 && used3p[i1-ib])) continue;
163 
164  int layer1 = layerIndex[i1-ib];
165  const CSCRecHit2D* h1 = *i1;
166 
167  for (ChamberHitContainerCIt i2 = ie-1; i2 != i1; --i2) {
168 
169  if(used[i2-ib] || recHits_per_layer[int(layerIndex[i2-ib])-1]>25 || (n_seg_min == 3 && used3p[i2-ib])) continue;
170 
171  int layer2 = layerIndex[i2-ib];
172  if((abs(layer2 - layer1) + 1) < int(n_seg_min)) break;//decrease n_seg_min
173 
174  const CSCRecHit2D* h2 = *i2;
175  if (areHitsCloseInR(h1, h2) && areHitsCloseInGlobalPhi(h1, h2)) {
176  proto_segment.clear();
177  if (!addHit(h1, layer1))continue;
178  if (!addHit(h2, layer2))continue;
179 
180  // Can only add hits if already have a segment
181  if ( sfit_ )tryAddingHitsToSegment(rechits, used, layerIndex, i1, i2);
182 
183  segok = isSegmentGood(rechits);
184  if (segok) {
185  if(proto_segment.size() > n_seg_min){
186  baseline(n_seg_min);
188  }
189 
190  if(sfit_->chi2()/sfit_->ndof() > chi2Norm_2D_*chi2D_iadd || proto_segment.size() < n_seg_min)proto_segment.clear();
191 
192  if (!proto_segment.empty()) {
193  updateParameters();
194 
195  //add same-size overcrossed protosegments to the collection
196  if(first_proto_segment){
197 
198  flagHitsAsUsed(rechits, common_used_it[0]);
199  min_chi[0] = sfit_->chi2()/sfit_->ndof();
200  best_proto_segment[0] = proto_segment;
201  first_proto_segment = false;
202  }else{ //for the rest of found proto_segments
203  common_it++;
204  flagHitsAsUsed(rechits, common_used_it[common_it]);
205 
206  min_chi[common_it] = sfit_->chi2()/sfit_->ndof();
207  best_proto_segment[common_it] = proto_segment;
208  ChamberHitContainerCIt hi, iu, ik;
209  int iter = common_it;
210  for(iu = ib; iu != ie; ++iu) {
211  for(hi = proto_segment.begin(); hi != proto_segment.end(); ++hi) {
212  if(*hi == *iu) {
213  int merge_nr = -1;
214  for(int k = 0; k < iter+1; k++){
215  if(common_used_it[k][iu-ib] == true){
216  if(merge_nr != -1){
217 
218  //merge the k and merge_nr collections of flaged hits into the merge_nr collection and unmark the k collection hits
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;
223  }
224  }
225  //change best_protoseg if min_chi_2 is smaller
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();
230  min_chi[k] = 9999;
231  }
232  common_it--;
233  }
234  else{
235  merge_nr = k;
236  }
237  }//if(common_used[k][iu-ib] == true)
238  }//for k
239  }//if
240  }//for proto_seg
241  }//for rec_hits
242  }//else
243  }//proto seg not empty
244  }
245  } // h1 & h2 close
246  if (segok)
247  break;
248  } // i2
249  } // i1
250 
251  //add the reconstructed segments
252  for(int j = 0;j < common_it+1; j++){
253 
254  proto_segment = best_proto_segment[j];
255  best_proto_segment[j].clear();
256 
257  //SKIP empty proto-segments
258  if(proto_segment.size() == 0) continue;
259 
260  updateParameters();
261 
262  // Create an actual CSCSegment - retrieve all info from the fit
263  CSCSegment temp(sfit_->hits(), sfit_->intercept(),
264  sfit_->localdir(), sfit_->covarianceMatrix(), sfit_->chi2());
265  sfit_ = 0;
266  segments.push_back(temp);
267  //if the segment has 3 hits flag them as used in a particular way
268  if(proto_segment.size() == 3){
269  flagHitsAsUsed(rechits, used3p);
270  }
271  else{
272  flagHitsAsUsed(rechits, used);
273  }
274  proto_segment.clear();
275  }
276  }//for n_seg_min
277 
278  if(!doCollisions && search_disp){
279 
280  //reset params and flags for the next chamber
281  search_disp = false;
282  doCollisions = true;
283  dRMax = 2.0;
284  dPhiMax = dPhiMax/2;
285  dRIntMax = dRIntMax/2;
288  chi2_str_ = 100;
289  chi2Max = chi2Max/2;
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 
302  if (good_segs) { // 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(!doCollisions && !search_disp) break;
309  windowScale = wideSeg;
310  } // ipass
311 
312  //get rid of enchansed 3p segments
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();
318  for (ChamberHitContainerCIt i1 = ib; i1 != ie; ++i1) {
319  if(used[i1-ib] && used3p[i1-ib]){
320  const CSCRecHit2D* sh1 = *i1;
321 
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 
331  //find segment hit coord
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();
336 
337  if(sh1layer == shlayer && SegRH_centerStrip == RH_centerStrip && SegRH_wg == RH_wg){
338  //remove the enchansed 3p segment
339  segments.erase(it,(it+1));
340  found_common = true;
341  break;
342  }
343  }//theserh
344  }
345  if(found_common) break;//current seg has already been erased
346  }//camber hits
347  if(!found_common)++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
356 
358  const BoolContainer& used, const LayerIndex& layerIndex,
359  const ChamberHitContainerCIt i1, const ChamberHitContainerCIt i2) {
360 
361  // Iterate over the layers with hits in the chamber
362  // Skip the layers containing the segment endpoints
363  // Test each hit on the other layers to see if it is near the segment
364  // If it is, see whether there is already a hit on the segment from the same layer
365  // - if so, and there are more than 2 hits on the segment, copy the segment,
366  // replace the old hit with the new hit. If the new segment chi2 is better
367  // then replace the original segment with the new one (by swap)
368  // - if not, copy the segment, add the hit. If the new chi2/dof is still satisfactory
369  // then replace the original segment with the new one (by swap)
370 
371  ChamberHitContainerCIt ib = rechits.begin();
372  ChamberHitContainerCIt ie = rechits.end();
373  for (ChamberHitContainerCIt i = ib; i != ie; ++i) {
374 
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])
378  continue;
379  }
380  }
381  else{
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])
384  continue;
385  }
386  }
387  int layer = layerIndex[i-ib];
388  const CSCRecHit2D* h = *i;
389  if (isHitNearSegment(h)) {
390  // Don't consider alternate hits on layers holding the two starting points
391  if (hasHitOnLayer(layer)) {
392  if (proto_segment.size() <= 2)continue;
393  compareProtoSegment(h, layer);
394  }
395  else{
396  increaseProtoSegment(h, layer, chi2D_iadd);
397  }
398  } // h & seg close
399  } // i
400 }
401 
402 bool CSCSegAlgoRU::areHitsCloseInR(const CSCRecHit2D* h1, const CSCRecHit2D* h2) const {
403  float maxWG_width[10] = {0, 0, 4.1, 5.69, 2.49, 5.06, 2.49, 5.06, 1.87, 5.06};
404  CSCDetId id = h1->cscDetId();
405  int iStn = id.iChamberType()-1;
406 
407  //find maxWG_width for ME11 (tilt = 29deg)
408  int wg_num = h2->hitWire();
409  if(iStn == 0 || iStn == 1){
410  if (wg_num == 1){
411  maxWG_width[0] = 9.25;
412  maxWG_width[1] = 9.25;
413  }
414  if (wg_num > 1 && wg_num < 48){
415  maxWG_width[0] = 3.14;
416  maxWG_width[1] = 3.14;
417  }
418  if (wg_num == 48){
419  maxWG_width[0] = 10.75;
420  maxWG_width[1] = 10.75;
421  }
422  }
423  const CSCLayer* l1 = theChamber->layer(h1->cscDetId().layer());
424  GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
425  const CSCLayer* l2 = theChamber->layer(h2->cscDetId().layer());
426  GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
427  //find z to understand the direction
428  float h1z = gp1.z();
429  float h2z = gp2.z();
430 
431  //switch off the IP check for non collisions case
432  if (!doCollisions){
433  h1z = 1;
434  h2z = 1;
435  }
436  if (gp2.perp() > ((gp1.perp() - dRMax*maxWG_width[iStn])*h2z)/h1z && gp2.perp() < ((gp1.perp() + dRMax*maxWG_width[iStn])*h2z)/h1z){
437  return true;
438  }
439  else{
440  return false;
441  }
442 }
443 
445 
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};//in rad
447 
448  const CSCLayer* l1 = theChamber->layer(h1->cscDetId().layer());
449  GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
450  const CSCLayer* l2 = theChamber->layer(h2->cscDetId().layer());
451  GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
452  float err_stpos_h1 = h1->errorWithinStrip();
453  float err_stpos_h2 = h2->errorWithinStrip();
454 
455  CSCDetId id = h1->cscDetId();
456  int iStn = id.iChamberType()-1;
457 
458 
459  float dphi_incr = 0;
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];
461  float dphi12 = deltaPhi(gp1.barePhi(),gp2.barePhi());
462 
463  return (fabs(dphi12) < (dPhiMax*strip_iadd+dphi_incr))? true:false; // +v
464 }
465 
467 
468  // Is hit near segment?
469  // Requires deltaphi and deltaR within ranges specified in parameter set.
470  // Note that to make intuitive cuts on delta(phi) one must work in
471  // phi range (-pi, +pi] not [0, 2pi)
472 
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};//in rad
474 
475  const CSCLayer* l1 = theChamber->layer((*(proto_segment.begin()))->cscDetId().layer());
476  GlobalPoint gp1 = l1->toGlobal((*(proto_segment.begin()))->localPosition());
477  const CSCLayer* l2 = theChamber->layer((*(proto_segment.begin()+1))->cscDetId().layer());
478  GlobalPoint gp2 = l2->toGlobal((*(proto_segment.begin()+1))->localPosition());
479  float err_stpos_h1 = (*(proto_segment.begin()))->errorWithinStrip();
480  float err_stpos_h2 = (*(proto_segment.begin()+1))->errorWithinStrip();
481 
482  const CSCLayer* l = theChamber->layer(h->cscDetId().layer());
483  GlobalPoint hp = l->toGlobal(h->localPosition());
484  float err_stpos_h = h->errorWithinStrip();
485 
486  float hphi = hp.phi(); // in (-pi, +pi]
487  if (hphi < 0.)
488  hphi += 2.*M_PI; // into range [0, 2pi)
489  float sphi = phiAtZ(hp.z()); // in [0, 2*pi)
490  float phidif = sphi-hphi;
491  if (phidif < 0.)
492  phidif += 2.*M_PI; // into range [0, 2pi)
493  if (phidif > M_PI)
494  phidif -= 2.*M_PI; // into range (-pi, pi]
495 
496  SVector6 r_glob;
497  CSCDetId id = h->cscDetId();
498  int iStn = id.iChamberType()-1;
499  float dphi_incr = 0;
500  float pos_str = 1;
501 
502  //increase dPhi cut if the hit is on the edge of the strip
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;
507  }
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];
510  }else{
511  if(centr_str) pos_str = 1.3;
512  }
513  r_glob((*(proto_segment.begin()))->cscDetId().layer()-1) = gp1.perp();
514  r_glob((*(proto_segment.begin()+1))->cscDetId().layer()-1) = gp2.perp();
515 
516  float R = hp.perp();
517  int layer = h->cscDetId().layer();
518  float r_interpolated = fit_r_phi(r_glob,layer);
519  float dr = fabs(r_interpolated - R);
520 
521  float maxWG_width[10] = {0, 0, 4.1, 5.69, 2.49, 5.06, 2.49, 5.06, 1.87, 5.06};
522 
523  //find maxWG_width for ME11 (tilt = 29deg)
524  int wg_num = h->hitWire();
525  if(iStn == 0 || iStn == 1){
526  if (wg_num == 1){
527  maxWG_width[0] = 9.25;
528  maxWG_width[1] = 9.25;
529  }
530  if (wg_num > 1 && wg_num < 48){
531  maxWG_width[0] = 3.14;
532  maxWG_width[1] = 3.14;
533  }
534  if (wg_num == 48){
535  maxWG_width[0] = 10.75;
536  maxWG_width[1] = 10.75;
537  }
538  }
539  return (fabs(phidif) < dPhiIntMax*strip_iadd*pos_str+dphi_incr && fabs(dr) < dRIntMax*maxWG_width[iStn])? true:false;
540 }
541 
542 float CSCSegAlgoRU::phiAtZ(float z) const {
543 
544  if ( !sfit_ ) return 0.;
545 
546  // Returns a phi in [ 0, 2*pi )
547  const CSCLayer* l1 = theChamber->layer((*(proto_segment.begin()))->cscDetId().layer());
548  GlobalPoint gp = l1->toGlobal(sfit_->intercept());
549  GlobalVector gv = l1->toGlobal(sfit_->localdir());
550 
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;
555 
556  return phi ;
557 }
558 
559 bool CSCSegAlgoRU::isSegmentGood(const ChamberHitContainer& rechitsInChamber) const {
560 
561  // If the chamber has 20 hits or fewer, require at least 3 hits on segment
562  // If the chamber has >20 hits require at least 4 hits
563  //@@ THESE VALUES SHOULD BECOME PARAMETERS?
564  bool ok = false;
565 
566  unsigned int iadd = ( rechitsInChamber.size() > 20)? 1 : 0;
567 
568  if (windowScale > 1.)
569  iadd = 1;
570 
571  if (proto_segment.size() >= 3+iadd)
572  ok = true;
573  return ok;
574 }
575 
576 void CSCSegAlgoRU::flagHitsAsUsed(const ChamberHitContainer& rechitsInChamber,
577  BoolContainer& used ) const {
578 
579  // Flag hits on segment as used
580  ChamberHitContainerCIt ib = rechitsInChamber.begin();
581  ChamberHitContainerCIt hi, iu;
582 
583  for(hi = proto_segment.begin(); hi != proto_segment.end(); ++hi) {
584  for(iu = ib; iu != rechitsInChamber.end(); ++iu) {
585  if(*hi == *iu)
586  used[iu-ib] = true;
587  }
588  }
589 }
590 
591 bool CSCSegAlgoRU::addHit(const CSCRecHit2D* aHit, int layer) {
592 
593  // Return true if hit was added successfully
594  // (and then parameters are updated).
595  // Return false if there is already a hit on the same layer, or insert failed.
596 
597  ChamberHitContainer::const_iterator it;
598 
599  for(it = proto_segment.begin(); it != proto_segment.end(); it++)
600  if (((*it)->cscDetId().layer() == layer) && (aHit != (*it)))
601  return false;
602 
603  proto_segment.push_back(aHit);
604  // make a fit
606  return true;
607 }
608 
610  // update the current CSCSegFit one and make the fit
611  sfit_.reset(new CSCSegFit( theChamber, proto_segment ));
612  sfit_->fit();
613 }
614 
615 float CSCSegAlgoRU::fit_r_phi(SVector6 points, int layer) const{
616  //find R or Phi on the given layer using the given points for the interpolation
617  float Sx = 0;
618  float Sy = 0;
619  float Sxx = 0;
620  float Sxy = 0;
621 
622  for (int i=1;i<7;i++){
623  if (points(i-1)== 0.) continue;
624  Sy = Sy + (points(i-1));
625  Sx = Sx + i;
626  Sxx = Sxx + (i*i);
627  Sxy = Sxy + ((i)*points(i-1));
628  }
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);
633 }
634 
635 void CSCSegAlgoRU::baseline(int n_seg_min){
636 
637  int nhits = proto_segment.size();
638  ChamberHitContainer::const_iterator iRH_worst;
639 
640  //initialise vectors for strip position and error within strip
641  SVector6 sp;
642  SVector6 se;
643 
644  unsigned int init_size = proto_segment.size();
645 
646  ChamberHitContainer buffer;
647  buffer.clear();
648  buffer.reserve(init_size);
649 
650  while (buffer.size()< init_size){
651  ChamberHitContainer::iterator min;
652  int min_layer = 10;
653  for(ChamberHitContainer::iterator k = proto_segment.begin(); k != proto_segment.end(); k++){
654  const CSCRecHit2D* iRHk = *k;
655  CSCDetId idRHk = iRHk->cscDetId();
656  int kLayer = idRHk.layer();
657  if(kLayer < min_layer){
658  min_layer = kLayer;
659  min = k;
660  }
661  }
662  buffer.push_back(*min);
663  proto_segment.erase(min);
664  }//while
665  proto_segment.clear();
666 
667  for (ChamberHitContainer::const_iterator cand = buffer.begin(); cand != buffer.end(); cand++) {
668  proto_segment.push_back(*cand);
669  }
670 
671  for(ChamberHitContainer::const_iterator iRH = proto_segment.begin(); iRH != proto_segment.end(); iRH++){
672  const CSCRecHit2D* iRHp = *iRH;
673  CSCDetId idRH = iRHp->cscDetId();
674  int kRing = idRH.ring();
675  int kStation = idRH.station();
676  int kLayer = idRH.layer();
677  // Find the strip containing this hit
678  int centerid = iRHp->nStrips()/2;
679  int centerStrip = iRHp->channels(centerid);
680  float stpos = (*iRHp).positionWithinStrip();
681  se(kLayer-1) = (*iRHp).errorWithinStrip();
682  // Take into account half-strip staggering of layers (ME1/1 has no staggering)
683  if (kStation == 1 && (kRing == 1 || kRing == 4)) sp(kLayer-1) = stpos + centerStrip;
684  else{
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;
687  }
688  }
689 
690  float chi2_str;
691  fitX(sp, se, -1, -1, chi2_str);
692 
693  //-----------------------------------------------------
694  // Optimal point rejection method
695  //-----------------------------------------------------
696 
697  float minSum = 1000;
698  int i1b = 0;
699  int i2b = 0;
700  int iworst = -1;
701  int bad_layer = -1;
702  ChamberHitContainer::const_iterator rh_to_be_deleted_1;
703  ChamberHitContainer::const_iterator rh_to_be_deleted_2;
704  if ( nhits > n_seg_min && (chi2_str/(nhits-2)) > chi2_str_*chi2D_iadd){
705  for (ChamberHitContainer::const_iterator i1 = proto_segment.begin(); i1 != proto_segment.end();++i1) {
706  ++i1b;
707  const CSCRecHit2D* i1_1 = *i1;
708  CSCDetId idRH1 = i1_1->cscDetId();
709  int z1 = idRH1.layer();
710  i2b = i1b;
711  for (ChamberHitContainer::const_iterator i2 = i1+1; i2 != proto_segment.end(); ++i2) {
712  ++i2b;
713  const CSCRecHit2D* i2_1 = *i2;
714  CSCDetId idRH2 = i2_1->cscDetId();
715  int z2 = idRH2.layer();
716  int irej = 0;
717 
718  for ( ChamberHitContainer::const_iterator ir = proto_segment.begin(); ir != proto_segment.end(); ++ir) {
719  ++irej;
720 
721  if (ir == i1 || ir == i2) continue;
722  float dsum = 0;
723  int hit_nr = 0;
724  const CSCRecHit2D* ir_1 = *ir;
725  CSCDetId idRH = ir_1->cscDetId();
726  int worst_layer = idRH.layer();
727  for (ChamberHitContainer::const_iterator i = proto_segment.begin(); i != proto_segment.end(); ++i) {
728  ++hit_nr;
729  const CSCRecHit2D* i_1 = *i;
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;
733  CSCDetId idRH = i_1->cscDetId();
734  int z = idRH.layer();
735  float di = fabs(sp(z-1) - intersept - slope*z);
736  dsum = dsum + di;
737  }//i
738  if (dsum < minSum){
739  minSum = dsum;
740  bad_layer = worst_layer;
741  iworst = irej;
742  rh_to_be_deleted_1 = ir;
743  }
744  }//ir
745  }//i2
746  }//i1
747  fitX(sp, se, bad_layer, -1, chi2_str);
748  }//if chi2prob<1.0e-4
749 
750  //find worst from n-1 hits
751  int iworst2 = -1;
752  int bad_layer2 = -1;
753  if (iworst > -1 && (nhits-1) > n_seg_min && (chi2_str/(nhits-3)) > chi2_str_*chi2D_iadd){
754  iworst = -1;
755  float minSum = 1000;
756  int i1b = 0;
757  int i2b = 0;
758  for (ChamberHitContainer::const_iterator i1 = proto_segment.begin(); i1 != proto_segment.end();++i1) {
759  ++i1b;
760  const CSCRecHit2D* i1_1 = *i1;
761  CSCDetId idRH1 = i1_1->cscDetId();
762  int z1 = idRH1.layer();
763  i2b = i1b;
764  for ( ChamberHitContainer::const_iterator i2 = i1+1; i2 != proto_segment.end(); ++i2) {
765  ++i2b;
766  const CSCRecHit2D* i2_1 = *i2;
767 
768  CSCDetId idRH2 = i2_1->cscDetId();
769  int z2 = idRH2.layer();
770  int irej = 0;
771 
772  for ( ChamberHitContainer::const_iterator ir = proto_segment.begin(); ir != proto_segment.end(); ++ir) {
773 
774  ++irej;
775  int irej2 = 0;
776  if (ir == i1 || ir == i2 ) continue;
777  const CSCRecHit2D* ir_1 = *ir;
778  CSCDetId idRH = ir_1->cscDetId();
779  int worst_layer = idRH.layer();
780  for ( ChamberHitContainer::const_iterator ir2 = proto_segment.begin(); ir2 != proto_segment.end(); ++ir2) {
781 
782  ++irej2;
783  if (ir2 == i1 || ir2 == i2 || ir2 ==ir ) continue;
784  float dsum = 0;
785  int hit_nr = 0;
786  const CSCRecHit2D* ir2_1 = *ir2;
787  CSCDetId idRH = ir2_1->cscDetId();
788  int worst_layer2 = idRH.layer();
789  for ( ChamberHitContainer::const_iterator i = proto_segment.begin(); i != proto_segment.end(); ++i) {
790  ++hit_nr;
791  const CSCRecHit2D* i_1 = *i;
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;
795  CSCDetId idRH = i_1->cscDetId();
796  int z = idRH.layer();
797  float di = fabs(sp(z-1) - intersept - slope*z);
798  dsum = dsum + di;
799  }//i
800  if (dsum < minSum){
801  minSum = dsum;
802  iworst2 = irej2;
803  iworst = irej;
804  bad_layer = worst_layer;
805  bad_layer2 = worst_layer2;
806  rh_to_be_deleted_1 = ir;
807  rh_to_be_deleted_2 = ir2;
808  }
809  }//ir2
810  }//ir
811  }//i2
812  }//i1
813 
814  fitX(sp, se, bad_layer ,bad_layer2, chi2_str);
815  }//if prob(n-1)<e-4
816 
817  //----------------------------------
818  //erase bad_hits
819  //----------------------------------
820 
821  if( iworst2-1 >= 0 && iworst2 <= int(proto_segment.size()) ) {
822  proto_segment.erase( rh_to_be_deleted_2);
823  }
824 
825  if( iworst-1 >= 0 && iworst <= int(proto_segment.size()) ){
826  proto_segment.erase(rh_to_be_deleted_1);
827  }
828 }
829 
830 float CSCSegAlgoRU::fitX(SVector6 points, SVector6 errors, int ir, int ir2, float &chi2_str){
831 
832  float S = 0;
833  float Sx = 0;
834  float Sy = 0;
835  float Sxx = 0;
836  float Sxy = 0;
837  float sigma2 = 0;
838 
839  for (int i=1;i<7;i++){
840  if (i == ir || i == ir2 || points(i-1) == 0.) continue;
841  sigma2 = errors(i-1)*errors(i-1);
842  float i1 = i - 3.5;
843  S = S + (1/sigma2);
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);
848  }
849 
850  float delta = S*Sxx - Sx*Sx;
851  float intercept = (Sxx*Sy - Sx*Sxy)/delta;
852  float slope = (S*Sxy - Sx*Sy)/delta;
853 
854  float chi_str = 0;
855  chi2_str = 0;
856 
857  // calculate chi2_str
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;
862  }
863 
864  return (intercept + slope*0);
865 }
866 
867 bool CSCSegAlgoRU::hasHitOnLayer(int layer) const {
868 
869  // Is there is already a hit on this layer?
871 
872  for(it = proto_segment.begin(); it != proto_segment.end(); it++)
873  if ((*it)->cscDetId().layer() == layer)
874  return true;
875 
876  return false;
877 }
878 
879 bool CSCSegAlgoRU::replaceHit(const CSCRecHit2D* h, int layer) {
880 
881  // replace a hit from a layer
882  ChamberHitContainer::const_iterator it;
883  for (it = proto_segment.begin(); it != proto_segment.end();) {
884  if ((*it)->cscDetId().layer() == layer)
885  it = proto_segment.erase(it);
886  else
887  ++it;
888  }
889 
890  return addHit(h, layer);
891 }
892 
894  // Copy the input CSCSegFit
895  std::unique_ptr<CSCSegFit> oldfit = std::make_unique<CSCSegFit>(theChamber, proto_segment);
896  oldfit->fit();
897 
898  // May create a new fit
899  bool ok = replaceHit(h, layer);
900 
901  if ( ( sfit_->chi2() >= oldfit->chi2() ) || !ok ) {
902  sfit_ = std::move(oldfit); // reset to the original input fit
903  }
904 }
905 
906 void CSCSegAlgoRU::increaseProtoSegment(const CSCRecHit2D* h, int layer, int chi2_factor) {
907  // Creates a new fit
908  std::unique_ptr<CSCSegFit> oldfit = std::make_unique<CSCSegFit>(theChamber, proto_segment);
909 
910  bool ok = addHit(h, layer);
911  //@@ TEST ON ndof<=0 IS JUST TO ACCEPT nhits=2 CASE??
912  if ( !ok || ( (sfit_->ndof() > 0) && (sfit_->chi2()/sfit_->ndof() >= chi2Max)) ) {
913  sfit_ = std::move(oldfit);
914  }
915 }
#define LogDebug(id)
dbl * delta
Definition: mlp_gen.cc:36
T getParameter(std::string const &) const
int i
Definition: DBlmapReader.cc:9
bool areHitsCloseInR(const CSCRecHit2D *h1, const CSCRecHit2D *h2) const
Utility functions.
CSCDetId cscDetId() const
Definition: CSCRecHit2D.h:52
float windowScale
Definition: CSCSegAlgoRU.h:134
int ib
Definition: cuy.py:660
T perp() const
Definition: PV3DBase.h:72
ROOT::Math::SVector< double, 6 > SVector6
Typedefs.
Definition: CSCSegAlgoRU.h:52
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
Definition: CSCSegAlgoRU.h:56
static const double slope[3]
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
T y() const
Definition: PV3DBase.h:63
bool isSegmentGood(const ChamberHitContainer &rechitsInChamber) const
const std::string myName
Definition: CSCSegAlgoRU.h:128
void updateParameters(void)
#define nullptr
int layer() const
Definition: CSCDetId.h:61
float dPhiIntMax
Definition: CSCSegAlgoRU.h:141
float float float z
std::vector< CSCSegment > run(const CSCChamber *aChamber, const ChamberHitContainer &rechits)
Definition: CSCSegAlgoRU.cc:57
T barePhi() const
Definition: PV3DBase.h:68
bool isHitNearSegment(const CSCRecHit2D *h) 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
bool replaceHit(const CSCRecHit2D *h, int layer)
ChamberHitContainer proto_segment
Definition: CSCSegAlgoRU.h:127
float fit_r_phi(SVector6 points, int layer) const
susybsm::HSCParticleRefProd hp
Definition: classes.h:27
std::vector< CSCSegment > buildSegments(const ChamberHitContainer &rechits)
Definition: CSCSegAlgoRU.cc:62
const CSCChamber * theChamber
Definition: CSCSegAlgoRU.h:126
T z() const
Definition: PV3DBase.h:64
def move
Definition: eostools.py:510
unsigned int nStrips() const
Definition: CSCRecHit2D.h:56
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void tryAddingHitsToSegment(const ChamberHitContainer &rechitsInChamber, const BoolContainer &used, const LayerIndex &layerIndex, const ChamberHitContainerCIt i1, const ChamberHitContainerCIt i2)
int j
Definition: DBlmapReader.cc:9
const CSCLayer * layer(CSCDetId id) const
Return the layer corresponding to the given id.
Definition: CSCChamber.cc:39
void compareProtoSegment(const CSCRecHit2D *h, int layer)
double f[11][100]
T min(T a, T b)
Definition: MathUtil.h:58
static const std::string kLayer("layer")
CSCSegAlgoRU(const edm::ParameterSet &ps)
Constructor.
Definition: CSCSegAlgoRU.cc:21
std::vector< const CSCRecHit2D * > ChamberHitContainer
Definition: CSCSegAlgoRU.h:55
std::vector< bool > BoolContainer
Definition: CSCSegAlgoRU.h:59
std::vector< int > LayerIndex
Definition: CSCSegAlgoRU.h:54
#define M_PI
float phiAtZ(float z) const
unsigned short iChamberType()
Definition: CSCDetId.h:107
bool areHitsCloseInGlobalPhi(const CSCRecHit2D *h1, const CSCRecHit2D *h2) const
int ring() const
Definition: CSCDetId.h:75
float chi2Norm_2D_
Definition: CSCSegAlgoRU.h:144
float fitX(SVector6 points, SVector6 errors, int ir, int ir2, float &chi2_str)
bool hasHitOnLayer(int layer) const
double S(const TLorentzVector &, const TLorentzVector &)
Definition: Particle.cc:99
short int hitWire() const
L1A.
Definition: CSCRecHit2D.h:62
void baseline(int n_seg_min)
const int init_size
float errorWithinStrip() const
The uncertainty of the estimated position within the strip.
Definition: CSCRecHit2D.h:79
int station() const
Definition: CSCDetId.h:86
Definition: sp.h:21
LocalPoint localPosition() const
Definition: CSCRecHit2D.h:50
bool addHit(const CSCRecHit2D *hit, int layer)
Utility functions.
std::unique_ptr< CSCSegFit > sfit_
Definition: CSCSegAlgoRU.h:149
T x() const
Definition: PV3DBase.h:62
void increaseProtoSegment(const CSCRecHit2D *h, int layer, int chi2_factor)