Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:24:41

0001 
0002 /**
0003  * \File CSCSegAlgoST.cc
0004  *
0005  *  \authors: S. Stoynev - NU
0006  *            I. Bloch    - FNAL
0007  *            E. James    - FNAL
0008  *            A. Sakharov - WSU
0009  *            T. Cox - UC Davis - segment fit factored out of entangled code - Jan 2015
0010  */
0011 
0012 #include "CSCSegAlgoST.h"
0013 #include "CSCCondSegFit.h"
0014 #include "CSCSegAlgoShowering.h"
0015 
0016 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0017 
0018 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0019 
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 
0022 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0023 
0024 #include <algorithm>
0025 #include <cmath>
0026 #include <iostream>
0027 #include <string>
0028 
0029 /* constructor
0030  *
0031  */
0032 CSCSegAlgoST::CSCSegAlgoST(const edm::ParameterSet& ps)
0033     : CSCSegmentAlgorithm(ps), myName_("CSCSegAlgoST"), ps_(ps), showering_(nullptr) {
0034   debug = ps.getUntrackedParameter<bool>("CSCDebug");
0035   //  minLayersApart         = ps.getParameter<int>("minLayersApart");
0036   //  nSigmaFromSegment      = ps.getParameter<double>("nSigmaFromSegment");
0037   minHitsPerSegment = ps.getParameter<int>("minHitsPerSegment");
0038   //  muonsPerChamberMax     = ps.getParameter<int>("CSCSegmentPerChamberMax");
0039   //  chi2Max                = ps.getParameter<double>("chi2Max");
0040   dXclusBoxMax = ps.getParameter<double>("dXclusBoxMax");
0041   dYclusBoxMax = ps.getParameter<double>("dYclusBoxMax");
0042   preClustering = ps.getParameter<bool>("preClustering");
0043   preClustering_useChaining = ps.getParameter<bool>("preClusteringUseChaining");
0044   Pruning = ps.getParameter<bool>("Pruning");
0045   BrutePruning = ps.getParameter<bool>("BrutePruning");
0046   BPMinImprovement = ps.getParameter<double>("BPMinImprovement");
0047   // maxRecHitsInCluster is the maximal number of hits in a precluster that is being processed
0048   // This cut is intended to remove messy events. Currently nothing is returned if there are
0049   // more that maxRecHitsInCluster hits. It could be useful to return an estimate of the
0050   // cluster position, which is available.
0051   maxRecHitsInCluster = ps.getParameter<int>("maxRecHitsInCluster");
0052   onlyBestSegment = ps.getParameter<bool>("onlyBestSegment");
0053 
0054   hitDropLimit4Hits = ps.getParameter<double>("hitDropLimit4Hits");
0055   hitDropLimit5Hits = ps.getParameter<double>("hitDropLimit5Hits");
0056   hitDropLimit6Hits = ps.getParameter<double>("hitDropLimit6Hits");
0057 
0058   yweightPenaltyThreshold = ps.getParameter<double>("yweightPenaltyThreshold");
0059   yweightPenalty = ps.getParameter<double>("yweightPenalty");
0060 
0061   curvePenaltyThreshold = ps.getParameter<double>("curvePenaltyThreshold");
0062   curvePenalty = ps.getParameter<double>("curvePenalty");
0063 
0064   useShowering = ps.getParameter<bool>("useShowering");
0065   if (useShowering)
0066     showering_ = new CSCSegAlgoShowering(ps);
0067 
0068   /// Improve the covariance matrix?
0069   adjustCovariance_ = ps.getParameter<bool>("CorrectTheErrors");
0070 
0071   chi2Norm_3D_ = ps.getParameter<double>("NormChi2Cut3D");
0072   prePrun_ = ps.getParameter<bool>("prePrun");
0073   prePrunLimit_ = ps.getParameter<double>("prePrunLimit");
0074 
0075   if (debug)
0076     edm::LogVerbatim("CSCSegment") << "CSCSegAlgoST: with factored conditioned segment fit";
0077 }
0078 
0079 /* Destructor
0080  *
0081  */
0082 CSCSegAlgoST::~CSCSegAlgoST() { delete showering_; }
0083 
0084 std::vector<CSCSegment> CSCSegAlgoST::run(const CSCChamber* aChamber, const ChamberHitContainer& rechits) {
0085   // Set member variable
0086   theChamber = aChamber;
0087 
0088   LogTrace("CSCSegAlgoST") << "[CSCSegAlgoST::run] Start building segments in chamber " << theChamber->id();
0089 
0090   // pre-cluster rechits and loop over all sub clusters seperately
0091   std::vector<CSCSegment> segments_temp;
0092   std::vector<CSCSegment> segments;
0093   std::vector<ChamberHitContainer> rechits_clusters;  // a collection of clusters of rechits
0094 
0095   // Define yweight penalty depending on chamber.
0096   // We fixed the relative ratios, but they can be scaled by parameters:
0097 
0098   for (int a = 0; a < 5; ++a) {
0099     for (int b = 0; b < 5; ++b) {
0100       a_yweightPenaltyThreshold[a][b] = 0.0;
0101     }
0102   }
0103 
0104   a_yweightPenaltyThreshold[1][1] = yweightPenaltyThreshold * 10.20;
0105   a_yweightPenaltyThreshold[1][2] = yweightPenaltyThreshold * 14.00;
0106   a_yweightPenaltyThreshold[1][3] = yweightPenaltyThreshold * 20.40;
0107   a_yweightPenaltyThreshold[1][4] = yweightPenaltyThreshold * 10.20;
0108   a_yweightPenaltyThreshold[2][1] = yweightPenaltyThreshold * 7.60;
0109   a_yweightPenaltyThreshold[2][2] = yweightPenaltyThreshold * 20.40;
0110   a_yweightPenaltyThreshold[3][1] = yweightPenaltyThreshold * 7.60;
0111   a_yweightPenaltyThreshold[3][2] = yweightPenaltyThreshold * 20.40;
0112   a_yweightPenaltyThreshold[4][1] = yweightPenaltyThreshold * 6.75;
0113   a_yweightPenaltyThreshold[4][2] = yweightPenaltyThreshold * 20.40;
0114 
0115   if (preClustering) {
0116     // run a pre-clusterer on the given rechits to split clearly-separated segment seeds:
0117     if (preClustering_useChaining) {
0118       // it uses X,Y,Z information; there are no configurable parameters used;
0119       // the X, Y, Z "cuts" are just (much) wider than the LCT readout ones
0120       // (which are actually not step functions); this new code could accomodate
0121       // the clusterHits one below but we leave it for security and backward
0122       // comparison reasons
0123       rechits_clusters = chainHits(theChamber, rechits);
0124     } else {
0125       // it uses X,Y information + configurable parameters
0126       rechits_clusters = clusterHits(theChamber, rechits);
0127     }
0128     // loop over the found clusters:
0129     for (std::vector<ChamberHitContainer>::iterator sub_rechits = rechits_clusters.begin();
0130          sub_rechits != rechits_clusters.end();
0131          ++sub_rechits) {
0132       // clear the buffer for the subset of segments:
0133       segments_temp.clear();
0134       // build the subset of segments:
0135       segments_temp = buildSegments((*sub_rechits));
0136       // add the found subset of segments to the collection of all segments in this chamber:
0137       segments.insert(segments.end(), segments_temp.begin(), segments_temp.end());
0138     }
0139     // Any pruning of hits?
0140     if (Pruning) {
0141       segments_temp.clear();  // segments_temp needed?!?!
0142       segments_temp = prune_bad_hits(theChamber, segments);
0143       segments.clear();              // segments_temp needed?!?!
0144       segments.swap(segments_temp);  // segments_temp needed?!?!
0145     }
0146 
0147     // Ganged strips in ME1/1A?
0148     if (("ME1/a" == aChamber->specs()->chamberTypeName()) && aChamber->specs()->gangedStrips()) {
0149       //  if ( aChamber->specs()->gangedStrips() ){
0150       findDuplicates(segments);
0151     }
0152     return segments;
0153   } else {
0154     segments = buildSegments(rechits);
0155     if (Pruning) {
0156       segments_temp.clear();  // segments_temp needed?!?!
0157       segments_temp = prune_bad_hits(theChamber, segments);
0158       segments.clear();              // segments_temp needed?!?!
0159       segments.swap(segments_temp);  // segments_temp needed?!?!
0160     }
0161 
0162     // Ganged strips in ME1/1A?
0163     if (("ME1/a" == aChamber->specs()->chamberTypeName()) && aChamber->specs()->gangedStrips()) {
0164       //  if ( aChamber->specs()->gangedStrips() ){
0165       findDuplicates(segments);
0166     }
0167     return segments;
0168     //return buildSegments(rechits);
0169   }
0170 }
0171 
0172 // ********************************************************************;
0173 // *** This method is meant to remove clearly bad hits, using as    ***;
0174 // *** much information from the chamber as possible (e.g. charge,  ***;
0175 // *** hit position, timing, etc.)                                  ***;
0176 // ********************************************************************;
0177 std::vector<CSCSegment> CSCSegAlgoST::prune_bad_hits(const CSCChamber* aChamber, std::vector<CSCSegment>& segments) {
0178   //   std::cout<<"*************************************************************"<<std::endl;
0179   //   std::cout<<"Called prune_bad_hits in Chamber "<< theChamber->specs()->chamberTypeName()<<std::endl;
0180   //   std::cout<<"*************************************************************"<<std::endl;
0181 
0182   std::vector<CSCSegment> segments_temp;
0183   std::vector<ChamberHitContainer> rechits_clusters;  // this is a collection of groups of rechits
0184 
0185   const float chi2ndfProbMin = 1.0e-4;
0186   bool use_brute_force = BrutePruning;
0187 
0188   int hit_nr = 0;
0189   int hit_nr_worst = -1;
0190   //int hit_nr_2ndworst = -1;
0191 
0192   for (std::vector<CSCSegment>::iterator it = segments.begin(); it != segments.end(); ++it) {
0193     // do nothing for nhit <= minHitPerSegment
0194     if ((*it).nRecHits() <= minHitsPerSegment)
0195       continue;
0196 
0197     if (!use_brute_force) {  // find worst hit
0198 
0199       float chisq = (*it).chi2();
0200       int nhits = (*it).nRecHits();
0201       LocalPoint localPos = (*it).localPosition();
0202       LocalVector segDir = (*it).localDirection();
0203       const CSCChamber* cscchamber = theChamber;
0204       float globZ;
0205 
0206       GlobalPoint globalPosition = cscchamber->toGlobal(localPos);
0207       globZ = globalPosition.z();
0208 
0209       if (ChiSquaredProbability((double)chisq, (double)(2 * nhits - 4)) < chi2ndfProbMin) {
0210         // find (rough) "residuals" (NOT excluding the hit from the fit - speed!) of hits on segment
0211         std::vector<CSCRecHit2D> theseRecHits = (*it).specificRecHits();
0212         std::vector<CSCRecHit2D>::const_iterator iRH_worst;
0213         //float xdist_local       = -99999.;
0214 
0215         float xdist_local_worst_sig = -99999.;
0216         float xdist_local_2ndworst_sig = -99999.;
0217         float xdist_local_sig = -99999.;
0218 
0219         hit_nr = 0;
0220         hit_nr_worst = -1;
0221         //hit_nr_2ndworst = -1;
0222 
0223         for (std::vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); ++iRH) {
0224           //mark "worst" hit:
0225 
0226           //float z_at_target ;
0227           //float radius      ;
0228           float loc_x_at_target;
0229           //float loc_y_at_target;
0230           //float loc_z_at_target;
0231 
0232           //z_at_target  = 0.;
0233           //radius       = 0.;
0234 
0235           // set the z target in CMS global coordinates:
0236           const CSCLayer* csclayerRH = theChamber->layer((*iRH).cscDetId().layer());
0237           LocalPoint localPositionRH = (*iRH).localPosition();
0238           GlobalPoint globalPositionRH = csclayerRH->toGlobal(localPositionRH);
0239 
0240           LocalError rerrlocal = (*iRH).localPositionError();
0241           float xxerr = rerrlocal.xx();
0242 
0243           float target_z = globalPositionRH.z();  // target z position in cm (z pos of the hit)
0244 
0245           if (target_z > 0.) {
0246             loc_x_at_target = localPos.x() + (segDir.x() / fabs(segDir.z()) * (target_z - globZ));
0247             //loc_y_at_target = localPos.y() + (segDir.y()/fabs(segDir.z())*( target_z - globZ ));
0248             //loc_z_at_target = target_z;
0249           } else {
0250             loc_x_at_target = localPos.x() + ((-1) * segDir.x() / fabs(segDir.z()) * (target_z - globZ));
0251             //loc_y_at_target = localPos.y() + ((-1)*segDir.y()/fabs(segDir.z())*( target_z - globZ ));
0252             //loc_z_at_target = target_z;
0253           }
0254           // have to transform the segments coordinates back to the local frame... how?!!!!!!!!!!!!
0255 
0256           //xdist_local  = fabs(localPositionRH.x() - loc_x_at_target);
0257           xdist_local_sig = fabs((localPositionRH.x() - loc_x_at_target) / (xxerr));
0258 
0259           if (xdist_local_sig > xdist_local_worst_sig) {
0260             xdist_local_2ndworst_sig = xdist_local_worst_sig;
0261             xdist_local_worst_sig = xdist_local_sig;
0262             iRH_worst = iRH;
0263             //hit_nr_2ndworst = hit_nr_worst;
0264             hit_nr_worst = hit_nr;
0265           } else if (xdist_local_sig > xdist_local_2ndworst_sig) {
0266             xdist_local_2ndworst_sig = xdist_local_sig;
0267             //hit_nr_2ndworst = hit_nr;
0268           }
0269           ++hit_nr;
0270         }
0271 
0272         // reset worst hit number if certain criteria apply.
0273         // Criteria: 2nd worst hit must be at least a factor of
0274         // 1.5 better than the worst in terms of sigma:
0275         if (xdist_local_worst_sig / xdist_local_2ndworst_sig < 1.5) {
0276           hit_nr_worst = -1;
0277           //hit_nr_2ndworst = -1;
0278         }
0279       }
0280     }
0281 
0282     // if worst hit was found, refit without worst hit and select if considerably better than original fit.
0283     // Can also use brute force: refit all n-1 hit segments and choose one over the n hit if considerably "better"
0284 
0285     std::vector<CSCRecHit2D> buffer;
0286     std::vector<std::vector<CSCRecHit2D> > reduced_segments;
0287     std::vector<CSCRecHit2D> theseRecHits = (*it).specificRecHits();
0288     float best_red_seg_prob = 0.0;
0289     // usefor chi2 1 diff   float best_red_seg_prob = 99999.;
0290     buffer.clear();
0291 
0292     if (ChiSquaredProbability((double)(*it).chi2(), (double)((2 * (*it).nRecHits()) - 4)) < chi2ndfProbMin) {
0293       buffer = theseRecHits;
0294 
0295       // Dirty switch: here one can select to refit all possible subsets or just the one without the
0296       // tagged worst hit:
0297       if (use_brute_force) {  // Brute force method: loop over all possible segments:
0298         for (size_t bi = 0; bi < buffer.size(); ++bi) {
0299           reduced_segments.push_back(buffer);
0300           reduced_segments[bi].erase(reduced_segments[bi].begin() + (bi), reduced_segments[bi].begin() + (bi + 1));
0301         }
0302       } else {  // More elegant but still biased: erase only worst hit
0303         // Comment: There is not a very strong correlation of the worst hit with the one that one should remove...
0304         if (hit_nr_worst >= 0 && hit_nr_worst <= int(buffer.size())) {
0305           // fill segment in buffer, delete worst hit
0306           buffer.erase(buffer.begin() + (hit_nr_worst), buffer.begin() + (hit_nr_worst + 1));
0307           reduced_segments.push_back(buffer);
0308         } else {
0309           // only fill segment in array, do not delete anything
0310           reduced_segments.push_back(buffer);
0311         }
0312       }
0313     }
0314 
0315     // Loop over the subsegments and fit (only one segment if "use_brute_force" is false):
0316     for (size_t iSegment = 0; iSegment < reduced_segments.size(); ++iSegment) {
0317       // loop over hits on given segment and push pointers to hits into protosegment
0318       protoSegment.clear();
0319       for (size_t m = 0; m < reduced_segments[iSegment].size(); ++m) {
0320         protoSegment.push_back(&reduced_segments[iSegment][m]);
0321       }
0322 
0323       // Create fitter object
0324       CSCCondSegFit* segfit = new CSCCondSegFit(pset(), chamber(), protoSegment);
0325       condpass1 = false;
0326       condpass2 = false;
0327       segfit->setScaleXError(1.0);
0328       segfit->fit(condpass1, condpass2);
0329 
0330       // Attempt to handle numerical instability of the fit;
0331       // The same as in the build method;
0332       // Preprune is not applied;
0333       if (adjustCovariance()) {
0334         if (segfit->chi2() / segfit->ndof() > chi2Norm_3D_) {
0335           condpass1 = true;
0336           segfit->fit(condpass1, condpass2);
0337         }
0338         if ((segfit->scaleXError() < 1.00005) && (segfit->chi2() / segfit->ndof() > chi2Norm_3D_)) {
0339           condpass2 = true;
0340           segfit->fit(condpass1, condpass2);
0341         }
0342       }
0343 
0344       // calculate error matrix
0345       //      AlgebraicSymMatrix temp2 = segfit->covarianceMatrix();
0346 
0347       // build an actual segment
0348       CSCSegment temp(
0349           protoSegment, segfit->intercept(), segfit->localdir(), segfit->covarianceMatrix(), segfit->chi2());
0350 
0351       // and finished with this fit
0352       delete segfit;
0353 
0354       // n-hit segment is (*it)
0355       // (n-1)-hit segment is temp
0356       // replace n-hit segment with (n-1)-hit segment if segment probability is BPMinImprovement better
0357       double oldchi2 = (*it).chi2();
0358       double olddof = 2 * (*it).nRecHits() - 4;
0359       double newchi2 = temp.chi2();
0360       double newdof = 2 * temp.nRecHits() - 4;
0361       if ((ChiSquaredProbability(oldchi2, olddof) < (1. / BPMinImprovement) * ChiSquaredProbability(newchi2, newdof)) &&
0362           (ChiSquaredProbability(newchi2, newdof) > best_red_seg_prob) &&
0363           (ChiSquaredProbability(newchi2, newdof) > 1e-10)) {
0364         best_red_seg_prob = ChiSquaredProbability(newchi2, newdof);
0365         // The (n-1)- segment is better than the n-hit segment.
0366         // If it has at least  minHitsPerSegment  hits replace the n-hit segment
0367         // with this  better (n-1)-hit segment:
0368         if (temp.nRecHits() >= minHitsPerSegment) {
0369           (*it) = temp;
0370         }
0371       }
0372     }  // end of loop over subsegments (iSegment)
0373 
0374   }  // end loop over segments (it)
0375 
0376   return segments;
0377 }
0378 
0379 // ********************************************************************;
0380 std::vector<std::vector<const CSCRecHit2D*> > CSCSegAlgoST::clusterHits(const CSCChamber* aChamber,
0381                                                                         const ChamberHitContainer& rechits) {
0382   theChamber = aChamber;
0383 
0384   std::vector<ChamberHitContainer> rechits_clusters;  // this is a collection of groups of rechits
0385   //   const float dXclus_box_cut       = 4.; // seems to work reasonably 070116
0386   //   const float dYclus_box_cut       = 8.; // seems to work reasonably 070116
0387 
0388   //float dXclus = 0.0;
0389   //float dYclus = 0.0;
0390   float dXclus_box = 0.0;
0391   float dYclus_box = 0.0;
0392 
0393   std::vector<const CSCRecHit2D*> temp;
0394 
0395   std::vector<ChamberHitContainer> seeds;
0396 
0397   std::vector<float> running_meanX;
0398   std::vector<float> running_meanY;
0399 
0400   std::vector<float> seed_minX;
0401   std::vector<float> seed_maxX;
0402   std::vector<float> seed_minY;
0403   std::vector<float> seed_maxY;
0404 
0405   //std::cout<<"*************************************************************"<<std::endl;
0406   //std::cout<<"Called clusterHits in Chamber "<< theChamber->specs()->chamberTypeName()<<std::endl;
0407   //std::cout<<"*************************************************************"<<std::endl;
0408 
0409   // split rechits into subvectors and return vector of vectors:
0410   // Loop over rechits
0411   // Create one seed per hit
0412   for (unsigned int i = 0; i < rechits.size(); ++i) {
0413     temp.clear();
0414 
0415     temp.push_back(rechits[i]);
0416 
0417     seeds.push_back(temp);
0418 
0419     // First added hit in seed defines the mean to which the next hit is compared
0420     // for this seed.
0421 
0422     running_meanX.push_back(rechits[i]->localPosition().x());
0423     running_meanY.push_back(rechits[i]->localPosition().y());
0424 
0425     // set min/max X and Y for box containing the hits in the precluster:
0426     seed_minX.push_back(rechits[i]->localPosition().x());
0427     seed_maxX.push_back(rechits[i]->localPosition().x());
0428     seed_minY.push_back(rechits[i]->localPosition().y());
0429     seed_maxY.push_back(rechits[i]->localPosition().y());
0430   }
0431 
0432   // merge clusters that are too close
0433   // measure distance between final "running mean"
0434   for (size_t NNN = 0; NNN < seeds.size(); ++NNN) {
0435     for (size_t MMM = NNN + 1; MMM < seeds.size(); ++MMM) {
0436       if (running_meanX[MMM] == 999999. || running_meanX[NNN] == 999999.) {
0437         LogTrace("CSCSegment|CSC")
0438             << "[CSCSegAlgoST::clusterHits] ALARM! Skipping used seeds, this should not happen - inform CSC DPG";
0439         //  std::cout<<"We should never see this line now!!!"<<std::endl;
0440         continue;  //skip seeds that have been used
0441       }
0442 
0443       // calculate cut criteria for simple running mean distance cut:
0444       //dXclus = fabs(running_meanX[NNN] - running_meanX[MMM]);
0445       //dYclus = fabs(running_meanY[NNN] - running_meanY[MMM]);
0446 
0447       // calculate minmal distance between precluster boxes containing the hits:
0448       if (running_meanX[NNN] > running_meanX[MMM])
0449         dXclus_box = seed_minX[NNN] - seed_maxX[MMM];
0450       else
0451         dXclus_box = seed_minX[MMM] - seed_maxX[NNN];
0452       if (running_meanY[NNN] > running_meanY[MMM])
0453         dYclus_box = seed_minY[NNN] - seed_maxY[MMM];
0454       else
0455         dYclus_box = seed_minY[MMM] - seed_maxY[NNN];
0456 
0457       if (dXclus_box < dXclusBoxMax && dYclus_box < dYclusBoxMax) {
0458         // merge clusters!
0459         // merge by adding seed NNN to seed MMM and erasing seed NNN
0460 
0461         // calculate running mean for the merged seed:
0462         running_meanX[MMM] = (running_meanX[NNN] * seeds[NNN].size() + running_meanX[MMM] * seeds[MMM].size()) /
0463                              (seeds[NNN].size() + seeds[MMM].size());
0464         running_meanY[MMM] = (running_meanY[NNN] * seeds[NNN].size() + running_meanY[MMM] * seeds[MMM].size()) /
0465                              (seeds[NNN].size() + seeds[MMM].size());
0466 
0467         // update min/max X and Y for box containing the hits in the merged cluster:
0468         if (seed_minX[NNN] <= seed_minX[MMM])
0469           seed_minX[MMM] = seed_minX[NNN];
0470         if (seed_maxX[NNN] > seed_maxX[MMM])
0471           seed_maxX[MMM] = seed_maxX[NNN];
0472         if (seed_minY[NNN] <= seed_minY[MMM])
0473           seed_minY[MMM] = seed_minY[NNN];
0474         if (seed_maxY[NNN] > seed_maxY[MMM])
0475           seed_maxY[MMM] = seed_maxY[NNN];
0476 
0477         // add seed NNN to MMM (lower to larger number)
0478         seeds[MMM].insert(seeds[MMM].end(), seeds[NNN].begin(), seeds[NNN].end());
0479 
0480         // mark seed NNN as used (at the moment just set running mean to 999999.)
0481         running_meanX[NNN] = 999999.;
0482         running_meanY[NNN] = 999999.;
0483         // we have merged a seed (NNN) to the highter seed (MMM) - need to contimue to
0484         // next seed (NNN+1)
0485         break;
0486       }
0487     }
0488   }
0489 
0490   // hand over the final seeds to the output
0491   // would be more elegant if we could do the above step with
0492   // erasing the merged ones, rather than the
0493   for (size_t NNN = 0; NNN < seeds.size(); ++NNN) {
0494     if (running_meanX[NNN] == 999999.)
0495       continue;  //skip seeds that have been marked as used up in merging
0496     rechits_clusters.push_back(seeds[NNN]);
0497   }
0498 
0499   //***************************************************************
0500 
0501   return rechits_clusters;
0502 }
0503 
0504 std::vector<std::vector<const CSCRecHit2D*> > CSCSegAlgoST::chainHits(const CSCChamber* aChamber,
0505                                                                       const ChamberHitContainer& rechits) {
0506   std::vector<ChamberHitContainer> rechits_chains;  // this is a collection of groups of rechits
0507 
0508   std::vector<const CSCRecHit2D*> temp;
0509 
0510   std::vector<ChamberHitContainer> seeds;
0511 
0512   std::vector<bool> usedCluster;
0513 
0514   // split rechits into subvectors and return vector of vectors:
0515   // Loop over rechits
0516   // Create one seed per hit
0517   //std::cout<<" rechits.size() = "<<rechits.size()<<std::endl;
0518   for (unsigned int i = 0; i < rechits.size(); ++i) {
0519     temp.clear();
0520     temp.push_back(rechits[i]);
0521     seeds.push_back(temp);
0522     usedCluster.push_back(false);
0523   }
0524   // Only ME1/1A can have ganged strips so no need to test name
0525   bool gangedME11a = false;
0526   if (("ME1/a" == aChamber->specs()->chamberTypeName()) && aChamber->specs()->gangedStrips()) {
0527     //  if ( aChamber->specs()->gangedStrips() ){
0528     gangedME11a = true;
0529   }
0530   // merge chains that are too close ("touch" each other)
0531   for (size_t NNN = 0; NNN < seeds.size(); ++NNN) {
0532     for (size_t MMM = NNN + 1; MMM < seeds.size(); ++MMM) {
0533       if (usedCluster[MMM] || usedCluster[NNN]) {
0534         continue;
0535       }
0536       // all is in the way we define "good";
0537       // try not to "cluster" the hits but to "chain" them;
0538       // it does the clustering but also does a better job
0539       // for inclined tracks (not clustering them together;
0540       // crossed tracks would be still clustered together)
0541       // 22.12.09: In fact it is not much more different
0542       // than the "clustering", we just introduce another
0543       // variable in the game - Z. And it makes sense
0544       // to re-introduce Y (or actually wire group mumber)
0545       // in a similar way as for the strip number - see
0546       // the code below.
0547       bool goodToMerge = isGoodToMerge(gangedME11a, seeds[NNN], seeds[MMM]);
0548       if (goodToMerge) {
0549         // merge chains!
0550         // merge by adding seed NNN to seed MMM and erasing seed NNN
0551 
0552         // add seed NNN to MMM (lower to larger number)
0553         seeds[MMM].insert(seeds[MMM].end(), seeds[NNN].begin(), seeds[NNN].end());
0554 
0555         // mark seed NNN as used
0556         usedCluster[NNN] = true;
0557         // we have merged a seed (NNN) to the highter seed (MMM) - need to contimue to
0558         // next seed (NNN+1)
0559         break;
0560       }
0561     }
0562   }
0563 
0564   // hand over the final seeds to the output
0565   // would be more elegant if we could do the above step with
0566   // erasing the merged ones, rather than the
0567 
0568   for (size_t NNN = 0; NNN < seeds.size(); ++NNN) {
0569     if (usedCluster[NNN])
0570       continue;  //skip seeds that have been marked as used up in merging
0571     rechits_chains.push_back(seeds[NNN]);
0572   }
0573 
0574   //***************************************************************
0575 
0576   return rechits_chains;
0577 }
0578 
0579 bool CSCSegAlgoST::isGoodToMerge(bool gangedME11a, ChamberHitContainer& newChain, ChamberHitContainer& oldChain) {
0580   for (size_t iRH_new = 0; iRH_new < newChain.size(); ++iRH_new) {
0581     int layer_new = newChain[iRH_new]->cscDetId().layer() - 1;
0582     int middleStrip_new = newChain[iRH_new]->nStrips() / 2;
0583     int centralStrip_new = newChain[iRH_new]->channels(middleStrip_new);
0584     int centralWire_new = newChain[iRH_new]->hitWire();
0585     bool layerRequirementOK = false;
0586     bool stripRequirementOK = false;
0587     bool wireRequirementOK = false;
0588     bool goodToMerge = false;
0589     for (size_t iRH_old = 0; iRH_old < oldChain.size(); ++iRH_old) {
0590       int layer_old = oldChain[iRH_old]->cscDetId().layer() - 1;
0591       int middleStrip_old = oldChain[iRH_old]->nStrips() / 2;
0592       int centralStrip_old = oldChain[iRH_old]->channels(middleStrip_old);
0593       int centralWire_old = oldChain[iRH_old]->hitWire();
0594 
0595       // to be chained, two hits need to be in neighbouring layers...
0596       // or better allow few missing layers (upto 3 to avoid inefficiencies);
0597       // however we'll not make an angle correction because it
0598       // worsen the situation in some of the "regular" cases
0599       // (not making the correction means that the conditions for
0600       // forming a cluster are different if we have missing layers -
0601       // this could affect events at the boundaries )
0602       if (layer_new == layer_old + 1 || layer_new == layer_old - 1 || layer_new == layer_old + 2 ||
0603           layer_new == layer_old - 2 || layer_new == layer_old + 3 || layer_new == layer_old - 3 ||
0604           layer_new == layer_old + 4 || layer_new == layer_old - 4) {
0605         layerRequirementOK = true;
0606       }
0607       int allStrips = 48;
0608       //to be chained, two hits need to be "close" in strip number (can do it in phi
0609       // but it doesn't really matter); let "close" means upto 2 strips (3?) -
0610       // this is more compared to what CLCT readout patterns allow
0611       if (centralStrip_new == centralStrip_old || centralStrip_new == centralStrip_old + 1 ||
0612           centralStrip_new == centralStrip_old - 1 || centralStrip_new == centralStrip_old + 2 ||
0613           centralStrip_new == centralStrip_old - 2) {
0614         stripRequirementOK = true;
0615       }
0616       // same for wires (and ALCT patterns)
0617       if (centralWire_new == centralWire_old || centralWire_new == centralWire_old + 1 ||
0618           centralWire_new == centralWire_old - 1 || centralWire_new == centralWire_old + 2 ||
0619           centralWire_new == centralWire_old - 2) {
0620         wireRequirementOK = true;
0621       }
0622 
0623       if (gangedME11a) {
0624         if (centralStrip_new == centralStrip_old + 1 - allStrips ||
0625             centralStrip_new == centralStrip_old - 1 - allStrips ||
0626             centralStrip_new == centralStrip_old + 2 - allStrips ||
0627             centralStrip_new == centralStrip_old - 2 - allStrips ||
0628             centralStrip_new == centralStrip_old + 1 + allStrips ||
0629             centralStrip_new == centralStrip_old - 1 + allStrips ||
0630             centralStrip_new == centralStrip_old + 2 + allStrips ||
0631             centralStrip_new == centralStrip_old - 2 + allStrips) {
0632           stripRequirementOK = true;
0633         }
0634       }
0635       if (layerRequirementOK && stripRequirementOK && wireRequirementOK) {
0636         goodToMerge = true;
0637         return goodToMerge;
0638       }
0639     }
0640   }
0641   return false;
0642 }
0643 
0644 double CSCSegAlgoST::theWeight(
0645     double coordinate_1, double coordinate_2, double coordinate_3, float layer_1, float layer_2, float layer_3) {
0646   double sub_weight = 0;
0647   sub_weight = fabs(((coordinate_2 - coordinate_3) / (layer_2 - layer_3)) -
0648                     ((coordinate_1 - coordinate_2) / (layer_1 - layer_2)));
0649   return sub_weight;
0650 }
0651 
0652 /* 
0653  * This algorithm uses a Minimum Spanning Tree (ST) approach to build
0654  * endcap muon track segments from the rechits in the 6 layers of a CSC.
0655  */
0656 std::vector<CSCSegment> CSCSegAlgoST::buildSegments(const ChamberHitContainer& rechits) {
0657   // Clear buffer for segment vector
0658   std::vector<CSCSegment> segmentInChamber;
0659   segmentInChamber.clear();  // list of final segments
0660 
0661   // CSC Ring;
0662   unsigned int thering = 999;
0663   unsigned int thestation = 999;
0664   //unsigned int thecham    = 999;
0665 
0666   std::vector<int> hits_onLayerNumber(6);
0667 
0668   unsigned int UpperLimit = maxRecHitsInCluster;
0669   if (int(rechits.size()) < minHitsPerSegment)
0670     return segmentInChamber;
0671 
0672   for (int iarray = 0; iarray < 6; ++iarray) {  // magic number 6: number of layers in CSC chamber - not gonna change :)
0673     PAhits_onLayer[iarray].clear();
0674     hits_onLayerNumber[iarray] = 0;
0675   }
0676 
0677   chosen_Psegments.clear();
0678   chosen_weight_A.clear();
0679 
0680   Psegments.clear();
0681   Psegments_noLx.clear();
0682   Psegments_noL1.clear();
0683   Psegments_noL2.clear();
0684   Psegments_noL3.clear();
0685   Psegments_noL4.clear();
0686   Psegments_noL5.clear();
0687   Psegments_noL6.clear();
0688 
0689   Psegments_hits.clear();
0690 
0691   weight_A.clear();
0692   weight_noLx_A.clear();
0693   weight_noL1_A.clear();
0694   weight_noL2_A.clear();
0695   weight_noL3_A.clear();
0696   weight_noL4_A.clear();
0697   weight_noL5_A.clear();
0698   weight_noL6_A.clear();
0699 
0700   weight_B.clear();
0701   weight_noL1_B.clear();
0702   weight_noL2_B.clear();
0703   weight_noL3_B.clear();
0704   weight_noL4_B.clear();
0705   weight_noL5_B.clear();
0706   weight_noL6_B.clear();
0707 
0708   curv_A.clear();
0709   curv_noL1_A.clear();
0710   curv_noL2_A.clear();
0711   curv_noL3_A.clear();
0712   curv_noL4_A.clear();
0713   curv_noL5_A.clear();
0714   curv_noL6_A.clear();
0715 
0716   // definition of middle layer for n-hit segment
0717   int midlayer_pointer[6] = {0, 0, 2, 3, 3, 4};
0718 
0719   // int n_layers_missed_tot = 0;
0720   int n_layers_occupied_tot = 0;
0721   int n_layers_processed = 0;
0722 
0723   float min_weight_A = 99999.9;
0724   float min_weight_noLx_A = 99999.9;
0725 
0726   //float best_weight_B = -1.;
0727   //float best_weight_noLx_B = -1.;
0728 
0729   //float best_curv_A = -1.;
0730   //float best_curv_noLx_A = -1.;
0731 
0732   int best_pseg = -1;
0733   int best_noLx_pseg = -1;
0734   int best_Layer_noLx = -1;
0735 
0736   //************************************************************************;
0737   //***   Start segment building   *****************************************;
0738   //************************************************************************;
0739 
0740   // Determine how many layers with hits we have
0741   // Fill all hits into the layer hit container:
0742 
0743   // Have 2 standard arrays: one giving the number of hits per layer.
0744   // The other the corresponding hits.
0745 
0746   // Loop all available hits, count hits per layer and fill the hits into array by layer
0747   for (size_t M = 0; M < rechits.size(); ++M) {
0748     // add hits to array per layer and count hits per layer:
0749     hits_onLayerNumber[rechits[M]->cscDetId().layer() - 1] += 1;
0750     if (hits_onLayerNumber[rechits[M]->cscDetId().layer() - 1] == 1)
0751       n_layers_occupied_tot += 1;
0752     // add hits to vector in array
0753     PAhits_onLayer[rechits[M]->cscDetId().layer() - 1].push_back(rechits[M]);
0754   }
0755 
0756   // We have now counted the hits per layer and filled pointers to the hits into an array
0757 
0758   int tothits = 0;
0759   int maxhits = 0;
0760   int nexthits = 0;
0761   int maxlayer = -1;
0762   int nextlayer = -1;
0763 
0764   for (size_t i = 0; i < hits_onLayerNumber.size(); ++i) {
0765     //std::cout<<"We have "<<hits_onLayerNumber[i]<<" hits on layer "<<i+1<<std::endl;
0766     tothits += hits_onLayerNumber[i];
0767     if (hits_onLayerNumber[i] > maxhits) {
0768       nextlayer = maxlayer;
0769       nexthits = maxhits;
0770       maxlayer = i;
0771       maxhits = hits_onLayerNumber[i];
0772     } else if (hits_onLayerNumber[i] > nexthits) {
0773       nextlayer = i;
0774       nexthits = hits_onLayerNumber[i];
0775     }
0776   }
0777 
0778   if (tothits > (int)UpperLimit) {
0779     if (n_layers_occupied_tot > 4) {
0780       tothits = tothits - hits_onLayerNumber[maxlayer];
0781       n_layers_occupied_tot = n_layers_occupied_tot - 1;
0782       PAhits_onLayer[maxlayer].clear();
0783       hits_onLayerNumber[maxlayer] = 0;
0784     }
0785   }
0786 
0787   if (tothits > (int)UpperLimit) {
0788     if (n_layers_occupied_tot > 4) {
0789       tothits = tothits - hits_onLayerNumber[nextlayer];
0790       n_layers_occupied_tot = n_layers_occupied_tot - 1;
0791       PAhits_onLayer[nextlayer].clear();
0792       hits_onLayerNumber[nextlayer] = 0;
0793     }
0794   }
0795 
0796   if (tothits > (int)UpperLimit) {
0797     //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0798     // Showering muon - returns nothing if chi2 == -1 (see comment in SegAlgoShowering)
0799     //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
0800     if (useShowering) {
0801       CSCSegment segShower = showering_->showerSeg(theChamber, rechits);
0802       if (debug)
0803         dumpSegment(segShower);
0804       // Make sure have at least 3 hits...
0805       if (segShower.nRecHits() < 3)
0806         return segmentInChamber;
0807       if (segShower.chi2() == -1)
0808         return segmentInChamber;
0809       segmentInChamber.push_back(segShower);
0810       return segmentInChamber;
0811     } else {
0812       //  LogTrace("CSCSegment|CSC") << "[CSCSegAlgoST::buildSegments] No. of rechits in the cluster/chamber > "
0813       //  << UpperLimit << " ... Segment finding in the cluster/chamber canceled!";
0814       CSCDetId id = rechits[0]->cscDetId();
0815       edm::LogVerbatim("CSCSegment|CSC") << "[CSCSegAlgoST::buildSegments] No. of rechits in the cluster/chamber > "
0816                                          << UpperLimit << " ... Segment finding canceled in CSC" << id;
0817       return segmentInChamber;
0818     }
0819   }
0820 
0821   // Find out which station, ring and chamber we are in
0822   // Used to choose station/ring dependant y-weight cuts
0823 
0824   if (!rechits.empty()) {
0825     thering = rechits[0]->cscDetId().ring();
0826     thestation = rechits[0]->cscDetId().station();
0827     //thecham = rechits[0]->cscDetId().chamber();
0828   }
0829 
0830   // std::cout<<"We are in Station/ring/chamber: "<<thestation <<" "<< thering<<" "<< thecham<<std::endl;
0831 
0832   // Cut-off parameter - don't reconstruct segments with less than X hits
0833   if (n_layers_occupied_tot < minHitsPerSegment) {
0834     return segmentInChamber;
0835   }
0836 
0837   // Start building all possible hit combinations:
0838 
0839   // loop over the six chamber layers and form segment candidates from the available hits:
0840 
0841   for (int layer = 0; layer < 6; ++layer) {
0842     // *****************************************************************
0843     // *** Set missed layer counter here (not currently implemented) ***
0844     // *****************************************************************
0845     // if( PAhits_onLayer[layer].size() == 0 ) {
0846     //   n_layers_missed_tot += 1;
0847     // }
0848 
0849     if (!PAhits_onLayer[layer].empty()) {
0850       n_layers_processed += 1;
0851     }
0852 
0853     // Save the size of the protosegment before hits were added on the current layer
0854     int orig_number_of_psegs = Psegments.size();
0855     int orig_number_of_noL1_psegs = Psegments_noL1.size();
0856     int orig_number_of_noL2_psegs = Psegments_noL2.size();
0857     int orig_number_of_noL3_psegs = Psegments_noL3.size();
0858     int orig_number_of_noL4_psegs = Psegments_noL4.size();
0859     int orig_number_of_noL5_psegs = Psegments_noL5.size();
0860     int orig_number_of_noL6_psegs = Psegments_noL6.size();
0861 
0862     // loop over the hits on the layer and initiate protosegments or add hits to protosegments
0863     for (int hit = 0; hit < int(PAhits_onLayer[layer].size()); ++hit) {  // loop all hits on the Layer number "layer"
0864 
0865       // create protosegments from all hits on the first layer with hits
0866       if (orig_number_of_psegs == 0) {  // would be faster to turn this around - ask for "orig_number_of_psegs != 0"
0867 
0868         Psegments_hits.push_back(PAhits_onLayer[layer][hit]);
0869 
0870         Psegments.push_back(Psegments_hits);
0871         Psegments_noL6.push_back(Psegments_hits);
0872         Psegments_noL5.push_back(Psegments_hits);
0873         Psegments_noL4.push_back(Psegments_hits);
0874         Psegments_noL3.push_back(Psegments_hits);
0875         Psegments_noL2.push_back(Psegments_hits);
0876 
0877         // Initialize weights corresponding to this segment for first hit (with 0)
0878 
0879         curv_A.push_back(0.0);
0880         curv_noL6_A.push_back(0.0);
0881         curv_noL5_A.push_back(0.0);
0882         curv_noL4_A.push_back(0.0);
0883         curv_noL3_A.push_back(0.0);
0884         curv_noL2_A.push_back(0.0);
0885 
0886         weight_A.push_back(0.0);
0887         weight_noL6_A.push_back(0.0);
0888         weight_noL5_A.push_back(0.0);
0889         weight_noL4_A.push_back(0.0);
0890         weight_noL3_A.push_back(0.0);
0891         weight_noL2_A.push_back(0.0);
0892 
0893         weight_B.push_back(0.0);
0894         weight_noL6_B.push_back(0.0);
0895         weight_noL5_B.push_back(0.0);
0896         weight_noL4_B.push_back(0.0);
0897         weight_noL3_B.push_back(0.0);
0898         weight_noL2_B.push_back(0.0);
0899 
0900         // reset array for next hit on next layer
0901         Psegments_hits.clear();
0902       } else {
0903         if (orig_number_of_noL1_psegs == 0) {
0904           Psegments_hits.push_back(PAhits_onLayer[layer][hit]);
0905 
0906           Psegments_noL1.push_back(Psegments_hits);
0907 
0908           // Initialize weight corresponding to this segment for first hit (with 0)
0909 
0910           curv_noL1_A.push_back(0.0);
0911 
0912           weight_noL1_A.push_back(0.0);
0913 
0914           weight_noL1_B.push_back(0.0);
0915 
0916           // reset array for next hit on next layer
0917           Psegments_hits.clear();
0918         }
0919 
0920         // loop over the protosegments and create a new protosegments for each hit-1 on this layer
0921 
0922         for (int pseg = 0; pseg < orig_number_of_psegs; ++pseg) {
0923           int pseg_pos = (pseg) + ((hit)*orig_number_of_psegs);
0924           int pseg_noL1_pos = (pseg) + ((hit)*orig_number_of_noL1_psegs);
0925           int pseg_noL2_pos = (pseg) + ((hit)*orig_number_of_noL2_psegs);
0926           int pseg_noL3_pos = (pseg) + ((hit)*orig_number_of_noL3_psegs);
0927           int pseg_noL4_pos = (pseg) + ((hit)*orig_number_of_noL4_psegs);
0928           int pseg_noL5_pos = (pseg) + ((hit)*orig_number_of_noL5_psegs);
0929           int pseg_noL6_pos = (pseg) + ((hit)*orig_number_of_noL6_psegs);
0930 
0931           // - Loop all psegs.
0932           // - If not last hit, clone  existing protosegments  (PAhits_onLayer[layer].size()-1) times
0933           // - then add the new hits
0934 
0935           if (!(hit == int(PAhits_onLayer[layer].size() -
0936                            1))) {  // not the last hit - prepare (copy) new protosegments for the following hits
0937             // clone psegs (to add next hits or last hit on layer):
0938 
0939             Psegments.push_back(Psegments[pseg_pos]);
0940             if (n_layers_processed != 2 && pseg < orig_number_of_noL1_psegs)
0941               Psegments_noL1.push_back(Psegments_noL1[pseg_noL1_pos]);
0942             if (n_layers_processed != 2 && pseg < orig_number_of_noL2_psegs)
0943               Psegments_noL2.push_back(Psegments_noL2[pseg_noL2_pos]);
0944             if (n_layers_processed != 3 && pseg < orig_number_of_noL3_psegs)
0945               Psegments_noL3.push_back(Psegments_noL3[pseg_noL3_pos]);
0946             if (n_layers_processed != 4 && pseg < orig_number_of_noL4_psegs)
0947               Psegments_noL4.push_back(Psegments_noL4[pseg_noL4_pos]);
0948             if (n_layers_processed != 5 && pseg < orig_number_of_noL5_psegs)
0949               Psegments_noL5.push_back(Psegments_noL5[pseg_noL5_pos]);
0950             if (n_layers_processed != 6 && pseg < orig_number_of_noL6_psegs)
0951               Psegments_noL6.push_back(Psegments_noL6[pseg_noL6_pos]);
0952             // clone weight corresponding to this segment too
0953             weight_A.push_back(weight_A[pseg_pos]);
0954             if (n_layers_processed != 2 && pseg < orig_number_of_noL1_psegs)
0955               weight_noL1_A.push_back(weight_noL1_A[pseg_noL1_pos]);
0956             if (n_layers_processed != 2 && pseg < orig_number_of_noL2_psegs)
0957               weight_noL2_A.push_back(weight_noL2_A[pseg_noL2_pos]);
0958             if (n_layers_processed != 3 && pseg < orig_number_of_noL3_psegs)
0959               weight_noL3_A.push_back(weight_noL3_A[pseg_noL3_pos]);
0960             if (n_layers_processed != 4 && pseg < orig_number_of_noL4_psegs)
0961               weight_noL4_A.push_back(weight_noL4_A[pseg_noL4_pos]);
0962             if (n_layers_processed != 5 && pseg < orig_number_of_noL5_psegs)
0963               weight_noL5_A.push_back(weight_noL5_A[pseg_noL5_pos]);
0964             if (n_layers_processed != 6 && pseg < orig_number_of_noL6_psegs)
0965               weight_noL6_A.push_back(weight_noL6_A[pseg_noL6_pos]);
0966             // clone curvature variable corresponding to this segment too
0967             curv_A.push_back(curv_A[pseg_pos]);
0968             if (n_layers_processed != 2 && pseg < orig_number_of_noL1_psegs)
0969               curv_noL1_A.push_back(curv_noL1_A[pseg_noL1_pos]);
0970             if (n_layers_processed != 2 && pseg < orig_number_of_noL2_psegs)
0971               curv_noL2_A.push_back(curv_noL2_A[pseg_noL2_pos]);
0972             if (n_layers_processed != 3 && pseg < orig_number_of_noL3_psegs)
0973               curv_noL3_A.push_back(curv_noL3_A[pseg_noL3_pos]);
0974             if (n_layers_processed != 4 && pseg < orig_number_of_noL4_psegs)
0975               curv_noL4_A.push_back(curv_noL4_A[pseg_noL4_pos]);
0976             if (n_layers_processed != 5 && pseg < orig_number_of_noL5_psegs)
0977               curv_noL5_A.push_back(curv_noL5_A[pseg_noL5_pos]);
0978             if (n_layers_processed != 6 && pseg < orig_number_of_noL6_psegs)
0979               curv_noL6_A.push_back(curv_noL6_A[pseg_noL6_pos]);
0980             // clone "y"-weight corresponding to this segment too
0981             weight_B.push_back(weight_B[pseg_pos]);
0982             if (n_layers_processed != 2 && pseg < orig_number_of_noL1_psegs)
0983               weight_noL1_B.push_back(weight_noL1_B[pseg_noL1_pos]);
0984             if (n_layers_processed != 2 && pseg < orig_number_of_noL2_psegs)
0985               weight_noL2_B.push_back(weight_noL2_B[pseg_noL2_pos]);
0986             if (n_layers_processed != 3 && pseg < orig_number_of_noL3_psegs)
0987               weight_noL3_B.push_back(weight_noL3_B[pseg_noL3_pos]);
0988             if (n_layers_processed != 4 && pseg < orig_number_of_noL4_psegs)
0989               weight_noL4_B.push_back(weight_noL4_B[pseg_noL4_pos]);
0990             if (n_layers_processed != 5 && pseg < orig_number_of_noL5_psegs)
0991               weight_noL5_B.push_back(weight_noL5_B[pseg_noL5_pos]);
0992             if (n_layers_processed != 6 && pseg < orig_number_of_noL6_psegs)
0993               weight_noL6_B.push_back(weight_noL6_B[pseg_noL6_pos]);
0994           }
0995           // add hits to original pseg:
0996           Psegments[pseg_pos].push_back(PAhits_onLayer[layer][hit]);
0997           if (n_layers_processed != 2 && pseg < orig_number_of_noL1_psegs)
0998             Psegments_noL1[pseg_noL1_pos].push_back(PAhits_onLayer[layer][hit]);
0999           if (n_layers_processed != 2 && pseg < orig_number_of_noL2_psegs)
1000             Psegments_noL2[pseg_noL2_pos].push_back(PAhits_onLayer[layer][hit]);
1001           if (n_layers_processed != 3 && pseg < orig_number_of_noL3_psegs)
1002             Psegments_noL3[pseg_noL3_pos].push_back(PAhits_onLayer[layer][hit]);
1003           if (n_layers_processed != 4 && pseg < orig_number_of_noL4_psegs)
1004             Psegments_noL4[pseg_noL4_pos].push_back(PAhits_onLayer[layer][hit]);
1005           if (n_layers_processed != 5 && pseg < orig_number_of_noL5_psegs)
1006             Psegments_noL5[pseg_noL5_pos].push_back(PAhits_onLayer[layer][hit]);
1007           if (n_layers_processed != 6 && pseg < orig_number_of_noL6_psegs)
1008             Psegments_noL6[pseg_noL6_pos].push_back(PAhits_onLayer[layer][hit]);
1009 
1010           // calculate/update the weight (only for >2 hits on psegment):
1011 
1012           if (Psegments[pseg_pos].size() > 2) {
1013             // looks more exciting than it is. Here the weight is calculated. It is the difference in x of the last two and one but the last two hits,
1014             // divided by the distance of the corresponding hits. Please refer to twiki page XXXX or CMS Note YYY (and use layer_distance)
1015 
1016             weight_A[pseg_pos] += theWeight((*(Psegments[pseg_pos].end() - 1))->localPosition().x(),
1017                                             (*(Psegments[pseg_pos].end() - 2))->localPosition().x(),
1018                                             (*(Psegments[pseg_pos].end() - 3))->localPosition().x(),
1019                                             float((*(Psegments[pseg_pos].end() - 1))->cscDetId().layer()),
1020                                             float((*(Psegments[pseg_pos].end() - 2))->cscDetId().layer()),
1021                                             float((*(Psegments[pseg_pos].end() - 3))->cscDetId().layer()));
1022 
1023             weight_B[pseg_pos] += theWeight((*(Psegments[pseg_pos].end() - 1))->localPosition().y(),
1024                                             (*(Psegments[pseg_pos].end() - 2))->localPosition().y(),
1025                                             (*(Psegments[pseg_pos].end() - 3))->localPosition().y(),
1026                                             float((*(Psegments[pseg_pos].end() - 1))->cscDetId().layer()),
1027                                             float((*(Psegments[pseg_pos].end() - 2))->cscDetId().layer()),
1028                                             float((*(Psegments[pseg_pos].end() - 3))->cscDetId().layer()));
1029 
1030             // if we have picked up the last hit go looking for pseg with the lowest (and second lowest?) weight
1031 
1032             if (int(Psegments[pseg_pos].size()) == n_layers_occupied_tot) {
1033               curv_A[pseg_pos] += theWeight(
1034                   (*(Psegments[pseg_pos].end() - 1))->localPosition().x(),
1035                   (*(Psegments[pseg_pos].end() - midlayer_pointer[n_layers_occupied_tot - 1]))->localPosition().x(),
1036                   (*(Psegments[pseg_pos].end() - n_layers_occupied_tot))->localPosition().x(),
1037                   float((*(Psegments[pseg_pos].end() - 1))->cscDetId().layer()),
1038                   float(
1039                       (*(Psegments[pseg_pos].end() - midlayer_pointer[n_layers_occupied_tot - 1]))->cscDetId().layer()),
1040                   float((*(Psegments[pseg_pos].end() - n_layers_occupied_tot))->cscDetId().layer()));
1041 
1042               if (curv_A[pseg_pos] > curvePenaltyThreshold)
1043                 weight_A[pseg_pos] = weight_A[pseg_pos] * curvePenalty;
1044 
1045               if (weight_B[pseg_pos] > a_yweightPenaltyThreshold[thestation][thering])
1046                 weight_A[pseg_pos] = weight_A[pseg_pos] * yweightPenalty;
1047 
1048               if (weight_A[pseg_pos] < min_weight_A) {
1049                 min_weight_A = weight_A[pseg_pos];
1050                 //best_weight_B = weight_B[ pseg_pos ];
1051                 //best_curv_A = curv_A[ pseg_pos ];
1052                 best_pseg = pseg_pos;
1053               }
1054             }
1055 
1056             // alternative: fill map with weight and pseg (which is already ordered)? Seems a very good tool to go looking for segments from.
1057             // As I understand, the segments would be inserted according to their weight, so the list would "automatically" be sorted.
1058           }
1059 
1060           if (n_layers_occupied_tot > 3) {
1061             if (pseg < orig_number_of_noL1_psegs && (n_layers_processed != 2)) {
1062               if ((Psegments_noL1[pseg_noL1_pos].size() > 2)) {
1063                 // looks more exciting than it is. Here the weight is calculated. It is the difference in x of the last two and one but the last two hits,
1064                 // divided by the distance of the corresponding hits. Please refer to twiki page XXXX or CMS Note YYY (and use layer_distance)
1065 
1066                 weight_noL1_A[pseg_noL1_pos] +=
1067                     theWeight((*(Psegments_noL1[pseg_noL1_pos].end() - 1))->localPosition().x(),
1068                               (*(Psegments_noL1[pseg_noL1_pos].end() - 2))->localPosition().x(),
1069                               (*(Psegments_noL1[pseg_noL1_pos].end() - 3))->localPosition().x(),
1070                               float((*(Psegments_noL1[pseg_noL1_pos].end() - 1))->cscDetId().layer()),
1071                               float((*(Psegments_noL1[pseg_noL1_pos].end() - 2))->cscDetId().layer()),
1072                               float((*(Psegments_noL1[pseg_noL1_pos].end() - 3))->cscDetId().layer()));
1073 
1074                 weight_noL1_B[pseg_noL1_pos] +=
1075                     theWeight((*(Psegments_noL1[pseg_noL1_pos].end() - 1))->localPosition().y(),
1076                               (*(Psegments_noL1[pseg_noL1_pos].end() - 2))->localPosition().y(),
1077                               (*(Psegments_noL1[pseg_noL1_pos].end() - 3))->localPosition().y(),
1078                               float((*(Psegments_noL1[pseg_noL1_pos].end() - 1))->cscDetId().layer()),
1079                               float((*(Psegments_noL1[pseg_noL1_pos].end() - 2))->cscDetId().layer()),
1080                               float((*(Psegments_noL1[pseg_noL1_pos].end() - 3))->cscDetId().layer()));
1081 
1082                 //if we have picked up the last hit go looking for pseg with the lowest (and second lowest?) weight
1083 
1084                 if (int(Psegments_noL1[pseg_noL1_pos].size()) == n_layers_occupied_tot - 1) {
1085                   curv_noL1_A[pseg_noL1_pos] += theWeight(
1086                       (*(Psegments_noL1[pseg_noL1_pos].end() - 1))->localPosition().x(),
1087                       (*(Psegments_noL1[pseg_noL1_pos].end() - midlayer_pointer[n_layers_occupied_tot - 2]))
1088                           ->localPosition()
1089                           .x(),
1090                       (*(Psegments_noL1[pseg_noL1_pos].end() - (n_layers_occupied_tot - 1)))->localPosition().x(),
1091                       float((*(Psegments_noL1[pseg_noL1_pos].end() - 1))->cscDetId().layer()),
1092                       float((*(Psegments_noL1[pseg_noL1_pos].end() - midlayer_pointer[n_layers_occupied_tot - 2]))
1093                                 ->cscDetId()
1094                                 .layer()),
1095                       float(
1096                           (*(Psegments_noL1[pseg_noL1_pos].end() - (n_layers_occupied_tot - 1)))->cscDetId().layer()));
1097 
1098                   if (curv_noL1_A[pseg_noL1_pos] > curvePenaltyThreshold)
1099                     weight_noL1_A[pseg_noL1_pos] = weight_noL1_A[pseg_noL1_pos] * curvePenalty;
1100 
1101                   if (weight_noL1_B[pseg_noL1_pos] > a_yweightPenaltyThreshold[thestation][thering])
1102                     weight_noL1_A[pseg_noL1_pos] = weight_noL1_A[pseg_noL1_pos] * yweightPenalty;
1103 
1104                   if (weight_noL1_A[pseg_noL1_pos] < min_weight_noLx_A) {
1105                     min_weight_noLx_A = weight_noL1_A[pseg_noL1_pos];
1106                     //best_weight_noLx_B = weight_noL1_B[ pseg_noL1_pos ];
1107                     //best_curv_noLx_A = curv_noL1_A[ pseg_noL1_pos ];
1108                     best_noLx_pseg = pseg_noL1_pos;
1109                     best_Layer_noLx = 1;
1110                   }
1111                 }
1112 
1113                 // alternative: fill map with weight and pseg (which is already ordered)? Seems a very good tool to go looking for segments from.
1114                 // As I understand, the segments would be inserted according to their weight, so the list would "automatically" be sorted.
1115               }
1116             }
1117           }
1118 
1119           if (n_layers_occupied_tot > 3) {
1120             if (pseg < orig_number_of_noL2_psegs && (n_layers_processed != 2)) {
1121               if ((Psegments_noL2[pseg_noL2_pos].size() > 2)) {
1122                 // looks more exciting than it is. Here the weight is calculated. It is the difference in x of the last two and one but the last two hits,
1123                 // divided by the distance of the corresponding hits. Please refer to twiki page XXXX or CMS Note YYY (and use layer_distance)
1124 
1125                 weight_noL2_A[pseg_noL2_pos] +=
1126                     theWeight((*(Psegments_noL2[pseg_noL2_pos].end() - 1))->localPosition().x(),
1127                               (*(Psegments_noL2[pseg_noL2_pos].end() - 2))->localPosition().x(),
1128                               (*(Psegments_noL2[pseg_noL2_pos].end() - 3))->localPosition().x(),
1129                               float((*(Psegments_noL2[pseg_noL2_pos].end() - 1))->cscDetId().layer()),
1130                               float((*(Psegments_noL2[pseg_noL2_pos].end() - 2))->cscDetId().layer()),
1131                               float((*(Psegments_noL2[pseg_noL2_pos].end() - 3))->cscDetId().layer()));
1132 
1133                 weight_noL2_B[pseg_noL2_pos] +=
1134                     theWeight((*(Psegments_noL2[pseg_noL2_pos].end() - 1))->localPosition().y(),
1135                               (*(Psegments_noL2[pseg_noL2_pos].end() - 2))->localPosition().y(),
1136                               (*(Psegments_noL2[pseg_noL2_pos].end() - 3))->localPosition().y(),
1137                               float((*(Psegments_noL2[pseg_noL2_pos].end() - 1))->cscDetId().layer()),
1138                               float((*(Psegments_noL2[pseg_noL2_pos].end() - 2))->cscDetId().layer()),
1139                               float((*(Psegments_noL2[pseg_noL2_pos].end() - 3))->cscDetId().layer()));
1140 
1141                 //if we have picked up the last hit go looking for pseg with the lowest (and second lowest?) weight
1142 
1143                 if (int(Psegments_noL2[pseg_noL2_pos].size()) == n_layers_occupied_tot - 1) {
1144                   curv_noL2_A[pseg_noL2_pos] += theWeight(
1145                       (*(Psegments_noL2[pseg_noL2_pos].end() - 1))->localPosition().x(),
1146                       (*(Psegments_noL2[pseg_noL2_pos].end() - midlayer_pointer[n_layers_occupied_tot - 2]))
1147                           ->localPosition()
1148                           .x(),
1149                       (*(Psegments_noL2[pseg_noL2_pos].end() - (n_layers_occupied_tot - 1)))->localPosition().x(),
1150                       float((*(Psegments_noL2[pseg_noL2_pos].end() - 1))->cscDetId().layer()),
1151                       float((*(Psegments_noL2[pseg_noL2_pos].end() - midlayer_pointer[n_layers_occupied_tot - 2]))
1152                                 ->cscDetId()
1153                                 .layer()),
1154                       float(
1155                           (*(Psegments_noL2[pseg_noL2_pos].end() - (n_layers_occupied_tot - 1)))->cscDetId().layer()));
1156 
1157                   if (curv_noL2_A[pseg_noL2_pos] > curvePenaltyThreshold)
1158                     weight_noL2_A[pseg_noL2_pos] = weight_noL2_A[pseg_noL2_pos] * curvePenalty;
1159 
1160                   if (weight_noL2_B[pseg_noL2_pos] > a_yweightPenaltyThreshold[thestation][thering])
1161                     weight_noL2_A[pseg_noL2_pos] = weight_noL2_A[pseg_noL2_pos] * yweightPenalty;
1162 
1163                   if (weight_noL2_A[pseg_noL2_pos] < min_weight_noLx_A) {
1164                     min_weight_noLx_A = weight_noL2_A[pseg_noL2_pos];
1165                     //best_weight_noLx_B = weight_noL2_B[ pseg_noL2_pos ];
1166                     //best_curv_noLx_A = curv_noL2_A[ pseg_noL2_pos ];
1167                     best_noLx_pseg = pseg_noL2_pos;
1168                     best_Layer_noLx = 2;
1169                   }
1170                 }
1171 
1172                 // alternative: fill map with weight and pseg (which is already ordered)? Seems a very good tool to go looking for segments from.
1173                 // As I understand, the segments would be inserted according to their weight, so the list would "automatically" be sorted.
1174               }
1175             }
1176           }
1177 
1178           if (n_layers_occupied_tot > 3) {
1179             if (pseg < orig_number_of_noL3_psegs && (n_layers_processed != 3)) {
1180               if ((Psegments_noL3[pseg_noL3_pos].size() > 2)) {
1181                 // looks more exciting than it is. Here the weight is calculated. It is the difference in x of the last two and one but the last two hits,
1182                 // divided by the distance of the corresponding hits. Please refer to twiki page XXXX or CMS Note YYY (and use layer_distance)
1183 
1184                 weight_noL3_A[pseg_noL3_pos] +=
1185                     theWeight((*(Psegments_noL3[pseg_noL3_pos].end() - 1))->localPosition().x(),
1186                               (*(Psegments_noL3[pseg_noL3_pos].end() - 2))->localPosition().x(),
1187                               (*(Psegments_noL3[pseg_noL3_pos].end() - 3))->localPosition().x(),
1188                               float((*(Psegments_noL3[pseg_noL3_pos].end() - 1))->cscDetId().layer()),
1189                               float((*(Psegments_noL3[pseg_noL3_pos].end() - 2))->cscDetId().layer()),
1190                               float((*(Psegments_noL3[pseg_noL3_pos].end() - 3))->cscDetId().layer()));
1191 
1192                 weight_noL3_B[pseg_noL3_pos] +=
1193                     theWeight((*(Psegments_noL3[pseg_noL3_pos].end() - 1))->localPosition().y(),
1194                               (*(Psegments_noL3[pseg_noL3_pos].end() - 2))->localPosition().y(),
1195                               (*(Psegments_noL3[pseg_noL3_pos].end() - 3))->localPosition().y(),
1196                               float((*(Psegments_noL3[pseg_noL3_pos].end() - 1))->cscDetId().layer()),
1197                               float((*(Psegments_noL3[pseg_noL3_pos].end() - 2))->cscDetId().layer()),
1198                               float((*(Psegments_noL3[pseg_noL3_pos].end() - 3))->cscDetId().layer()));
1199 
1200                 //if we have picked up the last hit go looking for pseg with the lowest (and second lowest?) weight
1201 
1202                 if (int(Psegments_noL3[pseg_noL3_pos].size()) == n_layers_occupied_tot - 1) {
1203                   curv_noL3_A[pseg_noL3_pos] += theWeight(
1204                       (*(Psegments_noL3[pseg_noL3_pos].end() - 1))->localPosition().x(),
1205                       (*(Psegments_noL3[pseg_noL3_pos].end() - midlayer_pointer[n_layers_occupied_tot - 2]))
1206                           ->localPosition()
1207                           .x(),
1208                       (*(Psegments_noL3[pseg_noL3_pos].end() - (n_layers_occupied_tot - 1)))->localPosition().x(),
1209                       float((*(Psegments_noL3[pseg_noL3_pos].end() - 1))->cscDetId().layer()),
1210                       float((*(Psegments_noL3[pseg_noL3_pos].end() - midlayer_pointer[n_layers_occupied_tot - 2]))
1211                                 ->cscDetId()
1212                                 .layer()),
1213                       float(
1214                           (*(Psegments_noL3[pseg_noL3_pos].end() - (n_layers_occupied_tot - 1)))->cscDetId().layer()));
1215 
1216                   if (curv_noL3_A[pseg_noL3_pos] > curvePenaltyThreshold)
1217                     weight_noL3_A[pseg_noL3_pos] = weight_noL3_A[pseg_noL3_pos] * curvePenalty;
1218 
1219                   if (weight_noL3_B[pseg_noL3_pos] > a_yweightPenaltyThreshold[thestation][thering])
1220                     weight_noL3_A[pseg_noL3_pos] = weight_noL3_A[pseg_noL3_pos] * yweightPenalty;
1221 
1222                   if (weight_noL3_A[pseg_noL3_pos] < min_weight_noLx_A) {
1223                     min_weight_noLx_A = weight_noL3_A[pseg_noL3_pos];
1224                     //best_weight_noLx_B = weight_noL3_B[ pseg_noL3_pos ];
1225                     //best_curv_noLx_A = curv_noL3_A[ pseg_noL3_pos ];
1226                     best_noLx_pseg = pseg_noL3_pos;
1227                     best_Layer_noLx = 3;
1228                   }
1229                 }
1230 
1231                 // alternative: fill map with weight and pseg (which is already ordered)? Seems a very good tool to go looking for segments from.
1232                 // As I understand, the segments would be inserted according to their weight, so the list would "automatically" be sorted.
1233               }
1234             }
1235           }
1236 
1237           if (n_layers_occupied_tot > 3) {
1238             if (pseg < orig_number_of_noL4_psegs && (n_layers_processed != 4)) {
1239               if ((Psegments_noL4[pseg_noL4_pos].size() > 2)) {
1240                 // looks more exciting than it is. Here the weight is calculated. It is the difference in x of the last two and one but the last two hits,
1241                 // divided by the distance of the corresponding hits. Please refer to twiki page XXXX or CMS Note YYY (and use layer_distance)
1242 
1243                 weight_noL4_A[pseg_noL4_pos] +=
1244                     theWeight((*(Psegments_noL4[pseg_noL4_pos].end() - 1))->localPosition().x(),
1245                               (*(Psegments_noL4[pseg_noL4_pos].end() - 2))->localPosition().x(),
1246                               (*(Psegments_noL4[pseg_noL4_pos].end() - 3))->localPosition().x(),
1247                               float((*(Psegments_noL4[pseg_noL4_pos].end() - 1))->cscDetId().layer()),
1248                               float((*(Psegments_noL4[pseg_noL4_pos].end() - 2))->cscDetId().layer()),
1249                               float((*(Psegments_noL4[pseg_noL4_pos].end() - 3))->cscDetId().layer()));
1250 
1251                 weight_noL4_B[pseg_noL4_pos] +=
1252                     theWeight((*(Psegments_noL4[pseg_noL4_pos].end() - 1))->localPosition().y(),
1253                               (*(Psegments_noL4[pseg_noL4_pos].end() - 2))->localPosition().y(),
1254                               (*(Psegments_noL4[pseg_noL4_pos].end() - 3))->localPosition().y(),
1255                               float((*(Psegments_noL4[pseg_noL4_pos].end() - 1))->cscDetId().layer()),
1256                               float((*(Psegments_noL4[pseg_noL4_pos].end() - 2))->cscDetId().layer()),
1257                               float((*(Psegments_noL4[pseg_noL4_pos].end() - 3))->cscDetId().layer()));
1258 
1259                 //if we have picked up the last hit go looking for pseg with the lowest (and second lowest?) weight
1260 
1261                 if (int(Psegments_noL4[pseg_noL4_pos].size()) == n_layers_occupied_tot - 1) {
1262                   curv_noL4_A[pseg_noL4_pos] += theWeight(
1263                       (*(Psegments_noL4[pseg_noL4_pos].end() - 1))->localPosition().x(),
1264                       (*(Psegments_noL4[pseg_noL4_pos].end() - midlayer_pointer[n_layers_occupied_tot - 2]))
1265                           ->localPosition()
1266                           .x(),
1267                       (*(Psegments_noL4[pseg_noL4_pos].end() - (n_layers_occupied_tot - 1)))->localPosition().x(),
1268                       float((*(Psegments_noL4[pseg_noL4_pos].end() - 1))->cscDetId().layer()),
1269                       float((*(Psegments_noL4[pseg_noL4_pos].end() - midlayer_pointer[n_layers_occupied_tot - 2]))
1270                                 ->cscDetId()
1271                                 .layer()),
1272                       float(
1273                           (*(Psegments_noL4[pseg_noL4_pos].end() - (n_layers_occupied_tot - 1)))->cscDetId().layer()));
1274 
1275                   if (curv_noL4_A[pseg_noL4_pos] > curvePenaltyThreshold)
1276                     weight_noL4_A[pseg_noL4_pos] = weight_noL4_A[pseg_noL4_pos] * curvePenalty;
1277 
1278                   if (weight_noL4_B[pseg_noL4_pos] > a_yweightPenaltyThreshold[thestation][thering])
1279                     weight_noL4_A[pseg_noL4_pos] = weight_noL4_A[pseg_noL4_pos] * yweightPenalty;
1280 
1281                   if (weight_noL4_A[pseg_noL4_pos] < min_weight_noLx_A) {
1282                     min_weight_noLx_A = weight_noL4_A[pseg_noL4_pos];
1283                     //best_weight_noLx_B = weight_noL4_B[ pseg_noL4_pos ];
1284                     //best_curv_noLx_A = curv_noL4_A[ pseg_noL4_pos ];
1285                     best_noLx_pseg = pseg_noL4_pos;
1286                     best_Layer_noLx = 4;
1287                   }
1288                 }
1289 
1290                 // alternative: fill map with weight and pseg (which is already ordered)? Seems a very good tool to go looking for segments from.
1291                 // As I understand, the segments would be inserted according to their weight, so the list would "automatically" be sorted.
1292               }
1293             }
1294           }
1295 
1296           if (n_layers_occupied_tot > 4) {
1297             if (pseg < orig_number_of_noL5_psegs && (n_layers_processed != 5)) {
1298               if ((Psegments_noL5[pseg_noL5_pos].size() > 2)) {
1299                 // looks more exciting than it is. Here the weight is calculated. It is the difference in x of the last two and one but the last two hits,
1300                 // divided by the distance of the corresponding hits. Please refer to twiki page XXXX or CMS Note YYY (and use layer_distance)
1301 
1302                 weight_noL5_A[pseg_noL5_pos] +=
1303                     theWeight((*(Psegments_noL5[pseg_noL5_pos].end() - 1))->localPosition().x(),
1304                               (*(Psegments_noL5[pseg_noL5_pos].end() - 2))->localPosition().x(),
1305                               (*(Psegments_noL5[pseg_noL5_pos].end() - 3))->localPosition().x(),
1306                               float((*(Psegments_noL5[pseg_noL5_pos].end() - 1))->cscDetId().layer()),
1307                               float((*(Psegments_noL5[pseg_noL5_pos].end() - 2))->cscDetId().layer()),
1308                               float((*(Psegments_noL5[pseg_noL5_pos].end() - 3))->cscDetId().layer()));
1309 
1310                 weight_noL5_B[pseg_noL5_pos] +=
1311                     theWeight((*(Psegments_noL5[pseg_noL5_pos].end() - 1))->localPosition().y(),
1312                               (*(Psegments_noL5[pseg_noL5_pos].end() - 2))->localPosition().y(),
1313                               (*(Psegments_noL5[pseg_noL5_pos].end() - 3))->localPosition().y(),
1314                               float((*(Psegments_noL5[pseg_noL5_pos].end() - 1))->cscDetId().layer()),
1315                               float((*(Psegments_noL5[pseg_noL5_pos].end() - 2))->cscDetId().layer()),
1316                               float((*(Psegments_noL5[pseg_noL5_pos].end() - 3))->cscDetId().layer()));
1317 
1318                 //if we have picked up the last hit go looking for pseg with the lowest (and second lowest?) weight
1319 
1320                 if (int(Psegments_noL5[pseg_noL5_pos].size()) == n_layers_occupied_tot - 1) {
1321                   curv_noL5_A[pseg_noL5_pos] += theWeight(
1322                       (*(Psegments_noL5[pseg_noL5_pos].end() - 1))->localPosition().x(),
1323                       (*(Psegments_noL5[pseg_noL5_pos].end() - midlayer_pointer[n_layers_occupied_tot - 2]))
1324                           ->localPosition()
1325                           .x(),
1326                       (*(Psegments_noL5[pseg_noL5_pos].end() - (n_layers_occupied_tot - 1)))->localPosition().x(),
1327                       float((*(Psegments_noL5[pseg_noL5_pos].end() - 1))->cscDetId().layer()),
1328                       float((*(Psegments_noL5[pseg_noL5_pos].end() - midlayer_pointer[n_layers_occupied_tot - 2]))
1329                                 ->cscDetId()
1330                                 .layer()),
1331                       float(
1332                           (*(Psegments_noL5[pseg_noL5_pos].end() - (n_layers_occupied_tot - 1)))->cscDetId().layer()));
1333 
1334                   if (curv_noL5_A[pseg_noL5_pos] > curvePenaltyThreshold)
1335                     weight_noL5_A[pseg_noL5_pos] = weight_noL5_A[pseg_noL5_pos] * curvePenalty;
1336 
1337                   if (weight_noL5_B[pseg_noL5_pos] > a_yweightPenaltyThreshold[thestation][thering])
1338                     weight_noL5_A[pseg_noL5_pos] = weight_noL5_A[pseg_noL5_pos] * yweightPenalty;
1339 
1340                   if (weight_noL5_A[pseg_noL5_pos] < min_weight_noLx_A) {
1341                     min_weight_noLx_A = weight_noL5_A[pseg_noL5_pos];
1342                     //best_weight_noLx_B = weight_noL5_B[ pseg_noL5_pos ];
1343                     //best_curv_noLx_A = curv_noL5_A[ pseg_noL5_pos ];
1344                     best_noLx_pseg = pseg_noL5_pos;
1345                     best_Layer_noLx = 5;
1346                   }
1347                 }
1348 
1349                 // alternative: fill map with weight and pseg (which is already ordered)? Seems a very good tool to go looking for segments from.
1350                 // As I understand, the segments would be inserted according to their weight, so the list would "automatically" be sorted.
1351               }
1352             }
1353           }
1354 
1355           if (n_layers_occupied_tot > 5) {
1356             if (pseg < orig_number_of_noL6_psegs && (n_layers_processed != 6)) {
1357               if ((Psegments_noL6[pseg_noL6_pos].size() > 2)) {
1358                 // looks more exciting than it is. Here the weight is calculated. It is the difference in x of the last two and one but the last two hits,
1359                 // divided by the distance of the corresponding hits. Please refer to twiki page XXXX or CMS Note YYY (and use layer_distance)
1360 
1361                 weight_noL6_A[pseg_noL6_pos] +=
1362                     theWeight((*(Psegments_noL6[pseg_noL6_pos].end() - 1))->localPosition().x(),
1363                               (*(Psegments_noL6[pseg_noL6_pos].end() - 2))->localPosition().x(),
1364                               (*(Psegments_noL6[pseg_noL6_pos].end() - 3))->localPosition().x(),
1365                               float((*(Psegments_noL6[pseg_noL6_pos].end() - 1))->cscDetId().layer()),
1366                               float((*(Psegments_noL6[pseg_noL6_pos].end() - 2))->cscDetId().layer()),
1367                               float((*(Psegments_noL6[pseg_noL6_pos].end() - 3))->cscDetId().layer()));
1368 
1369                 weight_noL6_B[pseg_noL6_pos] +=
1370                     theWeight((*(Psegments_noL6[pseg_noL6_pos].end() - 1))->localPosition().y(),
1371                               (*(Psegments_noL6[pseg_noL6_pos].end() - 2))->localPosition().y(),
1372                               (*(Psegments_noL6[pseg_noL6_pos].end() - 3))->localPosition().y(),
1373                               float((*(Psegments_noL6[pseg_noL6_pos].end() - 1))->cscDetId().layer()),
1374                               float((*(Psegments_noL6[pseg_noL6_pos].end() - 2))->cscDetId().layer()),
1375                               float((*(Psegments_noL6[pseg_noL6_pos].end() - 3))->cscDetId().layer()));
1376 
1377                 //if we have picked up the last hit go looking for pseg with the lowest (and second lowest?) weight
1378 
1379                 if (int(Psegments_noL6[pseg_noL6_pos].size()) == n_layers_occupied_tot - 1) {
1380                   curv_noL6_A[pseg_noL6_pos] += theWeight(
1381                       (*(Psegments_noL6[pseg_noL6_pos].end() - 1))->localPosition().x(),
1382                       (*(Psegments_noL6[pseg_noL6_pos].end() - midlayer_pointer[n_layers_occupied_tot - 2]))
1383                           ->localPosition()
1384                           .x(),
1385                       (*(Psegments_noL6[pseg_noL6_pos].end() - (n_layers_occupied_tot - 1)))->localPosition().x(),
1386                       float((*(Psegments_noL6[pseg_noL6_pos].end() - 1))->cscDetId().layer()),
1387                       float((*(Psegments_noL6[pseg_noL6_pos].end() - midlayer_pointer[n_layers_occupied_tot - 2]))
1388                                 ->cscDetId()
1389                                 .layer()),
1390                       float(
1391                           (*(Psegments_noL6[pseg_noL6_pos].end() - (n_layers_occupied_tot - 1)))->cscDetId().layer()));
1392 
1393                   if (curv_noL6_A[pseg_noL6_pos] > curvePenaltyThreshold)
1394                     weight_noL6_A[pseg_noL6_pos] = weight_noL6_A[pseg_noL6_pos] * curvePenalty;
1395 
1396                   if (weight_noL6_B[pseg_noL6_pos] > a_yweightPenaltyThreshold[thestation][thering])
1397                     weight_noL6_A[pseg_noL6_pos] = weight_noL6_A[pseg_noL6_pos] * yweightPenalty;
1398 
1399                   if (weight_noL6_A[pseg_noL6_pos] < min_weight_noLx_A) {
1400                     min_weight_noLx_A = weight_noL6_A[pseg_noL6_pos];
1401                     //best_weight_noLx_B = weight_noL6_B[ pseg_noL6_pos ];
1402                     //best_curv_noLx_A = curv_noL6_A[ pseg_noL6_pos ];
1403                     best_noLx_pseg = pseg_noL6_pos;
1404                     best_Layer_noLx = 6;
1405                   }
1406                 }
1407 
1408                 // alternative: fill map with weight and pseg (which is already ordered)? Seems a very good tool to go looking for segments from.
1409                 // As I understand, the segments would be inserted according to their weight, so the list would "automatically" be sorted.
1410               }
1411             }
1412           }
1413         }
1414       }
1415     }
1416   }
1417 
1418   //************************************************************************;
1419   //***   End segment building   *******************************************;
1420   //************************************************************************;
1421 
1422   // Important part! Here segment(s) are actually chosen. All the good segments
1423   // could be chosen or some (best) ones only (in order to save time).
1424 
1425   // Check if there is a segment with n-1 hits that has a signifcantly better
1426   // weight than the best n hit segment
1427 
1428   // IBL 070828: implicit assumption here is that there will always only be one segment per
1429   // cluster - if there are >1 we will need to find out which segment the alternative n-1 hit
1430   // protosegment belongs to!
1431 
1432   //float chosen_weight = min_weight_A;
1433   //float chosen_ywgt = best_weight_B;
1434   //float chosen_curv = best_curv_A;
1435   //int chosen_nlayers = n_layers_occupied_tot;
1436   int chosen_pseg = best_pseg;
1437   if (best_pseg < 0) {
1438     return segmentInChamber;
1439   }
1440   chosen_Psegments = (Psegments);
1441   chosen_weight_A = (weight_A);
1442 
1443   float hit_drop_limit = -999999.999;
1444 
1445   // define different weight improvement requirements depending on how many layers are in the segment candidate
1446   switch (n_layers_processed) {
1447     case 1:
1448       // do nothing;
1449       break;
1450     case 2:
1451       // do nothing;
1452       break;
1453     case 3:
1454       // do nothing;
1455       break;
1456     case 4:
1457       hit_drop_limit = hitDropLimit6Hits * (1. / 2.) * hitDropLimit4Hits;
1458       if ((best_Layer_noLx < 1) || (best_Layer_noLx > 4)) {
1459         //      std::cout<<"CSCSegAlgoST: For four layers, best_Layer_noLx = "<< best_Layer_noLx << std::endl;
1460       }
1461       if ((best_Layer_noLx == 2) || (best_Layer_noLx == 3))
1462         hit_drop_limit = hit_drop_limit * (1. / 2.);
1463       break;
1464     case 5:
1465       hit_drop_limit = hitDropLimit6Hits * (2. / 3.) * hitDropLimit5Hits;
1466       if ((best_Layer_noLx < 1) || (best_Layer_noLx > 5)) {
1467         //      std::cout<<"CSCSegAlgoST: For five layers, best_Layer_noLx = "<< best_Layer_noLx << std::endl;
1468       }
1469       if ((best_Layer_noLx == 2) || (best_Layer_noLx == 4))
1470         hit_drop_limit = hit_drop_limit * (1. / 2.);
1471       if (best_Layer_noLx == 3)
1472         hit_drop_limit = hit_drop_limit * (1. / 3.);
1473       break;
1474     case 6:
1475       hit_drop_limit = hitDropLimit6Hits * (3. / 4.);
1476       if ((best_Layer_noLx < 1) || (best_Layer_noLx > 6)) {
1477         //      std::cout<<"CSCSegAlgoST: For six layers, best_Layer_noLx = "<< best_Layer_noLx << std::endl;
1478       }
1479       if ((best_Layer_noLx == 2) || (best_Layer_noLx == 5))
1480         hit_drop_limit = hit_drop_limit * (1. / 2.);
1481       if ((best_Layer_noLx == 3) || (best_Layer_noLx == 4))
1482         hit_drop_limit = hit_drop_limit * (1. / 3.);
1483       break;
1484 
1485     default:
1486       // Fallback - should never occur.
1487       LogTrace("CSCSegment|CSC")
1488           << "[CSCSegAlgoST::buildSegments] Unexpected number of layers with hits - please inform CSC DPG.";
1489       hit_drop_limit = 0.1;
1490   }
1491 
1492   // choose the NoLx collection (the one that contains the best N-1 candidate)
1493   switch (best_Layer_noLx) {
1494     case 1:
1495       Psegments_noLx.clear();
1496       Psegments_noLx = Psegments_noL1;
1497       weight_noLx_A.clear();
1498       weight_noLx_A = weight_noL1_A;
1499       break;
1500     case 2:
1501       Psegments_noLx.clear();
1502       Psegments_noLx = Psegments_noL2;
1503       weight_noLx_A.clear();
1504       weight_noLx_A = weight_noL2_A;
1505       break;
1506     case 3:
1507       Psegments_noLx.clear();
1508       Psegments_noLx = Psegments_noL3;
1509       weight_noLx_A.clear();
1510       weight_noLx_A = weight_noL3_A;
1511       break;
1512     case 4:
1513       Psegments_noLx.clear();
1514       Psegments_noLx = Psegments_noL4;
1515       weight_noLx_A.clear();
1516       weight_noLx_A = weight_noL4_A;
1517       break;
1518     case 5:
1519       Psegments_noLx.clear();
1520       Psegments_noLx = Psegments_noL5;
1521       weight_noLx_A.clear();
1522       weight_noLx_A = weight_noL5_A;
1523       break;
1524     case 6:
1525       Psegments_noLx.clear();
1526       Psegments_noLx = Psegments_noL6;
1527       weight_noLx_A.clear();
1528       weight_noLx_A = weight_noL6_A;
1529       break;
1530 
1531     default:
1532       // Fallback - should occur only for preclusters with only 3 layers with hits.
1533       Psegments_noLx.clear();
1534       weight_noLx_A.clear();
1535   }
1536 
1537   if (min_weight_A > 0.) {
1538     if (min_weight_noLx_A / min_weight_A < hit_drop_limit) {
1539       //chosen_weight = min_weight_noLx_A;
1540       //chosen_ywgt = best_weight_noLx_B;
1541       //chosen_curv = best_curv_noLx_A;
1542       //chosen_nlayers = n_layers_occupied_tot-1;
1543       chosen_pseg = best_noLx_pseg;
1544       chosen_Psegments.clear();
1545       chosen_weight_A.clear();
1546       chosen_Psegments = (Psegments_noLx);
1547       chosen_weight_A = (weight_noLx_A);
1548     }
1549   }
1550 
1551   if (onlyBestSegment) {
1552     ChooseSegments2a(chosen_Psegments, chosen_pseg);
1553   } else {
1554     ChooseSegments3(chosen_Psegments, chosen_weight_A, chosen_pseg);
1555   }
1556 
1557   for (unsigned int iSegment = 0; iSegment < GoodSegments.size(); ++iSegment) {
1558     protoSegment = GoodSegments[iSegment];
1559 
1560     // Create new fitter object
1561     CSCCondSegFit* segfit = new CSCCondSegFit(pset(), chamber(), protoSegment);
1562     condpass1 = false;
1563     condpass2 = false;
1564     segfit->setScaleXError(1.0);
1565     segfit->fit(condpass1, condpass2);
1566 
1567     // Attempt to handle numerical instability of the fit.
1568     // (Any segment with chi2/dof > chi2Norm_3D_ is considered
1569     // as potentially suffering from numerical instability in fit.)
1570     if (adjustCovariance()) {
1571       // Call the fit with prefitting option:
1572       // First fit a straight line to X-Z coordinates and calculate chi2
1573       // This is done in CSCCondSegFit::correctTheCovX()
1574       // Scale up errors in X if this chi2 is too big (default 'too big' is >20);
1575       // Then refit XY-Z with the scaled-up X errors
1576       if (segfit->chi2() / segfit->ndof() > chi2Norm_3D_) {
1577         condpass1 = true;
1578         segfit->fit(condpass1, condpass2);
1579       }
1580       if (segfit->scaleXError() < 1.00005) {
1581         LogTrace("CSCWeirdSegment") << "[CSCSegAlgoST::buildSegments] Segment ErrXX scaled and refit " << std::endl;
1582         if (segfit->chi2() / segfit->ndof() > chi2Norm_3D_) {
1583           // Call the fit with direct adjustment of condition number;
1584           // If the attempt to improve fit by scaling up X error fails
1585           // call the procedure to make the condition number of M compatible with
1586           // the precision of X and Y measurements;
1587           // Achieved by decreasing abs value of the Covariance
1588           LogTrace("CSCWeirdSegment")
1589               << "[CSCSegAlgoST::buildSegments] Segment ErrXY changed to match cond. number and refit " << std::endl;
1590           condpass2 = true;
1591           segfit->fit(condpass1, condpass2);
1592         }
1593       }
1594       // Call the pre-pruning procedure;
1595       // If the attempt to improve fit by scaling up X error is successfull,
1596       // while scale factor for X errors is too big.
1597       // Prune the recHit inducing the biggest contribution into X-Z chi^2
1598       // and refit;
1599       if (prePrun_ && (sqrt(segfit->scaleXError()) > prePrunLimit_) && (segfit->nhits() > 3)) {
1600         LogTrace("CSCWeirdSegment")
1601             << "[CSCSegAlgoST::buildSegments] Scale factor chi2uCorrection too big, pre-Prune and refit " << std::endl;
1602         protoSegment.erase(protoSegment.begin() + segfit->worstHit(), protoSegment.begin() + segfit->worstHit() + 1);
1603 
1604         // Need to create new fitter object to repeat fit with fewer hits
1605         // Original code maintained current values of condpass1, condpass2, scaleXError - calc in CorrectTheCovX()
1606         //@@ DO THE SAME THING HERE, BUT IS THAT CORRECT?! It does make a difference.
1607         double tempcorr = segfit->scaleXError();  // save current value
1608         delete segfit;
1609         segfit = new CSCCondSegFit(pset(), chamber(), protoSegment);
1610         segfit->setScaleXError(tempcorr);  // reset to previous value (rather than init to 1)
1611         segfit->fit(condpass1, condpass2);
1612       }
1613     }
1614 
1615     // calculate covariance matrix
1616     //    AlgebraicSymMatrix temp2 = segfit->covarianceMatrix();
1617 
1618     // build an actual CSC segment
1619     CSCSegment temp(protoSegment, segfit->intercept(), segfit->localdir(), segfit->covarianceMatrix(), segfit->chi2());
1620     delete segfit;
1621 
1622     if (debug)
1623       dumpSegment(temp);
1624 
1625     segmentInChamber.push_back(temp);
1626   }
1627   return segmentInChamber;
1628 }
1629 
1630 void CSCSegAlgoST::ChooseSegments2a(std::vector<ChamberHitContainer>& chosen_segments, int chosen_seg) {
1631   // just return best segment
1632   GoodSegments.clear();
1633   GoodSegments.push_back(chosen_segments[chosen_seg]);
1634 }
1635 
1636 void CSCSegAlgoST::ChooseSegments3(std::vector<ChamberHitContainer>& chosen_segments,
1637                                    std::vector<float>& chosen_weight,
1638                                    int chosen_seg) {
1639   int SumCommonHits = 0;
1640   GoodSegments.clear();
1641   int nr_remaining_candidates;
1642   unsigned int nr_of_segment_candidates;
1643 
1644   nr_remaining_candidates = nr_of_segment_candidates = chosen_segments.size();
1645 
1646   // always select and return best protosegment:
1647   GoodSegments.push_back(chosen_segments[chosen_seg]);
1648 
1649   float chosen_weight_temp = 999999.;
1650   int chosen_seg_temp = -1;
1651 
1652   // try to find further segment candidates:
1653   while (nr_remaining_candidates > 0) {
1654     for (unsigned int iCand = 0; iCand < nr_of_segment_candidates; ++iCand) {
1655       //only compare current best to psegs that have not been marked bad:
1656       if (chosen_weight[iCand] < 0.)
1657         continue;
1658       SumCommonHits = 0;
1659 
1660       for (int ihits = 0; ihits < int(chosen_segments[iCand].size());
1661            ++ihits) {  // iCand and iiCand NEED to have same nr of hits! (always have by construction)
1662         if (chosen_segments[iCand][ihits] == chosen_segments[chosen_seg][ihits]) {
1663           ++SumCommonHits;
1664         }
1665       }
1666 
1667       //mark a pseg bad:
1668       if (SumCommonHits > 1) {  // needs to be a card; should be investigated first
1669         chosen_weight[iCand] = -1.;
1670         nr_remaining_candidates -= 1;
1671       } else {
1672         // save the protosegment with the smallest weight
1673         if (chosen_weight[iCand] < chosen_weight_temp) {
1674           chosen_weight_temp = chosen_weight[iCand];
1675           chosen_seg_temp = iCand;
1676         }
1677       }
1678     }
1679 
1680     if (chosen_seg_temp > -1)
1681       GoodSegments.push_back(chosen_segments[chosen_seg_temp]);
1682 
1683     chosen_seg = chosen_seg_temp;
1684     // re-initialze temporary best parameters
1685     chosen_weight_temp = 999999;
1686     chosen_seg_temp = -1;
1687   }
1688 }
1689 
1690 void CSCSegAlgoST::ChooseSegments2(int best_seg) {
1691   //  std::vector <int> CommonHits(6); // nice  concept :)
1692   std::vector<unsigned int> BadCandidate;
1693   int SumCommonHits = 0;
1694   GoodSegments.clear();
1695   BadCandidate.clear();
1696   for (unsigned int iCand = 0; iCand < Psegments.size(); ++iCand) {
1697     // skip here if segment was marked bad
1698     for (unsigned int iiCand = iCand + 1; iiCand < Psegments.size(); ++iiCand) {
1699       // skip here too if segment was marked bad
1700       SumCommonHits = 0;
1701       if (Psegments[iCand].size() != Psegments[iiCand].size()) {
1702         LogTrace("CSCSegment|CSC")
1703             << "[CSCSegAlgoST::ChooseSegments2] ALARM!! Should not happen - please inform CSC DPG";
1704       } else {
1705         for (int ihits = 0; ihits < int(Psegments[iCand].size());
1706              ++ihits) {  // iCand and iiCand NEED to have same nr of hits! (alsways have by construction)
1707           if (Psegments[iCand][ihits] == Psegments[iiCand][ihits]) {
1708             ++SumCommonHits;
1709           }
1710         }
1711       }
1712       if (SumCommonHits > 1) {
1713         if (weight_A[iCand] > weight_A[iiCand]) {  // use weight_A here
1714           BadCandidate.push_back(iCand);
1715           // rather mark segment bad by an array which is in sync with protosegments!! e.g. set weight = weight*1000 or have an addidional array or set it to weight *= -1
1716         } else {
1717           BadCandidate.push_back(iiCand);
1718           // rather mark segment bad by an array which is in sync with protosegments!! e.g. set weight = weight*1000 or have an addidional array or set it to weight *= -1
1719         }
1720       }
1721     }
1722   }
1723   bool discard;
1724   for (unsigned int isegm = 0; isegm < Psegments.size(); ++isegm) {
1725     // For best results another iteration/comparison over Psegments
1726     //should be applied here... It would make the program much slower.
1727     discard = false;
1728     for (unsigned int ibad = 0; ibad < BadCandidate.size(); ++ibad) {
1729       // can save this loop if we used an array in sync with Psegments!!!!
1730       if (isegm == BadCandidate[ibad]) {
1731         discard = true;
1732       }
1733     }
1734     if (!discard) {
1735       GoodSegments.push_back(Psegments[isegm]);
1736     }
1737   }
1738 }
1739 
1740 void CSCSegAlgoST::findDuplicates(std::vector<CSCSegment>& segments) {
1741   // this is intended for ME1/1a only - we have ghost segments because of the strips ganging
1742   // this function finds them (first the rechits by sharesInput() )
1743   // if a segment shares all the rechits with another segment it is a duplicate (even if
1744   // it has less rechits)
1745 
1746   for (std::vector<CSCSegment>::iterator it = segments.begin(); it != segments.end(); ++it) {
1747     std::vector<CSCSegment*> duplicateSegments;
1748     for (std::vector<CSCSegment>::iterator it2 = segments.begin(); it2 != segments.end(); ++it2) {
1749       //
1750       bool allShared = true;
1751       if (it != it2) {
1752         allShared = it->sharesRecHits(*it2);
1753       } else {
1754         allShared = false;
1755       }
1756       //
1757       if (allShared) {
1758         duplicateSegments.push_back(&(*it2));
1759       }
1760     }
1761     it->setDuplicateSegments(duplicateSegments);
1762   }
1763 }
1764 
1765 void CSCSegAlgoST::dumpSegment(const CSCSegment& seg) const {
1766   // Only called if pset value 'CSCDebug' is set in config
1767 
1768   edm::LogVerbatim("CSCSegment") << "CSCSegment in " << chamber()->id() << "\nlocal position = " << seg.localPosition()
1769                                  << "\nerror = " << seg.localPositionError()
1770                                  << "\nlocal direction = " << seg.localDirection()
1771                                  << "\nerror =" << seg.localDirectionError() << "\ncovariance matrix"
1772                                  << seg.parametersError() << "chi2/ndf = " << seg.chi2() << "/"
1773                                  << seg.degreesOfFreedom() << "\n#rechits = " << seg.specificRecHits().size()
1774                                  << "\ntime = " << seg.time();
1775 }