Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:39 

0001 /**
0002  * \file CSCSegAlgRU.cc
0003  *
0004  * \authors V.Palichik & N.Voytishin
0005  * \some functions and structure taken from SK algo by M.Sani and SegFit class by T.Cox
0006  */
0007 
0008 #include "CSCSegAlgoRU.h"
0009 #include "CSCSegFit.h"
0010 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0011 #include "DataFormats/Math/interface/deltaPhi.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0015 #include <algorithm>
0016 #include <cmath>
0017 #include <iostream>
0018 #include <memory>
0019 
0020 #include <string>
0021 
0022 CSCSegAlgoRU::CSCSegAlgoRU(const edm::ParameterSet& ps) : CSCSegmentAlgorithm(ps), myName("CSCSegAlgoRU") {
0023   doCollisions = ps.getParameter<bool>("doCollisions");
0024   chi2_str_ = ps.getParameter<double>("chi2_str");
0025   chi2Norm_2D_ = ps.getParameter<double>("chi2Norm_2D_");
0026   dRMax = ps.getParameter<double>("dRMax");
0027   dPhiMax = ps.getParameter<double>("dPhiMax");
0028   dRIntMax = ps.getParameter<double>("dRIntMax");
0029   dPhiIntMax = ps.getParameter<double>("dPhiIntMax");
0030   chi2Max = ps.getParameter<double>("chi2Max");
0031   wideSeg = ps.getParameter<double>("wideSeg");
0032   minLayersApart = ps.getParameter<int>("minLayersApart");
0033 
0034   LogDebug("CSC") << myName << " has algorithm cuts set to: \n"
0035                   << "--------------------------------------------------------------------\n"
0036                   << "dRMax = " << dRMax << '\n'
0037                   << "dPhiMax = " << dPhiMax << '\n'
0038                   << "dRIntMax = " << dRIntMax << '\n'
0039                   << "dPhiIntMax = " << dPhiIntMax << '\n'
0040                   << "chi2Max = " << chi2Max << '\n'
0041                   << "wideSeg = " << wideSeg << '\n'
0042                   << "minLayersApart = " << minLayersApart << std::endl;
0043 
0044   //reset the thresholds for non-collision data
0045   if (!doCollisions) {
0046     dRMax = 2.0;
0047     dPhiMax = 2 * dPhiMax;
0048     dRIntMax = 2 * dRIntMax;
0049     dPhiIntMax = 2 * dPhiIntMax;
0050     chi2Norm_2D_ = 5 * chi2Norm_2D_;
0051     chi2_str_ = 100;
0052     chi2Max = 2 * chi2Max;
0053   }
0054 }
0055 
0056 std::vector<CSCSegment> CSCSegAlgoRU::buildSegments(const CSCChamber* aChamber,
0057                                                     const ChamberHitContainer& urechits) const {
0058   ChamberHitContainer rechits = urechits;
0059   LayerIndex layerIndex(rechits.size());
0060   int recHits_per_layer[6] = {0, 0, 0, 0, 0, 0};
0061   //skip events with high multiplicity of hits
0062   if (rechits.size() > 150) {
0063     return std::vector<CSCSegment>();
0064   }
0065   int iadd = 0;
0066   for (unsigned int i = 0; i < rechits.size(); i++) {
0067     recHits_per_layer[rechits[i]->cscDetId().layer() - 1]++;  //count rh per chamber
0068     layerIndex[i] = rechits[i]->cscDetId().layer();
0069   }
0070   double z1 = aChamber->layer(1)->position().z();
0071   double z6 = aChamber->layer(6)->position().z();
0072   if (std::abs(z1) > std::abs(z6)) {
0073     reverse(layerIndex.begin(), layerIndex.end());
0074     reverse(rechits.begin(), rechits.end());
0075   }
0076   if (rechits.size() < 2) {
0077     return std::vector<CSCSegment>();
0078   }
0079   // We have at least 2 hits. We intend to try all possible pairs of hits to start
0080   // segment building. 'All possible' means each hit lies on different layers in the chamber.
0081   // after all same size segs are build we get rid of the overcrossed segments using the chi2 criteria
0082   // the hits from the segs that are left are marked as used and are not added to segs in future iterations
0083   // the hits from 3p segs are marked as used separately in order to try to assamble them in longer segments
0084   // in case there is a second pass
0085   // Choose first hit (as close to IP as possible) h1 and second hit
0086   // (as far from IP as possible) h2 To do this we iterate over hits
0087   // in the chamber by layer - pick two layers. Then we
0088   // iterate over hits within each of these layers and pick h1 and h2
0089   // these. If they are 'close enough' we build an empty
0090   // segment. Then try adding hits to this segment.
0091   // Initialize flags that a given hit has been allocated to a segment
0092   BoolContainer used(rechits.size(), false);
0093   BoolContainer used3p(rechits.size(), false);
0094   // This is going to point to fits to hits, and its content will be used to create a CSCSegment
0095   AlgoState aState;
0096   aState.aChamber = aChamber;
0097   aState.doCollisions = doCollisions;
0098   aState.dRMax = dRMax;
0099   aState.dPhiMax = dPhiMax;
0100   aState.dRIntMax = dRIntMax;
0101   aState.dPhiIntMax = dPhiIntMax;
0102   aState.chi2Norm_2D_ = chi2Norm_2D_;
0103   aState.chi2_str_ = chi2_str_;
0104   aState.chi2Max = chi2Max;
0105 
0106   // Define buffer for segments we build
0107   std::vector<CSCSegment> segments;
0108   ChamberHitContainerCIt ib = rechits.begin();
0109   ChamberHitContainerCIt ie = rechits.end();
0110   // Possibly allow 3 passes, second widening scale factor for cuts, third for segments from displaced vertices
0111   aState.windowScale = 1.;  // scale factor for cuts
0112   bool search_disp = false;
0113   aState.strip_iadd = 1;
0114   aState.chi2D_iadd = 1;
0115   int npass = (wideSeg > 1.) ? 3 : 2;
0116   for (int ipass = 0; ipass < npass; ++ipass) {
0117     if (aState.windowScale > 1.) {
0118       iadd = 1;
0119       aState.strip_iadd = 4;
0120       aState.chi2D_iadd = 4;
0121       aState.chi2Max = 2 * chi2Max;
0122       if (rechits.size() <= 12)
0123         iadd = 0;  //allow 3 hit segments for low hit multiplicity chambers
0124     }
0125 
0126     int used_rh = 0;
0127     for (ChamberHitContainerCIt i1 = ib; i1 != ie; ++i1) {
0128       if (used[i1 - ib])
0129         used_rh++;
0130     }
0131 
0132     //change the tresholds if it's time to look for displaced mu segments
0133     if (aState.doCollisions && search_disp &&
0134         int(rechits.size() - used_rh) >
0135             2) {  //check if there are enough recHits left to build a segment from displaced vertices
0136       aState.doCollisions = false;
0137       aState.windowScale = 1.;  // scale factor for cuts
0138       aState.dRMax = 4.0;
0139       aState.dPhiMax = 4 * aState.dPhiMax;
0140       aState.dRIntMax = 4 * aState.dRIntMax;
0141       aState.dPhiIntMax = 4 * aState.dPhiIntMax;
0142       aState.chi2Norm_2D_ = 10 * aState.chi2Norm_2D_;
0143       aState.chi2_str_ = 200;
0144       aState.chi2Max = 4 * aState.chi2Max;
0145     } else {
0146       search_disp = false;  //make sure the flag is off
0147     }
0148 
0149     for (unsigned int n_seg_min = 6u; n_seg_min > 2u + iadd; --n_seg_min) {
0150       BoolContainer common_used(rechits.size(), false);
0151       std::array<BoolContainer, 120> common_used_it = {};
0152       for (unsigned int i = 0; i < common_used_it.size(); i++) {
0153         common_used_it[i] = common_used;
0154       }
0155       ChamberHitContainer best_proto_segment[120];
0156       float min_chi[120] = {9999};
0157       int common_it = 0;
0158       bool first_proto_segment = true;
0159       for (ChamberHitContainerCIt i1 = ib; i1 != ie; ++i1) {
0160         bool segok = false;
0161         //skip if rh is used and the layer tat has big rh multiplicity(>25RHs)
0162         if (used[i1 - ib] || recHits_per_layer[int(layerIndex[i1 - ib]) - 1] > 25 ||
0163             (n_seg_min == 3 && used3p[i1 - ib]))
0164           continue;
0165         int layer1 = layerIndex[i1 - ib];
0166         const CSCRecHit2D* h1 = *i1;
0167         for (ChamberHitContainerCIt i2 = ie - 1; i2 != i1; --i2) {
0168           if (used[i2 - ib] || recHits_per_layer[int(layerIndex[i2 - ib]) - 1] > 25 ||
0169               (n_seg_min == 3 && used3p[i2 - ib]))
0170             continue;
0171           int layer2 = layerIndex[i2 - ib];
0172           if ((abs(layer2 - layer1) + 1) < int(n_seg_min))
0173             break;  //decrease n_seg_min
0174           const CSCRecHit2D* h2 = *i2;
0175           if (areHitsCloseInR(aState, h1, h2) && areHitsCloseInGlobalPhi(aState, h1, h2)) {
0176             aState.proto_segment.clear();
0177             if (!addHit(aState, h1, layer1))
0178               continue;
0179             if (!addHit(aState, h2, layer2))
0180               continue;
0181             // Can only add hits if already have a segment
0182             if (aState.sfit)
0183               tryAddingHitsToSegment(aState, rechits, used, layerIndex, i1, i2);
0184             segok = isSegmentGood(aState, rechits);
0185             if (segok) {
0186               if (aState.proto_segment.size() > n_seg_min) {
0187                 baseline(aState, n_seg_min);
0188                 updateParameters(aState);
0189               }
0190               if (aState.sfit->chi2() > aState.chi2Norm_2D_ * aState.chi2D_iadd ||
0191                   aState.proto_segment.size() < n_seg_min)
0192                 aState.proto_segment.clear();
0193               if (!aState.proto_segment.empty()) {
0194                 updateParameters(aState);
0195                 //add same-size overcrossed protosegments to the collection
0196                 if (first_proto_segment) {
0197                   flagHitsAsUsed(aState, rechits, common_used_it[0]);
0198                   min_chi[0] = aState.sfit->chi2();
0199                   best_proto_segment[0] = aState.proto_segment;
0200                   first_proto_segment = false;
0201                 } else {  //for the rest of found proto_segments
0202                   common_it++;
0203                   flagHitsAsUsed(aState, rechits, common_used_it[common_it]);
0204                   min_chi[common_it] = aState.sfit->chi2();
0205                   best_proto_segment[common_it] = aState.proto_segment;
0206                   ChamberHitContainerCIt hi, iu, ik;
0207                   int iter = common_it;
0208                   for (iu = ib; iu != ie; ++iu) {
0209                     for (hi = aState.proto_segment.begin(); hi != aState.proto_segment.end(); ++hi) {
0210                       if (*hi == *iu) {
0211                         int merge_nr = -1;
0212                         for (int k = 0; k < iter + 1; k++) {
0213                           if (common_used_it[k][iu - ib] == true) {
0214                             if (merge_nr != -1) {
0215                               //merge the k and merge_nr collections of flaged hits into the merge_nr collection and unmark the k collection hits
0216                               for (ik = ib; ik != ie; ++ik) {
0217                                 if (common_used_it[k][ik - ib] == true) {
0218                                   common_used_it[merge_nr][ik - ib] = true;
0219                                   common_used_it[k][ik - ib] = false;
0220                                 }
0221                               }
0222                               //change best_protoseg if min_chi_2 is smaller
0223                               if (min_chi[k] < min_chi[merge_nr]) {
0224                                 min_chi[merge_nr] = min_chi[k];
0225                                 best_proto_segment[merge_nr] = best_proto_segment[k];
0226                                 best_proto_segment[k].clear();
0227                                 min_chi[k] = 9999;
0228                               }
0229                               common_it--;
0230                             } else {
0231                               merge_nr = k;
0232                             }
0233                           }  //if(common_used[k][iu-ib] == true)
0234                         }  //for k
0235                       }  //if
0236                     }  //for proto_seg
0237                   }  //for rec_hits
0238                 }  //else
0239               }  //proto seg not empty
0240             }
0241           }  // h1 & h2 close
0242           if (segok)
0243             break;
0244         }  // i2
0245       }  // i1
0246 
0247       //add the reconstructed segments
0248       for (int j = 0; j < common_it + 1; j++) {
0249         aState.proto_segment = best_proto_segment[j];
0250         best_proto_segment[j].clear();
0251         //SKIP empty proto-segments
0252         if (aState.proto_segment.empty())
0253           continue;
0254         updateParameters(aState);
0255         // Create an actual CSCSegment - retrieve all info from the fit
0256         CSCSegment temp(aState.sfit->hits(),
0257                         aState.sfit->intercept(),
0258                         aState.sfit->localdir(),
0259                         aState.sfit->covarianceMatrix(),
0260                         aState.sfit->chi2());
0261         aState.sfit = nullptr;
0262         segments.push_back(temp);
0263         //if the segment has 3 hits flag them as used in a particular way
0264         if (aState.proto_segment.size() == 3) {
0265           flagHitsAsUsed(aState, rechits, used3p);
0266         } else {
0267           flagHitsAsUsed(aState, rechits, used);
0268         }
0269         aState.proto_segment.clear();
0270       }
0271     }  //for n_seg_min
0272 
0273     if (search_disp) {
0274       //reset params and flags for the next chamber
0275       search_disp = false;
0276       aState.doCollisions = true;
0277       aState.dRMax = 2.0;
0278       aState.chi2_str_ = 100;
0279       aState.dPhiMax = 0.25 * aState.dPhiMax;
0280       aState.dRIntMax = 0.25 * aState.dRIntMax;
0281       aState.dPhiIntMax = 0.25 * aState.dPhiIntMax;
0282       aState.chi2Norm_2D_ = 0.1 * aState.chi2Norm_2D_;
0283       aState.chi2Max = 0.25 * aState.chi2Max;
0284     }
0285 
0286     std::vector<CSCSegment>::iterator it = segments.begin();
0287     bool good_segs = false;
0288     while (it != segments.end()) {
0289       if ((*it).nRecHits() > 3) {
0290         good_segs = true;
0291         break;
0292       }
0293       ++it;
0294     }
0295     if (good_segs &&
0296         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)
0297       search_disp = true;
0298       continue;  //proceed to search the segs from displaced vertices
0299     }
0300 
0301     // Increase cut windows by factor of wideSeg only for collisions
0302     if (!aState.doCollisions && !search_disp)
0303       break;
0304     aState.windowScale = wideSeg;
0305   }  // ipass
0306 
0307   //get rid of enchansed 3p segments
0308   std::vector<CSCSegment>::iterator it = segments.begin();
0309   while (it != segments.end()) {
0310     if ((*it).nRecHits() == 3) {
0311       bool found_common = false;
0312       const std::vector<CSCRecHit2D>& theseRH = (*it).specificRecHits();
0313       for (ChamberHitContainerCIt i1 = ib; i1 != ie; ++i1) {
0314         if (used[i1 - ib] && used3p[i1 - ib]) {
0315           const CSCRecHit2D* sh1 = *i1;
0316           CSCDetId id = sh1->cscDetId();
0317           int sh1layer = id.layer();
0318           int RH_centerid = sh1->nStrips() / 2;
0319           int RH_centerStrip = sh1->channels(RH_centerid);
0320           int RH_wg = sh1->hitWire();
0321           std::vector<CSCRecHit2D>::const_iterator sh;
0322           for (sh = theseRH.begin(); sh != theseRH.end(); ++sh) {
0323             CSCDetId idRH = sh->cscDetId();
0324             //find segment hit coord
0325             int shlayer = idRH.layer();
0326             int SegRH_centerid = sh->nStrips() / 2;
0327             int SegRH_centerStrip = sh->channels(SegRH_centerid);
0328             int SegRH_wg = sh->hitWire();
0329             if (sh1layer == shlayer && SegRH_centerStrip == RH_centerStrip && SegRH_wg == RH_wg) {
0330               //remove the enchansed 3p segment
0331               segments.erase(it, (it + 1));
0332               found_common = true;
0333               break;
0334             }
0335           }  //theserh
0336         }
0337         if (found_common)
0338           break;  //current seg has already been erased
0339       }  //camber hits
0340       if (!found_common)
0341         ++it;
0342     }  //its a 3p seg
0343     else {
0344       ++it;  //go to the next seg
0345     }
0346   }  //while
0347   // Give the segments to the CSCChamber
0348   return segments;
0349 }  //build segments
0350 
0351 void CSCSegAlgoRU::tryAddingHitsToSegment(AlgoState& aState,
0352                                           const ChamberHitContainer& rechits,
0353                                           const BoolContainer& used,
0354                                           const LayerIndex& layerIndex,
0355                                           const ChamberHitContainerCIt i1,
0356                                           const ChamberHitContainerCIt i2) const {
0357   // Iterate over the layers with hits in the chamber
0358   // Skip the layers containing the segment endpoints
0359   // Test each hit on the other layers to see if it is near the segment
0360   // If it is, see whether there is already a hit on the segment from the same layer
0361   // - if so, and there are more than 2 hits on the segment, copy the segment,
0362   // replace the old hit with the new hit. If the new segment chi2 is better
0363   // then replace the original segment with the new one (by swap)
0364   // - if not, copy the segment, add the hit. If the new chi2/dof is still satisfactory
0365   // then replace the original segment with the new one (by swap)
0366   ChamberHitContainerCIt ib = rechits.begin();
0367   ChamberHitContainerCIt ie = rechits.end();
0368   for (ChamberHitContainerCIt i = ib; i != ie; ++i) {
0369     int layer = layerIndex[i - ib];
0370     if (hasHitOnLayer(aState, layer) && aState.proto_segment.size() <= 2)
0371       continue;
0372     if (layerIndex[i - ib] == layerIndex[i1 - ib] || layerIndex[i - ib] == layerIndex[i2 - ib] || used[i - ib])
0373       continue;
0374 
0375     const CSCRecHit2D* h = *i;
0376     if (isHitNearSegment(aState, h)) {
0377       // Don't consider alternate hits on layers holding the two starting points
0378       if (hasHitOnLayer(aState, layer)) {
0379         if (aState.proto_segment.size() <= 2)
0380           continue;
0381         compareProtoSegment(aState, h, layer);
0382       } else {
0383         increaseProtoSegment(aState, h, layer, aState.chi2D_iadd);
0384       }
0385     }  // h & seg close
0386   }  // i
0387 }
0388 
0389 bool CSCSegAlgoRU::areHitsCloseInR(const AlgoState& aState, const CSCRecHit2D* h1, const CSCRecHit2D* h2) const {
0390   float maxWG_width[10] = {0, 0, 4.1, 5.69, 2.49, 5.06, 2.49, 5.06, 1.87, 5.06};
0391   CSCDetId id = h1->cscDetId();
0392   int iStn = id.iChamberType() - 1;
0393   //find maxWG_width for ME11 (tilt = 29deg)
0394   int wg_num = h2->hitWire();
0395   if (iStn == 0 || iStn == 1) {
0396     if (wg_num == 1) {
0397       maxWG_width[0] = 9.25;
0398       maxWG_width[1] = 9.25;
0399     }
0400     if (wg_num > 1 && wg_num < 48) {
0401       maxWG_width[0] = 3.14;
0402       maxWG_width[1] = 3.14;
0403     }
0404     if (wg_num == 48) {
0405       maxWG_width[0] = 10.75;
0406       maxWG_width[1] = 10.75;
0407     }
0408   }
0409   const CSCLayer* l1 = aState.aChamber->layer(h1->cscDetId().layer());
0410   GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
0411   const CSCLayer* l2 = aState.aChamber->layer(h2->cscDetId().layer());
0412   GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
0413   //find z to understand the direction
0414   float h1z = gp1.z();
0415   float h2z = gp2.z();
0416   //switch off the IP check for non collisions case
0417   if (!aState.doCollisions) {
0418     h1z = 1;
0419     h2z = 1;
0420   }
0421   return (gp2.perp() > ((gp1.perp() - aState.dRMax * aState.strip_iadd * maxWG_width[iStn]) * h2z) / h1z &&
0422           gp2.perp() < ((gp1.perp() + aState.dRMax * aState.strip_iadd * maxWG_width[iStn]) * h2z) / h1z)
0423              ? true
0424              : false;
0425 }
0426 
0427 bool CSCSegAlgoRU::areHitsCloseInGlobalPhi(const AlgoState& aState,
0428                                            const CSCRecHit2D* h1,
0429                                            const CSCRecHit2D* h2) const {
0430   float strip_width[10] = {0.003878509,
0431                            0.002958185,
0432                            0.002327105,
0433                            0.00152552,
0434                            0.00465421,
0435                            0.002327105,
0436                            0.00465421,
0437                            0.002327105,
0438                            0.00465421,
0439                            0.002327105};  //in rad
0440   const CSCLayer* l1 = aState.aChamber->layer(h1->cscDetId().layer());
0441   GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
0442   const CSCLayer* l2 = aState.aChamber->layer(h2->cscDetId().layer());
0443   GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
0444   float err_stpos_h1 = h1->errorWithinStrip();
0445   float err_stpos_h2 = h2->errorWithinStrip();
0446   CSCDetId id = h1->cscDetId();
0447   int iStn = id.iChamberType() - 1;
0448   float dphi_incr = 0;
0449   if (err_stpos_h1 > 0.25 * strip_width[iStn] || err_stpos_h2 > 0.25 * strip_width[iStn])
0450     dphi_incr = 0.5 * strip_width[iStn];
0451   float dphi12 = deltaPhi(gp1.barePhi(), gp2.barePhi());
0452   return (fabs(dphi12) < (aState.dPhiMax * aState.strip_iadd + dphi_incr)) ? true : false;  // +v
0453 }
0454 
0455 bool CSCSegAlgoRU::isHitNearSegment(const AlgoState& aState, const CSCRecHit2D* h) const {
0456   // Is hit near segment?
0457   // Requires deltaphi and deltaR within ranges specified in parameter set.
0458   // Note that to make intuitive cuts on delta(phi) one must work in
0459   // phi range (-pi, +pi] not [0, 2pi)
0460   float strip_width[10] = {0.003878509,
0461                            0.002958185,
0462                            0.002327105,
0463                            0.00152552,
0464                            0.00465421,
0465                            0.002327105,
0466                            0.00465421,
0467                            0.002327105,
0468                            0.00465421,
0469                            0.002327105};  //in rad
0470   const CSCLayer* l1 = aState.aChamber->layer((*(aState.proto_segment.begin()))->cscDetId().layer());
0471   GlobalPoint gp1 = l1->toGlobal((*(aState.proto_segment.begin()))->localPosition());
0472   const CSCLayer* l2 = aState.aChamber->layer((*(aState.proto_segment.begin() + 1))->cscDetId().layer());
0473   GlobalPoint gp2 = l2->toGlobal((*(aState.proto_segment.begin() + 1))->localPosition());
0474   float err_stpos_h1 = (*(aState.proto_segment.begin()))->errorWithinStrip();
0475   float err_stpos_h2 = (*(aState.proto_segment.begin() + 1))->errorWithinStrip();
0476   const CSCLayer* l = aState.aChamber->layer(h->cscDetId().layer());
0477   GlobalPoint hp = l->toGlobal(h->localPosition());
0478   float err_stpos_h = h->errorWithinStrip();
0479   float hphi = hp.phi();  // in (-pi, +pi]
0480   if (hphi < 0.)
0481     hphi += 2. * M_PI;                  // into range [0, 2pi)
0482   float sphi = phiAtZ(aState, hp.z());  // in [0, 2*pi)
0483   float phidif = sphi - hphi;
0484   if (phidif < 0.)
0485     phidif += 2. * M_PI;  // into range [0, 2pi)
0486   if (phidif > M_PI)
0487     phidif -= 2. * M_PI;  // into range (-pi, pi]
0488   SVector6 r_glob;
0489   CSCDetId id = h->cscDetId();
0490   int iStn = id.iChamberType() - 1;
0491   float dphi_incr = 0;
0492   float pos_str = 1;
0493   //increase dPhi cut if the hit is on the edge of the strip
0494   float stpos = (*h).positionWithinStrip();
0495   bool centr_str = false;
0496   if (iStn != 0 && iStn != 1) {
0497     if (stpos > -0.25 && stpos < 0.25)
0498       centr_str = true;
0499   }
0500   if (err_stpos_h1 < 0.25 * strip_width[iStn] || err_stpos_h2 < 0.25 * strip_width[iStn] ||
0501       err_stpos_h < 0.25 * strip_width[iStn]) {
0502     dphi_incr = 0.5 * strip_width[iStn];
0503   } else {
0504     if (centr_str)
0505       pos_str = 1.3;
0506   }
0507   r_glob((*(aState.proto_segment.begin()))->cscDetId().layer() - 1) = gp1.perp();
0508   r_glob((*(aState.proto_segment.begin() + 1))->cscDetId().layer() - 1) = gp2.perp();
0509   float R = hp.perp();
0510   int layer = h->cscDetId().layer();
0511   float r_interpolated = fit_r_phi(aState, r_glob, layer);
0512   float dr = fabs(r_interpolated - R);
0513   float maxWG_width[10] = {0, 0, 4.1, 5.69, 2.49, 5.06, 2.49, 5.06, 1.87, 5.06};
0514   //find maxWG_width for ME11 (tilt = 29deg)
0515   int wg_num = h->hitWire();
0516   if (iStn == 0 || iStn == 1) {
0517     if (wg_num == 1) {
0518       maxWG_width[0] = 9.25;
0519       maxWG_width[1] = 9.25;
0520     }
0521     if (wg_num > 1 && wg_num < 48) {
0522       maxWG_width[0] = 3.14;
0523       maxWG_width[1] = 3.14;
0524     }
0525     if (wg_num == 48) {
0526       maxWG_width[0] = 10.75;
0527       maxWG_width[1] = 10.75;
0528     }
0529   }
0530   return (fabs(phidif) < aState.dPhiIntMax * aState.strip_iadd * pos_str + dphi_incr &&
0531           fabs(dr) < aState.dRIntMax * aState.strip_iadd * maxWG_width[iStn])
0532              ? true
0533              : false;
0534 }
0535 
0536 float CSCSegAlgoRU::phiAtZ(const AlgoState& aState, float z) const {
0537   if (!aState.sfit)
0538     return 0.;
0539   // Returns a phi in [ 0, 2*pi )
0540   const CSCLayer* l1 = aState.aChamber->layer((*(aState.proto_segment.begin()))->cscDetId().layer());
0541   GlobalPoint gp = l1->toGlobal(aState.sfit->intercept());
0542   GlobalVector gv = l1->toGlobal(aState.sfit->localdir());
0543   float x = gp.x() + (gv.x() / gv.z()) * (z - gp.z());
0544   float y = gp.y() + (gv.y() / gv.z()) * (z - gp.z());
0545   float phi = atan2(y, x);
0546   if (phi < 0.f)
0547     phi += 2. * M_PI;
0548   return phi;
0549 }
0550 
0551 bool CSCSegAlgoRU::isSegmentGood(const AlgoState& aState, const ChamberHitContainer& rechitsInChamber) const {
0552   // If the chamber has 20 hits or fewer, require at least 3 hits on segment
0553   // If the chamber has >20 hits require at least 4 hits
0554   // If it's the second cycle of the builder and there are <= 12 hits in chamber, require at least 3 hits on segment
0555   //@@ THESE VALUES SHOULD BECOME PARAMETERS?
0556   bool ok = false;
0557   unsigned int iadd = (rechitsInChamber.size() > 20) ? 1 : 0;
0558   if (aState.windowScale > 1.) {
0559     iadd = 1;
0560     if (rechitsInChamber.size() <= 12)
0561       iadd = 0;
0562   }
0563   if (aState.proto_segment.size() >= 3 + iadd)
0564     ok = true;
0565   return ok;
0566 }
0567 
0568 void CSCSegAlgoRU::flagHitsAsUsed(const AlgoState& aState,
0569                                   const ChamberHitContainer& rechitsInChamber,
0570                                   BoolContainer& used) const {
0571   // Flag hits on segment as used
0572   ChamberHitContainerCIt ib = rechitsInChamber.begin();
0573   ChamberHitContainerCIt hi, iu;
0574   for (hi = aState.proto_segment.begin(); hi != aState.proto_segment.end(); ++hi) {
0575     for (iu = ib; iu != rechitsInChamber.end(); ++iu) {
0576       if (*hi == *iu)
0577         used[iu - ib] = true;
0578     }
0579   }
0580 }
0581 
0582 bool CSCSegAlgoRU::addHit(AlgoState& aState, const CSCRecHit2D* aHit, int layer) const {
0583   // Return true if hit was added successfully
0584   // (and then parameters are updated).
0585   // Return false if there is already a hit on the same layer, or insert failed.
0586   ChamberHitContainer::const_iterator it;
0587   for (it = aState.proto_segment.begin(); it != aState.proto_segment.end(); it++)
0588     if (((*it)->cscDetId().layer() == layer) && (aHit != (*it)))
0589       return false;
0590   aState.proto_segment.push_back(aHit);
0591   // make a fit
0592   updateParameters(aState);
0593   return true;
0594 }
0595 
0596 void CSCSegAlgoRU::updateParameters(AlgoState& aState) const {
0597   // Delete input CSCSegFit, create a new one and make the fit
0598   // delete sfit;
0599   aState.sfit = std::make_unique<CSCSegFit>(aState.aChamber, aState.proto_segment);
0600   aState.sfit->fit();
0601 }
0602 
0603 float CSCSegAlgoRU::fit_r_phi(const AlgoState& aState, const SVector6& points, int layer) const {
0604   //find R or Phi on the given layer using the given points for the interpolation
0605   float Sx = 0;
0606   float Sy = 0;
0607   float Sxx = 0;
0608   float Sxy = 0;
0609   for (int i = 1; i < 7; i++) {
0610     if (points(i - 1) == 0.)
0611       continue;
0612     Sy = Sy + (points(i - 1));
0613     Sx = Sx + i;
0614     Sxx = Sxx + (i * i);
0615     Sxy = Sxy + ((i)*points(i - 1));
0616   }
0617   float delta = 2 * Sxx - Sx * Sx;
0618   float intercept = (Sxx * Sy - Sx * Sxy) / delta;
0619   float slope = (2 * Sxy - Sx * Sy) / delta;
0620   return (intercept + slope * layer);
0621 }
0622 
0623 void CSCSegAlgoRU::baseline(AlgoState& aState, int n_seg_min) const {
0624   int nhits = aState.proto_segment.size();
0625   //initialise vectors for strip position and error within strip
0626   SVector6 sp;
0627   SVector6 se;
0628   unsigned int init_size = aState.proto_segment.size();
0629   ChamberHitContainer buffer;
0630   buffer.clear();
0631   buffer.reserve(init_size);
0632   while (buffer.size() < init_size) {
0633     ChamberHitContainer::iterator min;
0634     int min_layer = 10;
0635     for (ChamberHitContainer::iterator k = aState.proto_segment.begin(); k != aState.proto_segment.end(); k++) {
0636       const CSCRecHit2D* iRHk = *k;
0637       CSCDetId idRHk = iRHk->cscDetId();
0638       int kLayer = idRHk.layer();
0639       if (kLayer < min_layer) {
0640         min_layer = kLayer;
0641         min = k;
0642       }
0643     }
0644     buffer.push_back(*min);
0645     aState.proto_segment.erase(min);
0646   }  //while
0647 
0648   aState.proto_segment.clear();
0649   for (ChamberHitContainer::const_iterator cand = buffer.begin(); cand != buffer.end(); cand++) {
0650     aState.proto_segment.push_back(*cand);
0651   }
0652 
0653   for (ChamberHitContainer::const_iterator iRH = aState.proto_segment.begin(); iRH != aState.proto_segment.end();
0654        iRH++) {
0655     const CSCRecHit2D* iRHp = *iRH;
0656     CSCDetId idRH = iRHp->cscDetId();
0657     int kRing = idRH.ring();
0658     int kStation = idRH.station();
0659     int kLayer = idRH.layer();
0660     // Find the strip containing this hit
0661     int centerid = iRHp->nStrips() / 2;
0662     int centerStrip = iRHp->channels(centerid);
0663     float stpos = (*iRHp).positionWithinStrip();
0664     se(kLayer - 1) = (*iRHp).errorWithinStrip();
0665     // Take into account half-strip staggering of layers (ME1/1 has no staggering)
0666     if (kStation == 1 && (kRing == 1 || kRing == 4))
0667       sp(kLayer - 1) = stpos + centerStrip;
0668     else {
0669       if (kLayer == 1 || kLayer == 3 || kLayer == 5)
0670         sp(kLayer - 1) = stpos + centerStrip;
0671       if (kLayer == 2 || kLayer == 4 || kLayer == 6)
0672         sp(kLayer - 1) = stpos - 0.5 + centerStrip;
0673     }
0674   }
0675   float chi2_str;
0676   fitX(aState, sp, se, -1, -1, chi2_str);
0677 
0678   //-----------------------------------------------------
0679   // Optimal point rejection method
0680   //-----------------------------------------------------
0681   float minSum = 1000;
0682   int iworst = -1;
0683   int bad_layer = -1;
0684   ChamberHitContainer::const_iterator rh_to_be_deleted_1;
0685   ChamberHitContainer::const_iterator rh_to_be_deleted_2;
0686   if ((chi2_str) > aState.chi2_str_ * aState.chi2D_iadd) {  ///(nhits-2)
0687     for (ChamberHitContainer::const_iterator i1 = aState.proto_segment.begin(); i1 != aState.proto_segment.end();
0688          ++i1) {
0689       const CSCRecHit2D* i1_1 = *i1;
0690       CSCDetId idRH1 = i1_1->cscDetId();
0691       int z1 = idRH1.layer();
0692       for (ChamberHitContainer::const_iterator i2 = i1 + 1; i2 != aState.proto_segment.end(); ++i2) {
0693         const CSCRecHit2D* i2_1 = *i2;
0694         CSCDetId idRH2 = i2_1->cscDetId();
0695         int z2 = idRH2.layer();
0696         int irej = 0;
0697         for (ChamberHitContainer::const_iterator ir = aState.proto_segment.begin(); ir != aState.proto_segment.end();
0698              ++ir) {
0699           ++irej;
0700           if (ir == i1 || ir == i2)
0701             continue;
0702           float dsum = 0;
0703           const CSCRecHit2D* ir_1 = *ir;
0704           CSCDetId idRH = ir_1->cscDetId();
0705           int worst_layer = idRH.layer();
0706           for (ChamberHitContainer::const_iterator i = aState.proto_segment.begin(); i != aState.proto_segment.end();
0707                ++i) {
0708             const CSCRecHit2D* i_1 = *i;
0709             if (i == i1 || i == i2 || i == ir)
0710               continue;
0711             float slope = (sp(z2 - 1) - sp(z1 - 1)) / (z2 - z1);
0712             float intersept = sp(z1 - 1) - slope * z1;
0713             CSCDetId idRH = i_1->cscDetId();
0714             int z = idRH.layer();
0715             float di = fabs(sp(z - 1) - intersept - slope * z);
0716             dsum = dsum + di;
0717           }  //i
0718           if (dsum < minSum) {
0719             minSum = dsum;
0720             bad_layer = worst_layer;
0721             iworst = irej;
0722             rh_to_be_deleted_1 = ir;
0723           }
0724         }  //ir
0725       }  //i2
0726     }  //i1
0727     fitX(aState, sp, se, bad_layer, -1, chi2_str);
0728   }  //if chi2prob<1.0e-4
0729 
0730   //find worst from n-1 hits
0731   int iworst2 = -1;
0732   int bad_layer2 = -1;
0733   if (iworst > -1 && (nhits - 1) > n_seg_min && (chi2_str) > aState.chi2_str_ * aState.chi2D_iadd) {  ///(nhits-3)
0734     iworst = -1;
0735     float minSum = 1000;
0736     for (ChamberHitContainer::const_iterator i1 = aState.proto_segment.begin(); i1 != aState.proto_segment.end();
0737          ++i1) {
0738       const CSCRecHit2D* i1_1 = *i1;
0739       CSCDetId idRH1 = i1_1->cscDetId();
0740       int z1 = idRH1.layer();
0741       for (ChamberHitContainer::const_iterator i2 = i1 + 1; i2 != aState.proto_segment.end(); ++i2) {
0742         const CSCRecHit2D* i2_1 = *i2;
0743         CSCDetId idRH2 = i2_1->cscDetId();
0744         int z2 = idRH2.layer();
0745         int irej = 0;
0746         for (ChamberHitContainer::const_iterator ir = aState.proto_segment.begin(); ir != aState.proto_segment.end();
0747              ++ir) {
0748           ++irej;
0749           int irej2 = 0;
0750           if (ir == i1 || ir == i2)
0751             continue;
0752           const CSCRecHit2D* ir_1 = *ir;
0753           CSCDetId idRH = ir_1->cscDetId();
0754           int worst_layer = idRH.layer();
0755           for (ChamberHitContainer::const_iterator ir2 = aState.proto_segment.begin();
0756                ir2 != aState.proto_segment.end();
0757                ++ir2) {
0758             ++irej2;
0759             if (ir2 == i1 || ir2 == i2 || ir2 == ir)
0760               continue;
0761             float dsum = 0;
0762             const CSCRecHit2D* ir2_1 = *ir2;
0763             CSCDetId idRH = ir2_1->cscDetId();
0764             int worst_layer2 = idRH.layer();
0765             for (ChamberHitContainer::const_iterator i = aState.proto_segment.begin(); i != aState.proto_segment.end();
0766                  ++i) {
0767               const CSCRecHit2D* i_1 = *i;
0768               if (i == i1 || i == i2 || i == ir || i == ir2)
0769                 continue;
0770               float slope = (sp(z2 - 1) - sp(z1 - 1)) / (z2 - z1);
0771               float intersept = sp(z1 - 1) - slope * z1;
0772               CSCDetId idRH = i_1->cscDetId();
0773               int z = idRH.layer();
0774               float di = fabs(sp(z - 1) - intersept - slope * z);
0775               dsum = dsum + di;
0776             }  //i
0777             if (dsum < minSum) {
0778               minSum = dsum;
0779               iworst2 = irej2;
0780               iworst = irej;
0781               bad_layer = worst_layer;
0782               bad_layer2 = worst_layer2;
0783               rh_to_be_deleted_1 = ir;
0784               rh_to_be_deleted_2 = ir2;
0785             }
0786           }  //ir2
0787         }  //ir
0788       }  //i2
0789     }  //i1
0790     fitX(aState, sp, se, bad_layer, bad_layer2, chi2_str);
0791   }  //if prob(n-1)<e-4
0792 
0793   //----------------------------------
0794   //erase bad_hits
0795   //----------------------------------
0796   if (iworst2 - 1 >= 0 && iworst2 <= int(aState.proto_segment.size())) {
0797     aState.proto_segment.erase(rh_to_be_deleted_2);
0798   }
0799   if (iworst - 1 >= 0 && iworst <= int(aState.proto_segment.size())) {
0800     aState.proto_segment.erase(rh_to_be_deleted_1);
0801   }
0802 }
0803 
0804 float CSCSegAlgoRU::fitX(
0805     const AlgoState& aState, SVector6 points, SVector6 errors, int ir, int ir2, float& chi2_str) const {
0806   float S = 0;
0807   float Sx = 0;
0808   float Sy = 0;
0809   float Sxx = 0;
0810   float Sxy = 0;
0811   float sigma2 = 0;
0812   for (int i = 1; i < 7; i++) {
0813     if (i == ir || i == ir2 || points(i - 1) == 0.)
0814       continue;
0815     sigma2 = errors(i - 1) * errors(i - 1);
0816     float i1 = i - 3.5;
0817     S = S + (1 / sigma2);
0818     Sy = Sy + (points(i - 1) / sigma2);
0819     Sx = Sx + ((i1) / sigma2);
0820     Sxx = Sxx + (i1 * i1) / sigma2;
0821     Sxy = Sxy + (((i1)*points(i - 1)) / sigma2);
0822   }
0823   float delta = S * Sxx - Sx * Sx;
0824   float intercept = (Sxx * Sy - Sx * Sxy) / delta;
0825   float slope = (S * Sxy - Sx * Sy) / delta;
0826   float chi_str = 0;
0827   chi2_str = 0;
0828   // calculate chi2_str
0829   for (int i = 1; i < 7; i++) {
0830     if (i == ir || i == ir2 || points(i - 1) == 0.)
0831       continue;
0832     chi_str = (points(i - 1) - intercept - slope * (i - 3.5)) / (errors(i - 1));
0833     chi2_str = chi2_str + chi_str * chi_str;
0834   }
0835   return (intercept + slope * 0);
0836 }
0837 
0838 bool CSCSegAlgoRU::hasHitOnLayer(const AlgoState& aState, int layer) const {
0839   // Is there is already a hit on this layer?
0840   ChamberHitContainerCIt it;
0841   for (it = aState.proto_segment.begin(); it != aState.proto_segment.end(); it++)
0842     if ((*it)->cscDetId().layer() == layer)
0843       return true;
0844   return false;
0845 }
0846 
0847 bool CSCSegAlgoRU::replaceHit(AlgoState& aState, const CSCRecHit2D* h, int layer) const {
0848   // replace a hit from a layer
0849   ChamberHitContainer::const_iterator it;
0850   for (it = aState.proto_segment.begin(); it != aState.proto_segment.end();) {
0851     if ((*it)->cscDetId().layer() == layer)
0852       it = aState.proto_segment.erase(it);
0853     else
0854       ++it;
0855   }
0856   return addHit(aState, h, layer);
0857 }
0858 
0859 void CSCSegAlgoRU::compareProtoSegment(AlgoState& aState, const CSCRecHit2D* h, int layer) const {
0860   // Copy the input CSCSegFit
0861   std::unique_ptr<CSCSegFit> oldfit;
0862   oldfit = std::make_unique<CSCSegFit>(aState.aChamber, aState.proto_segment);
0863   oldfit->fit();
0864   ChamberHitContainer oldproto = aState.proto_segment;
0865 
0866   // May create a new fit
0867   bool ok = replaceHit(aState, h, layer);
0868   if ((aState.sfit->chi2() >= oldfit->chi2()) || !ok) {
0869     // keep original fit
0870     aState.proto_segment = oldproto;
0871     aState.sfit = std::move(oldfit);  // reset to the original input fit
0872   }
0873 }
0874 
0875 void CSCSegAlgoRU::increaseProtoSegment(AlgoState& aState, const CSCRecHit2D* h, int layer, int chi2_factor) const {
0876   // Creates a new fit
0877   std::unique_ptr<CSCSegFit> oldfit;
0878   ChamberHitContainer oldproto = aState.proto_segment;
0879   oldfit = std::make_unique<CSCSegFit>(aState.aChamber, aState.proto_segment);
0880   oldfit->fit();
0881 
0882   bool ok = addHit(aState, h, layer);
0883   //@@ TEST ON ndof<=0 IS JUST TO ACCEPT nhits=2 CASE??
0884   if (!ok || ((aState.sfit->ndof() > 0) && (aState.sfit->chi2() / aState.sfit->ndof() >= aState.chi2Max))) {
0885     // reset to original fit
0886     aState.proto_segment = oldproto;
0887     aState.sfit = std::move(oldfit);
0888   }
0889 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!